Line data Source code
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 :
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 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
66 227358 : unsigned 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 153042 : auto base_elem = inf_elem->build_side_ptr(0);
72 :
73 : if (Dim > 1)
74 227358 : return FEInterface::n_dofs(fet, base_elem.get()) *
75 301674 : InfFERadial::n_dofs(fet.radial_order);
76 : else
77 0 : return InfFERadial::n_dofs(fet.radial_order);
78 74460 : }
79 :
80 :
81 :
82 : #ifdef LIBMESH_ENABLE_DEPRECATED
83 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
84 0 : unsigned 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 0 : const ElemType base_et (InfFEBase::get_elem_type(inf_elem_type));
91 :
92 : unsigned int n_base, n_radial;
93 0 : 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 0 : return FEInterface::n_dofs_at_node(Dim-1, fet, base_et, n_base)
104 0 : * InfFERadial::n_dofs_at_node(fet.radial_order, n_radial);
105 : else
106 0 : return InfFERadial::n_dofs_at_node(fet.radial_order, n_radial);
107 : }
108 : #endif // LIBMESH_ENABLE_DEPRECATED
109 :
110 :
111 :
112 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
113 258920 : unsigned 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 160640 : auto base_elem = inf_elem->build_side_ptr(0);
120 :
121 : unsigned int n_base, n_radial;
122 258920 : compute_node_indices(inf_elem->type(), n, n_base, n_radial);
123 :
124 : if (Dim > 1)
125 258920 : return FEInterface::n_dofs_at_node(fet, base_elem.get(), n_base)
126 258920 : * InfFERadial::n_dofs_at_node(fet.radial_order, n_radial);
127 : else
128 0 : return InfFERadial::n_dofs_at_node(fet.radial_order, n_radial);
129 62344 : }
130 :
131 :
132 :
133 : #ifdef LIBMESH_ENABLE_DEPRECATED
134 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
135 0 : unsigned 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 0 : const ElemType base_et (InfFEBase::get_elem_type(inf_elem_type));
141 :
142 : if (Dim > 1)
143 0 : return FEInterface::n_dofs_per_elem(Dim-1, fet, base_et)
144 0 : * InfFERadial::n_dofs_per_elem(fet.radial_order);
145 : else
146 0 : return InfFERadial::n_dofs_per_elem(fet.radial_order);
147 : }
148 : #endif // LIBMESH_ENABLE_DEPRECATED
149 :
150 :
151 :
152 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
153 22796 : unsigned 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 14256 : auto base_elem = inf_elem->build_side_ptr(0);
159 :
160 : if (Dim > 1)
161 22796 : return FEInterface::n_dofs_per_elem(fet, base_elem.get())
162 31336 : * InfFERadial::n_dofs_per_elem(fet.radial_order);
163 : else
164 0 : return InfFERadial::n_dofs_per_elem(fet.radial_order);
165 5714 : }
166 :
167 :
168 :
169 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
170 0 : void 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 0 : if (!_warned_for_nodal_soln)
177 : {
178 0 : libMesh::err << "WARNING: nodal_soln(...) does _not_ work for infinite elements." << std::endl
179 0 : << " Will return an empty nodal solution. Use " << std::endl
180 0 : << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!" << std::endl;
181 0 : _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 0 : nodal_soln.clear();
193 0 : libmesh_assert (nodal_soln.empty());
194 0 : return;
195 : }
196 :
197 :
198 :
199 :
200 :
201 :
202 :
203 :
204 : #ifdef LIBMESH_ENABLE_DEPRECATED
205 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
206 0 : Real 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 0 : if ((T_radial != INFINITE_MAP) && !_warned_for_shape)
218 : {
219 0 : libMesh::err << "WARNING: InfFE<Dim,T_radial,T_map>::shape(...) does _not_" << std::endl
220 0 : << " return the correct trial function! Use " << std::endl
221 0 : << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!"
222 0 : << std::endl;
223 0 : _warned_for_shape = true;
224 : }
225 : #endif
226 :
227 0 : const ElemType base_et (InfFEBase::get_elem_type(inf_elem_type));
228 0 : const Order o_radial (fet.radial_order);
229 0 : const Real v (p(Dim-1));
230 :
231 : unsigned int i_base, i_radial;
232 0 : 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 0 : return FEInterface::shape(Dim-1, fet, base_et, i_base, p)
239 0 : * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
240 0 : * InfFERadial::decay(Dim,v);
241 : else
242 0 : return InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
243 0 : * InfFERadial::decay(Dim,v);
244 : }
245 : #endif // LIBMESH_ENABLE_DEPRECATED
246 :
247 :
248 :
249 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
250 0 : Real 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 0 : 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 0 : if ((T_radial != INFINITE_MAP) && !_warned_for_shape)
261 : {
262 0 : libMesh::err << "WARNING: InfFE<Dim,T_radial,T_map>::shape(...) does _not_" << std::endl
263 0 : << " return the correct trial function! Use " << std::endl
264 0 : << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!"
265 0 : << std::endl;
266 0 : _warned_for_shape = true;
267 : }
268 : #endif
269 :
270 0 : const Order o_radial (fet.radial_order);
271 0 : const Real v (p(Dim-1));
272 0 : std::unique_ptr<const Elem> base_el (inf_elem->build_side_ptr(0));
273 :
274 : unsigned int i_base, i_radial;
275 0 : compute_shape_indices(fet, inf_elem, i, i_base, i_radial);
276 :
277 : if (Dim > 1)
278 0 : return FEInterface::shape(fet, base_el.get(), i_base, p)
279 0 : * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
280 0 : * InfFERadial::decay(Dim,v);
281 : else
282 0 : return InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
283 0 : * InfFERadial::decay(Dim,v);
284 0 : }
285 :
286 :
287 :
288 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
289 0 : Real 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 0 : if (add_p_level)
296 : {
297 0 : FEType tmp_fet=fet;
298 0 : tmp_fet = fet.order + add_p_level * inf_elem->p_level();
299 0 : return InfFE<Dim,T_radial,T_map>::shape(tmp_fet, inf_elem, i, p);
300 : }
301 0 : return InfFE<Dim,T_radial,T_map>::shape(fet, inf_elem, i, p);
302 : }
303 :
304 :
305 : #ifdef LIBMESH_ENABLE_DEPRECATED
306 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
307 0 : Real 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 0 : libmesh_assert_greater (Dim,j);
317 : #ifdef DEBUG
318 : // this makes only sense when used for mapping
319 0 : if ((T_radial != INFINITE_MAP) && !_warned_for_dshape)
320 : {
321 0 : libMesh::err << "WARNING: InfFE<Dim,T_radial,T_map>::shape_deriv(...) does _not_" << std::endl
322 0 : << " return the correct trial function gradients! Use " << std::endl
323 0 : << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!"
324 0 : << std::endl;
325 0 : _warned_for_dshape = true;
326 : }
327 : #endif
328 :
329 0 : const ElemType base_et (InfFEBase::get_elem_type(inf_elem_type));
330 0 : const Order o_radial (fe_t.radial_order);
331 0 : const Real v (p(Dim-1));
332 :
333 : unsigned int i_base, i_radial;
334 0 : compute_shape_indices(fe_t, inf_elem_type, i, i_base, i_radial);
335 :
336 0 : if (j== Dim -1)
337 : {
338 0 : Real RadialDeriv = InfFE<Dim,T_radial,T_map>::eval_deriv(v, o_radial, i_radial)
339 0 : * InfFERadial::decay(Dim,v)
340 0 : + InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
341 0 : * InfFERadial::decay_deriv(Dim, v);
342 :
343 0 : return FEInterface::shape(Dim-1, fe_t, base_et, i_base, p)*RadialDeriv;
344 : }
345 :
346 0 : return FEInterface::shape_deriv(Dim-1, fe_t, base_et, i_base, j, p)
347 0 : * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
348 0 : * InfFERadial::decay(Dim,v);
349 : }
350 : #endif // LIBMESH_ENABLE_DEPRECATED
351 :
352 :
353 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
354 0 : Real 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 0 : libmesh_assert_greater (Dim,j);
362 : #ifdef DEBUG
363 : // this makes only sense when used for mapping
364 0 : if ((T_radial != INFINITE_MAP) && !_warned_for_dshape)
365 : {
366 0 : libMesh::err << "WARNING: InfFE<Dim,T_radial,T_map>::shape_deriv(...) does _not_" << std::endl
367 0 : << " return the correct trial function gradients! Use " << std::endl
368 0 : << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!"
369 0 : << std::endl;
370 0 : _warned_for_dshape = true;
371 : }
372 : #endif
373 0 : const Order o_radial (fe_t.radial_order);
374 0 : const Real v (p(Dim-1));
375 :
376 0 : std::unique_ptr<const Elem> base_el (inf_elem->build_side_ptr(0));
377 :
378 : unsigned int i_base, i_radial;
379 :
380 0 : 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 0 : i_base=0;
385 : }
386 0 : compute_shape_indices(fe_t, inf_elem, i, i_base, i_radial);
387 :
388 0 : if (j== Dim -1)
389 : {
390 0 : Real RadialDeriv = InfFE<Dim,T_radial,T_map>::eval_deriv(v, o_radial, i_radial)
391 0 : * InfFERadial::decay(Dim,v)
392 0 : + InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
393 0 : * InfFERadial::decay_deriv(Dim,v);
394 :
395 0 : return FEInterface::shape(fe_t, base_el.get(), i_base, p)*RadialDeriv;
396 : }
397 0 : return FEInterface::shape_deriv(fe_t, base_el.get(), i_base, j, p)
398 0 : * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
399 0 : * InfFERadial::decay(Dim,v);
400 0 : }
401 :
402 :
403 :
404 :
405 :
406 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
407 0 : void InfFE<Dim,T_radial,T_map>::compute_data(const FEType & fet,
408 : const Elem * inf_elem,
409 : FEComputeData & data)
410 : {
411 0 : libmesh_assert(inf_elem);
412 : libmesh_assert_not_equal_to (Dim, 0);
413 :
414 0 : const Order o_radial (fet.radial_order);
415 0 : const Order radial_mapping_order (InfFERadial::mapping_order());
416 0 : const Point & p (data.p);
417 0 : const Real v (p(Dim-1));
418 0 : 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 0 : Real interpolated_dist = 0.;
427 : switch (Dim)
428 : {
429 : case 1:
430 : {
431 0 : libmesh_assert_equal_to (inf_elem->type(), INFEDGE2);
432 0 : interpolated_dist = Point(inf_elem->point(0) - inf_elem->point(1)).norm();
433 0 : break;
434 : }
435 :
436 : case 2:
437 : {
438 0 : const unsigned int n_base_nodes = base_el->n_nodes();
439 :
440 0 : const Point origin = inf_elem->origin();
441 0 : const Order base_mapping_order (base_el->default_order());
442 0 : const ElemType base_mapping_elem_type (base_el->type());
443 :
444 : // interpolate the base nodes' distances
445 0 : for (unsigned int n=0; n<n_base_nodes; n++)
446 0 : interpolated_dist += Point(base_el->point(n) - origin).norm()
447 0 : * FE<1,LAGRANGE>::shape (base_mapping_elem_type, base_mapping_order, n, p);
448 0 : break;
449 : }
450 :
451 : case 3:
452 : {
453 0 : const unsigned int n_base_nodes = base_el->n_nodes();
454 :
455 0 : const Point origin = inf_elem->origin();
456 0 : const Order base_mapping_order (base_el->default_order());
457 0 : const ElemType base_mapping_elem_type (base_el->type());
458 :
459 : // interpolate the base nodes' distances
460 0 : for (unsigned int n=0; n<n_base_nodes; n++)
461 0 : interpolated_dist += Point(base_el->point(n) - origin).norm()
462 0 : * FE<2,LAGRANGE>::shape (base_mapping_elem_type, base_mapping_order, n, p);
463 0 : break;
464 : }
465 :
466 : default:
467 : libmesh_error_msg("Unknown Dim = " << Dim);
468 : }
469 :
470 :
471 0 : 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 0 : data.phase = interpolated_dist /* together with next line: */
476 0 : * 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 0 : 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 0 : const unsigned int n_dof = n_dofs (fet, inf_elem);
502 0 : data.shape.resize(n_dof);
503 0 : if (data.need_derivative())
504 : {
505 0 : data.dshape.resize(n_dof);
506 0 : data.local_transform.resize(Dim);
507 :
508 0 : for (unsigned int d=0; d<Dim; d++)
509 0 : 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 0 : auto fe = FEBase::build_InfFE(Dim, fet);
515 0 : std::vector<Point> pt = {p};
516 0 : fe->get_dxidx(); // to compute the map
517 0 : fe->reinit(inf_elem, &pt);
518 :
519 : // compute the reference->physical map.
520 0 : data.local_transform[0][0] = fe->get_dxidx()[0];
521 0 : data.local_transform[1][0] = fe->get_detadx()[0];
522 0 : data.local_transform[1][1] = fe->get_detady()[0];
523 0 : data.local_transform[0][1] = fe->get_dxidy()[0];
524 : if (Dim > 2)
525 : {
526 0 : data.local_transform[2][0] = fe->get_dzetadx()[0];
527 0 : data.local_transform[2][1] = fe->get_dzetady()[0];
528 0 : data.local_transform[2][2] = fe->get_dzetadz()[0];
529 0 : data.local_transform[1][2] = fe->get_detadz()[0];
530 0 : data.local_transform[0][2] = fe->get_dxidz()[0];
531 : }
532 0 : } // endif data.need_derivative()
533 :
534 0 : 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 0 : compute_shape_indices(fet, inf_elem, i, i_base, i_radial);
539 :
540 0 : data.shape[i] = (InfFERadial::decay(Dim,v) /* (1.-v)/2. in 3D */
541 0 : * FEInterface::shape(fet, base_el.get(), i_base, p) /* S_n(s,t) */
542 0 : * 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 0 : if (data.need_derivative())
547 : {
548 0 : data.dshape[i](0) = (InfFERadial::decay(Dim,v)
549 0 : * FEInterface::shape_deriv(fet, base_el.get(), i_base, 0, p)
550 0 : * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial))
551 : * time_harmonic;
552 :
553 : if (Dim > 2)
554 : {
555 0 : data.dshape[i](1) = (InfFERadial::decay(Dim,v)
556 0 : * FEInterface::shape_deriv(fet, base_el.get(), i_base, 1, p)
557 0 : * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial))
558 : * time_harmonic;
559 :
560 : }
561 0 : data.dshape[i](Dim-1) = (InfFERadial::decay_deriv(Dim, v) * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
562 0 : +InfFERadial::decay(Dim,v) * InfFE<Dim,T_radial,T_map>::eval_deriv(v, o_radial, i_radial))
563 0 : * 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 0 : *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 0 : 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 0 : libmesh_error_msg("compute_data() for 1-dimensional InfFE not implemented.");
584 0 : }
585 :
586 :
587 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
588 0 : Real 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 0 : if (add_p_level)
596 : {
597 0 : FEType tmp_fet=fet;
598 0 : tmp_fet = fet.order + add_p_level * inf_elem->p_level();
599 0 : return InfFE<Dim,T_radial,T_map>::shape_deriv(tmp_fet, inf_elem, i, j, p);
600 : }
601 0 : return InfFE<Dim,T_radial,T_map>::shape_deriv(fet, inf_elem, i, j, p);
602 : }
603 :
604 :
605 :
606 :
607 :
608 :
609 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
610 258920 : void 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 258920 : switch (inf_elem_type)
616 : {
617 0 : case INFEDGE2:
618 : {
619 0 : libmesh_assert_less (outer_node_index, 2);
620 0 : base_node = 0;
621 0 : radial_node = outer_node_index;
622 0 : return;
623 : }
624 :
625 :
626 : // linear base approximation, easy to determine
627 0 : case INFQUAD4:
628 : {
629 0 : libmesh_assert_less (outer_node_index, 4);
630 0 : base_node = outer_node_index % 2;
631 0 : radial_node = outer_node_index / 2;
632 0 : return;
633 : }
634 :
635 0 : case INFPRISM6:
636 : {
637 0 : libmesh_assert_less (outer_node_index, 6);
638 0 : base_node = outer_node_index % 3;
639 0 : radial_node = outer_node_index / 3;
640 0 : return;
641 : }
642 :
643 110440 : case INFHEX8:
644 : {
645 42452 : libmesh_assert_less (outer_node_index, 8);
646 110440 : base_node = outer_node_index % 4;
647 110440 : radial_node = outer_node_index / 4;
648 110440 : return;
649 : }
650 :
651 :
652 : // higher order base approximation, more work necessary
653 0 : case INFQUAD6:
654 : {
655 0 : switch (outer_node_index)
656 : {
657 0 : case 0:
658 : case 1:
659 : {
660 0 : radial_node = 0;
661 0 : base_node = outer_node_index;
662 0 : return;
663 : }
664 :
665 0 : case 2:
666 : case 3:
667 : {
668 0 : radial_node = 1;
669 0 : base_node = outer_node_index-2;
670 0 : return;
671 : }
672 :
673 0 : case 4:
674 : {
675 0 : radial_node = 0;
676 0 : base_node = 2;
677 0 : return;
678 : }
679 :
680 0 : case 5:
681 : {
682 0 : radial_node = 1;
683 0 : base_node = 2;
684 0 : return;
685 : }
686 :
687 0 : default:
688 0 : libmesh_error_msg("Unrecognized outer_node_index = " << outer_node_index);
689 : }
690 : }
691 :
692 :
693 75260 : case INFHEX16:
694 : case INFHEX18:
695 : {
696 75260 : switch (outer_node_index)
697 : {
698 18416 : case 0:
699 : case 1:
700 : case 2:
701 : case 3:
702 : {
703 18416 : radial_node = 0;
704 18416 : base_node = outer_node_index;
705 18416 : return;
706 : }
707 :
708 18704 : case 4:
709 : case 5:
710 : case 6:
711 : case 7:
712 : {
713 18704 : radial_node = 1;
714 18704 : base_node = outer_node_index-4;
715 18704 : return;
716 : }
717 :
718 14826 : case 8:
719 : case 9:
720 : case 10:
721 : case 11:
722 : {
723 14826 : radial_node = 0;
724 14826 : base_node = outer_node_index-4;
725 14826 : return;
726 : }
727 :
728 16092 : case 12:
729 : case 13:
730 : case 14:
731 : case 15:
732 : {
733 16092 : radial_node = 1;
734 16092 : base_node = outer_node_index-8;
735 16092 : return;
736 : }
737 :
738 3716 : case 16:
739 : {
740 1283 : libmesh_assert_equal_to (inf_elem_type, INFHEX18);
741 3716 : radial_node = 0;
742 3716 : base_node = 8;
743 3716 : return;
744 : }
745 :
746 3506 : case 17:
747 : {
748 1203 : libmesh_assert_equal_to (inf_elem_type, INFHEX18);
749 3506 : radial_node = 1;
750 3506 : base_node = 8;
751 3506 : return;
752 : }
753 :
754 0 : default:
755 0 : libmesh_error_msg("Unrecognized outer_node_index = " << outer_node_index);
756 : }
757 : }
758 :
759 :
760 73220 : case INFPRISM12:
761 : {
762 73220 : switch (outer_node_index)
763 : {
764 20148 : case 0:
765 : case 1:
766 : case 2:
767 : {
768 20148 : radial_node = 0;
769 20148 : base_node = outer_node_index;
770 20148 : return;
771 : }
772 :
773 20120 : case 3:
774 : case 4:
775 : case 5:
776 : {
777 20120 : radial_node = 1;
778 20120 : base_node = outer_node_index-3;
779 20120 : return;
780 : }
781 :
782 16172 : case 6:
783 : case 7:
784 : case 8:
785 : {
786 16172 : radial_node = 0;
787 16172 : base_node = outer_node_index-3;
788 16172 : return;
789 : }
790 :
791 16780 : case 9:
792 : case 10:
793 : case 11:
794 : {
795 16780 : radial_node = 1;
796 16780 : base_node = outer_node_index-6;
797 16780 : return;
798 : }
799 :
800 0 : default:
801 0 : libmesh_error_msg("Unrecognized outer_node_index = " << outer_node_index);
802 : }
803 : }
804 :
805 :
806 0 : default:
807 0 : 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 :
816 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
817 : void 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 :
925 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
926 56871 : void 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 56871 : const unsigned int radial_order = static_cast<unsigned int>(fet.radial_order.get_order()); // 4
956 56871 : const unsigned int radial_order_p_one = radial_order+1; // 5
957 :
958 70391 : 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 13520 : unsigned int n_base_vertices = libMesh::invalid_uint; // 4
962 56871 : const unsigned int n_base_vertex_dof = FEInterface::n_dofs_at_node (fet, base_elem.get(), 0); // 2
963 :
964 13520 : unsigned int n_base_side_nodes = libMesh::invalid_uint; // 4
965 13520 : unsigned int n_base_side_dof = libMesh::invalid_uint; // 3
966 :
967 13520 : unsigned int n_base_face_nodes = libMesh::invalid_uint; // 1
968 13520 : unsigned int n_base_face_dof = libMesh::invalid_uint; // 5
969 :
970 56871 : const unsigned int n_base_elem_dof = FEInterface::n_dofs_per_elem (fet, base_elem.get()); // 9
971 :
972 56871 : const ElemType inf_elem_type = inf_elem->type();
973 :
974 56871 : switch (inf_elem_type)
975 : {
976 0 : case INFEDGE2:
977 : {
978 0 : n_base_vertices = 1;
979 0 : n_base_side_nodes = 0;
980 0 : n_base_face_nodes = 0;
981 0 : n_base_side_dof = 0;
982 0 : n_base_face_dof = 0;
983 0 : break;
984 : }
985 :
986 0 : case INFQUAD4:
987 : {
988 0 : n_base_vertices = 2;
989 0 : n_base_side_nodes = 0;
990 0 : n_base_face_nodes = 0;
991 0 : n_base_side_dof = 0;
992 0 : n_base_face_dof = 0;
993 0 : break;
994 : }
995 :
996 0 : case INFQUAD6:
997 : {
998 0 : n_base_vertices = 2;
999 0 : n_base_side_nodes = 1;
1000 0 : n_base_face_nodes = 0;
1001 0 : n_base_side_dof = FEInterface::n_dofs_at_node (fet, base_elem.get(), n_base_vertices);
1002 0 : n_base_face_dof = 0;
1003 0 : break;
1004 : }
1005 :
1006 15736 : case INFHEX8:
1007 : {
1008 5472 : n_base_vertices = 4;
1009 5472 : n_base_side_nodes = 0;
1010 5472 : n_base_face_nodes = 0;
1011 5472 : n_base_side_dof = 0;
1012 5472 : n_base_face_dof = 0;
1013 15736 : break;
1014 : }
1015 :
1016 0 : case INFHEX16:
1017 : {
1018 0 : n_base_vertices = 4;
1019 0 : n_base_side_nodes = 4;
1020 0 : n_base_face_nodes = 0;
1021 0 : n_base_side_dof = FEInterface::n_dofs_at_node (fet, base_elem.get(), n_base_vertices);
1022 0 : n_base_face_dof = 0;
1023 0 : break;
1024 : }
1025 :
1026 31205 : case INFHEX18:
1027 : {
1028 4712 : n_base_vertices = 4;
1029 4712 : n_base_side_nodes = 4;
1030 4712 : n_base_face_nodes = 1;
1031 31205 : n_base_side_dof = FEInterface::n_dofs_at_node (fet, base_elem.get(), n_base_vertices);
1032 31205 : n_base_face_dof = FEInterface::n_dofs_at_node (fet, base_elem.get(), 8);
1033 4712 : break;
1034 : }
1035 :
1036 :
1037 0 : case INFPRISM6:
1038 : {
1039 0 : n_base_vertices = 3;
1040 0 : n_base_side_nodes = 0;
1041 0 : n_base_face_nodes = 0;
1042 0 : n_base_side_dof = 0;
1043 0 : n_base_face_dof = 0;
1044 0 : break;
1045 : }
1046 :
1047 9930 : case INFPRISM12:
1048 : {
1049 3336 : n_base_vertices = 3;
1050 3336 : n_base_side_nodes = 3;
1051 3336 : n_base_face_nodes = 0;
1052 9930 : n_base_side_dof = FEInterface::n_dofs_at_node (fet, base_elem.get(), n_base_vertices);
1053 3336 : n_base_face_dof = 0;
1054 3336 : break;
1055 : }
1056 :
1057 0 : default:
1058 0 : 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 56871 : const unsigned int n_dof_at_base_vertices = n_base_vertices*n_base_vertex_dof; // 8
1065 56871 : const unsigned int n_dof_at_all_vertices = n_dof_at_base_vertices*radial_order_p_one; // 40
1066 :
1067 56871 : const unsigned int n_dof_at_base_sides = n_base_side_nodes*n_base_side_dof; // 12
1068 56871 : const unsigned int n_dof_at_all_sides = n_dof_at_base_sides*radial_order_p_one; // 60
1069 :
1070 56871 : const unsigned int n_dof_at_base_face = n_base_face_nodes*n_base_face_dof; // 5
1071 56871 : 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 56871 : if (i < n_dof_at_base_vertices) // range of i: 0..7
1076 : {
1077 : // belongs to vertex in the base
1078 7974 : radial_shape = 0;
1079 7974 : base_shape = i;
1080 : }
1081 :
1082 48897 : 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 36122 : 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 36122 : radial_shape = (i_offset % radial_order) + 1;
1093 36122 : base_shape = i_offset / radial_order;
1094 : }
1095 :
1096 12775 : 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 2251 : radial_shape = 0;
1100 2251 : base_shape = i - radial_order * n_dof_at_base_vertices; // 8..19
1101 : }
1102 :
1103 10524 : 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 11949 : const unsigned int i_offset = i - (n_dof_at_all_vertices
1107 2990 : + n_dof_at_base_sides); // 0..47
1108 8959 : radial_shape = (i_offset % radial_order) + 1;
1109 8959 : base_shape = (i_offset / radial_order) + n_dof_at_base_vertices;
1110 : }
1111 :
1112 1565 : 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 313 : radial_shape = 0;
1116 417 : base_shape = i - radial_order*(n_dof_at_base_vertices
1117 313 : + n_dof_at_base_sides); // 20..24
1118 : }
1119 :
1120 1252 : 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 1668 : const unsigned int i_offset = i - (n_dof_at_all_vertices
1124 416 : + n_dof_at_all_sides
1125 416 : + n_dof_at_base_face); // 0..19
1126 1252 : radial_shape = (i_offset % radial_order) + 1;
1127 1252 : base_shape = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides;
1128 : }
1129 :
1130 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
1131 : {
1132 : // belongs to the base and is an element associated shape
1133 0 : radial_shape = 0;
1134 0 : base_shape = i - (n_dof_at_all_vertices
1135 0 : + n_dof_at_all_sides
1136 0 : + n_dof_at_all_faces); // 0..8
1137 : }
1138 :
1139 : else // range of i: 134..169
1140 : {
1141 0 : libmesh_assert_less (i, n_dofs(fet, inf_elem));
1142 : // belongs to the outer shell and is an element associated shape
1143 0 : const unsigned int i_offset = i - (n_dof_at_all_vertices
1144 0 : + n_dof_at_all_sides
1145 0 : + n_dof_at_all_faces
1146 0 : + n_base_elem_dof); // 0..19
1147 0 : radial_shape = (i_offset % radial_order) + 1;
1148 0 : 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 70391 : return;
1153 29831 : }
1154 :
1155 :
1156 :
1157 : #ifdef LIBMESH_ENABLE_DEPRECATED
1158 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
1159 0 : void 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 0 : const unsigned int radial_order = static_cast<unsigned int>(fet.radial_order.get_order()); // 4
1170 0 : const unsigned int radial_order_p_one = radial_order+1; // 5
1171 :
1172 0 : 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 0 : unsigned int n_base_vertices = libMesh::invalid_uint; // 4
1176 0 : const unsigned int n_base_vertex_dof = FEInterface::n_dofs_at_node (Dim-1, fet, base_elem_type, 0);// 2
1177 :
1178 0 : unsigned int n_base_side_nodes = libMesh::invalid_uint; // 4
1179 0 : unsigned int n_base_side_dof = libMesh::invalid_uint; // 3
1180 :
1181 0 : unsigned int n_base_face_nodes = libMesh::invalid_uint; // 1
1182 0 : unsigned int n_base_face_dof = libMesh::invalid_uint; // 5
1183 :
1184 0 : const unsigned int n_base_elem_dof = FEInterface::n_dofs_per_elem (Dim-1, fet, base_elem_type);// 9
1185 :
1186 :
1187 0 : switch (inf_elem_type)
1188 : {
1189 0 : case INFEDGE2:
1190 : {
1191 0 : n_base_vertices = 1;
1192 0 : n_base_side_nodes = 0;
1193 0 : n_base_face_nodes = 0;
1194 0 : n_base_side_dof = 0;
1195 0 : n_base_face_dof = 0;
1196 0 : break;
1197 : }
1198 :
1199 0 : case INFQUAD4:
1200 : {
1201 0 : n_base_vertices = 2;
1202 0 : n_base_side_nodes = 0;
1203 0 : n_base_face_nodes = 0;
1204 0 : n_base_side_dof = 0;
1205 0 : n_base_face_dof = 0;
1206 0 : break;
1207 : }
1208 :
1209 0 : case INFQUAD6:
1210 : {
1211 0 : n_base_vertices = 2;
1212 0 : n_base_side_nodes = 1;
1213 0 : n_base_face_nodes = 0;
1214 0 : n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
1215 0 : n_base_face_dof = 0;
1216 0 : break;
1217 : }
1218 :
1219 0 : case INFHEX8:
1220 : {
1221 0 : n_base_vertices = 4;
1222 0 : n_base_side_nodes = 0;
1223 0 : n_base_face_nodes = 0;
1224 0 : n_base_side_dof = 0;
1225 0 : n_base_face_dof = 0;
1226 0 : break;
1227 : }
1228 :
1229 0 : case INFHEX16:
1230 : {
1231 0 : n_base_vertices = 4;
1232 0 : n_base_side_nodes = 4;
1233 0 : n_base_face_nodes = 0;
1234 0 : n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
1235 0 : n_base_face_dof = 0;
1236 0 : break;
1237 : }
1238 :
1239 0 : case INFHEX18:
1240 : {
1241 0 : n_base_vertices = 4;
1242 0 : n_base_side_nodes = 4;
1243 0 : n_base_face_nodes = 1;
1244 0 : n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
1245 0 : n_base_face_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, 8);
1246 0 : break;
1247 : }
1248 :
1249 :
1250 0 : case INFPRISM6:
1251 : {
1252 0 : n_base_vertices = 3;
1253 0 : n_base_side_nodes = 0;
1254 0 : n_base_face_nodes = 0;
1255 0 : n_base_side_dof = 0;
1256 0 : n_base_face_dof = 0;
1257 0 : break;
1258 : }
1259 :
1260 0 : case INFPRISM12:
1261 : {
1262 0 : n_base_vertices = 3;
1263 0 : n_base_side_nodes = 3;
1264 0 : n_base_face_nodes = 0;
1265 0 : n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
1266 0 : n_base_face_dof = 0;
1267 0 : break;
1268 : }
1269 :
1270 0 : default:
1271 0 : 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 0 : const unsigned int n_dof_at_base_vertices = n_base_vertices*n_base_vertex_dof; // 8
1278 0 : const unsigned int n_dof_at_all_vertices = n_dof_at_base_vertices*radial_order_p_one; // 40
1279 :
1280 0 : const unsigned int n_dof_at_base_sides = n_base_side_nodes*n_base_side_dof; // 12
1281 0 : const unsigned int n_dof_at_all_sides = n_dof_at_base_sides*radial_order_p_one; // 60
1282 :
1283 0 : const unsigned int n_dof_at_base_face = n_base_face_nodes*n_base_face_dof; // 5
1284 0 : 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 0 : if (i < n_dof_at_base_vertices) // range of i: 0..7
1289 : {
1290 : // belongs to vertex in the base
1291 0 : radial_shape = 0;
1292 0 : base_shape = i;
1293 : }
1294 :
1295 0 : 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 0 : 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 0 : radial_shape = (i_offset % radial_order) + 1;
1306 0 : base_shape = i_offset / radial_order;
1307 : }
1308 :
1309 0 : 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 0 : radial_shape = 0;
1313 0 : base_shape = i - radial_order * n_dof_at_base_vertices; // 8..19
1314 : }
1315 :
1316 0 : 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 0 : const unsigned int i_offset = i - (n_dof_at_all_vertices
1320 0 : + n_dof_at_base_sides); // 0..47
1321 0 : radial_shape = (i_offset % radial_order) + 1;
1322 0 : base_shape = (i_offset / radial_order) + n_dof_at_base_vertices;
1323 : }
1324 :
1325 0 : 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 0 : radial_shape = 0;
1329 0 : base_shape = i - radial_order*(n_dof_at_base_vertices
1330 0 : + n_dof_at_base_sides); // 20..24
1331 : }
1332 :
1333 0 : 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 0 : const unsigned int i_offset = i - (n_dof_at_all_vertices
1337 0 : + n_dof_at_all_sides
1338 0 : + n_dof_at_base_face); // 0..19
1339 0 : radial_shape = (i_offset % radial_order) + 1;
1340 0 : base_shape = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides;
1341 : }
1342 :
1343 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
1344 : {
1345 : // belongs to the base and is an element associated shape
1346 0 : radial_shape = 0;
1347 0 : base_shape = i - (n_dof_at_all_vertices
1348 0 : + n_dof_at_all_sides
1349 0 : + 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 0 : const unsigned int i_offset = i - (n_dof_at_all_vertices
1356 0 : + n_dof_at_all_sides
1357 0 : + n_dof_at_all_faces
1358 0 : + n_base_elem_dof); // 0..19
1359 0 : radial_shape = (i_offset % radial_order) + 1;
1360 0 : 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 0 : return;
1365 : }
1366 : #endif // LIBMESH_ENABLE_DEPRECATED
1367 :
1368 :
1369 : #ifdef LIBMESH_ENABLE_AMR
1370 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1371 :
1372 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
1373 2112 : void 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 0 : return;
1378 :
1379 1056 : libmesh_assert(elem);
1380 :
1381 : // only constrain active and ancestor elements
1382 2112 : if (elem->subactive())
1383 0 : 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 12072 : for (auto s : elem->side_index_range())
1403 : {
1404 19920 : if (elem->neighbor_ptr(s) != nullptr &&
1405 9960 : elem->neighbor_ptr(s) != remote_elem)
1406 9960 : if (elem->neighbor_ptr(s)->level() < elem->level())
1407 : {
1408 0 : libmesh_not_implemented();
1409 : }
1410 : }
1411 : }
1412 :
1413 : #endif //LIBMESH_ENABLE_NODE_CONSTRAINTS
1414 :
1415 :
1416 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
1417 3212 : void 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 3186 : return;
1426 :
1427 1056 : libmesh_assert(child_elem);
1428 :
1429 : // only constrain active and ancestor elements
1430 3212 : if (child_elem->subactive())
1431 0 : return;
1432 :
1433 : // Before we start to compute anything, lets check if any confinement is needed:
1434 1056 : bool need_constraints=false;
1435 18284 : for (auto child_neighbor : child_elem->neighbor_ptr_range())
1436 15098 : if (child_neighbor->level() < child_elem->level())
1437 : {
1438 0 : need_constraints = true;
1439 0 : break;
1440 : }
1441 3212 : if (!need_constraints)
1442 1056 : 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 26 : FEType fe_type = dof_map.variable_type(variable_number);
1472 :
1473 0 : libmesh_assert(fe_type.family == LAGRANGE);
1474 :
1475 0 : std::vector<dof_id_type> child_base_dof_indices, parent_base_dof_indices;
1476 0 : std::vector<dof_id_type> child_elem_dof_indices, parent_elem_dof_indices;
1477 :
1478 0 : 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 0 : libmesh_assert(parent_elem);
1484 :
1485 26 : dof_map.dof_indices (child_elem, child_elem_dof_indices,
1486 : variable_number);
1487 26 : dof_map.dof_indices (parent_elem, parent_elem_dof_indices,
1488 : variable_number);
1489 :
1490 26 : 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 26 : std::vector<unsigned int> radial_shape_index(n_total_dofs);
1494 26 : std::vector<unsigned int> base_shape_index(n_total_dofs);
1495 : // fill the shape index map
1496 : #ifdef DEBUG
1497 0 : unsigned int max_base_id=0;
1498 0 : unsigned int max_radial_id=0;
1499 : #endif
1500 1066 : for (unsigned int n=0; n<n_total_dofs; ++n)
1501 : {
1502 1040 : 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 0 : if (base_shape_index[n] > max_base_id)
1510 0 : max_base_id = base_shape_index[n];
1511 0 : if (radial_shape_index[n] > max_radial_id)
1512 0 : max_radial_id = radial_shape_index[n];
1513 : #endif
1514 : }
1515 :
1516 : #ifdef DEBUG
1517 0 : libmesh_assert_equal_to( (max_base_id+1)*(max_radial_id+1), n_total_dofs );
1518 : #endif
1519 :
1520 156 : for (auto s : child_elem->side_index_range())
1521 130 : if (child_elem->neighbor_ptr(s) != nullptr &&
1522 130 : child_elem->neighbor_ptr(s) != remote_elem)
1523 130 : 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 48 : std::unique_ptr<const Elem> child_base, parent_base;
1529 48 : child_elem->build_side_ptr(child_base, 0);
1530 48 : parent_elem->build_side_ptr(parent_base, 0);
1531 :
1532 : const unsigned int n_base_dofs =
1533 48 : FEInterface::n_dofs(fe_type, child_base.get());
1534 :
1535 : // We need global DOF indices for both base and 'full' elements
1536 48 : dof_map.dof_indices (child_base.get(), child_base_dof_indices,
1537 : variable_number);
1538 48 : 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 240 : for (unsigned int child_base_dof=0; child_base_dof != n_base_dofs; ++child_base_dof)
1545 : {
1546 0 : libmesh_assert_less (child_base_dof, child_base->n_nodes());
1547 :
1548 : // Childs global dof index.
1549 192 : 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 0 : bool self_constraint = false;
1554 844 : for (unsigned int parent_base_dof=0;
1555 844 : parent_base_dof != n_base_dofs; parent_base_dof++)
1556 : {
1557 0 : libmesh_assert_less (parent_base_dof, parent_base->n_nodes());
1558 :
1559 : // Their global dof index.
1560 700 : const dof_id_type parent_base_dof_g =
1561 : parent_base_dof_indices[parent_base_dof];
1562 :
1563 700 : if (parent_base_dof_g == child_base_dof_g)
1564 : {
1565 0 : self_constraint = true;
1566 0 : break;
1567 : }
1568 : }
1569 :
1570 192 : if (self_constraint)
1571 48 : 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 144 : unsigned int n_elem_dofs = FEInterface::n_dofs(fe_type, child_elem);
1577 0 : libmesh_assert_equal_to(n_elem_dofs, n_total_dofs);
1578 5904 : for(unsigned int child_elem_dof=0; child_elem_dof != n_elem_dofs; ++child_elem_dof)
1579 : {
1580 5760 : if (base_shape_index[child_elem_dof] != child_base_dof)
1581 5566 : 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 1440 : if (s==0)
1588 : {
1589 240 : if (radial_shape_index[child_elem_dof] > 0)
1590 216 : 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 1200 : if ( !child_elem->neighbor_ptr(s)->contains_point(child_base->point(child_base_dof)) )
1597 800 : continue;
1598 : }
1599 :
1600 :
1601 424 : 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 0 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1610 :
1611 424 : if (dof_map.is_constrained_dof(child_elem_dof_g))
1612 0 : continue;
1613 :
1614 194 : constraint_row = &(constraints[child_elem_dof_g]);
1615 0 : libmesh_assert(constraint_row->empty());
1616 : }
1617 :
1618 : // The support point of the DOF
1619 0 : 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 194 : 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 776 : for (unsigned int parent_base_dof=0;
1628 970 : parent_base_dof != n_base_dofs; parent_base_dof++)
1629 : {
1630 :
1631 776 : 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 31816 : for (unsigned int parent_elem_dof=0;
1647 31816 : parent_elem_dof != n_elem_dofs; parent_elem_dof++)
1648 : {
1649 31040 : if (base_shape_index[parent_elem_dof] != parent_base_dof)
1650 30264 : 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 7760 : if (radial_shape_index[parent_elem_dof] != radial_shape_index[child_elem_dof])
1655 6984 : continue;
1656 :
1657 : // Their global dof index.
1658 776 : 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 776 : if ((std::abs(parent_base_dof_value) > 1.e-5) &&
1664 0 : (std::abs(parent_base_dof_value) < .999))
1665 : {
1666 0 : 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 0 : else if (parent_base_dof_value >= .999)
1672 : {
1673 0 : libmesh_assert_equal_to (child_base_dof_g, parent_base_dof_indices[parent_base_dof]);
1674 0 : libmesh_assert_equal_to (child_elem_dof_g, parent_elem_dof_g);
1675 : }
1676 : #endif
1677 : }
1678 :
1679 : }
1680 : }
1681 :
1682 : }
1683 48 : }
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
1696 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const ElemType));
1697 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const ElemType));
1698 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const ElemType));
1699 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const ElemType, const unsigned int));
1700 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const ElemType, const unsigned int));
1701 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const ElemType, const unsigned int));
1702 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_shape_indices(const FEType &, const ElemType, const unsigned int, unsigned int &, unsigned int &));
1703 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_shape_indices(const FEType &, const ElemType, const unsigned int, unsigned int &, unsigned int &));
1704 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_shape_indices(const FEType &, const ElemType, const unsigned int, unsigned int &, unsigned int &));
1705 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_node_indices(const ElemType, const unsigned int, unsigned int &, unsigned int &));
1706 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_node_indices(const ElemType, const unsigned int, unsigned int &, unsigned int &));
1707 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_node_indices(const ElemType, const unsigned int, unsigned int &, unsigned int &));
1708 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
1709 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
1710 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
1711 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape_deriv(const FEType &, const ElemType, const unsigned int, const unsigned int, const Point &));
1712 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape_deriv(const FEType &, const ElemType, const unsigned int, const unsigned int, const Point &));
1713 : INSTANTIATE_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 :
1716 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs(const FEType &, const Elem*));
1717 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs(const FEType &, const Elem*));
1718 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs(const FEType &, const Elem*));
1719 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const Elem *));
1720 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const Elem *));
1721 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const Elem *));
1722 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const Elem *, const unsigned int));
1723 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const Elem *, const unsigned int));
1724 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const Elem *, const unsigned int));
1725 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_shape_indices(const FEType &, const Elem *, const unsigned int, unsigned int &, unsigned int &));
1726 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_shape_indices(const FEType &, const Elem *, const unsigned int, unsigned int &, unsigned int &));
1727 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_shape_indices(const FEType &, const Elem *, const unsigned int, unsigned int &, unsigned int &));
1728 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape(const FEType &, const Elem *, const unsigned int, const Point & p));
1729 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape(const FEType &, const Elem *, const unsigned int, const Point & p));
1730 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape(const FEType &, const Elem *, const unsigned int, const Point & p));
1731 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape(const FEType, const Elem *, const unsigned int, const Point &, const bool));
1732 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape(const FEType, const Elem *, const unsigned int, const Point &, const bool));
1733 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape(const FEType, const Elem *, const unsigned int, const Point &, const bool));
1734 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape_deriv(const FEType &, const Elem *, const unsigned int, const unsigned int, const Point &));
1735 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape_deriv(const FEType &, const Elem *, const unsigned int, const unsigned int, const Point &));
1736 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape_deriv(const FEType &, const Elem *, const unsigned int, const unsigned int, const Point &));
1737 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape_deriv(const FEType, const Elem *, const unsigned int, const unsigned int, const Point &, const bool));
1738 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape_deriv(const FEType, const Elem *, const unsigned int, const unsigned int, const Point &, const bool));
1739 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape_deriv(const FEType, const Elem *, const unsigned int, const unsigned int, const Point &, const bool));
1740 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_data(const FEType &, const Elem *, FEComputeData &));
1741 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_data(const FEType &, const Elem *, FEComputeData &));
1742 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_data(const FEType &, const Elem *, FEComputeData &));
1743 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, nodal_soln(const FEType &, const Elem *, const std::vector<Number> &, std::vector<Number> &));
1744 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, nodal_soln(const FEType &, const Elem *, const std::vector<Number> &, std::vector<Number> &));
1745 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, nodal_soln(const FEType &, const Elem *, const std::vector<Number> &, std::vector<Number> &));
1746 : #ifdef LIBMESH_ENABLE_AMR
1747 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, inf_compute_constraints(DofConstraints &, DofMap &, const unsigned int, const Elem *));
1748 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, inf_compute_constraints(DofConstraints &, DofMap &, const unsigned int, const Elem *));
1749 : INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, inf_compute_constraints(DofConstraints &, DofMap &, const unsigned int, const Elem *));
1750 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1751 : INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, inf_compute_node_constraints(NodeConstraints & constraints,const Elem * elem));
1752 : INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, inf_compute_node_constraints(NodeConstraints & constraints,const Elem * elem));
1753 : INSTANTIATE_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
|
