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 375 : void nonlagrange_dual_warning () {
37 : libmesh_warning("dual calculations have only been verified for the LAGRANGE family");
38 375 : }
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 2092 : void setup() { caching = new bool(!on_command_line("--disable-caching")); }
54 : public:
55 1624 : ~CachingSetup() { delete caching; caching = nullptr; }
56 : } caching_setup;
57 :
58 :
59 : // ------------------------------------------------------------
60 : // FE class members
61 : template <unsigned int Dim, FEFamily T>
62 2960044 : FE<Dim,T>::FE (const FEType & fet) :
63 : FEGenericBase<typename FEOutputType<T>::type> (Dim,fet),
64 1305937 : last_side(INVALID_ELEM),
65 2960044 : 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 827359 : libmesh_assert_equal_to (T, this->get_family());
71 2960044 : }
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 1868939 : void FE<Dim,T>::attach_quadrature_rule (QBase * q)
88 : {
89 548761 : libmesh_assert(q);
90 1868939 : this->qrule = q;
91 : // make sure we don't cache results from a previous quadrature rule
92 1868939 : this->_elem = nullptr;
93 1868939 : this->_elem_type = INVALID_ELEM;
94 1868939 : return;
95 : }
96 :
97 :
98 : template <unsigned int Dim, FEFamily T>
99 266257 : 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 76999 : libmesh_assert(elem);
106 76999 : libmesh_assert_less (s, elem->n_sides());
107 :
108 76999 : di.clear();
109 76999 : unsigned int nodenum = 0;
110 266257 : const unsigned int n_nodes = elem->n_nodes();
111 3108912 : for (unsigned int n = 0; n != n_nodes; ++n)
112 : {
113 : const unsigned int n_dofs =
114 4421277 : n_dofs_at_node(*elem, static_cast<Order>(o + add_p_level*elem->p_level()), n);
115 2842655 : if (elem->is_node_on_side(n, s))
116 2711307 : for (unsigned int i = 0; i != n_dofs; ++i)
117 1526880 : di.push_back(nodenum++);
118 : else
119 1658228 : nodenum += n_dofs;
120 : }
121 266257 : }
122 :
123 :
124 :
125 : template <unsigned int Dim, FEFamily T>
126 33900 : 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 33900 : const unsigned int n_nodes = elem->n_nodes();
138 720768 : for (unsigned int n = 0; n != n_nodes; ++n)
139 : {
140 : const unsigned int n_dofs =
141 1038132 : n_dofs_at_node(*elem, static_cast<Order>(o + add_p_level*elem->p_level()), n);
142 686868 : if (elem->is_node_on_edge(n, e))
143 232614 : for (unsigned int i = 0; i != n_dofs; ++i)
144 133230 : di.push_back(nodenum++);
145 : else
146 587484 : nodenum += n_dofs;
147 : }
148 33900 : }
149 :
150 :
151 :
152 : template <unsigned int Dim, FEFamily T>
153 1472622 : void FE<Dim,T>::cache(const Elem * elem)
154 : {
155 1472622 : cached_nodes.resize(elem->n_nodes());
156 17590194 : for (auto n : elem->node_index_range())
157 23969276 : cached_nodes[n] = elem->point(n);
158 :
159 1472622 : if (FEInterface::orientation_dependent(T))
160 : {
161 1179815 : cached_edges.resize(elem->n_edges());
162 6989940 : for (auto n : elem->edge_index_range())
163 5810125 : cached_edges[n] = elem->positive_edge_orientation(n);
164 :
165 1179815 : cached_faces.resize(elem->n_faces());
166 3621131 : for (auto n : elem->face_index_range())
167 2441316 : cached_faces[n] = elem->positive_face_orientation(n);
168 : }
169 1472622 : }
170 :
171 :
172 :
173 : template <unsigned int Dim, FEFamily T>
174 32653612 : bool FE<Dim,T>::matches_cache(const Elem * elem)
175 : {
176 42154171 : bool m = cached_nodes.size() == elem->n_nodes();
177 45839465 : for (unsigned n = 1; m && n < elem->n_nodes(); n++)
178 16735985 : m = (elem->point(n) - elem->point(0)).relative_fuzzy_equals(cached_nodes[n] - cached_nodes[0]);
179 :
180 32653612 : if (FEInterface::orientation_dependent(T))
181 : {
182 1213420 : m &= cached_edges.size() == elem->n_edges();
183 2183053 : for (unsigned n = 0; m && n < elem->n_edges(); n++)
184 969633 : m = elem->positive_edge_orientation(n) == cached_edges[n];
185 :
186 1213420 : m &= cached_faces.size() == elem->n_faces();
187 1321288 : for (unsigned n = 0; m && n < elem->n_faces(); n++)
188 107868 : m = elem->positive_face_orientation(n) == cached_faces[n];
189 : }
190 :
191 32653612 : return m;
192 : }
193 :
194 :
195 :
196 : template <unsigned int Dim, FEFamily T>
197 48463511 : 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 48463511 : this->determine_calculations();
207 :
208 : // Try to avoid calling init_shape_functions
209 : // even when shapes_need_reinit
210 14038357 : bool cached_elem_still_fits = false;
211 :
212 : // Most of the hard work happens when we have an actual element
213 48463511 : if (elem)
214 : {
215 : // Initialize the shape functions at the user-specified
216 : // points
217 48463511 : if (pts != nullptr)
218 : {
219 : // Set the type and p level for this element
220 15414822 : this->_elem = elem;
221 15414822 : this->_elem_type = elem->type();
222 15414822 : this->_elem_p_level = elem->p_level();
223 19575857 : this->_p_level = this->_add_p_level_in_reinit * elem->p_level();
224 :
225 : // Initialize the shape functions
226 8322070 : this->_fe_map->template init_reference_to_physical_map<Dim>
227 7092752 : (*pts, elem);
228 15414822 : this->init_shape_functions (*pts, elem);
229 :
230 : // The shape functions do not correspond to the qrule
231 15414822 : 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 9877322 : libmesh_assert(this->qrule);
244 33048689 : this->qrule->init(*elem);
245 :
246 33048689 : 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 56033078 : if (this->get_type() != elem->type() ||
252 32787926 : (elem->runtime_topology() &&
253 22984389 : this->_elem != elem) ||
254 32787926 : this->_elem_p_level != elem->p_level() ||
255 75637588 : !this->shapes_on_quadrature ||
256 19337218 : elem->mapping_type() != LAGRANGE_MAP)
257 : {
258 : // Set the type and p level for this element
259 395077 : this->_elem = elem;
260 395077 : this->_elem_type = elem->type();
261 395077 : this->_elem_p_level = elem->p_level();
262 505266 : this->_p_level = this->_add_p_level_in_reinit * elem->p_level();
263 :
264 : // Initialize the shape functions
265 331787 : this->_fe_map->template init_reference_to_physical_map<Dim>
266 395077 : (this->qrule->get_points(), elem);
267 505266 : this->init_shape_functions (this->qrule->get_points(), elem);
268 : }
269 : else
270 : {
271 32653612 : this->_elem = elem;
272 :
273 : // Check if cached element's nodes, edge and face orientations still fit
274 32653612 : cached_elem_still_fits = this->matches_cache(elem);
275 :
276 : // Initialize the shape functions if needed
277 32653612 : if (this->shapes_need_reinit() && !cached_elem_still_fits)
278 : {
279 913233 : this->_fe_map->template init_reference_to_physical_map<Dim>
280 1244474 : (this->qrule->get_points(), elem);
281 1548885 : 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 33048689 : if (this->shapes_need_reinit() && !cached_elem_still_fits && *caching)
287 1472622 : this->cache(elem);
288 :
289 : // The shape functions correspond to the qrule
290 33048689 : 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 48463511 : if (pts != nullptr)
326 : {
327 15414822 : if (weights != nullptr)
328 : {
329 42604 : this->_fe_map->compute_map (this->dim, *weights, elem, this->calculate_d2phi);
330 : }
331 : else
332 : {
333 23654042 : std::vector<Real> dummy_weights (pts->size(), 1.);
334 15372218 : this->_fe_map->compute_map (this->dim, dummy_weights, elem, this->calculate_d2phi);
335 : }
336 : }
337 : else
338 : {
339 42659437 : 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 48463511 : if (!cached_elem_still_fits)
345 : {
346 38399619 : if (pts != nullptr)
347 15414822 : this->compute_shape_functions (elem,*pts);
348 : else
349 29893892 : this->compute_shape_functions(elem,this->qrule->get_points());
350 38399619 : if (this->calculate_dual)
351 : {
352 : if (T != LAGRANGE)
353 375 : 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 1562 : if (elem && this->calculate_default_dual_coeff)
362 1562 : 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 1562 : this->compute_dual_shape_functions();
366 : }
367 : }
368 48463511 : }
369 :
370 : template <unsigned int Dim, FEFamily T>
371 1562 : 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 1562 : this->_elem = elem;
377 1562 : this->_elem_type = elem->type();
378 1562 : this->_elem_p_level = elem->p_level();
379 1990 : this->_p_level = this->_add_p_level_in_reinit * elem->p_level();
380 :
381 856 : const unsigned int n_shapes =
382 706 : this->n_dofs(elem, this->get_order());
383 :
384 1284 : std::vector<std::vector<OutputShape>> phi_vals;
385 1562 : phi_vals.resize(n_shapes);
386 26901 : for (const auto i : make_range(phi_vals.size()))
387 39255 : phi_vals[i].resize(pts.size());
388 :
389 1562 : all_shapes(elem, this->get_order(), pts, phi_vals);
390 1562 : this->compute_dual_shape_coeffs(JxW, phi_vals);
391 1562 : }
392 :
393 : template <unsigned int Dim, FEFamily T>
394 1562 : void FE<Dim,T>::reinit_default_dual_shape_coeffs (const Elem * elem)
395 : {
396 428 : libmesh_assert(elem);
397 :
398 428 : FEType default_fe_type(this->get_order(), T);
399 1990 : QGauss default_qrule(elem->dim(), default_fe_type.default_quadrature_order());
400 1562 : 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 1562 : 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 428 : this->set_calculate_default_dual_coeff(false);
407 1562 : }
408 :
409 :
410 : template <unsigned int Dim, FEFamily T>
411 1562 : void FE<Dim,T>::init_dual_shape_functions(const unsigned int n_shapes, const unsigned int n_qp)
412 : {
413 1562 : if (!this->calculate_dual)
414 0 : return;
415 :
416 428 : libmesh_assert_msg(this->calculate_phi,
417 : "dual shape function calculation relies on "
418 : "primal shape functions being calculated");
419 :
420 1562 : this->dual_phi.resize(n_shapes);
421 1562 : if (this->calculate_dphi)
422 1558 : this->dual_dphi.resize(n_shapes);
423 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
424 1562 : if (this->calculate_d2phi)
425 1477 : this->dual_d2phi.resize(n_shapes);
426 : #endif
427 :
428 26901 : for (auto i : index_range(this->dual_phi))
429 : {
430 32297 : this->dual_phi[i].resize(n_qp);
431 25339 : if (this->calculate_dphi)
432 32287 : this->dual_dphi[i].resize(n_qp);
433 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
434 25339 : if (this->calculate_d2phi)
435 31362 : this->dual_d2phi[i].resize(n_qp);
436 : #endif
437 : }
438 : }
439 :
440 : template <unsigned int Dim, FEFamily T>
441 16591293 : 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 8914604 : LOG_SCOPE("init_shape_functions()", "FE");
446 :
447 : // The number of quadrature points.
448 8913994 : const unsigned int n_qp = cast_int<unsigned int>(qp.size());
449 16591293 : this->_n_total_qp = n_qp;
450 :
451 : // Number of shape functions in the finite element approximation
452 : // space.
453 8913994 : const unsigned int n_approx_shape_functions =
454 7677299 : 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 4457302 : 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 16591293 : if (this->calculate_phi)
466 : {
467 18230540 : if (this->phi.size() == n_approx_shape_functions)
468 : {
469 13780772 : old_n_qp = n_approx_shape_functions ? this->phi[0].size() : 0;
470 3704481 : break;
471 : }
472 585891 : this->phi.resize (n_approx_shape_functions);
473 : }
474 2810521 : if (this->calculate_dphi)
475 : {
476 1843049 : if (this->dphi.size() == n_approx_shape_functions)
477 : {
478 842598 : old_n_qp = n_approx_shape_functions ? this->dphi[0].size() : 0;
479 237749 : break;
480 : }
481 604400 : this->dphi.resize (n_approx_shape_functions);
482 604400 : this->dphidx.resize (n_approx_shape_functions);
483 604400 : this->dphidy.resize (n_approx_shape_functions);
484 604400 : this->dphidz.resize (n_approx_shape_functions);
485 : }
486 :
487 1966886 : if (this->calculate_dphiref)
488 : {
489 : if (Dim > 0)
490 : {
491 1249305 : 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 981534 : this->dphidxi.resize (n_approx_shape_functions);
497 : }
498 :
499 : if (Dim > 1)
500 957721 : this->dphideta.resize (n_approx_shape_functions);
501 :
502 : if (Dim > 2)
503 784406 : this->dphidzeta.resize (n_approx_shape_functions);
504 : }
505 :
506 1955756 : if (this->calculate_curl_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
507 576 : this->curl_phi.resize(n_approx_shape_functions);
508 :
509 1955756 : 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 1955756 : if (this->calculate_d2phi)
514 : {
515 513546 : 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 410236 : this->d2phi.resize (n_approx_shape_functions);
522 410236 : this->d2phidx2.resize (n_approx_shape_functions);
523 410236 : this->d2phidxdy.resize (n_approx_shape_functions);
524 410236 : this->d2phidxdz.resize (n_approx_shape_functions);
525 410236 : this->d2phidy2.resize (n_approx_shape_functions);
526 410236 : this->d2phidydz.resize (n_approx_shape_functions);
527 410236 : this->d2phidz2.resize (n_approx_shape_functions);
528 :
529 : if (Dim > 0)
530 410236 : this->d2phidxi2.resize (n_approx_shape_functions);
531 :
532 : if (Dim > 1)
533 : {
534 388546 : this->d2phidxideta.resize (n_approx_shape_functions);
535 388546 : this->d2phideta2.resize (n_approx_shape_functions);
536 : }
537 : if (Dim > 2)
538 : {
539 382174 : this->d2phidxidzeta.resize (n_approx_shape_functions);
540 382174 : this->d2phidetadzeta.resize (n_approx_shape_functions);
541 382174 : this->d2phidzeta2.resize (n_approx_shape_functions);
542 : }
543 : }
544 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
545 : }
546 : while (false);
547 :
548 16591293 : if (old_n_qp != n_qp)
549 18453741 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
550 : {
551 16466607 : if (this->calculate_phi)
552 7798727 : this->phi[i].resize (n_qp);
553 :
554 16466607 : if (this->calculate_dphi)
555 : {
556 6510083 : this->dphi[i].resize (n_qp);
557 6510083 : this->dphidx[i].resize (n_qp);
558 6510083 : this->dphidy[i].resize (n_qp);
559 6510083 : this->dphidz[i].resize (n_qp);
560 : }
561 :
562 16465570 : if (this->calculate_dphiref)
563 : {
564 : if (Dim > 0)
565 10355403 : this->dphidxi[i].resize(n_qp);
566 :
567 : if (Dim > 1)
568 10231996 : this->dphideta[i].resize(n_qp);
569 :
570 : if (Dim > 2)
571 8632091 : this->dphidzeta[i].resize(n_qp);
572 : }
573 :
574 16466607 : if (this->calculate_curl_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
575 6252 : this->curl_phi[i].resize(n_qp);
576 :
577 16466607 : 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 16466607 : if (this->calculate_d2phi)
582 : {
583 4506007 : this->d2phi[i].resize (n_qp);
584 4506007 : this->d2phidx2[i].resize (n_qp);
585 4506007 : this->d2phidxdy[i].resize (n_qp);
586 4506007 : this->d2phidxdz[i].resize (n_qp);
587 4506007 : this->d2phidy2[i].resize (n_qp);
588 4506007 : this->d2phidydz[i].resize (n_qp);
589 4506007 : this->d2phidz2[i].resize (n_qp);
590 : if (Dim > 0)
591 4506007 : this->d2phidxi2[i].resize (n_qp);
592 : if (Dim > 1)
593 : {
594 4394840 : this->d2phidxideta[i].resize (n_qp);
595 4394840 : this->d2phideta2[i].resize (n_qp);
596 : }
597 : if (Dim > 2)
598 : {
599 4313553 : this->d2phidxidzeta[i].resize (n_qp);
600 4313553 : this->d2phidetadzeta[i].resize (n_qp);
601 4313553 : 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 12413444 : if (this->calculate_phi || this->calculate_dphi)
617 : {
618 11691707 : this->weight.resize (n_qp);
619 47424164 : for (unsigned int p=0; p<n_qp; p++)
620 48505489 : this->weight[p] = 1.;
621 : }
622 :
623 12413444 : if (this->calculate_dphi)
624 : {
625 4465260 : this->dweight.resize (n_qp);
626 4465260 : this->dphase.resize (n_qp);
627 15408772 : for (unsigned int p=0; p<n_qp; p++)
628 : {
629 10943512 : this->dweight[p].zero();
630 8066418 : 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 16522873 : if (this->calculate_dphiref && Dim > 0)
638 : {
639 13216367 : std::vector<std::vector<OutputShape>> * comps[3]
640 13216367 : { &this->dphidxi, &this->dphideta, &this->dphidzeta };
641 6340983 : 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 472036 : case 1:
657 : {
658 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
659 : // Compute the value of shape function i Hessians at quadrature point p
660 699433 : if (this->calculate_d2phi)
661 134113 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
662 492535 : for (unsigned int p=0; p<n_qp; p++)
663 385676 : this->d2phidxi2[i][p] = FE<Dim, T>::shape_second_deriv(
664 385676 : elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
665 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
666 :
667 227397 : break;
668 : }
669 :
670 :
671 :
672 : //------------------------------------------------------------
673 : // 2D
674 5997980 : case 2:
675 : {
676 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
677 : // Compute the value of shape function i Hessians at quadrature point p
678 8308950 : if (this->calculate_d2phi)
679 2922478 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
680 16531581 : for (unsigned int p=0; p<n_qp; p++)
681 : {
682 13877101 : this->d2phidxi2[i][p] = FE<Dim, T>::shape_second_deriv(
683 13877101 : elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
684 20698320 : this->d2phidxideta[i][p] = FE<Dim, T>::shape_second_deriv(
685 13877101 : elem, this->fe_type.order, i, 1, qp[p], this->_add_p_level_in_reinit);
686 13877101 : this->d2phideta2[i][p] = FE<Dim, T>::shape_second_deriv(
687 13877101 : 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 2310970 : break;
693 : }
694 :
695 :
696 :
697 : //------------------------------------------------------------
698 : // 3D
699 5614604 : case 3:
700 : {
701 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
702 : // Compute the value of shape function i Hessians at quadrature point p
703 7514490 : if (this->calculate_d2phi)
704 39482790 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
705 125913118 : for (unsigned int p=0; p<n_qp; p++)
706 : {
707 87970163 : this->d2phidxi2[i][p] = FE<Dim, T>::shape_second_deriv(
708 87970163 : elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
709 130944842 : this->d2phidxideta[i][p] = FE<Dim, T>::shape_second_deriv(
710 87970163 : elem, this->fe_type.order, i, 1, qp[p], this->_add_p_level_in_reinit);
711 130944842 : this->d2phideta2[i][p] = FE<Dim, T>::shape_second_deriv(
712 87970163 : elem, this->fe_type.order, i, 2, qp[p], this->_add_p_level_in_reinit);
713 130944842 : this->d2phidxidzeta[i][p] = FE<Dim, T>::shape_second_deriv(
714 87970163 : elem, this->fe_type.order, i, 3, qp[p], this->_add_p_level_in_reinit);
715 130944842 : this->d2phidetadzeta[i][p] = FE<Dim, T>::shape_second_deriv(
716 87970163 : elem, this->fe_type.order, i, 4, qp[p], this->_add_p_level_in_reinit);
717 87970163 : this->d2phidzeta2[i][p] = FE<Dim, T>::shape_second_deriv(
718 87970163 : 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 1899886 : break;
723 : }
724 :
725 :
726 : default:
727 : libmesh_error_msg("Invalid dimension Dim = " << Dim);
728 : }
729 :
730 16591293 : if (this->calculate_dual)
731 1319 : this->init_dual_shape_functions(n_approx_shape_functions, n_qp);
732 16591293 : }
733 :
734 : template <unsigned int Dim, FEFamily T>
735 : void
736 16117965 : 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 55410153 : for (unsigned int d=0; d != Dim; ++d)
743 : {
744 39292188 : auto & comps_d = *comps[d];
745 493756353 : for (auto i : index_range(comps_d))
746 : FE<Dim,T>::shape_derivs
747 573351979 : (elem,o,i,d,p,comps_d[i],add_p_level);
748 : }
749 16117965 : }
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 1138610465 : fdm_points(const unsigned int j, const Point & p)
796 : {
797 290683746 : libmesh_assert_less (j, LIBMESH_DIM);
798 :
799 : // cheat by using finite difference approximations:
800 290683746 : const Real eps = 1.e-6;
801 1138610465 : Point pp = p, pm = p;
802 :
803 1138610465 : switch (j)
804 : {
805 : // d()/dxi
806 134983212 : case 0:
807 : {
808 528987800 : pp(0) += eps;
809 528987800 : pm(0) -= eps;
810 528987800 : break;
811 : }
812 :
813 : // d()/deta
814 98592694 : case 1:
815 : {
816 385558730 : pp(1) += eps;
817 385558730 : pm(1) -= eps;
818 385558730 : break;
819 : }
820 :
821 : // d()/dzeta
822 57107840 : case 2:
823 : {
824 224063935 : pp(2) += eps;
825 224063935 : pm(2) -= eps;
826 224063935 : break;
827 : }
828 :
829 0 : default:
830 0 : libmesh_error_msg("Invalid derivative index j = " << j);
831 : }
832 :
833 1429294211 : return std::make_tuple(pp, pm, eps);
834 : }
835 :
836 : std::tuple<Point, Point, Real, unsigned int>
837 518975976 : 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 518975976 : Point pp = p, pm = p;
842 131324391 : unsigned int deriv_j = 0;
843 :
844 518975976 : switch (j)
845 : {
846 : // d^2() / dxi^2
847 22461599 : case 0:
848 : {
849 88771379 : pp(0) += eps;
850 88771379 : pm(0) -= eps;
851 22461599 : deriv_j = 0;
852 88771379 : break;
853 : }
854 :
855 : // d^2() / dxi deta
856 22461599 : case 1:
857 : {
858 88771379 : pp(1) += eps;
859 88771379 : pm(1) -= eps;
860 22461599 : deriv_j = 0;
861 88771379 : break;
862 : }
863 :
864 : // d^2() / deta^2
865 22461599 : case 2:
866 : {
867 88771379 : pp(1) += eps;
868 88771379 : pm(1) -= eps;
869 22461599 : deriv_j = 1;
870 88771379 : break;
871 : }
872 :
873 : // d^2()/dxidzeta
874 21313198 : case 3:
875 : {
876 84220613 : pp(2) += eps;
877 84220613 : pm(2) -= eps;
878 21313198 : deriv_j = 0;
879 84220613 : break;
880 : } // d^2()/deta^2
881 :
882 : // d^2()/detadzeta
883 21313198 : case 4:
884 : {
885 84220613 : pp(2) += eps;
886 84220613 : pm(2) -= eps;
887 21313198 : deriv_j = 1;
888 84220613 : break;
889 : }
890 :
891 : // d^2()/dzeta^2
892 21313198 : case 5:
893 : {
894 84220613 : pp(2) += eps;
895 84220613 : pm(2) -= eps;
896 21313198 : deriv_j = 2;
897 84220613 : break;
898 : }
899 :
900 0 : default:
901 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
902 : }
903 :
904 650300367 : return std::make_tuple(pp, pm, eps, deriv_j);
905 : }
906 :
907 : }
908 :
909 :
910 : template <typename OutputShape>
911 1105133921 : 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 281328192 : libmesh_assert(elem);
923 :
924 1105133921 : auto [pp, pm, eps] = fdm_points(j, p);
925 :
926 1105133921 : return (shape_func(elem, order, i, pp, add_p_level) -
927 1105133921 : 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 33476544 : 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 33476544 : auto [pp, pm, eps] = fdm_points(j, p);
961 :
962 33476544 : return (shape_func(type, order, elem, i, pp) -
963 33476544 : shape_func(type, order, elem, i, pm))/2./eps;
964 : }
965 :
966 :
967 : template <typename OutputShape>
968 : OutputShape
969 517945872 : 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 517945872 : auto [pp, pm, eps, deriv_j] = fdm_second_points(j, p);
981 :
982 517945872 : return (deriv_func(elem, order, i, deriv_j, pp, add_p_level) -
983 517945872 : 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 1030104 : 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 1030104 : auto [pp, pm, eps, deriv_j] = fdm_second_points(j, p);
1021 :
1022 1030104 : return (deriv_func(type, order, elem, i, deriv_j, pp) -
1023 1030104 : deriv_func(type, order, elem, i, deriv_j, pm))/2./eps;
1024 : }
1025 :
1026 :
1027 146433 : 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 146433 : const int extra_order = add_p_level * elem->p_level();
1034 :
1035 146433 : const int dim = elem->dim();
1036 :
1037 : const unsigned int n_sf =
1038 146433 : FEInterface::n_shape_functions(underlying_fe_type, extra_order,
1039 38905 : elem);
1040 :
1041 38905 : libmesh_assert_equal_to (n_sf, elem->n_nodes());
1042 :
1043 185338 : std::vector<Real> node_weights(n_sf);
1044 :
1045 77810 : const unsigned char datum_index = elem->mapping_data();
1046 1881100 : for (unsigned int n=0; n<n_sf; n++)
1047 2186020 : node_weights[n] =
1048 1734667 : elem->node_ref(n).get_extra_datum<Real>(datum_index);
1049 :
1050 77810 : const std::size_t n_p = p.size();
1051 :
1052 146433 : shapes.resize(n_sf);
1053 1881100 : for (unsigned int i=0; i != n_sf; ++i)
1054 : {
1055 1734667 : auto & shapes_i = shapes[i];
1056 1734667 : shapes_i.resize(n_p, 0);
1057 1734667 : FEInterface::shapes(dim, underlying_fe_type, elem, i, p,
1058 : shapes_i, add_p_level);
1059 9418774 : for (auto & s : shapes_i)
1060 9720576 : s *= node_weights[i];
1061 : }
1062 146433 : }
1063 :
1064 :
1065 189701 : 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 189701 : const int extra_order = add_p_level * elem->p_level();
1073 189701 : const unsigned int dim = elem->dim();
1074 :
1075 : const unsigned int n_sf =
1076 189701 : FEInterface::n_shape_functions(fe_type, extra_order, elem);
1077 :
1078 50947 : libmesh_assert_equal_to (n_sf, elem->n_nodes());
1079 :
1080 50947 : libmesh_assert_equal_to (dim, derivs.size());
1081 618400 : for (unsigned int d = 0; d != dim; ++d)
1082 542963 : derivs[d].resize(n_sf);
1083 :
1084 240648 : std::vector<Real> node_weights(n_sf);
1085 :
1086 101894 : const unsigned char datum_index = elem->mapping_data();
1087 2307509 : for (unsigned int n=0; n<n_sf; n++)
1088 2676828 : node_weights[n] =
1089 2117808 : elem->node_ref(n).get_extra_datum<Real>(datum_index);
1090 :
1091 101894 : const std::size_t n_p = p.size();
1092 :
1093 189701 : shapes.resize(n_sf);
1094 2307509 : for (unsigned int i=0; i != n_sf; ++i)
1095 2676828 : shapes[i].resize(n_p, 0);
1096 :
1097 189701 : FEInterface::all_shapes(dim, fe_type, elem, p, shapes, add_p_level);
1098 :
1099 2307509 : for (unsigned int i=0; i != n_sf; ++i)
1100 : {
1101 2117808 : auto & shapes_i = shapes[i];
1102 :
1103 15065975 : for (auto & s : shapes_i)
1104 16416291 : s *= node_weights[i];
1105 :
1106 7242068 : for (unsigned int d = 0; d != dim; ++d)
1107 : {
1108 6463174 : auto & derivs_di = derivs[d][i];
1109 5124260 : derivs_di.resize(n_p);
1110 5124260 : FEInterface::shape_derivs(fe_type, elem, i, d, p,
1111 : derivs_di, add_p_level);
1112 34090034 : for (auto & dip : derivs_di)
1113 36644086 : dip *= node_weights[i];
1114 : }
1115 : }
1116 189701 : }
1117 :
1118 :
1119 29952297 : 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 29952297 : int extra_order = add_p_level * elem.p_level();
1126 :
1127 : const unsigned int n_sf =
1128 29952297 : FEInterface::n_shape_functions(underlying_fe_type, extra_order, &elem);
1129 :
1130 8551638 : libmesh_assert_equal_to (n_sf, elem.n_nodes());
1131 :
1132 29952297 : std::vector<Real> node_weights(n_sf);
1133 :
1134 16976826 : const unsigned char datum_index = elem.mapping_data();
1135 :
1136 8551638 : Real weighted_shape_i = 0, weighted_sum = 0;
1137 :
1138 297152652 : for (unsigned int sf=0; sf<n_sf; sf++)
1139 : {
1140 : Real node_weight =
1141 267200355 : elem.node_ref(sf).get_extra_datum<Real>(datum_index);
1142 : Real weighted_shape = node_weight *
1143 267200355 : FEInterface::shape(underlying_fe_type, extra_order, &elem, sf, p);
1144 267200355 : weighted_sum += weighted_shape;
1145 267200355 : if (sf == i)
1146 8551638 : weighted_shape_i = weighted_shape;
1147 : }
1148 :
1149 46929123 : 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 4432345 : 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 4432345 : int extra_order = add_p_level * elem.p_level();
1244 :
1245 : const unsigned int n_sf =
1246 4432345 : FEInterface::n_shape_functions(underlying_fe_type, extra_order,
1247 1111222 : &elem);
1248 :
1249 4432345 : const unsigned int n_nodes = elem.n_nodes();
1250 1111222 : libmesh_assert_equal_to (n_sf, n_nodes);
1251 :
1252 4432345 : std::vector<Real> node_weights(n_nodes);
1253 :
1254 2222444 : const unsigned char datum_index = elem.mapping_data();
1255 102012276 : for (unsigned int n=0; n<n_nodes; n++)
1256 121986811 : node_weights[n] =
1257 97579931 : 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 102012276 : for (unsigned int sf=0; sf<n_sf; sf++)
1265 : {
1266 97579931 : Real weighted_shape = node_weights[sf] *
1267 97579931 : FEInterface::shape(underlying_fe_type, extra_order, &elem, sf,
1268 97579931 : p);
1269 97579931 : Real weighted_grada = node_weights[sf] *
1270 97579931 : FEInterface::shape_deriv(underlying_fe_type, extra_order,
1271 97579931 : &elem, sf, j1, p);
1272 97579931 : Real weighted_hess = node_weights[sf] *
1273 97579931 : FEInterface::shape_second_deriv(underlying_fe_type,
1274 97579931 : extra_order, &elem, sf, j, p);
1275 97579931 : weighted_sum += weighted_shape;
1276 97579931 : weighted_grada_sum += weighted_grada;
1277 24406880 : Real weighted_gradb = weighted_grada;
1278 97579931 : if (j1 != j2)
1279 : {
1280 48677657 : weighted_gradb = (j1 == j2) ? weighted_grada :
1281 48677657 : node_weights[sf] *
1282 48677657 : FEInterface::shape_deriv(underlying_fe_type, extra_order,
1283 : &elem, sf, j2, p);
1284 48677657 : weighted_grada_sum += weighted_grada;
1285 : }
1286 97579931 : weighted_hess_sum += weighted_hess;
1287 97579931 : 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 4432345 : if (j1 == j2)
1297 560012 : weighted_gradb_sum = weighted_grada_sum;
1298 :
1299 6654789 : return (weighted_sum * weighted_hess_i - weighted_grada_i * weighted_gradb_sum -
1300 8877233 : weighted_shape_i * weighted_hess_sum - weighted_gradb_i * weighted_grada_sum +
1301 6654789 : 2 * weighted_grada_sum * weighted_shape_i * weighted_gradb_sum / weighted_sum) /
1302 6654789 : weighted_sum / weighted_sum;
1303 : }
1304 :
1305 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1306 :
1307 :
1308 146433 : 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 116715 : std::vector<std::vector<Real>> shapes;
1315 :
1316 146433 : rational_fe_weighted_shapes(&elem, underlying_fe_type, shapes, p,
1317 : add_p_level);
1318 :
1319 224243 : std::vector<Real> shape_sums(p.size(), 0);
1320 :
1321 1881100 : for (auto i : index_range(v))
1322 : {
1323 451353 : libmesh_assert_equal_to ( p.size(), shapes[i].size() );
1324 9418774 : for (auto j : index_range(p))
1325 13793514 : shape_sums[j] += shapes[i][j];
1326 : }
1327 :
1328 1881100 : for (auto i : index_range(v))
1329 : {
1330 451353 : libmesh_assert_equal_to ( p.size(), v[i].size() );
1331 9418774 : for (auto j : index_range(v[i]))
1332 15829983 : v[i][j] = shapes[i][j] / shape_sums[j];
1333 : }
1334 146433 : }
1335 :
1336 :
1337 : template <typename OutputShape>
1338 189701 : 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 189701 : const int my_dim = elem.dim();
1345 :
1346 152841 : std::vector<std::vector<Real>> shapes;
1347 291595 : std::vector<std::vector<std::vector<Real>>> derivs(my_dim);
1348 :
1349 189701 : rational_fe_weighted_shapes_derivs(&elem, underlying_fe_type,
1350 : shapes, derivs, p, add_p_level);
1351 :
1352 291595 : std::vector<Real> shape_sums(p.size(), 0);
1353 291595 : std::vector<std::vector<Real>> shape_deriv_sums(my_dim);
1354 618400 : for (int d=0; d != my_dim; ++d)
1355 657227 : shape_deriv_sums[d].resize(p.size());
1356 :
1357 2307509 : for (auto i : index_range(shapes))
1358 : {
1359 559020 : libmesh_assert_equal_to ( p.size(), shapes[i].size() );
1360 15065975 : for (auto j : index_range(p))
1361 19884415 : shape_sums[j] += shapes[i][j];
1362 :
1363 7242068 : for (int d=0; d != my_dim; ++d)
1364 34090034 : for (auto j : index_range(p))
1365 67357334 : shape_deriv_sums[d][j] += derivs[d][i][j];
1366 : }
1367 :
1368 618400 : for (int d=0; d != my_dim; ++d)
1369 : {
1370 428699 : auto & comps_d = *comps[d];
1371 114264 : libmesh_assert_equal_to(comps_d.size(), elem.n_nodes());
1372 :
1373 5552959 : for (auto i : index_range(comps_d))
1374 : {
1375 1338914 : auto & comps_di = comps_d[i];
1376 4016742 : auto & derivs_di = derivs[d][i];
1377 :
1378 34090034 : for (auto j : index_range(comps_di))
1379 59679022 : comps_di[j] = (shape_sums[j] * derivs_di[j] -
1380 44322398 : shapes[i][j] * shape_deriv_sums[d][j]) /
1381 28965774 : shape_sums[j] / shape_sums[j];
1382 : }
1383 : }
1384 277508 : }
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
|
