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