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_QUAD9_HPP_
7 #define MINT_QUAD9_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 Quadratic Quadrilateral.
23  *
24  * \note The nodes are numbered in the following order: corners(4), edges(4)
25  *  and last the centroid (1) as depicted below.
26  *
27  * \verbatim
28  *
29  * quad_9:
30  *
31  *   3   6    2
32  *   +---+----+
33  *   |        |
34  * 7 +   + 8  + 5
35  *   |        |
36  *   +---+----+
37  *   0   4    1
38  *
39  * \endverbatim
40  *
41  * \see ShapeFunction
42  */
43 template <>
44 class Lagrange<mint::QUAD9> : public ShapeFunction<Lagrange<mint::QUAD9>>
45 {
46 public:
getCellType()47   static CellType getCellType() { return mint::QUAD9; }
48 
getType()49   static int getType() { return MINT_LAGRANGE_BASIS; }
50 
getNumDofs()51   static int getNumDofs() { return 9; }
52 
getMaxNewtonIters()53   static int getMaxNewtonIters() { return 16; }
54 
getDimension()55   static int getDimension() { return 2; }
56 
getMin()57   static double getMin() { return 0; }
58 
getMax()59   static double getMax() { return 1; }
60 
getCenter(double * center)61   static void getCenter(double* center)
62   {
63     SLIC_ASSERT(center != nullptr);
64 
65     center[0] = 0.5;
66     center[1] = 0.5;
67   }
68 
getCoords(double * coords)69   static void getCoords(double* coords)
70   {
71     SLIC_ASSERT(coords != nullptr);
72 
73     // corners:
74     // node 0
75     coords[0] = 0.0;
76     coords[1] = 0.0;
77 
78     // node 1
79     coords[2] = 1.0;
80     coords[3] = 0.0;
81 
82     // node 2
83     coords[4] = 1.0;
84     coords[5] = 1.0;
85 
86     // node 3
87     coords[6] = 0.0;
88     coords[7] = 1.0;
89 
90     // edges:
91     // node 4
92     coords[8] = 0.5;
93     coords[9] = 0.0;
94 
95     // node 5
96     coords[10] = 1.0;
97     coords[11] = 0.5;
98 
99     // node 6
100     coords[12] = 0.5;
101     coords[13] = 1.0;
102 
103     // node 7
104     coords[14] = 0.0;
105     coords[15] = 0.5;
106 
107     // centroid:
108     // node 8
109     coords[16] = 0.5;
110     coords[17] = 0.5;
111   }
112 
computeShape(const double * xr,double * phi)113   static void computeShape(const double* xr, double* phi)
114   {
115     SLIC_ASSERT(xr != nullptr);
116     SLIC_ASSERT(phi != nullptr);
117 
118     const double r = xr[0];
119     const double s = xr[1];
120 
121     // bar element along r
122     const double r1 = (r - 1.) * (2. * r - 1);
123     const double r2 = 4. * r * (1. - r);
124     const double r3 = r * (2. * r - 1.);
125 
126     // bar element along s
127     const double s1 = (s - 1.) * (2. * s - 1);
128     const double s2 = 4. * s * (1. - s);
129     const double s3 = s * (2. * s - 1.);
130 
131     phi[0] = r1 * s1;
132     phi[1] = r3 * s1;
133     phi[2] = r3 * s3;
134     phi[3] = r1 * s3;
135     phi[4] = r2 * s1;
136     phi[5] = r3 * s2;
137     phi[6] = r2 * s3;
138     phi[7] = r1 * s2;
139     phi[8] = r2 * s2;
140   }
141 
computeDerivatives(const double * xr,double * phidot)142   static void computeDerivatives(const double* xr, double* phidot)
143   {
144     SLIC_ASSERT(xr != nullptr);
145     SLIC_ASSERT(phidot != nullptr);
146 
147     const double r = xr[0];
148     const double s = xr[1];
149 
150     // bar element along r
151     const double r1 = (r - 1.) * (2. * r - 1);
152     const double r2 = 4. * r * (1. - r);
153     const double r3 = r * (2. * r - 1.);
154 
155     // bar element along s
156     const double s1 = (s - 1.) * (2. * s - 1);
157     const double s2 = 4. * s * (1. - s);
158     const double s3 = s * (2. * s - 1.);
159 
160     // bar element derivatives along r
161     const double dr1 = 4. * r - 3.;
162     const double dr2 = 4. - 8. * r;
163     const double dr3 = 4. * r - 1.;
164 
165     // bar element derivatives along s
166     const double ds1 = 4. * s - 3.;
167     const double ds2 = 4. - 8. * s;
168     const double ds3 = 4. * s - 1.;
169 
170     // r derivatives
171     phidot[0] = dr1 * s1;
172     phidot[1] = dr3 * s1;
173     phidot[2] = dr3 * s3;
174     phidot[3] = dr1 * s3;
175     phidot[4] = dr2 * s1;
176     phidot[5] = dr3 * s2;
177     phidot[6] = dr2 * s3;
178     phidot[7] = dr1 * s2;
179     phidot[8] = dr2 * s2;
180 
181     // s derivatives
182     phidot[9] = r1 * ds1;
183     phidot[10] = r3 * ds1;
184     phidot[11] = r3 * ds3;
185     phidot[12] = r1 * ds3;
186     phidot[13] = r2 * ds1;
187     phidot[14] = r3 * ds2;
188     phidot[15] = r2 * ds3;
189     phidot[16] = r1 * ds2;
190     phidot[17] = r2 * ds2;
191   }
192 };
193 
194 } /* namespace mint */
195 }  // namespace axom
196 
197 #endif /* MINT_QUAD_9_HPP_ */
198