Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/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 372 : void nonlagrange_dual_warning () {
37 : libmesh_warning("dual calculations have only been verified for the LAGRANGE family");
38 372 : }
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 16843 : void setup() { caching = new bool(!on_command_line("--disable-caching")); }
54 : public:
55 16389 : ~CachingSetup() { delete caching; caching = nullptr; }
56 : } caching_setup;
57 :
58 :
59 : // ------------------------------------------------------------
60 : // FE class members
61 : template <unsigned int Dim, FEFamily T>
62 15325910 : FE<Dim,T>::FE (const FEType & fet) :
63 : FEGenericBase<typename FEOutputType<T>::type> (Dim,fet),
64 13587597 : last_side(INVALID_ELEM),
65 15325910 : 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 869451 : libmesh_assert_equal_to (T, this->get_family());
71 15325910 : }
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 12884294 : void FE<Dim,T>::attach_quadrature_rule (QBase * q)
88 : {
89 589951 : libmesh_assert(q);
90 12884294 : this->qrule = q;
91 : // make sure we don't cache results from a previous quadrature rule
92 12884294 : this->_elem = nullptr;
93 12884294 : this->_elem_type = INVALID_ELEM;
94 12884294 : return;
95 : }
96 :
97 :
98 : template <unsigned int Dim, FEFamily T>
99 476885 : 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 40563 : libmesh_assert(elem);
106 40563 : libmesh_assert_less (s, elem->n_sides());
107 :
108 40563 : di.clear();
109 40563 : unsigned int nodenum = 0;
110 476885 : const unsigned int n_nodes = elem->n_nodes();
111 6848420 : for (unsigned int n = 0; n != n_nodes; ++n)
112 : {
113 : const unsigned int n_dofs =
114 7392037 : n_dofs_at_node(*elem, static_cast<Order>(o + add_p_level*elem->p_level()), n);
115 6371535 : if (elem->is_node_on_side(n, s))
116 6171404 : for (unsigned int i = 0; i != n_dofs; ++i)
117 3456875 : di.push_back(nodenum++);
118 : else
119 3657006 : nodenum += n_dofs;
120 : }
121 476885 : }
122 :
123 :
124 :
125 : template <unsigned int Dim, FEFamily T>
126 119562 : 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 8709 : libmesh_assert(elem);
133 8709 : libmesh_assert_less (e, elem->n_edges());
134 :
135 8709 : di.clear();
136 8709 : unsigned int nodenum = 0;
137 119562 : const unsigned int n_nodes = elem->n_nodes();
138 2505660 : for (unsigned int n = 0; n != n_nodes; ++n)
139 : {
140 : const unsigned int n_dofs =
141 2737362 : n_dofs_at_node(*elem, static_cast<Order>(o + add_p_level*elem->p_level()), n);
142 2386098 : if (elem->is_node_on_edge(n, e))
143 818043 : for (unsigned int i = 0; i != n_dofs; ++i)
144 466305 : di.push_back(nodenum++);
145 : else
146 2034360 : nodenum += n_dofs;
147 : }
148 119562 : }
149 :
150 :
151 :
152 : template <unsigned int Dim, FEFamily T>
153 4774234 : void FE<Dim,T>::cache(const Elem * elem)
154 : {
155 4774234 : cached_nodes.resize(elem->n_nodes());
156 57296552 : for (auto n : elem->node_index_range())
157 61212954 : cached_nodes[n] = elem->point(n);
158 :
159 4774234 : if (FEInterface::orientation_dependent(T))
160 : {
161 3655740 : cached_edges.resize(elem->n_edges());
162 22631232 : for (auto n : elem->edge_index_range())
163 18975492 : cached_edges[n] = elem->positive_edge_orientation(n);
164 :
165 3655740 : cached_faces.resize(elem->n_faces());
166 12026854 : for (auto n : elem->face_index_range())
167 8371114 : cached_faces[n] = elem->positive_face_orientation(n);
168 : }
169 4774234 : }
170 :
171 :
172 :
173 : template <unsigned int Dim, FEFamily T>
174 117481366 : bool FE<Dim,T>::matches_cache(const Elem * elem)
175 : {
176 126897259 : bool m = cached_nodes.size() == elem->n_nodes();
177 166719434 : for (unsigned n = 1; m && n < elem->n_nodes(); n++)
178 52459826 : m = (elem->point(n) - elem->point(0)).relative_fuzzy_equals(cached_nodes[n] - cached_nodes[0]);
179 :
180 117481366 : if (FEInterface::orientation_dependent(T))
181 : {
182 3332063 : m &= cached_edges.size() == elem->n_edges();
183 4551583 : for (unsigned n = 0; m && n < elem->n_edges(); n++)
184 1219520 : m = elem->positive_edge_orientation(n) == cached_edges[n];
185 :
186 3332063 : m &= cached_faces.size() == elem->n_faces();
187 3687811 : for (unsigned n = 0; m && n < elem->n_faces(); n++)
188 355748 : m = elem->positive_face_orientation(n) == cached_faces[n];
189 : }
190 :
191 117481366 : return m;
192 : }
193 :
194 :
195 :
196 : template <unsigned int Dim, FEFamily T>
197 160343635 : 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 160343635 : this->determine_calculations();
207 :
208 : // Try to avoid calling init_shape_functions
209 : // even when shapes_need_reinit
210 13490119 : bool cached_elem_still_fits = false;
211 :
212 : // Most of the hard work happens when we have an actual element
213 160343635 : if (elem)
214 : {
215 : // Initialize the shape functions at the user-specified
216 : // points
217 160343635 : if (pts != nullptr)
218 : {
219 : // Set the type and p level for this element
220 41226766 : this->_elem = elem;
221 41226766 : this->_elem_type = elem->type();
222 41226766 : this->_elem_p_level = elem->p_level();
223 44925168 : this->_p_level = this->_add_p_level_in_reinit * elem->p_level();
224 :
225 : // Initialize the shape functions
226 7396774 : this->_fe_map->template init_reference_to_physical_map<Dim>
227 33829992 : (*pts, elem);
228 41226766 : this->init_shape_functions (*pts, elem);
229 :
230 : // The shape functions do not correspond to the qrule
231 41226766 : 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 9791747 : libmesh_assert(this->qrule);
244 119116869 : this->qrule->init(*elem);
245 :
246 119116869 : 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 227262888 : if (this->get_type() != elem->type() ||
252 117864976 : (elem->runtime_topology() &&
253 108146019 : this->_elem != elem) ||
254 117864976 : this->_elem_p_level != elem->p_level() ||
255 246696927 : !this->shapes_on_quadrature ||
256 19168045 : elem->mapping_type() != LAGRANGE_MAP)
257 : {
258 : // Set the type and p level for this element
259 1635503 : this->_elem = elem;
260 1635503 : this->_elem_type = elem->type();
261 1635503 : this->_elem_p_level = elem->p_level();
262 1744647 : this->_p_level = this->_add_p_level_in_reinit * elem->p_level();
263 :
264 : // Initialize the shape functions
265 328590 : this->_fe_map->template init_reference_to_physical_map<Dim>
266 1635503 : (this->qrule->get_points(), elem);
267 1744647 : this->init_shape_functions (this->qrule->get_points(), elem);
268 : }
269 : else
270 : {
271 117481366 : this->_elem = elem;
272 :
273 : // Check if cached element's nodes, edge and face orientations still fit
274 117481366 : cached_elem_still_fits = this->matches_cache(elem);
275 :
276 : // Initialize the shape functions if needed
277 117481366 : if (this->shapes_need_reinit() && !cached_elem_still_fits)
278 : {
279 978427 : this->_fe_map->template init_reference_to_physical_map<Dim>
280 3898713 : (this->qrule->get_points(), elem);
281 4224864 : 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 119116869 : if (this->shapes_need_reinit() && !cached_elem_still_fits && *caching)
287 4774234 : this->cache(elem);
288 :
289 : // The shape functions correspond to the qrule
290 119116869 : 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 160343635 : if (pts != nullptr)
326 : {
327 41226766 : if (weights != nullptr)
328 : {
329 52000 : this->_fe_map->compute_map (this->dim, *weights, elem, this->calculate_d2phi);
330 : }
331 : else
332 : {
333 48531300 : std::vector<Real> dummy_weights (pts->size(), 1.);
334 41174766 : this->_fe_map->compute_map (this->dim, dummy_weights, elem, this->calculate_d2phi);
335 : }
336 : }
337 : else
338 : {
339 128641906 : 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 160343635 : if (!cached_elem_still_fits)
345 : {
346 117792099 : if (pts != nullptr)
347 41226766 : this->compute_shape_functions (elem,*pts);
348 : else
349 83440609 : this->compute_shape_functions(elem,this->qrule->get_points());
350 117792099 : if (this->calculate_dual)
351 : {
352 : if (T != LAGRANGE)
353 372 : 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 6348 : if (elem && this->calculate_default_dual_coeff)
362 6348 : 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 6348 : this->compute_dual_shape_functions();
366 : }
367 : }
368 160343635 : }
369 :
370 : template <unsigned int Dim, FEFamily T>
371 6348 : 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 6348 : this->_elem = elem;
377 6348 : this->_elem_type = elem->type();
378 6348 : this->_elem_p_level = elem->p_level();
379 6773 : this->_p_level = this->_add_p_level_in_reinit * elem->p_level();
380 :
381 850 : const unsigned int n_shapes =
382 5498 : this->n_dofs(elem, this->get_order());
383 :
384 1275 : std::vector<std::vector<OutputShape>> phi_vals;
385 6348 : phi_vals.resize(n_shapes);
386 108062 : for (const auto i : make_range(phi_vals.size()))
387 115594 : phi_vals[i].resize(pts.size());
388 :
389 6348 : all_shapes(elem, this->get_order(), pts, phi_vals);
390 6348 : this->compute_dual_shape_coeffs(JxW, phi_vals);
391 6348 : }
392 :
393 : template <unsigned int Dim, FEFamily T>
394 6348 : void FE<Dim,T>::reinit_default_dual_shape_coeffs (const Elem * elem)
395 : {
396 425 : libmesh_assert(elem);
397 :
398 425 : FEType default_fe_type(this->get_order(), T);
399 6773 : QGauss default_qrule(elem->dim(), default_fe_type.default_quadrature_order());
400 6348 : 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 6348 : 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 425 : this->set_calculate_default_dual_coeff(false);
407 6348 : }
408 :
409 :
410 : template <unsigned int Dim, FEFamily T>
411 6348 : void FE<Dim,T>::init_dual_shape_functions(const unsigned int n_shapes, const unsigned int n_qp)
412 : {
413 6348 : if (!this->calculate_dual)
414 0 : return;
415 :
416 425 : libmesh_assert_msg(this->calculate_phi,
417 : "dual shape function calculation relies on "
418 : "primal shape functions being calculated");
419 :
420 6348 : this->dual_phi.resize(n_shapes);
421 6348 : if (this->calculate_dphi)
422 6336 : this->dual_dphi.resize(n_shapes);
423 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
424 6348 : if (this->calculate_d2phi)
425 6086 : this->dual_d2phi.resize(n_shapes);
426 : #endif
427 :
428 108062 : for (auto i : index_range(this->dual_phi))
429 : {
430 108654 : this->dual_phi[i].resize(n_qp);
431 101714 : if (this->calculate_dphi)
432 108628 : this->dual_dphi[i].resize(n_qp);
433 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
434 101714 : if (this->calculate_d2phi)
435 106184 : this->dual_d2phi[i].resize(n_qp);
436 : #endif
437 : }
438 : }
439 :
440 : template <unsigned int Dim, FEFamily T>
441 45401603 : 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 8043712 : LOG_SCOPE("init_shape_functions()", "FE");
446 :
447 : // The number of quadrature points.
448 8043176 : const unsigned int n_qp = cast_int<unsigned int>(qp.size());
449 45401603 : this->_n_total_qp = n_qp;
450 :
451 : // Number of shape functions in the finite element approximation
452 : // space.
453 8043176 : const unsigned int n_approx_shape_functions =
454 37358427 : 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 4021856 : 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 45401603 : if (this->calculate_phi)
466 : {
467 42846257 : if (this->phi.size() == n_approx_shape_functions)
468 : {
469 37254985 : old_n_qp = n_approx_shape_functions ? this->phi[0].size() : 0;
470 3301750 : break;
471 : }
472 2142374 : this->phi.resize (n_approx_shape_functions);
473 : }
474 8146618 : if (this->calculate_dphi)
475 : {
476 4432748 : if (this->dphi.size() == n_approx_shape_functions)
477 : {
478 1852770 : old_n_qp = n_approx_shape_functions ? this->dphi[0].size() : 0;
479 213164 : break;
480 : }
481 2209487 : this->dphi.resize (n_approx_shape_functions);
482 2209487 : this->dphidx.resize (n_approx_shape_functions);
483 2209487 : this->dphidy.resize (n_approx_shape_functions);
484 2209487 : this->dphidz.resize (n_approx_shape_functions);
485 : }
486 :
487 6287471 : if (this->calculate_dphiref)
488 : {
489 : if (Dim > 0)
490 : {
491 3638141 : if (this->dphidxi.size() == n_approx_shape_functions)
492 : {
493 34935 : old_n_qp = n_approx_shape_functions ? this->dphidxi[0].size() : 0;
494 2985 : break;
495 : }
496 3348533 : this->dphidxi.resize (n_approx_shape_functions);
497 : }
498 :
499 : if (Dim > 1)
500 3052920 : this->dphideta.resize (n_approx_shape_functions);
501 :
502 : if (Dim > 2)
503 2384054 : this->dphidzeta.resize (n_approx_shape_functions);
504 : }
505 :
506 6258913 : if (this->calculate_curl_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
507 7023 : this->curl_phi.resize(n_approx_shape_functions);
508 :
509 6258913 : if (this->calculate_div_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
510 11185 : this->div_phi.resize(n_approx_shape_functions);
511 :
512 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
513 6258913 : if (this->calculate_d2phi)
514 : {
515 1586721 : 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 1483387 : this->d2phi.resize (n_approx_shape_functions);
522 1483387 : this->d2phidx2.resize (n_approx_shape_functions);
523 1483387 : this->d2phidxdy.resize (n_approx_shape_functions);
524 1483387 : this->d2phidxdz.resize (n_approx_shape_functions);
525 1483387 : this->d2phidy2.resize (n_approx_shape_functions);
526 1483387 : this->d2phidydz.resize (n_approx_shape_functions);
527 1483387 : this->d2phidz2.resize (n_approx_shape_functions);
528 :
529 : if (Dim > 0)
530 1483387 : this->d2phidxi2.resize (n_approx_shape_functions);
531 :
532 : if (Dim > 1)
533 : {
534 1197829 : this->d2phidxideta.resize (n_approx_shape_functions);
535 1197829 : this->d2phideta2.resize (n_approx_shape_functions);
536 : }
537 : if (Dim > 2)
538 : {
539 1164163 : this->d2phidxidzeta.resize (n_approx_shape_functions);
540 1164163 : this->d2phidetadzeta.resize (n_approx_shape_functions);
541 1164163 : this->d2phidzeta2.resize (n_approx_shape_functions);
542 : }
543 : }
544 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
545 : }
546 : while (false);
547 :
548 45401603 : if (old_n_qp != n_qp)
549 55582335 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
550 : {
551 49246488 : if (this->calculate_phi)
552 20012351 : this->phi[i].resize (n_qp);
553 :
554 49246488 : if (this->calculate_dphi)
555 : {
556 18327160 : this->dphi[i].resize (n_qp);
557 18327160 : this->dphidx[i].resize (n_qp);
558 18327160 : this->dphidy[i].resize (n_qp);
559 18327160 : this->dphidz[i].resize (n_qp);
560 : }
561 :
562 49240111 : if (this->calculate_dphiref)
563 : {
564 : if (Dim > 0)
565 28423488 : this->dphidxi[i].resize(n_qp);
566 :
567 : if (Dim > 1)
568 27208577 : this->dphideta[i].resize(n_qp);
569 :
570 : if (Dim > 2)
571 22650448 : this->dphidzeta[i].resize(n_qp);
572 : }
573 :
574 49246488 : if (this->calculate_curl_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
575 58194 : this->curl_phi[i].resize(n_qp);
576 :
577 49246488 : if (this->calculate_div_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
578 92816 : this->div_phi[i].resize(n_qp);
579 :
580 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
581 49246488 : if (this->calculate_d2phi)
582 : {
583 13134965 : this->d2phi[i].resize (n_qp);
584 13134965 : this->d2phidx2[i].resize (n_qp);
585 13134965 : this->d2phidxdy[i].resize (n_qp);
586 13134965 : this->d2phidxdz[i].resize (n_qp);
587 13134965 : this->d2phidy2[i].resize (n_qp);
588 13134965 : this->d2phidydz[i].resize (n_qp);
589 13134965 : this->d2phidz2[i].resize (n_qp);
590 : if (Dim > 0)
591 13134965 : this->d2phidxi2[i].resize (n_qp);
592 : if (Dim > 1)
593 : {
594 11966378 : this->d2phidxideta[i].resize (n_qp);
595 11966378 : this->d2phideta2[i].resize (n_qp);
596 : }
597 : if (Dim > 2)
598 : {
599 11606097 : this->d2phidxidzeta[i].resize (n_qp);
600 11606097 : this->d2phidetadzeta[i].resize (n_qp);
601 11606097 : 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 11590807 : if (this->calculate_phi || this->calculate_dphi)
617 : {
618 10846151 : this->weight.resize (n_qp);
619 46769968 : for (unsigned int p=0; p<n_qp; p++)
620 48562168 : this->weight[p] = 1.;
621 : }
622 :
623 11590807 : if (this->calculate_dphi)
624 : {
625 4257362 : this->dweight.resize (n_qp);
626 4257362 : this->dphase.resize (n_qp);
627 15841385 : for (unsigned int p=0; p<n_qp; p++)
628 : {
629 11584023 : this->dweight[p].zero();
630 8524301 : 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 45202823 : if (this->calculate_dphiref && Dim > 0)
638 : {
639 23502630 : std::vector<std::vector<OutputShape>> * comps[3]
640 23502630 : { &this->dphidxi, &this->dphideta, &this->dphidzeta };
641 17004056 : 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 2509323 : case 1:
657 : {
658 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
659 : // Compute the value of shape function i Hessians at quadrature point p
660 2736774 : if (this->calculate_d2phi)
661 1504803 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
662 5891882 : for (unsigned int p=0; p<n_qp; p++)
663 4689280 : this->d2phidxi2[i][p] = FE<Dim, T>::shape_second_deriv(
664 4689280 : elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
665 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
666 :
667 227451 : break;
668 : }
669 :
670 :
671 :
672 : //------------------------------------------------------------
673 : // 2D
674 17483924 : case 2:
675 : {
676 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
677 : // Compute the value of shape function i Hessians at quadrature point p
678 19327293 : if (this->calculate_d2phi)
679 9246335 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
680 55429666 : for (unsigned int p=0; p<n_qp; p++)
681 : {
682 47034584 : this->d2phidxi2[i][p] = FE<Dim, T>::shape_second_deriv(
683 47034584 : elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
684 87013326 : this->d2phidxideta[i][p] = FE<Dim, T>::shape_second_deriv(
685 47034584 : elem, this->fe_type.order, i, 1, qp[p], this->_add_p_level_in_reinit);
686 47034584 : this->d2phideta2[i][p] = FE<Dim, T>::shape_second_deriv(
687 47034584 : 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 1843369 : break;
693 : }
694 :
695 :
696 :
697 : //------------------------------------------------------------
698 : // 3D
699 21206769 : case 3:
700 : {
701 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
702 : // Compute the value of shape function i Hessians at quadrature point p
703 23138756 : if (this->calculate_d2phi)
704 123459180 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
705 391609763 : for (unsigned int p=0; p<n_qp; p++)
706 : {
707 273029274 : this->d2phidxi2[i][p] = FE<Dim, T>::shape_second_deriv(
708 273029274 : elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
709 501090664 : this->d2phidxideta[i][p] = FE<Dim, T>::shape_second_deriv(
710 273029274 : elem, this->fe_type.order, i, 1, qp[p], this->_add_p_level_in_reinit);
711 501090664 : this->d2phideta2[i][p] = FE<Dim, T>::shape_second_deriv(
712 273029274 : elem, this->fe_type.order, i, 2, qp[p], this->_add_p_level_in_reinit);
713 501090664 : this->d2phidxidzeta[i][p] = FE<Dim, T>::shape_second_deriv(
714 273029274 : elem, this->fe_type.order, i, 3, qp[p], this->_add_p_level_in_reinit);
715 501090664 : this->d2phidetadzeta[i][p] = FE<Dim, T>::shape_second_deriv(
716 273029274 : elem, this->fe_type.order, i, 4, qp[p], this->_add_p_level_in_reinit);
717 273029274 : this->d2phidzeta2[i][p] = FE<Dim, T>::shape_second_deriv(
718 273029274 : 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 1931987 : break;
723 : }
724 :
725 :
726 : default:
727 : libmesh_error_msg("Invalid dimension Dim = " << Dim);
728 : }
729 :
730 45401603 : if (this->calculate_dual)
731 5394 : this->init_dual_shape_functions(n_approx_shape_functions, n_qp);
732 45401603 : }
733 :
734 : template <unsigned int Dim, FEFamily T>
735 : void
736 43677192 : 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 152552174 : for (unsigned int d=0; d != Dim; ++d)
743 : {
744 108874982 : auto & comps_d = *comps[d];
745 1425361080 : for (auto i : index_range(comps_d))
746 : FE<Dim,T>::shape_derivs
747 1430485086 : (elem,o,i,d,p,comps_d[i],add_p_level);
748 : }
749 43677192 : }
750 :
751 :
752 : template <unsigned int Dim, FEFamily T>
753 : void
754 3216 : 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 3216 : nodal_soln(elem, o, elem_soln, full_nodal_soln, add_p_level, vdim);
763 3216 : const std::vector<unsigned int> side_nodes =
764 3216 : elem->nodes_on_side(side);
765 :
766 536 : std::size_t n_side_nodes = side_nodes.size();
767 3216 : nodal_soln_on_side.resize(n_side_nodes);
768 22944 : for (auto n : make_range(n_side_nodes))
769 24660 : nodal_soln_on_side[n] = full_nodal_soln[side_nodes[n]];
770 3216 : }
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 3483800262 : fdm_points(const unsigned int j, const Point & p)
796 : {
797 281635446 : libmesh_assert_less (j, LIBMESH_DIM);
798 :
799 : // cheat by using finite difference approximations:
800 281635446 : const Real eps = 1.e-6;
801 3483800262 : Point pp = p, pm = p;
802 :
803 3483800262 : switch (j)
804 : {
805 : // d()/dxi
806 131199192 : case 0:
807 : {
808 1612007960 : pp(0) += eps;
809 1612007960 : pm(0) -= eps;
810 1612007960 : break;
811 : }
812 :
813 : // d()/deta
814 94808674 : case 1:
815 : {
816 1172365656 : pp(1) += eps;
817 1172365656 : pm(1) -= eps;
818 1172365656 : break;
819 : }
820 :
821 : // d()/dzeta
822 55627580 : case 2:
823 : {
824 699426646 : pp(2) += eps;
825 699426646 : pm(2) -= eps;
826 699426646 : break;
827 : }
828 :
829 0 : default:
830 0 : libmesh_error_msg("Invalid derivative index j = " << j);
831 : }
832 :
833 3765435708 : return std::make_tuple(pp, pm, eps);
834 : }
835 :
836 : std::tuple<Point, Point, Real, unsigned int>
837 1585743489 : fdm_second_points(const unsigned int j, const Point & p)
838 : {
839 : // cheat by using finite difference approximations:
840 131324391 : const Real eps = 1.e-5;
841 1585743489 : Point pp = p, pm = p;
842 131324391 : unsigned int deriv_j = 0;
843 :
844 1585743489 : switch (j)
845 : {
846 : // d^2() / dxi^2
847 22461599 : case 0:
848 : {
849 271227443 : pp(0) += eps;
850 271227443 : pm(0) -= eps;
851 22461599 : deriv_j = 0;
852 271227443 : break;
853 : }
854 :
855 : // d^2() / dxi deta
856 22461599 : case 1:
857 : {
858 271227443 : pp(1) += eps;
859 271227443 : pm(1) -= eps;
860 22461599 : deriv_j = 0;
861 271227443 : break;
862 : }
863 :
864 : // d^2() / deta^2
865 22461599 : case 2:
866 : {
867 271227443 : pp(1) += eps;
868 271227443 : pm(1) -= eps;
869 22461599 : deriv_j = 1;
870 271227443 : break;
871 : }
872 :
873 : // d^2()/dxidzeta
874 21313198 : case 3:
875 : {
876 257353720 : pp(2) += eps;
877 257353720 : pm(2) -= eps;
878 21313198 : deriv_j = 0;
879 257353720 : break;
880 : } // d^2()/deta^2
881 :
882 : // d^2()/detadzeta
883 21313198 : case 4:
884 : {
885 257353720 : pp(2) += eps;
886 257353720 : pm(2) -= eps;
887 21313198 : deriv_j = 1;
888 257353720 : break;
889 : }
890 :
891 : // d^2()/dzeta^2
892 21313198 : case 5:
893 : {
894 257353720 : pp(2) += eps;
895 257353720 : pm(2) -= eps;
896 21313198 : deriv_j = 2;
897 257353720 : break;
898 : }
899 :
900 0 : default:
901 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
902 : }
903 :
904 1717067880 : return std::make_tuple(pp, pm, eps, deriv_j);
905 : }
906 :
907 : }
908 :
909 :
910 : template <typename OutputShape>
911 3335220990 : 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 276680352 : libmesh_assert(elem);
923 :
924 3335220990 : auto [pp, pm, eps] = fdm_points(j, p);
925 :
926 3335220990 : return (shape_func(elem, order, i, pp, add_p_level) -
927 3335220990 : 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 148579272 : 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 148579272 : auto [pp, pm, eps] = fdm_points(j, p);
961 :
962 148579272 : return (shape_func(type, order, elem, i, pp) -
963 148579272 : shape_func(type, order, elem, i, pm))/2./eps;
964 : }
965 :
966 :
967 : template <typename OutputShape>
968 : OutputShape
969 1578358017 : 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 1578358017 : auto [pp, pm, eps, deriv_j] = fdm_second_points(j, p);
981 :
982 1578358017 : return (deriv_func(elem, order, i, deriv_j, pp, add_p_level) -
983 1578358017 : 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 7385472 : 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 7385472 : auto [pp, pm, eps, deriv_j] = fdm_second_points(j, p);
1021 :
1022 7385472 : return (deriv_func(type, order, elem, i, deriv_j, pp) -
1023 7385472 : deriv_func(type, order, elem, i, deriv_j, pm))/2./eps;
1024 : }
1025 :
1026 :
1027 447222 : 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 447222 : const int extra_order = add_p_level * elem->p_level();
1034 :
1035 447222 : const int dim = elem->dim();
1036 :
1037 : const unsigned int n_sf =
1038 447222 : FEInterface::n_shape_functions(underlying_fe_type, extra_order,
1039 38757 : elem);
1040 :
1041 38757 : libmesh_assert_equal_to (n_sf, elem->n_nodes());
1042 :
1043 485979 : std::vector<Real> node_weights(n_sf);
1044 :
1045 77514 : const unsigned char datum_index = elem->mapping_data();
1046 5615830 : for (unsigned int n=0; n<n_sf; n++)
1047 5615965 : node_weights[n] =
1048 5168608 : elem->node_ref(n).get_extra_datum<Real>(datum_index);
1049 :
1050 77514 : const std::size_t n_p = p.size();
1051 :
1052 447222 : shapes.resize(n_sf);
1053 5615830 : for (unsigned int i=0; i != n_sf; ++i)
1054 : {
1055 5168608 : auto & shapes_i = shapes[i];
1056 5168608 : shapes_i.resize(n_p, 0);
1057 5168608 : FEInterface::shapes(dim, underlying_fe_type, elem, i, p,
1058 : shapes_i, add_p_level);
1059 27693992 : for (auto & s : shapes_i)
1060 24557857 : s *= node_weights[i];
1061 : }
1062 447222 : }
1063 :
1064 :
1065 561870 : 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 561870 : const int extra_order = add_p_level * elem->p_level();
1073 561870 : const unsigned int dim = elem->dim();
1074 :
1075 : const unsigned int n_sf =
1076 561870 : FEInterface::n_shape_functions(fe_type, extra_order, elem);
1077 :
1078 49896 : libmesh_assert_equal_to (n_sf, elem->n_nodes());
1079 :
1080 49896 : libmesh_assert_equal_to (dim, derivs.size());
1081 1846502 : for (unsigned int d = 0; d != dim; ++d)
1082 1396520 : derivs[d].resize(n_sf);
1083 :
1084 611766 : std::vector<Real> node_weights(n_sf);
1085 :
1086 99792 : const unsigned char datum_index = elem->mapping_data();
1087 6745720 : for (unsigned int n=0; n<n_sf; n++)
1088 6731265 : node_weights[n] =
1089 6183850 : elem->node_ref(n).get_extra_datum<Real>(datum_index);
1090 :
1091 99792 : const std::size_t n_p = p.size();
1092 :
1093 561870 : shapes.resize(n_sf);
1094 6745720 : for (unsigned int i=0; i != n_sf; ++i)
1095 6731265 : shapes[i].resize(n_p, 0);
1096 :
1097 561870 : FEInterface::all_shapes(dim, fe_type, elem, p, shapes, add_p_level);
1098 :
1099 6745720 : for (unsigned int i=0; i != n_sf; ++i)
1100 : {
1101 6183850 : auto & shapes_i = shapes[i];
1102 :
1103 42928406 : for (auto & s : shapes_i)
1104 40114047 : s *= node_weights[i];
1105 :
1106 21210010 : for (unsigned int d = 0; d != dim; ++d)
1107 : {
1108 16335434 : auto & derivs_di = derivs[d][i];
1109 15026160 : derivs_di.resize(n_p);
1110 15026160 : FEInterface::shape_derivs(fe_type, elem, i, d, p,
1111 : derivs_di, add_p_level);
1112 96961110 : for (auto & dip : derivs_di)
1113 89357066 : dip *= node_weights[i];
1114 : }
1115 : }
1116 561870 : }
1117 :
1118 :
1119 23320440 : 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 23320440 : int extra_order = add_p_level * elem.p_level();
1126 :
1127 : const unsigned int n_sf =
1128 23320440 : FEInterface::n_shape_functions(underlying_fe_type, extra_order, &elem);
1129 :
1130 2116962 : libmesh_assert_equal_to (n_sf, elem.n_nodes());
1131 :
1132 23320440 : std::vector<Real> node_weights(n_sf);
1133 :
1134 4107474 : const unsigned char datum_index = elem.mapping_data();
1135 :
1136 2116962 : Real weighted_shape_i = 0, weighted_sum = 0;
1137 :
1138 318152052 : for (unsigned int sf=0; sf<n_sf; sf++)
1139 : {
1140 : Real node_weight =
1141 294831612 : elem.node_ref(sf).get_extra_datum<Real>(datum_index);
1142 : Real weighted_shape = node_weight *
1143 294831612 : FEInterface::shape(underlying_fe_type, extra_order, &elem, sf, p);
1144 294831612 : weighted_sum += weighted_shape;
1145 294831612 : if (sf == i)
1146 2116962 : weighted_shape_i = weighted_shape;
1147 : }
1148 :
1149 27427914 : return weighted_shape_i / weighted_sum;
1150 : }
1151 :
1152 :
1153 49542084 : 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 4181172 : libmesh_assert_less(j, elem.dim());
1161 :
1162 49542084 : int extra_order = add_p_level * elem.p_level();
1163 :
1164 : const unsigned int n_sf =
1165 49542084 : FEInterface::n_shape_functions(underlying_fe_type, extra_order, &elem);
1166 :
1167 49542084 : const unsigned int n_nodes = elem.n_nodes();
1168 4181172 : libmesh_assert_equal_to (n_sf, n_nodes);
1169 :
1170 49542084 : std::vector<Real> node_weights(n_nodes);
1171 :
1172 8362344 : const unsigned char datum_index = elem.mapping_data();
1173 754581492 : for (unsigned int n=0; n<n_nodes; n++)
1174 761551656 : node_weights[n] =
1175 705039408 : elem.node_ref(n).get_extra_datum<Real>(datum_index);
1176 :
1177 4181172 : Real weighted_shape_i = 0, weighted_sum = 0,
1178 4181172 : weighted_grad_i = 0, weighted_grad_sum = 0;
1179 :
1180 754581492 : for (unsigned int sf=0; sf<n_sf; sf++)
1181 : {
1182 705039408 : Real weighted_shape = node_weights[sf] *
1183 705039408 : FEInterface::shape(underlying_fe_type, extra_order, &elem, sf, p);
1184 705039408 : Real weighted_grad = node_weights[sf] *
1185 705039408 : FEInterface::shape_deriv(underlying_fe_type, extra_order, &elem, sf, j, p);
1186 705039408 : weighted_sum += weighted_shape;
1187 705039408 : weighted_grad_sum += weighted_grad;
1188 705039408 : if (sf == i)
1189 : {
1190 4181172 : weighted_shape_i = weighted_shape;
1191 4181172 : weighted_grad_i = weighted_grad;
1192 : }
1193 : }
1194 :
1195 49542084 : return (weighted_sum * weighted_grad_i - weighted_shape_i * weighted_grad_sum) /
1196 57904428 : weighted_sum / weighted_sum;
1197 : }
1198 :
1199 :
1200 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1201 :
1202 13696218 : 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 1111222 : switch (j)
1211 : {
1212 189701 : case 0:
1213 : // j = 0 ==> d^2 phi / dxi^2
1214 189701 : j1 = j2 = 0;
1215 189701 : break;
1216 189412 : case 1:
1217 : // j = 1 ==> d^2 phi / dxi deta
1218 189412 : j1 = 0;
1219 189412 : j2 = 1;
1220 189412 : break;
1221 189412 : case 2:
1222 : // j = 2 ==> d^2 phi / deta^2
1223 189412 : j1 = j2 = 1;
1224 189412 : break;
1225 180899 : case 3:
1226 : // j = 3 ==> d^2 phi / dxi dzeta
1227 180899 : j1 = 0;
1228 180899 : j2 = 2;
1229 180899 : break;
1230 180899 : case 4:
1231 : // j = 4 ==> d^2 phi / deta dzeta
1232 180899 : j1 = 1;
1233 180899 : j2 = 2;
1234 180899 : break;
1235 180899 : case 5:
1236 : // j = 5 ==> d^2 phi / dzeta^2
1237 180899 : j1 = j2 = 2;
1238 180899 : break;
1239 0 : default:
1240 0 : libmesh_error();
1241 : }
1242 :
1243 13696218 : int extra_order = add_p_level * elem.p_level();
1244 :
1245 : const unsigned int n_sf =
1246 13696218 : FEInterface::n_shape_functions(underlying_fe_type, extra_order,
1247 1111222 : &elem);
1248 :
1249 13696218 : const unsigned int n_nodes = elem.n_nodes();
1250 1111222 : libmesh_assert_equal_to (n_sf, n_nodes);
1251 :
1252 13696218 : std::vector<Real> node_weights(n_nodes);
1253 :
1254 2222444 : const unsigned char datum_index = elem.mapping_data();
1255 308030160 : for (unsigned int n=0; n<n_nodes; n++)
1256 318740822 : node_weights[n] =
1257 294333942 : elem.node_ref(n).get_extra_datum<Real>(datum_index);
1258 :
1259 1111222 : Real weighted_shape_i = 0, weighted_sum = 0,
1260 1111222 : weighted_grada_i = 0, weighted_grada_sum = 0,
1261 1111222 : weighted_gradb_i = 0, weighted_gradb_sum = 0,
1262 1111222 : weighted_hess_i = 0, weighted_hess_sum = 0;
1263 :
1264 308030160 : for (unsigned int sf=0; sf<n_sf; sf++)
1265 : {
1266 294333942 : Real weighted_shape = node_weights[sf] *
1267 294333942 : FEInterface::shape(underlying_fe_type, extra_order, &elem, sf,
1268 294333942 : p);
1269 294333942 : Real weighted_grada = node_weights[sf] *
1270 294333942 : FEInterface::shape_deriv(underlying_fe_type, extra_order,
1271 294333942 : &elem, sf, j1, p);
1272 294333942 : Real weighted_hess = node_weights[sf] *
1273 294333942 : FEInterface::shape_second_deriv(underlying_fe_type,
1274 294333942 : extra_order, &elem, sf, j, p);
1275 294333942 : weighted_sum += weighted_shape;
1276 294333942 : weighted_grada_sum += weighted_grada;
1277 24406880 : Real weighted_gradb = weighted_grada;
1278 294333942 : if (j1 != j2)
1279 : {
1280 146828754 : weighted_gradb = (j1 == j2) ? weighted_grada :
1281 146828754 : node_weights[sf] *
1282 146828754 : FEInterface::shape_deriv(underlying_fe_type, extra_order,
1283 : &elem, sf, j2, p);
1284 146828754 : weighted_grada_sum += weighted_grada;
1285 : }
1286 294333942 : weighted_hess_sum += weighted_hess;
1287 294333942 : if (sf == i)
1288 : {
1289 1111222 : weighted_shape_i = weighted_shape;
1290 1111222 : weighted_grada_i = weighted_grada;
1291 1111222 : weighted_gradb_i = weighted_gradb;
1292 1111222 : weighted_hess_i = weighted_hess;
1293 : }
1294 : }
1295 :
1296 13696218 : if (j1 == j2)
1297 560012 : weighted_gradb_sum = weighted_grada_sum;
1298 :
1299 15918662 : return (weighted_sum * weighted_hess_i - weighted_grada_i * weighted_gradb_sum -
1300 18141106 : weighted_shape_i * weighted_hess_sum - weighted_gradb_i * weighted_grada_sum +
1301 15918662 : 2 * weighted_grada_sum * weighted_shape_i * weighted_gradb_sum / weighted_sum) /
1302 15918662 : weighted_sum / weighted_sum;
1303 : }
1304 :
1305 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1306 :
1307 :
1308 447222 : 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 116271 : std::vector<std::vector<Real>> shapes;
1315 :
1316 447222 : rational_fe_weighted_shapes(&elem, underlying_fe_type, shapes, p,
1317 : add_p_level);
1318 :
1319 524736 : std::vector<Real> shape_sums(p.size(), 0);
1320 :
1321 5615830 : for (auto i : index_range(v))
1322 : {
1323 447357 : libmesh_assert_equal_to ( p.size(), shapes[i].size() );
1324 27693992 : for (auto j : index_range(p))
1325 28622803 : shape_sums[j] += shapes[i][j];
1326 : }
1327 :
1328 5615830 : for (auto i : index_range(v))
1329 : {
1330 447357 : libmesh_assert_equal_to ( p.size(), v[i].size() );
1331 27693992 : for (auto j : index_range(v[i]))
1332 30655276 : v[i][j] = shapes[i][j] / shape_sums[j];
1333 : }
1334 447222 : }
1335 :
1336 :
1337 : template <typename OutputShape>
1338 561870 : 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 561870 : const int my_dim = elem.dim();
1345 :
1346 149688 : std::vector<std::vector<Real>> shapes;
1347 661662 : std::vector<std::vector<std::vector<Real>>> derivs(my_dim);
1348 :
1349 561870 : rational_fe_weighted_shapes_derivs(&elem, underlying_fe_type,
1350 : shapes, derivs, p, add_p_level);
1351 :
1352 661662 : std::vector<Real> shape_sums(p.size(), 0);
1353 661662 : std::vector<std::vector<Real>> shape_deriv_sums(my_dim);
1354 1846502 : for (int d=0; d != my_dim; ++d)
1355 1508408 : shape_deriv_sums[d].resize(p.size());
1356 :
1357 6745720 : for (auto i : index_range(shapes))
1358 : {
1359 547415 : libmesh_assert_equal_to ( p.size(), shapes[i].size() );
1360 42928406 : for (auto j : index_range(p))
1361 43483538 : shape_sums[j] += shapes[i][j];
1362 :
1363 21210010 : for (int d=0; d != my_dim; ++d)
1364 96961110 : for (auto j : index_range(p))
1365 119045530 : shape_deriv_sums[d][j] += derivs[d][i][j];
1366 : }
1367 :
1368 1846502 : for (int d=0; d != my_dim; ++d)
1369 : {
1370 1284632 : auto & comps_d = *comps[d];
1371 111888 : libmesh_assert_equal_to(comps_d.size(), elem.n_nodes());
1372 :
1373 16310792 : for (auto i : index_range(comps_d))
1374 : {
1375 1309274 : auto & comps_di = comps_d[i];
1376 3927822 : auto & derivs_di = derivs[d][i];
1377 :
1378 96961110 : for (auto j : index_range(comps_di))
1379 111623414 : comps_di[j] = (shape_sums[j] * derivs_di[j] -
1380 96779182 : shapes[i][j] * shape_deriv_sums[d][j]) /
1381 81934950 : shape_sums[j] / shape_sums[j];
1382 : }
1383 : }
1384 1023948 : }
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
|
