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 // ------ P3dc 40 class TypeOfFE_P3dcLagrange : 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 static const R2 G; 53 static const R cshrink; 54 static const R cshrink1; Shrink(const R2 & P)55 static R2 Shrink(const R2 &P) { return (P - G) * cshrink + G; } 56 Shrink1(const R2 & P)57 static R2 Shrink1(const R2 &P) { return (P - G) * cshrink1 + G; } 58 TypeOfFE_P3dcLagrange()59 TypeOfFE_P3dcLagrange( ) : TypeOfFE(3 + 2 * 3 + 1, 1, Data, 4, 1, 10, 10, Pi_h_coef) { 60 static const R2 Pt[10] = {R2(0 / 3., 0 / 3.), R2(3 / 3., 0 / 3.), R2(0 / 3., 3 / 3.), 61 R2(2 / 3., 1 / 3.), R2(1 / 3., 2 / 3.), R2(0 / 3., 2 / 3.), 62 R2(0 / 3., 1 / 3.), R2(1 / 3., 0 / 3.), R2(2 / 3., 0 / 3.), 63 R2(1 / 3., 1 / 3.)}; 64 65 // 3,4,5,6,7,8 66 for (int i = 0; i < NbDoF; i++) { 67 pij_alpha[i] = IPJ(i, i, 0); 68 P_Pi_h[i] = Shrink(Pt[i]); 69 } 70 } 71 72 void FB(const bool *whatd, const Mesh &Th, const Triangle &K, const RdHat &P1, 73 RNMK_ &val) const; 74 }; 75 76 const R2 TypeOfFE_P3dcLagrange::G(1. / 3., 1. / 3.); 77 const R TypeOfFE_P3dcLagrange::cshrink = 1 - 1e-2; 78 const R TypeOfFE_P3dcLagrange::cshrink1 = 1. / TypeOfFE_P3dcLagrange::cshrink; 79 80 // on what nu df on node node of df 81 int TypeOfFE_P3dcLagrange::Data[] = { 82 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, // the support number of the node of the df 83 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, // the number of the df on the node 84 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // the node of the df 85 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // the df come from which FE (generaly 0) 86 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, // which are de df on sub FE 87 0, // for each compontant $j=0,N-1$ it give the sub FE associated 88 0, 10}; 89 double TypeOfFE_P3dcLagrange::Pi_h_coef[] = {1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}; FB(const bool * whatd,const Mesh &,const Triangle & K,const RdHat & P1,RNMK_ & val) const90 void TypeOfFE_P3dcLagrange::FB(const bool *whatd, const Mesh &, const Triangle &K, 91 const RdHat &P1, RNMK_ &val) const { 92 R2 P = Shrink1(P1); 93 R2 A(K[0]), B(K[1]), C(K[2]); 94 R l0 = 1 - P.x - P.y, l1 = P.x, l2 = P.y; 95 R L[3] = {l0 * k, l1 * k, l2 * k}; 96 97 throwassert(val.N( ) >= 10); 98 throwassert(val.M( ) == 1); 99 // Attention il faut renumeroter les fonction de bases 100 // car dans freefem++, il y a un node par sommet, arete or element 101 // et la numerotation naturelle mais 2 noud pas arete 102 // donc p est la perumation 103 // echange de numerotation si les arete sont dans le mauvais sens 104 int p[10]; 105 106 for (int i = 0; i < 10; ++i) { 107 p[i] = i; 108 } 109 110 /* 111 * if(K.EdgeOrientation(0) <0) Exchange(p[3],p[4]);// 3,4 112 * if(K.EdgeOrientation(1) <0) Exchange(p[5],p[6]);// 5,6 113 * if(K.EdgeOrientation(2) <0) Exchange(p[7],p[8]);// 7,8 114 */ 115 // cout << KN_<int>(p,10) <<endl; 116 val = 0; 117 /* 118 * // les fonction de base du Pk Lagrange sont 119 * // 120 * // 121 */ 122 // -- 123 124 if (whatd[op_id]) { 125 RN_ f0(val('.', 0, op_id)); 126 127 for (int df = 0; df < ndf; df++) { 128 int pdf = p[df]; 129 R f = 1. / ff[df]; 130 131 for (int i = 0; i < k; ++i) { 132 f *= L[nn[df][i]] - aa[df][i]; 133 } 134 135 f0[pdf] = f; 136 } 137 } 138 139 if (whatd[op_dx] || whatd[op_dy] || whatd[op_dxx] || whatd[op_dyy] || whatd[op_dxy]) { 140 R ks = k * cshrink1; 141 R2 D[] = {K.H(0) * ks, K.H(1) * ks, K.H(2) * ks}; 142 if (whatd[op_dx] || whatd[op_dy]) { 143 for (int df = 0; df < ndf; df++) { 144 int pdf = p[df]; 145 R fx = 0., fy = 0., f = 1. / ff[df]; 146 147 for (int i = 0; i < k; ++i) { 148 int n = nn[df][i]; 149 R Ln = L[n] - aa[df][i]; 150 fx = fx * Ln + f * D[n].x; 151 fy = fy * Ln + f * D[n].y; 152 f = f * Ln; 153 } 154 155 if (whatd[op_dx]) { 156 val(pdf, 0, op_dx) = fx; 157 } 158 159 if (whatd[op_dy]) { 160 val(pdf, 0, op_dy) = fy; 161 } 162 } 163 } 164 165 if (whatd[op_dyy] || whatd[op_dxy] || whatd[op_dxx]) { 166 for (int df = 0; df < ndf; df++) { 167 int pdf = p[df]; 168 R fx = 0., fy = 0., f = 1. / ff[df]; 169 R fxx = 0., fyy = 0., fxy = 0.; 170 171 for (int i = 0; i < k; ++i) { 172 int n = nn[df][i]; 173 R Ln = L[n] - aa[df][i]; 174 fxx = fxx * Ln + 2. * fx * D[n].x; 175 fyy = fyy * Ln + 2. * fy * D[n].y; 176 fxy = fxy * Ln + fx * D[n].y + fy * D[n].x; 177 fx = fx * Ln + f * D[n].x; 178 fy = fy * Ln + f * D[n].y; 179 f = f * Ln; 180 } 181 182 if (whatd[op_dxx]) { 183 val(pdf, 0, op_dxx) = fxx; 184 } 185 186 if (whatd[op_dyy]) { 187 val(pdf, 0, op_dyy) = fyy; 188 } 189 190 if (whatd[op_dxy]) { 191 val(pdf, 0, op_dxy) = fxy; 192 } 193 } 194 } 195 } 196 } 197 198 #include "Element_P3dc.hpp" 199 200 // link with FreeFem++ 201 static TypeOfFE_P3dcLagrange P3dcLagrangeP3dc; 202 // a static variable to add the finite element to freefem++ 203 static AddNewFE P3dcLagrange("P3dc", &P3dcLagrangeP3dc); 204 } // namespace Fem2D 205 206 // --- fin -- 207