1 /****************************************************************************/ 2 /* This file is part of FreeFEM. */ 3 /* */ 4 /* FreeFEM is free software: you can redistribute it and/or modify */ 5 /* it under the terms of the GNU Lesser General Public License as */ 6 /* published by the Free Software Foundation, either version 3 of */ 7 /* the License, or (at your option) any later version. */ 8 /* */ 9 /* FreeFEM 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 */ 12 /* GNU Lesser General Public License for more details. */ 13 /* */ 14 /* You should have received a copy of the GNU Lesser General Public License */ 15 /* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>. */ 16 /****************************************************************************/ 17 // SUMMARY : ... 18 // LICENSE : LGPLv3 19 // ORG : LJLL Universite Pierre et Marie Curie, Paris, FRANCE 20 // AUTHORS : ... 21 // E-MAIL : ... 22 23 /* clang-format off */ 24 //ff-c++-LIBRARY-dep: 25 //ff-c++-cpp-dep: 26 /* clang-format on */ 27 28 #include "ff++.hpp" 29 #include "AddNewFE.h" 30 // Attention probleme de numerotation des inconnues 31 // ------------------------------------------------- 32 // dans freefem, il y a un noeud par objets sommet, arete, element. 33 // et donc la numerotation des dl dans l'element depend 34 // de l'orientation des aretes 35 // 36 /// --------------------------------------------------------------- 37 namespace Fem2D { 38 // ------ P4 Hierarchical (just remove P1 node of the P2 finite element) -------- 39 class TypeOfFE_P4Lagrange : public TypeOfFE { 40 public: 41 static const int k = 4; 42 static const int ndf = (k + 2) * (k + 1) / 2; 43 static int Data[]; 44 static double Pi_h_coef[]; 45 static const int nn[15][4]; 46 static const int aa[15][4]; 47 static const int ff[15]; 48 static const int il[15]; 49 static const int jl[15]; 50 static const int kl[15]; 51 TypeOfFE_P4Lagrange()52 TypeOfFE_P4Lagrange( ) : TypeOfFE(3 + 3 * 3 + 3, 1, Data, 4, 1, 15 + 6, 15, 0) { 53 static const R2 Pt[15] = {R2(0 / 4., 0 / 4.), R2(4 / 4., 0 / 4.), R2(0 / 4., 4 / 4.), 54 R2(3 / 4., 1 / 4.), R2(2 / 4., 2 / 4.), R2(1 / 4., 3 / 4.), 55 R2(0 / 4., 3 / 4.), R2(0 / 4., 2 / 4.), R2(0 / 4., 1 / 4.), 56 R2(1 / 4., 0 / 4.), R2(2 / 4., 0 / 4.), R2(3 / 4., 0 / 4.), 57 R2(1 / 4., 2 / 4.), R2(2 / 4., 1 / 4.), R2(1 / 4., 1 / 4.)}; 58 59 // 3,4,5, 6,7,8, 9,10,11, 60 int other[15] = {0, 1, 2, 5, 4, 3, 8, 7, 6, 11, 10, 9, 12, 13, 14}; 61 int kk = 0; 62 63 for (int i = 0; i < NbDoF; i++) { 64 pij_alpha[kk++] = IPJ(i, i, 0); 65 if (other[i] != i) { 66 pij_alpha[kk++] = IPJ(i, other[i], 0); 67 } 68 69 P_Pi_h[i] = Pt[i]; 70 } 71 72 assert(P_Pi_h.N( ) == NbDoF); 73 assert(pij_alpha.N( ) == kk); 74 } 75 76 void FB(const bool *whatd, const Mesh &Th, const Triangle &K, const RdHat &PHat, 77 RNMK_ &val) const; Pi_h_alpha(const baseFElement & K,KN_<double> & v) const78 void Pi_h_alpha(const baseFElement &K, KN_< double > &v) const { 79 for (int i = 0; i < 15 + 6; ++i) { 80 v[i] = 1; 81 } 82 83 int e0 = K.EdgeOrientation(0); 84 int e1 = K.EdgeOrientation(1); 85 int e2 = K.EdgeOrientation(2); 86 int ooo[6] = {e0, e0, e1, e1, e2, e2}; 87 /* 3,4 88 * 5, 89 * 6,7 90 * 8,9, 91 * 10, 92 * 11,12, 93 * 13,14, 94 * 15 95 * 16,17 96 */ 97 int iii[6] = {3, 6, 8, 11, 13, 16}; 98 int jjj[6] = {}; 99 100 for (int i = 0; i < 6; ++i) { 101 jjj[i] = iii[i] + 1; // si orient = -1 102 } 103 104 for (int i = 0; i < 6; ++i) { 105 if (ooo[i] == 1) { 106 v[jjj[i]] = 0; 107 } else { 108 v[iii[i]] = 0; 109 } 110 } 111 } 112 }; 113 114 // on what nu df on node node of df 115 int TypeOfFE_P4Lagrange::Data[] = { 116 0, 1, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, // the support number of the node of the df 117 0, 0, 0, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, // the number of the df on the node 118 0, 1, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, // the node of the df 119 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // the df come from which FE (generaly 0) 120 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, // which are de df on sub FE 121 0, // for each compontant $j=0,N-1$ it give the sub FE associated 122 0, 15}; FB(const bool * whatd,const Mesh &,const Triangle & K,const RdHat & PHat,RNMK_ & val) const123 void TypeOfFE_P4Lagrange::FB(const bool *whatd, const Mesh &, const Triangle &K, 124 const RdHat &PHat, RNMK_ &val) const { 125 R2 A(K[0]), B(K[1]), C(K[2]); 126 R l0 = 1. - PHat.x - PHat.y, l1 = PHat.x, l2 = PHat.y; 127 R L[3] = {l0 * k, l1 * k, l2 * k}; 128 129 throwassert(val.N( ) >= 10); 130 throwassert(val.M( ) == 1); 131 // Attention il faut renumeroter les fonction de bases 132 // car dans freefem++, il y a un node par sommet, arete or element 133 // et la numerotation naturelle mais 2 noud pas arete 134 // donc p est la perumation 135 // echange de numerotation si les arete sont dans le mauvais sens 136 int p[15] = {}; 137 138 for (int i = 0; i < 15; ++i) { 139 p[i] = i; 140 } 141 142 if (K.EdgeOrientation(0) < 0) { 143 Exchange(p[3], p[5]); // 3,4 144 } 145 146 if (K.EdgeOrientation(1) < 0) { 147 Exchange(p[6], p[8]); // 5,6 148 } 149 150 if (K.EdgeOrientation(2) < 0) { 151 Exchange(p[9], p[11]); // 7,8 152 } 153 154 val = 0; 155 /* 156 * // les fonction de base du Pk Lagrange sont 157 */ 158 159 if (whatd[op_id]) { 160 RN_ f0(val('.', 0, op_id)); 161 162 for (int df = 0; df < ndf; df++) { 163 int pdf = p[df]; 164 R f = 1. / ff[df]; 165 166 for (int i = 0; i < k; ++i) { 167 f *= L[nn[df][i]] - aa[df][i]; 168 } 169 170 f0[pdf] = f; 171 } 172 } 173 174 if (whatd[op_dx] || whatd[op_dy] || whatd[op_dxx] || whatd[op_dyy] || whatd[op_dxy]) { 175 R2 D[] = {K.H(0) * k, K.H(1) * k, K.H(2) * k}; 176 if (whatd[op_dx] || whatd[op_dy]) { 177 for (int df = 0; df < ndf; df++) { 178 int pdf = p[df]; 179 R fx = 0., fy = 0., f = 1. / ff[df]; 180 181 for (int i = 0; i < k; ++i) { 182 int n = nn[df][i]; 183 R Ln = L[n] - aa[df][i]; 184 fx = fx * Ln + f * D[n].x; 185 fy = fy * Ln + f * D[n].y; 186 f = f * Ln; 187 } 188 189 if (whatd[op_dx]) { 190 val(pdf, 0, op_dx) = fx; 191 } 192 193 if (whatd[op_dy]) { 194 val(pdf, 0, op_dy) = fy; 195 } 196 } 197 } 198 199 if (whatd[op_dyy] || whatd[op_dxy] || whatd[op_dxx]) { 200 for (int df = 0; df < ndf; df++) { 201 int pdf = p[df]; 202 R fx = 0., fy = 0., f = 1. / ff[df]; 203 R fxx = 0., fyy = 0., fxy = 0.; 204 205 for (int i = 0; i < k; ++i) { 206 int n = nn[df][i]; 207 R Ln = L[n] - aa[df][i]; 208 fxx = fxx * Ln + 2. * fx * D[n].x; 209 fyy = fyy * Ln + 2. * fy * D[n].y; 210 fxy = fxy * Ln + fx * D[n].y + fy * D[n].x; 211 fx = fx * Ln + f * D[n].x; 212 fy = fy * Ln + f * D[n].y; 213 f = f * Ln; 214 } 215 216 if (whatd[op_dxx]) { 217 val(pdf, 0, op_dxx) = fxx; 218 } 219 220 if (whatd[op_dyy]) { 221 val(pdf, 0, op_dyy) = fyy; 222 } 223 224 if (whatd[op_dxy]) { 225 val(pdf, 0, op_dxy) = fxy; 226 } 227 } 228 } 229 } 230 } 231 232 #include "Element_P4.hpp" 233 234 // link with FreeFem++ 235 static TypeOfFE_P4Lagrange P4LagrangeP4; 236 // a static variable to add the finite element to freefem++ 237 static AddNewFE P4Lagrange("P4", &P4LagrangeP4); 238 } // namespace Fem2D 239 240 // --- fin -- 241