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