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 :
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 474 : InfFE<Dim,T_radial,T_map>::InfFE (const FEType & fet) :
44 : FEBase (Dim, fet),
45 :
46 107 : calculate_map_scaled(false),
47 107 : calculate_phi_scaled(false),
48 107 : calculate_dphi_scaled(false),
49 107 : calculate_xyz(false),
50 107 : calculate_jxw(false),
51 107 : _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 367 : current_fe_type (FEType(fet.order,
61 474 : fet.family,
62 : INVALID_ORDER,
63 474 : fet.radial_family,
64 474 : fet.inf_map))
65 :
66 : {
67 : // Sanity checks
68 180 : libmesh_assert_equal_to (T_radial, fe_type.radial_family);
69 180 : libmesh_assert_equal_to (T_map, fe_type.inf_map);
70 :
71 : // build the base_fe object
72 : if (Dim != 1)
73 581 : base_fe = FEBase::build(Dim-1, fet);
74 474 : }
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 did not explicitly pre-request something (or nothing)
334 : // to be computed, then we throw an error here.
335 1540 : bool requested_ok =
336 5339 : this->calculate_nothing || this->calculate_phi ||
337 2592 : this->calculate_dphi || this->calculate_dphiref ||
338 2592 : this->calculate_phi_scaled || this->calculate_dphi_scaled ||
339 0 : this->calculate_xyz || this->calculate_jxw ||
340 0 : this->calculate_map_scaled || this->calculate_map ||
341 6879 : this->calculate_curl_phi || this->calculate_div_phi;
342 :
343 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
344 1540 : requested_ok = requested_ok || this->calculate_d2phi;
345 : #endif
346 :
347 1540 : libmesh_error_msg_if(
348 : !requested_ok,
349 : "You must call one or more of the FE accessors "
350 : "(e.g. get_phi(), get_dphi(), get_nothing()) "
351 : "_before_ calling reinit()!");
352 :
353 : // set further terms necessary to do the requested task
354 5339 : if (calculate_jxw)
355 0 : this->calculate_map = true;
356 5339 : if (this->calculate_dphi)
357 85 : this->calculate_map = true;
358 5339 : if (this->calculate_dphi_scaled)
359 2612 : this->calculate_map_scaled = true;
360 : // if Cartesian positions were requested but the calculation of map
361 : // was not triggered, we'll opt for the 'scaled' variant.
362 5339 : if (calculate_xyz && !calculate_map)
363 0 : this->calculate_map_scaled = true;
364 4132 : base_fe->calculate_phi = this->calculate_phi || this->calculate_phi_scaled
365 6203 : || this->calculate_dphi || this->calculate_dphi_scaled;
366 5339 : base_fe->calculate_dphi = this->calculate_dphi || this->calculate_dphi_scaled;
367 5339 : if (this->calculate_map || this->calculate_map_scaled
368 2652 : || this->calculate_dphiref)
369 : {
370 2687 : base_fe->calculate_dphiref = true;
371 2687 : base_fe->get_xyz(); // trigger base_fe->fe_map to 'calculate_xyz'
372 2687 : base_fe->get_JxW(); // trigger base_fe->fe_map to 'calculate_dxyz'
373 : }
374 5339 : base_fe->determine_calculations();
375 :
376 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
377 5339 : if (this->calculate_d2phi)
378 0 : libmesh_not_implemented();
379 : #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
380 5339 : }
381 :
382 :
383 :
384 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
385 : void
386 1431 : InfFE<Dim, T_radial, T_map>::
387 : init_radial_shape_functions(const Elem * libmesh_dbg_var(inf_elem),
388 : const std::vector<Point> * radial_pts)
389 : {
390 361 : libmesh_assert(radial_qrule.get() || radial_pts);
391 361 : libmesh_assert(inf_elem);
392 :
393 : // Start logging the radial shape function initialization
394 722 : LOG_SCOPE("init_radial_shape_functions()", "InfFE");
395 :
396 : // initialize most of the things related to physical approximation
397 722 : const Order radial_approx_order = fe_type.radial_order;
398 : const unsigned int n_radial_approx_shape_functions =
399 361 : InfFERadial::n_dofs(radial_approx_order);
400 :
401 722 : const std::size_t n_radial_qp =
402 1395 : radial_pts ? radial_pts->size() : radial_qrule->n_points();
403 722 : const std::vector<Point> & radial_qp =
404 36 : radial_pts ? *radial_pts : radial_qrule->get_points();
405 :
406 : // the radial polynomials (eval)
407 1431 : if (calculate_phi || calculate_dphi || calculate_phi_scaled || calculate_dphi_scaled)
408 : {
409 1431 : mode.resize (n_radial_approx_shape_functions);
410 9490 : for (unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
411 9477 : mode[i].resize (n_radial_qp);
412 :
413 : // evaluate the mode shapes in radial direction at radial quadrature points
414 9490 : for (unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
415 17848 : for (std::size_t p=0; p<n_radial_qp; ++p)
416 11899 : mode[i][p] = InfFE<Dim,T_radial,T_map>::eval (radial_qp[p](0), radial_approx_order, i);
417 : }
418 :
419 1431 : if (calculate_dphi || calculate_dphi_scaled)
420 : {
421 95 : dmodedv.resize (n_radial_approx_shape_functions);
422 430 : for (unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
423 469 : dmodedv[i].resize (n_radial_qp);
424 :
425 : // evaluate the mode shapes in radial direction at radial quadrature points
426 430 : for (unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
427 2400 : for (std::size_t p=0; p<n_radial_qp; ++p)
428 2891 : dmodedv[i][p] = InfFE<Dim,T_radial,T_map>::eval_deriv (radial_qp[p](0), radial_approx_order, i);
429 : }
430 :
431 : // the (1-v)/2 weight.
432 1431 : if (calculate_phi || calculate_phi_scaled || calculate_dphi || calculate_dphi_scaled)
433 : {
434 1431 : som.resize (n_radial_qp);
435 : // compute scalar values at radial quadrature points
436 3327 : for (std::size_t p=0; p<n_radial_qp; ++p)
437 2443 : som[p] = InfFERadial::decay (Dim, radial_qp[p](0));
438 : }
439 1431 : if (calculate_dphi || calculate_dphi_scaled)
440 : {
441 95 : dsomdv.resize (n_radial_qp);
442 : // compute scalar values at radial quadrature points
443 655 : for (std::size_t p=0; p<n_radial_qp; ++p)
444 784 : dsomdv[p] = InfFERadial::decay_deriv (Dim, radial_qp[p](0));
445 : }
446 1431 : }
447 :
448 :
449 :
450 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
451 2052 : void InfFE<Dim,T_radial,T_map>::init_shape_functions(const std::vector<Point> & radial_qp,
452 : const std::vector<Point> & base_qp,
453 : const Elem * inf_elem)
454 : {
455 567 : libmesh_assert(inf_elem);
456 :
457 : // Start logging the radial shape function initialization
458 1134 : LOG_SCOPE("init_shape_functions()", "InfFE");
459 :
460 : // fast access to some const ints for the radial data
461 1134 : const unsigned int n_radial_approx_sf = InfFERadial::n_dofs(fe_type.radial_order);
462 1134 : const std::size_t n_radial_qp = radial_qp.size();
463 : #ifdef DEBUG
464 567 : if (calculate_phi || calculate_dphi || calculate_phi_scaled || calculate_dphi_scaled)
465 567 : libmesh_assert_equal_to(n_radial_approx_sf, mode.size());
466 567 : if (calculate_phi || calculate_dphi || calculate_dphi_scaled)
467 567 : libmesh_assert_equal_to(som.size(), n_radial_qp);
468 : #endif
469 :
470 :
471 : // initialize most of the quantities related to mapping
472 :
473 : // The element type and order to use in the base map
474 : //const Order base_mapping_order = base_elem->default_order();
475 :
476 : // the number of base shape functions used to construct the map
477 : // (Lagrange shape functions are used for mapping in the base)
478 : //unsigned int n_base_mapping_shape_functions =
479 : // InfFEBase::n_base_mapping_sf(*base_elem,
480 : // base_mapping_order);
481 :
482 : // initialize most of the things related to physical approximation
483 : unsigned int n_base_approx_shape_functions;
484 : if (Dim > 1)
485 567 : n_base_approx_shape_functions =
486 2052 : FEInterface::n_dofs(base_fe->get_fe_type(), base_elem.get());
487 : else
488 0 : n_base_approx_shape_functions = 1;
489 :
490 :
491 : // update class member field
492 2052 : _n_total_approx_sf =
493 2052 : n_radial_approx_sf * n_base_approx_shape_functions;
494 :
495 :
496 : // The number of the base quadrature points.
497 1134 : const unsigned int n_base_qp = cast_int<unsigned int>(base_qp.size());
498 :
499 : // The total number of quadrature points.
500 2052 : _n_total_qp = n_radial_qp * n_base_qp;
501 :
502 :
503 : // initialize the node and shape numbering maps
504 : {
505 : // similar for the shapes: the i-th entry stores
506 : // the associated base/radial shape number
507 2052 : _radial_shape_index.resize(_n_total_approx_sf);
508 2052 : _base_shape_index.resize(_n_total_approx_sf);
509 :
510 : // fill the shape index map
511 57883 : for (unsigned int n=0; n<_n_total_approx_sf; ++n)
512 : {
513 82871 : compute_shape_indices (this->fe_type,
514 : inf_elem,
515 : n,
516 : _base_shape_index[n],
517 : _radial_shape_index[n]);
518 13520 : libmesh_assert_less (_base_shape_index[n], n_base_approx_shape_functions);
519 13520 : libmesh_assert_less (_radial_shape_index[n], n_radial_approx_sf);
520 : }
521 : }
522 :
523 : // resize the base data fields
524 : //dist.resize(n_base_mapping_shape_functions);
525 :
526 : // resize the total data fields
527 :
528 : // the phase term varies with xi, eta and zeta(v): store it for _all_ qp
529 : //
530 : // when computing the phase, we need the base approximations
531 : // therefore, initialize the phase here, but evaluate it
532 : // in compute_shape_functions().
533 : //
534 : // the weight, though, is only needed at the radial quadrature points, n_radial_qp.
535 : // but for a uniform interface to the protected data fields
536 : // the weight data field (which are accessible from the outside) are expanded to _n_total_qp.
537 2052 : if (calculate_phi || calculate_dphi)
538 1406 : weight.resize (_n_total_qp);
539 2052 : if (calculate_phi_scaled || calculate_dphi_scaled)
540 656 : weightxr_sq.resize (_n_total_qp);
541 2052 : if (calculate_dphi || calculate_dphi_scaled)
542 721 : dweightdv.resize (n_radial_qp);
543 2052 : if (calculate_dphi)
544 75 : dweight.resize (_n_total_qp);
545 2052 : if (calculate_dphi_scaled)
546 656 : dweightxr_sq.resize(_n_total_qp);
547 :
548 2052 : if (calculate_dphi || calculate_dphi_scaled)
549 721 : dphase.resize (_n_total_qp);
550 :
551 : // this vector contains the integration weights for the combined quadrature rule
552 : // if no quadrature rules are given, use only ones.
553 2052 : _total_qrule_weights.resize(_n_total_qp,1);
554 :
555 : // InfFE's data fields phi, dphi, dphidx, phi_map etc hold the _total_
556 : // shape and mapping functions, respectively
557 : {
558 2052 : if (calculate_map_scaled)
559 661 : JxWxdecay.resize(_n_total_qp);
560 2052 : if (calculate_jxw)
561 0 : JxW.resize(_n_total_qp);
562 2052 : if (calculate_map_scaled || calculate_map)
563 : {
564 726 : xyz.resize(_n_total_qp);
565 726 : dxidx_map_scaled.resize(_n_total_qp);
566 726 : dxidy_map_scaled.resize(_n_total_qp);
567 726 : dxidz_map_scaled.resize(_n_total_qp);
568 726 : detadx_map_scaled.resize(_n_total_qp);
569 726 : detady_map_scaled.resize(_n_total_qp);
570 726 : detadz_map_scaled.resize(_n_total_qp);
571 726 : dzetadx_map_scaled.resize(_n_total_qp);
572 726 : dzetady_map_scaled.resize(_n_total_qp);
573 726 : dzetadz_map_scaled.resize(_n_total_qp);
574 : }
575 2052 : if (calculate_map)
576 : {
577 80 : dxidx_map.resize(_n_total_qp);
578 80 : dxidy_map.resize(_n_total_qp);
579 80 : dxidz_map.resize(_n_total_qp);
580 80 : detadx_map.resize(_n_total_qp);
581 80 : detady_map.resize(_n_total_qp);
582 80 : detadz_map.resize(_n_total_qp);
583 80 : dzetadx_map.resize(_n_total_qp);
584 80 : dzetady_map.resize(_n_total_qp);
585 80 : dzetadz_map.resize(_n_total_qp);
586 : }
587 2052 : if (calculate_phi)
588 1406 : phi.resize (_n_total_approx_sf);
589 2052 : if (calculate_phi_scaled)
590 656 : phixr.resize (_n_total_approx_sf);
591 2052 : if (calculate_dphi)
592 : {
593 75 : dphi.resize (_n_total_approx_sf);
594 75 : dphidx.resize (_n_total_approx_sf);
595 75 : dphidy.resize (_n_total_approx_sf);
596 75 : dphidz.resize (_n_total_approx_sf);
597 : }
598 :
599 2052 : if (calculate_dphi_scaled)
600 : {
601 656 : dphixr.resize (_n_total_approx_sf);
602 656 : dphixr_sq.resize(_n_total_approx_sf);
603 : }
604 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
605 :
606 2052 : if (calculate_d2phi)
607 : {
608 0 : libmesh_not_implemented();
609 : d2phi.resize (_n_total_approx_sf);
610 : d2phidx2.resize (_n_total_approx_sf);
611 : d2phidxdy.resize (_n_total_approx_sf);
612 : d2phidxdz.resize (_n_total_approx_sf);
613 : d2phidy2.resize (_n_total_approx_sf);
614 : d2phidydz.resize (_n_total_approx_sf);
615 : d2phidz2.resize (_n_total_approx_sf);
616 : d2phidxi2.resize (_n_total_approx_sf);
617 :
618 : if (Dim > 1)
619 : {
620 : d2phidxideta.resize(_n_total_approx_sf);
621 : d2phideta2.resize(_n_total_approx_sf);
622 : }
623 :
624 : if (Dim > 2)
625 : {
626 : d2phidetadzeta.resize(_n_total_approx_sf);
627 : d2phidxidzeta.resize(_n_total_approx_sf);
628 : d2phidzeta2.resize(_n_total_approx_sf);
629 : }
630 : }
631 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
632 :
633 2052 : if (calculate_dphi || calculate_dphi_scaled)
634 : {
635 721 : dphidxi.resize (_n_total_approx_sf);
636 :
637 : if (Dim > 1)
638 721 : dphideta.resize(_n_total_approx_sf);
639 :
640 : if (Dim == 3)
641 721 : dphidzeta.resize(_n_total_approx_sf);
642 : }
643 :
644 : }
645 :
646 : // collect all the for loops, where inner vectors are
647 : // resized to the appropriate number of quadrature points
648 : {
649 2052 : if (calculate_phi)
650 33232 : for (unsigned int i=0; i<_n_total_approx_sf; ++i)
651 37334 : phi[i].resize (_n_total_qp);
652 :
653 2052 : if (calculate_dphi)
654 1025 : for (unsigned int i=0; i<_n_total_approx_sf; ++i)
655 : {
656 1330 : dphi[i].resize (_n_total_qp);
657 1330 : dphidx[i].resize (_n_total_qp);
658 1330 : dphidy[i].resize (_n_total_qp);
659 1330 : dphidz[i].resize (_n_total_qp);
660 : }
661 :
662 2052 : if (calculate_phi_scaled)
663 24741 : for (unsigned int i=0; i<_n_total_approx_sf; ++i)
664 : {
665 32129 : phixr[i].resize (_n_total_qp);
666 : }
667 2052 : if (calculate_dphi_scaled)
668 24741 : for (unsigned int i=0; i<_n_total_approx_sf; ++i)
669 : {
670 32129 : dphixr[i].resize(_n_total_qp);
671 32129 : dphixr_sq[i].resize(_n_total_qp);
672 : }
673 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
674 2052 : if (calculate_d2phi)
675 0 : for (unsigned int i=0; i<_n_total_approx_sf; ++i)
676 : {
677 0 : d2phi[i].resize (_n_total_qp);
678 0 : d2phidx2[i].resize (_n_total_qp);
679 0 : d2phidxdy[i].resize (_n_total_qp);
680 0 : d2phidxdz[i].resize (_n_total_qp);
681 0 : d2phidy2[i].resize (_n_total_qp);
682 0 : d2phidydz[i].resize (_n_total_qp);
683 0 : d2phidy2[i].resize (_n_total_qp);
684 0 : d2phidxi2[i].resize (_n_total_qp);
685 :
686 : if (Dim > 1)
687 : {
688 0 : d2phidxideta[i].resize (_n_total_qp);
689 0 : d2phideta2[i].resize (_n_total_qp);
690 : }
691 : if (Dim > 2)
692 : {
693 0 : d2phidxidzeta[i].resize (_n_total_qp);
694 0 : d2phidetadzeta[i].resize (_n_total_qp);
695 0 : d2phidzeta2[i].resize (_n_total_qp);
696 : }
697 : }
698 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
699 :
700 2052 : if (calculate_dphi || calculate_dphi_scaled)
701 25676 : for (unsigned int i=0; i<_n_total_approx_sf; ++i)
702 : {
703 33347 : dphidxi[i].resize (_n_total_qp);
704 :
705 : if (Dim > 1)
706 33347 : dphideta[i].resize (_n_total_qp);
707 :
708 : if (Dim == 3)
709 33347 : dphidzeta[i].resize (_n_total_qp);
710 :
711 : }
712 :
713 : }
714 : {
715 : // (a) compute scalar values at _all_ quadrature points -- for uniform
716 : // access from the outside to these fields
717 : // (b) form a std::vector<Real> which contains the appropriate weights
718 : // of the combined quadrature rule!
719 567 : libmesh_assert_equal_to (radial_qp.size(), n_radial_qp);
720 :
721 2052 : if (radial_qrule && base_qrule)
722 : {
723 244 : const std::vector<Real> & radial_qw = radial_qrule->get_weights();
724 244 : const std::vector<Real> & base_qw = base_qrule->get_weights();
725 244 : libmesh_assert_equal_to (radial_qw.size(), n_radial_qp);
726 244 : libmesh_assert_equal_to (base_qw.size(), n_base_qp);
727 :
728 5588 : for (unsigned int rp=0; rp<n_radial_qp; ++rp)
729 97938 : for (unsigned int bp=0; bp<n_base_qp; ++bp)
730 186492 : _total_qrule_weights[bp + rp*n_base_qp] = radial_qw[rp] * base_qw[bp];
731 : }
732 :
733 :
734 8325 : for (unsigned int rp=0; rp<n_radial_qp; ++rp)
735 : {
736 6273 : if (calculate_phi || calculate_dphi)
737 4717 : for (unsigned int bp=0; bp<n_base_qp; ++bp)
738 5880 : weight[bp + rp*n_base_qp] = InfFERadial::D(radial_qp[rp](0));
739 :
740 6273 : if ( calculate_phi_scaled)
741 96423 : for (unsigned int bp=0; bp<n_base_qp; ++bp)
742 122479 : weightxr_sq[bp + rp*n_base_qp] = InfFERadial::Dxr_sq(radial_qp[rp](0));
743 :
744 6273 : if ( calculate_dphi || calculate_dphi_scaled)
745 9982 : dweightdv[rp] = InfFERadial::D_deriv(radial_qp[rp](0));
746 : }
747 : }
748 2052 : }
749 :
750 :
751 :
752 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
753 3998 : void InfFE<Dim,T_radial,T_map>::compute_shape_functions(const Elem * inf_elem,
754 : const std::vector<Point> & base_qp,
755 : const std::vector<Point> & radial_qp
756 : )
757 : {
758 1215 : libmesh_assert(inf_elem);
759 : // at least check whether the base element type is correct.
760 : // otherwise this version of computing dist would give problems
761 1215 : libmesh_assert_equal_to (base_elem->type(),
762 : InfFEBase::get_elem_type(inf_elem->type()));
763 :
764 : // Start logging the overall computation of shape functions
765 2430 : LOG_SCOPE("compute_shape_functions()", "InfFE");
766 :
767 : //const unsigned int n_radial_qp = cast_int<unsigned int>(som.size());
768 : //const unsigned int n_base_qp = cast_int<unsigned int>(S_map[0].size());
769 2430 : const std::size_t n_radial_qp = radial_qp.size();
770 3998 : const unsigned int n_base_qp = base_qp.size();
771 :
772 1215 : libmesh_assert_equal_to (_n_total_qp, n_radial_qp*n_base_qp);
773 1215 : libmesh_assert_equal_to (_n_total_qp, _total_qrule_weights.size());
774 : #ifdef DEBUG
775 1215 : if (som.size() > 0)
776 1215 : libmesh_assert_equal_to(n_radial_qp, som.size());
777 :
778 1215 : if (this->calculate_map || this->calculate_map_scaled)
779 : {
780 : // these vectors are needed later; initialize here already to have access to
781 : // n_base_qp etc.
782 896 : const std::vector<std::vector<Real>> & S_map = (base_fe->get_fe_map()).get_phi_map();
783 896 : if (S_map[0].size() > 0)
784 896 : libmesh_assert_equal_to(n_base_qp, S_map[0].size());
785 : }
786 1215 : if (radial_qrule)
787 892 : libmesh_assert_equal_to(n_radial_qp, radial_qrule->n_points());
788 1215 : if (base_qrule)
789 892 : libmesh_assert_equal_to(n_base_qp, base_qrule->n_points());
790 1215 : libmesh_assert_equal_to(_n_total_qp % n_radial_qp, 0); // "Error in the structure of quadrature points!");
791 : #endif
792 :
793 :
794 3998 : _n_total_approx_sf = InfFERadial::n_dofs(fe_type.radial_order) *
795 3998 : FEInterface::n_dofs(base_fe->get_fe_type(), base_elem.get());
796 :
797 :
798 :
799 3998 : const Point origin = inf_elem->origin();
800 :
801 : // Compute the shape function values (and derivatives)
802 : // at the Quadrature points. Note that the actual values
803 : // have already been computed via init_shape_functions
804 :
805 3998 : unsigned int elem_dim = inf_elem->dim();
806 : // Compute the value of the derivative shape function i at quadrature point p
807 3998 : switch (elem_dim)
808 : {
809 0 : case 1:
810 : case 2:
811 : {
812 0 : libmesh_not_implemented();
813 : break;
814 : }
815 1215 : case 3:
816 : {
817 6428 : std::vector<std::vector<Real>> S (0);
818 6428 : std::vector<std::vector<Real>> Ss(0);
819 6428 : std::vector<std::vector<Real>> St(0);
820 :
821 5213 : std::vector<Real> base_dxidx (0);
822 5213 : std::vector<Real> base_dxidy (0);
823 5213 : std::vector<Real> base_dxidz (0);
824 5213 : std::vector<Real> base_detadx(0);
825 5213 : std::vector<Real> base_detady(0);
826 5213 : std::vector<Real> base_detadz(0);
827 :
828 5213 : std::vector<Point> base_xyz (0);
829 :
830 3998 : if (calculate_phi || calculate_phi_scaled ||
831 0 : calculate_dphi || calculate_dphi_scaled)
832 3998 : S=base_fe->phi;
833 :
834 : // fast access to the approximation and mapping shapes of base_fe
835 3998 : if (calculate_map || calculate_map_scaled)
836 : {
837 2672 : Ss = base_fe->dphidxi;
838 2672 : St = base_fe->dphideta;
839 :
840 2672 : base_dxidx = base_fe->get_dxidx();
841 2672 : base_dxidy = base_fe->get_dxidy();
842 2672 : base_dxidz = base_fe->get_dxidz();
843 2672 : base_detadx = base_fe->get_detadx();
844 2672 : base_detady = base_fe->get_detady();
845 2672 : base_detadz = base_fe->get_detadz();
846 :
847 2672 : base_xyz = base_fe->get_xyz();
848 : }
849 :
850 3998 : ElemType base_type= base_elem->type();
851 :
852 : #ifdef DEBUG
853 1215 : if (calculate_phi)
854 351 : libmesh_assert_equal_to (phi.size(), _n_total_approx_sf);
855 1215 : if (calculate_dphi)
856 : {
857 30 : libmesh_assert_equal_to (dphidxi.size(), _n_total_approx_sf);
858 30 : libmesh_assert_equal_to (dphideta.size(), _n_total_approx_sf);
859 30 : libmesh_assert_equal_to (dphidzeta.size(), _n_total_approx_sf);
860 : }
861 : #endif
862 :
863 1215 : unsigned int tp=0; // total qp
864 22746 : for (unsigned int rp=0; rp<n_radial_qp; ++rp) // over radial qps
865 244657 : for (unsigned int bp=0; bp<n_base_qp; ++bp) // over base qps
866 :
867 : { // First compute the map from base element quantities to physical space:
868 :
869 : // initialize them with invalid value to not use them
870 : // without setting them to the correct value before.
871 75289 : Point unit_r(NAN);
872 75289 : RealGradient grad_a_scaled(NAN);
873 75289 : Real a(NAN);
874 75289 : Real r_norm(NAN);
875 225909 : if (calculate_map || calculate_map_scaled)
876 : {
877 449493 : xyz[tp] = InfFEMap::map(elem_dim, inf_elem, Point(base_qp[bp](0),base_qp[bp](1),radial_qp[rp](0)));
878 :
879 74970 : const Point r(xyz[tp]-origin);
880 224583 : a=(base_xyz[bp]-origin).norm();
881 149613 : r_norm = r.norm();
882 :
883 : // check that 'som' == a/r.
884 : #ifndef NDEVEL
885 74970 : if (som.size())
886 74970 : libmesh_assert_less(std::abs(som[rp] -a/r_norm) , 1e-7);
887 : #endif
888 74970 : unit_r=(r/r_norm);
889 :
890 : // They are used for computing the normal and do not correspond to the direction of eta and xi in this element:
891 : // Due to the stretch of these axes in radial direction, they are deformed.
892 374523 : Point e_xi(base_dxidx[bp],
893 : base_dxidy[bp],
894 : base_dxidz[bp]);
895 449493 : Point e_eta(base_detadx[bp],
896 : base_detady[bp],
897 : base_detadz[bp]);
898 :
899 224583 : const RealGradient normal=e_eta.cross(e_xi).unit();
900 :
901 : // grad a = a/r.norm() * grad_a_scaled
902 74970 : grad_a_scaled=unit_r - normal/(normal*unit_r);
903 :
904 374523 : const Real dxi_er=base_dxidx[bp]* unit_r(0) + base_dxidy[bp] *unit_r(1) + base_dxidz[bp] *unit_r(2);
905 449493 : const Real deta_er=base_detadx[bp]*unit_r(0) + base_detady[bp]*unit_r(1) + base_detadz[bp]*unit_r(2);
906 :
907 : // in case of non-affine map, further terms need to be taken into account,
908 : // involving \p e_eta and \p e_xi and thus recursive computation is needed
909 224583 : if (!base_elem->has_affine_map())
910 : {
911 : /**
912 : * The full form for 'a' is
913 : * a = (r0*normal)/(normal*unit_r);
914 : * where r0 is some point on the base plane(!)
915 : * when the base element is not a plane, r0 and normal are functions of space.
916 : * Here, some approximation is used:
917 : */
918 :
919 200 : const unsigned int n_sf = base_elem->n_nodes();
920 80 : RealGradient tmp(0.,0.,0.);
921 2000 : for (unsigned int i=0; i< n_sf; ++i)
922 : {
923 3240 : RealGradient dL_da_i = (FE<2,LAGRANGE>::shape_deriv(base_type,
924 1800 : base_elem->default_order(),
925 720 : i, 0, base_qp[bp]) * e_xi
926 3240 : +FE<2,LAGRANGE>::shape_deriv(base_type,
927 1800 : base_elem->default_order(),
928 720 : i, 1, base_qp[bp]) * e_eta);
929 :
930 1080 : tmp += (base_elem->node_ref(i) -origin).norm()* dL_da_i;
931 :
932 : }
933 80 : libmesh_assert(tmp*unit_r < .95 ); // in a proper setup, tmp should have only a small radial component.
934 200 : grad_a_scaled = ( tmp - (tmp*unit_r)*unit_r ) / ( 1. - tmp*unit_r);
935 :
936 : }
937 :
938 : // 'scale' = r/a
939 299553 : dxidx_map_scaled[tp] = (grad_a_scaled(0) - unit_r(0))*dxi_er +base_dxidx[bp];
940 299553 : dxidy_map_scaled[tp] = (grad_a_scaled(1) - unit_r(1))*dxi_er +base_dxidy[bp];
941 299553 : dxidz_map_scaled[tp] = (grad_a_scaled(2) - unit_r(2))*dxi_er +base_dxidz[bp];
942 :
943 : // 'scale' = r/a
944 299553 : detadx_map_scaled[tp] = (grad_a_scaled(0) - unit_r(0))*deta_er + base_detadx[bp];
945 299553 : detady_map_scaled[tp] = (grad_a_scaled(1) - unit_r(1))*deta_er + base_detady[bp];
946 299553 : detadz_map_scaled[tp] = (grad_a_scaled(2) - unit_r(2))*deta_er + base_detadz[bp];
947 :
948 : // 'scale' = (r/a)**2
949 224583 : dzetadx_map_scaled[tp] =-2./a*(grad_a_scaled(0) - unit_r(0));
950 224583 : dzetady_map_scaled[tp] =-2./a*(grad_a_scaled(1) - unit_r(1));
951 299553 : dzetadz_map_scaled[tp] =-2./a*(grad_a_scaled(2) - unit_r(2));
952 :
953 : }
954 :
955 225909 : if (calculate_map)
956 : {
957 2289 : dxidx_map[tp] = a/r_norm * dxidx_map_scaled[tp];
958 2289 : dxidy_map[tp] = a/r_norm * dxidy_map_scaled[tp];
959 2289 : dxidz_map[tp] = a/r_norm * dxidz_map_scaled[tp];
960 :
961 2289 : detadx_map[tp] = a/r_norm * detadx_map_scaled[tp];
962 2289 : detady_map[tp] = a/r_norm * detady_map_scaled[tp];
963 2289 : detadz_map[tp] = a/r_norm * detadz_map_scaled[tp];
964 :
965 : // dzetadx = dzetadr*dr/dx - 2/r * grad_a
966 : // = dzetadr*dr/dx - 2*a/r^2 * grad_a_scaled
967 1635 : dzetadx_map[tp] =-2.*a/(r_norm*r_norm)*(grad_a_scaled(0) - unit_r(0));
968 1635 : dzetady_map[tp] =-2.*a/(r_norm*r_norm)*(grad_a_scaled(1) - unit_r(1));
969 1635 : dzetadz_map[tp] =-2.*a/(r_norm*r_norm)*(grad_a_scaled(2) - unit_r(2));
970 :
971 1635 : if (calculate_jxw)
972 : {
973 0 : Real inv_jac = (dxidx_map[tp]*( detady_map[tp]*dzetadz_map[tp]- dzetady_map[tp]*detadz_map[tp]) +
974 0 : detadx_map[tp]*(dzetady_map[tp]* dxidz_map[tp]- dxidy_map[tp]*dzetadz_map[tp]) +
975 0 : dzetadx_map[tp]*( dxidy_map[tp]*detadz_map[tp]- detady_map[tp]* dxidz_map[tp]));
976 :
977 0 : if (inv_jac <= 1e-10)
978 : {
979 0 : libmesh_error_msg("ERROR: negative inverse Jacobian " \
980 : << inv_jac \
981 : << " at point " \
982 : << xyz[tp] \
983 : << " in element " \
984 : << inf_elem->id());
985 : }
986 :
987 :
988 0 : JxW[tp] = _total_qrule_weights[tp]/inv_jac;
989 : }
990 :
991 : }
992 225909 : if (calculate_map_scaled)
993 : {
994 372003 : Real inv_jacxR_pow4 = (dxidx_map_scaled[tp] *( detady_map_scaled[tp]*dzetadz_map_scaled[tp]
995 297593 : - dzetady_map_scaled[tp]* detadz_map_scaled[tp]) +
996 297593 : detadx_map_scaled[tp] *( dzetady_map_scaled[tp]* dxidz_map_scaled[tp]
997 223183 : - dxidy_map_scaled[tp]*dzetadz_map_scaled[tp]) +
998 223183 : dzetadx_map_scaled[tp]*( dxidy_map_scaled[tp]* detadz_map_scaled[tp]
999 223183 : -detady_map_scaled[tp]* dxidz_map_scaled[tp]));
1000 223183 : if (inv_jacxR_pow4 <= 1e-7)
1001 : {
1002 0 : libmesh_error_msg("ERROR: negative weighted inverse Jacobian " \
1003 : << inv_jacxR_pow4 \
1004 : << " at point " \
1005 : << xyz[tp] \
1006 : << " in element " \
1007 : << inf_elem->id());
1008 : }
1009 :
1010 372003 : JxWxdecay[tp] = _total_qrule_weights[tp]/inv_jacxR_pow4;
1011 : }
1012 :
1013 : // phase term mu(r)=i*k*(r-a).
1014 : // skip i*k: it is added separately during matrix assembly.
1015 :
1016 225909 : if (calculate_dphi || calculate_dphi_scaled)
1017 299504 : dphase[tp] = unit_r - grad_a_scaled*a/r_norm;
1018 :
1019 225909 : if (calculate_dphi)
1020 : {
1021 2880 : dweight[tp](0) = dweightdv[rp] * dzetadx_map[tp];
1022 1600 : dweight[tp](1) = dweightdv[rp] * dzetady_map[tp];
1023 2240 : dweight[tp](2) = dweightdv[rp] * dzetadz_map[tp];
1024 : }
1025 225909 : if (calculate_dphi_scaled)
1026 : {
1027 371940 : dweightxr_sq[tp](0) = dweightdv[rp] * dzetadx_map_scaled[tp];
1028 223148 : dweightxr_sq[tp](1) = dweightdv[rp] * dzetady_map_scaled[tp];
1029 297544 : dweightxr_sq[tp](2) = dweightdv[rp] * dzetadz_map_scaled[tp];
1030 : }
1031 :
1032 225909 : if (calculate_phi || calculate_phi_scaled || calculate_dphi || calculate_dphi_scaled)
1033 : // compute the shape-functions and derivative quantities:
1034 7997773 : for (unsigned int i=0; i <_n_total_approx_sf ; ++i)
1035 : {
1036 : // let the index vectors take care of selecting the appropriate base/radial shape
1037 7771864 : unsigned int bi = _base_shape_index [i];
1038 7771864 : unsigned int ri = _radial_shape_index[i];
1039 7771864 : if (calculate_phi)
1040 150828 : phi [i][tp] = S [bi][bp] * mode[ri][rp] * som[rp];
1041 :
1042 7771864 : if (calculate_phi_scaled)
1043 23162044 : phixr [i][tp] = S [bi][bp] * mode[ri][rp];
1044 :
1045 7771864 : if (calculate_dphi || calculate_dphi_scaled)
1046 : {
1047 23230588 : dphidxi [i][tp] = Ss[bi][bp] * mode[ri][rp] * som[rp];
1048 15485608 : dphideta [i][tp] = St[bi][bp] * mode[ri][rp] * som[rp];
1049 10322288 : dphidzeta[i][tp] = S [bi][bp]
1050 20648928 : * (dmodedv[ri][rp] * som[rp] + mode[ri][rp] * dsomdv[rp]);
1051 : }
1052 :
1053 7771864 : if (calculate_dphi)
1054 : {
1055 :
1056 : // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
1057 30464 : dphi[i][tp](0) =
1058 56576 : dphidx[i][tp] = (dphidxi[i][tp]*dxidx_map[tp] +
1059 39168 : dphideta[i][tp]*detadx_map[tp] +
1060 47872 : dphidzeta[i][tp]*dzetadx_map[tp]);
1061 :
1062 : // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
1063 21760 : dphi[i][tp](1) =
1064 39168 : dphidy[i][tp] = (dphidxi[i][tp]*dxidy_map[tp] +
1065 21760 : dphideta[i][tp]*detady_map[tp] +
1066 30464 : dphidzeta[i][tp]*dzetady_map[tp]);
1067 :
1068 : // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
1069 21760 : dphi[i][tp](2) =
1070 39168 : dphidz[i][tp] = (dphidxi[i][tp]*dxidz_map[tp] +
1071 21760 : dphideta[i][tp]*detadz_map[tp] +
1072 30464 : dphidzeta[i][tp]*dzetadz_map[tp]);
1073 :
1074 : }
1075 7771864 : if (calculate_dphi_scaled)
1076 : { // we don't distinguish between the different levels of scaling here...
1077 :
1078 18014852 : dphixr[i][tp](0)= (dphidxi[i][tp]*dxidx_map_scaled[tp] +
1079 12867660 : dphideta[i][tp]*detadx_map_scaled[tp] +
1080 18014852 : dphidzeta[i][tp]*dzetadx_map_scaled[tp]*som[rp]);
1081 :
1082 10294064 : dphixr[i][tp](1) = (dphidxi[i][tp]*dxidy_map_scaled[tp] +
1083 7720468 : dphideta[i][tp]*detady_map_scaled[tp] +
1084 7720468 : dphidzeta[i][tp]*dzetady_map_scaled[tp]*som[rp]);
1085 :
1086 10294064 : dphixr[i][tp](2) = (dphidxi[i][tp]*dxidz_map_scaled[tp] +
1087 7720468 : dphideta[i][tp]*detadz_map_scaled[tp] +
1088 7720468 : dphidzeta[i][tp]*dzetadz_map_scaled[tp]*som[rp]);
1089 :
1090 15441256 : const Real dphidxixr = Ss[bi][bp] * mode[ri][rp];
1091 10294064 : const Real dphidetaxr= St[bi][bp] * mode[ri][rp];
1092 :
1093 10294064 : dphixr_sq[i][tp](0)= (dphidxixr*dxidx_map_scaled[tp] +
1094 10294064 : dphidetaxr*detadx_map_scaled[tp] +
1095 7720468 : dphidzeta[i][tp]*dzetadx_map_scaled[tp]);
1096 :
1097 15441256 : dphixr_sq[i][tp](1) = (dphidxixr*dxidy_map_scaled[tp] +
1098 10294064 : dphidetaxr*detady_map_scaled[tp] +
1099 7720468 : dphidzeta[i][tp]*dzetady_map_scaled[tp]);
1100 :
1101 15441256 : dphixr_sq[i][tp](2) = (dphidxixr*dxidz_map_scaled[tp] +
1102 10294064 : dphidetaxr*detadz_map_scaled[tp] +
1103 7720468 : dphidzeta[i][tp]*dzetadz_map_scaled[tp]);
1104 : }
1105 :
1106 : }
1107 225909 : tp++;
1108 : }
1109 :
1110 1215 : break;
1111 1568 : }
1112 :
1113 0 : default:
1114 0 : libmesh_error_msg("Unsupported dim = " << dim);
1115 : }
1116 3998 : }
1117 :
1118 :
1119 :
1120 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
1121 0 : bool InfFE<Dim,T_radial,T_map>::shapes_need_reinit() const
1122 : {
1123 : // We never call this.
1124 0 : libmesh_not_implemented();
1125 : return false;
1126 : }
1127 :
1128 : } // namespace libMesh
1129 :
1130 :
1131 : // Explicit instantiations
1132 : #include "libmesh/inf_fe_instantiate_1D.h"
1133 : #include "libmesh/inf_fe_instantiate_2D.h"
1134 : #include "libmesh/inf_fe_instantiate_3D.h"
1135 :
1136 :
1137 :
1138 : #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
|
