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 31 // Attention probleme de numerotation des inconnues 32 // ------------------------------------------------- 33 // dans freefem, il y a un noeud par objets sommet, arete, element. 34 // et donc la numerotation des dl dans l'element depend 35 // de l'orientation des aretes 36 // 37 /// --------------------------------------------------------------- 38 namespace Fem2D { 39 // ------ P3 Hierarchical (just remove P1 node of the P2 finite element) -------- 40 class TypeOfFE_P3Lagrange : public TypeOfFE { 41 public: 42 static const int k = 3; 43 static const int ndf = (k + 2) * (k + 1) / 2; 44 static int Data[]; 45 static double Pi_h_coef[]; 46 static const int nn[10][3]; 47 static const int aa[10][3]; 48 static const int ff[10]; 49 static const int il[10]; 50 static const int jl[10]; 51 static const int kl[10]; 52 TypeOfFE_P3Lagrange()53 TypeOfFE_P3Lagrange( ) : TypeOfFE(3 + 2 * 3 + 1, 1, Data, 4, 1, 16, 10, 0) { 54 static const R2 Pt[10] = {R2(0 / 3., 0 / 3.), R2(3 / 3., 0 / 3.), R2(0 / 3., 3 / 3.), 55 R2(2 / 3., 1 / 3.), R2(1 / 3., 2 / 3.), R2(0 / 3., 2 / 3.), 56 R2(0 / 3., 1 / 3.), R2(1 / 3., 0 / 3.), R2(2 / 3., 0 / 3.), 57 R2(1 / 3., 1 / 3.)}; 58 // 3,4,5,6,7,8 59 int other[10] = {-1, -1, -1, 4, 3, 6, 5, 8, 7, -1}; 60 int kk = 0; 61 62 for (int i = 0; i < NbDoF; i++) { 63 pij_alpha[kk++] = IPJ(i, i, 0); 64 if (other[i] >= 0) { 65 pij_alpha[kk++] = IPJ(i, other[i], 0); 66 } 67 68 P_Pi_h[i] = Pt[i]; 69 } 70 71 assert(P_Pi_h.N( ) == NbDoF); 72 assert(pij_alpha.N( ) == kk); 73 } 74 75 void FB(const bool *whatd, const Mesh &Th, const Triangle &K, const RdHat &PHat, 76 RNMK_ &val) const; Pi_h_alpha(const baseFElement & K,KN_<double> & v) const77 void Pi_h_alpha(const baseFElement &K, KN_< double > &v) const { 78 for (int i = 0; i < 16; ++i) { 79 v[i] = 1; 80 } 81 82 int e0 = K.EdgeOrientation(0); 83 int e1 = K.EdgeOrientation(1); 84 int e2 = K.EdgeOrientation(2); 85 int ooo[6] = {e0, e0, e1, e1, e2, e2}; 86 int iii[6] = {}; 87 int jjj[6] = {}; 88 89 for (int i = 0; i < 6; ++i) { 90 iii[i] = 3 + 2 * i; // si orient = 1 91 jjj[i] = 4 + 2 * i; // si orient = -1 92 } 93 94 for (int i = 0; i < 6; ++i) { 95 if (ooo[i] == 1) { 96 v[jjj[i]] = 0; 97 } else { 98 v[iii[i]] = 0; 99 } 100 } 101 } 102 }; 103 104 // on what nu df on node node of df 105 int TypeOfFE_P3Lagrange::Data[] = { 106 0, 1, 2, 3, 3, 4, 4, 5, 5, 6, // the support number of the node of the df 107 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, // the number of the df on the node 108 0, 1, 2, 3, 3, 4, 4, 5, 5, 6, // the node of the df 109 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // the df come from which FE (generaly 0) 110 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, // which are de df on sub FE 111 0, // for each compontant $j=0,N-1$ it give the sub FE associated 112 0, 10}; 113 double TypeOfFE_P3Lagrange::Pi_h_coef[] = {1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}; FB(const bool * whatd,const Mesh &,const Triangle & K,const RdHat & PHat,RNMK_ & val) const114 void TypeOfFE_P3Lagrange::FB(const bool *whatd, const Mesh &, const Triangle &K, 115 const RdHat &PHat, RNMK_ &val) const { 116 R2 A(K[0]), B(K[1]), C(K[2]); 117 R l0 = 1 - PHat.x - PHat.y, l1 = PHat.x, l2 = PHat.y; 118 R L[3] = {l0 * k, l1 * k, l2 * k}; 119 120 throwassert(val.N( ) >= 10); 121 throwassert(val.M( ) == 1); 122 // Attention il faut renumeroter les fonction de bases 123 // car dans freefem++, il y a un node par sommet, arete or element 124 // et la numerotation naturelle mais 2 noud pas arete 125 // donc p est la perumation 126 // echange de numerotation si les arete sont dans le mauvais sens 127 int p[10] = {}; 128 129 for (int i = 0; i < 10; ++i) { 130 p[i] = i; 131 } 132 133 if (K.EdgeOrientation(0) < 0) { 134 Exchange(p[3], p[4]); // 3,4 135 } 136 137 if (K.EdgeOrientation(1) < 0) { 138 Exchange(p[5], p[6]); // 5,6 139 } 140 141 if (K.EdgeOrientation(2) < 0) { 142 Exchange(p[7], p[8]); // 7,8 143 } 144 145 val = 0; 146 147 if (whatd[op_id]) { 148 RN_ f0(val('.', 0, op_id)); 149 150 for (int df = 0; df < ndf; df++) { 151 int pdf = p[df]; 152 R f = 1. / ff[df]; 153 154 for (int i = 0; i < k; ++i) { 155 f *= L[nn[df][i]] - aa[df][i]; 156 } 157 158 f0[pdf] = f; 159 } 160 } 161 162 if (whatd[op_dx] || whatd[op_dy] || whatd[op_dxx] || whatd[op_dyy] || whatd[op_dxy]) { 163 R2 D[] = {K.H(0) * k, K.H(1) * k, K.H(2) * k}; 164 if (whatd[op_dx] || whatd[op_dy]) { 165 for (int df = 0; df < ndf; df++) { 166 int pdf = p[df]; 167 R fx = 0., fy = 0., f = 1. / ff[df]; 168 169 for (int i = 0; i < k; ++i) { 170 int n = nn[df][i]; 171 R Ln = L[n] - aa[df][i]; 172 fx = fx * Ln + f * D[n].x; 173 fy = fy * Ln + f * D[n].y; 174 f = f * Ln; 175 } 176 177 if (whatd[op_dx]) { 178 val(pdf, 0, op_dx) = fx; 179 } 180 181 if (whatd[op_dy]) { 182 val(pdf, 0, op_dy) = fy; 183 } 184 } 185 } 186 187 if (whatd[op_dyy] || whatd[op_dxy] || whatd[op_dxx]) { 188 for (int df = 0; df < ndf; df++) { 189 int pdf = p[df]; 190 R fx = 0., fy = 0., f = 1. / ff[df]; 191 R fxx = 0., fyy = 0., fxy = 0.; 192 193 for (int i = 0; i < k; ++i) { 194 int n = nn[df][i]; 195 R Ln = L[n] - aa[df][i]; 196 fxx = fxx * Ln + 2. * fx * D[n].x; 197 fyy = fyy * Ln + 2. * fy * D[n].y; 198 fxy = fxy * Ln + fx * D[n].y + fy * D[n].x; 199 fx = fx * Ln + f * D[n].x; 200 fy = fy * Ln + f * D[n].y; 201 f = f * Ln; 202 } 203 204 if (whatd[op_dxx]) { 205 val(pdf, 0, op_dxx) = fxx; 206 } 207 208 if (whatd[op_dyy]) { 209 val(pdf, 0, op_dyy) = fyy; 210 } 211 212 if (whatd[op_dxy]) { 213 val(pdf, 0, op_dxy) = fxy; 214 } 215 } 216 } 217 } 218 } 219 220 #include "Element_P3.hpp" 221 222 // Author: F. Hecht , P-H Tournier, Jet Hoe Tang jethoe.tang@googlemail.com 223 // Jan 2017 224 // in tets 225 class TypeOfFE_P3_3d : public GTypeOfFE< Mesh3 > { 226 public: 227 typedef Mesh3 Mesh; 228 typedef Mesh3::Element Element; 229 typedef GFElement< Mesh3 > FElement; 230 static const int kp = 3; // P3 231 static const int ndof = (kp + 3) * (kp + 2) * (kp + 1) / 6; 232 static int dfon[]; 233 static int nl[20][3]; 234 static int cl[20][3]; 235 static int cp[20]; 236 static int pp[20][4]; 237 static const int d = Mesh::Rd::d; 238 239 TypeOfFE_P3_3d( ); // constructor 240 void FB(const What_d whatd, const Mesh &Th, const Mesh3::Element &K, const RdHat &PHat, 241 RNMK_ &val) const; 242 void set(const Mesh &Th, const Element &K, InterpolationMatrix< RdHat > &M, int ocoef, int odf, 243 int *nump) const; 244 }; 245 246 int TypeOfFE_P3_3d::nl[20][3] = { 247 {0, 0, 0} /* 0 */, {1, 1, 1} /* 1 */, {2, 2, 2} /* 2 */, {3, 3, 3} /* 3 */, 248 {0, 0, 1} /* 4 */, {0, 1, 1} /* 5 */, {0, 0, 2} /* 6 */, {0, 2, 2} /* 7 */, 249 {0, 0, 3} /* 8 */, {0, 3, 3} /* 9 */, {1, 1, 2} /* 10 */, {1, 2, 2} /* 11 */, 250 {1, 1, 3} /* 12 */, {1, 3, 3} /* 13 */, {2, 2, 3} /* 14 */, {2, 3, 3} /* 15 */, 251 {1, 2, 3} /* 16 */, {0, 2, 3} /* 17 */, {0, 1, 3} /* 18 */, {0, 1, 2} /* 19 */}; 252 int TypeOfFE_P3_3d::cl[20][3] = { 253 {0, 1, 2} /* 0 */, {0, 1, 2} /* 1 */, {0, 1, 2} /* 2 */, {0, 1, 2} /* 3 */, 254 {0, 1, 0} /* 4 */, {0, 0, 1} /* 5 */, {0, 1, 0} /* 6 */, {0, 0, 1} /* 7 */, 255 {0, 1, 0} /* 8 */, {0, 0, 1} /* 9 */, {0, 1, 0} /* 10 */, {0, 0, 1} /* 11 */, 256 {0, 1, 0} /* 12 */, {0, 0, 1} /* 13 */, {0, 1, 0} /* 14 */, {0, 0, 1} /* 15 */, 257 {0, 0, 0} /* 16 */, {0, 0, 0} /* 17 */, {0, 0, 0} /* 18 */, {0, 0, 0} /* 19 */}; 258 int TypeOfFE_P3_3d::cp[20] = {6, 6, 6, 6, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1}; 259 int TypeOfFE_P3_3d::pp[20][4] = { 260 {3, 0, 0, 0} /* 0 */, {0, 3, 0, 0} /* 1 */, {0, 0, 3, 0} /* 2 */, {0, 0, 0, 3} /* 3 */, 261 {2, 1, 0, 0} /* 4 */, {1, 2, 0, 0} /* 5 */, {2, 0, 1, 0} /* 6 */, {1, 0, 2, 0} /* 7 */, 262 {2, 0, 0, 1} /* 8 */, {1, 0, 0, 2} /* 9 */, {0, 2, 1, 0} /* 10 */, {0, 1, 2, 0} /* 11 */, 263 {0, 2, 0, 1} /* 12 */, {0, 1, 0, 2} /* 13 */, {0, 0, 2, 1} /* 14 */, {0, 0, 1, 2} /* 15 */, 264 {0, 1, 1, 1} /* 16 */, {1, 0, 1, 1} /* 17 */, {1, 1, 0, 1} /* 18 */, {1, 1, 1, 0} /* 19 */}; 265 int TypeOfFE_P3_3d::dfon[] = {1, 2, 1, 0}; // 2 dofs on each edge, 2 dofs on each face 266 TypeOfFE_P3_3d()267 TypeOfFE_P3_3d::TypeOfFE_P3_3d( ) : GTypeOfFE< Mesh >(TypeOfFE_P3_3d::dfon, 1, 3, false, false) { 268 typedef Element E; 269 int n = this->NbDoF; 270 bool dd = verbosity > 5; 271 if (dd) { 272 cout << "\n +++ P3 : ndof : " << n << " " << this->PtInterpolation.N( ) << endl; 273 } 274 275 R3 *Pt = this->PtInterpolation; 276 // construction of interpolation ppoint 277 278 { 279 double cc = 1. / 3.; 280 281 for (int i = 0; i < ndof; ++i) { 282 Pt[i] = R3::KHat[0] * cc * pp[i][0] + R3::KHat[1] * cc * pp[i][1] + 283 R3::KHat[2] * cc * pp[i][2] + R3::KHat[3] * cc * pp[i][3]; 284 } 285 286 if (dd) { 287 cout << this->PtInterpolation << endl; 288 } 289 } 290 291 for (int i = 0; i < n; i++) { 292 this->pInterpolation[i] = i; 293 this->cInterpolation[i] = 0; 294 this->dofInterpolation[i] = i; 295 this->coefInterpolation[i] = 1.; 296 } 297 } 298 set(const Mesh & Th,const Element & K,InterpolationMatrix<RdHat> & M,int ocoef,int odf,int * nump) const299 void TypeOfFE_P3_3d::set(const Mesh &Th, const Element &K, InterpolationMatrix< RdHat > &M, 300 int ocoef, int odf, int *nump) const { 301 int n = this->NbDoF; 302 int *p = M.p; 303 304 for (int i = 0; i < n; ++i) { 305 M.p[i] = i; 306 } 307 308 if (verbosity > 9) { 309 cout << " P3 set:"; 310 } 311 312 int dof = 4; 313 314 for (int e = 0; e < 6; ++e) { 315 int oe = K.EdgeOrientation(e); 316 if (oe < 0) { 317 swap(p[dof], p[dof + 1]); 318 } 319 320 dof += 2; 321 } 322 } 323 FB(const What_d whatd,const Mesh & Th,const Mesh3::Element & K,const RdHat & PHat,RNMK_ & val) const324 void TypeOfFE_P3_3d::FB(const What_d whatd, const Mesh &Th, const Mesh3::Element &K, 325 const RdHat &PHat, RNMK_ &val) const { 326 assert(val.N( ) >= 20); // 23 degrees of freedom 327 assert(val.M( ) == 1); // 3 components 328 // int n = this->NbDoF; 329 // ------------- 330 // perm: the permutation for which the 4 tetrahedron vertices are listed with increasing GLOBAL 331 // number (i.e. perm[0] is the local number of the vertex with the smallest global number, ... 332 // perm[3] is the local number of the vertex with the biggest global number.) 333 // ------------- 334 R ld[4]; 335 PHat.toBary(ld); 336 ld[0] *= 3.; 337 ld[1] *= 3.; 338 ld[2] *= 3.; 339 ld[3] *= 3.; 340 341 int p[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19}; 342 343 { 344 int dof = 4; 345 346 for (int e = 0; e < 6; ++e) { 347 int oe = K.EdgeOrientation(e); 348 if (oe < 0) { 349 swap(p[dof], p[dof + 1]); 350 } 351 352 dof += 2; 353 } 354 } 355 static int ddd = 100; 356 ddd++; 357 val = 0.; 358 RN_ f0(val('.', 0, op_id)); 359 if (ddd < 20) { 360 cout << ld[0] << " " << ld[1] << " " << ld[2] << " " << ld[3] << " ::"; 361 } 362 363 if (whatd & Fop_D0) { 364 for (int i = 0; i < 20; ++i) { 365 R fi = 1. / cp[i]; 366 367 for (int l = 0; l < 3; ++l) { 368 fi *= ld[nl[i][l]] - cl[i][l]; 369 } 370 371 if (ddd < 20) { 372 cout << " " << fi; 373 } 374 375 f0[p[i]] = fi; 376 } 377 378 if (ddd < 20) { 379 cout << endl; 380 } 381 } 382 383 if (whatd & (Fop_D1 | Fop_D2)) { 384 R3 Dld[4], Df[20]; 385 K.Gradlambda(Dld); 386 Dld[0] *= 3.; 387 Dld[1] *= 3.; 388 Dld[2] *= 3.; 389 Dld[3] *= 3.; 390 391 for (int i = 0; i < 20; ++i) { 392 R fi = 1. / cp[i]; 393 R3 &dfi = Df[p[i]]; 394 395 for (int l = 0; l < 3; ++l) { 396 double ci = ld[nl[i][l]] - cl[i][l]; 397 dfi *= ci; 398 dfi += fi * Dld[nl[i][l]]; 399 fi *= ci; 400 } 401 402 RN_ f0x(val('.', 0, op_dx)); 403 RN_ f0y(val('.', 0, op_dy)); 404 RN_ f0z(val('.', 0, op_dz)); 405 if (whatd & Fop_dx) { 406 for (int i = 0; i < 20; ++i) { 407 f0x[i] = Df[i].x; 408 } 409 } 410 411 if (whatd & Fop_dy) { 412 for (int i = 0; i < 20; ++i) { 413 f0y[i] = Df[i].y; 414 } 415 } 416 417 if (whatd & Fop_dz) { 418 for (int i = 0; i < 20; ++i) { 419 f0z[i] = Df[i].z; 420 } 421 } 422 423 ffassert(!(whatd & Fop_D2)); // no D2 to do !!! 424 } 425 } 426 } 427 428 // link with FreeFem++ 429 static TypeOfFE_P3Lagrange P3LagrangeP3; 430 // a static variable to add the finite element to freefem++ 431 static AddNewFE P3Lagrange("P3", &P3LagrangeP3); 432 static TypeOfFE_P3_3d P3_3d; 433 GTypeOfFE< Mesh3 > &Elm_P3_3d(P3_3d); 434 435 static AddNewFE3 TFE_P3_3d("P33d", &Elm_P3_3d); init()436 static void init( ) { 437 TEF2dto3d[&P3LagrangeP3] = &Elm_P3_3d; // P3 -> P33d 438 } 439 } // namespace Fem2D 440 LOADFUNC(Fem2D::init); 441 442 // --- fin -- 443