libMesh
inf_fe.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 // Local includes
21 #include "libmesh/libmesh_config.h"
22 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
23 #include "libmesh/inf_fe.h"
24 #include "libmesh/quadrature_gauss.h"
25 #include "libmesh/elem.h"
26 #include "libmesh/libmesh_logging.h"
27 
28 namespace libMesh
29 {
30 
31 
32 
33 // Constructor
34 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
36  FEBase (Dim, fet),
37 
38  _n_total_approx_sf (0),
39  _n_total_qp (0),
40 
41  // initialize the current_fe_type to all the same
42  // values as \p fet (since the FE families and coordinate
43  // map type should not change), but use an invalid order
44  // for the radial part (since this is the only order
45  // that may change!).
46  // the data structures like \p phi etc are not initialized
47  // through the constructor, but through reinit()
48  current_fe_type (FEType(fet.order,
49  fet.family,
51  fet.radial_family,
52  fet.inf_map))
53 
54 {
55  // Sanity checks
56  libmesh_assert_equal_to (T_radial, fe_type.radial_family);
57  libmesh_assert_equal_to (T_map, fe_type.inf_map);
58 
59  // build the base_fe object
60  if (Dim != 1)
61  base_fe.reset(FEBase::build(Dim-1, fet).release());
62 }
63 
64 
65 
66 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
68 {
69  libmesh_assert(q);
70  libmesh_assert(base_fe.get());
71 
72  const Order base_int_order = q->get_order();
73  const Order radial_int_order = static_cast<Order>(2 * (static_cast<unsigned int>(fe_type.radial_order.get_order()) + 1) +2);
74  const unsigned int qrule_dim = q->get_dim();
75 
76  if (Dim != 1)
77  {
78  // build a Dim-1 quadrature rule of the type that we received
79  base_qrule.reset(QBase::build(q->type(), qrule_dim-1, base_int_order).release());
80  base_fe->attach_quadrature_rule(base_qrule.get());
81  }
82 
83  // in radial direction, always use Gauss quadrature
84  radial_qrule.reset(new QGauss(1, radial_int_order));
85 
86  // currently not used. But maybe helpful to store the QBase *
87  // with which we initialized our own quadrature rules
88  qrule = q;
89 }
90 
91 
92 
93 
94 
95 template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
97 {
98  base_elem.reset(Base::build_elem(inf_elem));
99 }
100 
101 
102 
103 
104 
105 
106 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
108  const std::vector<Point> * const pts,
109  const std::vector<Real> * const weights)
110 {
111  libmesh_assert(base_fe.get());
112  libmesh_assert(inf_elem);
113 
114  // I don't understand infinite elements well enough to risk
115  // calculating too little. :-( RHS
116  this->calculate_phi = this->calculate_dphi = this->calculate_dphiref = true;
117  this->get_xyz();
118  this->determine_calculations();
119  base_fe->calculate_phi = base_fe->calculate_dphi = base_fe->calculate_dphiref = true;
120  base_fe->get_xyz();
121  base_fe->determine_calculations();
122 
123  if (pts == libmesh_nullptr)
124  {
125  libmesh_assert(base_fe->qrule);
126  libmesh_assert_equal_to (base_fe->qrule, base_qrule.get());
128 
129  bool init_shape_functions_required = false;
130 
131  // init the radial data fields only when the radial order changes
133  {
135 
136  // Watch out: this call to QBase->init() only works for
137  // current_fe_type = const! To allow variable Order,
138  // the init() of QBase has to be modified...
139  radial_qrule->init(EDGE2);
140 
141  // initialize the radial shape functions
142  this->init_radial_shape_functions(inf_elem);
143 
144  init_shape_functions_required=true;
145  }
146 
147 
148  bool update_base_elem_required=true;
149 
150  // update the type in accordance to the current cell
151  // and reinit if the cell type has changed or (as in
152  // the case of the hierarchics) the shape functions
153  // depend on the particular element and need a reinit
154  if ((Dim != 1) &&
155  ((this->get_type() != inf_elem->type()) ||
156  (base_fe->shapes_need_reinit())))
157  {
158  // store the new element type, update base_elem
159  // here. Through \p update_base_elem_required,
160  // remember whether it has to be updated (see below).
161  elem_type = inf_elem->type();
162  this->update_base_elem(inf_elem);
163  update_base_elem_required=false;
164 
165  // initialize the base quadrature rule for the new element
166  base_qrule->init(base_elem->type());
167 
168  // initialize the shape functions in the base
169  base_fe->init_base_shape_functions(base_fe->qrule->get_points(),
170  base_elem.get());
171 
172  // compute the shape functions and map functions of base_fe
173  // before using them later in combine_base_radial.
174  base_fe->_fe_map->compute_map (base_fe->dim, base_fe->qrule->get_weights(),
175  base_elem.get(), base_fe->calculate_d2phi);
176  base_fe->compute_shape_functions(base_elem.get(), base_fe->qrule->get_points());
177 
178  init_shape_functions_required=true;
179  }
180 
181 
182  // when either the radial or base part change,
183  // we have to init the whole fields
184  if (init_shape_functions_required)
185  this->init_shape_functions (radial_qrule->get_points(),
186  base_fe->qrule->get_points(),
187  inf_elem);
188 
189  // computing the distance only works when we have the current
190  // base_elem stored. This happens when fe_type is const,
191  // the inf_elem->type remains the same. Then we have to
192  // update the base elem _here_.
193  if (update_base_elem_required)
194  this->update_base_elem(inf_elem);
195 
196  // compute dist (depends on geometry, therefore has to be updated for
197  // each and every new element), throw radial and base part together
198  this->combine_base_radial (inf_elem);
199 
200  this->_fe_map->compute_map (this->dim, _total_qrule_weights, inf_elem, this->calculate_d2phi);
201 
202  // Compute the shape functions and the derivatives
203  // at all quadrature points.
204  this->compute_shape_functions (inf_elem,base_fe->qrule->get_points());
205  }
206 
207  else // if pts != libmesh_nullptr
208  {
209  // update the elem_type
210  elem_type = inf_elem->type();
211 
212  // We'll assume that pts is a tensor product mesh of points.
213  // That will handle the pts.size()==1 case that we care about
214  // right now, and it will generalize a bit, and it won't break
215  // the assumptions elsewhere in InfFE.
216  std::vector<Point> radial_pts;
217  for (std::size_t p=0; p != pts->size(); ++p)
218  {
219  Real radius = (*pts)[p](Dim-1);
220  if (radial_pts.size() && radial_pts[0](0) == radius)
221  break;
222  radial_pts.push_back(Point(radius));
223  }
224  const unsigned int radial_pts_size = radial_pts.size();
225  const unsigned int base_pts_size = pts->size() / radial_pts_size;
226  // If we're a tensor product we should have no remainder
227  libmesh_assert_equal_to
228  (base_pts_size * radial_pts_size, pts->size());
229 
230  std::vector<Point> base_pts;
231  base_pts.reserve(base_pts_size);
232  for (std::size_t p=0; p != pts->size(); p += radial_pts_size)
233  {
234  Point pt = (*pts)[p];
235  pt(Dim-1) = 0;
236  base_pts.push_back(pt);
237  }
238 
239  // init radial shapes
240  this->init_radial_shape_functions(inf_elem, &radial_pts);
241 
242  // update the base
243  this->update_base_elem(inf_elem);
244 
245  // the finite element on the ifem base
246  base_fe.reset(FEBase::build(Dim-1, this->fe_type).release());
247 
248  base_fe->calculate_phi = base_fe->calculate_dphi = base_fe->calculate_dphiref = true;
249  base_fe->get_xyz();
250  base_fe->determine_calculations();
251 
252  // init base shapes
253  base_fe->init_base_shape_functions(base_pts,
254  base_elem.get());
255 
256  // compute the shape functions and map functions of base_fe
257  // before using them later in combine_base_radial.
258 
259  if (weights)
260  {
261  base_fe->_fe_map->compute_map (base_fe->dim, *weights,
262  base_elem.get(), base_fe->calculate_d2phi);
263  }
264  else
265  {
266  std::vector<Real> dummy_weights (pts->size(), 1.);
267  base_fe->_fe_map->compute_map (base_fe->dim, dummy_weights,
268  base_elem.get(), base_fe->calculate_d2phi);
269  }
270 
271  base_fe->compute_shape_functions(base_elem.get(), *pts);
272 
273  this->init_shape_functions (radial_pts, base_pts, inf_elem);
274 
275  // combine the base and radial shapes
276  this->combine_base_radial (inf_elem);
277 
278  // weights
279  if (weights != libmesh_nullptr)
280  {
281  this->_fe_map->compute_map (this->dim, *weights, inf_elem, this->calculate_d2phi);
282  }
283  else
284  {
285  std::vector<Real> dummy_weights (pts->size(), 1.);
286  this->_fe_map->compute_map (this->dim, dummy_weights, inf_elem, this->calculate_d2phi);
287  }
288 
289  // finally compute the ifem shapes
290  this->compute_shape_functions (inf_elem,*pts);
291  }
292 
293 }
294 
295 
296 
297 
298 
299 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
300 void
303  const std::vector<Point> * radial_pts)
304 {
305  libmesh_assert(radial_qrule.get() || radial_pts);
306  libmesh_assert(inf_elem);
307 
308  // Start logging the radial shape function initialization
309  LOG_SCOPE("init_radial_shape_functions()", "InfFE");
310 
311  // initialize most of the things related to mapping
312 
313  // The order to use in the radial map (currently independent of the element type)
314  const Order radial_mapping_order = Radial::mapping_order();
315  const unsigned int n_radial_mapping_shape_functions = Radial::n_dofs(radial_mapping_order);
316 
317  // initialize most of the things related to physical approximation
318  const Order radial_approx_order = fe_type.radial_order;
319  const unsigned int n_radial_approx_shape_functions = Radial::n_dofs(radial_approx_order);
320 
321  const unsigned int n_radial_qp =
322  radial_pts ? radial_pts->size() : radial_qrule->n_points();
323  const std::vector<Point> & radial_qp =
324  radial_pts ? *radial_pts : radial_qrule->get_points();
325 
326  // resize the radial data fields
327 
328  // the radial polynomials (eval)
329  mode.resize (n_radial_approx_shape_functions);
330  dmodedv.resize (n_radial_approx_shape_functions);
331 
332  // the (1-v)/2 weight
333  som.resize (n_radial_qp);
334  dsomdv.resize (n_radial_qp);
335 
336  // the radial map
337  radial_map.resize (n_radial_mapping_shape_functions);
338  dradialdv_map.resize (n_radial_mapping_shape_functions);
339 
340 
341  for (unsigned int i=0; i<n_radial_mapping_shape_functions; i++)
342  {
343  radial_map[i].resize (n_radial_qp);
344  dradialdv_map[i].resize (n_radial_qp);
345  }
346 
347 
348  for (unsigned int i=0; i<n_radial_approx_shape_functions; i++)
349  {
350  mode[i].resize (n_radial_qp);
351  dmodedv[i].resize (n_radial_qp);
352  }
353 
354 
355  // compute scalar values at radial quadrature points
356  for (unsigned int p=0; p<n_radial_qp; p++)
357  {
358  som[p] = Radial::decay (radial_qp[p](0));
359  dsomdv[p] = Radial::decay_deriv (radial_qp[p](0));
360  }
361 
362 
363  // evaluate the mode shapes in radial direction at radial quadrature points
364  for (unsigned int i=0; i<n_radial_approx_shape_functions; i++)
365  for (unsigned int p=0; p<n_radial_qp; p++)
366  {
367  mode[i][p] = InfFE<Dim,T_radial,T_map>::eval (radial_qp[p](0), radial_approx_order, i);
368  dmodedv[i][p] = InfFE<Dim,T_radial,T_map>::eval_deriv (radial_qp[p](0), radial_approx_order, i);
369  }
370 
371 
372  // evaluate the mapping functions in radial direction at radial quadrature points
373  for (unsigned int i=0; i<n_radial_mapping_shape_functions; i++)
374  for (unsigned int p=0; p<n_radial_qp; p++)
375  {
376  radial_map[i][p] = InfFE<Dim,INFINITE_MAP,T_map>::eval (radial_qp[p](0), radial_mapping_order, i);
377  dradialdv_map[i][p] = InfFE<Dim,INFINITE_MAP,T_map>::eval_deriv (radial_qp[p](0), radial_mapping_order, i);
378  }
379 }
380 
381 
382 
383 
384 
385 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
386 void InfFE<Dim,T_radial,T_map>::init_shape_functions(const std::vector<Point> & radial_qp,
387  const std::vector<Point> & base_qp,
388  const Elem * inf_elem)
389 {
390  libmesh_assert(inf_elem);
391 
392  // Start logging the radial shape function initialization
393  LOG_SCOPE("init_shape_functions()", "InfFE");
394 
395  // fast access to some const ints for the radial data
396  const unsigned int n_radial_mapping_sf = cast_int<unsigned int>(radial_map.size());
397  const unsigned int n_radial_approx_sf = cast_int<unsigned int>(mode.size());
398  const unsigned int n_radial_qp = cast_int<unsigned int>(som.size());
399 
400 
401  // initialize most of the things related to mapping
402 
403  // The element type and order to use in the base map
404  const Order base_mapping_order = base_elem->default_order();
405  const ElemType base_mapping_elem_type = base_elem->type();
406 
407  // the number of base shape functions used to construct the map
408  // (Lagrange shape functions are used for mapping in the base)
409  unsigned int n_base_mapping_shape_functions = Base::n_base_mapping_sf(base_mapping_elem_type,
410  base_mapping_order);
411 
412  const unsigned int n_total_mapping_shape_functions =
413  n_radial_mapping_sf * n_base_mapping_shape_functions;
414 
415  // initialize most of the things related to physical approximation
416  unsigned int n_base_approx_shape_functions;
417  if (Dim > 1)
418  n_base_approx_shape_functions = base_fe->n_shape_functions();
419  else
420  n_base_approx_shape_functions = 1;
421 
422 
423  const unsigned int n_total_approx_shape_functions =
424  n_radial_approx_sf * n_base_approx_shape_functions;
425 
426  // update class member field
427  _n_total_approx_sf = n_total_approx_shape_functions;
428 
429 
430  // The number of the base quadrature points.
431  const unsigned int n_base_qp = base_qp.size();
432 
433  // The total number of quadrature points.
434  const unsigned int n_total_qp = n_radial_qp * n_base_qp;
435 
436 
437  // update class member field
438  _n_total_qp = n_total_qp;
439 
440 
441 
442  // initialize the node and shape numbering maps
443  {
444  // these vectors work as follows: the i-th entry stores
445  // the associated base/radial node number
446  _radial_node_index.resize(n_total_mapping_shape_functions);
447  _base_node_index.resize(n_total_mapping_shape_functions);
448 
449  // similar for the shapes: the i-th entry stores
450  // the associated base/radial shape number
451  _radial_shape_index.resize(n_total_approx_shape_functions);
452  _base_shape_index.resize(n_total_approx_shape_functions);
453 
454  const ElemType inf_elem_type = inf_elem->type();
455 
456  // fill the node index map
457  for (unsigned int n=0; n<n_total_mapping_shape_functions; n++)
458  {
459  compute_node_indices (inf_elem_type,
460  n,
461  _base_node_index[n],
462  _radial_node_index[n]);
463  libmesh_assert_less (_base_node_index[n], n_base_mapping_shape_functions);
464  libmesh_assert_less (_radial_node_index[n], n_radial_mapping_sf);
465  }
466 
467  // fill the shape index map
468  for (unsigned int n=0; n<n_total_approx_shape_functions; n++)
469  {
471  inf_elem_type,
472  n,
475  libmesh_assert_less (_base_shape_index[n], n_base_approx_shape_functions);
476  libmesh_assert_less (_radial_shape_index[n], n_radial_approx_sf);
477  }
478  }
479 
480  // resize the base data fields
481  dist.resize(n_base_mapping_shape_functions);
482 
483  // resize the total data fields
484 
485  // the phase term varies with xi, eta and zeta(v): store it for _all_ qp
486  //
487  // when computing the phase, we need the base approximations
488  // therefore, initialize the phase here, but evaluate it
489  // in combine_base_radial().
490  //
491  // the weight, though, is only needed at the radial quadrature points, n_radial_qp.
492  // but for a uniform interface to the protected data fields
493  // the weight data field (which are accessible from the outside) are expanded to n_total_qp.
494  weight.resize (n_total_qp);
495  dweightdv.resize (n_total_qp);
496  dweight.resize (n_total_qp);
497 
498  dphase.resize (n_total_qp);
499  dphasedxi.resize (n_total_qp);
500  dphasedeta.resize (n_total_qp);
501  dphasedzeta.resize (n_total_qp);
502 
503  // this vector contains the integration weights for the combined quadrature rule
504  _total_qrule_weights.resize(n_total_qp);
505 
506  // InfFE's data fields phi, dphi, dphidx, phi_map etc hold the _total_
507  // shape and mapping functions, respectively
508  {
509  phi.resize (n_total_approx_shape_functions);
510  dphi.resize (n_total_approx_shape_functions);
511  dphidx.resize (n_total_approx_shape_functions);
512  dphidy.resize (n_total_approx_shape_functions);
513  dphidz.resize (n_total_approx_shape_functions);
514  dphidxi.resize (n_total_approx_shape_functions);
515 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
516  libmesh_do_once(libMesh::err << "Second derivatives for Infinite elements"
517  << " are not yet implemented!"
518  << std::endl);
519 
520  d2phi.resize (n_total_approx_shape_functions);
521  d2phidx2.resize (n_total_approx_shape_functions);
522  d2phidxdy.resize (n_total_approx_shape_functions);
523  d2phidxdz.resize (n_total_approx_shape_functions);
524  d2phidy2.resize (n_total_approx_shape_functions);
525  d2phidydz.resize (n_total_approx_shape_functions);
526  d2phidz2.resize (n_total_approx_shape_functions);
527  d2phidxi2.resize (n_total_approx_shape_functions);
528 
529  if (Dim > 1)
530  {
531  d2phidxideta.resize(n_total_approx_shape_functions);
532  d2phideta2.resize(n_total_approx_shape_functions);
533  }
534 
535  if (Dim > 2)
536  {
537  d2phidetadzeta.resize(n_total_approx_shape_functions);
538  d2phidxidzeta.resize(n_total_approx_shape_functions);
539  d2phidzeta2.resize(n_total_approx_shape_functions);
540  }
541 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
542 
543  if (Dim > 1)
544  dphideta.resize(n_total_approx_shape_functions);
545 
546  if (Dim == 3)
547  dphidzeta.resize(n_total_approx_shape_functions);
548 
549  std::vector<std::vector<Real>> & phi_map = this->_fe_map->get_phi_map();
550  std::vector<std::vector<Real>> & dphidxi_map = this->_fe_map->get_dphidxi_map();
551 
552  phi_map.resize(n_total_mapping_shape_functions);
553  dphidxi_map.resize(n_total_mapping_shape_functions);
554 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
555  std::vector<std::vector<Real>> & d2phidxi2_map = this->_fe_map->get_d2phidxi2_map();
556  d2phidxi2_map.resize(n_total_mapping_shape_functions);
557 
558  if (Dim > 1)
559  {
560  std::vector<std::vector<Real>> & d2phidxideta_map = this->_fe_map->get_d2phidxideta_map();
561  std::vector<std::vector<Real>> & d2phideta2_map = this->_fe_map->get_d2phideta2_map();
562  d2phidxideta_map.resize(n_total_mapping_shape_functions);
563  d2phideta2_map.resize(n_total_mapping_shape_functions);
564  }
565 
566  if (Dim == 3)
567  {
568  std::vector<std::vector<Real>> & d2phidxidzeta_map = this->_fe_map->get_d2phidxidzeta_map();
569  std::vector<std::vector<Real>> & d2phidetadzeta_map = this->_fe_map->get_d2phidetadzeta_map();
570  std::vector<std::vector<Real>> & d2phidzeta2_map = this->_fe_map->get_d2phidzeta2_map();
571  d2phidxidzeta_map.resize(n_total_mapping_shape_functions);
572  d2phidetadzeta_map.resize(n_total_mapping_shape_functions);
573  d2phidzeta2_map.resize(n_total_mapping_shape_functions);
574  }
575 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
576 
577  if (Dim > 1)
578  {
579  std::vector<std::vector<Real>> & dphideta_map = this->_fe_map->get_dphideta_map();
580  dphideta_map.resize(n_total_mapping_shape_functions);
581  }
582 
583  if (Dim == 3)
584  {
585  std::vector<std::vector<Real>> & dphidzeta_map = this->_fe_map->get_dphidzeta_map();
586  dphidzeta_map.resize(n_total_mapping_shape_functions);
587  }
588  }
589 
590  // collect all the for loops, where inner vectors are
591  // resized to the appropriate number of quadrature points
592  {
593  for (unsigned int i=0; i<n_total_approx_shape_functions; i++)
594  {
595  phi[i].resize (n_total_qp);
596  dphi[i].resize (n_total_qp);
597  dphidx[i].resize (n_total_qp);
598  dphidy[i].resize (n_total_qp);
599  dphidz[i].resize (n_total_qp);
600  dphidxi[i].resize (n_total_qp);
601 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
602  d2phi[i].resize (n_total_qp);
603  d2phidx2[i].resize (n_total_qp);
604  d2phidxdy[i].resize (n_total_qp);
605  d2phidxdz[i].resize (n_total_qp);
606  d2phidy2[i].resize (n_total_qp);
607  d2phidydz[i].resize (n_total_qp);
608  d2phidy2[i].resize (n_total_qp);
609  d2phidxi2[i].resize (n_total_qp);
610 
611  if (Dim > 1)
612  {
613  d2phidxideta[i].resize (n_total_qp);
614  d2phideta2[i].resize (n_total_qp);
615  }
616  if (Dim > 2)
617  {
618  d2phidxidzeta[i].resize (n_total_qp);
619  d2phidetadzeta[i].resize (n_total_qp);
620  d2phidzeta2[i].resize (n_total_qp);
621  }
622 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
623 
624  if (Dim > 1)
625  dphideta[i].resize (n_total_qp);
626 
627  if (Dim == 3)
628  dphidzeta[i].resize (n_total_qp);
629 
630  }
631 
632  for (unsigned int i=0; i<n_total_mapping_shape_functions; i++)
633  {
634  std::vector<std::vector<Real>> & phi_map = this->_fe_map->get_phi_map();
635  std::vector<std::vector<Real>> & dphidxi_map = this->_fe_map->get_dphidxi_map();
636  phi_map[i].resize (n_total_qp);
637  dphidxi_map[i].resize (n_total_qp);
638 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
639  std::vector<std::vector<Real>> & d2phidxi2_map = this->_fe_map->get_d2phidxi2_map();
640  d2phidxi2_map[i].resize (n_total_qp);
641  if (Dim > 1)
642  {
643  std::vector<std::vector<Real>> & d2phidxideta_map = this->_fe_map->get_d2phidxideta_map();
644  std::vector<std::vector<Real>> & d2phideta2_map = this->_fe_map->get_d2phideta2_map();
645  d2phidxideta_map[i].resize (n_total_qp);
646  d2phideta2_map[i].resize (n_total_qp);
647  }
648 
649  if (Dim > 2)
650  {
651  std::vector<std::vector<Real>> & d2phidxidzeta_map = this->_fe_map->get_d2phidxidzeta_map();
652  std::vector<std::vector<Real>> & d2phidetadzeta_map = this->_fe_map->get_d2phidetadzeta_map();
653  std::vector<std::vector<Real>> & d2phidzeta2_map = this->_fe_map->get_d2phidzeta2_map();
654  d2phidxidzeta_map[i].resize (n_total_qp);
655  d2phidetadzeta_map[i].resize (n_total_qp);
656  d2phidzeta2_map[i].resize (n_total_qp);
657  }
658 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
659 
660  if (Dim > 1)
661  {
662  std::vector<std::vector<Real>> & dphideta_map = this->_fe_map->get_dphideta_map();
663  dphideta_map[i].resize (n_total_qp);
664  }
665 
666  if (Dim == 3)
667  {
668  std::vector<std::vector<Real>> & dphidzeta_map = this->_fe_map->get_dphidzeta_map();
669  dphidzeta_map[i].resize (n_total_qp);
670  }
671  }
672  }
673 
674 
675 
676  {
677  // (a) compute scalar values at _all_ quadrature points -- for uniform
678  // access from the outside to these fields
679  // (b) form a std::vector<Real> which contains the appropriate weights
680  // of the combined quadrature rule!
681  libmesh_assert_equal_to (radial_qp.size(), n_radial_qp);
682 
683  if (radial_qrule && base_qrule)
684  {
685  const std::vector<Real> & radial_qw = radial_qrule->get_weights();
686  const std::vector<Real> & base_qw = base_qrule->get_weights();
687  libmesh_assert_equal_to (radial_qw.size(), n_radial_qp);
688  libmesh_assert_equal_to (base_qw.size(), n_base_qp);
689 
690  for (unsigned int rp=0; rp<n_radial_qp; rp++)
691  for (unsigned int bp=0; bp<n_base_qp; bp++)
692  {
693  weight[bp + rp*n_base_qp] = Radial::D(radial_qp[rp](0));
694  dweightdv[bp + rp*n_base_qp] = Radial::D_deriv(radial_qp[rp](0));
695  _total_qrule_weights[bp + rp*n_base_qp] = radial_qw[rp] * base_qw[bp];
696  }
697  }
698  else
699  {
700  for (unsigned int rp=0; rp<n_radial_qp; rp++)
701  for (unsigned int bp=0; bp<n_base_qp; bp++)
702  {
703  weight[bp + rp*n_base_qp] = Radial::D(radial_qp[rp](0));
704  dweightdv[bp + rp*n_base_qp] = Radial::D_deriv(radial_qp[rp](0));
705  }
706  }
707  }
708 }
709 
710 
711 
712 
713 
714 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
716 {
717  libmesh_assert(inf_elem);
718  // at least check whether the base element type is correct.
719  // otherwise this version of computing dist would give problems
720  libmesh_assert_equal_to (base_elem->type(), Base::get_elem_type(inf_elem->type()));
721 
722  // Start logging the combination of radial and base parts
723  LOG_SCOPE("combine_base_radial()", "InfFE");
724 
725  // zero the phase, since it is to be summed up
726  std::fill (dphasedxi.begin(), dphasedxi.end(), 0.);
727  std::fill (dphasedeta.begin(), dphasedeta.end(), 0.);
728  std::fill (dphasedzeta.begin(), dphasedzeta.end(), 0.);
729 
730 
731  const unsigned int n_base_mapping_sf = cast_int<unsigned int>(dist.size());
732  const Point origin = inf_elem->origin();
733 
734  // for each new infinite element, compute the radial distances
735  for (unsigned int n=0; n<n_base_mapping_sf; n++)
736  dist[n] = Point(base_elem->point(n) - origin).norm();
737 
738 
739  switch (Dim)
740  {
741  // 1D
742  case 1:
743  {
744  libmesh_not_implemented();
745  break;
746  }
747 
748  // 2D
749  case 2:
750  {
751  libmesh_not_implemented();
752  break;
753  }
754 
755  // 3D
756  case 3:
757  {
758  // fast access to the approximation and mapping shapes of base_fe
759  const std::vector<std::vector<Real>> & S = base_fe->phi;
760  const std::vector<std::vector<Real>> & Ss = base_fe->dphidxi;
761  const std::vector<std::vector<Real>> & St = base_fe->dphideta;
762  const std::vector<std::vector<Real>> & S_map = (base_fe->get_fe_map()).get_phi_map();
763  const std::vector<std::vector<Real>> & Ss_map = (base_fe->get_fe_map()).get_dphidxi_map();
764  const std::vector<std::vector<Real>> & St_map = (base_fe->get_fe_map()).get_dphideta_map();
765 
766  const unsigned int n_radial_qp = som.size();
767  if (radial_qrule)
768  libmesh_assert_equal_to(n_radial_qp, radial_qrule->n_points());
769  const unsigned int n_base_qp = S_map[0].size();
770  if (base_qrule)
771  libmesh_assert_equal_to(n_base_qp, base_qrule->n_points());
772 
773  const unsigned int n_total_mapping_sf =
774  cast_int<unsigned int>(radial_map.size()) * n_base_mapping_sf;
775 
776  const unsigned int n_total_approx_sf = Radial::n_dofs(fe_type.radial_order) * base_fe->n_shape_functions();
777 
778 
779  // compute the phase term derivatives
780  {
781  unsigned int tp=0;
782  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qps
783  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qps
784  {
785  // sum over all base shapes, to get the average distance
786  for (unsigned int i=0; i<n_base_mapping_sf; i++)
787  {
788  dphasedxi[tp] += Ss_map[i][bp] * dist[i] * radial_map [1][rp];
789  dphasedeta[tp] += St_map[i][bp] * dist[i] * radial_map [1][rp];
790  dphasedzeta[tp] += S_map [i][bp] * dist[i] * dradialdv_map[1][rp];
791  }
792 
793  tp++;
794 
795  } // loop radial and base qps
796  }
797 
798  libmesh_assert_equal_to (phi.size(), n_total_approx_sf);
799  libmesh_assert_equal_to (dphidxi.size(), n_total_approx_sf);
800  libmesh_assert_equal_to (dphideta.size(), n_total_approx_sf);
801  libmesh_assert_equal_to (dphidzeta.size(), n_total_approx_sf);
802 
803  // compute the overall approximation shape functions,
804  // pick the appropriate radial and base shapes through using
805  // _base_shape_index and _radial_shape_index
806  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qps
807  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qps
808  for (unsigned int ti=0; ti<n_total_approx_sf; ti++) // over _all_ approx_sf
809  {
810  // let the index vectors take care of selecting the appropriate base/radial shape
811  const unsigned int bi = _base_shape_index [ti];
812  const unsigned int ri = _radial_shape_index[ti];
813  phi [ti][bp+rp*n_base_qp] = S [bi][bp] * mode[ri][rp] * som[rp];
814  dphidxi [ti][bp+rp*n_base_qp] = Ss[bi][bp] * mode[ri][rp] * som[rp];
815  dphideta [ti][bp+rp*n_base_qp] = St[bi][bp] * mode[ri][rp] * som[rp];
816  dphidzeta[ti][bp+rp*n_base_qp] = S [bi][bp]
817  * (dmodedv[ri][rp] * som[rp] + mode[ri][rp] * dsomdv[rp]);
818  }
819 
820  std::vector<std::vector<Real>> & phi_map = this->_fe_map->get_phi_map();
821  std::vector<std::vector<Real>> & dphidxi_map = this->_fe_map->get_dphidxi_map();
822  std::vector<std::vector<Real>> & dphideta_map = this->_fe_map->get_dphideta_map();
823  std::vector<std::vector<Real>> & dphidzeta_map = this->_fe_map->get_dphidzeta_map();
824 
825  libmesh_assert_equal_to (phi_map.size(), n_total_mapping_sf);
826  libmesh_assert_equal_to (dphidxi_map.size(), n_total_mapping_sf);
827  libmesh_assert_equal_to (dphideta_map.size(), n_total_mapping_sf);
828  libmesh_assert_equal_to (dphidzeta_map.size(), n_total_mapping_sf);
829 
830  // compute the overall mapping functions,
831  // pick the appropriate radial and base entries through using
832  // _base_node_index and _radial_node_index
833  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qps
834  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qps
835  for (unsigned int ti=0; ti<n_total_mapping_sf; ti++) // over all mapping shapes
836  {
837  // let the index vectors take care of selecting the appropriate base/radial mapping shape
838  const unsigned int bi = _base_node_index [ti];
839  const unsigned int ri = _radial_node_index[ti];
840  phi_map [ti][bp+rp*n_base_qp] = S_map [bi][bp] * radial_map [ri][rp];
841  dphidxi_map [ti][bp+rp*n_base_qp] = Ss_map[bi][bp] * radial_map [ri][rp];
842  dphideta_map [ti][bp+rp*n_base_qp] = St_map[bi][bp] * radial_map [ri][rp];
843  dphidzeta_map[ti][bp+rp*n_base_qp] = S_map [bi][bp] * dradialdv_map[ri][rp];
844  }
845 
846  break;
847  }
848 
849  default:
850  libmesh_error_msg("Unsupported Dim = " << Dim);
851  }
852 }
853 
854 
855 
856 
857 
858 
859 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
861  const std::vector<Point> &)
862 {
863  // Start logging the overall computation of shape functions
864  LOG_SCOPE("compute_shape_functions()", "InfFE");
865 
866  const unsigned int n_total_qp = _n_total_qp;
867 
868  // Compute the shape function values (and derivatives)
869  // at the Quadrature points. Note that the actual values
870  // have already been computed via init_shape_functions
871 
872  // Compute the value of the derivative shape function i at quadrature point p
873  switch (dim)
874  {
875 
876  case 1:
877  {
878  libmesh_not_implemented();
879  break;
880  }
881 
882  case 2:
883  {
884  libmesh_not_implemented();
885  break;
886  }
887 
888  case 3:
889  {
890  const std::vector<Real> & dxidx_map = this->_fe_map->get_dxidx();
891  const std::vector<Real> & dxidy_map = this->_fe_map->get_dxidy();
892  const std::vector<Real> & dxidz_map = this->_fe_map->get_dxidz();
893 
894  const std::vector<Real> & detadx_map = this->_fe_map->get_detadx();
895  const std::vector<Real> & detady_map = this->_fe_map->get_detady();
896  const std::vector<Real> & detadz_map = this->_fe_map->get_detadz();
897 
898  const std::vector<Real> & dzetadx_map = this->_fe_map->get_dzetadx();
899  const std::vector<Real> & dzetady_map = this->_fe_map->get_dzetady();
900  const std::vector<Real> & dzetadz_map = this->_fe_map->get_dzetadz();
901 
902  // These are _all_ shape functions of this infinite element
903  for (std::size_t i=0; i<phi.size(); i++)
904  for (unsigned int p=0; p<n_total_qp; p++)
905  {
906  // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
907  dphi[i][p](0) =
908  dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
909  dphideta[i][p]*detadx_map[p] +
910  dphidzeta[i][p]*dzetadx_map[p]);
911 
912  // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
913  dphi[i][p](1) =
914  dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
915  dphideta[i][p]*detady_map[p] +
916  dphidzeta[i][p]*dzetady_map[p]);
917 
918  // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
919  dphi[i][p](2) =
920  dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
921  dphideta[i][p]*detadz_map[p] +
922  dphidzeta[i][p]*dzetadz_map[p]);
923  }
924 
925 
926  // This is the derivative of the phase term of this infinite element
927  for (unsigned int p=0; p<n_total_qp; p++)
928  {
929  // the derivative of the phase term
930  dphase[p](0) = (dphasedxi[p] * dxidx_map[p] +
931  dphasedeta[p] * detadx_map[p] +
932  dphasedzeta[p] * dzetadx_map[p]);
933 
934  dphase[p](1) = (dphasedxi[p] * dxidy_map[p] +
935  dphasedeta[p] * detady_map[p] +
936  dphasedzeta[p] * dzetady_map[p]);
937 
938  dphase[p](2) = (dphasedxi[p] * dxidz_map[p] +
939  dphasedeta[p] * detadz_map[p] +
940  dphasedzeta[p] * dzetadz_map[p]);
941 
942  // the derivative of the radial weight - varies only in radial direction,
943  // therefore dweightdxi = dweightdeta = 0.
944  dweight[p](0) = dweightdv[p] * dzetadx_map[p];
945 
946  dweight[p](1) = dweightdv[p] * dzetady_map[p];
947 
948  dweight[p](2) = dweightdv[p] * dzetadz_map[p];
949  }
950 
951  break;
952  }
953 
954  default:
955  libmesh_error_msg("Unsupported dim = " << dim);
956  }
957 }
958 
959 
960 
961 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
963 {
964  return false;
965 }
966 
967 } // namespace libMesh
968 
969 
970 // Explicit instantiations
971 #include "libmesh/inf_fe_instantiate_1D.h"
972 #include "libmesh/inf_fe_instantiate_2D.h"
973 #include "libmesh/inf_fe_instantiate_3D.h"
974 
975 
976 
977 #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
bool calculate_d2phi
Should we calculate shape function hessians?
Definition: fe_abstract.h:544
std::vector< std::vector< OutputTensor > > d2phi
Shape function second derivative values.
Definition: fe_base.h:552
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
OStreamProxy err
std::vector< std::vector< OutputShape > > dphidxi
Shape function derivatives in the xi direction.
Definition: fe_base.h:519
int get_order() const
Explicitly request the order as an int.
Definition: fe_type.h:77
std::vector< std::vector< OutputShape > > d2phidxdz
Shape function second derivatives in the x-z direction.
Definition: fe_base.h:597
Order get_order() const
Definition: quadrature.h:189
std::vector< std::vector< OutputShape > > dphidzeta
Shape function derivatives in the zeta direction.
Definition: fe_base.h:529
UniquePtr< QBase > radial_qrule
The quadrature rule for the base element associated with the current infinite element.
Definition: inf_fe.h:757
static Real decay_deriv(const Real)
Definition: inf_fe.h:114
std::vector< Real > dsomdv
the first local derivative of the radial decay in local coordinates.
Definition: inf_fe.h:633
std::vector< std::vector< Real > > dradialdv_map
the first local derivative of the radial mapping shapes
Definition: inf_fe.h:656
std::vector< std::vector< OutputShape > > d2phidydz
Shape function second derivatives in the y-z direction.
Definition: fe_base.h:607
bool calculate_phi
Should we calculate shape functions?
Definition: fe_abstract.h:534
std::vector< std::vector< Real > > radial_map
the radial mapping shapes in local coordinates
Definition: inf_fe.h:650
static Real D(const Real v)
Definition: inf_fe.h:120
const Real radius
std::vector< Real > dphasedxi
the first local derivative (for 3D, the first in the base) of the phase term in local coordinates...
Definition: inf_fe.h:663
virtual ElemType type() const =0
void combine_base_radial(const Elem *inf_elem)
Combines the shape functions, which were formed in init_shape_functions(Elem *), with geometric data...
Definition: inf_fe.C:715
virtual bool shapes_need_reinit() const libmesh_override
Definition: inf_fe.C:962
static Real eval_deriv(Real v, Order o_radial, unsigned int i)
std::vector< Real > dweightdv
the additional radial weight in local coordinates, over all quadrature points.
Definition: inf_fe.h:619
std::vector< std::vector< OutputShape > > d2phidxideta
Shape function second derivatives in the xi-eta direction.
Definition: fe_base.h:562
std::vector< unsigned int > _base_node_index
The internal structure of the InfFE – tensor product of base element times radial nodes – has to be...
Definition: inf_fe.h:702
ElemType
Defines an enum for geometric element types.
OrderWrapper radial_order
The approximation order in the base of the infinite element.
Definition: fe_type.h:236
void init_shape_functions(const std::vector< Point > &radial_qp, const std::vector< Point > &base_qp, const Elem *inf_elem)
Initialize all the data fields like weight, mode, phi, dphidxi, dphideta, dphidzeta, etc.
Definition: inf_fe.C:386
unsigned int _n_total_qp
The total number of quadrature points for the current configuration.
Definition: inf_fe.h:739
static Elem * build_elem(const Elem *inf_elem)
Build the base element of an infinite element.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
void init_radial_shape_functions(const Elem *inf_elem, const std::vector< Point > *radial_pts=libmesh_nullptr)
Some of the member data only depend on the radial part of the infinite element.
Definition: inf_fe.C:302
const class libmesh_nullptr_t libmesh_nullptr
std::vector< Real > weight
Used for certain infinite element families: the additional radial weight in local coordinates...
Definition: fe_base.h:644
std::vector< std::vector< OutputShape > > d2phidx2
Shape function second derivatives in the x direction.
Definition: fe_base.h:587
The libMesh namespace provides an interface to certain functionality in the library.
static ElemType get_elem_type(const ElemType type)
static void compute_node_indices(const ElemType inf_elem_type, const unsigned int outer_node_index, unsigned int &base_node, unsigned int &radial_node)
Computes the indices in the base base_node and in radial direction radial_node (either 0 or 1) associ...
std::vector< std::vector< OutputShape > > dphidy
Shape function derivatives in the y direction.
Definition: fe_base.h:539
std::vector< std::vector< OutputShape > > d2phidy2
Shape function second derivatives in the y direction.
Definition: fe_base.h:602
std::vector< std::vector< Real > > mode
the radial approximation shapes in local coordinates Needed when setting up the overall shape functio...
Definition: inf_fe.h:639
UniquePtr< QBase > base_qrule
The quadrature rule for the base element associated with the current infinite element.
Definition: inf_fe.h:751
UniquePtr< Elem > base_elem
The base element associated with the current infinite element.
Definition: inf_fe.h:763
libmesh_assert(j)
std::vector< std::vector< OutputShape > > d2phidetadzeta
Shape function second derivatives in the eta-zeta direction.
Definition: fe_base.h:577
ElemType get_type() const
Definition: fe_abstract.h:410
std::vector< std::vector< OutputShape > > d2phidxidzeta
Shape function second derivatives in the xi-zeta direction.
Definition: fe_base.h:567
std::vector< std::vector< OutputShape > > d2phidxdy
Shape function second derivatives in the x-y direction.
Definition: fe_base.h:592
std::vector< std::vector< OutputShape > > dphidx
Shape function derivatives in the x direction.
Definition: fe_base.h:534
friend class InfFE
Make all InfFE<Dim,T_radial,T_map> classes friends of each other, so that the protected eval() may be...
Definition: inf_fe.h:817
std::vector< std::vector< OutputShape > > phi
Shape function values.
Definition: fe_base.h:499
std::vector< std::vector< OutputShape > > d2phideta2
Shape function second derivatives in the eta direction.
Definition: fe_base.h:572
const unsigned int dim
The dimensionality of the object.
Definition: fe_abstract.h:523
static unsigned int n_base_mapping_sf(const ElemType base_elem_type, const Order base_mapping_order)
static Real decay(const Real v)
Definition: inf_fe.h:826
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr) libmesh_override
This is at the core of this class.
Definition: inf_fe.C:107
std::vector< OutputGradient > dphase
Used for certain infinite element families: the first derivatives of the phase term in global coordin...
Definition: fe_base.h:630
std::vector< Real > dphasedzeta
the third local derivative (for 3D, the derivative in radial direction) of the phase term in local co...
Definition: inf_fe.h:677
static Real eval(Real v, Order o_radial, unsigned int i)
InfMapType inf_map
The coordinate mapping type of the infinite element.
Definition: fe_type.h:257
FEType current_fe_type
This FEType stores the characteristics for which the data structures phi, phi_map etc are currently i...
Definition: inf_fe.h:781
virtual void compute_shape_functions(const Elem *, const std::vector< Point > &) libmesh_override
After having updated the jacobian and the transformation from local to global coordinates in FEAbstra...
Definition: inf_fe.C:860
QBase * qrule
A pointer to the quadrature rule employed.
Definition: fe_abstract.h:584
PetscErrorCode Vec Mat libmesh_dbg_var(j)
std::vector< Real > dphasedeta
the second local derivative (for 3D, the second in the base) of the phase term in local coordinates...
Definition: inf_fe.h:670
static Real D_deriv(const Real v)
Definition: inf_fe.h:125
static unsigned int n_dofs(const Order o_radial)
Definition: inf_fe.h:147
virtual QuadratureType type() const =0
void determine_calculations()
Determine which values are to be calculated, for both the FE itself and for the FEMap.
Definition: fe_base.C:739
bool calculate_dphiref
Should we calculate reference shape function gradients?
Definition: fe_abstract.h:559
std::vector< std::vector< OutputGradient > > dphi
Shape function derivative values.
Definition: fe_base.h:504
unsigned int _n_total_approx_sf
The number of total approximation shape functions for the current configuration.
Definition: inf_fe.h:733
FEFamily radial_family
For InfFE, family contains the radial shape family, while base_family contains the approximation type...
Definition: fe_type.h:249
virtual Point origin() const
Definition: elem.h:1480
std::vector< unsigned int > _radial_node_index
The internal structure of the InfFE – tensor product of base element times radial nodes – has to be...
Definition: inf_fe.h:692
std::vector< Real > som
the radial decay in local coordinates.
Definition: inf_fe.h:628
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:517
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool calculate_dphi
Should we calculate shape function gradients?
Definition: fe_abstract.h:539
unsigned int get_dim() const
Definition: quadrature.h:122
static void compute_shape_indices(const FEType &fet, const ElemType inf_elem_type, const unsigned int i, unsigned int &base_shape, unsigned int &radial_shape)
Computes the indices of shape functions in the base base_shape and in radial direction radial_shape (...
std::vector< unsigned int > _radial_shape_index
The internal structure of the InfFE – tensor product of base element shapes times radial shapes – h...
Definition: inf_fe.h:712
std::vector< std::vector< OutputShape > > d2phidz2
Shape function second derivatives in the z direction.
Definition: fe_base.h:612
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
This class implements specific orders of Gauss quadrature.
UniquePtr< FEBase > base_fe
Have a FE<Dim-1,T_base> handy for base approximation.
Definition: inf_fe.h:771
static UniquePtr< 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.
std::vector< std::vector< OutputShape > > d2phidxi2
Shape function second derivatives in the xi direction.
Definition: fe_base.h:557
std::vector< unsigned int > _base_shape_index
The internal structure of the InfFE – tensor product of base element shapes times radial shapes – h...
Definition: inf_fe.h:722
FEType fe_type
The finite element type for this object.
Definition: fe_abstract.h:567
std::vector< Real > dist
the radial distance of the base nodes from the origin
Definition: inf_fe.h:611
static Order mapping_order()
Definition: inf_fe.h:131
ElemType elem_type
The element type the current data structures are set up for.
Definition: fe_abstract.h:573
std::vector< std::vector< Real > > dmodedv
the first local derivative of the radial approximation shapes.
Definition: inf_fe.h:645
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
The QBase class provides the basic functionality from which various quadrature rules can be derived...
Definition: quadrature.h:53
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:230
virtual void attach_quadrature_rule(QBase *q) libmesh_override
The use of quadrature rules with the InfFE class is somewhat different from the approach of the FE cl...
Definition: inf_fe.C:67
This class forms the foundation from which generic finite elements may be derived.
std::vector< Real > _total_qrule_weights
this vector contains the combined integration weights, so that FEAbstract::compute_map() can still be...
Definition: inf_fe.h:745
std::vector< std::vector< OutputShape > > d2phidzeta2
Shape function second derivatives in the zeta direction.
Definition: fe_base.h:582
std::vector< std::vector< OutputShape > > dphidz
Shape function derivatives in the z direction.
Definition: fe_base.h:544
void update_base_elem(const Elem *inf_elem)
Updates the protected member base_elem to the appropriate base element for the given inf_elem...
Definition: inf_fe.C:96
std::vector< std::vector< OutputShape > > dphideta
Shape function derivatives in the eta direction.
Definition: fe_base.h:524
std::vector< RealGradient > dweight
Used for certain infinite element families: the global derivative of the additional radial weight ...
Definition: fe_base.h:637
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.