libMesh
inf_fe.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 
23 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
24 
25 #include "libmesh/inf_fe.h"
26 #include "libmesh/fe.h"
27 #include "libmesh/quadrature_gauss.h"
28 #include "libmesh/elem.h"
29 #include "libmesh/libmesh_logging.h"
30 #include "libmesh/int_range.h"
31 #include "libmesh/type_tensor.h"
32 #include "libmesh/fe_interface.h"
33 
34 #include <memory>
35 
36 namespace libMesh
37 {
38 
39 
40 
41 // Constructor
42 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
44  FEBase (Dim, fet),
45 
46  calculate_map_scaled(false),
47  calculate_phi_scaled(false),
48  calculate_dphi_scaled(false),
49  calculate_xyz(false),
50  calculate_jxw(false),
51  _n_total_approx_sf (0),
52 
53  // initialize the current_fe_type to all the same
54  // values as \p fet (since the FE families and coordinate
55  // map type should not change), but use an invalid order
56  // for the radial part (since this is the only order
57  // that may change!).
58  // the data structures like \p phi etc are not initialized
59  // through the constructor, but through reinit()
60  current_fe_type (FEType(fet.order,
61  fet.family,
63  fet.radial_family,
64  fet.inf_map))
65 
66 {
67  // Sanity checks
68  libmesh_assert_equal_to (T_radial, fe_type.radial_family);
69  libmesh_assert_equal_to (T_map, fe_type.inf_map);
70 
71  // build the base_fe object
72  if (Dim != 1)
73  base_fe = FEBase::build(Dim-1, fet);
74 }
75 
76 
77 
78 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
80 {
81  libmesh_assert(q);
82  libmesh_assert(base_fe);
83 
84  const Order base_int_order = q->get_order();
85  const Order radial_int_order = static_cast<Order>(2 * (static_cast<unsigned int>(fe_type.radial_order.get_order()) + 1) +2);
86  const unsigned int qrule_dim = q->get_dim();
87 
88  if (Dim != 1)
89  {
90  // build a Dim-1 quadrature rule of the type that we received
91  base_qrule = QBase::build(q->type(), qrule_dim-1, base_int_order);
92  base_fe->attach_quadrature_rule(base_qrule.get());
93  }
94 
95  // in radial direction, always use Gauss quadrature
96  radial_qrule = std::make_unique<QGauss>(1, radial_int_order);
97 
98  // Maybe helpful to store the QBase *
99  // with which we initialized our own quadrature rules.
100  // Used e.g. in \p InfFE::reinit(elem,side)
101  qrule = q;
102 }
103 
104 
105 
106 
107 
108 template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
110 {
111  base_elem = InfFEBase::build_elem(inf_elem);
112 }
113 
114 
115 
116 
117 
118 
119 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
121  const std::vector<Point> * const pts,
122  const std::vector<Real> * const weights)
123 {
124  libmesh_assert(base_fe.get());
125  libmesh_assert(inf_elem);
126 
127  // checks for consistency of requested calculations,
128  // adds further quantities as needed.
129  this->determine_calculations();
130 
131  if (pts == nullptr)
132  {
133  libmesh_assert(base_fe->qrule);
134  libmesh_assert_equal_to (base_fe->qrule, base_qrule.get());
135  libmesh_assert(radial_qrule.get());
136 
137  bool init_shape_functions_required = false;
138 
139  // init the radial data fields only when the radial order changes
140  if (current_fe_type.radial_order != fe_type.radial_order)
141  {
142  current_fe_type.radial_order = fe_type.radial_order;
143 
144  // Watch out: this call to QBase->init() only works for
145  // current_fe_type = const! To allow variable Order,
146  // the init() of QBase has to be modified...
147  radial_qrule->init(EDGE2);
148 
149  // initialize the radial shape functions
150  this->init_radial_shape_functions(inf_elem);
151 
152  init_shape_functions_required=true;
153  }
154 
155 
156  bool update_base_elem_required=true;
157 
158  // update the type in accordance to the current cell
159  // and reinit if the cell type has changed or (as in
160  // the case of the hierarchics) the shape functions
161  // depend on the particular element and need a reinit
162  if ((Dim != 1) &&
163  ((this->get_type() != inf_elem->type()) ||
164  (base_fe->shapes_need_reinit())))
165  {
166  // store the new element type, update base_elem
167  // here. Through \p update_base_elem_required,
168  // remember whether it has to be updated (see below).
169  this->_elem = inf_elem;
170  this->_elem_type = inf_elem->type();
171  this->update_base_elem(inf_elem);
172  update_base_elem_required=false;
173 
174  // initialize the base quadrature rule for the new element
175  base_qrule->init(*base_elem);
176  init_shape_functions_required=true;
177 
178  }
179  else
180  this->_elem = inf_elem;
181 
182  // computing the reference-to-physical map and coordinates works
183  // only, if we have the current base_elem stored.
184  // This happens when fe_type is const,
185  // the inf_elem->type remains the same. Then we have to
186  // update the base elem _here_.
187  if (update_base_elem_required)
188  this->update_base_elem(inf_elem);
189 
190  if (calculate_phi_scaled || calculate_dphi_scaled || calculate_phi || calculate_dphi)
191  // initialize the shape functions in the base
192  base_fe->init_base_shape_functions(base_fe->qrule->get_points(),
193  base_elem.get());
194 
195  // compute the shape functions and map functions of base_fe
196  // before using them later in compute_shape_functions.
197  base_fe->_fe_map->compute_map (base_fe->dim, base_fe->qrule->get_weights(),
198  base_elem.get(), base_fe->calculate_d2phi);
199  base_fe->compute_shape_functions(base_elem.get(), base_fe->qrule->get_points());
200 
201  // when either the radial or base part change,
202  // we have to init the whole fields
203  if (init_shape_functions_required)
204  this->init_shape_functions (radial_qrule->get_points(),
205  base_fe->qrule->get_points(),
206  inf_elem);
207 
208  // Compute the shape functions and the derivatives
209  // at all quadrature points.
210  this->compute_shape_functions (inf_elem,
211  base_fe->qrule->get_points(),
212  radial_qrule->get_points()
213  /* weights are computed inside the function*/
214  );
215  }
216 
217  else // if pts != nullptr
218  {
219  // update the elem
220  this->_elem = inf_elem;
221  this->_elem_type = inf_elem->type();
222 
223  // We'll assume that pts is a tensor product mesh of points.
224  // pts[i] = pts[ angular_index + n_angular_pts * radial_index]
225  // That will handle the pts.size()==1 case that we care about
226  // right now, and it will generalize a bit, and it won't break
227  // the assumptions elsewhere in InfFE.
228  std::vector<Point> radial_pts;
229  if (pts->size() > 0)
230  {
231  Real radius = (*pts)[0](Dim-1);
232  radial_pts.push_back(radius);
233  unsigned int n_radial_pts=1;
234  unsigned int n_angular_pts=1;
235  for (auto p : IntRange<std::size_t>(1, pts->size()))
236  {
237  radius = (*pts)[p](Dim-1);
238  // check for changes of radius: The max. allowed distance is somewhat arbitrary
239  // but the given value should not produce false positives...
240  if (std::abs(radial_pts[n_radial_pts-1](0) - radius) > 1e-4)
241  {
242  // it may change only every n_angular_pts:
243  if (p == (n_radial_pts)*n_angular_pts)
244  {
245  radial_pts.push_back(radius);
246  ++n_radial_pts;
247  }
248  else
249  {
250  libmesh_error_msg("We assumed that the "<<pts->size()
251  <<" points are of tensor-product type with "
252  <<n_radial_pts<<" radial points and "
253  <<n_angular_pts<< " angular points."<<std::endl
254  <<"But apparently point "<<p+1
255  <<" does not fit that scheme: Its radius is "
256  <<radius <<"but should have "
257  <<radial_pts[n_radial_pts*n_angular_pts-p]<<".");
258  //<<radial_pts[p-n_radial_pts*n_angular_pts]<<".");
259  }
260  }
261  // if we are still at the first radial segment,
262  // we consider another angular point
263  else if (n_radial_pts == 1)
264  {
265  ++n_angular_pts;
266  }
267  // if there was repetition but this does not, the assumed
268  // format does not work:
269  }
270  }
271  else
272  {
273  // I don't see any reason to call this function with no points.
274  libmesh_error_msg("Calling reinit() with an empty point list is prohibited.\n");
275  }
276 
277  const std::size_t radial_pts_size = radial_pts.size();
278  const std::size_t base_pts_size = pts->size() / radial_pts_size;
279  // If we're a tensor product we should have no remainder
280  libmesh_assert_equal_to
281  (base_pts_size * radial_pts_size, pts->size());
282 
283 
284  std::vector<Point> base_pts;
285  base_pts.reserve(base_pts_size);
286  for (std::size_t p=0; p != base_pts_size; ++p)
287  {
288  Point pt = (*pts)[p];
289  pt(Dim-1) = 0;
290  base_pts.push_back(pt);
291  }
292 
293  // init radial shapes
294  this->init_radial_shape_functions(inf_elem, &radial_pts);
295 
296  // update the base
297  this->update_base_elem(inf_elem);
298 
299  // the finite element on the ifem base
300  base_fe = FEBase::build(Dim-1, this->fe_type);
301 
302  // having a new base_fe, we need to redetermine the tasks...
303  this->determine_calculations();
304 
305  base_fe->reinit( base_elem.get(), &base_pts);
306 
307  this->init_shape_functions (radial_pts, base_pts, inf_elem);
308 
309  // finally compute the ifem shapes
310  if (weights != nullptr)
311  {
312  this->compute_shape_functions (inf_elem,base_pts,radial_pts);
313  }
314  else
315  {
316  this->compute_shape_functions (inf_elem, base_pts, radial_pts);
317  }
318 
319  }
320 
321  if (this->calculate_dual)
322  libmesh_not_implemented_msg("Dual shape support for infinite elements is "
323  "not currently implemented");
324 }
325 
326 
327 
328 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
330 {
331  this->calculations_started = true;
332 
333  // If the user forgot to request anything, but we're enabling
334  // deprecated backwards compatibility, then we'll be safe and
335  // calculate everything. If we haven't enable deprecated backwards
336  // compatibility then we'll scream and die.
337 #ifdef LIBMESH_ENABLE_DEPRECATED
338  if (!this->calculate_nothing &&
339  !this->calculate_phi && !this->calculate_dphi &&
340  !this->calculate_dphiref &&
341  !this->calculate_phi_scaled && !this->calculate_dphi_scaled &&
342  !this->calculate_xyz && !this->calculate_jxw &&
343  !this->calculate_map_scaled && !this->calculate_map &&
344 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
345  !this->calculate_d2phi &&
346 #endif
347  !this->calculate_curl_phi && !this->calculate_div_phi)
348  {
349  libmesh_deprecated();
350  this->calculate_phi = this->calculate_dphi = this->calculate_jxw = true;
351  this->calculate_dphiref = true;
352 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
353  this->calculate_d2phi = true;
354 #endif
355  this->calculate_phi_scaled = this->calculate_dphi_scaled = this->calculate_xyz = true;
356  if (FEInterface::field_type(fe_type.family) == TYPE_VECTOR)
357  {
358  this->calculate_curl_phi = true;
359  this->calculate_div_phi = true;
360  }
361  }
362 #else //LIBMESH_ENABLE_DEPRECATED
363  // ANSI C does not allow to embed the preprocessor-statement into the assert, so we
364  // make two statements, just different by 'calculate_d2phi'.
365 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
366  libmesh_assert (this->calculate_nothing || this->calculate_d2phi ||
367  this->calculate_phi || this->calculate_dphi ||
368  this->calculate_dphiref ||
369  this->calculate_phi_scaled || this->calculate_dphi_scaled ||
370  this->calculate_xyz || this->calculate_jxw ||
371  this->calculate_map_scaled || this->calculate_map ||
372  this->calculate_curl_phi || this->calculate_div_phi);
373 #else
374  libmesh_assert (this->calculate_nothing ||
375  this->calculate_phi || this->calculate_dphi ||
376  this->calculate_dphiref ||
377  this->calculate_phi_scaled || this->calculate_dphi_scaled ||
378  this->calculate_xyz || this->calculate_jxw ||
379  this->calculate_map_scaled || this->calculate_map ||
380  this->calculate_curl_phi || this->calculate_div_phi);
381 #endif
382 #endif // LIBMESH_ENABLE_DEPRECATED
383 
384  // set further terms necessary to do the requested task
385  if (calculate_jxw)
386  this->calculate_map = true;
387  if (this->calculate_dphi)
388  this->calculate_map = true;
389  if (this->calculate_dphi_scaled)
390  this->calculate_map_scaled = true;
391  if (calculate_xyz && !calculate_map)
392  // if Cartesian positions were requested but the calculation of map
393  // was not triggered, we'll opt for the 'scaled' variant.
394  this->calculate_map_scaled = true;
395  base_fe->calculate_phi = this->calculate_phi || this->calculate_phi_scaled
396  || this->calculate_dphi || this->calculate_dphi_scaled;
397  base_fe->calculate_dphi = this->calculate_dphi || this->calculate_dphi_scaled;
398  if (this->calculate_map || this->calculate_map_scaled
399  || this->calculate_dphiref)
400  {
401  base_fe->calculate_dphiref = true;
402  base_fe->get_xyz(); // trigger base_fe->fe_map to 'calculate_xyz'
403  base_fe->get_JxW(); // trigger base_fe->fe_map to 'calculate_dxyz'
404  }
405  base_fe->determine_calculations();
406 
407 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
408  if (this->calculate_d2phi)
409  libmesh_not_implemented();
410 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
411 }
412 
413 
414 
415 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
416 void
418 init_radial_shape_functions(const Elem * libmesh_dbg_var(inf_elem),
419  const std::vector<Point> * radial_pts)
420 {
421  libmesh_assert(radial_qrule.get() || radial_pts);
422  libmesh_assert(inf_elem);
423 
424  // Start logging the radial shape function initialization
425  LOG_SCOPE("init_radial_shape_functions()", "InfFE");
426 
427  // initialize most of the things related to physical approximation
428  const Order radial_approx_order = fe_type.radial_order;
429  const unsigned int n_radial_approx_shape_functions =
430  InfFERadial::n_dofs(radial_approx_order);
431 
432  const std::size_t n_radial_qp =
433  radial_pts ? radial_pts->size() : radial_qrule->n_points();
434  const std::vector<Point> & radial_qp =
435  radial_pts ? *radial_pts : radial_qrule->get_points();
436 
437  // the radial polynomials (eval)
438  if (calculate_phi || calculate_dphi || calculate_phi_scaled || calculate_dphi_scaled)
439  {
440  mode.resize (n_radial_approx_shape_functions);
441  for (unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
442  mode[i].resize (n_radial_qp);
443 
444  // evaluate the mode shapes in radial direction at radial quadrature points
445  for (unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
446  for (std::size_t p=0; p<n_radial_qp; ++p)
447  mode[i][p] = InfFE<Dim,T_radial,T_map>::eval (radial_qp[p](0), radial_approx_order, i);
448  }
449 
450  if (calculate_dphi || calculate_dphi_scaled)
451  {
452  dmodedv.resize (n_radial_approx_shape_functions);
453  for (unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
454  dmodedv[i].resize (n_radial_qp);
455 
456  // evaluate the mode shapes in radial direction at radial quadrature points
457  for (unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
458  for (std::size_t p=0; p<n_radial_qp; ++p)
459  dmodedv[i][p] = InfFE<Dim,T_radial,T_map>::eval_deriv (radial_qp[p](0), radial_approx_order, i);
460  }
461 
462  // the (1-v)/2 weight.
463  if (calculate_phi || calculate_phi_scaled || calculate_dphi || calculate_dphi_scaled)
464  {
465  som.resize (n_radial_qp);
466  // compute scalar values at radial quadrature points
467  for (std::size_t p=0; p<n_radial_qp; ++p)
468  som[p] = InfFERadial::decay (Dim, radial_qp[p](0));
469  }
470  if (calculate_dphi || calculate_dphi_scaled)
471  {
472  dsomdv.resize (n_radial_qp);
473  // compute scalar values at radial quadrature points
474  for (std::size_t p=0; p<n_radial_qp; ++p)
475  dsomdv[p] = InfFERadial::decay_deriv (Dim, radial_qp[p](0));
476  }
477 }
478 
479 
480 
481 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
482 void InfFE<Dim,T_radial,T_map>::init_shape_functions(const std::vector<Point> & radial_qp,
483  const std::vector<Point> & base_qp,
484  const Elem * inf_elem)
485 {
486  libmesh_assert(inf_elem);
487 
488  // Start logging the radial shape function initialization
489  LOG_SCOPE("init_shape_functions()", "InfFE");
490 
491  // fast access to some const ints for the radial data
492  const unsigned int n_radial_approx_sf = InfFERadial::n_dofs(fe_type.radial_order);
493  const std::size_t n_radial_qp = radial_qp.size();
494 #ifdef DEBUG
495  if (calculate_phi || calculate_dphi || calculate_phi_scaled || calculate_dphi_scaled)
496  libmesh_assert_equal_to(n_radial_approx_sf, mode.size());
497  if (calculate_phi || calculate_dphi || calculate_dphi_scaled)
498  libmesh_assert_equal_to(som.size(), n_radial_qp);
499 #endif
500 
501 
502  // initialize most of the quantities related to mapping
503 
504  // The element type and order to use in the base map
505  //const Order base_mapping_order = base_elem->default_order();
506 
507  // the number of base shape functions used to construct the map
508  // (Lagrange shape functions are used for mapping in the base)
509  //unsigned int n_base_mapping_shape_functions =
510  // InfFEBase::n_base_mapping_sf(*base_elem,
511  // base_mapping_order);
512 
513  // initialize most of the things related to physical approximation
514  unsigned int n_base_approx_shape_functions;
515  if (Dim > 1)
516  n_base_approx_shape_functions =
517  FEInterface::n_dofs(base_fe->get_fe_type(), base_elem.get());
518  else
519  n_base_approx_shape_functions = 1;
520 
521 
522  // update class member field
523  _n_total_approx_sf =
524  n_radial_approx_sf * n_base_approx_shape_functions;
525 
526 
527  // The number of the base quadrature points.
528  const unsigned int n_base_qp = cast_int<unsigned int>(base_qp.size());
529 
530  // The total number of quadrature points.
531  _n_total_qp = n_radial_qp * n_base_qp;
532 
533 
534  // initialize the node and shape numbering maps
535  {
536  // similar for the shapes: the i-th entry stores
537  // the associated base/radial shape number
538  _radial_shape_index.resize(_n_total_approx_sf);
539  _base_shape_index.resize(_n_total_approx_sf);
540 
541  // fill the shape index map
542  for (unsigned int n=0; n<_n_total_approx_sf; ++n)
543  {
544  compute_shape_indices (this->fe_type,
545  inf_elem,
546  n,
547  _base_shape_index[n],
548  _radial_shape_index[n]);
549  libmesh_assert_less (_base_shape_index[n], n_base_approx_shape_functions);
550  libmesh_assert_less (_radial_shape_index[n], n_radial_approx_sf);
551  }
552  }
553 
554  // resize the base data fields
555  //dist.resize(n_base_mapping_shape_functions);
556 
557  // resize the total data fields
558 
559  // the phase term varies with xi, eta and zeta(v): store it for _all_ qp
560  //
561  // when computing the phase, we need the base approximations
562  // therefore, initialize the phase here, but evaluate it
563  // in compute_shape_functions().
564  //
565  // the weight, though, is only needed at the radial quadrature points, n_radial_qp.
566  // but for a uniform interface to the protected data fields
567  // the weight data field (which are accessible from the outside) are expanded to _n_total_qp.
568  if (calculate_phi || calculate_dphi)
569  weight.resize (_n_total_qp);
570  if (calculate_phi_scaled || calculate_dphi_scaled)
571  weightxr_sq.resize (_n_total_qp);
572  if (calculate_dphi || calculate_dphi_scaled)
573  dweightdv.resize (n_radial_qp);
574  if (calculate_dphi)
575  dweight.resize (_n_total_qp);
576  if (calculate_dphi_scaled)
577  dweightxr_sq.resize(_n_total_qp);
578 
579  if (calculate_dphi || calculate_dphi_scaled)
580  dphase.resize (_n_total_qp);
581 
582  // this vector contains the integration weights for the combined quadrature rule
583  // if no quadrature rules are given, use only ones.
584  _total_qrule_weights.resize(_n_total_qp,1);
585 
586  // InfFE's data fields phi, dphi, dphidx, phi_map etc hold the _total_
587  // shape and mapping functions, respectively
588  {
589  if (calculate_map_scaled)
590  JxWxdecay.resize(_n_total_qp);
591  if (calculate_jxw)
592  JxW.resize(_n_total_qp);
593  if (calculate_map_scaled || calculate_map)
594  {
595  xyz.resize(_n_total_qp);
596  dxidx_map_scaled.resize(_n_total_qp);
597  dxidy_map_scaled.resize(_n_total_qp);
598  dxidz_map_scaled.resize(_n_total_qp);
599  detadx_map_scaled.resize(_n_total_qp);
600  detady_map_scaled.resize(_n_total_qp);
601  detadz_map_scaled.resize(_n_total_qp);
602  dzetadx_map_scaled.resize(_n_total_qp);
603  dzetady_map_scaled.resize(_n_total_qp);
604  dzetadz_map_scaled.resize(_n_total_qp);
605  }
606  if (calculate_map)
607  {
608  dxidx_map.resize(_n_total_qp);
609  dxidy_map.resize(_n_total_qp);
610  dxidz_map.resize(_n_total_qp);
611  detadx_map.resize(_n_total_qp);
612  detady_map.resize(_n_total_qp);
613  detadz_map.resize(_n_total_qp);
614  dzetadx_map.resize(_n_total_qp);
615  dzetady_map.resize(_n_total_qp);
616  dzetadz_map.resize(_n_total_qp);
617  }
618  if (calculate_phi)
619  phi.resize (_n_total_approx_sf);
620  if (calculate_phi_scaled)
621  phixr.resize (_n_total_approx_sf);
622  if (calculate_dphi)
623  {
624  dphi.resize (_n_total_approx_sf);
625  dphidx.resize (_n_total_approx_sf);
626  dphidy.resize (_n_total_approx_sf);
627  dphidz.resize (_n_total_approx_sf);
628  }
629 
630  if (calculate_dphi_scaled)
631  {
632  dphixr.resize (_n_total_approx_sf);
633  dphixr_sq.resize(_n_total_approx_sf);
634  }
635 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
636 
637  if (calculate_d2phi)
638  {
639  libmesh_not_implemented();
640  d2phi.resize (_n_total_approx_sf);
641  d2phidx2.resize (_n_total_approx_sf);
642  d2phidxdy.resize (_n_total_approx_sf);
643  d2phidxdz.resize (_n_total_approx_sf);
644  d2phidy2.resize (_n_total_approx_sf);
645  d2phidydz.resize (_n_total_approx_sf);
646  d2phidz2.resize (_n_total_approx_sf);
647  d2phidxi2.resize (_n_total_approx_sf);
648 
649  if (Dim > 1)
650  {
651  d2phidxideta.resize(_n_total_approx_sf);
652  d2phideta2.resize(_n_total_approx_sf);
653  }
654 
655  if (Dim > 2)
656  {
657  d2phidetadzeta.resize(_n_total_approx_sf);
658  d2phidxidzeta.resize(_n_total_approx_sf);
659  d2phidzeta2.resize(_n_total_approx_sf);
660  }
661  }
662 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
663 
664  if (calculate_dphi || calculate_dphi_scaled)
665  {
666  dphidxi.resize (_n_total_approx_sf);
667 
668  if (Dim > 1)
669  dphideta.resize(_n_total_approx_sf);
670 
671  if (Dim == 3)
672  dphidzeta.resize(_n_total_approx_sf);
673  }
674 
675  }
676 
677  // collect all the for loops, where inner vectors are
678  // resized to the appropriate number of quadrature points
679  {
680  if (calculate_phi)
681  for (unsigned int i=0; i<_n_total_approx_sf; ++i)
682  phi[i].resize (_n_total_qp);
683 
684  if (calculate_dphi)
685  for (unsigned int i=0; i<_n_total_approx_sf; ++i)
686  {
687  dphi[i].resize (_n_total_qp);
688  dphidx[i].resize (_n_total_qp);
689  dphidy[i].resize (_n_total_qp);
690  dphidz[i].resize (_n_total_qp);
691  }
692 
693  if (calculate_phi_scaled)
694  for (unsigned int i=0; i<_n_total_approx_sf; ++i)
695  {
696  phixr[i].resize (_n_total_qp);
697  }
698  if (calculate_dphi_scaled)
699  for (unsigned int i=0; i<_n_total_approx_sf; ++i)
700  {
701  dphixr[i].resize(_n_total_qp);
702  dphixr_sq[i].resize(_n_total_qp);
703  }
704 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
705  if (calculate_d2phi)
706  for (unsigned int i=0; i<_n_total_approx_sf; ++i)
707  {
708  d2phi[i].resize (_n_total_qp);
709  d2phidx2[i].resize (_n_total_qp);
710  d2phidxdy[i].resize (_n_total_qp);
711  d2phidxdz[i].resize (_n_total_qp);
712  d2phidy2[i].resize (_n_total_qp);
713  d2phidydz[i].resize (_n_total_qp);
714  d2phidy2[i].resize (_n_total_qp);
715  d2phidxi2[i].resize (_n_total_qp);
716 
717  if (Dim > 1)
718  {
719  d2phidxideta[i].resize (_n_total_qp);
720  d2phideta2[i].resize (_n_total_qp);
721  }
722  if (Dim > 2)
723  {
724  d2phidxidzeta[i].resize (_n_total_qp);
725  d2phidetadzeta[i].resize (_n_total_qp);
726  d2phidzeta2[i].resize (_n_total_qp);
727  }
728  }
729 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
730 
731  if (calculate_dphi || calculate_dphi_scaled)
732  for (unsigned int i=0; i<_n_total_approx_sf; ++i)
733  {
734  dphidxi[i].resize (_n_total_qp);
735 
736  if (Dim > 1)
737  dphideta[i].resize (_n_total_qp);
738 
739  if (Dim == 3)
740  dphidzeta[i].resize (_n_total_qp);
741 
742  }
743 
744  }
745  {
746  // (a) compute scalar values at _all_ quadrature points -- for uniform
747  // access from the outside to these fields
748  // (b) form a std::vector<Real> which contains the appropriate weights
749  // of the combined quadrature rule!
750  libmesh_assert_equal_to (radial_qp.size(), n_radial_qp);
751 
752  if (radial_qrule && base_qrule)
753  {
754  const std::vector<Real> & radial_qw = radial_qrule->get_weights();
755  const std::vector<Real> & base_qw = base_qrule->get_weights();
756  libmesh_assert_equal_to (radial_qw.size(), n_radial_qp);
757  libmesh_assert_equal_to (base_qw.size(), n_base_qp);
758 
759  for (unsigned int rp=0; rp<n_radial_qp; ++rp)
760  for (unsigned int bp=0; bp<n_base_qp; ++bp)
761  _total_qrule_weights[bp + rp*n_base_qp] = radial_qw[rp] * base_qw[bp];
762  }
763 
764 
765  for (unsigned int rp=0; rp<n_radial_qp; ++rp)
766  {
767  if (calculate_phi || calculate_dphi)
768  for (unsigned int bp=0; bp<n_base_qp; ++bp)
769  weight[bp + rp*n_base_qp] = InfFERadial::D(radial_qp[rp](0));
770 
771  if ( calculate_phi_scaled)
772  for (unsigned int bp=0; bp<n_base_qp; ++bp)
773  weightxr_sq[bp + rp*n_base_qp] = InfFERadial::Dxr_sq(radial_qp[rp](0));
774 
775  if ( calculate_dphi || calculate_dphi_scaled)
776  dweightdv[rp] = InfFERadial::D_deriv(radial_qp[rp](0));
777  }
778  }
779 }
780 
781 
782 
783 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
785  const std::vector<Point> & base_qp,
786  const std::vector<Point> & radial_qp
787  )
788 {
789  libmesh_assert(inf_elem);
790  // at least check whether the base element type is correct.
791  // otherwise this version of computing dist would give problems
792  libmesh_assert_equal_to (base_elem->type(),
793  InfFEBase::get_elem_type(inf_elem->type()));
794 
795  // Start logging the overall computation of shape functions
796  LOG_SCOPE("compute_shape_functions()", "InfFE");
797 
798  //const unsigned int n_radial_qp = cast_int<unsigned int>(som.size());
799  //const unsigned int n_base_qp = cast_int<unsigned int>(S_map[0].size());
800  const std::size_t n_radial_qp = radial_qp.size();
801  const unsigned int n_base_qp = base_qp.size();
802 
803  libmesh_assert_equal_to (_n_total_qp, n_radial_qp*n_base_qp);
804  libmesh_assert_equal_to (_n_total_qp, _total_qrule_weights.size());
805 #ifdef DEBUG
806  if (som.size() > 0)
807  libmesh_assert_equal_to(n_radial_qp, som.size());
808 
809  if (this->calculate_map || this->calculate_map_scaled)
810  {
811  // these vectors are needed later; initialize here already to have access to
812  // n_base_qp etc.
813  const std::vector<std::vector<Real>> & S_map = (base_fe->get_fe_map()).get_phi_map();
814  if (S_map[0].size() > 0)
815  libmesh_assert_equal_to(n_base_qp, S_map[0].size());
816  }
817  if (radial_qrule)
818  libmesh_assert_equal_to(n_radial_qp, radial_qrule->n_points());
819  if (base_qrule)
820  libmesh_assert_equal_to(n_base_qp, base_qrule->n_points());
821  libmesh_assert_equal_to(_n_total_qp % n_radial_qp, 0); // "Error in the structure of quadrature points!");
822 #endif
823 
824 
825  _n_total_approx_sf = InfFERadial::n_dofs(fe_type.radial_order) *
826  FEInterface::n_dofs(base_fe->get_fe_type(), base_elem.get());
827 
828 
829 
830  const Point origin = inf_elem->origin();
831 
832  // Compute the shape function values (and derivatives)
833  // at the Quadrature points. Note that the actual values
834  // have already been computed via init_shape_functions
835 
836  unsigned int elem_dim = inf_elem->dim();
837  // Compute the value of the derivative shape function i at quadrature point p
838  switch (elem_dim)
839  {
840  case 1:
841  case 2:
842  {
843  libmesh_not_implemented();
844  break;
845  }
846  case 3:
847  {
848  std::vector<std::vector<Real>> S (0);
849  std::vector<std::vector<Real>> Ss(0);
850  std::vector<std::vector<Real>> St(0);
851 
852  std::vector<Real> base_dxidx (0);
853  std::vector<Real> base_dxidy (0);
854  std::vector<Real> base_dxidz (0);
855  std::vector<Real> base_detadx(0);
856  std::vector<Real> base_detady(0);
857  std::vector<Real> base_detadz(0);
858 
859  std::vector<Point> base_xyz (0);
860 
861  if (calculate_phi || calculate_phi_scaled ||
862  calculate_dphi || calculate_dphi_scaled)
863  S=base_fe->phi;
864 
865  // fast access to the approximation and mapping shapes of base_fe
866  if (calculate_map || calculate_map_scaled)
867  {
868  Ss = base_fe->dphidxi;
869  St = base_fe->dphideta;
870 
871  base_dxidx = base_fe->get_dxidx();
872  base_dxidy = base_fe->get_dxidy();
873  base_dxidz = base_fe->get_dxidz();
874  base_detadx = base_fe->get_detadx();
875  base_detady = base_fe->get_detady();
876  base_detadz = base_fe->get_detadz();
877 
878  base_xyz = base_fe->get_xyz();
879  }
880 
881  ElemType base_type= base_elem->type();
882 
883 #ifdef DEBUG
884  if (calculate_phi)
885  libmesh_assert_equal_to (phi.size(), _n_total_approx_sf);
886  if (calculate_dphi)
887  {
888  libmesh_assert_equal_to (dphidxi.size(), _n_total_approx_sf);
889  libmesh_assert_equal_to (dphideta.size(), _n_total_approx_sf);
890  libmesh_assert_equal_to (dphidzeta.size(), _n_total_approx_sf);
891  }
892 #endif
893 
894  unsigned int tp=0; // total qp
895  for (unsigned int rp=0; rp<n_radial_qp; ++rp) // over radial qps
896  for (unsigned int bp=0; bp<n_base_qp; ++bp) // over base qps
897 
898  { // First compute the map from base element quantities to physical space:
899 
900  // initialize them with invalid value to not use them
901  // without setting them to the correct value before.
902  Point unit_r(NAN);
903  RealGradient grad_a_scaled(NAN);
904  Real a(NAN);
905  Real r_norm(NAN);
906  if (calculate_map || calculate_map_scaled)
907  {
908  xyz[tp] = InfFEMap::map(elem_dim, inf_elem, Point(base_qp[bp](0),base_qp[bp](1),radial_qp[rp](0)));
909 
910  const Point r(xyz[tp]-origin);
911  a=(base_xyz[bp]-origin).norm();
912  r_norm = r.norm();
913 
914  // check that 'som' == a/r.
915 #ifndef NDEVEL
916  if (som.size())
917  libmesh_assert_less(std::abs(som[rp] -a/r_norm) , 1e-7);
918 #endif
919  unit_r=(r/r_norm);
920 
921  // They are used for computing the normal and do not correspond to the direction of eta and xi in this element:
922  // Due to the stretch of these axes in radial direction, they are deformed.
923  Point e_xi(base_dxidx[bp],
924  base_dxidy[bp],
925  base_dxidz[bp]);
926  Point e_eta(base_detadx[bp],
927  base_detady[bp],
928  base_detadz[bp]);
929 
930  const RealGradient normal=e_eta.cross(e_xi).unit();
931 
932  // grad a = a/r.norm() * grad_a_scaled
933  grad_a_scaled=unit_r - normal/(normal*unit_r);
934 
935  const Real dxi_er=base_dxidx[bp]* unit_r(0) + base_dxidy[bp] *unit_r(1) + base_dxidz[bp] *unit_r(2);
936  const Real deta_er=base_detadx[bp]*unit_r(0) + base_detady[bp]*unit_r(1) + base_detadz[bp]*unit_r(2);
937 
938  // in case of non-affine map, further terms need to be taken into account,
939  // involving \p e_eta and \p e_xi and thus recursive computation is needed
940  if (!base_elem->has_affine_map())
941  {
950  const unsigned int n_sf = base_elem->n_nodes();
951  RealGradient tmp(0.,0.,0.);
952  for (unsigned int i=0; i< n_sf; ++i)
953  {
954  RealGradient dL_da_i = (FE<2,LAGRANGE>::shape_deriv(base_type,
955  base_elem->default_order(),
956  i, 0, base_qp[bp]) * e_xi
957  +FE<2,LAGRANGE>::shape_deriv(base_type,
958  base_elem->default_order(),
959  i, 1, base_qp[bp]) * e_eta);
960 
961  tmp += (base_elem->node_ref(i) -origin).norm()* dL_da_i;
962 
963  }
964  libmesh_assert(tmp*unit_r < .95 ); // in a proper setup, tmp should have only a small radial component.
965  grad_a_scaled = ( tmp - (tmp*unit_r)*unit_r ) / ( 1. - tmp*unit_r);
966 
967  }
968 
969  // 'scale' = r/a
970  dxidx_map_scaled[tp] = (grad_a_scaled(0) - unit_r(0))*dxi_er +base_dxidx[bp];
971  dxidy_map_scaled[tp] = (grad_a_scaled(1) - unit_r(1))*dxi_er +base_dxidy[bp];
972  dxidz_map_scaled[tp] = (grad_a_scaled(2) - unit_r(2))*dxi_er +base_dxidz[bp];
973 
974  // 'scale' = r/a
975  detadx_map_scaled[tp] = (grad_a_scaled(0) - unit_r(0))*deta_er + base_detadx[bp];
976  detady_map_scaled[tp] = (grad_a_scaled(1) - unit_r(1))*deta_er + base_detady[bp];
977  detadz_map_scaled[tp] = (grad_a_scaled(2) - unit_r(2))*deta_er + base_detadz[bp];
978 
979  // 'scale' = (r/a)**2
980  dzetadx_map_scaled[tp] =-2./a*(grad_a_scaled(0) - unit_r(0));
981  dzetady_map_scaled[tp] =-2./a*(grad_a_scaled(1) - unit_r(1));
982  dzetadz_map_scaled[tp] =-2./a*(grad_a_scaled(2) - unit_r(2));
983 
984  }
985 
986  if (calculate_map)
987  {
988  dxidx_map[tp] = a/r_norm * dxidx_map_scaled[tp];
989  dxidy_map[tp] = a/r_norm * dxidy_map_scaled[tp];
990  dxidz_map[tp] = a/r_norm * dxidz_map_scaled[tp];
991 
992  detadx_map[tp] = a/r_norm * detadx_map_scaled[tp];
993  detady_map[tp] = a/r_norm * detady_map_scaled[tp];
994  detadz_map[tp] = a/r_norm * detadz_map_scaled[tp];
995 
996  // dzetadx = dzetadr*dr/dx - 2/r * grad_a
997  // = dzetadr*dr/dx - 2*a/r^2 * grad_a_scaled
998  dzetadx_map[tp] =-2.*a/(r_norm*r_norm)*(grad_a_scaled(0) - unit_r(0));
999  dzetady_map[tp] =-2.*a/(r_norm*r_norm)*(grad_a_scaled(1) - unit_r(1));
1000  dzetadz_map[tp] =-2.*a/(r_norm*r_norm)*(grad_a_scaled(2) - unit_r(2));
1001 
1002  if (calculate_jxw)
1003  {
1004  Real inv_jac = (dxidx_map[tp]*( detady_map[tp]*dzetadz_map[tp]- dzetady_map[tp]*detadz_map[tp]) +
1005  detadx_map[tp]*(dzetady_map[tp]* dxidz_map[tp]- dxidy_map[tp]*dzetadz_map[tp]) +
1006  dzetadx_map[tp]*( dxidy_map[tp]*detadz_map[tp]- detady_map[tp]* dxidz_map[tp]));
1007 
1008  if (inv_jac <= 1e-10)
1009  {
1010  libmesh_error_msg("ERROR: negative inverse Jacobian " \
1011  << inv_jac \
1012  << " at point " \
1013  << xyz[tp] \
1014  << " in element " \
1015  << inf_elem->id());
1016  }
1017 
1018 
1019  JxW[tp] = _total_qrule_weights[tp]/inv_jac;
1020  }
1021 
1022  }
1023  if (calculate_map_scaled)
1024  {
1025  Real inv_jacxR_pow4 = (dxidx_map_scaled[tp] *( detady_map_scaled[tp]*dzetadz_map_scaled[tp]
1026  - dzetady_map_scaled[tp]* detadz_map_scaled[tp]) +
1027  detadx_map_scaled[tp] *( dzetady_map_scaled[tp]* dxidz_map_scaled[tp]
1028  - dxidy_map_scaled[tp]*dzetadz_map_scaled[tp]) +
1029  dzetadx_map_scaled[tp]*( dxidy_map_scaled[tp]* detadz_map_scaled[tp]
1030  -detady_map_scaled[tp]* dxidz_map_scaled[tp]));
1031  if (inv_jacxR_pow4 <= 1e-7)
1032  {
1033  libmesh_error_msg("ERROR: negative weighted inverse Jacobian " \
1034  << inv_jacxR_pow4 \
1035  << " at point " \
1036  << xyz[tp] \
1037  << " in element " \
1038  << inf_elem->id());
1039  }
1040 
1041  JxWxdecay[tp] = _total_qrule_weights[tp]/inv_jacxR_pow4;
1042  }
1043 
1044  // phase term mu(r)=i*k*(r-a).
1045  // skip i*k: it is added separately during matrix assembly.
1046 
1047  if (calculate_dphi || calculate_dphi_scaled)
1048  dphase[tp] = unit_r - grad_a_scaled*a/r_norm;
1049 
1050  if (calculate_dphi)
1051  {
1052  dweight[tp](0) = dweightdv[rp] * dzetadx_map[tp];
1053  dweight[tp](1) = dweightdv[rp] * dzetady_map[tp];
1054  dweight[tp](2) = dweightdv[rp] * dzetadz_map[tp];
1055  }
1056  if (calculate_dphi_scaled)
1057  {
1058  dweightxr_sq[tp](0) = dweightdv[rp] * dzetadx_map_scaled[tp];
1059  dweightxr_sq[tp](1) = dweightdv[rp] * dzetady_map_scaled[tp];
1060  dweightxr_sq[tp](2) = dweightdv[rp] * dzetadz_map_scaled[tp];
1061  }
1062 
1063  if (calculate_phi || calculate_phi_scaled || calculate_dphi || calculate_dphi_scaled)
1064  // compute the shape-functions and derivative quantities:
1065  for (unsigned int i=0; i <_n_total_approx_sf ; ++i)
1066  {
1067  // let the index vectors take care of selecting the appropriate base/radial shape
1068  unsigned int bi = _base_shape_index [i];
1069  unsigned int ri = _radial_shape_index[i];
1070  if (calculate_phi)
1071  phi [i][tp] = S [bi][bp] * mode[ri][rp] * som[rp];
1072 
1073  if (calculate_phi_scaled)
1074  phixr [i][tp] = S [bi][bp] * mode[ri][rp];
1075 
1076  if (calculate_dphi || calculate_dphi_scaled)
1077  {
1078  dphidxi [i][tp] = Ss[bi][bp] * mode[ri][rp] * som[rp];
1079  dphideta [i][tp] = St[bi][bp] * mode[ri][rp] * som[rp];
1080  dphidzeta[i][tp] = S [bi][bp]
1081  * (dmodedv[ri][rp] * som[rp] + mode[ri][rp] * dsomdv[rp]);
1082  }
1083 
1084  if (calculate_dphi)
1085  {
1086 
1087  // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
1088  dphi[i][tp](0) =
1089  dphidx[i][tp] = (dphidxi[i][tp]*dxidx_map[tp] +
1090  dphideta[i][tp]*detadx_map[tp] +
1091  dphidzeta[i][tp]*dzetadx_map[tp]);
1092 
1093  // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
1094  dphi[i][tp](1) =
1095  dphidy[i][tp] = (dphidxi[i][tp]*dxidy_map[tp] +
1096  dphideta[i][tp]*detady_map[tp] +
1097  dphidzeta[i][tp]*dzetady_map[tp]);
1098 
1099  // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
1100  dphi[i][tp](2) =
1101  dphidz[i][tp] = (dphidxi[i][tp]*dxidz_map[tp] +
1102  dphideta[i][tp]*detadz_map[tp] +
1103  dphidzeta[i][tp]*dzetadz_map[tp]);
1104 
1105  }
1106  if (calculate_dphi_scaled)
1107  { // we don't distinguish between the different levels of scaling here...
1108 
1109  dphixr[i][tp](0)= (dphidxi[i][tp]*dxidx_map_scaled[tp] +
1110  dphideta[i][tp]*detadx_map_scaled[tp] +
1111  dphidzeta[i][tp]*dzetadx_map_scaled[tp]*som[rp]);
1112 
1113  dphixr[i][tp](1) = (dphidxi[i][tp]*dxidy_map_scaled[tp] +
1114  dphideta[i][tp]*detady_map_scaled[tp] +
1115  dphidzeta[i][tp]*dzetady_map_scaled[tp]*som[rp]);
1116 
1117  dphixr[i][tp](2) = (dphidxi[i][tp]*dxidz_map_scaled[tp] +
1118  dphideta[i][tp]*detadz_map_scaled[tp] +
1119  dphidzeta[i][tp]*dzetadz_map_scaled[tp]*som[rp]);
1120 
1121  const Real dphidxixr = Ss[bi][bp] * mode[ri][rp];
1122  const Real dphidetaxr= St[bi][bp] * mode[ri][rp];
1123 
1124  dphixr_sq[i][tp](0)= (dphidxixr*dxidx_map_scaled[tp] +
1125  dphidetaxr*detadx_map_scaled[tp] +
1126  dphidzeta[i][tp]*dzetadx_map_scaled[tp]);
1127 
1128  dphixr_sq[i][tp](1) = (dphidxixr*dxidy_map_scaled[tp] +
1129  dphidetaxr*detady_map_scaled[tp] +
1130  dphidzeta[i][tp]*dzetady_map_scaled[tp]);
1131 
1132  dphixr_sq[i][tp](2) = (dphidxixr*dxidz_map_scaled[tp] +
1133  dphidetaxr*detadz_map_scaled[tp] +
1134  dphidzeta[i][tp]*dzetadz_map_scaled[tp]);
1135  }
1136 
1137  }
1138  tp++;
1139  }
1140 
1141  break;
1142  }
1143 
1144  default:
1145  libmesh_error_msg("Unsupported dim = " << dim);
1146  }
1147 }
1148 
1149 
1150 
1151 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
1153 {
1154  // We never call this.
1155  libmesh_not_implemented();
1156  return false;
1157 }
1158 
1159 } // namespace libMesh
1160 
1161 
1162 // Explicit instantiations
1163 #include "libmesh/inf_fe_instantiate_1D.h"
1164 #include "libmesh/inf_fe_instantiate_2D.h"
1165 #include "libmesh/inf_fe_instantiate_3D.h"
1166 
1167 
1168 
1169 #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
virtual void determine_calculations() override
Determine which values are to be calculated, for both the FE itself and for the FEMap.
Definition: inf_fe.C:329
ElemType
Defines an enum for geometric element types.
static ElemType get_elem_type(const ElemType type)
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
A specific instantiation of the FEBase class.
Definition: fe.h:42
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:355
static unsigned int n_dofs(const Order o_radial)
Definition: inf_fe.h:113
const Real radius
auto norm() const -> decltype(std::norm(T()))
Definition: type_vector.h:907
virtual Point origin() const
Definition: elem.h:1914
static Point map(const unsigned int dim, const Elem *inf_elem, const Point &reference_point)
Definition: inf_fe_map.C:40
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
Definition: int_range.h:53
static Real decay_deriv(const unsigned int dim, const Real)
Definition: inf_fe.h:1309
unsigned int dim
static Real Dxr_sq(const Real)
Definition: inf_fe.h:84
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:482
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
static FEFieldType field_type(const FEType &fe_type)
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The libMesh namespace provides an interface to certain functionality in the library.
static std::unique_ptr< const Elem > build_elem(const Elem *inf_elem)
Build the base element of an infinite element.
virtual QuadratureType type() const =0
static Real D(const Real v)
Definition: inf_fe.h:82
A specific instantiation of the FEBase class.
Definition: fe.h:127
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:1270
dof_id_type id() const
Definition: dof_object.h:828
std::unique_ptr< FEBase > base_fe
Have a FE<Dim-1,T_base> handy for base approximation.
Definition: inf_fe.h:1223
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:437
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libmesh_assert(ctx)
Order get_order() const
Definition: quadrature.h:249
static Real decay(const unsigned int dim, const Real v)
Definition: inf_fe.h:1283
unsigned int get_dim() const
Definition: quadrature.h:150
InfMapType inf_map
The coordinate mapping type of the infinite element.
Definition: fe_type.h:275
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
Definition: type_vector.h:884
auto norm(const T &a) -> decltype(std::abs(a))
Definition: tensor_tools.h:74
virtual bool shapes_need_reinit() const override
Definition: inf_fe.C:1152
void init_radial_shape_functions(const Elem *inf_elem, const std::vector< Point > *radial_pts=nullptr)
Some of the member data only depend on the radial part of the infinite element.
Definition: inf_fe.C:418
FEFamily radial_family
The type of approximation in radial direction.
Definition: fe_type.h:267
void compute_shape_functions(const Elem *inf_elem, const std::vector< Point > &base_qp, const std::vector< Point > &radial_qp)
After having updated the jacobian and the transformation from local to global coordinates in FEAbstra...
Definition: inf_fe.C:784
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned short dim() const =0
static Real D_deriv(const Real v)
Definition: inf_fe.h:90
static std::unique_ptr< QBase > build(std::string_view name, const unsigned int dim, const Order order=INVALID_ORDER)
Builds a specific quadrature rule based on the name string.
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:120
virtual void attach_quadrature_rule(QBase *q) override
The use of quadrature rules with the InfFE class is somewhat different from the approach of the FE cl...
Definition: inf_fe.C:79
FEType fe_type
The finite element type for this object.
Definition: fe_abstract.h:730
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
The QBase class provides the basic functionality from which various quadrature rules can be derived...
Definition: quadrature.h:61
This class forms the foundation from which generic finite elements may be derived.
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:109