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 : HCT finite element 18 // LICENSE : LGPLv3 19 // ORG : LJLL Universite Pierre et Marie Curie, Paris, FRANCE 20 // AUTHORS : Frederic Hecht 21 // Hanen Narje 22 // E-MAIL : frederic.hecht@sorbonne-universite.fr 23 // ferchichihanen@gmail.com 24 25 /* clang-format off */ 26 //ff-c++-LIBRARY-dep: 27 //ff-c++-cpp-dep: 28 /* clang-format on */ 29 30 // Reference 31 // Sur l'implementation des elements finis de Hsieh-Clough-Tocher complet et reduit 32 // https://hal.inria.fr/file/index/docid/76557/filename/RR-0004.pdf 33 34 // Related files: 35 // to check and validate: testFEHCT.edp 36 // to get a real example: bilapHCT.edp 37 38 #include "ff++.hpp" 39 #include "AddNewFE.h" 40 41 namespace Fem2D { 42 class TypeOfFE_HCT : public TypeOfFE { 43 public: 44 static int Data[]; 45 TypeOfFE_HCT()46 TypeOfFE_HCT( ) 47 : TypeOfFE(12, 48 3, // hack u, u_x, u_y for interpolation 49 Data, 2, 1, 50 9 + 6, // nb coef to build interpolation 51 6, // np point to build interpolation 52 0) { 53 const R2 Pt[] = {R2(0, 0), R2(1, 0), R2(0, 1), R2(0.5, 0.5), R2(0, 0.5), R2(0.5, 0)}; 54 // for the 3 vertices 3 coef => 9 coef .. 55 int kk = 0; 56 57 for (int p = 0; p < 3; p++) { 58 P_Pi_h[p] = Pt[p]; 59 pij_alpha[kk] = IPJ(kk, p, 0); 60 kk++; // VALUE 61 pij_alpha[kk] = IPJ(kk, p, 1); 62 kk++; // DX 63 pij_alpha[kk] = IPJ(kk, p, 2); 64 kk++; // DY 65 } 66 67 int p = 3; 68 69 for (int e = 0; e < 3; ++e) { // point d'integration sur l'arete e 70 P_Pi_h[p] = Pt[p]; 71 pij_alpha[kk++] = IPJ(9 + e, p, 1); // coef = ne_x * sge 72 pij_alpha[kk++] = IPJ(9 + e, p, 2); // coef = ne_y * sge 73 p++; 74 } 75 76 assert(P_Pi_h.N( ) == p); 77 assert(pij_alpha.N( ) == kk); 78 } 79 80 void FB(const bool *whatd, const Mesh &Th, const Triangle &K, const R2 &P, RNMK_ &val) const; 81 void Pi_h_alpha(const baseFElement &K, KN_< double > &v) const; 82 }; 83 // on what nu df on node node of df 84 int TypeOfFE_HCT::Data[] = {0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 4, 5, // support on what 85 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 0, 0, // df on node 86 0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 4, 5, // th node 87 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // df of prevoius FE 88 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, // ++ which df on prevoiu 89 0, 0, 0, // ??? 90 0, 0, 0, 12, 12, 12}; Pi_h_alpha(const baseFElement & K,KN_<double> & v) const91 void TypeOfFE_HCT::Pi_h_alpha(const baseFElement &K, KN_< double > &v) const { 92 const Triangle &T(K.T); 93 int k = 0; 94 95 // coef pour les 3 sommets fois le 3 composantes 96 for (int i = 0; i < 9; i++) { 97 v[k++] = 1; 98 } 99 100 // integration sur les aretes 101 for (int i = 0; i < 3; i++) { 102 R2 N(T.Edge(i).perp( )); 103 N *= T.EdgeOrientation(i) / N.norme( ); 104 v[k++] = N.x; 105 v[k++] = N.y; 106 } 107 108 ffassert(v.N( ) == k); 109 } 110 set2zero(double * p,int k)111 void set2zero(double *p, int k) { 112 for (int i = 0; i < k; ++i) { 113 p[i] = 0.; 114 } 115 } 116 117 #define P3(a, b, c) a *b *c 118 #define P3abcx(x, a, b, c) a##x *b *c + a *b##x *c + a *b *c##x 119 #define P3abcxy(x, y, a, b, c) \ 120 a##x *b##y *c + a##y *b##x *c + a##y *b *c##x + a##x *b *c##y + a *b##x *c##y + a *b##y *c##x 121 #define P3abcxyz(x, y, z, a, b, c) \ 122 a##x *b##y *c##z + a##y *b##x *c##z + a##y *b##z *c##x + a##x *b##z *c##y + a##z *b##x *c##y + \ 123 a##z *b##y *c##x 124 125 #define P3X(a, b, c) P3abcx(x, a, b, c) 126 #define P3Y(a, b, c) P3abcx(y, a, b, c) 127 #define P3XY(a, b, c) P3abcxy(x, y, a, b, c) 128 #define P3XX(a, b, c) P3abcxy(x, x, a, b, c) 129 #define P3YY(a, b, c) P3abcxy(y, y, a, b, c) 130 #define P3XXX(a, b, c) P3abcxyz(x, x, x, a, b, c) 131 #define P3YYY(a, b, c) P3abcxyz(y, y, y, a, b, c) 132 #define P3XXY(a, b, c) P3abcxyz(x, x, y, a, b, c) 133 #define P3XYY(a, b, c) P3abcxyz(x, y, y, a, b, c) 134 135 #define LL10(P3) \ 136 { \ 137 P3(li, li, li), P3(li1, li1, li1), P3(li2, li2, li2), P3(li, li, li2), P3(li, li, li1), \ 138 P3(li1, li1, li), P3(li1, li1, li2), P3(li2, li2, li1), P3(li2, li2, li), P3(li, li1, li2) \ 139 } 140 FB(const bool * whatd,const Mesh &,const Triangle & K,const R2 & P,RNMK_ & val) const141 void TypeOfFE_HCT::FB(const bool *whatd, const Mesh &, const Triangle &K, const R2 &P, 142 RNMK_ &val) const { 143 typedef double R; 144 double area = K.area; 145 int Nop = val.K( ); 146 R2 A(K[0]), B(K[1]), C(K[2]); 147 R l[3] = {1 - P.x - P.y, P.x, P.y}; 148 R2 Dl[3] = {K.H(0), K.H(1), K.H(2)}; 149 R2 E[3] = {K.Edge(0), K.Edge(1), K.Edge(2)}; 150 R lg2[3] = {E[0].norme2( ), E[1].norme2( ), E[2].norme2( )}; 151 R lg[3] = {sqrt(lg2[0]), sqrt(lg2[1]), sqrt(lg2[2])}; 152 R eta[3] = {(lg2[2] - lg2[1]) / lg2[0], (lg2[0] - lg2[2]) / lg2[1], (lg2[1] - lg2[0]) / lg2[2]}; 153 double sgE[3] = {K.EdgeOrientation(0), K.EdgeOrientation(1), K.EdgeOrientation(2)}; 154 val = 0; 155 156 throwassert(val.N( ) >= 6); 157 throwassert(val.M( ) == 3); 158 159 int i0 = 0; 160 if (l[1] < l[i0]) { 161 i0 = 1; 162 } 163 164 if (l[2] < l[i0]) { 165 i0 = 2; 166 } 167 168 int i1 = (i0 + 1) % 3, i2 = (i0 + 2) % 3; 169 double etai = eta[i0], etai1 = eta[i1], etai2 = eta[i2]; 170 double li = l[i0], li1 = l[i1], li2 = l[i2]; 171 double lix = Dl[i0].x, li1x = Dl[i1].x, li2x = Dl[i2].x; 172 double liy = Dl[i0].y, li1y = Dl[i1].y, li2y = Dl[i2].y; 173 // i0,i1,i2, 174 int p12[12] = { 175 3 * i0, 3 * i1, 3 * i2, 3 * i0 + 2, 3 * i0 + 1, 3 * i1 + 2, 3 * i1 + 1, 176 3 * i2 + 2, 3 * i2 + 1, 9 + i0, 9 + i1, 9 + i2}; // renumerotation DL .. ff-> paper 177 178 /* 179 * double ll[10]={ P3(li,li,li), P3(li1,li1,li1), P3(li2,li2,li2), 180 * P3(li,li,li2), P3(li,li,li1), P3(li1,li1,li), 181 * P3(li1,li1,li2), P3(li2,li2,li1), P3(li2,li2,li), 182 * P3(li,li1,li2) }; 183 */ 184 // paper DOF fig 1.5.1 corresponding array Ai: 185 // ( P(q_ij), Dp(q_ij-1 - q_ij), Dp(q_ij+1 - q_ij), (j=0,1,2) (9 DOF) 186 // Dp(b_ij)(h_ij) ou bij = middle of edge ij (j=0,1,2) (3 DOF) 187 // Warning 188 // double ccc[] = { 1./(Dl[0],Ne[0]), 1./(Dl[1],Ne[1]), 1./(Dl[2],Ne[2]) }; 189 double c12 = 1. / 12.; 190 double Ai[12][10] = { 191 {(-0.5) * (etai1 - etai2), 0, 0, (1.5) * (3 + etai1), (1.5) * (3 - etai2), 0, 0, 0, 0, 0}, 192 {(0.5) * (1 - 2 * etai - etai2), 1, 0, (-1.5) * (1 - etai), (1.5) * (etai + etai2), 3, 3, 0, 193 0, 3 * (1 - etai)}, 194 {(0.5) * (1 + 2 * etai + etai1), 0, 1, (-1.5) * (etai + etai1), (-1.5) * (1 + etai), 0, 0, 3, 195 3, 3 * (1 + etai)}, 196 {(-c12) * (1 + etai1), 0, 0, (0.25) * (7 + etai1), (-0.5), 0, 0, 0, 0, 0}, 197 {(-c12) * (1 - etai2), 0, 0, (-0.5), (0.25) * (7 - etai2), 0, 0, 0, 0, 0}, 198 {(-c12) * (7 + etai2), 0, 0, (0.5), (0.25) * (5 + etai2), 1, 0, 0, 0, -1}, 199 {(1. / 6.) * (4 - etai), 0, 0, (-0.25) * (3 - etai), (-0.25) * (5 - etai), 0, 1, 0, 0, 200 (0.5) * (3 - etai)}, 201 {(1. / 6.) * (4 + etai), 0, 0, (-0.25) * (5 + etai), (-0.25) * (3 + etai), 0, 0, 1, 0, 202 (0.5) * (3 + etai)}, 203 {(-c12) * (7 - etai1), 0, 0, (0.25) * (5 - etai1), (0.5), 0, 0, 0, 1, -1}, 204 {(4. / 3.), 0, 0, -2, -2, 0, 0, 0, 0, 4}, 205 {(-2. / 3.), 0, 0, 2, 0, 0, 0, 0, 0, 0}, 206 {(-2. / 3.), 0, 0, 0, 2, 0, 0, 0, 0, 0}}; 207 const int nnzdd = 6 * 3; // nb coef.. 208 double add[] = {1., 1., 0., 0., 1., 1., 1., 0., 0., 1., 1., 1., 0., 0., 1., 1., 1., 1.}; 209 int idd[] = {0, 1, 1, 2, 2, 3, 4, 4, 5, 5, 6, 7, 7, 8, 8, 9, 10, 11}; 210 int jdd[] = {0, 1, 2, 1, 2, 3, 4, 5, 4, 5, 6, 7, 8, 7, 8, 9, 10, 11}; 211 212 for (int i = 0, kb = 0, kc = 15; i < 3; ++i, kb += 5, ++kc) { 213 int ii = i * 5 + 1; 214 int ip = (i + 1) % 3; 215 int is = (i + 2) % 3; 216 R2 Es = E[is]; 217 R2 Ep = -E[ip]; 218 double dd[4] = {Es.x, Es.y, Ep.x, Ep.y}; 219 220 add[ii] = dd[0]; 221 add[ii + 1] = dd[2]; 222 add[ii + 2] = dd[1]; 223 add[ii + 3] = dd[3]; 224 225 add[kc] = (sgE[i] * 2 * area / lg[i]); // hauteur 226 } 227 228 double AAA[12][10]; 229 set2zero(&AAA[0][0], 120); 230 double AA[12][10]; 231 set2zero(&AA[0][0], 120); 232 233 for (int jj = 0; jj < 10; ++jj) { 234 for (int i = 0; i < 12; ++i) { 235 AAA[p12[i]][jj] += Ai[i][jj]; 236 } 237 } 238 239 for (int k = 0; k < nnzdd; ++k) { 240 int i = idd[k]; 241 int j = jdd[k]; 242 double dij = add[k]; 243 244 for (int jj = 0; jj < 10; ++jj) { 245 AA[i][jj] += dij * AAA[j][jj]; 246 } 247 } 248 249 if (whatd[op_id] || whatd[op_dx] || whatd[op_dy]) { 250 251 double ll[10] = LL10(P3); 252 double llx[10] = LL10(P3X); 253 double lly[10] = LL10(P3Y); 254 RN_ f(val('.', 0, op_id)); 255 RN_ fx(val('.', 1, op_id)); 256 RN_ fy(val('.', 2, op_id)); 257 258 for (int i = 0; i < 12; ++i) { 259 for (int j = 0; j < 10; ++j) { 260 f[i] += AA[i][j] * ll[j]; 261 fx[i] += AA[i][j] * llx[j]; 262 fy[i] += AA[i][j] * lly[j]; 263 } 264 } 265 266 if (whatd[op_dx]) { 267 val('.', 0, op_dx) = fx; 268 } 269 270 if (whatd[op_dy]) { 271 val('.', 0, op_dy) = fy; 272 } 273 } 274 275 if (whatd[op_dx] || whatd[op_dxx] || whatd[op_dxy]) { 276 // cout << "dx dxx dxy"<< endl; 277 double ll[10] = LL10(P3X); 278 double llx[10] = LL10(P3XX); 279 double lly[10] = LL10(P3XY); 280 RN_ f(val('.', 0, op_dx)); 281 RN_ fx(val('.', 1, op_dx)); 282 RN_ fy(val('.', 2, op_dx)); 283 f = 0.; 284 285 for (int i = 0; i < 12; ++i) { 286 for (int j = 0; j < 10; ++j) { 287 f[i] += AA[i][j] * ll[j]; 288 fx[i] += AA[i][j] * llx[j]; 289 fy[i] += AA[i][j] * lly[j]; 290 } 291 } 292 293 if (whatd[op_dxx]) { 294 val('.', 0, op_dxx) = fx; 295 } 296 297 if (whatd[op_dxy]) { 298 val('.', 0, op_dxy) = fy; 299 } 300 } 301 302 if (whatd[op_dy] || whatd[op_dyy]) { 303 double ll[10] = LL10(P3Y); 304 double llx[10] = LL10(P3XY); 305 double lly[10] = LL10(P3YY); 306 RN_ f(val('.', 0, op_dy)); 307 RN_ fx(val('.', 1, op_dy)); 308 RN_ fy(val('.', 2, op_dy)); 309 f = 0.; 310 311 for (int i = 0; i < 12; ++i) { 312 for (int j = 0; j < 10; ++j) { 313 f[i] += AA[i][j] * ll[j]; 314 fx[i] += AA[i][j] * llx[j]; 315 fy[i] += AA[i][j] * lly[j]; 316 } 317 } 318 319 if (whatd[op_dyy]) { 320 val('.', 0, op_dyy) = fy; 321 } 322 } 323 324 if (Nop > op_dxx) { 325 val('.', 1, op_dxx) = NAN; 326 val('.', 2, op_dxx) = NAN; 327 } 328 329 if (Nop > op_dyy) { 330 val('.', 1, op_dyy) = NAN; 331 val('.', 2, op_dyy) = NAN; 332 } 333 334 if (Nop > op_dxy) { 335 val('.', 1, op_dxy) = NAN; 336 val('.', 2, op_dxy) = NAN; 337 } 338 } 339 340 // a static variable to add the finite element to freefem++ 341 static TypeOfFE_HCT Lagrange_HCT; 342 static AddNewFE FE_HCT("HCT", &Lagrange_HCT); 343 } // namespace Fem2D 344 // --- fin -- 345