NGSolve
4.9
|
00001 namespace ngfem 00002 { 00003 00004 class TetShapesInnerLegendre 00005 { 00006 public: 00007 template <typename Sx, typename Sy, typename Sz, typename T> 00008 static int Calc (int n, Sx x, Sy y, Sz z, T & values) 00009 { 00010 int ii = 0, i, j, k; 00011 ArrayMem<Sx, 20> polx(n+1), poly(n+1), polz(n+1); 00012 00013 ScaledLegendrePolynomial (n-4, x, 1-y-z, polx); 00014 ScaledLegendrePolynomial (n-4, 2*y - (1-z) , (1-z),poly); //SZ (1-z)-1, poly); 00015 LegendrePolynomial (n-4, 2*z-1, polz); 00016 Sx bub = (1-x-y-z) * (1+x-y-z) * y * z ; 00017 00018 for (i = 0; i <= n-4; i++) 00019 for (j = 0; j <= n-4-i; j++) 00020 for (k = 0; k <= n-4-i-j; k++) 00021 values[ii++] = bub * polx[i] * poly[j] * polz[k]; 00022 00023 return ii; 00024 } 00025 00026 template <typename Sx, typename Sy, typename Sz, typename T1, typename T2, typename T3> 00027 static void CalcSplitted (int n, Sx x, Sy y, Sz z, T1 & val1, T2 & val2, T3 & val3) 00028 { 00029 Sx bub1 = sqr (1-y-z) - sqr(x); 00030 ScaledLegendrePolynomialMult (n-4, x, 1-y-z, bub1, val1); 00031 ScaledLegendrePolynomialMult (n-4, 2*y - (1-z) , (1-z), y, val2); 00032 00033 // LegendrePolynomialMult (n-4, 2*z-1, z, val3); 00034 LegendrePolynomial leg; 00035 leg.EvalMult (n-4, 2*z-1, z, val3); 00036 00037 /* 00038 ScaledLegendrePolynomial (n-4, x, 1-y-z, val1); 00039 // Sx bub1 = (1-x-y-z) * (1+x-y-z); 00040 Sx bub1 = sqr (1-y-z) - sqr(x); 00041 for (int i = 0; i <= n-4; i++) 00042 val1[i] *= bub1; 00043 00044 ScaledLegendrePolynomial (n-4, 2*y - (1-z) , (1-z), val2); 00045 for (int i = 0; i <= n-4; i++) 00046 val2[i] *= y; 00047 00048 LegendrePolynomial (n-4, 2*z-1, val3); 00049 for (int i = 0; i <= n-4; i++) 00050 val3[i] *= z; 00051 */ 00052 } 00053 }; 00054 00055 00056 00057 class TetShapesInnerJacobi 00058 { 00059 public: 00060 template <typename Sx, typename Sy, typename Sz, typename T> 00061 static int Calc (int n, Sx x, Sy y, Sz z, T & values) 00062 { 00063 int ii = 0, j, k; 00064 ArrayMem<Sx, 20> polx(n+1), poly(n+1), polz(n+1); 00065 00066 Sx bub = y * z * (1-x-y-z) * (1+x-y-z); 00067 ScaledJacobiPolynomial (n-4, x, (1-y-z), 2, 2, polx); 00068 00069 for (int ix = 0; ix <= n-4; ix++) 00070 { 00071 ScaledJacobiPolynomial (n-4, (2*y-1+z),(1-z), 2*ix+5, 2, poly); 00072 JacobiPolynomial (n-4, 2*z-1, 2*ix+5, 2, polz); 00073 00074 for (j = 0; j <= n-4-ix; j++) 00075 for (k = 0; k <= n-4-ix-j; k++) 00076 values[ii++] = bub * polx[ix] * poly[j] * polz[k]; 00077 } 00078 return ii; 00079 } 00080 }; 00081 00082 class TetShapesFaceLegendre 00083 { 00084 public: 00085 template <typename Sx, typename Sy, typename Sz, typename T> 00086 static int Calc (int n, Sx x, Sy y, Sz z, T & values) 00087 { 00088 ArrayMem<Sx, 20> polx(n+1), poly(n+1); 00089 00090 ScaledLegendrePolynomial (n-3, x, 1-y-z, polx); 00091 ScaledLegendrePolynomial (n-3, 2*y-(1-z),(1-z), poly); 00092 Sx bub = (1-x-y-z) * (1+x-y-z)*y; 00093 00094 int ii = 0; 00095 for (int i = 0; i <= n-3; i++) 00096 for (int j = 0; j <= n-3-i; j++) 00097 values[ii++] = bub * polx[i] * poly[j]; 00098 00099 return ii; 00100 } 00101 template <typename Sx, typename Sy, typename Sz, typename T> 00102 static void CalcSplitted (int n, Sx x, Sy y, Sz z, T & val1, T & val2) 00103 { 00104 Sx bub1 = sqr (1-y-z) - sqr (x); 00105 ScaledLegendrePolynomialMult (n-3, x, 1-y-z, bub1, val1); 00106 ScaledLegendrePolynomialMult (n-3, 2*y-(1-z),(1-z), y, val2); 00107 /* 00108 ScaledLegendrePolynomial (n-3, x, 1-y-z, val1); 00109 00110 // Sx bub1 = (1-x-y-z) * (1+x-y-z); 00111 Sx bub1 = sqr (1-y-z) - sqr (x); 00112 for (int i = 0; i <= n-3; i++) 00113 val1[i] *= bub1; 00114 00115 ScaledLegendrePolynomial (n-3, 2*y-(1-z),(1-z), val2); 00116 for (int i = 0; i <= n-3; i++) 00117 val2[i] *= y; 00118 */ 00119 } 00120 }; 00121 00122 00123 00124 00125 class TetShapesFaceJacobi 00126 { 00127 public: 00128 template <typename Sx, typename Sy, typename Sz, typename T> 00129 static int Calc (int n, Sx x, Sy y, Sz z, T & values) 00130 { 00131 int ii = 0; 00132 ArrayMem<Sx, 20> polx(n+1), poly(n+1); 00133 00134 Sx bub = y * (1-x-y-z) * (1+x-y-z); 00135 ScaledJacobiPolynomial (n-3, x, 1-y-z, 2, 2, polx); 00136 00137 for (int ix = 0; ix <= n-3; ix++) 00138 { 00139 ScaledJacobiPolynomial (n-3, 2*y-1+z, 1-z, 2*ix+5, 2, poly); 00140 for (int j = 0; j <= n-3-ix; j++) 00141 values[ii++] = bub * polx[ix] * poly[j]; 00142 } 00143 return ii; 00144 } 00145 }; 00146 00147 00148 class TetShapesFaceOpt1 00149 { 00150 public: 00151 template <typename Sx, typename Sy, typename Sz, typename T> 00152 static int Calc (int n, Sx x, Sy y, Sz z, T & values) 00153 { 00154 int nd = Calc1 (n, x, y, z, values); 00155 ArrayMem<Sx, 100> hvalues(nd); 00156 00157 Sx lam1 = 0.5 * (1 + x - y - z); 00158 Sx lam2 = 0.5 * (1 - x - y - z); 00159 Sx lam3 = y; 00160 Sx lam4 = z; 00161 00162 Sx hlam1 = 0; 00163 Sx hlam2 = lam2; 00164 Sx hlam3 = lam3; 00165 Sx hlam4 = lam4 + lam1; 00166 00167 Sx hx = hlam1 - hlam2; 00168 Sx hy = hlam3; 00169 Sx hz = hlam4; 00170 00171 Sx frac; 00172 if (hlam4 < 1e-12) 00173 frac = 0.0; 00174 else 00175 frac = lam4 / hlam4; 00176 00177 Calc1 (n, hx, hy, hz, hvalues); 00178 for (int i = 0; i < nd; i++) 00179 values[i] -= frac * hvalues[i]; 00180 00181 return nd; 00182 } 00183 00184 template <typename Sx, typename Sy, typename Sz, typename T> 00185 static int Calc1 (int n, Sx x, Sy y, Sz z, T & values) 00186 { 00187 int nd = Calc2 (n, x, y, z, values); 00188 ArrayMem<Sx, 100> hvalues(nd); 00189 00190 Sx lam1 = 0.5 * (1 + x - y - z); 00191 Sx lam2 = 0.5 * (1 - x - y - z); 00192 Sx lam3 = y; 00193 Sx lam4 = z; 00194 00195 Sx hlam1 = lam1; 00196 Sx hlam2 = 0; 00197 Sx hlam3 = lam3; 00198 Sx hlam4 = lam4 + lam2; 00199 00200 Sx hx = hlam1 - hlam2; 00201 Sx hy = hlam3; 00202 Sx hz = hlam4; 00203 00204 Sx frac; 00205 if (hlam4 < 1e-12) 00206 frac = 0.0; 00207 else 00208 frac = lam4 / hlam4; 00209 00210 Calc2 (n, hx, hy, hz, hvalues); 00211 for (int i = 0; i < nd; i++) 00212 values[i] -= frac * hvalues[i]; 00213 00214 return nd; 00215 } 00216 00217 template <typename Sx, typename Sy, typename Sz, typename T> 00218 static int Calc2 (int n, Sx x, Sy y, Sz z, T & values) 00219 { 00220 Sx * hp = &values[0]; 00221 int nd = Calc3 (n, x, y, z, hp); 00222 ArrayMem<Sx, 100> hvalues(nd); 00223 00224 Sx lam1 = 0.5 * (1 + x - y - z); 00225 Sx lam2 = 0.5 * (1 - x - y - z); 00226 Sx lam3 = y; 00227 Sx lam4 = z; 00228 00229 Sx hlam1 = lam1; 00230 Sx hlam2 = lam2; 00231 Sx hlam3 = 0; 00232 Sx hlam4 = lam4 + lam3; 00233 00234 Sx hx = hlam1 - hlam2; 00235 Sx hy = hlam3; 00236 Sx hz = hlam4; 00237 00238 Sx frac; 00239 if (hlam4 < 1e-12) 00240 frac = 0.0; 00241 else 00242 frac = lam4 / hlam4; 00243 00244 hp = &hvalues[0]; 00245 Calc3 (n, hx, hy, hz, hp); 00246 for (int i = 0; i < nd; i++) 00247 values[i] -= frac * hvalues[i]; 00248 00249 return nd; 00250 } 00251 00252 00253 00254 00255 00256 00257 template <typename Sx, typename Sy, typename Sz, typename T> 00258 static int Calc3 (int n, Sx x, Sy y, Sz z, T & values) 00259 { 00260 int ii = 0; 00261 ArrayMem<Sx, 20> polx(n+1), poly(n+1); 00262 00263 const IntegrationRule & rule = SelectIntegrationRule (ET_TRIG, n+2); 00264 00265 for (int ix = 0; ix <= n-3; ix++) 00266 for (int j = 0; j <= n-3-ix; j++) 00267 values[ii++] = 0; 00268 for (int i = 0; i < rule.GetNIP(); i++) 00269 { 00270 ii = 0; 00271 const IntegrationPoint & ip = rule[i]; 00272 00273 Sx hx = x + z * (-1+2*ip(0)+ip(1)); 00274 Sy hy = y + z * ip(1); 00275 00276 Sx bub = hy * (1-hx-hy) * (1+hx-hy); 00277 ScaledJacobiPolynomial (n-3, hx, 1-hy, 2, 2, polx); 00278 00279 Sx fac = 2 * bub * ip.Weight(); // / (z*z); 00280 00281 for (int ix = 0; ix <= n-3; ix++) 00282 { 00283 ScaledJacobiPolynomial (n-3, 2*hy-1, 1, 2*ix+5, 2, poly); 00284 for (int j = 0; j <= n-3-ix; j++) 00285 values[ii++] += fac * polx[ix] * poly[j]; 00286 } 00287 } 00288 return ii; 00289 } 00290 }; 00291 00292 00293 00294 00295 class TetShapesFaceOpt2 00296 { 00297 public: 00298 template <typename Sx, typename Sy, typename Sz, typename T> 00299 static int Calc (int n, Sx x, Sy y, Sz z, T & values) 00300 { 00301 int ii = 0, i, j; 00302 ArrayMem<Sx, 20> polx(n+1), poly(n+1); 00303 00304 const IntegrationRule & rule = SelectIntegrationRule (ET_TRIG, n+2); 00305 00306 for (int ix = 0; ix <= n-3; ix++) 00307 for (j = 0; j <= n-3-ix; j++) 00308 values[ii++] = 0; 00309 for (i = 0; i < rule.GetNIP(); i++) 00310 { 00311 ii = 0; 00312 const IntegrationPoint & ip = rule[i]; 00313 00314 Sx hx = x + z * (-1+2*ip(0)+ip(1)); 00315 Sy hy = y + z * ip(1); 00316 00317 //Sx bub = hy * (1-hx-hy) * (1+hx-hy); 00318 ScaledJacobiPolynomial (n-3, hx, 1-hy, 2, 2, polx); 00319 00320 Sx fac = 2 * ip.Weight(); // / (z*z); 00321 00322 for (int ix = 0; ix <= n-3; ix++) 00323 { 00324 ScaledJacobiPolynomial (n-3, 2*hy-1, 1, 2*ix+5, 2, poly); 00325 for (j = 0; j <= n-3-ix; j++) 00326 values[ii++] += fac * polx[ix] * poly[j]; 00327 } 00328 } 00329 int jj = ii; 00330 ArrayMem<Sx, 100> hvalues(jj); 00331 for (int i = 0; i < ii; i++) 00332 hvalues[i] = values[i]; 00333 ii = 0; jj = 0; 00334 for (int k = 0; k <= n-3; k++) 00335 { 00336 for (int m = 0; m <= n-3-k; m++) 00337 { 00338 values[ii++] = y*(1-x-y-z)*(1+x-y-z)*hvalues[jj++]; 00339 } 00340 } 00341 00342 return ii; 00343 } 00344 }; 00345 00346 00347 00348 00349 00350 00351 00352 00353 00354 00355 /* 00357 virtual void CalcShape (const IntegrationPoint & ip, 00358 FlatVector<> shape) const 00359 { 00360 double lam1 = ip(0); 00361 double lam2 = 1-ip(0)-ip(1); 00362 double lam3 = ip(1); 00363 double bub = lam1*lam2*lam3; 00364 00365 ArrayMem<double, 50> polx(order-2); 00366 ArrayMem<double, 50> poly(order-2); 00367 00368 double fac = 1; 00369 for (int i = 0; i <= order-3; i++) 00370 { 00371 polx[i] *= fac; 00372 fac *= (1-lam3); 00373 } 00374 00375 int ii = 0; 00376 for (int ix = 0; ix <= order-3; ix++) 00377 { 00378 JacobiPolynomial (order-3, lam3-(lam1+lam2), 2*ix+5, 2, poly); 00379 for (int iy = 0; iy <= order-3-ix; iy++) 00380 shape(ii++) = bub*polx[ix] * poly[iy]; 00381 } 00382 } 00383 00384 */ 00385 00386 00387 }