1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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/quadrature_gauss.h"
27 #include "libmesh/elem.h"
28 #include "libmesh/libmesh_logging.h"
29 #include "libmesh/int_range.h"
30 #include "libmesh/auto_ptr.h"
31
32
33 namespace {
34 // Put this outside a templated class, so we only get 1 warning
35 // during our unit tests, not 1 warning for each of the zillion FE
36 // specializations we test.
inffe_hessian_warning()37 void inffe_hessian_warning () {
38 libmesh_warning("Warning: Second derivatives for Infinite elements"
39 << " are not yet implemented!"
40 << std::endl);
41 }
42 }
43
44
45 namespace libMesh
46 {
47
48
49
50 // Constructor
51 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
InfFE(const FEType & fet)52 InfFE<Dim,T_radial,T_map>::InfFE (const FEType & fet) :
53 FEBase (Dim, fet),
54
55 _n_total_approx_sf (0),
56 _n_total_qp (0),
57
58 // initialize the current_fe_type to all the same
59 // values as \p fet (since the FE families and coordinate
60 // map type should not change), but use an invalid order
61 // for the radial part (since this is the only order
62 // that may change!).
63 // the data structures like \p phi etc are not initialized
64 // through the constructor, but through reinit()
65 current_fe_type (FEType(fet.order,
66 fet.family,
67 INVALID_ORDER,
68 fet.radial_family,
69 fet.inf_map))
70
71 {
72 // Sanity checks
73 libmesh_assert_equal_to (T_radial, fe_type.radial_family);
74 libmesh_assert_equal_to (T_map, fe_type.inf_map);
75
76 // build the base_fe object
77 if (Dim != 1)
78 base_fe = FEBase::build(Dim-1, fet);
79 }
80
81
82
83 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
attach_quadrature_rule(QBase * q)84 void InfFE<Dim,T_radial,T_map>::attach_quadrature_rule (QBase * q)
85 {
86 libmesh_assert(q);
87 libmesh_assert(base_fe);
88
89 const Order base_int_order = q->get_order();
90 const Order radial_int_order = static_cast<Order>(2 * (static_cast<unsigned int>(fe_type.radial_order.get_order()) + 1) +2);
91 const unsigned int qrule_dim = q->get_dim();
92
93 if (Dim != 1)
94 {
95 // build a Dim-1 quadrature rule of the type that we received
96 base_qrule = QBase::build(q->type(), qrule_dim-1, base_int_order);
97 base_fe->attach_quadrature_rule(base_qrule.get());
98 }
99
100 // in radial direction, always use Gauss quadrature
101 radial_qrule = libmesh_make_unique<QGauss>(1, radial_int_order);
102
103 // Maybe helpful to store the QBase *
104 // with which we initialized our own quadrature rules.
105 // Used e.g. in \p InfFE::reinit(elem,side)
106 qrule = q;
107 }
108
109
110
111
112
113 template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
update_base_elem(const Elem * inf_elem)114 void InfFE<Dim,T_radial,T_base>::update_base_elem (const Elem * inf_elem)
115 {
116 base_elem.reset(InfFEBase::build_elem(inf_elem));
117 }
118
119
120
121
122
123
124 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
reinit(const Elem * inf_elem,const std::vector<Point> * const pts,const std::vector<Real> * const weights)125 void InfFE<Dim,T_radial,T_map>::reinit(const Elem * inf_elem,
126 const std::vector<Point> * const pts,
127 const std::vector<Real> * const weights)
128 {
129 libmesh_assert(base_fe.get());
130 libmesh_assert(inf_elem);
131
132 // I don't understand infinite elements well enough to risk
133 // calculating too little. :-( RHS
134 this->calculate_phi = this->calculate_dphi = this->calculate_dphiref = true;
135 this->get_xyz();
136 this->determine_calculations();
137 base_fe->calculate_phi = base_fe->calculate_dphi = base_fe->calculate_dphiref = true;
138 base_fe->get_xyz();
139 base_fe->determine_calculations();
140
141 if (pts == nullptr)
142 {
143 libmesh_assert(base_fe->qrule);
144 libmesh_assert_equal_to (base_fe->qrule, base_qrule.get());
145 libmesh_assert(radial_qrule.get());
146
147 bool init_shape_functions_required = false;
148
149 // init the radial data fields only when the radial order changes
150 if (current_fe_type.radial_order != fe_type.radial_order)
151 {
152 current_fe_type.radial_order = fe_type.radial_order;
153
154 // Watch out: this call to QBase->init() only works for
155 // current_fe_type = const! To allow variable Order,
156 // the init() of QBase has to be modified...
157 radial_qrule->init(EDGE2);
158
159 // initialize the radial shape functions
160 this->init_radial_shape_functions(inf_elem);
161
162 init_shape_functions_required=true;
163 }
164
165
166 bool update_base_elem_required=true;
167
168 // update the type in accordance to the current cell
169 // and reinit if the cell type has changed or (as in
170 // the case of the hierarchics) the shape functions
171 // depend on the particular element and need a reinit
172 if ((Dim != 1) &&
173 ((this->get_type() != inf_elem->type()) ||
174 (base_fe->shapes_need_reinit())))
175 {
176 // store the new element type, update base_elem
177 // here. Through \p update_base_elem_required,
178 // remember whether it has to be updated (see below).
179 elem_type = inf_elem->type();
180 this->update_base_elem(inf_elem);
181 update_base_elem_required=false;
182
183 // initialize the base quadrature rule for the new element
184 base_qrule->init(base_elem->type());
185
186 // initialize the shape functions in the base
187 base_fe->init_base_shape_functions(base_fe->qrule->get_points(),
188 base_elem.get());
189
190 // compute the shape functions and map functions of base_fe
191 // before using them later in combine_base_radial.
192 base_fe->_fe_map->compute_map (base_fe->dim, base_fe->qrule->get_weights(),
193 base_elem.get(), base_fe->calculate_d2phi);
194 base_fe->compute_shape_functions(base_elem.get(), base_fe->qrule->get_points());
195
196 init_shape_functions_required=true;
197 }
198
199
200 // when either the radial or base part change,
201 // we have to init the whole fields
202 if (init_shape_functions_required)
203 this->init_shape_functions (radial_qrule->get_points(),
204 base_fe->qrule->get_points(),
205 inf_elem);
206
207 // computing the distance only works when we have the current
208 // base_elem stored. This happens when fe_type is const,
209 // the inf_elem->type remains the same. Then we have to
210 // update the base elem _here_.
211 if (update_base_elem_required)
212 this->update_base_elem(inf_elem);
213
214 // compute dist (depends on geometry, therefore has to be updated for
215 // each and every new element), throw radial and base part together
216 this->combine_base_radial (inf_elem);
217
218 this->_fe_map->compute_map (this->dim, _total_qrule_weights, inf_elem, this->calculate_d2phi);
219
220 // Compute the shape functions and the derivatives
221 // at all quadrature points.
222 this->compute_shape_functions (inf_elem,base_fe->qrule->get_points());
223 }
224
225 else // if pts != nullptr
226 {
227 // update the elem_type
228 elem_type = inf_elem->type();
229
230 // We'll assume that pts is a tensor product mesh of points.
231 // That will handle the pts.size()==1 case that we care about
232 // right now, and it will generalize a bit, and it won't break
233 // the assumptions elsewhere in InfFE.
234 std::vector<Point> radial_pts;
235 if (pts->size() > 0)
236 {
237 Real radius = (*pts)[0](Dim-1);
238 radial_pts.push_back(Point(radius));
239 unsigned int n_radial_pts=1;
240 unsigned int n_angular_pts=1;
241 for (auto p : IntRange<std::size_t>(1, pts->size()))
242 {
243 radius = (*pts)[p](Dim-1);
244 // check for repetition of radius: The max. allowed distance is somewhat arbitrary
245 // but the given value should not produce false positives...
246 if (abs(radial_pts[p-n_radial_pts*n_angular_pts](0) - radius)< 1e-4)
247 {
248 // should the next one be in the next segment?
249 if (p+1 == n_radial_pts*(n_angular_pts+1))
250 n_angular_pts++;
251 }
252 // if none yet repeated:
253 else if (n_angular_pts == 1)
254 radial_pts.push_back(Point(radius));
255 // if there was repetition but this does not, the assumed
256 // format does not work:
257 else
258 {
259 libmesh_error_msg("We assumed that the "<<pts->size()
260 <<" points are of tensor-product type with "
261 <<n_radial_pts<<" radial points and "
262 <<n_angular_pts<< " angular points."<<std::endl
263 <<"But apparently point "<<p
264 <<" does not fit that scheme: Its radius is "
265 <<radius <<"but should have "
266 <<radial_pts[p-n_radial_pts*n_angular_pts]<<".");
267 }
268 }
269 }
270 const std::size_t radial_pts_size = radial_pts.size();
271 const std::size_t base_pts_size = pts->size() / radial_pts_size;
272 // If we're a tensor product we should have no remainder
273 libmesh_assert_equal_to
274 (base_pts_size * radial_pts_size, pts->size());
275
276
277 std::vector<Point> base_pts;
278 base_pts.reserve(base_pts_size);
279 for (std::size_t p=0, ps=pts->size(); p != ps; p += radial_pts_size)
280 {
281 Point pt = (*pts)[p];
282 pt(Dim-1) = 0;
283 base_pts.push_back(pt);
284 }
285
286 // init radial shapes
287 this->init_radial_shape_functions(inf_elem, &radial_pts);
288
289 // update the base
290 this->update_base_elem(inf_elem);
291
292 // the finite element on the ifem base
293 base_fe = FEBase::build(Dim-1, this->fe_type);
294
295 base_fe->calculate_phi = base_fe->calculate_dphi = base_fe->calculate_dphiref = true;
296 base_fe->get_xyz();
297 base_fe->determine_calculations();
298
299 // init base shapes
300 base_fe->init_base_shape_functions(base_pts,
301 base_elem.get());
302
303 // compute the shape functions and map functions of base_fe
304 // before using them later in combine_base_radial.
305
306 if (weights)
307 {
308 base_fe->_fe_map->compute_map (base_fe->dim, *weights,
309 base_elem.get(), base_fe->calculate_d2phi);
310 }
311 else
312 {
313 std::vector<Real> dummy_weights (pts->size(), 1.);
314 base_fe->_fe_map->compute_map (base_fe->dim, dummy_weights,
315 base_elem.get(), base_fe->calculate_d2phi);
316 }
317
318 base_fe->compute_shape_functions(base_elem.get(), *pts);
319
320 this->init_shape_functions (radial_pts, base_pts, inf_elem);
321
322 // combine the base and radial shapes
323 this->combine_base_radial (inf_elem);
324
325 // weights
326 if (weights != nullptr)
327 {
328 this->_fe_map->compute_map (this->dim, *weights, inf_elem, this->calculate_d2phi);
329 }
330 else
331 {
332 std::vector<Real> dummy_weights (pts->size(), 1.);
333 this->_fe_map->compute_map (this->dim, dummy_weights, inf_elem, this->calculate_d2phi);
334 }
335
336 // finally compute the ifem shapes
337 this->compute_shape_functions (inf_elem,*pts);
338 }
339
340 if (this->calculate_dual)
341 libmesh_not_implemented_msg("Dual shape support for infinite elements is "
342 "not currently implemented");
343 }
344
345
346
347
348
349 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
350 void
351 InfFE<Dim, T_radial, T_map>::
init_radial_shape_functions(const Elem * libmesh_dbg_var (inf_elem),const std::vector<Point> * radial_pts)352 init_radial_shape_functions(const Elem * libmesh_dbg_var(inf_elem),
353 const std::vector<Point> * radial_pts)
354 {
355 libmesh_assert(radial_qrule.get() || radial_pts);
356 libmesh_assert(inf_elem);
357
358 // Start logging the radial shape function initialization
359 LOG_SCOPE("init_radial_shape_functions()", "InfFE");
360
361 // initialize most of the things related to mapping
362
363 // The order to use in the radial map (currently independent of the element type)
364 const Order radial_mapping_order = InfFERadial::mapping_order();
365 const unsigned int n_radial_mapping_shape_functions =
366 InfFERadial::n_dofs(radial_mapping_order);
367
368 // initialize most of the things related to physical approximation
369 const Order radial_approx_order = fe_type.radial_order;
370 const unsigned int n_radial_approx_shape_functions =
371 InfFERadial::n_dofs(radial_approx_order);
372
373 const std::size_t n_radial_qp =
374 radial_pts ? radial_pts->size() : radial_qrule->n_points();
375 const std::vector<Point> & radial_qp =
376 radial_pts ? *radial_pts : radial_qrule->get_points();
377
378 // resize the radial data fields
379
380 // the radial polynomials (eval)
381 mode.resize (n_radial_approx_shape_functions);
382 dmodedv.resize (n_radial_approx_shape_functions);
383
384 // the (1-v)/2 weight
385 som.resize (n_radial_qp);
386 dsomdv.resize (n_radial_qp);
387
388 // the radial map
389 radial_map.resize (n_radial_mapping_shape_functions);
390 dradialdv_map.resize (n_radial_mapping_shape_functions);
391
392
393 for (unsigned int i=0; i<n_radial_mapping_shape_functions; i++)
394 {
395 radial_map[i].resize (n_radial_qp);
396 dradialdv_map[i].resize (n_radial_qp);
397 }
398
399
400 for (unsigned int i=0; i<n_radial_approx_shape_functions; i++)
401 {
402 mode[i].resize (n_radial_qp);
403 dmodedv[i].resize (n_radial_qp);
404 }
405
406
407 // compute scalar values at radial quadrature points
408 for (std::size_t p=0; p<n_radial_qp; p++)
409 {
410 som[p] = InfFERadial::decay (Dim, radial_qp[p](0));
411 dsomdv[p] = InfFERadial::decay_deriv (radial_qp[p](0));
412 }
413
414
415 // evaluate the mode shapes in radial direction at radial quadrature points
416 for (unsigned int i=0; i<n_radial_approx_shape_functions; i++)
417 for (std::size_t p=0; p<n_radial_qp; p++)
418 {
419 mode[i][p] = InfFE<Dim,T_radial,T_map>::eval (radial_qp[p](0), radial_approx_order, i);
420 dmodedv[i][p] = InfFE<Dim,T_radial,T_map>::eval_deriv (radial_qp[p](0), radial_approx_order, i);
421 }
422
423
424 // evaluate the mapping functions in radial direction at radial quadrature points
425 for (unsigned int i=0; i<n_radial_mapping_shape_functions; i++)
426 for (std::size_t p=0; p<n_radial_qp; p++)
427 {
428 radial_map[i][p] = InfFEMap::eval (radial_qp[p](0), radial_mapping_order, i);
429 dradialdv_map[i][p] = InfFEMap::eval_deriv (radial_qp[p](0), radial_mapping_order, i);
430 }
431 }
432
433
434
435
436
437 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
init_shape_functions(const std::vector<Point> & radial_qp,const std::vector<Point> & base_qp,const Elem * inf_elem)438 void InfFE<Dim,T_radial,T_map>::init_shape_functions(const std::vector<Point> & radial_qp,
439 const std::vector<Point> & base_qp,
440 const Elem * inf_elem)
441 {
442 libmesh_assert(inf_elem);
443
444 // Start logging the radial shape function initialization
445 LOG_SCOPE("init_shape_functions()", "InfFE");
446
447 // fast access to some const ints for the radial data
448 const unsigned int n_radial_mapping_sf = cast_int<unsigned int>(radial_map.size());
449 const unsigned int n_radial_approx_sf = cast_int<unsigned int>(mode.size());
450 const unsigned int n_radial_qp = cast_int<unsigned int>(som.size());
451
452
453 // initialize most of the things related to mapping
454
455 // The element type and order to use in the base map
456 const Order base_mapping_order = base_elem->default_order();
457
458 // the number of base shape functions used to construct the map
459 // (Lagrange shape functions are used for mapping in the base)
460 unsigned int n_base_mapping_shape_functions =
461 InfFEBase::n_base_mapping_sf(*base_elem,
462 base_mapping_order);
463
464 const unsigned int n_total_mapping_shape_functions =
465 n_radial_mapping_sf * n_base_mapping_shape_functions;
466
467 // initialize most of the things related to physical approximation
468 unsigned int n_base_approx_shape_functions;
469 if (Dim > 1)
470 n_base_approx_shape_functions = base_fe->n_shape_functions();
471 else
472 n_base_approx_shape_functions = 1;
473
474
475 const unsigned int n_total_approx_shape_functions =
476 n_radial_approx_sf * n_base_approx_shape_functions;
477
478 // update class member field
479 _n_total_approx_sf = n_total_approx_shape_functions;
480
481
482 // The number of the base quadrature points.
483 const unsigned int n_base_qp = cast_int<unsigned int>(base_qp.size());
484
485 // The total number of quadrature points.
486 const unsigned int n_total_qp = n_radial_qp * n_base_qp;
487
488
489 // update class member field
490 _n_total_qp = n_total_qp;
491
492
493
494 // initialize the node and shape numbering maps
495 {
496 // these vectors work as follows: the i-th entry stores
497 // the associated base/radial node number
498 _radial_node_index.resize(n_total_mapping_shape_functions);
499 _base_node_index.resize(n_total_mapping_shape_functions);
500
501 // similar for the shapes: the i-th entry stores
502 // the associated base/radial shape number
503 _radial_shape_index.resize(n_total_approx_shape_functions);
504 _base_shape_index.resize(n_total_approx_shape_functions);
505
506 const ElemType inf_elem_type = inf_elem->type();
507
508 // fill the node index map
509 for (unsigned int n=0; n<n_total_mapping_shape_functions; n++)
510 {
511 compute_node_indices (inf_elem_type,
512 n,
513 _base_node_index[n],
514 _radial_node_index[n]);
515 libmesh_assert_less (_base_node_index[n], n_base_mapping_shape_functions);
516 libmesh_assert_less (_radial_node_index[n], n_radial_mapping_sf);
517 }
518
519 // fill the shape index map
520 for (unsigned int n=0; n<n_total_approx_shape_functions; n++)
521 {
522 compute_shape_indices (this->fe_type,
523 inf_elem,
524 n,
525 _base_shape_index[n],
526 _radial_shape_index[n]);
527 libmesh_assert_less (_base_shape_index[n], n_base_approx_shape_functions);
528 libmesh_assert_less (_radial_shape_index[n], n_radial_approx_sf);
529 }
530 }
531
532 // resize the base data fields
533 dist.resize(n_base_mapping_shape_functions);
534
535 // resize the total data fields
536
537 // the phase term varies with xi, eta and zeta(v): store it for _all_ qp
538 //
539 // when computing the phase, we need the base approximations
540 // therefore, initialize the phase here, but evaluate it
541 // in combine_base_radial().
542 //
543 // the weight, though, is only needed at the radial quadrature points, n_radial_qp.
544 // but for a uniform interface to the protected data fields
545 // the weight data field (which are accessible from the outside) are expanded to n_total_qp.
546 weight.resize (n_total_qp);
547 dweightdv.resize (n_total_qp);
548 dweight.resize (n_total_qp);
549
550 dphase.resize (n_total_qp);
551 dphasedxi.resize (n_total_qp);
552 dphasedeta.resize (n_total_qp);
553 dphasedzeta.resize (n_total_qp);
554
555 // this vector contains the integration weights for the combined quadrature rule
556 _total_qrule_weights.resize(n_total_qp);
557
558 // InfFE's data fields phi, dphi, dphidx, phi_map etc hold the _total_
559 // shape and mapping functions, respectively
560 {
561 phi.resize (n_total_approx_shape_functions);
562 dphi.resize (n_total_approx_shape_functions);
563 dphidx.resize (n_total_approx_shape_functions);
564 dphidy.resize (n_total_approx_shape_functions);
565 dphidz.resize (n_total_approx_shape_functions);
566 dphidxi.resize (n_total_approx_shape_functions);
567 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
568 inffe_hessian_warning();
569
570 d2phi.resize (n_total_approx_shape_functions);
571 d2phidx2.resize (n_total_approx_shape_functions);
572 d2phidxdy.resize (n_total_approx_shape_functions);
573 d2phidxdz.resize (n_total_approx_shape_functions);
574 d2phidy2.resize (n_total_approx_shape_functions);
575 d2phidydz.resize (n_total_approx_shape_functions);
576 d2phidz2.resize (n_total_approx_shape_functions);
577 d2phidxi2.resize (n_total_approx_shape_functions);
578
579 if (Dim > 1)
580 {
581 d2phidxideta.resize(n_total_approx_shape_functions);
582 d2phideta2.resize(n_total_approx_shape_functions);
583 }
584
585 if (Dim > 2)
586 {
587 d2phidetadzeta.resize(n_total_approx_shape_functions);
588 d2phidxidzeta.resize(n_total_approx_shape_functions);
589 d2phidzeta2.resize(n_total_approx_shape_functions);
590 }
591 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
592
593 if (Dim > 1)
594 dphideta.resize(n_total_approx_shape_functions);
595
596 if (Dim == 3)
597 dphidzeta.resize(n_total_approx_shape_functions);
598
599 std::vector<std::vector<Real>> & phi_map = this->_fe_map->get_phi_map();
600 std::vector<std::vector<Real>> & dphidxi_map = this->_fe_map->get_dphidxi_map();
601
602 phi_map.resize(n_total_mapping_shape_functions);
603 dphidxi_map.resize(n_total_mapping_shape_functions);
604 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
605 std::vector<std::vector<Real>> & d2phidxi2_map = this->_fe_map->get_d2phidxi2_map();
606 d2phidxi2_map.resize(n_total_mapping_shape_functions);
607
608 if (Dim > 1)
609 {
610 std::vector<std::vector<Real>> & d2phidxideta_map = this->_fe_map->get_d2phidxideta_map();
611 std::vector<std::vector<Real>> & d2phideta2_map = this->_fe_map->get_d2phideta2_map();
612 d2phidxideta_map.resize(n_total_mapping_shape_functions);
613 d2phideta2_map.resize(n_total_mapping_shape_functions);
614 }
615
616 if (Dim == 3)
617 {
618 std::vector<std::vector<Real>> & d2phidxidzeta_map = this->_fe_map->get_d2phidxidzeta_map();
619 std::vector<std::vector<Real>> & d2phidetadzeta_map = this->_fe_map->get_d2phidetadzeta_map();
620 std::vector<std::vector<Real>> & d2phidzeta2_map = this->_fe_map->get_d2phidzeta2_map();
621 d2phidxidzeta_map.resize(n_total_mapping_shape_functions);
622 d2phidetadzeta_map.resize(n_total_mapping_shape_functions);
623 d2phidzeta2_map.resize(n_total_mapping_shape_functions);
624 }
625 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
626
627 if (Dim > 1)
628 {
629 std::vector<std::vector<Real>> & dphideta_map = this->_fe_map->get_dphideta_map();
630 dphideta_map.resize(n_total_mapping_shape_functions);
631 }
632
633 if (Dim == 3)
634 {
635 std::vector<std::vector<Real>> & dphidzeta_map = this->_fe_map->get_dphidzeta_map();
636 dphidzeta_map.resize(n_total_mapping_shape_functions);
637 }
638 }
639
640 // collect all the for loops, where inner vectors are
641 // resized to the appropriate number of quadrature points
642 {
643 for (unsigned int i=0; i<n_total_approx_shape_functions; i++)
644 {
645 phi[i].resize (n_total_qp);
646 dphi[i].resize (n_total_qp);
647 dphidx[i].resize (n_total_qp);
648 dphidy[i].resize (n_total_qp);
649 dphidz[i].resize (n_total_qp);
650 dphidxi[i].resize (n_total_qp);
651 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
652 d2phi[i].resize (n_total_qp);
653 d2phidx2[i].resize (n_total_qp);
654 d2phidxdy[i].resize (n_total_qp);
655 d2phidxdz[i].resize (n_total_qp);
656 d2phidy2[i].resize (n_total_qp);
657 d2phidydz[i].resize (n_total_qp);
658 d2phidy2[i].resize (n_total_qp);
659 d2phidxi2[i].resize (n_total_qp);
660
661 if (Dim > 1)
662 {
663 d2phidxideta[i].resize (n_total_qp);
664 d2phideta2[i].resize (n_total_qp);
665 }
666 if (Dim > 2)
667 {
668 d2phidxidzeta[i].resize (n_total_qp);
669 d2phidetadzeta[i].resize (n_total_qp);
670 d2phidzeta2[i].resize (n_total_qp);
671 }
672 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
673
674 if (Dim > 1)
675 dphideta[i].resize (n_total_qp);
676
677 if (Dim == 3)
678 dphidzeta[i].resize (n_total_qp);
679
680 }
681
682 for (unsigned int i=0; i<n_total_mapping_shape_functions; i++)
683 {
684 std::vector<std::vector<Real>> & phi_map = this->_fe_map->get_phi_map();
685 std::vector<std::vector<Real>> & dphidxi_map = this->_fe_map->get_dphidxi_map();
686 phi_map[i].resize (n_total_qp);
687 dphidxi_map[i].resize (n_total_qp);
688 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
689 std::vector<std::vector<Real>> & d2phidxi2_map = this->_fe_map->get_d2phidxi2_map();
690 d2phidxi2_map[i].resize (n_total_qp);
691 if (Dim > 1)
692 {
693 std::vector<std::vector<Real>> & d2phidxideta_map = this->_fe_map->get_d2phidxideta_map();
694 std::vector<std::vector<Real>> & d2phideta2_map = this->_fe_map->get_d2phideta2_map();
695 d2phidxideta_map[i].resize (n_total_qp);
696 d2phideta2_map[i].resize (n_total_qp);
697 }
698
699 if (Dim > 2)
700 {
701 std::vector<std::vector<Real>> & d2phidxidzeta_map = this->_fe_map->get_d2phidxidzeta_map();
702 std::vector<std::vector<Real>> & d2phidetadzeta_map = this->_fe_map->get_d2phidetadzeta_map();
703 std::vector<std::vector<Real>> & d2phidzeta2_map = this->_fe_map->get_d2phidzeta2_map();
704 d2phidxidzeta_map[i].resize (n_total_qp);
705 d2phidetadzeta_map[i].resize (n_total_qp);
706 d2phidzeta2_map[i].resize (n_total_qp);
707 }
708 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
709
710 if (Dim > 1)
711 {
712 std::vector<std::vector<Real>> & dphideta_map = this->_fe_map->get_dphideta_map();
713 dphideta_map[i].resize (n_total_qp);
714 }
715
716 if (Dim == 3)
717 {
718 std::vector<std::vector<Real>> & dphidzeta_map = this->_fe_map->get_dphidzeta_map();
719 dphidzeta_map[i].resize (n_total_qp);
720 }
721 }
722 }
723
724
725
726 {
727 // (a) compute scalar values at _all_ quadrature points -- for uniform
728 // access from the outside to these fields
729 // (b) form a std::vector<Real> which contains the appropriate weights
730 // of the combined quadrature rule!
731 libmesh_assert_equal_to (radial_qp.size(), n_radial_qp);
732
733 if (radial_qrule && base_qrule)
734 {
735 const std::vector<Real> & radial_qw = radial_qrule->get_weights();
736 const std::vector<Real> & base_qw = base_qrule->get_weights();
737 libmesh_assert_equal_to (radial_qw.size(), n_radial_qp);
738 libmesh_assert_equal_to (base_qw.size(), n_base_qp);
739
740 for (unsigned int rp=0; rp<n_radial_qp; rp++)
741 for (unsigned int bp=0; bp<n_base_qp; bp++)
742 {
743 weight[bp + rp*n_base_qp] = InfFERadial::D(radial_qp[rp](0));
744 dweightdv[bp + rp*n_base_qp] = InfFERadial::D_deriv(radial_qp[rp](0));
745 _total_qrule_weights[bp + rp*n_base_qp] = radial_qw[rp] * base_qw[bp];
746 }
747 }
748 else
749 {
750 for (unsigned int rp=0; rp<n_radial_qp; rp++)
751 for (unsigned int bp=0; bp<n_base_qp; bp++)
752 {
753 weight[bp + rp*n_base_qp] = InfFERadial::D(radial_qp[rp](0));
754 dweightdv[bp + rp*n_base_qp] = InfFERadial::D_deriv(radial_qp[rp](0));
755 }
756 }
757 }
758 }
759
760
761
762
763
764 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
combine_base_radial(const Elem * inf_elem)765 void InfFE<Dim,T_radial,T_map>::combine_base_radial(const Elem * inf_elem)
766 {
767 libmesh_assert(inf_elem);
768 // at least check whether the base element type is correct.
769 // otherwise this version of computing dist would give problems
770 libmesh_assert_equal_to (base_elem->type(),
771 InfFEBase::get_elem_type(inf_elem->type()));
772
773 // Start logging the combination of radial and base parts
774 LOG_SCOPE("combine_base_radial()", "InfFE");
775
776 // zero the phase, since it is to be summed up
777 std::fill (dphasedxi.begin(), dphasedxi.end(), 0.);
778 std::fill (dphasedeta.begin(), dphasedeta.end(), 0.);
779 std::fill (dphasedzeta.begin(), dphasedzeta.end(), 0.);
780
781
782 const unsigned int n_base_mapping_sf = cast_int<unsigned int>(dist.size());
783 const Point origin = inf_elem->origin();
784
785 // for each new infinite element, compute the radial distances
786 for (unsigned int n=0; n<n_base_mapping_sf; n++)
787 dist[n] = Point(base_elem->point(n) - origin).norm();
788
789
790 switch (Dim)
791 {
792 // 1D
793 case 1:
794 {
795 libmesh_not_implemented();
796 break;
797 }
798
799 // 2D
800 case 2:
801 {
802 libmesh_not_implemented();
803 break;
804 }
805
806 // 3D
807 case 3:
808 {
809 // fast access to the approximation and mapping shapes of base_fe
810 const std::vector<std::vector<Real>> & S = base_fe->phi;
811 const std::vector<std::vector<Real>> & Ss = base_fe->dphidxi;
812 const std::vector<std::vector<Real>> & St = base_fe->dphideta;
813 const std::vector<std::vector<Real>> & S_map = (base_fe->get_fe_map()).get_phi_map();
814 const std::vector<std::vector<Real>> & Ss_map = (base_fe->get_fe_map()).get_dphidxi_map();
815 const std::vector<std::vector<Real>> & St_map = (base_fe->get_fe_map()).get_dphideta_map();
816
817 const unsigned int n_radial_qp = cast_int<unsigned int>(som.size());
818 if (radial_qrule)
819 libmesh_assert_equal_to(n_radial_qp, radial_qrule->n_points());
820 const unsigned int n_base_qp = cast_int<unsigned int>(S_map[0].size());
821 if (base_qrule)
822 libmesh_assert_equal_to(n_base_qp, base_qrule->n_points());
823
824 const unsigned int n_total_mapping_sf =
825 cast_int<unsigned int>(radial_map.size()) * n_base_mapping_sf;
826
827 const unsigned int n_total_approx_sf =
828 InfFERadial::n_dofs(fe_type.radial_order) *
829 base_fe->n_shape_functions();
830
831
832 // compute the phase term derivatives
833 {
834 unsigned int tp=0;
835 for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qps
836 for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qps
837 {
838 // sum over all base shapes, to get the average distance
839 for (unsigned int i=0; i<n_base_mapping_sf; i++)
840 {
841 dphasedxi[tp] += Ss_map[i][bp] * dist[i] * radial_map [1][rp];
842 dphasedeta[tp] += St_map[i][bp] * dist[i] * radial_map [1][rp];
843 dphasedzeta[tp] += S_map [i][bp] * dist[i] * dradialdv_map[1][rp];
844 }
845
846 tp++;
847
848 } // loop radial and base qps
849 }
850
851 libmesh_assert_equal_to (phi.size(), n_total_approx_sf);
852 libmesh_assert_equal_to (dphidxi.size(), n_total_approx_sf);
853 libmesh_assert_equal_to (dphideta.size(), n_total_approx_sf);
854 libmesh_assert_equal_to (dphidzeta.size(), n_total_approx_sf);
855
856 // compute the overall approximation shape functions,
857 // pick the appropriate radial and base shapes through using
858 // _base_shape_index and _radial_shape_index
859 for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qps
860 for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qps
861 for (unsigned int ti=0; ti<n_total_approx_sf; ti++) // over _all_ approx_sf
862 {
863 // let the index vectors take care of selecting the appropriate base/radial shape
864 const unsigned int bi = _base_shape_index [ti];
865 const unsigned int ri = _radial_shape_index[ti];
866 phi [ti][bp+rp*n_base_qp] = S [bi][bp] * mode[ri][rp] * som[rp];
867 dphidxi [ti][bp+rp*n_base_qp] = Ss[bi][bp] * mode[ri][rp] * som[rp];
868 dphideta [ti][bp+rp*n_base_qp] = St[bi][bp] * mode[ri][rp] * som[rp];
869 dphidzeta[ti][bp+rp*n_base_qp] = S [bi][bp]
870 * (dmodedv[ri][rp] * som[rp] + mode[ri][rp] * dsomdv[rp]);
871 }
872
873 std::vector<std::vector<Real>> & phi_map = this->_fe_map->get_phi_map();
874 std::vector<std::vector<Real>> & dphidxi_map = this->_fe_map->get_dphidxi_map();
875 std::vector<std::vector<Real>> & dphideta_map = this->_fe_map->get_dphideta_map();
876 std::vector<std::vector<Real>> & dphidzeta_map = this->_fe_map->get_dphidzeta_map();
877
878 libmesh_assert_equal_to (phi_map.size(), n_total_mapping_sf);
879 libmesh_assert_equal_to (dphidxi_map.size(), n_total_mapping_sf);
880 libmesh_assert_equal_to (dphideta_map.size(), n_total_mapping_sf);
881 libmesh_assert_equal_to (dphidzeta_map.size(), n_total_mapping_sf);
882
883 // compute the overall mapping functions,
884 // pick the appropriate radial and base entries through using
885 // _base_node_index and _radial_node_index
886 for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qps
887 for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qps
888 for (unsigned int ti=0; ti<n_total_mapping_sf; ti++) // over all mapping shapes
889 {
890 // let the index vectors take care of selecting the appropriate base/radial mapping shape
891 const unsigned int bi = _base_node_index [ti];
892 const unsigned int ri = _radial_node_index[ti];
893 phi_map [ti][bp+rp*n_base_qp] = S_map [bi][bp] * radial_map [ri][rp];
894 dphidxi_map [ti][bp+rp*n_base_qp] = Ss_map[bi][bp] * radial_map [ri][rp];
895 dphideta_map [ti][bp+rp*n_base_qp] = St_map[bi][bp] * radial_map [ri][rp];
896 dphidzeta_map[ti][bp+rp*n_base_qp] = S_map [bi][bp] * dradialdv_map[ri][rp];
897 }
898
899 break;
900 }
901
902 default:
903 libmesh_error_msg("Unsupported Dim = " << Dim);
904 }
905 }
906
907
908
909
910
911
912 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
compute_shape_functions(const Elem *,const std::vector<Point> &)913 void InfFE<Dim,T_radial,T_map>::compute_shape_functions(const Elem *,
914 const std::vector<Point> &)
915 {
916 // Start logging the overall computation of shape functions
917 LOG_SCOPE("compute_shape_functions()", "InfFE");
918
919 const unsigned int n_total_qp = _n_total_qp;
920
921 // Compute the shape function values (and derivatives)
922 // at the Quadrature points. Note that the actual values
923 // have already been computed via init_shape_functions
924
925 // Compute the value of the derivative shape function i at quadrature point p
926 switch (dim)
927 {
928
929 case 1:
930 {
931 libmesh_not_implemented();
932 break;
933 }
934
935 case 2:
936 {
937 libmesh_not_implemented();
938 break;
939 }
940
941 case 3:
942 {
943 const std::vector<Real> & dxidx_map = this->_fe_map->get_dxidx();
944 const std::vector<Real> & dxidy_map = this->_fe_map->get_dxidy();
945 const std::vector<Real> & dxidz_map = this->_fe_map->get_dxidz();
946
947 const std::vector<Real> & detadx_map = this->_fe_map->get_detadx();
948 const std::vector<Real> & detady_map = this->_fe_map->get_detady();
949 const std::vector<Real> & detadz_map = this->_fe_map->get_detadz();
950
951 const std::vector<Real> & dzetadx_map = this->_fe_map->get_dzetadx();
952 const std::vector<Real> & dzetady_map = this->_fe_map->get_dzetady();
953 const std::vector<Real> & dzetadz_map = this->_fe_map->get_dzetadz();
954
955 // These are _all_ shape functions of this infinite element
956 for (auto i : index_range(phi))
957 for (unsigned int p=0; p<n_total_qp; p++)
958 {
959 // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
960 dphi[i][p](0) =
961 dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
962 dphideta[i][p]*detadx_map[p] +
963 dphidzeta[i][p]*dzetadx_map[p]);
964
965 // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
966 dphi[i][p](1) =
967 dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
968 dphideta[i][p]*detady_map[p] +
969 dphidzeta[i][p]*dzetady_map[p]);
970
971 // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
972 dphi[i][p](2) =
973 dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
974 dphideta[i][p]*detadz_map[p] +
975 dphidzeta[i][p]*dzetadz_map[p]);
976 }
977
978
979 // This is the derivative of the phase term of this infinite element
980 for (unsigned int p=0; p<n_total_qp; p++)
981 {
982 // the derivative of the phase term
983 dphase[p](0) = (dphasedxi[p] * dxidx_map[p] +
984 dphasedeta[p] * detadx_map[p] +
985 dphasedzeta[p] * dzetadx_map[p]);
986
987 dphase[p](1) = (dphasedxi[p] * dxidy_map[p] +
988 dphasedeta[p] * detady_map[p] +
989 dphasedzeta[p] * dzetady_map[p]);
990
991 dphase[p](2) = (dphasedxi[p] * dxidz_map[p] +
992 dphasedeta[p] * detadz_map[p] +
993 dphasedzeta[p] * dzetadz_map[p]);
994
995 // the derivative of the radial weight - varies only in radial direction,
996 // therefore dweightdxi = dweightdeta = 0.
997 dweight[p](0) = dweightdv[p] * dzetadx_map[p];
998
999 dweight[p](1) = dweightdv[p] * dzetady_map[p];
1000
1001 dweight[p](2) = dweightdv[p] * dzetadz_map[p];
1002 }
1003
1004 break;
1005 }
1006
1007 default:
1008 libmesh_error_msg("Unsupported dim = " << dim);
1009 }
1010 }
1011
1012
1013
1014 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
shapes_need_reinit()1015 bool InfFE<Dim,T_radial,T_map>::shapes_need_reinit() const
1016 {
1017 return false;
1018 }
1019
1020 } // namespace libMesh
1021
1022
1023 // Explicit instantiations
1024 #include "libmesh/inf_fe_instantiate_1D.h"
1025 #include "libmesh/inf_fe_instantiate_2D.h"
1026 #include "libmesh/inf_fe_instantiate_3D.h"
1027
1028
1029
1030 #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1031