NGSolve
4.9
|
00001 namespace ngfem 00002 { 00003 00008 class TrigShapesInnerLegendre 00009 { 00010 public: 00012 template <typename Sx, typename Sy, typename T> 00013 static void CalcSplitted (int n, Sx x, Sy y, T & val1, T & val2) 00014 { 00015 LegendrePolynomial leg; 00016 00017 Sx bub1 = (1-x-y)*(1+x-y); 00018 // ScaledLegendrePolynomialMult (n-3, x, 1-y, bub1, val1); 00019 leg.EvalScaledMult (n-3, x, 1-y, bub1, val1); 00020 00021 //LegendrePolynomialMult (n-3, 2*y-1, y, val2); 00022 leg.EvalMult (n-3, 2*y-1, y, val2); 00023 00024 /* 00025 ScaledLegendrePolynomial (n-3, x, 1-y, val1); 00026 for (int i = 0; i <= n-3; i++) 00027 val1[i] *= bub1; 00028 00029 LegendrePolynomial (n-3, 2*y-1, val2); 00030 for (int i = 0; i <= n-3; i++) 00031 val2[i] *= y; 00032 */ 00033 } 00034 00035 00036 template <int n, typename Sx, typename Sy, typename T> 00037 static void CalcSplitted (Sx x, Sy y, T & val1, T & val2) 00038 { 00039 Sx bub1 = (1-x-y)*(1+x-y); 00040 // ScaledLegendrePolynomialMult (n-3, x, 1-y, bub1, val1); 00041 // LegendrePolynomialFO<n-3>::EvalScaledMult (x, 1-y, bub1, val1); 00042 LegendrePolynomial leg; 00043 leg.EvalScaledMultFO<n-3> (x, 1-y, bub1, val1); 00044 00045 // LegendrePolynomialMult (n-3, 2*y-1, y, val2); 00046 // LegendrePolynomialFO<n-3>::EvalMult (2*y-1, y, val2); 00047 leg.EvalMultFO<n-3> (2*y-1, y, val2); 00048 } 00049 00050 00052 template < typename Sx, typename Sy, typename T> 00053 static int Calc (int n, Sx x, Sy y, T & values) 00054 { 00055 ArrayMem<Sx, 20> polx(n-2), poly(n-2); 00056 00057 ScaledLegendrePolynomial (n-3, x, 1-y, polx); 00058 LegendrePolynomial (n-3, 2*y-1, poly); 00059 Sx bub = y * (1-x-y) * (1+x-y); 00060 00061 int ii = 0; 00062 for (int i = 0; i <= n-3; i++) 00063 for (int j = 0; j <= n-3-i; j++) 00064 values[ii++] = bub * polx[i] * poly[j]; 00065 00066 return ii; 00067 } 00068 00069 template <int n, typename Sx, typename Sy, typename T> 00070 static int Calc (Sx x, Sy y, T & values) 00071 { 00072 // ArrayMem<Sx, 20> polx(n-2), poly(n-2); 00073 Sx polx[n], poly[n]; 00074 00075 /* 00076 ScaledLegendrePolynomial (n-3, x, 1-y, polx); 00077 LegendrePolynomial (n-3, 2*y-1, poly); 00078 Sx bub = y * (1-x-y) * (1+x-y); 00079 */ 00080 CalcSplitted<n> (x, y, polx, poly); 00081 int ii = 0; 00082 for (int i = 0; i <= n-3; i++) 00083 for (int j = 0; j <= n-3-i; j++) 00084 values[ii++] = /* bub * */ polx[i] * poly[j]; 00085 00086 return ii; 00087 } 00088 00089 00090 00091 00092 00094 template <typename Sx, typename Sy, typename Sf, typename T> 00095 static int CalcMult (int n, Sx x, Sy y, Sf & fac, T & values) 00096 { 00097 ArrayMem<Sx, 20> polx(n-2), poly(n-2); 00098 00099 ScaledLegendrePolynomial (n-3, x, 1-y, polx); 00100 LegendrePolynomial (n-3, 2*y-1, poly); 00101 00102 Sx bub = fac * y * (1-x-y) * (1+x-y); 00103 00104 int ii = 0; 00105 for (int i = 0; i <= n-3; i++) 00106 for (int j = 0; j <= n-3-i; j++) 00107 values[ii++] = bub * polx[i] * poly[j]; 00108 00109 return ii; 00110 } 00111 00112 00113 }; 00114 00115 00120 class TrigShapesInnerJacobi 00121 { 00122 public: 00124 template <typename Sx, typename Sy, typename T> 00125 static int Calc (int n, Sx x, Sy y, T & values) 00126 { 00127 int ii = 0; 00128 ArrayMem<Sx, 20> polx(n+1), poly(n+1); 00129 00130 Sx bub = y * (1-x-y) * (1+x-y); 00131 ScaledJacobiPolynomial (n-3, x, 1-y, 2, 2, polx); 00132 00133 for (int ix = 0; ix <= n-3; ix++) 00134 { 00135 JacobiPolynomial (n-3, 2*y-1, 2*ix+5, 2, poly); 00136 for (int j = 0; j <= n-3-ix; j++) 00137 values[ii++] = bub * polx[ix] * poly[j]; 00138 } 00139 return ii; 00140 } 00141 00142 00143 }; 00144 00145 00146 //*************************************MONOMIAL EXTENSION*********************************************************** 00147 00161 class TrigExtensionMonomial 00162 { 00163 public: 00165 template <class Sx, class Sy, class T> 00166 inline static int CalcTrigExt (int n, Sx x, Sy y, T & values) 00167 { 00168 Sy fy = (1-y)*(1-y); 00169 Sx p3 = 0; 00170 Sx p2 = -1; 00171 Sx p1 = x; 00172 00173 for (int j=2; j<=n; j++) 00174 { 00175 p3=p2; p2=p1; 00176 p1=( (2*j-3) * x * p2 - (j-3) * fy * p3) / j; 00177 values[j-2] = p1; 00178 } 00179 return n-1; 00180 } 00181 00182 00184 template <class T> 00185 inline static int CalcTrigExtDeriv (int n, double x, double y, T & values) 00186 { 00187 double fy = (1-y)*(1-y); 00188 double p3 = 0, p3x = 0, p3y = 0; 00189 double p2 = -1, p2x = 0, p2y = 0; 00190 double p1 = x, p1x = 1, p1y = 0; 00191 00192 for (int j=2; j<=n; j++) 00193 { 00194 p3=p2; p3x = p2x; p3y = p2y; 00195 p2=p1; p2x = p1x; p2y = p1y; 00196 double c1 = (2.0*j-3) / j; 00197 double c2 = (j-3.0) / j; 00198 00199 p1 = c1 * x * p2 - c2 * fy * p3; 00200 p1x = c1 * p2 + c1 * x * p2x - c2 * fy * p3x; 00201 p1y = c1 * x * p2y - (c2 * 2 * (y-1) * p3 + c2 * fy * p3y); 00202 values (j-2, 0) = p1x; 00203 values (j-2, 1) = p1y; 00204 } 00205 return n-1; 00206 } 00207 00209 template <class Sx, class T> 00210 inline static int Calc (int n, Sx x, T & values) 00211 { 00212 Sx p3 = 0; 00213 Sx p2 = -1; 00214 Sx p1 = x; 00215 00216 for (int j=2; j<=n; j++) 00217 { 00218 p3=p2; p2=p1; 00219 p1=( (2*j-3) * x * p2 - (j-3) * p3) / j; 00220 values[j-2] = p1; 00221 } 00222 return n-1; 00223 00224 } 00225 00227 template <class T> 00228 inline static int CalcDeriv (int n, double x, T & values) 00229 { 00230 double p1 = 1.0, p2 = 0.0, p3; 00231 00232 for (int j=1; j<=n-1; j++) 00233 { 00234 p3 = p2; p2 = p1; 00235 p1 = ((2.0*j-1.0)*x*p2 - (j-1.0)*p3) / j; 00236 values[j-1] = p1; 00237 } 00238 return n-1; 00239 00240 } 00241 }; 00242 00243 00244 // *********************** OPTIMAL EXTENSION ***************************************** 00258 class TrigExtensionOptimal 00259 { 00260 enum { SIZE = 1000 }; 00261 static double coefs[SIZE][5]; 00262 static bool initialized; 00263 public: 00264 00265 TrigExtensionOptimal (); 00266 00267 template <class Sx, class Sy, class T> 00268 inline static int CalcTrigExt (int p, Sx x, Sy y, T & values) 00269 { 00270 Sx p1(0.0), p2(0.0), p3(0.0), p4(0.0), p5(0.0); 00271 Sx bub = (1.+x-y)*(x+y-1); 00272 switch(p) 00273 { 00274 default: 00275 case 5: 00276 00277 { 00278 p5 = -1./24. * bub*x * (-9.+21.*x*x-14.*y+35.*y*y); 00279 values[3] = p5; 00280 } 00281 00282 case 4: 00283 00284 { 00285 p4 = -0.125 * bub * (-1.+5.*x*x-2.*y+3.*y*y); 00286 values[2] = p4; 00287 } 00288 00289 case 3: 00290 00291 { 00292 p3 = -0.5 * bub * x; 00293 values[1] = p3; 00294 } 00295 00296 case 2: 00297 00298 { 00299 p2 = -0.5 * bub; 00300 values[0] = p2; 00301 } 00302 case 1: 00303 case 0: 00304 ; 00305 } 00306 00307 for (int j = 6; j <= p; j++) 00308 { 00309 p1 = 00310 -coefs[j-6][0] * p2 00311 + coefs[j-6][1] * x * p3 00312 - (coefs[j-6][2]+coefs[j-6][3]*(x*x-y*y)) * p4 00313 + coefs[j-6][4] * x * p5; 00314 p2 = p3; p3 = p4; p4 = p5; p5 = p1; 00315 values[j-2] = p1; 00316 } 00317 00318 return p-1; 00319 } 00320 00321 template <class T> 00322 inline static int CalcTrigExtDeriv (int n, double x, double y, T & values) 00323 { 00324 ArrayMem<AutoDiff<2>,10> ad_values(n-1); 00325 AutoDiff<2> ad_x(x, 0); 00326 AutoDiff<2> ad_y(y, 1); 00327 00328 CalcTrigExt (n, ad_x, ad_y, ad_values); 00329 00330 for (int i = 0; i < n-1; i++) 00331 for (int j = 0; j < 2; j++) 00332 values(i,j) = ad_values[i].DValue(j); 00333 return n-1; 00334 } 00335 00336 00337 template <class Sx, class T> 00338 inline static int Calc (int n, Sx x, T & values) 00339 { 00340 Sx y = 0.0; 00341 return CalcTrigExt (n, x, y, values); 00342 } 00343 00344 template <class T> 00345 inline static int CalcDeriv (int n, double x, T & values) 00346 { 00347 return CalcTrigExtDeriv (n, x, 0.0, values); 00348 } 00349 }; 00350 00351 00352 00353 // ************************ MINIMAL EXTENSION ************************* 00354 00355 00359 class TrigExtensionMin 00360 { 00361 public: 00362 template <class Sx, class Sy, class T> 00363 inline static int CalcTrigExt (int n, Sx x, Sy y, T & values) 00364 { 00365 #ifdef ABC 00366 // TrigExtensionOptimal::CalcTrigExt (n, x, y, values); 00367 00368 ArrayMem<Sx, 20> polx(n+1); 00369 00370 const IntegrationRule & rule = 00371 GetIntegrationRules().SelectIntegrationRule (ET_SEGM, n+2); 00372 00373 Sx bub = (1-x-y) * (1+x-y); 00374 for (int i = 0; i <= n-2; i++) 00375 values[i] = 0; 00376 00377 for (int i = 0; i < rule.GetNIP(); i++) 00378 { 00379 const IntegrationPoint & ip = rule.GetIP(i); 00380 00381 Sx hx = x + y * (2.0*ip(0)-1.0); 00382 00383 JacobiPolynomial (n-2, hx, 1, 1, polx); 00384 // LegendrePolynomial (n-2, hx, polx); 00385 00386 Sx fac = bub * ip.Weight(); // / (z*z); 00387 00388 for (int j = 0; j <= n-2; j++) 00389 values[j] += fac * polx[j]; 00390 } 00391 00392 00393 00394 // IntegratedLegendreMonomialExt::CalcTrigExt (n, x, y, values); 00395 00396 /* 00397 Sy fy = (1-y)*(1-y); 00398 Sx p3 = 0; 00399 Sx p2 = -1; 00400 Sx p1 = x; 00401 00402 for (int j=2; j<=n; j++) 00403 { 00404 p3=p2; p2=p1; 00405 p1=( (2*j-3) * x * p2 - (j-3) * fy * p3) / j; 00406 values[j-2] = p1; 00407 } 00408 */ 00409 00410 /* 00411 Sx fac = 1; 00412 for (int j = n-2; j >= 0; j--) 00413 { 00414 values[j] = fac * values[j]; 00415 fac = fac * (1-y); 00416 } 00417 */ 00418 00419 for (int j = 0; j < n-3; j++) 00420 values[j] *= VertexExtensionOptimal<2>::Calc (n-2-j, 1-y); 00421 00422 return n-1; 00423 #endif 00424 00425 if (y == 0.0) y += 1e-8; 00426 ArrayMem<Sx, 20> pol1(n+2), pol2(n+2); 00427 JacobiPolynomial (n+1, x-y, 1, 1, pol1); 00428 JacobiPolynomial (n+1, x+y, 1, 1, pol2); 00429 for (int i = 0; i <= n-2; i++) 00430 values[i] = (pol2[i+1]-pol1[i+1]) / (2*y); 00431 00432 /* 00433 ScaledJacobiPolynomial (n-2, x, 1-y, 2, 2, values); 00434 for (int i = 0; i < n-1; i++) 00435 values[i] *= 0.5 * (i+4) * (1-x*x); 00436 */ 00437 00438 00439 00440 LowEnergyVertexPolynomials2D (n, 1-2*y, pol1); 00441 for (int i = 0; i <= n-2; i++) 00442 values[i] *= pol1[n-2-i] * (1-x-y) * (1+x-y); 00443 00444 return n-1; 00445 } 00446 00447 template <class T> 00448 inline static int CalcTrigExtDeriv (int n, double x, double y, T & values) 00449 { 00450 ArrayMem<AutoDiff<2>,10> ad_values(n-1); 00451 AutoDiff<2> ad_x(x, 0); 00452 AutoDiff<2> ad_y(y, 1); 00453 00454 CalcTrigExt (n, ad_x, ad_y, ad_values); 00455 00456 for (int i = 0; i < n-1; i++) 00457 for (int j = 0; j < 2; j++) 00458 values(i,j) = ad_values[i].DValue(j); 00459 return n-1; 00460 } 00461 00462 00463 00464 template <class Sx, class T> 00465 inline static int Calc (int n, Sx x, T & values) 00466 { 00467 /* 00468 Sx y = 0.0; 00469 CalcTrigExt (n, x, y, values); 00470 */ 00471 00472 JacobiPolynomial (n-2, x, 2, 2, values); 00473 for (int i = 0; i < n-1; i++) 00474 values[i] *= 0.5 * (i+4) * (1-x*x); 00475 00476 return n-1; 00477 00478 } 00479 00480 template <class T> 00481 inline static int CalcDeriv (int n, double x, T & values) 00482 { 00483 ArrayMem<AutoDiff<1>,10> ad_values(n-1); 00484 AutoDiff<1> ad_x(x, 0); 00485 00486 Calc (n, ad_x, ad_values); 00487 00488 for (int i = 0; i < n-1; i++) 00489 values[i] = ad_values[i].DValue(0); 00490 return n-1; 00491 00492 // return CalcTrigExtDeriv (n, x, 0.0, values); 00493 } 00494 }; 00495 00496 00497 }