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 : #ifndef LIBMESH_FE_BASE_H
21 : #define LIBMESH_FE_BASE_H
22 :
23 : // Local includes
24 : #include "libmesh/libmesh_common.h"
25 : #include "libmesh/compare_types.h"
26 : #include "libmesh/fe_abstract.h"
27 : #include "libmesh/fe_transformation_base.h"
28 : #include "libmesh/point.h"
29 : #include "libmesh/reference_counted_object.h"
30 : #include "libmesh/tensor_tools.h"
31 : #include "libmesh/type_n_tensor.h"
32 : #include "libmesh/vector_value.h"
33 : #include "libmesh/dense_matrix.h"
34 :
35 : // C++ includes
36 : #include <cstddef>
37 : #include <vector>
38 : #include <memory>
39 :
40 : namespace libMesh
41 : {
42 :
43 :
44 : // forward declarations
45 : template <typename T> class DenseMatrix;
46 : template <typename T> class DenseVector;
47 : class BoundaryInfo;
48 : class DofConstraints;
49 : class DofMap;
50 : class Elem;
51 : class MeshBase;
52 : template <typename T> class NumericVector;
53 : class QBase;
54 : template <typename T> class FETransformationBase;
55 : class FEType;
56 :
57 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
58 : class NodeConstraints;
59 : #endif
60 :
61 : #ifdef LIBMESH_ENABLE_PERIODIC
62 : class PeriodicBoundaries;
63 : class PointLocatorBase;
64 : #endif
65 :
66 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
67 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
68 : class InfFE;
69 : #endif
70 :
71 : /**
72 : * This class forms the foundation from which generic finite
73 : * elements may be derived. In the current implementation the
74 : * templated derived class \p FE offers a wide variety of commonly
75 : * used finite element concepts. Check there for details.
76 : *
77 : * Use the \p FEGenericBase<OutputType>::build() method to create an
78 : * object of any of the derived classes which is compatible with
79 : * OutputType.
80 : *
81 : * \author Benjamin S. Kirk
82 : * \date 2002
83 : */
84 : template <typename OutputType>
85 : class FEGenericBase : public FEAbstract
86 : {
87 : protected:
88 :
89 : /**
90 : * Constructor. Optionally initializes required data
91 : * structures. Protected so that this base class
92 : * cannot be explicitly instantiated.
93 : */
94 : FEGenericBase (const unsigned int dim,
95 : const FEType & fet);
96 :
97 : public:
98 :
99 : /**
100 : * Destructor.
101 : */
102 : virtual ~FEGenericBase();
103 :
104 : /**
105 : * Builds a specific finite element type. A \p
106 : * std::unique_ptr<FEGenericBase> is returned to prevent a memory leak. This
107 : * way the user need not remember to delete the object.
108 : *
109 : * The build call will fail if the OutputType of this class is not
110 : * compatible with the output required for the requested \p type
111 : */
112 : static std::unique_ptr<FEGenericBase> build (const unsigned int dim,
113 : const FEType & type);
114 :
115 : /**
116 : * Convenient typedefs for gradients of output, hessians of output,
117 : * and potentially-complex-valued versions of same.
118 : */
119 : typedef OutputType OutputShape;
120 : typedef typename TensorTools::IncrementRank<OutputShape>::type OutputGradient;
121 : typedef typename TensorTools::IncrementRank<OutputGradient>::type OutputTensor;
122 : typedef typename TensorTools::DecrementRank<OutputShape>::type OutputDivergence;
123 : typedef typename TensorTools::MakeNumber<OutputShape>::type OutputNumber;
124 : typedef typename TensorTools::IncrementRank<OutputNumber>::type OutputNumberGradient;
125 : typedef typename TensorTools::IncrementRank<OutputNumberGradient>::type OutputNumberTensor;
126 : typedef typename TensorTools::DecrementRank<OutputNumber>::type OutputNumberDivergence;
127 :
128 :
129 :
130 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
131 :
132 : /**
133 : * Builds a specific infinite element type. A \p
134 : * std::unique_ptr<FEGenericBase> is returned to prevent a memory leak. This
135 : * way the user need not remember to delete the object.
136 : *
137 : * The build call will fail if the OutputShape of this class is not
138 : * compatible with the output required for the requested \p type
139 : */
140 : static std::unique_ptr<FEGenericBase> build_InfFE (const unsigned int dim,
141 : const FEType & type);
142 :
143 : #endif
144 :
145 : #ifdef LIBMESH_ENABLE_AMR
146 :
147 : /**
148 : * Computes the constraint matrix contributions (for
149 : * non-conforming adapted meshes) corresponding to
150 : * variable number \p var_number, using generic
151 : * projections.
152 : */
153 : static void compute_proj_constraints (DofConstraints & constraints,
154 : DofMap & dof_map,
155 : const unsigned int variable_number,
156 : const Elem * elem);
157 :
158 : /**
159 : * Creates a local projection on \p coarse_elem, based on the
160 : * DoF values in \p global_vector for it's children. Computes a
161 : * vector of coefficients corresponding to dof_indices for only the
162 : * single given \p var
163 : */
164 :
165 : static void coarsened_dof_values(const NumericVector<Number> & global_vector,
166 : const DofMap & dof_map,
167 : const Elem * coarse_elem,
168 : DenseVector<Number> & coarse_dofs,
169 : const unsigned int var,
170 : const bool use_old_dof_indices = false);
171 :
172 : /**
173 : * Creates a local projection on \p coarse_elem, based on the
174 : * DoF values in \p global_vector for it's children. Computes a
175 : * vector of coefficients corresponding to all dof_indices.
176 : */
177 :
178 : static void coarsened_dof_values(const NumericVector<Number> & global_vector,
179 : const DofMap & dof_map,
180 : const Elem * coarse_elem,
181 : DenseVector<Number> & coarse_dofs,
182 : const bool use_old_dof_indices = false);
183 :
184 : #endif // #ifdef LIBMESH_ENABLE_AMR
185 :
186 : #ifdef LIBMESH_ENABLE_PERIODIC
187 :
188 : /**
189 : * Computes the constraint matrix contributions (for
190 : * meshes with periodic boundary conditions) corresponding to
191 : * variable number \p var_number, using generic projections.
192 : */
193 : static void compute_periodic_constraints (DofConstraints & constraints,
194 : DofMap & dof_map,
195 : const PeriodicBoundaries & boundaries,
196 : const MeshBase & mesh,
197 : const PointLocatorBase * point_locator,
198 : const unsigned int variable_number,
199 : const Elem * elem);
200 :
201 : #endif // LIBMESH_ENABLE_PERIODIC
202 :
203 : /**
204 : * \returns The shape function values at the quadrature points
205 : * on the element.
206 : */
207 6520 : const std::vector<std::vector<OutputShape>> & get_phi() const
208 6520 : { libmesh_assert(!calculations_started || calculate_phi);
209 324062184 : calculate_phi = true; return phi; }
210 :
211 : /**
212 : * \returns The dual (biorthogonal) shape function values at the quadrature
213 : * points on the element.
214 : *
215 : * The dual basis is constructed to satisfy
216 : * \f$ \int \Phi_j N_k \, d\gamma = \delta_{jk} \int N_k \, d\gamma \f$
217 : * against the primal basis \f$ N \f$. That biorthogonality is what makes a
218 : * mortar coupling matrix diagonal, so a Lagrange multiplier can be condensed
219 : * with a diagonal inverse.
220 : *
221 : * \note For SECOND order LAGRANGE and L2_LAGRANGE on QUAD8, TRI6 and TRI7
222 : * these functions are biorthogonal to a locally
223 : * transformed primal basis rather than to \f$ N \f$ itself: \f$ \int N_k \f$
224 : * is not positive on those faces (exactly 0 at a TRI6 vertex, -1/3 at a QUAD8
225 : * corner), so the usual construction degenerates. For those element types
226 : * \f$ \int \Phi_j N_k \f$ is therefore not diagonal; it factors as
227 : * \f$ D = \tilde{D} T^{-1} \f$ with \f$ \tilde{D} \f$ diagonal, leaving
228 : * \f$ D^{-1} = T \tilde{D}^{-1} \f$ sparse and cheap to apply. The functions
229 : * returned here are still expressed in the primal basis, so the type and
230 : * meaning of this return value are unchanged. See Popp, Wohlmuth, Gee and
231 : * Wall, SIAM J. Sci. Comput. 34(4):B421-B446, 2012, Sec. 4.4.1 and Eq. (4.12).
232 : */
233 0 : const std::vector<std::vector<OutputShape>> & get_dual_phi() const
234 : {
235 0 : libmesh_assert(!calculations_started || calculate_dual);
236 0 : calculate_dual = true;
237 : // Dual phi computation relies on primal phi computation
238 0 : this->request_phi();
239 0 : return dual_phi;
240 : }
241 :
242 982518 : virtual void request_phi() const override
243 982518 : { get_phi(); }
244 :
245 0 : virtual void request_dual_phi() const override
246 0 : { get_dual_phi(); }
247 :
248 : /**
249 : * \returns The shape function derivatives at the quadrature
250 : * points.
251 : */
252 107648 : const std::vector<std::vector<OutputGradient>> & get_dphi() const
253 107648 : { libmesh_assert(!calculations_started || calculate_dphi);
254 276154572 : calculate_dphi = calculate_dphiref = true; return dphi; }
255 :
256 : /**
257 : * \returns The dual shape function derivatives at the quadrature points.
258 : *
259 : * \note See get_dual_phi() for the QUAD8/TRI6 caveat.
260 : */
261 0 : const std::vector<std::vector<OutputGradient>> & get_dual_dphi() const
262 0 : { libmesh_assert(!calculations_started || calculate_dphi);
263 0 : calculate_dphi = calculate_dual = calculate_dphiref = true; return dual_dphi; }
264 :
265 115497 : virtual void request_dphi() const override
266 115497 : { get_dphi(); }
267 :
268 0 : virtual void request_dual_dphi() const override
269 0 : { get_dual_dphi(); }
270 :
271 : /**
272 : * \returns The coefficients expressing the dual basis in the primal basis.
273 : *
274 : * \note See get_dual_phi() for the QUAD8/TRI6 caveat.
275 : */
276 0 : const DenseMatrix<Real> & get_dual_coeff() const
277 0 : { return dual_coeff; }
278 :
279 : /**
280 : * \returns The curl of the shape function at the quadrature
281 : * points.
282 : */
283 : virtual_for_inffe
284 525637 : const std::vector<std::vector<OutputShape>> & get_curl_phi() const
285 202175 : { libmesh_assert(!calculations_started || calculate_curl_phi);
286 2859877 : calculate_curl_phi = calculate_dphiref = true; return curl_phi; }
287 :
288 : /**
289 : * \returns The divergence of the shape function at the quadrature
290 : * points.
291 : */
292 : virtual_for_inffe
293 140745 : const std::vector<std::vector<OutputDivergence>> & get_div_phi() const
294 47027 : { libmesh_assert(!calculations_started || calculate_div_phi);
295 687657 : calculate_div_phi = calculate_dphiref = true; return div_phi; }
296 :
297 : /**
298 : * \returns The shape function x-derivative at the quadrature
299 : * points.
300 : */
301 0 : const std::vector<std::vector<OutputShape>> & get_dphidx() const
302 0 : { libmesh_assert(!calculations_started || calculate_dphi);
303 0 : calculate_dphi = calculate_dphiref = true; return dphidx; }
304 :
305 : /**
306 : * \returns The shape function y-derivative at the quadrature
307 : * points.
308 : */
309 0 : const std::vector<std::vector<OutputShape>> & get_dphidy() const
310 0 : { libmesh_assert(!calculations_started || calculate_dphi);
311 0 : calculate_dphi = calculate_dphiref = true; return dphidy; }
312 :
313 : /**
314 : * \returns The shape function z-derivative at the quadrature
315 : * points.
316 : */
317 0 : const std::vector<std::vector<OutputShape>> & get_dphidz() const
318 0 : { libmesh_assert(!calculations_started || calculate_dphi);
319 0 : calculate_dphi = calculate_dphiref = true; return dphidz; }
320 :
321 : /**
322 : * \returns The shape function xi-derivative at the quadrature
323 : * points.
324 : */
325 8344158 : const std::vector<std::vector<OutputShape>> & get_dphidxi() const
326 8344158 : { libmesh_assert(!calculations_started || calculate_dphiref);
327 91903193 : calculate_dphiref = true; return dphidxi; }
328 :
329 : /**
330 : * \returns The shape function eta-derivative at the quadrature
331 : * points.
332 : */
333 8321108 : const std::vector<std::vector<OutputShape>> & get_dphideta() const
334 8321108 : { libmesh_assert(!calculations_started || calculate_dphiref);
335 8321108 : calculate_dphiref = true; return dphideta; }
336 :
337 : /**
338 : * \returns The shape function zeta-derivative at the quadrature
339 : * points.
340 : */
341 1419146 : const std::vector<std::vector<OutputShape>> & get_dphidzeta() const
342 1419146 : { libmesh_assert(!calculations_started || calculate_dphiref);
343 1419146 : calculate_dphiref = true; return dphidzeta; }
344 :
345 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
346 :
347 : /**
348 : * \returns The shape function second derivatives at the quadrature
349 : * points.
350 : */
351 176499 : const std::vector<std::vector<OutputTensor>> & get_d2phi() const
352 176499 : { libmesh_assert(!calculations_started || calculate_d2phi);
353 2229655 : calculate_d2phi = calculate_dphiref = true; return d2phi; }
354 :
355 : /**
356 : * \returns The dual shape function second derivatives at the quadrature
357 : * points.
358 : *
359 : * \note See get_dual_phi() for the QUAD8/TRI6 caveat.
360 : */
361 0 : const std::vector<std::vector<OutputTensor>> & get_dual_d2phi() const
362 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
363 0 : calculate_d2phi = calculate_dual = calculate_dphiref = true; return dual_d2phi; }
364 :
365 : /**
366 : * \returns The shape function second derivatives at the quadrature
367 : * points.
368 : */
369 0 : const std::vector<std::vector<OutputShape>> & get_d2phidx2() const
370 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
371 0 : calculate_d2phi = calculate_dphiref = true; return d2phidx2; }
372 :
373 : /**
374 : * \returns The shape function second derivatives at the quadrature
375 : * points.
376 : */
377 0 : const std::vector<std::vector<OutputShape>> & get_d2phidxdy() const
378 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
379 0 : calculate_d2phi = calculate_dphiref = true; return d2phidxdy; }
380 :
381 : /**
382 : * \returns The shape function second derivatives at the quadrature
383 : * points.
384 : */
385 0 : const std::vector<std::vector<OutputShape>> & get_d2phidxdz() const
386 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
387 0 : calculate_d2phi = calculate_dphiref = true; return d2phidxdz; }
388 :
389 : /**
390 : * \returns The shape function second derivatives at the quadrature
391 : * points.
392 : */
393 0 : const std::vector<std::vector<OutputShape>> & get_d2phidy2() const
394 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
395 0 : calculate_d2phi = calculate_dphiref = true; return d2phidy2; }
396 :
397 : /**
398 : * \returns The shape function second derivatives at the quadrature
399 : * points.
400 : */
401 0 : const std::vector<std::vector<OutputShape>> & get_d2phidydz() const
402 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
403 0 : calculate_d2phi = calculate_dphiref = true; return d2phidydz; }
404 :
405 : /**
406 : * \returns The shape function second derivatives at the quadrature
407 : * points.
408 : */
409 0 : const std::vector<std::vector<OutputShape>> & get_d2phidz2() const
410 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
411 0 : calculate_d2phi = calculate_dphiref = true; return d2phidz2; }
412 :
413 : /**
414 : * \returns The shape function second derivatives at the quadrature
415 : * points, in reference coordinates
416 : */
417 619362 : const std::vector<std::vector<OutputShape>> & get_d2phidxi2() const
418 619362 : { libmesh_assert(!calculations_started || calculate_d2phi);
419 7861374 : calculate_d2phi = calculate_dphiref = true; return d2phidxi2; }
420 :
421 : /**
422 : * \returns The shape function second derivatives at the quadrature
423 : * points, in reference coordinates
424 : */
425 610754 : const std::vector<std::vector<OutputShape>> & get_d2phidxideta() const
426 610754 : { libmesh_assert(!calculations_started || calculate_d2phi);
427 610754 : calculate_d2phi = calculate_dphiref = true; return d2phidxideta; }
428 :
429 : /**
430 : * \returns The shape function second derivatives at the quadrature
431 : * points, in reference coordinates
432 : */
433 517559 : const std::vector<std::vector<OutputShape>> & get_d2phidxidzeta() const
434 517559 : { libmesh_assert(!calculations_started || calculate_d2phi);
435 517559 : calculate_d2phi = calculate_dphiref = true; return d2phidxidzeta; }
436 :
437 : /**
438 : * \returns The shape function second derivatives at the quadrature
439 : * points, in reference coordinates
440 : */
441 610754 : const std::vector<std::vector<OutputShape>> & get_d2phideta2() const
442 610754 : { libmesh_assert(!calculations_started || calculate_d2phi);
443 610754 : calculate_d2phi = calculate_dphiref = true; return d2phideta2; }
444 :
445 : /**
446 : * \returns The shape function second derivatives at the quadrature
447 : * points, in reference coordinates
448 : */
449 517559 : const std::vector<std::vector<OutputShape>> & get_d2phidetadzeta() const
450 517559 : { libmesh_assert(!calculations_started || calculate_d2phi);
451 517559 : calculate_d2phi = calculate_dphiref = true; return d2phidetadzeta; }
452 :
453 : /**
454 : * \returns The shape function second derivatives at the quadrature
455 : * points, in reference coordinates
456 : */
457 517559 : const std::vector<std::vector<OutputShape>> & get_d2phidzeta2() const
458 517559 : { libmesh_assert(!calculations_started || calculate_d2phi);
459 517559 : calculate_d2phi = calculate_dphiref = true; return d2phidzeta2; }
460 :
461 : #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
462 :
463 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
464 :
465 : /**
466 : * \returns The global first derivative of the phase term
467 : * which is used in infinite elements, evaluated at the
468 : * quadrature points.
469 : *
470 : * In case of the general finite element class \p FE this
471 : * field is initialized to all zero, so that the variational
472 : * formulation for an infinite element produces correct element
473 : * matrices for a mesh using both finite and infinite elements.
474 : */
475 0 : const std::vector<OutputGradient> & get_dphase() const
476 0 : { return dphase; }
477 :
478 :
479 : /**
480 : * \returns The multiplicative weight at each quadrature point.
481 : * This weight is used for certain infinite element weak
482 : * formulations, so that weighted Sobolev spaces are
483 : * used for the trial function space. This renders the
484 : * variational form easily computable.
485 : *
486 : * In case of the general finite element class \p FE this
487 : * field is initialized to all ones, so that the variational
488 : * formulation for an infinite element produces correct element
489 : * matrices for a mesh using both finite and infinite elements.
490 : */
491 2 : virtual const std::vector<Real> & get_Sobolev_weight() const
492 2 : { return weight; }
493 :
494 : /**
495 : * \returns The first global derivative of the multiplicative
496 : * weight at each quadrature point. See \p get_Sobolev_weight()
497 : * for details. In case of \p FE initialized to all zero.
498 : */
499 0 : virtual const std::vector<RealGradient> & get_Sobolev_dweight() const
500 0 : { return dweight; }
501 :
502 : /**
503 : * \returns The multiplicative weight (see \p get_Sobolev_weight)
504 : * but weighted with the radial coordinate square.
505 : *
506 : * In finite elements, this gives just 1, similar to \p get_Sobolev_Weight()
507 : */
508 1340 : virtual const std::vector<Real> & get_Sobolev_weightxR_sq() const
509 1340 : { return weight; }
510 :
511 : /**
512 : * \returns The first global derivative of the multiplicative weight
513 : * (see \p dget_Sobolev_weight) but weighted with the square of the
514 : * radial coordinate.
515 : *
516 : * In finite elements, this is 0.
517 : */
518 1340 : virtual const std::vector<RealGradient> & get_Sobolev_dweightxR_sq() const
519 1340 : { return dweight; }
520 :
521 : /**
522 : * \returns The shape function \p phi (for FE) and \p phi weighted by r/decay
523 : * for InfFE.
524 : *
525 : * To compensate for the decay function applied to the Jacobian (see \p get_JxWxdecay_sq),
526 : * the wave function \p phi should be divided by this function.
527 : *
528 : * The factor r must be compensated for by the Sobolev \p weight.
529 : * (i.e. by using \p get_Sobolev_weightxR_sq())
530 : **/
531 1340 : virtual const std::vector<std::vector<OutputShape>> & get_phi_over_decayxR () const
532 1340 : { return get_phi();}
533 :
534 : /**
535 : * \returns the gradient of the shape function (see \p get_dphi()),
536 : * but in case of \p InfFE, weighted with r/decay.
537 : * See \p get_phi_over_decayxR() for details.
538 : */
539 1340 : virtual const std::vector<std::vector<OutputGradient>> & get_dphi_over_decayxR () const
540 1340 : { return get_dphi();}
541 :
542 : /**
543 : * \returns the gradient of the shape function (see \p get_dphi()),
544 : * but in case of \p InfFE, weighted with 1/decay.
545 : *
546 : * In contrast to the shape function, its gradient stays finite
547 : * when divided by the decay function.
548 : */
549 0 : virtual const std::vector<std::vector<OutputGradient>> & get_dphi_over_decay () const
550 0 : { return get_dphi();}
551 :
552 : #endif
553 :
554 : /**
555 : * Prints the value of each shape function at each quadrature point.
556 : */
557 : virtual void print_phi(std::ostream & os) const override;
558 : virtual void print_dual_phi(std::ostream & os) const override;
559 :
560 : /**
561 : * Prints the value of each shape function's derivative
562 : * at each quadrature point.
563 : */
564 : virtual void print_dphi(std::ostream & os) const override;
565 : virtual void print_dual_dphi(std::ostream & os) const override;
566 :
567 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
568 :
569 : /**
570 : * Prints the value of each shape function's second derivatives
571 : * at each quadrature point.
572 : */
573 : virtual void print_d2phi(std::ostream & os) const override;
574 : virtual void print_dual_d2phi(std::ostream & os) const override;
575 :
576 : #endif
577 :
578 :
579 : protected:
580 :
581 :
582 :
583 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
584 :
585 : /**
586 : * Initialize the data fields for the base of an
587 : * an infinite element. Implement this in the derived
588 : * class \p FE<Dim,T>.
589 : */
590 : virtual void init_base_shape_functions(const std::vector<Point> & qp,
591 : const Elem * e) = 0;
592 :
593 : #endif
594 :
595 : /**
596 : * Determine which values are to be calculated, for both the FE
597 : * itself and for the FEMap.
598 : */
599 : virtual_for_inffe
600 : void determine_calculations();
601 :
602 : /**
603 : * \returns true iff no calculations have been requested of this
604 : * FE object or of its associated FEMap
605 : */
606 27918101 : bool calculating_nothing() const
607 : {
608 29809658 : return calculate_nothing &&
609 1891557 : !this->calculate_phi && !this->calculate_dphi &&
610 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
611 1891557 : !this->calculate_d2phi &&
612 : #endif
613 29975371 : !this->calculate_curl_phi && !this->calculate_div_phi &&
614 4400698 : !this->calculate_map;
615 : }
616 :
617 : /**
618 : * After having updated the jacobian and the transformation
619 : * from local to global coordinates in \p FEMap::compute_map(),
620 : * the first derivatives of the shape functions are
621 : * transformed to global coordinates, giving \p dphi,
622 : * \p dphidx, \p dphidy, and \p dphidz. This method
623 : * should rarely be re-defined in derived classes, but
624 : * still should be usable for children. Therefore, keep
625 : * it protected.
626 : */
627 : virtual void compute_shape_functions(const Elem * elem, const std::vector<Point> & qp) override;
628 :
629 : /**
630 : * Compute the dual basis coefficients \p dual_coeff
631 : * we rely on the \p JxW (or weights) and the \p phi values,
632 : * which can come from default or customized qrule
633 : */
634 : void compute_dual_shape_coeffs(const std::vector<Real> & JxW, const std::vector<std::vector<OutputShape>> & phi);
635 :
636 : /**
637 : * Compute \p dual_phi, \p dual_dphi, \p dual_d2phi
638 : * It is only valid for this to be called after reinit has occurred with a
639 : * quadrature rule
640 : */
641 : void compute_dual_shape_functions();
642 :
643 : /**
644 : * Object that handles computing shape function values, gradients, etc
645 : * in the physical domain.
646 : */
647 : std::unique_ptr<FETransformationBase<OutputType>> _fe_trans;
648 :
649 : /**
650 : * Shape function values.
651 : */
652 : std::vector<std::vector<OutputShape>> phi;
653 : std::vector<std::vector<OutputShape>> dual_phi;
654 :
655 : /**
656 : * Shape function derivative values.
657 : */
658 : std::vector<std::vector<OutputGradient>> dphi;
659 : std::vector<std::vector<OutputGradient>> dual_dphi;
660 :
661 : /**
662 : * Coefficient matrix for the dual basis.
663 : */
664 : mutable DenseMatrix<Real> dual_coeff;
665 :
666 : /**
667 : * Shape function curl values. Only defined for vector types.
668 : */
669 : std::vector<std::vector<OutputShape>> curl_phi;
670 :
671 : /**
672 : * Shape function divergence values. Only defined for vector types.
673 : */
674 : std::vector<std::vector<OutputDivergence>> div_phi;
675 :
676 : /**
677 : * Shape function derivatives in the xi direction.
678 : */
679 : std::vector<std::vector<OutputShape>> dphidxi;
680 :
681 : /**
682 : * Shape function derivatives in the eta direction.
683 : */
684 : std::vector<std::vector<OutputShape>> dphideta;
685 :
686 : /**
687 : * Shape function derivatives in the zeta direction.
688 : */
689 : std::vector<std::vector<OutputShape>> dphidzeta;
690 :
691 : /**
692 : * Shape function derivatives in the x direction.
693 : */
694 : std::vector<std::vector<OutputShape>> dphidx;
695 :
696 : /**
697 : * Shape function derivatives in the y direction.
698 : */
699 : std::vector<std::vector<OutputShape>> dphidy;
700 :
701 : /**
702 : * Shape function derivatives in the z direction.
703 : */
704 : std::vector<std::vector<OutputShape>> dphidz;
705 :
706 :
707 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
708 :
709 : /**
710 : * Shape function second derivative values.
711 : */
712 : std::vector<std::vector<OutputTensor>> d2phi;
713 : std::vector<std::vector<OutputTensor>> dual_d2phi;
714 :
715 : /**
716 : * Shape function second derivatives in the xi direction.
717 : */
718 : std::vector<std::vector<OutputShape>> d2phidxi2;
719 :
720 : /**
721 : * Shape function second derivatives in the xi-eta direction.
722 : */
723 : std::vector<std::vector<OutputShape>> d2phidxideta;
724 :
725 : /**
726 : * Shape function second derivatives in the xi-zeta direction.
727 : */
728 : std::vector<std::vector<OutputShape>> d2phidxidzeta;
729 :
730 : /**
731 : * Shape function second derivatives in the eta direction.
732 : */
733 : std::vector<std::vector<OutputShape>> d2phideta2;
734 :
735 : /**
736 : * Shape function second derivatives in the eta-zeta direction.
737 : */
738 : std::vector<std::vector<OutputShape>> d2phidetadzeta;
739 :
740 : /**
741 : * Shape function second derivatives in the zeta direction.
742 : */
743 : std::vector<std::vector<OutputShape>> d2phidzeta2;
744 :
745 : /**
746 : * Shape function second derivatives in the x direction.
747 : */
748 : std::vector<std::vector<OutputShape>> d2phidx2;
749 :
750 : /**
751 : * Shape function second derivatives in the x-y direction.
752 : */
753 : std::vector<std::vector<OutputShape>> d2phidxdy;
754 :
755 : /**
756 : * Shape function second derivatives in the x-z direction.
757 : */
758 : std::vector<std::vector<OutputShape>> d2phidxdz;
759 :
760 : /**
761 : * Shape function second derivatives in the y direction.
762 : */
763 : std::vector<std::vector<OutputShape>> d2phidy2;
764 :
765 : /**
766 : * Shape function second derivatives in the y-z direction.
767 : */
768 : std::vector<std::vector<OutputShape>> d2phidydz;
769 :
770 : /**
771 : * Shape function second derivatives in the z direction.
772 : */
773 : std::vector<std::vector<OutputShape>> d2phidz2;
774 :
775 : #endif
776 :
777 :
778 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
779 :
780 : //--------------------------------------------------------------
781 : /* protected members for infinite elements, which are accessed
782 : * from the outside through some inline functions
783 : */
784 :
785 :
786 : /**
787 : * Used for certain infinite element families:
788 : * the first derivatives of the phase term in global coordinates,
789 : * over all quadrature points.
790 : */
791 : std::vector<OutputGradient> dphase;
792 :
793 : /**
794 : * Used for certain infinite element families:
795 : * the global derivative of the additional radial weight \f$ 1/{r^2} \f$,
796 : * over all quadrature points.
797 : */
798 : std::vector<RealGradient> dweight;
799 :
800 : /**
801 : * Used for certain infinite element families:
802 : * the additional radial weight \f$ 1/{r^2} \f$ in local coordinates,
803 : * over all quadrature points.
804 : */
805 : std::vector<Real> weight;
806 :
807 : #endif
808 :
809 : private:
810 :
811 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
812 :
813 : /**
814 : * Make all \p InfFE<Dim,T_radial,T_map> classes friends
815 : * so that they can safely used \p FE<Dim-1,T_base> through
816 : * a \p FEGenericBase * as base approximation.
817 : */
818 : template <unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
819 : friend class InfFE;
820 :
821 : #endif
822 :
823 :
824 : };
825 :
826 : // --------------------------------------------------------------------
827 : // Generic templates. We specialize for OutputType = Real, so these are
828 : // only used for OutputType = RealVectorValue
829 : template <typename OutputType>
830 0 : void FEGenericBase<OutputType>::compute_dual_shape_functions ()
831 : {
832 0 : libmesh_error_msg(
833 : "Computation of dual shape functions for vector finite element "
834 : "families is not currently implemented");
835 : }
836 :
837 : template <typename OutputType>
838 0 : void FEGenericBase<OutputType>::compute_dual_shape_coeffs(const std::vector<Real> & /*JxW*/, const std::vector<std::vector<OutputShape>> & /*phi_vals*/)
839 : {
840 0 : libmesh_error_msg(
841 : "Computation of dual shape functions for vector finite element "
842 : "families is not currently implemented");
843 : }
844 :
845 : // -----------------------------------------------------------
846 : // Forward declaration of specialization
847 : template <>
848 : void FEGenericBase<Real>::compute_dual_shape_functions();
849 :
850 : template <>
851 : void FEGenericBase<Real>::compute_dual_shape_coeffs(const std::vector<Real> & /*JxW*/, const std::vector<std::vector<OutputShape>> & /*phi_vals*/);
852 :
853 :
854 : // Typedefs for convenience and backwards compatibility
855 : typedef FEGenericBase<Real> FEBase;
856 : typedef FEGenericBase<RealGradient> FEVectorBase;
857 :
858 :
859 :
860 :
861 : // ------------------------------------------------------------
862 : // FEGenericBase class inline members
863 : template <typename OutputType>
864 : inline
865 12861431 : FEGenericBase<OutputType>::FEGenericBase(const unsigned int d,
866 : const FEType & fet) :
867 : FEAbstract(d,fet),
868 11418504 : _fe_trans( FETransformationBase<OutputType>::build(fet) ),
869 11418504 : phi(),
870 11418504 : dual_phi(),
871 11418504 : dphi(),
872 11418504 : dual_dphi(),
873 : curl_phi(),
874 : div_phi(),
875 : dphidxi(),
876 : dphideta(),
877 : dphidzeta(),
878 : dphidx(),
879 : dphidy(),
880 : dphidz()
881 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
882 : ,d2phi(),
883 : dual_d2phi(),
884 : d2phidxi2(),
885 : d2phidxideta(),
886 : d2phidxidzeta(),
887 : d2phideta2(),
888 : d2phidetadzeta(),
889 : d2phidzeta2(),
890 : d2phidx2(),
891 : d2phidxdy(),
892 : d2phidxdz(),
893 : d2phidy2(),
894 : d2phidydz(),
895 10888189 : d2phidz2()
896 : #endif
897 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
898 : ,dphase(),
899 : dweight(),
900 1973242 : weight()
901 : #endif
902 : {
903 12861431 : }
904 :
905 :
906 :
907 : template <typename OutputType>
908 : inline
909 12861431 : FEGenericBase<OutputType>::~FEGenericBase()
910 : {
911 14828560 : }
912 :
913 : } // namespace libMesh
914 :
915 : #endif // LIBMESH_FE_BASE_H
|
