1 // Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
2 // other Axom Project Developers. See the top-level LICENSE file for details.
3 //
4 // SPDX-License-Identifier: (BSD-3-Clause)
5 
6 #ifndef MINT_LAGRANGE_PYRA_5_HPP_
7 #define MINT_LAGRANGE_PYRA_5_HPP_
8 
9 // Mint includes
10 #include "axom/mint/mesh/CellTypes.hpp"
11 #include "axom/mint/fem/FEBasisTypes.hpp"
12 #include "axom/mint/fem/shape_functions/ShapeFunction.hpp"
13 
14 // Slic includes
15 #include "axom/slic/interface/slic.hpp"
16 
17 namespace axom
18 {
19 namespace mint
20 {
21 /*!
22  * \brief Lagrange Finite Element definition for the Linear Pyramnid
23  *
24  * \verbatim
25  *
26  * pyra_5:
27  *
28  *              4
29  *            / |\
30  *           /  | \
31  *          /   |  \
32  *        0/_ __|_ 3\
33  *         \    |    \
34  *          \ _ |_  _ \
35  *            1       2
36  *
37  * \endverbatim
38  *
39  * \warning The jacobian for the pyramid elements may become singular near the
40  *  apex in which case the isoparametric mapping will fail.
41  *
42  * \see ShapeFunction
43  */
44 template <>
45 class Lagrange<mint::PYRAMID> : public ShapeFunction<Lagrange<mint::PYRAMID>>
46 {
47 public:
getCellType()48   static CellType getCellType() { return mint::PYRAMID; }
49 
getType()50   static int getType() { return MINT_LAGRANGE_BASIS; }
51 
getNumDofs()52   static int getNumDofs() { return 5; }
53 
getMaxNewtonIters()54   static int getMaxNewtonIters() { return 16; }
55 
getDimension()56   static int getDimension() { return 3; }
57 
getMin()58   static double getMin() { return 0; }
59 
getMax()60   static double getMax() { return 1; }
61 
getCenter(double * center)62   static void getCenter(double* center)
63   {
64     SLIC_ASSERT(center != nullptr);
65     center[0] = center[1] = 0.4;
66     center[2] = 0.2;
67   }
68 
getCoords(double * coords)69   static void getCoords(double* coords)
70   {
71     SLIC_ASSERT(coords != nullptr);
72 
73     // node 0
74     coords[0] = 0.0;
75     coords[1] = 0.0;
76     coords[2] = 0.0;
77 
78     // node 1
79     coords[3] = 1.0;
80     coords[4] = 0.0;
81     coords[5] = 0.0;
82 
83     // node 2
84     coords[6] = 1.0;
85     coords[7] = 1.0;
86     coords[8] = 0.0;
87 
88     // node 3
89     coords[9] = 0.0;
90     coords[10] = 1.0;
91     coords[11] = 0.0;
92 
93     // node 4
94     coords[12] = 0.0;
95     coords[13] = 0.0;
96     coords[14] = 1.0;
97   }
98 
computeShape(const double * xr,double * phi)99   static void computeShape(const double* xr, double* phi)
100   {
101     SLIC_ASSERT(xr != nullptr);
102     SLIC_ASSERT(phi != nullptr);
103 
104     const double r = xr[0];
105     const double s = xr[1];
106     const double t = xr[2];
107     const double rm = 1. - r;
108     const double sm = 1. - s;
109     const double tm = 1. - t;
110 
111     phi[0] = rm * sm * tm;
112     phi[1] = r * sm * tm;
113     phi[2] = r * s * tm;
114     phi[3] = rm * s * tm;
115     phi[4] = t;
116   }
117 
computeDerivatives(const double * xr,double * phidot)118   static void computeDerivatives(const double* xr, double* phidot)
119   {
120     SLIC_ASSERT(xr != nullptr);
121     SLIC_ASSERT(phidot != nullptr);
122 
123     const double r = xr[0];
124     const double s = xr[1];
125     const double t = xr[2];
126     const double rm = 1. - r;
127     const double sm = 1. - s;
128     const double tm = 1. - t;
129 
130     // r-derivatives
131     phidot[0] = -sm * tm;
132     phidot[1] = sm * tm;
133     phidot[2] = s * tm;
134     phidot[3] = -s * tm;
135     phidot[4] = 0.0;
136 
137     // s-derivatives
138     phidot[5] = -rm * tm;
139     phidot[6] = -r * tm;
140     phidot[7] = r * tm;
141     phidot[8] = rm * tm;
142     phidot[9] = 0.0;
143 
144     // t-derivatives
145     phidot[10] = -rm * sm;
146     phidot[11] = -r * sm;
147     phidot[12] = -r * s;
148     phidot[13] = -rm * s;
149     phidot[14] = 1.0;
150   }
151 };
152 
153 } /* namespace mint */
154 } /* namespace axom */
155 
156 #endif /* MINT_LAGRANGE_PYRA_5_HPP_ */
157