libMesh
inf_fe_boundary.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 
22 
23 // Local includes
24 #include "libmesh/libmesh_config.h"
25 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
26 #include "libmesh/inf_fe.h"
27 #include "libmesh/inf_fe_macro.h"
28 #include "libmesh/quadrature.h"
29 #include "libmesh/enum_quadrature_type.h"
30 #include "libmesh/elem.h"
31 
32 namespace libMesh
33 {
34 
35 // Method for 2D, 3D -- see inf_fe_1D.C for a 1D version of this
36 template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
37 void InfFE<Dim,T_radial,T_base>::reinit(const Elem * inf_elem,
38  const unsigned int s,
39  const Real tolerance,
40  const std::vector<Point> * const pts,
41  const std::vector<Real> * const weights)
42 {
43  if (weights != nullptr)
44  libmesh_not_implemented_msg("ERROR: User-specified weights for infinite elements are not implemented!");
45 
46  if (pts != nullptr)
47  libmesh_not_implemented_msg("ERROR: User-specified points for infinite elements are not implemented!");
48 
49  // We don't do this for 1D elements!
50  libmesh_assert_not_equal_to (Dim, 1);
51 
52  libmesh_assert(inf_elem);
53  libmesh_assert(qrule);
54 
55  // We need to allow this as well...
56  // libmesh_assert_not_equal_to (s, 0);
57 
58 
59  // Build the side of interest
60  const std::unique_ptr<const Elem> side(inf_elem->build_side_ptr(s));
61 
62  // set the element type
63  elem_type = inf_elem->type();
64 
65  // eventually initialize radial quadrature rule
66  bool radial_qrule_initialized = false;
67 
68  // if we are working on the base-side, the radial function is constant.
69  // With this, we ensure that at least for base elements we reinitialize all quantities
70  // when we enter for the first time.
71  if (s == 0)
72  current_fe_type.radial_order = 0;
73 
74  if (current_fe_type.radial_order != fe_type.radial_order)
75  {
76  if (s > 0)
77  {
78  current_fe_type.radial_order = fe_type.radial_order;
79  radial_qrule->init(EDGE2, inf_elem->p_level());
80  }
81  else
82  {
83  // build a new 0-dimensional quadrature-rule:
84  radial_qrule=QBase::build(QGAUSS, 0, fe_type.radial_order);
85  radial_qrule->init(NODEELEM, 0);
86 
87  //the base_qrule is set up with dim-1, but apparently we need dim, so we replace it:
88  base_qrule=QBase::build(qrule->type(), side->dim(), qrule->get_order());
89 
90  //FIXME: Do I have to care about the order of my neighbours element?
91  //unsigned int side_p_level = elem->p_level();
92  //if (elem->neighbor_ptr(s) != nullptr)
93  // side_p_level = std::max(side_p_level, elem->neighbor_ptr(s)->p_level());
94  base_qrule->init(side->type(), side->p_level());
95  }
96  radial_qrule_initialized = true;
97  }
98 
99  // Initialize the face shape functions
100  if (this->get_type() != inf_elem->type() ||
101  base_fe->shapes_need_reinit() ||
102  radial_qrule_initialized)
103  this->init_face_shape_functions (qrule->get_points(), side.get());
104 
105 
106  // compute the face map
107  this->_fe_map->compute_face_map(this->dim, _total_qrule_weights, side.get());
108 
109  // make a copy of the Jacobian for integration
110  const std::vector<Real> JxW_int(this->_fe_map->get_JxW());
111 
112  // Find where the integration points are located on the
113  // full element.
114  std::vector<Point> qp;
115  this->inverse_map (inf_elem, this->_fe_map->get_xyz(), qp, tolerance);
116 
117  // just to ensure that we are working on the base and not, for numeric reasons,
118  // somewhere else...
119  if (s==0)
120  {
121  for (auto & p : qp)
122  {
123  p(Dim-1)=-1.;
124  }
125  }
126 
127  // compute the shape function and derivative values
128  // at the points qp
129  this->reinit (inf_elem, &qp);
130 
131  // copy back old data
132  this->_fe_map->get_JxW() = JxW_int;
133 
134 }
135 
136 
137 
138 // Method for 2D, 3D -- see inf_fe_1D.C for a 1D version of this
139 template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
140 void InfFE<Dim,T_radial,T_base>::edge_reinit(const Elem *,
141  const unsigned int,
142  const Real,
143  const std::vector<Point> * const pts,
144  const std::vector<Real> * const /*weights*/)
145 {
146  // We don't do this for 1D elements!
147  //libmesh_assert_not_equal_to (Dim, 1);
148  libmesh_not_implemented_msg("ERROR: Edge conditions for infinite elements not implemented!");
149 
150  if (pts != nullptr)
151  libmesh_not_implemented_msg("ERROR: User-specified points for infinite elements not implemented!");
152 }
153 
154 
155 
156 
157 template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
158 void InfFE<Dim,T_radial,T_base>::init_face_shape_functions(const std::vector<Point> &,
159  const Elem * inf_side)
160 {
161  libmesh_assert(inf_side);
162 
163  // Currently, this makes only sense in 3-D!
164  libmesh_assert_equal_to (Dim, 3);
165 
166  // Initialize the radial shape functions
167  this->init_radial_shape_functions(inf_side);
168 
169  // Initialize the base shape functions
170  if (inf_side->infinite())
171  this->update_base_elem(inf_side);
172  else
173  // in this case, I need the 2D base
174  this->update_base_elem(inf_side->parent());
175 
176  // Initialize the base quadrature rule
177  base_qrule->init(base_elem->type(), inf_side->p_level());
178 
179  // base_fe still corresponds to the (dim-1)-dimensional base of the InfFE object,
180  // so update the fe_base.
181  if (inf_side->infinite())
182  {
183  libmesh_assert_equal_to (Dim, 3);
184  base_fe = FEBase::build(Dim-2, this->fe_type);
185  base_fe->attach_quadrature_rule(base_qrule.get());
186  }
187  else
188  {
189  base_fe = FEBase::build(Dim-1, this->fe_type);
190  base_fe->attach_quadrature_rule(base_qrule.get());
191  }
192 
193  //before initializing, we should say what to compute:
194  base_fe->_fe_map->get_xyz();
195  base_fe->_fe_map->get_JxW();
196 
197  // initialize the shape functions on the base
198  base_fe->init_base_shape_functions(base_fe->qrule->get_points(),
199  base_elem.get());
200 
201  // the number of quadrature points
202  const unsigned int n_radial_qp =
203  cast_int<unsigned int>(som.size());
204  const unsigned int n_base_qp = base_qrule->n_points();
205  const unsigned int n_total_qp = n_radial_qp * n_base_qp;
206 
207  // the quadrature weights
208  _total_qrule_weights.resize(n_total_qp);
209 
210  // now init the shapes for boundary work
211  {
212 
213  // The element type and order to use in the base map
214  const Order base_mapping_order ( base_elem->default_order() );
215 
216  // the number of mapping shape functions. For base side it is 1.
217  // (Lagrange shape functions are used for mapping in the base)
218  const unsigned int n_radial_mapping_sf =
219  inf_side->infinite() ? cast_int<unsigned int>(radial_map.size()) : 1;
220 
221  const unsigned int n_base_mapping_shape_functions =
222  InfFEBase::n_base_mapping_sf(*base_elem, base_mapping_order);
223 
224  const unsigned int n_total_mapping_shape_functions =
225  n_radial_mapping_sf * n_base_mapping_shape_functions;
226 
227 
228  // initialize the node and shape numbering maps
229  _radial_node_index.resize (n_total_mapping_shape_functions);
230  _base_node_index.resize (n_total_mapping_shape_functions);
231 
232  if (inf_side->infinite())
233  {
234  const ElemType inf_face_elem_type (inf_side->type());
235 
236  // fill the node index map
237  for (unsigned int n=0; n<n_total_mapping_shape_functions; n++)
238  {
239  compute_node_indices (inf_face_elem_type,
240  n,
241  _base_node_index[n],
242  _radial_node_index[n]);
243 
244  libmesh_assert_less (_base_node_index[n], n_base_mapping_shape_functions);
245  libmesh_assert_less (_radial_node_index[n], n_radial_mapping_sf);
246  }
247  }
248  else
249  {
250  for (unsigned int n=0; n<n_total_mapping_shape_functions; n++)
251  {
252  _base_node_index[n] = n;
253  _radial_node_index[n] = 0;
254  }
255  }
256 
257  // resize map data fields
258  std::vector<std::vector<Real>> & psi_map = this->_fe_map->get_psi();
259  std::vector<std::vector<Real>> & dpsidxi_map = this->_fe_map->get_dpsidxi();
260  std::vector<std::vector<Real>> & dpsideta_map = this->_fe_map->get_dpsideta();
261  psi_map.resize (n_total_mapping_shape_functions);
262  dpsidxi_map.resize (n_total_mapping_shape_functions);
263  dpsideta_map.resize (n_total_mapping_shape_functions);
264 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
265  std::vector<std::vector<Real>> & d2psidxi2_map = this->_fe_map->get_d2psidxi2();
266  std::vector<std::vector<Real>> & d2psidxideta_map = this->_fe_map->get_d2psidxideta();
267  std::vector<std::vector<Real>> & d2psideta2_map = this->_fe_map->get_d2psideta2();
268  d2psidxi2_map.resize (n_total_mapping_shape_functions);
269  d2psidxideta_map.resize (n_total_mapping_shape_functions);
270  d2psideta2_map.resize (n_total_mapping_shape_functions);
271 #endif
272 
273  for (unsigned int i=0; i<n_total_mapping_shape_functions; i++)
274  {
275  psi_map[i].resize (n_total_qp);
276  dpsidxi_map[i].resize (n_total_qp);
277  dpsideta_map[i].resize (n_total_qp);
278 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
279  d2psidxi2_map[i].resize (n_total_qp);
280  d2psidxideta_map[i].resize(n_total_qp);
281  d2psideta2_map[i].resize (n_total_qp);
282 #endif
283  }
284 
285  // compute shape maps
286  if (inf_side->infinite())
287  {
288  const std::vector<std::vector<Real>> & S_map = (base_fe->get_fe_map()).get_phi_map();
289  const std::vector<std::vector<Real>> & Ss_map = (base_fe->get_fe_map()).get_dphidxi_map();
290 
291  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qps
292  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qps
293  for (unsigned int ti=0; ti<n_total_mapping_shape_functions; ti++) // over all mapping shapes
294  {
295  // let the index vectors take care of selecting the appropriate base/radial mapping shape
296  const unsigned int bi = _base_node_index [ti];
297  const unsigned int ri = _radial_node_index[ti];
298  psi_map [ti][bp+rp*n_base_qp] = S_map [bi][bp] * radial_map [ri][rp];
299  dpsidxi_map [ti][bp+rp*n_base_qp] = Ss_map[bi][bp] * radial_map [ri][rp];
300  dpsideta_map [ti][bp+rp*n_base_qp] = S_map [bi][bp] * dradialdv_map[ri][rp];
301 
302  // second derivatives are not implemented for infinite elements
303 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
304  d2psidxi2_map [ti][bp+rp*n_base_qp] = 0.;
305  d2psidxideta_map [ti][bp+rp*n_base_qp] = 0.;
306  d2psideta2_map [ti][bp+rp*n_base_qp] = 0.;
307 #endif
308  }
309 
310  }
311  else
312  {
313  const std::vector<std::vector<Real>> & S_map = (base_fe->get_fe_map()).get_phi_map();
314  const std::vector<std::vector<Real>> & Ss_map = (base_fe->get_fe_map()).get_dphidxi_map();
315  const std::vector<std::vector<Real>> & St_map = (base_fe->get_fe_map()).get_dphideta_map();
316  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qps
317  for (unsigned int ti=0; ti<n_total_mapping_shape_functions; ti++) // over all mapping shapes
318  {
319  psi_map [ti][bp] = S_map[ti][bp] ;
320  dpsidxi_map [ti][bp] = Ss_map[ti][bp] ;
321  dpsideta_map [ti][bp] = St_map[ti][bp] ;
322 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
323  d2psidxi2_map [ti][bp] = 0.;
324  d2psidxideta_map [ti][bp] = 0.;
325  d2psideta2_map [ti][bp] = 0.;
326 #endif
327  }
328  }
329  }
330 
331  // quadrature rule weights
332  {
333  const std::vector<Real> & radial_qw = radial_qrule->get_weights();
334  const std::vector<Real> & base_qw = base_qrule->get_weights();
335 
336  libmesh_assert_equal_to (radial_qw.size(), n_radial_qp);
337  libmesh_assert_equal_to (base_qw.size(), n_base_qp);
338 
339  for (unsigned int rp=0; rp<n_radial_qp; rp++)
340  for (unsigned int bp=0; bp<n_base_qp; bp++)
341  {
342  _total_qrule_weights[bp + rp*n_base_qp] = radial_qw[rp] * base_qw[bp];
343  }
344  }
345 
346 }
347 
348 
349 
350 
351 // Explicit instantiations - doesn't make sense in 1D, but as
352 // long as we only return errors, we are fine... ;-)
353 //#include "libmesh/inf_fe_instantiate_1D.h"
354 //#include "libmesh/inf_fe_instantiate_2D.h"
355 //#include "libmesh/inf_fe_instantiate_3D.h"
356 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
357 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
358 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
359 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, edge_reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
360 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, edge_reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
361 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, edge_reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
362 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, init_face_shape_functions(const std::vector<Point> &, const Elem *));
363 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, init_face_shape_functions(const std::vector<Point> &, const Elem *));
364 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, init_face_shape_functions(const std::vector<Point> &, const Elem *));
365 
366 } // namespace libMesh
367 
368 #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
libMesh::Elem::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy=true)=0
libMesh::InfFE::init_face_shape_functions
void init_face_shape_functions(const std::vector< Point > &, const Elem *inf_side)
Initialize all the data fields like weight, phi, etc for the side s.
Definition: inf_fe_boundary.C:158
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
libMesh::Order
Order
Definition: enum_order.h:40
libMesh::Elem::infinite
virtual bool infinite() const =0
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2512
libMesh::QGAUSS
Definition: enum_quadrature_type.h:34
dim
unsigned int dim
Definition: adaptivity_ex3.C:113
libMesh::QBase::build
static std::unique_ptr< QBase > build(const std::string &name, const unsigned int dim, const Order order=INVALID_ORDER)
Builds a specific quadrature rule based on the name string.
Definition: quadrature_build.C:41
libMesh::CARTESIAN
Definition: enum_inf_map_type.h:35
libMesh::INSTANTIATE_INF_FE_MBRF
INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, reinit(const Elem *, const unsigned int, const Real, const std::vector< Point > *const, const std::vector< Real > *const))
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::FEGenericBase::build
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libMesh::InfFE::reinit
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr) override
This is at the core of this class.
Definition: inf_fe.C:112
libMesh::Elem::parent
const Elem * parent() const
Definition: elem.h:2434
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:100
libMesh::NODEELEM
Definition: enum_elem_type.h:66
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::InfFE::edge_reinit
virtual void edge_reinit(const Elem *elem, const unsigned int edge, const Real tolerance=TOLERANCE, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr) override
Not implemented yet.
Definition: inf_fe_boundary.C:140
libMesh::EDGE2
Definition: enum_elem_type.h:35
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::InfFEBase::n_base_mapping_sf
static unsigned int n_base_mapping_sf(const Elem &base_elem, const Order base_mapping_order)
Definition: inf_fe_base_radial.C:95
Generated on Sat Jan 25 2020 12:06:59 for libMesh by   doxygen 1.8.16