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 323978286 : 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 On QUAD8 and TRI6 these functions are biorthogonal to a locally
222 : * transformed primal basis rather than to \f$ N \f$ itself: \f$ \int N_k \f$
223 : * is not positive on those faces (exactly 0 at a TRI6 vertex, -1/3 at a QUAD8
224 : * corner), so the usual construction degenerates. For those two element types
225 : * \f$ \int \Phi_j N_k \f$ is therefore not diagonal; it factors as
226 : * \f$ D = \tilde{D} T^{-1} \f$ with \f$ \tilde{D} \f$ diagonal, leaving
227 : * \f$ D^{-1} = T \tilde{D}^{-1} \f$ sparse and cheap to apply. The functions
228 : * returned here are still expressed in the primal basis, so the type and
229 : * meaning of this return value are unchanged. See Popp, Wohlmuth, Gee and
230 : * Wall, SIAM J. Sci. Comput. 34(4):B421-B446, 2012, Sec. 4.4.1 and Eq. (4.12).
231 : */
232 0 : const std::vector<std::vector<OutputShape>> & get_dual_phi() const
233 : {
234 0 : libmesh_assert(!calculations_started || calculate_dual);
235 0 : calculate_dual = true;
236 : // Dual phi computation relies on primal phi computation
237 0 : this->request_phi();
238 0 : return dual_phi;
239 : }
240 :
241 953109 : virtual void request_phi() const override
242 953109 : { get_phi(); }
243 :
244 0 : virtual void request_dual_phi() const override
245 0 : { get_dual_phi(); }
246 :
247 : /**
248 : * \returns The shape function derivatives at the quadrature
249 : * points.
250 : */
251 107648 : const std::vector<std::vector<OutputGradient>> & get_dphi() const
252 107648 : { libmesh_assert(!calculations_started || calculate_dphi);
253 276147866 : calculate_dphi = calculate_dphiref = true; return dphi; }
254 :
255 : /**
256 : * \returns The dual shape function derivatives at the quadrature points.
257 : *
258 : * \note See get_dual_phi() for the QUAD8/TRI6 caveat.
259 : */
260 0 : const std::vector<std::vector<OutputGradient>> & get_dual_dphi() const
261 0 : { libmesh_assert(!calculations_started || calculate_dphi);
262 0 : calculate_dphi = calculate_dual = calculate_dphiref = true; return dual_dphi; }
263 :
264 101652 : virtual void request_dphi() const override
265 101652 : { get_dphi(); }
266 :
267 0 : virtual void request_dual_dphi() const override
268 0 : { get_dual_dphi(); }
269 :
270 : /**
271 : * \returns The coefficients expressing the dual basis in the primal basis.
272 : *
273 : * \note See get_dual_phi() for the QUAD8/TRI6 caveat.
274 : */
275 0 : const DenseMatrix<Real> & get_dual_coeff() const
276 0 : { return dual_coeff; }
277 :
278 : /**
279 : * \returns The curl of the shape function at the quadrature
280 : * points.
281 : */
282 : virtual_for_inffe
283 525637 : const std::vector<std::vector<OutputShape>> & get_curl_phi() const
284 202175 : { libmesh_assert(!calculations_started || calculate_curl_phi);
285 2859877 : calculate_curl_phi = calculate_dphiref = true; return curl_phi; }
286 :
287 : /**
288 : * \returns The divergence of the shape function at the quadrature
289 : * points.
290 : */
291 : virtual_for_inffe
292 140745 : const std::vector<std::vector<OutputDivergence>> & get_div_phi() const
293 47027 : { libmesh_assert(!calculations_started || calculate_div_phi);
294 687657 : calculate_div_phi = calculate_dphiref = true; return div_phi; }
295 :
296 : /**
297 : * \returns The shape function x-derivative at the quadrature
298 : * points.
299 : */
300 0 : const std::vector<std::vector<OutputShape>> & get_dphidx() const
301 0 : { libmesh_assert(!calculations_started || calculate_dphi);
302 0 : calculate_dphi = calculate_dphiref = true; return dphidx; }
303 :
304 : /**
305 : * \returns The shape function y-derivative at the quadrature
306 : * points.
307 : */
308 0 : const std::vector<std::vector<OutputShape>> & get_dphidy() const
309 0 : { libmesh_assert(!calculations_started || calculate_dphi);
310 0 : calculate_dphi = calculate_dphiref = true; return dphidy; }
311 :
312 : /**
313 : * \returns The shape function z-derivative at the quadrature
314 : * points.
315 : */
316 0 : const std::vector<std::vector<OutputShape>> & get_dphidz() const
317 0 : { libmesh_assert(!calculations_started || calculate_dphi);
318 0 : calculate_dphi = calculate_dphiref = true; return dphidz; }
319 :
320 : /**
321 : * \returns The shape function xi-derivative at the quadrature
322 : * points.
323 : */
324 8149204 : const std::vector<std::vector<OutputShape>> & get_dphidxi() const
325 8149204 : { libmesh_assert(!calculations_started || calculate_dphiref);
326 89327523 : calculate_dphiref = true; return dphidxi; }
327 :
328 : /**
329 : * \returns The shape function eta-derivative at the quadrature
330 : * points.
331 : */
332 8128085 : const std::vector<std::vector<OutputShape>> & get_dphideta() const
333 8128085 : { libmesh_assert(!calculations_started || calculate_dphiref);
334 8128085 : calculate_dphiref = true; return dphideta; }
335 :
336 : /**
337 : * \returns The shape function zeta-derivative at the quadrature
338 : * points.
339 : */
340 1254814 : const std::vector<std::vector<OutputShape>> & get_dphidzeta() const
341 1254814 : { libmesh_assert(!calculations_started || calculate_dphiref);
342 1254814 : calculate_dphiref = true; return dphidzeta; }
343 :
344 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
345 :
346 : /**
347 : * \returns The shape function second derivatives at the quadrature
348 : * points.
349 : */
350 176227 : const std::vector<std::vector<OutputTensor>> & get_d2phi() const
351 176227 : { libmesh_assert(!calculations_started || calculate_d2phi);
352 2226397 : calculate_d2phi = calculate_dphiref = true; return d2phi; }
353 :
354 : /**
355 : * \returns The dual shape function second derivatives at the quadrature
356 : * points.
357 : *
358 : * \note See get_dual_phi() for the QUAD8/TRI6 caveat.
359 : */
360 0 : const std::vector<std::vector<OutputTensor>> & get_dual_d2phi() const
361 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
362 0 : calculate_d2phi = calculate_dual = calculate_dphiref = true; return dual_d2phi; }
363 :
364 : /**
365 : * \returns The shape function second derivatives at the quadrature
366 : * points.
367 : */
368 0 : const std::vector<std::vector<OutputShape>> & get_d2phidx2() const
369 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
370 0 : calculate_d2phi = calculate_dphiref = true; return d2phidx2; }
371 :
372 : /**
373 : * \returns The shape function second derivatives at the quadrature
374 : * points.
375 : */
376 0 : const std::vector<std::vector<OutputShape>> & get_d2phidxdy() const
377 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
378 0 : calculate_d2phi = calculate_dphiref = true; return d2phidxdy; }
379 :
380 : /**
381 : * \returns The shape function second derivatives at the quadrature
382 : * points.
383 : */
384 0 : const std::vector<std::vector<OutputShape>> & get_d2phidxdz() const
385 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
386 0 : calculate_d2phi = calculate_dphiref = true; return d2phidxdz; }
387 :
388 : /**
389 : * \returns The shape function second derivatives at the quadrature
390 : * points.
391 : */
392 0 : const std::vector<std::vector<OutputShape>> & get_d2phidy2() const
393 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
394 0 : calculate_d2phi = calculate_dphiref = true; return d2phidy2; }
395 :
396 : /**
397 : * \returns The shape function second derivatives at the quadrature
398 : * points.
399 : */
400 0 : const std::vector<std::vector<OutputShape>> & get_d2phidydz() const
401 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
402 0 : calculate_d2phi = calculate_dphiref = true; return d2phidydz; }
403 :
404 : /**
405 : * \returns The shape function second derivatives at the quadrature
406 : * points.
407 : */
408 0 : const std::vector<std::vector<OutputShape>> & get_d2phidz2() const
409 0 : { libmesh_assert(!calculations_started || calculate_d2phi);
410 0 : calculate_d2phi = calculate_dphiref = true; return d2phidz2; }
411 :
412 : /**
413 : * \returns The shape function second derivatives at the quadrature
414 : * points, in reference coordinates
415 : */
416 524820 : const std::vector<std::vector<OutputShape>> & get_d2phidxi2() const
417 524820 : { libmesh_assert(!calculations_started || calculate_d2phi);
418 6608235 : calculate_d2phi = calculate_dphiref = true; return d2phidxi2; }
419 :
420 : /**
421 : * \returns The shape function second derivatives at the quadrature
422 : * points, in reference coordinates
423 : */
424 517160 : const std::vector<std::vector<OutputShape>> & get_d2phidxideta() const
425 517160 : { libmesh_assert(!calculations_started || calculate_d2phi);
426 517160 : calculate_d2phi = calculate_dphiref = true; return d2phidxideta; }
427 :
428 : /**
429 : * \returns The shape function second derivatives at the quadrature
430 : * points, in reference coordinates
431 : */
432 435393 : const std::vector<std::vector<OutputShape>> & get_d2phidxidzeta() const
433 435393 : { libmesh_assert(!calculations_started || calculate_d2phi);
434 435393 : calculate_d2phi = calculate_dphiref = true; return d2phidxidzeta; }
435 :
436 : /**
437 : * \returns The shape function second derivatives at the quadrature
438 : * points, in reference coordinates
439 : */
440 517160 : const std::vector<std::vector<OutputShape>> & get_d2phideta2() const
441 517160 : { libmesh_assert(!calculations_started || calculate_d2phi);
442 517160 : calculate_d2phi = calculate_dphiref = true; return d2phideta2; }
443 :
444 : /**
445 : * \returns The shape function second derivatives at the quadrature
446 : * points, in reference coordinates
447 : */
448 435393 : const std::vector<std::vector<OutputShape>> & get_d2phidetadzeta() const
449 435393 : { libmesh_assert(!calculations_started || calculate_d2phi);
450 435393 : calculate_d2phi = calculate_dphiref = true; return d2phidetadzeta; }
451 :
452 : /**
453 : * \returns The shape function second derivatives at the quadrature
454 : * points, in reference coordinates
455 : */
456 435393 : const std::vector<std::vector<OutputShape>> & get_d2phidzeta2() const
457 435393 : { libmesh_assert(!calculations_started || calculate_d2phi);
458 435393 : calculate_d2phi = calculate_dphiref = true; return d2phidzeta2; }
459 :
460 : #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
461 :
462 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
463 :
464 : /**
465 : * \returns The global first derivative of the phase term
466 : * which is used in infinite elements, evaluated at the
467 : * quadrature points.
468 : *
469 : * In case of the general finite element class \p FE this
470 : * field is initialized to all zero, so that the variational
471 : * formulation for an infinite element produces correct element
472 : * matrices for a mesh using both finite and infinite elements.
473 : */
474 0 : const std::vector<OutputGradient> & get_dphase() const
475 0 : { return dphase; }
476 :
477 :
478 : /**
479 : * \returns The multiplicative weight at each quadrature point.
480 : * This weight is used for certain infinite element weak
481 : * formulations, so that weighted Sobolev spaces are
482 : * used for the trial function space. This renders the
483 : * variational form easily computable.
484 : *
485 : * In case of the general finite element class \p FE this
486 : * field is initialized to all ones, so that the variational
487 : * formulation for an infinite element produces correct element
488 : * matrices for a mesh using both finite and infinite elements.
489 : */
490 2 : virtual const std::vector<Real> & get_Sobolev_weight() const
491 2 : { return weight; }
492 :
493 : /**
494 : * \returns The first global derivative of the multiplicative
495 : * weight at each quadrature point. See \p get_Sobolev_weight()
496 : * for details. In case of \p FE initialized to all zero.
497 : */
498 0 : virtual const std::vector<RealGradient> & get_Sobolev_dweight() const
499 0 : { return dweight; }
500 :
501 : /**
502 : * \returns The multiplicative weight (see \p get_Sobolev_weight)
503 : * but weighted with the radial coordinate square.
504 : *
505 : * In finite elements, this gives just 1, similar to \p get_Sobolev_Weight()
506 : */
507 1340 : virtual const std::vector<Real> & get_Sobolev_weightxR_sq() const
508 1340 : { return weight; }
509 :
510 : /**
511 : * \returns The first global derivative of the multiplicative weight
512 : * (see \p dget_Sobolev_weight) but weighted with the square of the
513 : * radial coordinate.
514 : *
515 : * In finite elements, this is 0.
516 : */
517 1340 : virtual const std::vector<RealGradient> & get_Sobolev_dweightxR_sq() const
518 1340 : { return dweight; }
519 :
520 : /**
521 : * \returns The shape function \p phi (for FE) and \p phi weighted by r/decay
522 : * for InfFE.
523 : *
524 : * To compensate for the decay function applied to the Jacobian (see \p get_JxWxdecay_sq),
525 : * the wave function \p phi should be divided by this function.
526 : *
527 : * The factor r must be compensated for by the Sobolev \p weight.
528 : * (i.e. by using \p get_Sobolev_weightxR_sq())
529 : **/
530 1340 : virtual const std::vector<std::vector<OutputShape>> & get_phi_over_decayxR () const
531 1340 : { return get_phi();}
532 :
533 : /**
534 : * \returns the gradient of the shape function (see \p get_dphi()),
535 : * but in case of \p InfFE, weighted with r/decay.
536 : * See \p get_phi_over_decayxR() for details.
537 : */
538 1340 : virtual const std::vector<std::vector<OutputGradient>> & get_dphi_over_decayxR () const
539 1340 : { return get_dphi();}
540 :
541 : /**
542 : * \returns the gradient of the shape function (see \p get_dphi()),
543 : * but in case of \p InfFE, weighted with 1/decay.
544 : *
545 : * In contrast to the shape function, its gradient stays finite
546 : * when divided by the decay function.
547 : */
548 0 : virtual const std::vector<std::vector<OutputGradient>> & get_dphi_over_decay () const
549 0 : { return get_dphi();}
550 :
551 : #endif
552 :
553 : /**
554 : * Prints the value of each shape function at each quadrature point.
555 : */
556 : virtual void print_phi(std::ostream & os) const override;
557 : virtual void print_dual_phi(std::ostream & os) const override;
558 :
559 : /**
560 : * Prints the value of each shape function's derivative
561 : * at each quadrature point.
562 : */
563 : virtual void print_dphi(std::ostream & os) const override;
564 : virtual void print_dual_dphi(std::ostream & os) const override;
565 :
566 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
567 :
568 : /**
569 : * Prints the value of each shape function's second derivatives
570 : * at each quadrature point.
571 : */
572 : virtual void print_d2phi(std::ostream & os) const override;
573 : virtual void print_dual_d2phi(std::ostream & os) const override;
574 :
575 : #endif
576 :
577 :
578 : protected:
579 :
580 :
581 :
582 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
583 :
584 : /**
585 : * Initialize the data fields for the base of an
586 : * an infinite element. Implement this in the derived
587 : * class \p FE<Dim,T>.
588 : */
589 : virtual void init_base_shape_functions(const std::vector<Point> & qp,
590 : const Elem * e) = 0;
591 :
592 : #endif
593 :
594 : /**
595 : * Determine which values are to be calculated, for both the FE
596 : * itself and for the FEMap.
597 : */
598 : virtual_for_inffe
599 : void determine_calculations();
600 :
601 : /**
602 : * \returns true iff no calculations have been requested of this
603 : * FE object or of its associated FEMap
604 : */
605 27878742 : bool calculating_nothing() const
606 : {
607 29745627 : return calculate_nothing &&
608 1866885 : !this->calculate_phi && !this->calculate_dphi &&
609 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
610 1866885 : !this->calculate_d2phi &&
611 : #endif
612 29909284 : !this->calculate_curl_phi && !this->calculate_div_phi &&
613 4372758 : !this->calculate_map;
614 : }
615 :
616 : /**
617 : * After having updated the jacobian and the transformation
618 : * from local to global coordinates in \p FEMap::compute_map(),
619 : * the first derivatives of the shape functions are
620 : * transformed to global coordinates, giving \p dphi,
621 : * \p dphidx, \p dphidy, and \p dphidz. This method
622 : * should rarely be re-defined in derived classes, but
623 : * still should be usable for children. Therefore, keep
624 : * it protected.
625 : */
626 : virtual void compute_shape_functions(const Elem * elem, const std::vector<Point> & qp) override;
627 :
628 : /**
629 : * Compute the dual basis coefficients \p dual_coeff
630 : * we rely on the \p JxW (or weights) and the \p phi values,
631 : * which can come from default or customized qrule
632 : */
633 : void compute_dual_shape_coeffs(const std::vector<Real> & JxW, const std::vector<std::vector<OutputShape>> & phi);
634 :
635 : /**
636 : * Compute \p dual_phi, \p dual_dphi, \p dual_d2phi
637 : * It is only valid for this to be called after reinit has occurred with a
638 : * quadrature rule
639 : */
640 : void compute_dual_shape_functions();
641 :
642 : /**
643 : * Object that handles computing shape function values, gradients, etc
644 : * in the physical domain.
645 : */
646 : std::unique_ptr<FETransformationBase<OutputType>> _fe_trans;
647 :
648 : /**
649 : * Shape function values.
650 : */
651 : std::vector<std::vector<OutputShape>> phi;
652 : std::vector<std::vector<OutputShape>> dual_phi;
653 :
654 : /**
655 : * Shape function derivative values.
656 : */
657 : std::vector<std::vector<OutputGradient>> dphi;
658 : std::vector<std::vector<OutputGradient>> dual_dphi;
659 :
660 : /**
661 : * Coefficient matrix for the dual basis.
662 : */
663 : mutable DenseMatrix<Real> dual_coeff;
664 :
665 : /**
666 : * Shape function curl values. Only defined for vector types.
667 : */
668 : std::vector<std::vector<OutputShape>> curl_phi;
669 :
670 : /**
671 : * Shape function divergence values. Only defined for vector types.
672 : */
673 : std::vector<std::vector<OutputDivergence>> div_phi;
674 :
675 : /**
676 : * Shape function derivatives in the xi direction.
677 : */
678 : std::vector<std::vector<OutputShape>> dphidxi;
679 :
680 : /**
681 : * Shape function derivatives in the eta direction.
682 : */
683 : std::vector<std::vector<OutputShape>> dphideta;
684 :
685 : /**
686 : * Shape function derivatives in the zeta direction.
687 : */
688 : std::vector<std::vector<OutputShape>> dphidzeta;
689 :
690 : /**
691 : * Shape function derivatives in the x direction.
692 : */
693 : std::vector<std::vector<OutputShape>> dphidx;
694 :
695 : /**
696 : * Shape function derivatives in the y direction.
697 : */
698 : std::vector<std::vector<OutputShape>> dphidy;
699 :
700 : /**
701 : * Shape function derivatives in the z direction.
702 : */
703 : std::vector<std::vector<OutputShape>> dphidz;
704 :
705 :
706 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
707 :
708 : /**
709 : * Shape function second derivative values.
710 : */
711 : std::vector<std::vector<OutputTensor>> d2phi;
712 : std::vector<std::vector<OutputTensor>> dual_d2phi;
713 :
714 : /**
715 : * Shape function second derivatives in the xi direction.
716 : */
717 : std::vector<std::vector<OutputShape>> d2phidxi2;
718 :
719 : /**
720 : * Shape function second derivatives in the xi-eta direction.
721 : */
722 : std::vector<std::vector<OutputShape>> d2phidxideta;
723 :
724 : /**
725 : * Shape function second derivatives in the xi-zeta direction.
726 : */
727 : std::vector<std::vector<OutputShape>> d2phidxidzeta;
728 :
729 : /**
730 : * Shape function second derivatives in the eta direction.
731 : */
732 : std::vector<std::vector<OutputShape>> d2phideta2;
733 :
734 : /**
735 : * Shape function second derivatives in the eta-zeta direction.
736 : */
737 : std::vector<std::vector<OutputShape>> d2phidetadzeta;
738 :
739 : /**
740 : * Shape function second derivatives in the zeta direction.
741 : */
742 : std::vector<std::vector<OutputShape>> d2phidzeta2;
743 :
744 : /**
745 : * Shape function second derivatives in the x direction.
746 : */
747 : std::vector<std::vector<OutputShape>> d2phidx2;
748 :
749 : /**
750 : * Shape function second derivatives in the x-y direction.
751 : */
752 : std::vector<std::vector<OutputShape>> d2phidxdy;
753 :
754 : /**
755 : * Shape function second derivatives in the x-z direction.
756 : */
757 : std::vector<std::vector<OutputShape>> d2phidxdz;
758 :
759 : /**
760 : * Shape function second derivatives in the y direction.
761 : */
762 : std::vector<std::vector<OutputShape>> d2phidy2;
763 :
764 : /**
765 : * Shape function second derivatives in the y-z direction.
766 : */
767 : std::vector<std::vector<OutputShape>> d2phidydz;
768 :
769 : /**
770 : * Shape function second derivatives in the z direction.
771 : */
772 : std::vector<std::vector<OutputShape>> d2phidz2;
773 :
774 : #endif
775 :
776 :
777 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
778 :
779 : //--------------------------------------------------------------
780 : /* protected members for infinite elements, which are accessed
781 : * from the outside through some inline functions
782 : */
783 :
784 :
785 : /**
786 : * Used for certain infinite element families:
787 : * the first derivatives of the phase term in global coordinates,
788 : * over all quadrature points.
789 : */
790 : std::vector<OutputGradient> dphase;
791 :
792 : /**
793 : * Used for certain infinite element families:
794 : * the global derivative of the additional radial weight \f$ 1/{r^2} \f$,
795 : * over all quadrature points.
796 : */
797 : std::vector<RealGradient> dweight;
798 :
799 : /**
800 : * Used for certain infinite element families:
801 : * the additional radial weight \f$ 1/{r^2} \f$ in local coordinates,
802 : * over all quadrature points.
803 : */
804 : std::vector<Real> weight;
805 :
806 : #endif
807 :
808 : private:
809 :
810 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
811 :
812 : /**
813 : * Make all \p InfFE<Dim,T_radial,T_map> classes friends
814 : * so that they can safely used \p FE<Dim-1,T_base> through
815 : * a \p FEGenericBase * as base approximation.
816 : */
817 : template <unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
818 : friend class InfFE;
819 :
820 : #endif
821 :
822 :
823 : };
824 :
825 : // --------------------------------------------------------------------
826 : // Generic templates. We specialize for OutputType = Real, so these are
827 : // only used for OutputType = RealVectorValue
828 : template <typename OutputType>
829 0 : void FEGenericBase<OutputType>::compute_dual_shape_functions ()
830 : {
831 0 : libmesh_error_msg(
832 : "Computation of dual shape functions for vector finite element "
833 : "families is not currently implemented");
834 : }
835 :
836 : template <typename OutputType>
837 0 : void FEGenericBase<OutputType>::compute_dual_shape_coeffs(const std::vector<Real> & /*JxW*/, const std::vector<std::vector<OutputShape>> & /*phi_vals*/)
838 : {
839 0 : libmesh_error_msg(
840 : "Computation of dual shape functions for vector finite element "
841 : "families is not currently implemented");
842 : }
843 :
844 : // -----------------------------------------------------------
845 : // Forward declaration of specialization
846 : template <>
847 : void FEGenericBase<Real>::compute_dual_shape_functions();
848 :
849 : template <>
850 : void FEGenericBase<Real>::compute_dual_shape_coeffs(const std::vector<Real> & /*JxW*/, const std::vector<std::vector<OutputShape>> & /*phi_vals*/);
851 :
852 :
853 : // Typedefs for convenience and backwards compatibility
854 : typedef FEGenericBase<Real> FEBase;
855 : typedef FEGenericBase<RealGradient> FEVectorBase;
856 :
857 :
858 :
859 :
860 : // ------------------------------------------------------------
861 : // FEGenericBase class inline members
862 : template <typename OutputType>
863 : inline
864 16955391 : FEGenericBase<OutputType>::FEGenericBase(const unsigned int d,
865 : const FEType & fet) :
866 : FEAbstract(d,fet),
867 15267342 : _fe_trans( FETransformationBase<OutputType>::build(fet) ),
868 15267342 : phi(),
869 15267342 : dual_phi(),
870 15267342 : dphi(),
871 15267342 : dual_dphi(),
872 : curl_phi(),
873 : div_phi(),
874 : dphidxi(),
875 : dphideta(),
876 : dphidzeta(),
877 : dphidx(),
878 : dphidy(),
879 : dphidz()
880 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
881 : ,d2phi(),
882 : dual_d2phi(),
883 : d2phidxi2(),
884 : d2phidxideta(),
885 : d2phidxidzeta(),
886 : d2phideta2(),
887 : d2phidetadzeta(),
888 : d2phidzeta2(),
889 : d2phidx2(),
890 : d2phidxdy(),
891 : d2phidxdz(),
892 : d2phidy2(),
893 : d2phidydz(),
894 14694900 : d2phidz2()
895 : #endif
896 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
897 : ,dphase(),
898 : dweight(),
899 2260491 : weight()
900 : #endif
901 : {
902 16955391 : }
903 :
904 :
905 :
906 : template <typename OutputType>
907 : inline
908 16955391 : FEGenericBase<OutputType>::~FEGenericBase()
909 : {
910 19215277 : }
911 :
912 : } // namespace libMesh
913 :
914 : #endif // LIBMESH_FE_BASE_H
|
