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