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/elem.h"
22 : #include "libmesh/fe.h"
23 : #include "libmesh/fe_interface.h"
24 : #include "libmesh/fe_macro.h"
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/quadrature.h"
27 : #include "libmesh/tensor_value.h"
28 : #include "libmesh/enum_elem_type.h"
29 : #include "libmesh/quadrature_gauss.h"
30 : #include "libmesh/libmesh_singleton.h"
31 :
32 : namespace {
33 : // Put this outside a templated class, so we only get 1 warning
34 : // during our unit tests, not 1 warning for each of the zillion FE
35 : // specializations we test.
36 417 : void nonlagrange_dual_warning () {
37 : libmesh_warning("dual calculations have only been verified for the LAGRANGE family");
38 417 : }
39 : }
40 :
41 :
42 : namespace libMesh
43 : {
44 : // ------------------------------------------------------------
45 : // Whether we cache the node locations, edge and face orientations of the last
46 : // element we computed on as needed to avoid calling init_shape_functions and
47 : // compute_shape_functions
48 : static const bool * caching = nullptr;
49 :
50 : class CachingSetup: public Singleton::Setup
51 : {
52 : private:
53 2119 : void setup() { caching = new bool(!on_command_line("--disable-caching")); }
54 : public:
55 1645 : ~CachingSetup() { delete caching; caching = nullptr; }
56 : } caching_setup;
57 :
58 :
59 : // ------------------------------------------------------------
60 : // FE class members
61 : template <unsigned int Dim, FEFamily T>
62 2619946 : FE<Dim,T>::FE (const FEType & fet) :
63 : FEGenericBase<typename FEOutputType<T>::type> (Dim,fet),
64 1175500 : last_side(INVALID_ELEM),
65 2619946 : last_edge(INVALID_ELEM)
66 : {
67 : // Sanity check. Make sure the
68 : // Family specified in the template instantiation
69 : // matches the one in the FEType object
70 725253 : libmesh_assert_equal_to (T, this->get_family());
71 2619946 : }
72 :
73 :
74 : template <unsigned int Dim, FEFamily T>
75 0 : unsigned int FE<Dim,T>::n_shape_functions () const
76 : {
77 0 : if (this->_elem)
78 0 : return this->n_dofs (this->_elem,
79 0 : this->fe_type.order + this->_p_level);
80 :
81 0 : return this->n_dofs (this->get_type(),
82 0 : this->fe_type.order + this->_p_level);
83 : }
84 :
85 :
86 : template <unsigned int Dim, FEFamily T>
87 1454742 : void FE<Dim,T>::attach_quadrature_rule (QBase * q)
88 : {
89 424952 : libmesh_assert(q);
90 1454742 : this->qrule = q;
91 : // make sure we don't cache results from a previous quadrature rule
92 1454742 : this->_elem = nullptr;
93 1454742 : this->_elem_type = INVALID_ELEM;
94 1454742 : return;
95 : }
96 :
97 :
98 : template <unsigned int Dim, FEFamily T>
99 272559 : void FE<Dim,T>::dofs_on_side(const Elem * const elem,
100 : const Order o,
101 : unsigned int s,
102 : std::vector<unsigned int> & di,
103 : const bool add_p_level)
104 : {
105 78602 : libmesh_assert(elem);
106 78602 : libmesh_assert_less (s, elem->n_sides());
107 :
108 78602 : di.clear();
109 78602 : unsigned int nodenum = 0;
110 272559 : const unsigned int n_nodes = elem->n_nodes();
111 3202189 : for (unsigned int n = 0; n != n_nodes; ++n)
112 : {
113 : const unsigned int n_dofs =
114 4552612 : n_dofs_at_node(*elem, static_cast<Order>(o + add_p_level*elem->p_level()), n);
115 2929630 : if (elem->is_node_on_side(n, s))
116 2808816 : for (unsigned int i = 0; i != n_dofs; ++i)
117 1587843 : di.push_back(nodenum++);
118 : else
119 1708657 : nodenum += n_dofs;
120 : }
121 272559 : }
122 :
123 :
124 :
125 : template <unsigned int Dim, FEFamily T>
126 41868 : void FE<Dim,T>::dofs_on_edge(const Elem * const elem,
127 : const Order o,
128 : unsigned int e,
129 : std::vector<unsigned int> & di,
130 : const bool add_p_level)
131 : {
132 10789 : libmesh_assert(elem);
133 10789 : libmesh_assert_less (e, elem->n_edges());
134 :
135 10789 : di.clear();
136 10789 : unsigned int nodenum = 0;
137 41868 : const unsigned int n_nodes = elem->n_nodes();
138 861904 : for (unsigned int n = 0; n != n_nodes; ++n)
139 : {
140 : const unsigned int n_dofs =
141 1240500 : n_dofs_at_node(*elem, static_cast<Order>(o + add_p_level*elem->p_level()), n);
142 820036 : if (elem->is_node_on_edge(n, e))
143 286038 : for (unsigned int i = 0; i != n_dofs; ++i)
144 162750 : di.push_back(nodenum++);
145 : else
146 696748 : nodenum += n_dofs;
147 : }
148 41868 : }
149 :
150 :
151 :
152 : template <unsigned int Dim, FEFamily T>
153 1479690 : void FE<Dim,T>::cache(const Elem * elem)
154 : {
155 1479690 : cached_nodes.resize(elem->n_nodes());
156 17669320 : for (auto n : elem->node_index_range())
157 24077328 : cached_nodes[n] = elem->point(n);
158 :
159 1479690 : if (FEInterface::orientation_dependent(T))
160 : {
161 1184025 : cached_edges.resize(elem->n_edges());
162 7017660 : for (auto n : elem->edge_index_range())
163 5833635 : cached_edges[n] = elem->positive_edge_orientation(n);
164 :
165 1184025 : cached_faces.resize(elem->n_faces());
166 3637893 : for (auto n : elem->face_index_range())
167 2453868 : cached_faces[n] = elem->positive_face_orientation(n);
168 : }
169 1479690 : }
170 :
171 :
172 :
173 : template <unsigned int Dim, FEFamily T>
174 32666967 : bool FE<Dim,T>::matches_cache(const Elem * elem)
175 : {
176 42168858 : bool m = cached_nodes.size() == elem->n_nodes();
177 45864051 : for (unsigned n = 1; m && n < elem->n_nodes(); n++)
178 16750466 : m = (elem->point(n) - elem->point(0)).relative_fuzzy_equals(cached_nodes[n] - cached_nodes[0]);
179 :
180 32666967 : if (FEInterface::orientation_dependent(T))
181 : {
182 1216260 : m &= cached_edges.size() == elem->n_edges();
183 2188689 : for (unsigned n = 0; m && n < elem->n_edges(); n++)
184 972429 : m = elem->positive_edge_orientation(n) == cached_edges[n];
185 :
186 1216260 : m &= cached_faces.size() == elem->n_faces();
187 1324128 : for (unsigned n = 0; m && n < elem->n_faces(); n++)
188 107868 : m = elem->positive_face_orientation(n) == cached_faces[n];
189 : }
190 :
191 32666967 : return m;
192 : }
193 :
194 :
195 :
196 : template <unsigned int Dim, FEFamily T>
197 48942028 : void FE<Dim,T>::reinit(const Elem * elem,
198 : const std::vector<Point> * const pts,
199 : const std::vector<Real> * const weights)
200 : {
201 : // We can be called with no element. If we're evaluating SCALAR
202 : // dofs we'll still have work to do.
203 : // libmesh_assert(elem);
204 :
205 : // We're calculating now! Time to determine what.
206 48942028 : this->determine_calculations();
207 :
208 : // Try to avoid calling init_shape_functions
209 : // even when shapes_need_reinit
210 14164176 : bool cached_elem_still_fits = false;
211 :
212 : // Most of the hard work happens when we have an actual element
213 48942028 : if (elem)
214 : {
215 : // Initialize the shape functions at the user-specified
216 : // points
217 48942028 : if (pts != nullptr)
218 : {
219 : // Set the type and p level for this element
220 15822308 : this->_elem = elem;
221 15822308 : this->_elem_type = elem->type();
222 15822308 : this->_elem_p_level = elem->p_level();
223 20089150 : this->_p_level = this->_add_p_level_in_reinit * elem->p_level();
224 :
225 : // Initialize the shape functions
226 8533684 : this->_fe_map->template init_reference_to_physical_map<Dim>
227 7288624 : (*pts, elem);
228 15822308 : this->init_shape_functions (*pts, elem);
229 :
230 : // The shape functions do not correspond to the qrule
231 15822308 : this->shapes_on_quadrature = false;
232 : }
233 :
234 : // If there are no user specified points, we use the
235 : // quadrature rule
236 :
237 : // update the type in accordance to the current cell
238 : // and reinit if the cell type has changed or (as in
239 : // the case of the hierarchics) the shape functions need
240 : // reinit, since they depend on the particular element shape
241 : else
242 : {
243 9897334 : libmesh_assert(this->qrule);
244 33119720 : this->qrule->init(*elem);
245 :
246 33119720 : if (this->qrule->shapes_need_reinit())
247 0 : this->shapes_on_quadrature = false;
248 :
249 : // We're not going to bother trying to cache nodal
250 : // points *and* weights for fancier mapping types.
251 56116123 : if (this->get_type() != elem->type() ||
252 32801269 : (elem->runtime_topology() &&
253 22996403 : this->_elem != elem) ||
254 32801269 : this->_elem_p_level != elem->p_level() ||
255 75723294 : !this->shapes_on_quadrature ||
256 19339882 : elem->mapping_type() != LAGRANGE_MAP)
257 : {
258 : // Set the type and p level for this element
259 452753 : this->_elem = elem;
260 452753 : this->_elem_type = elem->type();
261 452753 : this->_elem_p_level = elem->p_level();
262 576166 : this->_p_level = this->_add_p_level_in_reinit * elem->p_level();
263 :
264 : // Initialize the shape functions
265 382371 : this->_fe_map->template init_reference_to_physical_map<Dim>
266 452753 : (this->qrule->get_points(), elem);
267 576166 : this->init_shape_functions (this->qrule->get_points(), elem);
268 : }
269 : else
270 : {
271 32666967 : this->_elem = elem;
272 :
273 : // Check if cached element's nodes, edge and face orientations still fit
274 32666967 : cached_elem_still_fits = this->matches_cache(elem);
275 :
276 : // Initialize the shape functions if needed
277 32666967 : if (this->shapes_need_reinit() && !cached_elem_still_fits)
278 : {
279 916023 : this->_fe_map->template init_reference_to_physical_map<Dim>
280 1248490 : (this->qrule->get_points(), elem);
281 1553831 : this->init_shape_functions (this->qrule->get_points(), elem);
282 : }
283 : }
284 :
285 : // Replace cached nodes, edge and face orientations if no longer fitting
286 33119720 : if (this->shapes_need_reinit() && !cached_elem_still_fits && *caching)
287 1479690 : this->cache(elem);
288 :
289 : // The shape functions correspond to the qrule
290 33119720 : this->shapes_on_quadrature = true;
291 : }
292 : }
293 : else // With no defined elem, so mapping or caching to
294 : // be done, and our "quadrature rule" is one point for nonlocal
295 : // (SCALAR) variables and zero points for local variables.
296 : {
297 0 : this->_elem = nullptr;
298 0 : this->_elem_type = INVALID_ELEM;
299 0 : this->_elem_p_level = 0;
300 0 : this->_p_level = 0;
301 :
302 0 : if (!pts)
303 : {
304 : if (T == SCALAR)
305 : {
306 0 : this->qrule->get_points() =
307 0 : std::vector<Point>(1,Point(0));
308 :
309 0 : this->qrule->get_weights() =
310 0 : std::vector<Real>(1,1);
311 : }
312 : else
313 : {
314 0 : this->qrule->get_points().clear();
315 0 : this->qrule->get_weights().clear();
316 : }
317 :
318 0 : this->init_shape_functions (this->qrule->get_points(), elem);
319 : }
320 : else
321 0 : this->init_shape_functions (*pts, elem);
322 : }
323 :
324 : // Compute the map for this element.
325 48942028 : if (pts != nullptr)
326 : {
327 15822308 : if (weights != nullptr)
328 : {
329 42776 : this->_fe_map->compute_map (this->dim, *weights, elem, this->calculate_d2phi);
330 : }
331 : else
332 : {
333 24272884 : std::vector<Real> dummy_weights (pts->size(), 1.);
334 15779532 : this->_fe_map->compute_map (this->dim, dummy_weights, elem, this->calculate_d2phi);
335 : }
336 : }
337 : else
338 : {
339 42745024 : this->_fe_map->compute_map (this->dim, this->qrule->get_weights(), elem, this->calculate_d2phi);
340 : }
341 :
342 : // Compute the shape functions and the derivatives at all of the
343 : // quadrature points.
344 48942028 : if (!cached_elem_still_fits)
345 : {
346 38877064 : if (pts != nullptr)
347 15822308 : this->compute_shape_functions (elem,*pts);
348 : else
349 29978083 : this->compute_shape_functions(elem,this->qrule->get_points());
350 38877064 : if (this->calculate_dual)
351 : {
352 : if (T != LAGRANGE)
353 417 : nonlagrange_dual_warning();
354 : // Check if we need to calculate the dual coefficients based on the default QRule
355 : // We keep the default dual coeff calculation for the initial stage of the simulation
356 : // and in the middel of the simulation when a customized QRule is not provided.
357 : // This is used in MOOSE mortar-based contact. Currently, we re-compute dual_coeff
358 : // for all the elements on the mortar segment mesh by setting `calculate_default_dual_coeff' = false
359 : // in MOOSE (in `Assembly::reinitDual`) and use the customized QRule for calculating the dual shape coefficients
360 : // This is to be improved in the future
361 1753 : if (elem && this->calculate_default_dual_coeff)
362 1753 : this->reinit_default_dual_shape_coeffs(elem);
363 : // The dual shape functions relies on the customized shape functions
364 : // and the coefficient matrix, \p dual_coeff
365 1753 : this->compute_dual_shape_functions();
366 : }
367 : }
368 48942028 : }
369 :
370 : template <unsigned int Dim, FEFamily T>
371 1753 : void FE<Dim,T>::reinit_dual_shape_coeffs(const Elem * elem,
372 : const std::vector<Point> & pts,
373 : const std::vector<Real> & JxW)
374 : {
375 : // Set the type and p level for this element
376 1753 : this->_elem = elem;
377 1753 : this->_elem_type = elem->type();
378 1753 : this->_elem_p_level = elem->p_level();
379 2230 : this->_p_level = this->_add_p_level_in_reinit * elem->p_level();
380 :
381 954 : const unsigned int n_shapes =
382 799 : this->n_dofs(elem, this->get_order());
383 :
384 1431 : std::vector<std::vector<OutputShape>> phi_vals;
385 1753 : phi_vals.resize(n_shapes);
386 30942 : for (const auto i : make_range(phi_vals.size()))
387 45049 : phi_vals[i].resize(pts.size());
388 :
389 1753 : all_shapes(elem, this->get_order(), pts, phi_vals);
390 1753 : this->compute_dual_shape_coeffs(JxW, phi_vals);
391 1753 : }
392 :
393 : template <unsigned int Dim, FEFamily T>
394 1753 : void FE<Dim,T>::reinit_default_dual_shape_coeffs (const Elem * elem)
395 : {
396 477 : libmesh_assert(elem);
397 :
398 477 : FEType default_fe_type(this->get_order(), T);
399 2230 : QGauss default_qrule(elem->dim(), default_fe_type.default_quadrature_order());
400 1753 : default_qrule.init(*elem);
401 : // In preparation of computing dual_coeff, we compute the default shape
402 : // function values and use these to compute the dual shape coefficients.
403 : // The TRUE dual_phi values are computed in compute_dual_shape_functions()
404 1753 : this->reinit_dual_shape_coeffs(elem, default_qrule.get_points(), default_qrule.get_weights());
405 : // we do not compute default dual coeff many times as this can be expensive
406 477 : this->set_calculate_default_dual_coeff(false);
407 1753 : }
408 :
409 :
410 : template <unsigned int Dim, FEFamily T>
411 1753 : void FE<Dim,T>::init_dual_shape_functions(const unsigned int n_shapes, const unsigned int n_qp)
412 : {
413 1753 : if (!this->calculate_dual)
414 0 : return;
415 :
416 477 : libmesh_assert_msg(this->calculate_phi,
417 : "dual shape function calculation relies on "
418 : "primal shape functions being calculated");
419 :
420 1753 : this->dual_phi.resize(n_shapes);
421 1753 : if (this->calculate_dphi)
422 1713 : this->dual_dphi.resize(n_shapes);
423 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
424 1753 : if (this->calculate_d2phi)
425 1632 : this->dual_d2phi.resize(n_shapes);
426 : #endif
427 :
428 30942 : for (auto i : index_range(this->dual_phi))
429 : {
430 37119 : this->dual_phi[i].resize(n_qp);
431 29189 : if (this->calculate_dphi)
432 36813 : this->dual_dphi[i].resize(n_qp);
433 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
434 29189 : if (this->calculate_d2phi)
435 35888 : this->dual_d2phi[i].resize(n_qp);
436 : #endif
437 : }
438 : }
439 :
440 : template <unsigned int Dim, FEFamily T>
441 17026223 : void FE<Dim,T>::init_shape_functions(const std::vector<Point> & qp,
442 : const Elem * elem)
443 : {
444 : // Start logging the shape function initialization
445 13714395 : LOG_SCOPE("init_shape_functions()", "FE");
446 :
447 : // The number of quadrature points.
448 9140908 : const unsigned int n_qp = cast_int<unsigned int>(qp.size());
449 17026223 : this->_n_total_qp = n_qp;
450 :
451 : // Number of shape functions in the finite element approximation
452 : // space.
453 9140908 : const unsigned int n_approx_shape_functions =
454 7885315 : this->n_dofs(elem, this->get_order());
455 :
456 : // Maybe we already have correctly-sized data? Check data sizes,
457 : // and get ready to break out of a "loop" if all these resize()
458 : // calls are redundant.
459 4573487 : unsigned int old_n_qp = 0;
460 : do
461 : {
462 : // resize the vectors to hold current data
463 : // Phi are the shape functions used for the FE approximation
464 : // Phi_map are the shape functions used for the FE mapping
465 17026223 : if (this->calculate_phi)
466 : {
467 18715203 : if (this->phi.size() == n_approx_shape_functions)
468 : {
469 14153863 : old_n_qp = n_approx_shape_functions ? this->phi[0].size() : 0;
470 3801211 : break;
471 : }
472 597852 : this->phi.resize (n_approx_shape_functions);
473 : }
474 2872360 : if (this->calculate_dphi)
475 : {
476 1849314 : if (this->dphi.size() == n_approx_shape_functions)
477 : {
478 842592 : old_n_qp = n_approx_shape_functions ? this->dphi[0].size() : 0;
479 237739 : break;
480 : }
481 609339 : this->dphi.resize (n_approx_shape_functions);
482 609339 : this->dphidx.resize (n_approx_shape_functions);
483 609339 : this->dphidy.resize (n_approx_shape_functions);
484 609339 : this->dphidz.resize (n_approx_shape_functions);
485 : }
486 :
487 2028731 : if (this->calculate_dphiref)
488 : {
489 : if (Dim > 0)
490 : {
491 1255572 : if (this->dphidxi.size() == n_approx_shape_functions)
492 : {
493 12167 : old_n_qp = n_approx_shape_functions ? this->dphidxi[0].size() : 0;
494 2985 : break;
495 : }
496 986473 : this->dphidxi.resize (n_approx_shape_functions);
497 : }
498 :
499 : if (Dim > 1)
500 961860 : this->dphideta.resize (n_approx_shape_functions);
501 :
502 : if (Dim > 2)
503 786270 : this->dphidzeta.resize (n_approx_shape_functions);
504 : }
505 :
506 2017601 : if (this->calculate_curl_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
507 576 : this->curl_phi.resize(n_approx_shape_functions);
508 :
509 2017601 : if (this->calculate_div_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
510 1022 : this->div_phi.resize(n_approx_shape_functions);
511 :
512 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
513 2017601 : if (this->calculate_d2phi)
514 : {
515 518600 : if (this->d2phi.size() == n_approx_shape_functions)
516 : {
517 0 : old_n_qp = n_approx_shape_functions ? this->d2phi[0].size() : 0;
518 0 : break;
519 : }
520 :
521 414228 : this->d2phi.resize (n_approx_shape_functions);
522 414228 : this->d2phidx2.resize (n_approx_shape_functions);
523 414228 : this->d2phidxdy.resize (n_approx_shape_functions);
524 414228 : this->d2phidxdz.resize (n_approx_shape_functions);
525 414228 : this->d2phidy2.resize (n_approx_shape_functions);
526 414228 : this->d2phidydz.resize (n_approx_shape_functions);
527 414228 : this->d2phidz2.resize (n_approx_shape_functions);
528 :
529 : if (Dim > 0)
530 414228 : this->d2phidxi2.resize (n_approx_shape_functions);
531 :
532 : if (Dim > 1)
533 : {
534 391802 : this->d2phidxideta.resize (n_approx_shape_functions);
535 391802 : this->d2phideta2.resize (n_approx_shape_functions);
536 : }
537 : if (Dim > 2)
538 : {
539 384022 : this->d2phidxidzeta.resize (n_approx_shape_functions);
540 384022 : this->d2phidetadzeta.resize (n_approx_shape_functions);
541 384022 : this->d2phidzeta2.resize (n_approx_shape_functions);
542 : }
543 : }
544 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
545 : }
546 : while (false);
547 :
548 17026223 : if (old_n_qp != n_qp)
549 19473660 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
550 : {
551 17424149 : if (this->calculate_phi)
552 8135109 : this->phi[i].resize (n_qp);
553 :
554 17424149 : if (this->calculate_dphi)
555 : {
556 6621486 : this->dphi[i].resize (n_qp);
557 6621486 : this->dphidx[i].resize (n_qp);
558 6621486 : this->dphidy[i].resize (n_qp);
559 6621486 : this->dphidz[i].resize (n_qp);
560 : }
561 :
562 17423112 : if (this->calculate_dphiref)
563 : {
564 : if (Dim > 0)
565 10466806 : this->dphidxi[i].resize(n_qp);
566 :
567 : if (Dim > 1)
568 10339879 : this->dphideta[i].resize(n_qp);
569 :
570 : if (Dim > 2)
571 8709177 : this->dphidzeta[i].resize(n_qp);
572 : }
573 :
574 17424149 : if (this->calculate_curl_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
575 6252 : this->curl_phi[i].resize(n_qp);
576 :
577 17424149 : if (this->calculate_div_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
578 10966 : this->div_phi[i].resize(n_qp);
579 :
580 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
581 17424149 : if (this->calculate_d2phi)
582 : {
583 4605239 : this->d2phi[i].resize (n_qp);
584 4605239 : this->d2phidx2[i].resize (n_qp);
585 4605239 : this->d2phidxdy[i].resize (n_qp);
586 4605239 : this->d2phidxdz[i].resize (n_qp);
587 4605239 : this->d2phidy2[i].resize (n_qp);
588 4605239 : this->d2phidydz[i].resize (n_qp);
589 4605239 : this->d2phidz2[i].resize (n_qp);
590 : if (Dim > 0)
591 4605239 : this->d2phidxi2[i].resize (n_qp);
592 : if (Dim > 1)
593 : {
594 4490852 : this->d2phidxideta[i].resize (n_qp);
595 4490852 : this->d2phideta2[i].resize (n_qp);
596 : }
597 : if (Dim > 2)
598 : {
599 4390407 : this->d2phidxidzeta[i].resize (n_qp);
600 4390407 : this->d2phidetadzeta[i].resize (n_qp);
601 4390407 : this->d2phidzeta2[i].resize (n_qp);
602 : }
603 : }
604 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
605 : }
606 :
607 :
608 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
609 : //------------------------------------------------------------
610 : // Initialize the data fields, which should only be used for infinite
611 : // elements, to some sensible values, so that using a FE with the
612 : // variational formulation of an InfFE, correct element matrices are
613 : // returned
614 :
615 : {
616 12737002 : if (this->calculate_phi || this->calculate_dphi)
617 : {
618 11976510 : this->weight.resize (n_qp);
619 48246902 : for (unsigned int p=0; p<n_qp; p++)
620 49230953 : this->weight[p] = 1.;
621 : }
622 :
623 12737002 : if (this->calculate_dphi)
624 : {
625 4737088 : this->dweight.resize (n_qp);
626 4737088 : this->dphase.resize (n_qp);
627 15965547 : for (unsigned int p=0; p<n_qp; p++)
628 : {
629 11228459 : this->dweight[p].zero();
630 8267578 : this->dphase[p].zero();
631 : }
632 : }
633 : }
634 : #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
635 :
636 : // Compute the values of the shape function derivatives
637 16957803 : if (this->calculate_dphiref && Dim > 0)
638 : {
639 13966371 : std::vector<std::vector<OutputShape>> * comps[3]
640 13966371 : { &this->dphidxi, &this->dphideta, &this->dphidzeta };
641 6708507 : FE<Dim,T>::all_shape_derivs(elem, this->fe_type.order, qp, comps, this->_add_p_level_in_reinit);
642 : }
643 :
644 : switch (Dim)
645 : {
646 :
647 : //------------------------------------------------------------
648 : // 0D
649 : case 0:
650 : {
651 19049 : break;
652 : }
653 :
654 : //------------------------------------------------------------
655 : // 1D
656 476031 : case 1:
657 : {
658 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
659 : // Compute the value of shape function i Hessians at quadrature point p
660 704763 : if (this->calculate_d2phi)
661 152317 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
662 522455 : for (unsigned int p=0; p<n_qp; p++)
663 401428 : this->d2phidxi2[i][p] = FE<Dim, T>::shape_second_deriv(
664 401428 : elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
665 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
666 :
667 228732 : break;
668 : }
669 :
670 :
671 :
672 : //------------------------------------------------------------
673 : // 2D
674 6034152 : case 2:
675 : {
676 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
677 : // Compute the value of shape function i Hessians at quadrature point p
678 8357237 : if (this->calculate_d2phi)
679 3470612 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
680 17574603 : for (unsigned int p=0; p<n_qp; p++)
681 : {
682 14414929 : this->d2phidxi2[i][p] = FE<Dim, T>::shape_second_deriv(
683 14414929 : elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
684 21498204 : this->d2phidxideta[i][p] = FE<Dim, T>::shape_second_deriv(
685 14414929 : elem, this->fe_type.order, i, 1, qp[p], this->_add_p_level_in_reinit);
686 14414929 : this->d2phideta2[i][p] = FE<Dim, T>::shape_second_deriv(
687 14414929 : elem, this->fe_type.order, i, 2, qp[p], this->_add_p_level_in_reinit);
688 : }
689 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
690 :
691 :
692 2323085 : break;
693 : }
694 :
695 :
696 :
697 : //------------------------------------------------------------
698 : // 3D
699 5893182 : case 3:
700 : {
701 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
702 : // Compute the value of shape function i Hessians at quadrature point p
703 7895803 : if (this->calculate_d2phi)
704 53660616 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
705 154375354 : for (unsigned int p=0; p<n_qp; p++)
706 : {
707 102570323 : this->d2phidxi2[i][p] = FE<Dim, T>::shape_second_deriv(
708 102570323 : elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
709 152686110 : this->d2phidxideta[i][p] = FE<Dim, T>::shape_second_deriv(
710 102570323 : elem, this->fe_type.order, i, 1, qp[p], this->_add_p_level_in_reinit);
711 152686110 : this->d2phideta2[i][p] = FE<Dim, T>::shape_second_deriv(
712 102570323 : elem, this->fe_type.order, i, 2, qp[p], this->_add_p_level_in_reinit);
713 152686110 : this->d2phidxidzeta[i][p] = FE<Dim, T>::shape_second_deriv(
714 102570323 : elem, this->fe_type.order, i, 3, qp[p], this->_add_p_level_in_reinit);
715 152686110 : this->d2phidetadzeta[i][p] = FE<Dim, T>::shape_second_deriv(
716 102570323 : elem, this->fe_type.order, i, 4, qp[p], this->_add_p_level_in_reinit);
717 102570323 : this->d2phidzeta2[i][p] = FE<Dim, T>::shape_second_deriv(
718 102570323 : elem, this->fe_type.order, i, 5, qp[p], this->_add_p_level_in_reinit);
719 : }
720 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
721 :
722 2002621 : break;
723 : }
724 :
725 :
726 : default:
727 : libmesh_error_msg("Invalid dimension Dim = " << Dim);
728 : }
729 :
730 17026223 : if (this->calculate_dual)
731 1510 : this->init_dual_shape_functions(n_approx_shape_functions, n_qp);
732 17026223 : }
733 :
734 : template <unsigned int Dim, FEFamily T>
735 : void
736 16771853 : FE<Dim,T>::default_all_shape_derivs (const Elem * elem,
737 : const Order o,
738 : const std::vector<Point> & p,
739 : std::vector<std::vector<OutputShape>> * comps[3],
740 : const bool add_p_level)
741 : {
742 57910761 : for (unsigned int d=0; d != Dim; ++d)
743 : {
744 41138908 : auto & comps_d = *comps[d];
745 550545891 : for (auto i : index_range(comps_d))
746 : FE<Dim,T>::shape_derivs
747 642402647 : (elem,o,i,d,p,comps_d[i],add_p_level);
748 : }
749 16771853 : }
750 :
751 :
752 : template <unsigned int Dim, FEFamily T>
753 : void
754 1072 : FE<Dim,T>::default_side_nodal_soln(const Elem * elem, const Order o,
755 : const unsigned int side,
756 : const std::vector<Number> & elem_soln,
757 : std::vector<Number> & nodal_soln_on_side,
758 : const bool add_p_level,
759 : const unsigned vdim)
760 : {
761 536 : std::vector<Number> full_nodal_soln;
762 1072 : nodal_soln(elem, o, elem_soln, full_nodal_soln, add_p_level, vdim);
763 1072 : const std::vector<unsigned int> side_nodes =
764 1072 : elem->nodes_on_side(side);
765 :
766 536 : std::size_t n_side_nodes = side_nodes.size();
767 1072 : nodal_soln_on_side.resize(n_side_nodes);
768 7648 : for (auto n : make_range(n_side_nodes))
769 11508 : nodal_soln_on_side[n] = full_nodal_soln[side_nodes[n]];
770 1072 : }
771 :
772 :
773 :
774 :
775 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
776 :
777 : template <unsigned int Dim, FEFamily T>
778 2662 : void FE<Dim,T>::init_base_shape_functions(const std::vector<Point> & qp,
779 : const Elem * e)
780 : {
781 2662 : this->_elem = e;
782 2662 : this->_elem_type = e->type();
783 2662 : this->_fe_map->template init_reference_to_physical_map<Dim>(qp, e);
784 2662 : init_shape_functions(qp, e);
785 2662 : }
786 :
787 : #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS
788 :
789 :
790 : // Helper for FDM methods
791 : namespace {
792 : using namespace libMesh;
793 :
794 : std::tuple<Point, Point, Real>
795 1176980239 : fdm_points(const unsigned int j, const Point & p)
796 : {
797 301067220 : libmesh_assert_less (j, LIBMESH_DIM);
798 :
799 : // cheat by using finite difference approximations:
800 301067220 : const Real eps = 1.e-6;
801 1176980239 : Point pp = p, pm = p;
802 :
803 1176980239 : switch (j)
804 : {
805 : // d()/dxi
806 139687818 : case 0:
807 : {
808 546396762 : pp(0) += eps;
809 546396762 : pm(0) -= eps;
810 546396762 : break;
811 : }
812 :
813 : // d()/deta
814 102058300 : case 1:
815 : {
816 398366316 : pp(1) += eps;
817 398366316 : pm(1) -= eps;
818 398366316 : break;
819 : }
820 :
821 : // d()/dzeta
822 59321102 : case 2:
823 : {
824 232217161 : pp(2) += eps;
825 232217161 : pm(2) -= eps;
826 232217161 : break;
827 : }
828 :
829 0 : default:
830 0 : libmesh_error_msg("Invalid derivative index j = " << j);
831 : }
832 :
833 1478047459 : return std::make_tuple(pp, pm, eps);
834 : }
835 :
836 : std::tuple<Point, Point, Real, unsigned int>
837 540342048 : fdm_second_points(const unsigned int j, const Point & p)
838 : {
839 : // cheat by using finite difference approximations:
840 136926207 : const Real eps = 1.e-5;
841 540342048 : Point pp = p, pm = p;
842 136926207 : unsigned int deriv_j = 0;
843 :
844 540342048 : switch (j)
845 : {
846 : // d^2() / dxi^2
847 23406539 : case 0:
848 : {
849 92373827 : pp(0) += eps;
850 92373827 : pm(0) -= eps;
851 23406539 : deriv_j = 0;
852 92373827 : break;
853 : }
854 :
855 : // d^2() / dxi deta
856 23406539 : case 1:
857 : {
858 92373827 : pp(1) += eps;
859 92373827 : pm(1) -= eps;
860 23406539 : deriv_j = 0;
861 92373827 : break;
862 : }
863 :
864 : // d^2() / deta^2
865 23406539 : case 2:
866 : {
867 92373827 : pp(1) += eps;
868 92373827 : pm(1) -= eps;
869 23406539 : deriv_j = 1;
870 92373827 : break;
871 : }
872 :
873 : // d^2()/dxidzeta
874 22235530 : case 3:
875 : {
876 87740189 : pp(2) += eps;
877 87740189 : pm(2) -= eps;
878 22235530 : deriv_j = 0;
879 87740189 : break;
880 : } // d^2()/deta^2
881 :
882 : // d^2()/detadzeta
883 22235530 : case 4:
884 : {
885 87740189 : pp(2) += eps;
886 87740189 : pm(2) -= eps;
887 22235530 : deriv_j = 1;
888 87740189 : break;
889 : }
890 :
891 : // d^2()/dzeta^2
892 22235530 : case 5:
893 : {
894 87740189 : pp(2) += eps;
895 87740189 : pm(2) -= eps;
896 22235530 : deriv_j = 2;
897 87740189 : break;
898 : }
899 :
900 0 : default:
901 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
902 : }
903 :
904 677268255 : return std::make_tuple(pp, pm, eps, deriv_j);
905 : }
906 :
907 : }
908 :
909 :
910 : template <typename OutputShape>
911 1137542905 : OutputShape fe_fdm_deriv(const Elem * elem,
912 : const Order order,
913 : const unsigned int i,
914 : const unsigned int j,
915 : const Point & p,
916 : const bool add_p_level,
917 : OutputShape(*shape_func)
918 : (const Elem *, const Order,
919 : const unsigned int, const Point &,
920 : const bool))
921 : {
922 290037936 : libmesh_assert(elem);
923 :
924 1137542905 : auto [pp, pm, eps] = fdm_points(j, p);
925 :
926 1137542905 : return (shape_func(elem, order, i, pp, add_p_level) -
927 1137542905 : shape_func(elem, order, i, pm, add_p_level))/2./eps;
928 : }
929 :
930 :
931 : template <typename OutputShape>
932 0 : OutputShape fe_fdm_deriv(const ElemType type,
933 : const Order order,
934 : const unsigned int i,
935 : const unsigned int j,
936 : const Point & p,
937 : OutputShape(*shape_func)
938 : (const ElemType, const Order,
939 : const unsigned int, const Point &))
940 : {
941 0 : auto [pp, pm, eps] = fdm_points(j, p);
942 :
943 0 : return (shape_func(type, order, i, pp) -
944 0 : shape_func(type, order, i, pm))/2./eps;
945 : }
946 :
947 :
948 : template <typename OutputShape>
949 39437334 : OutputShape fe_fdm_deriv(const ElemType type,
950 : const Order order,
951 : const Elem * elem,
952 : const unsigned int i,
953 : const unsigned int j,
954 : const Point & p,
955 : OutputShape(*shape_func)
956 : (const ElemType, const Order,
957 : const Elem *,
958 : const unsigned int, const Point &))
959 : {
960 39437334 : auto [pp, pm, eps] = fdm_points(j, p);
961 :
962 39437334 : return (shape_func(type, order, elem, i, pp) -
963 39437334 : shape_func(type, order, elem, i, pm))/2./eps;
964 : }
965 :
966 :
967 : template <typename OutputShape>
968 : OutputShape
969 538590504 : fe_fdm_second_deriv(const Elem * elem,
970 : const Order order,
971 : const unsigned int i,
972 : const unsigned int j,
973 : const Point & p,
974 : const bool add_p_level,
975 : OutputShape(*deriv_func)
976 : (const Elem *, const Order,
977 : const unsigned int, const unsigned int,
978 : const Point &, const bool))
979 : {
980 538590504 : auto [pp, pm, eps, deriv_j] = fdm_second_points(j, p);
981 :
982 538590504 : return (deriv_func(elem, order, i, deriv_j, pp, add_p_level) -
983 538590504 : deriv_func(elem, order, i, deriv_j, pm, add_p_level))/2./eps;
984 : }
985 :
986 :
987 : template <typename OutputShape>
988 : OutputShape
989 0 : fe_fdm_second_deriv(const ElemType type,
990 : const Order order,
991 : const unsigned int i,
992 : const unsigned int j,
993 : const Point & p,
994 : OutputShape(*deriv_func)
995 : (const ElemType, const Order,
996 : const unsigned int, const unsigned int,
997 : const Point &))
998 : {
999 0 : auto [pp, pm, eps, deriv_j] = fdm_second_points(j, p);
1000 :
1001 0 : return (deriv_func(type, order, i, deriv_j, pp) -
1002 0 : deriv_func(type, order, i, deriv_j, pm))/2./eps;
1003 : }
1004 :
1005 :
1006 : template <typename OutputShape>
1007 : OutputShape
1008 1751544 : fe_fdm_second_deriv(const ElemType type,
1009 : const Order order,
1010 : const Elem * elem,
1011 : const unsigned int i,
1012 : const unsigned int j,
1013 : const Point & p,
1014 : OutputShape(*deriv_func)
1015 : (const ElemType, const Order,
1016 : const Elem *,
1017 : const unsigned int, const unsigned int,
1018 : const Point &))
1019 : {
1020 1751544 : auto [pp, pm, eps, deriv_j] = fdm_second_points(j, p);
1021 :
1022 1751544 : return (deriv_func(type, order, elem, i, deriv_j, pp) -
1023 1751544 : deriv_func(type, order, elem, i, deriv_j, pm))/2./eps;
1024 : }
1025 :
1026 :
1027 151405 : void rational_fe_weighted_shapes(const Elem * elem,
1028 : const FEType underlying_fe_type,
1029 : std::vector<std::vector<Real>> & shapes,
1030 : const std::vector<Point> & p,
1031 : const bool add_p_level)
1032 : {
1033 151405 : const int extra_order = add_p_level * elem->p_level();
1034 :
1035 151405 : const int dim = elem->dim();
1036 :
1037 : const unsigned int n_sf =
1038 151405 : FEInterface::n_shape_functions(underlying_fe_type, extra_order,
1039 40181 : elem);
1040 :
1041 40181 : libmesh_assert_equal_to (n_sf, elem->n_nodes());
1042 :
1043 191586 : std::vector<Real> node_weights(n_sf);
1044 :
1045 80362 : const unsigned char datum_index = elem->mapping_data();
1046 1955740 : for (unsigned int n=0; n<n_sf; n++)
1047 2273228 : node_weights[n] =
1048 1804335 : elem->node_ref(n).get_extra_datum<Real>(datum_index);
1049 :
1050 80362 : const std::size_t n_p = p.size();
1051 :
1052 151405 : shapes.resize(n_sf);
1053 1955740 : for (unsigned int i=0; i != n_sf; ++i)
1054 : {
1055 1804335 : auto & shapes_i = shapes[i];
1056 1804335 : shapes_i.resize(n_p, 0);
1057 1804335 : FEInterface::shapes(dim, underlying_fe_type, elem, i, p,
1058 : shapes_i, add_p_level);
1059 9602846 : for (auto & s : shapes_i)
1060 9863704 : s *= node_weights[i];
1061 : }
1062 151405 : }
1063 :
1064 :
1065 194273 : void rational_fe_weighted_shapes_derivs(const Elem * elem,
1066 : const FEType fe_type,
1067 : std::vector<std::vector<Real>> & shapes,
1068 : std::vector<std::vector<std::vector<Real>>> & derivs,
1069 : const std::vector<Point> & p,
1070 : const bool add_p_level)
1071 : {
1072 194273 : const int extra_order = add_p_level * elem->p_level();
1073 194273 : const unsigned int dim = elem->dim();
1074 :
1075 : const unsigned int n_sf =
1076 194273 : FEInterface::n_shape_functions(fe_type, extra_order, elem);
1077 :
1078 52123 : libmesh_assert_equal_to (n_sf, elem->n_nodes());
1079 :
1080 52123 : libmesh_assert_equal_to (dim, derivs.size());
1081 635796 : for (unsigned int d = 0; d != dim; ++d)
1082 559083 : derivs[d].resize(n_sf);
1083 :
1084 246396 : std::vector<Real> node_weights(n_sf);
1085 :
1086 104246 : const unsigned char datum_index = elem->mapping_data();
1087 2372773 : for (unsigned int n=0; n<n_sf; n++)
1088 2752816 : node_weights[n] =
1089 2178500 : elem->node_ref(n).get_extra_datum<Real>(datum_index);
1090 :
1091 104246 : const std::size_t n_p = p.size();
1092 :
1093 194273 : shapes.resize(n_sf);
1094 2372773 : for (unsigned int i=0; i != n_sf; ++i)
1095 2752816 : shapes[i].resize(n_p, 0);
1096 :
1097 194273 : FEInterface::all_shapes(dim, fe_type, elem, p, shapes, add_p_level);
1098 :
1099 2372773 : for (unsigned int i=0; i != n_sf; ++i)
1100 : {
1101 2178500 : auto & shapes_i = shapes[i];
1102 :
1103 15187359 : for (auto & s : shapes_i)
1104 16492279 : s *= node_weights[i];
1105 :
1106 7480432 : for (unsigned int d = 0; d != dim; ++d)
1107 : {
1108 6685606 : auto & derivs_di = derivs[d][i];
1109 5301932 : derivs_di.resize(n_p);
1110 5301932 : FEInterface::shape_derivs(fe_type, elem, i, d, p,
1111 : derivs_di, add_p_level);
1112 34445378 : for (auto & dip : derivs_di)
1113 36866518 : dip *= node_weights[i];
1114 : }
1115 : }
1116 194273 : }
1117 :
1118 :
1119 7852905 : Real rational_fe_shape(const Elem & elem,
1120 : const FEType underlying_fe_type,
1121 : const unsigned int i,
1122 : const Point & p,
1123 : const bool add_p_level)
1124 : {
1125 7852905 : int extra_order = add_p_level * elem.p_level();
1126 :
1127 : const unsigned int n_sf =
1128 7852905 : FEInterface::n_shape_functions(underlying_fe_type, extra_order, &elem);
1129 :
1130 2237526 : libmesh_assert_equal_to (n_sf, elem.n_nodes());
1131 :
1132 7852905 : std::vector<Real> node_weights(n_sf);
1133 :
1134 4348602 : const unsigned char datum_index = elem.mapping_data();
1135 :
1136 2237526 : Real weighted_shape_i = 0, weighted_sum = 0;
1137 :
1138 103020924 : for (unsigned int sf=0; sf<n_sf; sf++)
1139 : {
1140 : Real node_weight =
1141 95168019 : elem.node_ref(sf).get_extra_datum<Real>(datum_index);
1142 : Real weighted_shape = node_weight *
1143 95168019 : FEInterface::shape(underlying_fe_type, extra_order, &elem, sf, p);
1144 95168019 : weighted_sum += weighted_shape;
1145 95168019 : if (sf == i)
1146 2237526 : weighted_shape_i = weighted_shape;
1147 : }
1148 :
1149 12201507 : return weighted_shape_i / weighted_sum;
1150 : }
1151 :
1152 :
1153 15580257 : Real rational_fe_shape_deriv(const Elem & elem,
1154 : const FEType underlying_fe_type,
1155 : const unsigned int i,
1156 : const unsigned int j,
1157 : const Point & p,
1158 : const bool add_p_level)
1159 : {
1160 4186500 : libmesh_assert_less(j, elem.dim());
1161 :
1162 15580257 : int extra_order = add_p_level * elem.p_level();
1163 :
1164 : const unsigned int n_sf =
1165 15580257 : FEInterface::n_shape_functions(underlying_fe_type, extra_order, &elem);
1166 :
1167 15580257 : const unsigned int n_nodes = elem.n_nodes();
1168 4186500 : libmesh_assert_equal_to (n_sf, n_nodes);
1169 :
1170 15580257 : std::vector<Real> node_weights(n_nodes);
1171 :
1172 8373000 : const unsigned char datum_index = elem.mapping_data();
1173 232431996 : for (unsigned int n=0; n<n_nodes; n++)
1174 273411939 : node_weights[n] =
1175 216851739 : elem.node_ref(n).get_extra_datum<Real>(datum_index);
1176 :
1177 4186500 : Real weighted_shape_i = 0, weighted_sum = 0,
1178 4186500 : weighted_grad_i = 0, weighted_grad_sum = 0;
1179 :
1180 232431996 : for (unsigned int sf=0; sf<n_sf; sf++)
1181 : {
1182 216851739 : Real weighted_shape = node_weights[sf] *
1183 216851739 : FEInterface::shape(underlying_fe_type, extra_order, &elem, sf, p);
1184 216851739 : Real weighted_grad = node_weights[sf] *
1185 216851739 : FEInterface::shape_deriv(underlying_fe_type, extra_order, &elem, sf, j, p);
1186 216851739 : weighted_sum += weighted_shape;
1187 216851739 : weighted_grad_sum += weighted_grad;
1188 216851739 : if (sf == i)
1189 : {
1190 4186500 : weighted_shape_i = weighted_shape;
1191 4186500 : weighted_grad_i = weighted_grad;
1192 : }
1193 : }
1194 :
1195 15580257 : return (weighted_sum * weighted_grad_i - weighted_shape_i * weighted_grad_sum) /
1196 23953257 : weighted_sum / weighted_sum;
1197 : }
1198 :
1199 :
1200 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1201 :
1202 4783445 : Real rational_fe_shape_second_deriv(const Elem & elem,
1203 : const FEType underlying_fe_type,
1204 : const unsigned int i,
1205 : const unsigned int j,
1206 : const Point & p,
1207 : const bool add_p_level)
1208 : {
1209 : unsigned int j1, j2;
1210 1199654 : switch (j)
1211 : {
1212 204997 : case 0:
1213 : // j = 0 ==> d^2 phi / dxi^2
1214 204997 : j1 = j2 = 0;
1215 204997 : break;
1216 204668 : case 1:
1217 : // j = 1 ==> d^2 phi / dxi deta
1218 204668 : j1 = 0;
1219 204668 : j2 = 1;
1220 204668 : break;
1221 204668 : case 2:
1222 : // j = 2 ==> d^2 phi / deta^2
1223 204668 : j1 = j2 = 1;
1224 204668 : break;
1225 195107 : case 3:
1226 : // j = 3 ==> d^2 phi / dxi dzeta
1227 195107 : j1 = 0;
1228 195107 : j2 = 2;
1229 195107 : break;
1230 195107 : case 4:
1231 : // j = 4 ==> d^2 phi / deta dzeta
1232 195107 : j1 = 1;
1233 195107 : j2 = 2;
1234 195107 : break;
1235 195107 : case 5:
1236 : // j = 5 ==> d^2 phi / dzeta^2
1237 195107 : j1 = j2 = 2;
1238 195107 : break;
1239 0 : default:
1240 0 : libmesh_error();
1241 : }
1242 :
1243 4783445 : int extra_order = add_p_level * elem.p_level();
1244 :
1245 : const unsigned int n_sf =
1246 4783445 : FEInterface::n_shape_functions(underlying_fe_type, extra_order,
1247 1199654 : &elem);
1248 :
1249 4783445 : const unsigned int n_nodes = elem.n_nodes();
1250 1199654 : libmesh_assert_equal_to (n_sf, n_nodes);
1251 :
1252 4783445 : std::vector<Real> node_weights(n_nodes);
1253 :
1254 2399308 : const unsigned char datum_index = elem.mapping_data();
1255 109817892 : for (unsigned int n=0; n<n_nodes; n++)
1256 131307503 : node_weights[n] =
1257 105034447 : elem.node_ref(n).get_extra_datum<Real>(datum_index);
1258 :
1259 1199654 : Real weighted_shape_i = 0, weighted_sum = 0,
1260 1199654 : weighted_grada_i = 0, weighted_grada_sum = 0,
1261 1199654 : weighted_gradb_i = 0, weighted_gradb_sum = 0,
1262 1199654 : weighted_hess_i = 0, weighted_hess_sum = 0;
1263 :
1264 109817892 : for (unsigned int sf=0; sf<n_sf; sf++)
1265 : {
1266 105034447 : Real weighted_shape = node_weights[sf] *
1267 105034447 : FEInterface::shape(underlying_fe_type, extra_order, &elem, sf,
1268 105034447 : p);
1269 105034447 : Real weighted_grada = node_weights[sf] *
1270 105034447 : FEInterface::shape_deriv(underlying_fe_type, extra_order,
1271 105034447 : &elem, sf, j1, p);
1272 105034447 : Real weighted_hess = node_weights[sf] *
1273 105034447 : FEInterface::shape_second_deriv(underlying_fe_type,
1274 105034447 : extra_order, &elem, sf, j, p);
1275 105034447 : weighted_sum += weighted_shape;
1276 105034447 : weighted_grada_sum += weighted_grada;
1277 26273056 : Real weighted_gradb = weighted_grada;
1278 105034447 : if (j1 != j2)
1279 : {
1280 52391157 : weighted_gradb = (j1 == j2) ? weighted_grada :
1281 52391157 : node_weights[sf] *
1282 52391157 : FEInterface::shape_deriv(underlying_fe_type, extra_order,
1283 : &elem, sf, j2, p);
1284 52391157 : weighted_grada_sum += weighted_grada;
1285 : }
1286 105034447 : weighted_hess_sum += weighted_hess;
1287 105034447 : if (sf == i)
1288 : {
1289 1199654 : weighted_shape_i = weighted_shape;
1290 1199654 : weighted_grada_i = weighted_grada;
1291 1199654 : weighted_gradb_i = weighted_gradb;
1292 1199654 : weighted_hess_i = weighted_hess;
1293 : }
1294 : }
1295 :
1296 4783445 : if (j1 == j2)
1297 604772 : weighted_gradb_sum = weighted_grada_sum;
1298 :
1299 7182753 : return (weighted_sum * weighted_hess_i - weighted_grada_i * weighted_gradb_sum -
1300 9582061 : weighted_shape_i * weighted_hess_sum - weighted_gradb_i * weighted_grada_sum +
1301 7182753 : 2 * weighted_grada_sum * weighted_shape_i * weighted_gradb_sum / weighted_sum) /
1302 7182753 : weighted_sum / weighted_sum;
1303 : }
1304 :
1305 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1306 :
1307 :
1308 151405 : void rational_all_shapes (const Elem & elem,
1309 : const FEType underlying_fe_type,
1310 : const std::vector<Point> & p,
1311 : std::vector<std::vector<Real>> & v,
1312 : const bool add_p_level)
1313 : {
1314 120543 : std::vector<std::vector<Real>> shapes;
1315 :
1316 151405 : rational_fe_weighted_shapes(&elem, underlying_fe_type, shapes, p,
1317 : add_p_level);
1318 :
1319 231767 : std::vector<Real> shape_sums(p.size(), 0);
1320 :
1321 1955740 : for (auto i : index_range(v))
1322 : {
1323 468893 : libmesh_assert_equal_to ( p.size(), shapes[i].size() );
1324 9602846 : for (auto j : index_range(p))
1325 13994090 : shape_sums[j] += shapes[i][j];
1326 : }
1327 :
1328 1955740 : for (auto i : index_range(v))
1329 : {
1330 468893 : libmesh_assert_equal_to ( p.size(), v[i].size() );
1331 9602846 : for (auto j : index_range(v[i]))
1332 16059283 : v[i][j] = shapes[i][j] / shape_sums[j];
1333 : }
1334 151405 : }
1335 :
1336 :
1337 : template <typename OutputShape>
1338 194273 : void rational_all_shape_derivs (const Elem & elem,
1339 : const FEType underlying_fe_type,
1340 : const std::vector<Point> & p,
1341 : std::vector<std::vector<OutputShape>> * comps[3],
1342 : const bool add_p_level)
1343 : {
1344 194273 : const int my_dim = elem.dim();
1345 :
1346 156369 : std::vector<std::vector<Real>> shapes;
1347 298519 : std::vector<std::vector<std::vector<Real>>> derivs(my_dim);
1348 :
1349 194273 : rational_fe_weighted_shapes_derivs(&elem, underlying_fe_type,
1350 : shapes, derivs, p, add_p_level);
1351 :
1352 298519 : std::vector<Real> shape_sums(p.size(), 0);
1353 298519 : std::vector<std::vector<Real>> shape_deriv_sums(my_dim);
1354 635796 : for (int d=0; d != my_dim; ++d)
1355 676643 : shape_deriv_sums[d].resize(p.size());
1356 :
1357 2372773 : for (auto i : index_range(shapes))
1358 : {
1359 574316 : libmesh_assert_equal_to ( p.size(), shapes[i].size() );
1360 15187359 : for (auto j : index_range(p))
1361 19975699 : shape_sums[j] += shapes[i][j];
1362 :
1363 7480432 : for (int d=0; d != my_dim; ++d)
1364 34445378 : for (auto j : index_range(p))
1365 67758806 : shape_deriv_sums[d][j] += derivs[d][i][j];
1366 : }
1367 :
1368 635796 : for (int d=0; d != my_dim; ++d)
1369 : {
1370 441523 : auto & comps_d = *comps[d];
1371 117560 : libmesh_assert_equal_to(comps_d.size(), elem.n_nodes());
1372 :
1373 5743455 : for (auto i : index_range(comps_d))
1374 : {
1375 1383674 : auto & comps_di = comps_d[i];
1376 4151022 : auto & derivs_di = derivs[d][i];
1377 :
1378 34445378 : for (auto j : index_range(comps_di))
1379 60035734 : comps_di[j] = (shape_sums[j] * derivs_di[j] -
1380 44589590 : shapes[i][j] * shape_deriv_sums[d][j]) /
1381 29143446 : shape_sums[j] / shape_sums[j];
1382 : }
1383 : }
1384 284300 : }
1385 :
1386 :
1387 :
1388 : template
1389 : Real fe_fdm_deriv<Real>(const Elem *, const Order, const unsigned int,
1390 : const unsigned int, const Point &, const bool,
1391 : Real(*shape_func)
1392 : (const Elem *, const Order, const unsigned int,
1393 : const Point &, const bool));
1394 :
1395 : template
1396 : Real fe_fdm_deriv<Real>(const ElemType, const Order, const unsigned int,
1397 : const unsigned int, const Point &,
1398 : Real(*shape_func)
1399 : (const ElemType, const Order, const unsigned int,
1400 : const Point &));
1401 :
1402 : template
1403 : Real fe_fdm_deriv<Real>(const ElemType, const Order, const Elem *,
1404 : const unsigned int, const unsigned int, const Point &,
1405 : Real(*shape_func)
1406 : (const ElemType, const Order, const Elem *,
1407 : const unsigned int, const Point &));
1408 :
1409 : template
1410 : RealGradient
1411 : fe_fdm_deriv<RealGradient>(const Elem *, const Order, const unsigned int,
1412 : const unsigned int, const Point &, const bool,
1413 : RealGradient(*shape_func)
1414 : (const Elem *, const Order, const unsigned int,
1415 : const Point &, const bool));
1416 :
1417 : template
1418 : Real
1419 : fe_fdm_second_deriv<Real>(const ElemType, const Order, const unsigned int,
1420 : const unsigned int, const Point &,
1421 : Real(*shape_func)
1422 : (const ElemType, const Order, const unsigned int,
1423 : const unsigned int, const Point &));
1424 :
1425 : template
1426 : Real
1427 : fe_fdm_second_deriv<Real>(const Elem *, const Order, const unsigned int,
1428 : const unsigned int, const Point &, const bool,
1429 : Real(*shape_func)
1430 : (const Elem *, const Order, const unsigned int,
1431 : const unsigned int, const Point &, const bool));
1432 :
1433 : template
1434 : Real
1435 : fe_fdm_second_deriv<Real>(const ElemType, const Order, const Elem *,
1436 : const unsigned int, const unsigned int, const Point &,
1437 : Real(*shape_func)
1438 : (const ElemType, const Order, const Elem *,
1439 : const unsigned int, const unsigned int, const Point &));
1440 :
1441 : template
1442 : RealGradient
1443 : fe_fdm_second_deriv<RealGradient>(const Elem *, const Order, const unsigned int,
1444 : const unsigned int, const Point &, const bool,
1445 : RealGradient(*shape_func)
1446 : (const Elem *, const Order, const unsigned int,
1447 : const unsigned int, const Point &, const bool));
1448 :
1449 :
1450 :
1451 :
1452 : //--------------------------------------------------------------
1453 : // Explicit instantiations using macro from fe_macro.h
1454 :
1455 : INSTANTIATE_FE(0);
1456 :
1457 : INSTANTIATE_FE(1);
1458 :
1459 : INSTANTIATE_FE(2);
1460 :
1461 : INSTANTIATE_FE(3);
1462 :
1463 : INSTANTIATE_SUBDIVISION_FE;
1464 :
1465 : template LIBMESH_EXPORT void rational_all_shape_derivs<Real> (const Elem & elem, const FEType underlying_fe_type, const std::vector<Point> & p, std::vector<std::vector<Real>> * comps[3], const bool add_p_level);
1466 : } // namespace libMesh
|
