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