libMesh
Loading...
Searching...
No Matches
inf_fe_static.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 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#include "libmesh/libmesh_config.h"
20
21#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
22
23// libmesh includes
24#include "libmesh/inf_fe.h"
25#include "libmesh/inf_fe_macro.h"
26#include "libmesh/fe.h"
27#include "libmesh/fe_interface.h"
28#include "libmesh/fe_compute_data.h"
29#include "libmesh/elem.h"
30#include "libmesh/enum_to_string.h"
31
32#include "libmesh/remote_elem.h"
33
34// to fetch NodeConstraintRow:
35#include "libmesh/dof_map.h"
36
37namespace libMesh
38{
39
40
41// ------------------------------------------------------------
42// InfFE class static member initialization
43template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
44ElemType InfFE<Dim,T_radial,T_map>::_compute_node_indices_fast_current_elem_type = INVALID_ELEM;
45
46#ifdef DEBUG
47
48template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
49bool InfFE<Dim,T_radial,T_map>::_warned_for_nodal_soln = false;
50
51
52template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
53bool InfFE<Dim,T_radial,T_map>::_warned_for_shape = false;
54
55template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
56bool InfFE<Dim,T_radial,T_map>::_warned_for_dshape = false;
57
58#endif
59
60
61
62
63// ------------------------------------------------------------
64// InfFE static class members
65template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
66unsigned int InfFE<Dim,T_radial,T_map>::n_dofs(const FEType & fet,
67 const Elem * inf_elem)
68{
69 // The "base" Elem is a non-infinite Elem corresponding to side 0 of
70 // the InfElem.
71 auto base_elem = inf_elem->build_side_ptr(0);
72
73 if (Dim > 1)
74 return FEInterface::n_dofs(fet, base_elem.get()) *
75 InfFERadial::n_dofs(fet.radial_order);
76 else
77 return InfFERadial::n_dofs(fet.radial_order);
78}
79
80
81
82#ifdef LIBMESH_ENABLE_DEPRECATED
83template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
84unsigned int InfFE<Dim,T_radial,T_map>::n_dofs_at_node (const FEType & fet,
85 const ElemType inf_elem_type,
86 const unsigned int n)
87{
88 libmesh_deprecated();
89
90 const ElemType base_et (InfFEBase::get_elem_type(inf_elem_type));
91
92 unsigned int n_base, n_radial;
93 compute_node_indices(inf_elem_type, n, n_base, n_radial);
94
95 // libMesh::out << "elem_type=" << inf_elem_type
96 // << ", fet.radial_order=" << fet.radial_order
97 // << ", n=" << n
98 // << ", n_radial=" << n_radial
99 // << ", n_base=" << n_base
100 // << std::endl;
101
102 if (Dim > 1)
103 return FEInterface::n_dofs_at_node(Dim-1, fet, base_et, n_base)
104 * InfFERadial::n_dofs_at_node(fet.radial_order, n_radial);
105 else
106 return InfFERadial::n_dofs_at_node(fet.radial_order, n_radial);
107}
108#endif // LIBMESH_ENABLE_DEPRECATED
109
110
111
112template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
113unsigned int InfFE<Dim,T_radial,T_map>::n_dofs_at_node (const FEType & fet,
114 const Elem * inf_elem,
115 const unsigned int n)
116{
117 // The "base" Elem is a non-infinite Elem corresponding to side 0 of
118 // the InfElem.
119 auto base_elem = inf_elem->build_side_ptr(0);
120
121 unsigned int n_base, n_radial;
122 compute_node_indices(inf_elem->type(), n, n_base, n_radial);
123
124 if (Dim > 1)
125 return FEInterface::n_dofs_at_node(fet, base_elem.get(), n_base)
126 * InfFERadial::n_dofs_at_node(fet.radial_order, n_radial);
127 else
128 return InfFERadial::n_dofs_at_node(fet.radial_order, n_radial);
129}
130
131
132
133#ifdef LIBMESH_ENABLE_DEPRECATED
134template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
135unsigned int InfFE<Dim,T_radial,T_map>::n_dofs_per_elem (const FEType & fet,
136 const ElemType inf_elem_type)
137{
138 libmesh_deprecated();
139
140 const ElemType base_et (InfFEBase::get_elem_type(inf_elem_type));
141
142 if (Dim > 1)
143 return FEInterface::n_dofs_per_elem(Dim-1, fet, base_et)
144 * InfFERadial::n_dofs_per_elem(fet.radial_order);
145 else
146 return InfFERadial::n_dofs_per_elem(fet.radial_order);
147}
148#endif // LIBMESH_ENABLE_DEPRECATED
149
150
151
152template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
153unsigned int InfFE<Dim,T_radial,T_map>::n_dofs_per_elem (const FEType & fet,
154 const Elem * inf_elem)
155{
156 // The "base" Elem is a non-infinite Elem corresponding to side 0 of
157 // the InfElem.
158 auto base_elem = inf_elem->build_side_ptr(0);
159
160 if (Dim > 1)
161 return FEInterface::n_dofs_per_elem(fet, base_elem.get())
162 * InfFERadial::n_dofs_per_elem(fet.radial_order);
163 else
164 return InfFERadial::n_dofs_per_elem(fet.radial_order);
165}
166
167
168
169template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
170void InfFE<Dim,T_radial,T_map>::nodal_soln(const FEType & /* fet */,
171 const Elem * /* elem */,
172 const std::vector<Number> & /* elem_soln */,
173 std::vector<Number> & nodal_soln)
174{
175#ifdef DEBUG
176 if (!_warned_for_nodal_soln)
177 {
178 libMesh::err << "WARNING: nodal_soln(...) does _not_ work for infinite elements." << std::endl
179 << " Will return an empty nodal solution. Use " << std::endl
180 << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!" << std::endl;
181 _warned_for_nodal_soln = true;
182 }
183#endif
184
185 /*
186 * In the base the infinite element couples to
187 * conventional finite elements. To not destroy
188 * the values there, clear \p nodal_soln. This
189 * indicates to the user of \p nodal_soln to
190 * not use this result.
191 */
192 nodal_soln.clear();
193 libmesh_assert (nodal_soln.empty());
194 return;
195}
196
197
198
199
200
201
202
203
204#ifdef LIBMESH_ENABLE_DEPRECATED
205template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
206Real InfFE<Dim,T_radial,T_map>::shape(const FEType & fet,
207 const ElemType inf_elem_type,
208 const unsigned int i,
209 const Point & p)
210{
211 libmesh_deprecated();
212
213 libmesh_assert_not_equal_to (Dim, 0);
214
215#ifdef DEBUG
216 // this makes only sense when used for mapping
217 if ((T_radial != INFINITE_MAP) && !_warned_for_shape)
218 {
219 libMesh::err << "WARNING: InfFE<Dim,T_radial,T_map>::shape(...) does _not_" << std::endl
220 << " return the correct trial function! Use " << std::endl
221 << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!"
222 << std::endl;
223 _warned_for_shape = true;
224 }
225#endif
226
227 const ElemType base_et (InfFEBase::get_elem_type(inf_elem_type));
228 const Order o_radial (fet.radial_order);
229 const Real v (p(Dim-1));
230
231 unsigned int i_base, i_radial;
232 compute_shape_indices(fet, inf_elem_type, i, i_base, i_radial);
233
234 //TODO:[SP/DD] exp(ikr) is still missing here!
235 // but is it intended? It would be probably somehow nice, but than it would be Number, not Real !
236 // --> thus it would destroy the interface...
237 if (Dim > 1)
238 return FEInterface::shape(Dim-1, fet, base_et, i_base, p)
239 * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
240 * InfFERadial::decay(Dim,v);
241 else
242 return InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
243 * InfFERadial::decay(Dim,v);
244}
245#endif // LIBMESH_ENABLE_DEPRECATED
246
247
248
249template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
250Real InfFE<Dim,T_radial,T_map>::shape(const FEType & fet,
251 const Elem * inf_elem,
252 const unsigned int i,
253 const Point & p)
254{
255 libmesh_assert(inf_elem);
256 libmesh_assert_not_equal_to (Dim, 0);
257
258#ifdef DEBUG
259 // this makes only sense when used for mapping
260 if ((T_radial != INFINITE_MAP) && !_warned_for_shape)
261 {
262 libMesh::err << "WARNING: InfFE<Dim,T_radial,T_map>::shape(...) does _not_" << std::endl
263 << " return the correct trial function! Use " << std::endl
264 << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!"
265 << std::endl;
266 _warned_for_shape = true;
267 }
268#endif
269
270 const Order o_radial (fet.radial_order);
271 const Real v (p(Dim-1));
272 std::unique_ptr<const Elem> base_el (inf_elem->build_side_ptr(0));
273
274 unsigned int i_base, i_radial;
275 compute_shape_indices(fet, inf_elem, i, i_base, i_radial);
276
277 if (Dim > 1)
278 return FEInterface::shape(fet, base_el.get(), i_base, p)
279 * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
280 * InfFERadial::decay(Dim,v);
281 else
282 return InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
283 * InfFERadial::decay(Dim,v);
284}
285
286
287
288template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
289Real InfFE<Dim,T_radial,T_map>::shape(const FEType fet,
290 const Elem * inf_elem,
291 const unsigned int i,
292 const Point & p,
293 const bool add_p_level)
294{
295 if (add_p_level)
296 {
297 FEType tmp_fet=fet;
298 tmp_fet = fet.order + add_p_level * inf_elem->p_level();
299 return InfFE<Dim,T_radial,T_map>::shape(tmp_fet, inf_elem, i, p);
300 }
301 return InfFE<Dim,T_radial,T_map>::shape(fet, inf_elem, i, p);
302}
303
304
305#ifdef LIBMESH_ENABLE_DEPRECATED
306template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
307Real InfFE<Dim,T_radial,T_map>::shape_deriv (const FEType & fe_t,
308 const ElemType inf_elem_type,
309 const unsigned int i,
310 const unsigned int j,
311 const Point & p)
312{
313 libmesh_deprecated();
314
315 libmesh_assert_not_equal_to (Dim, 0);
316 libmesh_assert_greater (Dim,j);
317#ifdef DEBUG
318 // this makes only sense when used for mapping
319 if ((T_radial != INFINITE_MAP) && !_warned_for_dshape)
320 {
321 libMesh::err << "WARNING: InfFE<Dim,T_radial,T_map>::shape_deriv(...) does _not_" << std::endl
322 << " return the correct trial function gradients! Use " << std::endl
323 << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!"
324 << std::endl;
325 _warned_for_dshape = true;
326 }
327#endif
328
329 const ElemType base_et (InfFEBase::get_elem_type(inf_elem_type));
330 const Order o_radial (fe_t.radial_order);
331 const Real v (p(Dim-1));
332
333 unsigned int i_base, i_radial;
334 compute_shape_indices(fe_t, inf_elem_type, i, i_base, i_radial);
335
336 if (j== Dim -1)
337 {
338 Real RadialDeriv = InfFE<Dim,T_radial,T_map>::eval_deriv(v, o_radial, i_radial)
339 * InfFERadial::decay(Dim,v)
340 + InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
341 * InfFERadial::decay_deriv(Dim, v);
342
343 return FEInterface::shape(Dim-1, fe_t, base_et, i_base, p)*RadialDeriv;
344 }
345
346 return FEInterface::shape_deriv(Dim-1, fe_t, base_et, i_base, j, p)
347 * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
348 * InfFERadial::decay(Dim,v);
349}
350#endif // LIBMESH_ENABLE_DEPRECATED
351
352
353template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
354Real InfFE<Dim,T_radial,T_map>::shape_deriv (const FEType & fe_t,
355 const Elem * inf_elem,
356 const unsigned int i,
357 const unsigned int j,
358 const Point & p)
359{
360 libmesh_assert_not_equal_to (Dim, 0);
361 libmesh_assert_greater (Dim,j);
362#ifdef DEBUG
363 // this makes only sense when used for mapping
364 if ((T_radial != INFINITE_MAP) && !_warned_for_dshape)
365 {
366 libMesh::err << "WARNING: InfFE<Dim,T_radial,T_map>::shape_deriv(...) does _not_" << std::endl
367 << " return the correct trial function gradients! Use " << std::endl
368 << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!"
369 << std::endl;
370 _warned_for_dshape = true;
371 }
372#endif
373 const Order o_radial (fe_t.radial_order);
374 const Real v (p(Dim-1));
375
376 std::unique_ptr<const Elem> base_el (inf_elem->build_side_ptr(0));
377
378 unsigned int i_base, i_radial;
379
380 if ((-1. > v ) || (v > 1.))
381 {
382 //TODO: This is for debugging. We should never come here.
383 // Therefore we can do very useless things then:
384 i_base=0;
385 }
386 compute_shape_indices(fe_t, inf_elem, i, i_base, i_radial);
387
388 if (j== Dim -1)
389 {
390 Real RadialDeriv = InfFE<Dim,T_radial,T_map>::eval_deriv(v, o_radial, i_radial)
391 * InfFERadial::decay(Dim,v)
392 + InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
393 * InfFERadial::decay_deriv(Dim,v);
394
395 return FEInterface::shape(fe_t, base_el.get(), i_base, p)*RadialDeriv;
396 }
397 return FEInterface::shape_deriv(fe_t, base_el.get(), i_base, j, p)
398 * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
399 * InfFERadial::decay(Dim,v);
400}
401
402
403
404
405
406template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
407void InfFE<Dim,T_radial,T_map>::compute_data(const FEType & fet,
408 const Elem * inf_elem,
409 FEComputeData & data)
410{
411 libmesh_assert(inf_elem);
412 libmesh_assert_not_equal_to (Dim, 0);
413
414 const Order o_radial (fet.radial_order);
415 const Order radial_mapping_order (InfFERadial::mapping_order());
416 const Point & p (data.p);
417 const Real v (p(Dim-1));
418 std::unique_ptr<const Elem> base_el (inf_elem->build_side_ptr(0));
419
420 /*
421 * compute \p interpolated_dist containing the mapping-interpolated
422 * distance of the base point to the origin. This is the same
423 * for all shape functions. Set \p interpolated_dist to 0, it
424 * is added to.
425 */
426 Real interpolated_dist = 0.;
427 switch (Dim)
428 {
429 case 1:
430 {
431 libmesh_assert_equal_to (inf_elem->type(), INFEDGE2);
432 interpolated_dist = Point(inf_elem->point(0) - inf_elem->point(1)).norm();
433 break;
434 }
435
436 case 2:
437 {
438 const unsigned int n_base_nodes = base_el->n_nodes();
439
440 const Point origin = inf_elem->origin();
441 const Order base_mapping_order (base_el->default_order());
442 const ElemType base_mapping_elem_type (base_el->type());
443
444 // interpolate the base nodes' distances
445 for (unsigned int n=0; n<n_base_nodes; n++)
446 interpolated_dist += Point(base_el->point(n) - origin).norm()
447 * FE<1,LAGRANGE>::shape (base_mapping_elem_type, base_mapping_order, n, p);
448 break;
449 }
450
451 case 3:
452 {
453 const unsigned int n_base_nodes = base_el->n_nodes();
454
455 const Point origin = inf_elem->origin();
456 const Order base_mapping_order (base_el->default_order());
457 const ElemType base_mapping_elem_type (base_el->type());
458
459 // interpolate the base nodes' distances
460 for (unsigned int n=0; n<n_base_nodes; n++)
461 interpolated_dist += Point(base_el->point(n) - origin).norm()
462 * FE<2,LAGRANGE>::shape (base_mapping_elem_type, base_mapping_order, n, p);
463 break;
464 }
465
466 default:
467 libmesh_error_msg("Unknown Dim = " << Dim);
468 }
469
470
471 const Real speed = data.speed;
472
473 //TODO: I find it inconvenient to have a quantity phase which is phase/speed.
474 // But it might be better than redefining a quantities meaning.
475 data.phase = interpolated_dist /* together with next line: */
476 * InfFE<Dim,INFINITE_MAP,T_map>::eval(v, radial_mapping_order, 1)/speed; /* phase(s,t,v)/c */
477
478 // We assume time-harmonic behavior in this function!
479
480#ifdef LIBMESH_USE_COMPLEX_NUMBERS
481 // the wave number
482 const Number wavenumber = 2. * libMesh::pi * data.frequency / speed;
483
484 // the exponent for time-harmonic behavior
485 // \note: this form is much less general than the implementation of dphase, which can be easily extended to
486 // other forms than e^{i kr}.
487 const Number exponent = imaginary /* imaginary unit */
488 * wavenumber /* k (can be complex) */
489 * data.phase*speed;
490
491 const Number time_harmonic = exp(exponent); /* e^(i*k*phase(s,t,v)) */
492#else
493 const Number time_harmonic = 1;
494#endif //LIBMESH_USE_COMPLEX_NUMBERS
495
496 /*
497 * compute \p shape for all dof in the element
498 */
499 if (Dim > 1)
500 {
501 const unsigned int n_dof = n_dofs (fet, inf_elem);
502 data.shape.resize(n_dof);
503 if (data.need_derivative())
504 {
505 data.dshape.resize(n_dof);
506 data.local_transform.resize(Dim);
507
508 for (unsigned int d=0; d<Dim; d++)
509 data.local_transform[d].resize(Dim);
510
511 // compute the reference->physical map at the point \p p.
512 // Use another fe_map to avoid interference with \p this->_fe_map
513 // which is initialized at the quadrature points...
514 auto fe = FEBase::build_InfFE(Dim, fet);
515 std::vector<Point> pt = {p};
516 fe->get_dxidx(); // to compute the map
517 fe->reinit(inf_elem, &pt);
518
519 // compute the reference->physical map.
520 data.local_transform[0][0] = fe->get_dxidx()[0];
521 data.local_transform[1][0] = fe->get_detadx()[0];
522 data.local_transform[1][1] = fe->get_detady()[0];
523 data.local_transform[0][1] = fe->get_dxidy()[0];
524 if (Dim > 2)
525 {
526 data.local_transform[2][0] = fe->get_dzetadx()[0];
527 data.local_transform[2][1] = fe->get_dzetady()[0];
528 data.local_transform[2][2] = fe->get_dzetadz()[0];
529 data.local_transform[1][2] = fe->get_detadz()[0];
530 data.local_transform[0][2] = fe->get_dxidz()[0];
531 }
532 } // endif data.need_derivative()
533
534 for (unsigned int i=0; i<n_dof; i++)
535 {
536 // compute base and radial shape indices
537 unsigned int i_base, i_radial;
538 compute_shape_indices(fet, inf_elem, i, i_base, i_radial);
539
540 data.shape[i] = (InfFERadial::decay(Dim,v) /* (1.-v)/2. in 3D */
541 * FEInterface::shape(fet, base_el.get(), i_base, p) /* S_n(s,t) */
542 * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)) /* L_n(v) */
543 * time_harmonic; /* e^(i*k*phase(s,t,v) */
544
545 // use differentiation of the above equation
546 if (data.need_derivative())
547 {
548 data.dshape[i](0) = (InfFERadial::decay(Dim,v)
549 * FEInterface::shape_deriv(fet, base_el.get(), i_base, 0, p)
550 * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial))
551 * time_harmonic;
552
553 if (Dim > 2)
554 {
555 data.dshape[i](1) = (InfFERadial::decay(Dim,v)
556 * FEInterface::shape_deriv(fet, base_el.get(), i_base, 1, p)
557 * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial))
558 * time_harmonic;
559
560 }
561 data.dshape[i](Dim-1) = (InfFERadial::decay_deriv(Dim, v) * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
562 +InfFERadial::decay(Dim,v) * InfFE<Dim,T_radial,T_map>::eval_deriv(v, o_radial, i_radial))
563 * FEInterface::shape(fet, base_el.get(), i_base, p) * time_harmonic;
564
565#ifdef LIBMESH_USE_COMPLEX_NUMBERS
566 // derivative of time_harmonic (works for harmonic behavior only):
567 data.dshape[i](Dim-1)+= data.shape[i]*imaginary*wavenumber
568 *interpolated_dist*InfFE<Dim,INFINITE_MAP,T_map>::eval_deriv(v, radial_mapping_order, 1);
569
570#else
571 /*
572 * The gradient in infinite elements is dominated by the contribution due to the oscillating phase.
573 * Since this term is imaginary, I think there is no means to look at it without having complex numbers.
574 */
575 libmesh_not_implemented();
576 // Maybe we can solve it with a warning as well, but I think one really should not do this...
577#endif
578 }
579 }
580 }
581
582 else
583 libmesh_error_msg("compute_data() for 1-dimensional InfFE not implemented.");
584}
585
586
587template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
588Real InfFE<Dim,T_radial,T_map>::shape_deriv(const FEType fet,
589 const Elem * inf_elem,
590 const unsigned int i,
591 const unsigned int j,
592 const Point & p,
593 const bool add_p_level)
594{
595 if (add_p_level)
596 {
597 FEType tmp_fet=fet;
598 tmp_fet = fet.order + add_p_level * inf_elem->p_level();
599 return InfFE<Dim,T_radial,T_map>::shape_deriv(tmp_fet, inf_elem, i, j, p);
600 }
601 return InfFE<Dim,T_radial,T_map>::shape_deriv(fet, inf_elem, i, j, p);
602}
603
604
605
606
607
608
609template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
610void InfFE<Dim,T_radial,T_map>::compute_node_indices (const ElemType inf_elem_type,
611 const unsigned int outer_node_index,
612 unsigned int & base_node,
613 unsigned int & radial_node)
614{
615 switch (inf_elem_type)
616 {
617 case INFEDGE2:
618 {
619 libmesh_assert_less (outer_node_index, 2);
620 base_node = 0;
621 radial_node = outer_node_index;
622 return;
623 }
624
625
626 // linear base approximation, easy to determine
627 case INFQUAD4:
628 {
629 libmesh_assert_less (outer_node_index, 4);
630 base_node = outer_node_index % 2;
631 radial_node = outer_node_index / 2;
632 return;
633 }
634
635 case INFPRISM6:
636 {
637 libmesh_assert_less (outer_node_index, 6);
638 base_node = outer_node_index % 3;
639 radial_node = outer_node_index / 3;
640 return;
641 }
642
643 case INFHEX8:
644 {
645 libmesh_assert_less (outer_node_index, 8);
646 base_node = outer_node_index % 4;
647 radial_node = outer_node_index / 4;
648 return;
649 }
650
651
652 // higher order base approximation, more work necessary
653 case INFQUAD6:
654 {
655 switch (outer_node_index)
656 {
657 case 0:
658 case 1:
659 {
660 radial_node = 0;
661 base_node = outer_node_index;
662 return;
663 }
664
665 case 2:
666 case 3:
667 {
668 radial_node = 1;
669 base_node = outer_node_index-2;
670 return;
671 }
672
673 case 4:
674 {
675 radial_node = 0;
676 base_node = 2;
677 return;
678 }
679
680 case 5:
681 {
682 radial_node = 1;
683 base_node = 2;
684 return;
685 }
686
687 default:
688 libmesh_error_msg("Unrecognized outer_node_index = " << outer_node_index);
689 }
690 }
691
692
693 case INFHEX16:
694 case INFHEX18:
695 {
696 switch (outer_node_index)
697 {
698 case 0:
699 case 1:
700 case 2:
701 case 3:
702 {
703 radial_node = 0;
704 base_node = outer_node_index;
705 return;
706 }
707
708 case 4:
709 case 5:
710 case 6:
711 case 7:
712 {
713 radial_node = 1;
714 base_node = outer_node_index-4;
715 return;
716 }
717
718 case 8:
719 case 9:
720 case 10:
721 case 11:
722 {
723 radial_node = 0;
724 base_node = outer_node_index-4;
725 return;
726 }
727
728 case 12:
729 case 13:
730 case 14:
731 case 15:
732 {
733 radial_node = 1;
734 base_node = outer_node_index-8;
735 return;
736 }
737
738 case 16:
739 {
740 libmesh_assert_equal_to (inf_elem_type, INFHEX18);
741 radial_node = 0;
742 base_node = 8;
743 return;
744 }
745
746 case 17:
747 {
748 libmesh_assert_equal_to (inf_elem_type, INFHEX18);
749 radial_node = 1;
750 base_node = 8;
751 return;
752 }
753
754 default:
755 libmesh_error_msg("Unrecognized outer_node_index = " << outer_node_index);
756 }
757 }
758
759
760 case INFPRISM12:
761 {
762 switch (outer_node_index)
763 {
764 case 0:
765 case 1:
766 case 2:
767 {
768 radial_node = 0;
769 base_node = outer_node_index;
770 return;
771 }
772
773 case 3:
774 case 4:
775 case 5:
776 {
777 radial_node = 1;
778 base_node = outer_node_index-3;
779 return;
780 }
781
782 case 6:
783 case 7:
784 case 8:
785 {
786 radial_node = 0;
787 base_node = outer_node_index-3;
788 return;
789 }
790
791 case 9:
792 case 10:
793 case 11:
794 {
795 radial_node = 1;
796 base_node = outer_node_index-6;
797 return;
798 }
799
800 default:
801 libmesh_error_msg("Unrecognized outer_node_index = " << outer_node_index);
802 }
803 }
804
805
806 default:
807 libmesh_error_msg("ERROR: Bad infinite element type=" << Utility::enum_to_string(inf_elem_type) << ", node=" << outer_node_index);
808 }
809}
810
811
812
813
814
815
816template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
817void InfFE<Dim,T_radial,T_map>::compute_node_indices_fast (const ElemType inf_elem_type,
818 const unsigned int outer_node_index,
819 unsigned int & base_node,
820 unsigned int & radial_node)
821{
822 libmesh_assert_not_equal_to (inf_elem_type, INVALID_ELEM);
823
824 static std::vector<unsigned int> _static_base_node_index;
825 static std::vector<unsigned int> _static_radial_node_index;
826
827 /*
828 * fast counterpart to compute_node_indices(), uses local static buffers
829 * to store the index maps. The class member
830 * \p _compute_node_indices_fast_current_elem_type remembers
831 * the current element type.
832 *
833 * Note that there exist non-static members storing the
834 * same data. However, you never know what element type
835 * is currently used by the \p InfFE object, and what
836 * request is currently directed to the static \p InfFE
837 * members (which use \p compute_node_indices_fast()).
838 * So separate these.
839 *
840 * check whether the work for this elemtype has already
841 * been done. If so, use this index. Otherwise, refresh
842 * the buffer to this element type.
843 */
844 if (inf_elem_type==_compute_node_indices_fast_current_elem_type)
845 {
846 base_node = _static_base_node_index [outer_node_index];
847 radial_node = _static_radial_node_index[outer_node_index];
848 return;
849 }
850 else
851 {
852 // store the map for _all_ nodes for this element type
853 _compute_node_indices_fast_current_elem_type = inf_elem_type;
854
855 unsigned int n_nodes = libMesh::invalid_uint;
856
857 switch (inf_elem_type)
858 {
859 case INFEDGE2:
860 {
861 n_nodes = 2;
862 break;
863 }
864 case INFQUAD4:
865 {
866 n_nodes = 4;
867 break;
868 }
869 case INFQUAD6:
870 {
871 n_nodes = 6;
872 break;
873 }
874 case INFHEX8:
875 {
876 n_nodes = 8;
877 break;
878 }
879 case INFHEX16:
880 {
881 n_nodes = 16;
882 break;
883 }
884 case INFHEX18:
885 {
886 n_nodes = 18;
887 break;
888 }
889 case INFPRISM6:
890 {
891 n_nodes = 6;
892 break;
893 }
894 case INFPRISM12:
895 {
896 n_nodes = 12;
897 break;
898 }
899 default:
900 libmesh_error_msg("ERROR: Bad infinite element type=" << Utility::enum_to_string(inf_elem_type) << ", node=" << outer_node_index);
901 }
902
903
904 _static_base_node_index.resize (n_nodes);
905 _static_radial_node_index.resize(n_nodes);
906
907 for (unsigned int n=0; n<n_nodes; n++)
908 compute_node_indices (inf_elem_type,
909 n,
910 _static_base_node_index [outer_node_index],
911 _static_radial_node_index[outer_node_index]);
912
913 // and return for the specified node
914 base_node = _static_base_node_index [outer_node_index];
915 radial_node = _static_radial_node_index[outer_node_index];
916 return;
917 }
918}
919
920
921
922
923
924
925template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
926void InfFE<Dim,T_radial,T_map>::compute_shape_indices (const FEType & fet,
927 const Elem * inf_elem,
928 const unsigned int i,
929 unsigned int & base_shape,
930 unsigned int & radial_shape)
931{
932 // An example is provided: the numbers in comments refer to
933 // a fictitious InfHex18. The numbers are chosen as exemplary
934 // values. There is currently no base approximation that
935 // requires this many dof's at nodes, sides, faces and in the element.
936 //
937 // the order of the shape functions is heavily related with the
938 // order the dofs are assigned in \p DofMap::distributed_dofs().
939 // Due to the infinite elements with higher-order base approximation,
940 // some more effort is necessary.
941 //
942 // numbering scheme:
943 // 1. all vertices in the base, assign node->n_comp() dofs to each vertex
944 // 2. all vertices further out: innermost loop: radial shapes,
945 // then the base approximation shapes
946 // 3. all side nodes in the base, assign node->n_comp() dofs to each side node
947 // 4. all side nodes further out: innermost loop: radial shapes,
948 // then the base approximation shapes
949 // 5. (all) face nodes in the base, assign node->n_comp() dofs to each face node
950 // 6. (all) face nodes further out: innermost loop: radial shapes,
951 // then the base approximation shapes
952 // 7. element-associated dof in the base
953 // 8. element-associated dof further out
954
955 const unsigned int radial_order = static_cast<unsigned int>(fet.radial_order.get_order()); // 4
956 const unsigned int radial_order_p_one = radial_order+1; // 5
957
958 std::unique_ptr<const Elem> base_elem = InfFEBase::build_elem(inf_elem); // QUAD9
959
960 // assume that the number of dof is the same for all vertices
961 unsigned int n_base_vertices = libMesh::invalid_uint; // 4
962 const unsigned int n_base_vertex_dof = FEInterface::n_dofs_at_node (fet, base_elem.get(), 0); // 2
963
964 unsigned int n_base_side_nodes = libMesh::invalid_uint; // 4
965 unsigned int n_base_side_dof = libMesh::invalid_uint; // 3
966
967 unsigned int n_base_face_nodes = libMesh::invalid_uint; // 1
968 unsigned int n_base_face_dof = libMesh::invalid_uint; // 5
969
970 const unsigned int n_base_elem_dof = FEInterface::n_dofs_per_elem (fet, base_elem.get()); // 9
971
972 const ElemType inf_elem_type = inf_elem->type();
973
974 switch (inf_elem_type)
975 {
976 case INFEDGE2:
977 {
978 n_base_vertices = 1;
979 n_base_side_nodes = 0;
980 n_base_face_nodes = 0;
981 n_base_side_dof = 0;
982 n_base_face_dof = 0;
983 break;
984 }
985
986 case INFQUAD4:
987 {
988 n_base_vertices = 2;
989 n_base_side_nodes = 0;
990 n_base_face_nodes = 0;
991 n_base_side_dof = 0;
992 n_base_face_dof = 0;
993 break;
994 }
995
996 case INFQUAD6:
997 {
998 n_base_vertices = 2;
999 n_base_side_nodes = 1;
1000 n_base_face_nodes = 0;
1001 n_base_side_dof = FEInterface::n_dofs_at_node (fet, base_elem.get(), n_base_vertices);
1002 n_base_face_dof = 0;
1003 break;
1004 }
1005
1006 case INFHEX8:
1007 {
1008 n_base_vertices = 4;
1009 n_base_side_nodes = 0;
1010 n_base_face_nodes = 0;
1011 n_base_side_dof = 0;
1012 n_base_face_dof = 0;
1013 break;
1014 }
1015
1016 case INFHEX16:
1017 {
1018 n_base_vertices = 4;
1019 n_base_side_nodes = 4;
1020 n_base_face_nodes = 0;
1021 n_base_side_dof = FEInterface::n_dofs_at_node (fet, base_elem.get(), n_base_vertices);
1022 n_base_face_dof = 0;
1023 break;
1024 }
1025
1026 case INFHEX18:
1027 {
1028 n_base_vertices = 4;
1029 n_base_side_nodes = 4;
1030 n_base_face_nodes = 1;
1031 n_base_side_dof = FEInterface::n_dofs_at_node (fet, base_elem.get(), n_base_vertices);
1032 n_base_face_dof = FEInterface::n_dofs_at_node (fet, base_elem.get(), 8);
1033 break;
1034 }
1035
1036
1037 case INFPRISM6:
1038 {
1039 n_base_vertices = 3;
1040 n_base_side_nodes = 0;
1041 n_base_face_nodes = 0;
1042 n_base_side_dof = 0;
1043 n_base_face_dof = 0;
1044 break;
1045 }
1046
1047 case INFPRISM12:
1048 {
1049 n_base_vertices = 3;
1050 n_base_side_nodes = 3;
1051 n_base_face_nodes = 0;
1052 n_base_side_dof = FEInterface::n_dofs_at_node (fet, base_elem.get(), n_base_vertices);
1053 n_base_face_dof = 0;
1054 break;
1055 }
1056
1057 default:
1058 libmesh_error_msg("Unrecognized inf_elem type = " << Utility::enum_to_string(inf_elem_type));
1059 }
1060
1061
1062 {
1063 // these are the limits describing the intervals where the shape function lies
1064 const unsigned int n_dof_at_base_vertices = n_base_vertices*n_base_vertex_dof; // 8
1065 const unsigned int n_dof_at_all_vertices = n_dof_at_base_vertices*radial_order_p_one; // 40
1066
1067 const unsigned int n_dof_at_base_sides = n_base_side_nodes*n_base_side_dof; // 12
1068 const unsigned int n_dof_at_all_sides = n_dof_at_base_sides*radial_order_p_one; // 60
1069
1070 const unsigned int n_dof_at_base_face = n_base_face_nodes*n_base_face_dof; // 5
1071 const unsigned int n_dof_at_all_faces = n_dof_at_base_face*radial_order_p_one; // 25
1072
1073
1074 // start locating the shape function
1075 if (i < n_dof_at_base_vertices) // range of i: 0..7
1076 {
1077 // belongs to vertex in the base
1078 radial_shape = 0;
1079 base_shape = i;
1080 }
1081
1082 else if (i < n_dof_at_all_vertices) // range of i: 8..39
1083 {
1084 /* belongs to vertex in the outer shell
1085 *
1086 * subtract the number of dof already counted,
1087 * so that i_offset contains only the offset for the base
1088 */
1089 const unsigned int i_offset = i - n_dof_at_base_vertices; // 0..31
1090
1091 // first the radial dof are counted, then the base dof
1092 radial_shape = (i_offset % radial_order) + 1;
1093 base_shape = i_offset / radial_order;
1094 }
1095
1096 else if (i < n_dof_at_all_vertices+n_dof_at_base_sides) // range of i: 40..51
1097 {
1098 // belongs to base, is a side node
1099 radial_shape = 0;
1100 base_shape = i - radial_order * n_dof_at_base_vertices; // 8..19
1101 }
1102
1103 else if (i < n_dof_at_all_vertices+n_dof_at_all_sides) // range of i: 52..99
1104 {
1105 // belongs to side node in the outer shell
1106 const unsigned int i_offset = i - (n_dof_at_all_vertices
1107 + n_dof_at_base_sides); // 0..47
1108 radial_shape = (i_offset % radial_order) + 1;
1109 base_shape = (i_offset / radial_order) + n_dof_at_base_vertices;
1110 }
1111
1112 else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_base_face) // range of i: 100..104
1113 {
1114 // belongs to the node in the base face
1115 radial_shape = 0;
1116 base_shape = i - radial_order*(n_dof_at_base_vertices
1117 + n_dof_at_base_sides); // 20..24
1118 }
1119
1120 else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_all_faces) // range of i: 105..124
1121 {
1122 // belongs to the node in the outer face
1123 const unsigned int i_offset = i - (n_dof_at_all_vertices
1124 + n_dof_at_all_sides
1125 + n_dof_at_base_face); // 0..19
1126 radial_shape = (i_offset % radial_order) + 1;
1127 base_shape = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides;
1128 }
1129
1130 else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_all_faces+n_base_elem_dof) // range of i: 125..133
1131 {
1132 // belongs to the base and is an element associated shape
1133 radial_shape = 0;
1134 base_shape = i - (n_dof_at_all_vertices
1135 + n_dof_at_all_sides
1136 + n_dof_at_all_faces); // 0..8
1137 }
1138
1139 else // range of i: 134..169
1140 {
1141 libmesh_assert_less (i, n_dofs(fet, inf_elem));
1142 // belongs to the outer shell and is an element associated shape
1143 const unsigned int i_offset = i - (n_dof_at_all_vertices
1144 + n_dof_at_all_sides
1145 + n_dof_at_all_faces
1146 + n_base_elem_dof); // 0..19
1147 radial_shape = (i_offset % radial_order) + 1;
1148 base_shape = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides + n_dof_at_base_face;
1149 }
1150 }
1151
1152 return;
1153}
1154
1155
1156
1157#ifdef LIBMESH_ENABLE_DEPRECATED
1158template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
1159void InfFE<Dim,T_radial,T_map>::compute_shape_indices (const FEType & fet,
1160 const ElemType inf_elem_type,
1161 const unsigned int i,
1162 unsigned int & base_shape,
1163 unsigned int & radial_shape)
1164{
1165 // This is basically cut-and-paste copy of the non-deprecated code
1166 // above, and should be removed eventually.
1167 libmesh_deprecated();
1168
1169 const unsigned int radial_order = static_cast<unsigned int>(fet.radial_order.get_order()); // 4
1170 const unsigned int radial_order_p_one = radial_order+1; // 5
1171
1172 const ElemType base_elem_type (InfFEBase::get_elem_type(inf_elem_type)); // QUAD9
1173
1174 // assume that the number of dof is the same for all vertices
1175 unsigned int n_base_vertices = libMesh::invalid_uint; // 4
1176 const unsigned int n_base_vertex_dof = FEInterface::n_dofs_at_node (Dim-1, fet, base_elem_type, 0);// 2
1177
1178 unsigned int n_base_side_nodes = libMesh::invalid_uint; // 4
1179 unsigned int n_base_side_dof = libMesh::invalid_uint; // 3
1180
1181 unsigned int n_base_face_nodes = libMesh::invalid_uint; // 1
1182 unsigned int n_base_face_dof = libMesh::invalid_uint; // 5
1183
1184 const unsigned int n_base_elem_dof = FEInterface::n_dofs_per_elem (Dim-1, fet, base_elem_type);// 9
1185
1186
1187 switch (inf_elem_type)
1188 {
1189 case INFEDGE2:
1190 {
1191 n_base_vertices = 1;
1192 n_base_side_nodes = 0;
1193 n_base_face_nodes = 0;
1194 n_base_side_dof = 0;
1195 n_base_face_dof = 0;
1196 break;
1197 }
1198
1199 case INFQUAD4:
1200 {
1201 n_base_vertices = 2;
1202 n_base_side_nodes = 0;
1203 n_base_face_nodes = 0;
1204 n_base_side_dof = 0;
1205 n_base_face_dof = 0;
1206 break;
1207 }
1208
1209 case INFQUAD6:
1210 {
1211 n_base_vertices = 2;
1212 n_base_side_nodes = 1;
1213 n_base_face_nodes = 0;
1214 n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
1215 n_base_face_dof = 0;
1216 break;
1217 }
1218
1219 case INFHEX8:
1220 {
1221 n_base_vertices = 4;
1222 n_base_side_nodes = 0;
1223 n_base_face_nodes = 0;
1224 n_base_side_dof = 0;
1225 n_base_face_dof = 0;
1226 break;
1227 }
1228
1229 case INFHEX16:
1230 {
1231 n_base_vertices = 4;
1232 n_base_side_nodes = 4;
1233 n_base_face_nodes = 0;
1234 n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
1235 n_base_face_dof = 0;
1236 break;
1237 }
1238
1239 case INFHEX18:
1240 {
1241 n_base_vertices = 4;
1242 n_base_side_nodes = 4;
1243 n_base_face_nodes = 1;
1244 n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
1245 n_base_face_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, 8);
1246 break;
1247 }
1248
1249
1250 case INFPRISM6:
1251 {
1252 n_base_vertices = 3;
1253 n_base_side_nodes = 0;
1254 n_base_face_nodes = 0;
1255 n_base_side_dof = 0;
1256 n_base_face_dof = 0;
1257 break;
1258 }
1259
1260 case INFPRISM12:
1261 {
1262 n_base_vertices = 3;
1263 n_base_side_nodes = 3;
1264 n_base_face_nodes = 0;
1265 n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
1266 n_base_face_dof = 0;
1267 break;
1268 }
1269
1270 default:
1271 libmesh_error_msg("Unrecognized inf_elem_type = " << Utility::enum_to_string(inf_elem_type));
1272 }
1273
1274
1275 {
1276 // these are the limits describing the intervals where the shape function lies
1277 const unsigned int n_dof_at_base_vertices = n_base_vertices*n_base_vertex_dof; // 8
1278 const unsigned int n_dof_at_all_vertices = n_dof_at_base_vertices*radial_order_p_one; // 40
1279
1280 const unsigned int n_dof_at_base_sides = n_base_side_nodes*n_base_side_dof; // 12
1281 const unsigned int n_dof_at_all_sides = n_dof_at_base_sides*radial_order_p_one; // 60
1282
1283 const unsigned int n_dof_at_base_face = n_base_face_nodes*n_base_face_dof; // 5
1284 const unsigned int n_dof_at_all_faces = n_dof_at_base_face*radial_order_p_one; // 25
1285
1286
1287 // start locating the shape function
1288 if (i < n_dof_at_base_vertices) // range of i: 0..7
1289 {
1290 // belongs to vertex in the base
1291 radial_shape = 0;
1292 base_shape = i;
1293 }
1294
1295 else if (i < n_dof_at_all_vertices) // range of i: 8..39
1296 {
1297 /* belongs to vertex in the outer shell
1298 *
1299 * subtract the number of dof already counted,
1300 * so that i_offset contains only the offset for the base
1301 */
1302 const unsigned int i_offset = i - n_dof_at_base_vertices; // 0..31
1303
1304 // first the radial dof are counted, then the base dof
1305 radial_shape = (i_offset % radial_order) + 1;
1306 base_shape = i_offset / radial_order;
1307 }
1308
1309 else if (i < n_dof_at_all_vertices+n_dof_at_base_sides) // range of i: 40..51
1310 {
1311 // belongs to base, is a side node
1312 radial_shape = 0;
1313 base_shape = i - radial_order * n_dof_at_base_vertices; // 8..19
1314 }
1315
1316 else if (i < n_dof_at_all_vertices+n_dof_at_all_sides) // range of i: 52..99
1317 {
1318 // belongs to side node in the outer shell
1319 const unsigned int i_offset = i - (n_dof_at_all_vertices
1320 + n_dof_at_base_sides); // 0..47
1321 radial_shape = (i_offset % radial_order) + 1;
1322 base_shape = (i_offset / radial_order) + n_dof_at_base_vertices;
1323 }
1324
1325 else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_base_face) // range of i: 100..104
1326 {
1327 // belongs to the node in the base face
1328 radial_shape = 0;
1329 base_shape = i - radial_order*(n_dof_at_base_vertices
1330 + n_dof_at_base_sides); // 20..24
1331 }
1332
1333 else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_all_faces) // range of i: 105..124
1334 {
1335 // belongs to the node in the outer face
1336 const unsigned int i_offset = i - (n_dof_at_all_vertices
1337 + n_dof_at_all_sides
1338 + n_dof_at_base_face); // 0..19
1339 radial_shape = (i_offset % radial_order) + 1;
1340 base_shape = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides;
1341 }
1342
1343 else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_all_faces+n_base_elem_dof) // range of i: 125..133
1344 {
1345 // belongs to the base and is an element associated shape
1346 radial_shape = 0;
1347 base_shape = i - (n_dof_at_all_vertices
1348 + n_dof_at_all_sides
1349 + n_dof_at_all_faces); // 0..8
1350 }
1351
1352 else // range of i: 134..169
1353 {
1354 // belongs to the outer shell and is an element associated shape
1355 const unsigned int i_offset = i - (n_dof_at_all_vertices
1356 + n_dof_at_all_sides
1357 + n_dof_at_all_faces
1358 + n_base_elem_dof); // 0..19
1359 radial_shape = (i_offset % radial_order) + 1;
1360 base_shape = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides + n_dof_at_base_face;
1361 }
1362 }
1363
1364 return;
1365}
1366#endif // LIBMESH_ENABLE_DEPRECATED
1367
1368
1369#ifdef LIBMESH_ENABLE_AMR
1370#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1371
1372template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
1373void InfFE<Dim, T_radial, T_map>::inf_compute_node_constraints (NodeConstraints & /* constraints */,const Elem * elem)
1374{
1375 // only constrain elements in 2,3d.
1376 if (Dim == 1)
1377 return;
1378
1379 libmesh_assert(elem);
1380
1381 // only constrain active and ancestor elements
1382 if (elem->subactive())
1383 return;
1384
1385 // for infinite elements, the computation of constraints is somewhat different
1386 // than for Lagrange elements:
1387 // 1) Only the base element (i.e. side(0) ) may be refined.
1388 // Thus, in radial direction no constraints must be considered.
1389 // 2) Due to the tensorial structure of shape functions (base_shape * radial_function),
1390 // it must be ensured that all element DOFs inherit that constraint.
1391 // Consequently, the constraints are computed on the base (baseh_shape) but must
1392 // be applied to all DOFs with the respective base_shape index (i.e. for all radial_functions).
1393 //
1394 // FIXME: In the current form, this function does not work for infinite elements
1395 // because constraining the non-base points requires knowledge of the T_map and T_radial
1396 // parameters; but they are not accessible via the element and may differ between variables.
1397 //
1398 // For the moment being, we just check if this element can be skipped and fail otherwise.
1399
1400 // if one of the sides needs a constraint, an error is thrown.
1401 // In other cases, we leave the function regularly.
1402 for (auto s : elem->side_index_range())
1403 {
1404 if (elem->neighbor_ptr(s) != nullptr &&
1405 elem->neighbor_ptr(s) != remote_elem)
1406 if (elem->neighbor_ptr(s)->level() < elem->level())
1407 {
1408 libmesh_not_implemented();
1409 }
1410 }
1411}
1412
1413#endif //LIBMESH_ENABLE_NODE_CONSTRAINTS
1414
1415
1416template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
1417void InfFE<Dim, T_radial, T_map>::inf_compute_constraints (DofConstraints & constraints,
1418 DofMap & dof_map,
1419 const unsigned int variable_number,
1420 const Elem * child_elem)
1421{
1422
1423 // only constrain elements in 2,3d.
1424 if (Dim == 1)
1425 return;
1426
1427 libmesh_assert(child_elem);
1428
1429 // only constrain active and ancestor elements
1430 if (child_elem->subactive())
1431 return;
1432
1433 // Before we start to compute anything, lets check if any confinement is needed:
1434 bool need_constraints=false;
1435 for (auto child_neighbor : child_elem->neighbor_ptr_range())
1436 if (child_neighbor->level() < child_elem->level())
1437 {
1438 need_constraints = true;
1439 break;
1440 }
1441 if (!need_constraints)
1442 return;
1443
1444 // For infinite elements, the computation of constraints is somewhat different
1445 // than for Lagrange elements:
1446 // 1) When an infinite element is refined, only the base element (i.e. side(0) ) is refined.
1447 //
1448 // 2) Due to the tensorial structure of shape functions (base_shape * radial_function),
1449 // it must be ensured that all element DOFs inherit that constraint.
1450 // It is important here to distinguish the (total) DOF from base DOF and radial DOF contributions.
1451 //
1452 // 3) Due to the generality of radial polynomial (of type fe_type.radial_family and with order fe_type.radial_order)
1453 // here basis functions cannot be mapped to nodes: Independent from the radial polynomial,
1454 // infinite elements have one set of nodes at the base (side(0)) and a second set at twice the distance to their origin.
1455 //
1456 // Independent from the polynomial and degree used, the first radial DOF is 1 at the base while all others are 0 there
1457 //
1458 //Constraining of DOFs is only needed when a DOF is nonzero at the elements face shared with a coarser element.
1459 // Thus, the following scheme is used here:
1460 //
1461 // -If the coarser element is the neighbor(0) (i.e. we share only the base), we must constrain
1462 // all DOFs that correspond to the first order radial contribution.
1463 // -if an infinite neighbor is coarser (than 'child_elem'), all those DOFs must be constrained
1464 // whose contribution from the base is non-zero at the interface.
1465 // In this case, we lack a point-assignement between DOFs and nodes, but since there is no refinement in radial direction,
1466 // the radial polynomials coincide on neighboring elements.
1467 // Thus, if one constraines these DOFs at one (arbitrary) point correctly, they match for each point along the radial direction.
1468 // Hence, we constrain them with the same values as those DOFs belonging to the first order polynomial, obtaining consistent
1469 // constraints that mimic constraints that are computed at the support points for each radial polynomial contribution.
1470
1471 FEType fe_type = dof_map.variable_type(variable_number);
1472
1473 libmesh_assert(fe_type.family == LAGRANGE);
1474
1475 std::vector<dof_id_type> child_base_dof_indices, parent_base_dof_indices;
1476 std::vector<dof_id_type> child_elem_dof_indices, parent_elem_dof_indices;
1477
1478 const Elem * parent_elem = child_elem->parent();
1479
1480 // This can't happen... Only level-0 elements have nullptr
1481 // parents, and no level-0 elements can be at a higher
1482 // level than their neighbors!
1483 libmesh_assert(parent_elem);
1484
1485 dof_map.dof_indices (child_elem, child_elem_dof_indices,
1486 variable_number);
1487 dof_map.dof_indices (parent_elem, parent_elem_dof_indices,
1488 variable_number);
1489
1490 const unsigned int n_total_dofs = child_elem_dof_indices.size();
1491 // fill the elements shape index map: we will have to use it later
1492 // to find the elements dofs that correspond to certain base_elem_dofs.
1493 std::vector<unsigned int> radial_shape_index(n_total_dofs);
1494 std::vector<unsigned int> base_shape_index(n_total_dofs);
1495 // fill the shape index map
1496#ifdef DEBUG
1497 unsigned int max_base_id=0;
1498 unsigned int max_radial_id=0;
1499#endif
1500 for (unsigned int n=0; n<n_total_dofs; ++n)
1501 {
1502 compute_shape_indices (fe_type,
1503 child_elem,
1504 n,
1505 base_shape_index[n],
1506 radial_shape_index[n]);
1507
1508#ifdef DEBUG
1509 if (base_shape_index[n] > max_base_id)
1510 max_base_id = base_shape_index[n];
1511 if (radial_shape_index[n] > max_radial_id)
1512 max_radial_id = radial_shape_index[n];
1513#endif
1514 }
1515
1516#ifdef DEBUG
1517 libmesh_assert_equal_to( (max_base_id+1)*(max_radial_id+1), n_total_dofs );
1518#endif
1519
1520 for (auto s : child_elem->side_index_range())
1521 if (child_elem->neighbor_ptr(s) != nullptr &&
1522 child_elem->neighbor_ptr(s) != remote_elem)
1523 if (child_elem->neighbor_ptr(s)->level() < child_elem->level())
1524 {
1525 // we ALWAYS take the base element for reference:
1526 // - For s=0, we refine all dofs with `radial_shape_index == 0
1527 // - for s>0, we refine all dofs whose corresponding base_shape has its support point shared with neighbor(s)
1528 std::unique_ptr<const Elem> child_base, parent_base;
1529 child_elem->build_side_ptr(child_base, 0);
1530 parent_elem->build_side_ptr(parent_base, 0);
1531
1532 const unsigned int n_base_dofs =
1533 FEInterface::n_dofs(fe_type, child_base.get());
1534
1535 // We need global DOF indices for both base and 'full' elements
1536 dof_map.dof_indices (child_base.get(), child_base_dof_indices,
1537 variable_number);
1538 dof_map.dof_indices (parent_base.get(), parent_base_dof_indices,
1539 variable_number);
1540
1541
1542 // First we loop over the childs base DOFs (nodes) and check which of them needs constraint
1543 // and which can be skipped.
1544 for (unsigned int child_base_dof=0; child_base_dof != n_base_dofs; ++child_base_dof)
1545 {
1546 libmesh_assert_less (child_base_dof, child_base->n_nodes());
1547
1548 // Childs global dof index.
1549 const dof_id_type child_base_dof_g = child_base_dof_indices[child_base_dof];
1550
1551 // Hunt for "constraining against myself" cases before
1552 // we bother creating a constraint row
1553 bool self_constraint = false;
1554 for (unsigned int parent_base_dof=0;
1555 parent_base_dof != n_base_dofs; parent_base_dof++)
1556 {
1557 libmesh_assert_less (parent_base_dof, parent_base->n_nodes());
1558
1559 // Their global dof index.
1560 const dof_id_type parent_base_dof_g =
1561 parent_base_dof_indices[parent_base_dof];
1562
1563 if (parent_base_dof_g == child_base_dof_g)
1564 {
1565 self_constraint = true;
1566 break;
1567 }
1568 }
1569
1570 if (self_constraint)
1571 continue;
1572
1573 // now we need to constrain all __child_elem__ DOFs whose base corresponds to
1574 // child_base_dof.
1575 // --> loop over all child_elem dofs whose base_shape_index == child_base_dof
1576 unsigned int n_elem_dofs = FEInterface::n_dofs(fe_type, child_elem);
1577 libmesh_assert_equal_to(n_elem_dofs, n_total_dofs);
1578 for(unsigned int child_elem_dof=0; child_elem_dof != n_elem_dofs; ++child_elem_dof)
1579 {
1580 if (base_shape_index[child_elem_dof] != child_base_dof)
1581 continue;
1582
1583 // independent from the radial description, the first radial DOF is 1 at the base
1584 // while all others start with 0.
1585 // Thus, to confine for the bases neighbor, we only need to refine DOFs that correspond
1586 // to the first radial DOF
1587 if (s==0)
1588 {
1589 if (radial_shape_index[child_elem_dof] > 0)
1590 continue;
1591 }
1592 else
1593 {
1594 // If the neighbor is not the base, we must check now if the support point of the dof
1595 // is actually shared with that neighbor:
1596 if ( !child_elem->neighbor_ptr(s)->contains_point(child_base->point(child_base_dof)) )
1597 continue;
1598 }
1599
1600
1601 const dof_id_type child_elem_dof_g = child_elem_dof_indices[child_elem_dof];
1602
1603 DofConstraintRow * constraint_row;
1604
1605 // we may be running constraint methods concurrently
1606 // on multiple threads, so we need a lock to
1607 // ensure that this constraint is "ours"
1608 {
1609 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1610
1611 if (dof_map.is_constrained_dof(child_elem_dof_g))
1612 continue;
1613
1614 constraint_row = &(constraints[child_elem_dof_g]);
1615 libmesh_assert(constraint_row->empty());
1616 }
1617
1618 // The support point of the DOF
1619 const Point & support_point = child_base->point(child_base_dof);
1620
1621 // Figure out where my (base) node lies on the parents reference element.
1622 const Point mapped_point = FEMap::inverse_map(Dim-1,
1623 parent_base.get(),
1624 support_point);
1625
1626 // now we need the parents base DOFs, evaluated at the mapped_point for refinement:
1627 for (unsigned int parent_base_dof=0;
1628 parent_base_dof != n_base_dofs; parent_base_dof++)
1629 {
1630
1631 const Real parent_base_dof_value = FEInterface::shape(fe_type,
1632 parent_base.get(),
1633 parent_base_dof,
1634 mapped_point);
1635
1636
1637 // all parent elements DOFs whose base_index corresponds to parent_base_dof
1638 // must be constrained with the parent_base_dof_value.
1639
1640 // The value of the radial function does not play a role here:
1641 // 1) only the function with radial_shape_index[] == 0 are 1 at the base,
1642 // the others are 0.
1643 // 2) The radial basis is (usually) not a Lagrange polynomial.
1644 // Thus, constraining according to a support point doesn't work.
1645 // However, they reach '1' at a certain (radial) distance which is the same for parent and child.
1646 for (unsigned int parent_elem_dof=0;
1647 parent_elem_dof != n_elem_dofs; parent_elem_dof++)
1648 {
1649 if (base_shape_index[parent_elem_dof] != parent_base_dof)
1650 continue;
1651
1652 // only constrain with coinciding radial DOFs.
1653 // Otherwise, we start coupling all DOFs with each other and end up in a mess.
1654 if (radial_shape_index[parent_elem_dof] != radial_shape_index[child_elem_dof])
1655 continue;
1656
1657 // Their global dof index.
1658 const dof_id_type parent_elem_dof_g =
1659 parent_elem_dof_indices[parent_elem_dof];
1660
1661 // Only add non-zero and non-identity values
1662 // for Lagrange basis functions. (parent_base is assumed to be of Lagrange-type).
1663 if ((std::abs(parent_base_dof_value) > 1.e-5) &&
1664 (std::abs(parent_base_dof_value) < .999))
1665 {
1666 constraint_row->emplace(parent_elem_dof_g, parent_base_dof_value);
1667 }
1668#ifdef DEBUG
1669 // Protect for the case u_i = 0.999 u_j,
1670 // in which case i better equal j.
1671 else if (parent_base_dof_value >= .999)
1672 {
1673 libmesh_assert_equal_to (child_base_dof_g, parent_base_dof_indices[parent_base_dof]);
1674 libmesh_assert_equal_to (child_elem_dof_g, parent_elem_dof_g);
1675 }
1676#endif
1677 }
1678
1679 }
1680 }
1681
1682 }
1683 }
1684}
1685
1686#endif // LIBMESH_ENABLE_AMR
1687
1688
1689//--------------------------------------------------------------
1690// Explicit instantiations
1691//#include "libmesh/inf_fe_instantiate_1D.h"
1692//#include "libmesh/inf_fe_instantiate_2D.h"
1693//#include "libmesh/inf_fe_instantiate_3D.h"
1694
1695#ifdef LIBMESH_ENABLE_DEPRECATED
1696INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const ElemType));
1697INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const ElemType));
1698INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const ElemType));
1699INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const ElemType, const unsigned int));
1700INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const ElemType, const unsigned int));
1701INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const ElemType, const unsigned int));
1702INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_shape_indices(const FEType &, const ElemType, const unsigned int, unsigned int &, unsigned int &));
1703INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_shape_indices(const FEType &, const ElemType, const unsigned int, unsigned int &, unsigned int &));
1704INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_shape_indices(const FEType &, const ElemType, const unsigned int, unsigned int &, unsigned int &));
1705INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_node_indices(const ElemType, const unsigned int, unsigned int &, unsigned int &));
1706INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_node_indices(const ElemType, const unsigned int, unsigned int &, unsigned int &));
1707INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_node_indices(const ElemType, const unsigned int, unsigned int &, unsigned int &));
1708INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
1709INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
1710INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
1711INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape_deriv(const FEType &, const ElemType, const unsigned int, const unsigned int, const Point &));
1712INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape_deriv(const FEType &, const ElemType, const unsigned int, const unsigned int, const Point &));
1713INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape_deriv(const FEType &, const ElemType, const unsigned int, const unsigned int, const Point &));
1714#endif // LIBMESH_ENABLE_DEPRECATED
1715
1716INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs(const FEType &, const Elem*));
1717INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs(const FEType &, const Elem*));
1718INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs(const FEType &, const Elem*));
1719INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const Elem *));
1720INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const Elem *));
1721INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const Elem *));
1722INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const Elem *, const unsigned int));
1723INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const Elem *, const unsigned int));
1724INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const Elem *, const unsigned int));
1725INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_shape_indices(const FEType &, const Elem *, const unsigned int, unsigned int &, unsigned int &));
1726INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_shape_indices(const FEType &, const Elem *, const unsigned int, unsigned int &, unsigned int &));
1727INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_shape_indices(const FEType &, const Elem *, const unsigned int, unsigned int &, unsigned int &));
1728INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape(const FEType &, const Elem *, const unsigned int, const Point & p));
1729INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape(const FEType &, const Elem *, const unsigned int, const Point & p));
1730INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape(const FEType &, const Elem *, const unsigned int, const Point & p));
1731INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape(const FEType, const Elem *, const unsigned int, const Point &, const bool));
1732INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape(const FEType, const Elem *, const unsigned int, const Point &, const bool));
1733INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape(const FEType, const Elem *, const unsigned int, const Point &, const bool));
1734INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape_deriv(const FEType &, const Elem *, const unsigned int, const unsigned int, const Point &));
1735INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape_deriv(const FEType &, const Elem *, const unsigned int, const unsigned int, const Point &));
1736INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape_deriv(const FEType &, const Elem *, const unsigned int, const unsigned int, const Point &));
1737INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape_deriv(const FEType, const Elem *, const unsigned int, const unsigned int, const Point &, const bool));
1738INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape_deriv(const FEType, const Elem *, const unsigned int, const unsigned int, const Point &, const bool));
1739INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape_deriv(const FEType, const Elem *, const unsigned int, const unsigned int, const Point &, const bool));
1740INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_data(const FEType &, const Elem *, FEComputeData &));
1741INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_data(const FEType &, const Elem *, FEComputeData &));
1742INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_data(const FEType &, const Elem *, FEComputeData &));
1743INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, nodal_soln(const FEType &, const Elem *, const std::vector<Number> &, std::vector<Number> &));
1744INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, nodal_soln(const FEType &, const Elem *, const std::vector<Number> &, std::vector<Number> &));
1745INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, nodal_soln(const FEType &, const Elem *, const std::vector<Number> &, std::vector<Number> &));
1746#ifdef LIBMESH_ENABLE_AMR
1747INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, inf_compute_constraints(DofConstraints &, DofMap &, const unsigned int, const Elem *));
1748INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, inf_compute_constraints(DofConstraints &, DofMap &, const unsigned int, const Elem *));
1749INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, inf_compute_constraints(DofConstraints &, DofMap &, const unsigned int, const Elem *));
1750#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1751INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, inf_compute_node_constraints(NodeConstraints & constraints,const Elem * elem));
1752INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, inf_compute_node_constraints(NodeConstraints & constraints,const Elem * elem));
1753INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, inf_compute_node_constraints(NodeConstraints & constraints,const Elem * elem));
1754#endif
1755#endif
1756
1757} // namespace libMesh
1758
1759#endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
libMesh::DofConstraints
The constraint matrix storage format.
Definition dof_map.h:112
libMesh::DofMap
This class handles the numbering of degrees of freedom on a mesh.
Definition dof_map.h:181
libMesh::DofMap::dof_indices
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition dof_map.C:2201
libMesh::DofMap::variable_type
const FEType & variable_type(const unsigned int i) const
Definition dof_map.h:2388
libMesh::DofMap::is_constrained_dof
bool is_constrained_dof(const dof_id_type dof) const
Definition dof_map.h:2426
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::neighbor_ptr_range
SimpleRange< NeighborPtrIter > neighbor_ptr_range()
Returns a range with all neighbors of an element, usable in range-based for loops.
Definition elem.h:3520
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition elem.h:2462
libMesh::Elem::origin
virtual Point origin() const
Definition elem.h:1923
libMesh::Elem::parent
const Elem * parent() const
Definition elem.h:3047
libMesh::Elem::contains_point
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const
Definition elem.C:2784
libMesh::Elem::level
unsigned int level() const
Definition elem.h:3091
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Elem::subactive
bool subactive() const
Definition elem.h:2976
libMesh::Elem::p_level
unsigned int p_level() const
Definition elem.h:3125
libMesh::Elem::neighbor_ptr
const Elem * neighbor_ptr(unsigned int i) const
Definition elem.h:2615
libMesh::Elem::side_index_range
IntRange< unsigned short > side_index_range() const
Definition elem.h:2727
libMesh::Elem::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i)=0
libMesh::FEComputeData
class FEComputeData hides arbitrary data to be passed to and from children of FEBase through the FEIn...
Definition fe_compute_data.h:52
libMesh::FEComputeData::p
const Point & p
Holds the point where the data are to be computed.
Definition fe_compute_data.h:76
libMesh::FEComputeData::dshape
std::vector< Gradient > dshape
Storage for the computed shape derivative values.
Definition fe_compute_data.h:86
libMesh::FEComputeData::need_derivative
bool need_derivative()
Check whether derivatives should be computed or not.
Definition fe_compute_data.h:143
libMesh::FEComputeData::local_transform
std::vector< std::vector< Real > > local_transform
Storage for local to global mapping at p.
Definition fe_compute_data.h:96
libMesh::FEComputeData::shape
std::vector< Number > shape
Storage for the computed shape function values.
Definition fe_compute_data.h:81
libMesh::FEComputeData::phase
Real phase
Storage for the computed phase lag.
Definition fe_compute_data.h:102
libMesh::FEComputeData::speed
Real speed
The wave speed.
Definition fe_compute_data.h:107
libMesh::FEComputeData::frequency
Number frequency
The frequency to evaluate shape functions including the wave number depending terms.
Definition fe_compute_data.h:114
libMesh::FEGenericBase::build_InfFE
static std::unique_ptr< FEGenericBase > build_InfFE(const unsigned int dim, const FEType &type)
Builds a specific infinite element type.
libMesh::FEInterface::shape
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition fe_interface.C:761
libMesh::FEInterface::shape_deriv
static Real shape_deriv(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const unsigned int j, const Point &p)
Definition fe_interface.C:1381
libMesh::FEInterface::n_dofs
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition fe_interface.C:353
libMesh::FEInterface::n_dofs_at_node
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
Definition fe_interface.C:437
libMesh::FEInterface::n_dofs_per_elem
static unsigned int n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition fe_interface.C:531
libMesh::FEMap::inverse_map
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
Definition fe_map.C:1512
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
libMesh::FEType::radial_order
OrderWrapper radial_order
The approximation order in radial direction of the infinite element.
Definition fe_type.h:263
libMesh::FEType::order
OrderWrapper order
The approximation order of the element (at 0 p-refinement level).
Definition fe_type.h:203
libMesh::FEType::family
FEFamily family
The type of finite element.
Definition fe_type.h:228
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::InfFEBase::get_elem_type
static ElemType get_elem_type(const ElemType type)
Definition inf_fe_base_radial.C:45
libMesh::InfFEBase::build_elem
static std::unique_ptr< const Elem > build_elem(const Elem *inf_elem)
Build the base element of an infinite element.
Definition inf_fe_base_radial.C:38
libMesh::InfFERadial::decay_deriv
static Real decay_deriv(const unsigned int dim, const Real)
Definition inf_fe.h:1292
libMesh::InfFERadial::n_dofs_per_elem
static unsigned int n_dofs_per_elem(const Order o_radial)
Definition inf_fe.h:136
libMesh::InfFERadial::n_dofs_at_node
static unsigned int n_dofs_at_node(const Order o_radial, const unsigned int n_onion)
Definition inf_fe_base_radial.C:115
libMesh::InfFERadial::mapping_order
static Order mapping_order()
Definition inf_fe.h:97
libMesh::InfFERadial::n_dofs
static unsigned int n_dofs(const Order o_radial)
Definition inf_fe.h:113
libMesh::InfFERadial::decay
static Real decay(const unsigned int dim, const Real v)
Definition inf_fe.h:1266
libMesh::InfFE
A specific instantiation of the FEBase class.
Definition inf_fe.h:225
libMesh::InfFE::n_dofs
static unsigned int n_dofs(const FEType &fet, const Elem *inf_elem)
Definition inf_fe_static.C:66
libMesh::InfFE::n_dofs_at_node
static unsigned int n_dofs_at_node(const FEType &fet, const ElemType inf_elem_type, const unsigned int n)
Definition inf_fe_static.C:84
libMesh::InfFE::eval
static Real eval(Real v, Order o_radial, unsigned int i)
libMesh::InfFE::_warned_for_nodal_soln
static bool _warned_for_nodal_soln
static members that are used to issue warning messages only once.
Definition inf_fe.h:1241
libMesh::InfFE::inf_compute_constraints
static void inf_compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *child_elem)
Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to var...
Definition inf_fe_static.C:1417
libMesh::InfFE::compute_node_indices_fast
static void compute_node_indices_fast(const ElemType inf_elem_type, const unsigned int outer_node_index, unsigned int &base_node, unsigned int &radial_node)
Does the same as compute_node_indices(), but stores the maps for the current element type.
Definition inf_fe_static.C:817
libMesh::InfFE::eval_deriv
static Real eval_deriv(Real v, Order o_radial, unsigned int i)
libMesh::InfFE::inf_compute_node_constraints
static void inf_compute_node_constraints(NodeConstraints &constraints, const Elem *elem)
Definition inf_fe_static.C:1373
libMesh::InfFE::_warned_for_shape
static bool _warned_for_shape
Definition inf_fe.h:1242
libMesh::InfFE::shape_deriv
static Real shape_deriv(const FEType &fet, const Elem *inf_elem, const unsigned int i, const unsigned int j, const Point &p)
Definition inf_fe_static.C:354
libMesh::InfFE::compute_shape_indices
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)
Definition inf_fe_static.C:1159
libMesh::InfFE::shape
static Real shape(const FEType &fet, const ElemType t, const unsigned int i, const Point &p)
Definition inf_fe_static.C:206
libMesh::InfFE::_warned_for_dshape
static bool _warned_for_dshape
Definition inf_fe.h:1243
libMesh::InfFE::compute_data
static void compute_data(const FEType &fe_t, const Elem *inf_elem, FEComputeData &data)
Generalized version of shape(), takes an Elem *.
Definition inf_fe_static.C:407
libMesh::InfFE::compute_node_indices
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...
Definition inf_fe_static.C:610
libMesh::InfFE::n_dofs_per_elem
static unsigned int n_dofs_per_elem(const FEType &fet, const ElemType inf_elem_type)
Definition inf_fe_static.C:135
libMesh::InfFE::_compute_node_indices_fast_current_elem_type
static ElemType _compute_node_indices_fast_current_elem_type
When compute_node_indices_fast() is used, this static variable remembers the element type for which t...
Definition inf_fe.h:1233
libMesh::InfFE::nodal_soln
static void nodal_soln(const FEType &fet, const Elem *elem, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
Usually, this method would build the nodal soln from the element soln.
Definition inf_fe_static.C:170
libMesh::NodeConstraints
The Node constraint storage format.
Definition dof_map.h:160
libMesh::OrderWrapper::get_order
int get_order() const
Explicitly request the order as an int.
Definition fe_type.h:80
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::TypeVector::norm
auto norm() const
Definition type_vector.h:908
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition threads.C:30
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::err
OStreamProxy err
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition libmesh.h:303
libMesh::pi
const Real pi
.
Definition libmesh.h:292
libMesh::remote_elem
const RemoteElem * remote_elem
Definition remote_elem.C:57
libMesh::DofConstraintRow
std::map< dof_id_type, Real, std::less< dof_id_type >, Threads::scalable_allocator< std::pair< const dof_id_type, Real > > > DofConstraintRow
A row of the Dof constraint matrix.
Definition dof_map.h:100
libMesh::imaginary
const Number imaginary
The imaginary unit, .
libMesh::dof_id_type
uint8_t dof_id_type
Definition id_types.h:67
libMesh::INSTANTIATE_INF_FE_MBRF
INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, reinit(const Elem *, const unsigned int, const Real, const std::vector< Point > *const, const std::vector< Real > *const))
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
n_nodes
const dof_id_type n_nodes
Definition tecplot_io.C:67
Generated on Fri Aug 21 2026 13:32:29 for libMesh by doxygen 1.9.8