1 #include <mystdlib.h> 2 3 #include "meshing.hpp" 4 #ifndef CURVEDELEMS_NEW 5 6 namespace netgen 7 { 8 9 10 // computes Gaussean integration formula on (0,1) with n points 11 // in: Numerical algs in C (or, was it the Fortran book ?) ComputeGaussRule(int n,ARRAY<double> & xi,ARRAY<double> & wi)12 void ComputeGaussRule (int n, ARRAY<double> & xi, ARRAY<double> & wi) 13 { 14 xi.SetSize (n); 15 wi.SetSize (n); 16 17 int m = (n+1)/2; 18 double p1, p2, p3; 19 double pp, z, z1; 20 for (int i = 1; i <= m; i++) 21 { 22 z = cos ( M_PI * (i - 0.25) / (n + 0.5)); 23 while(1) 24 { 25 p1 = 1; p2 = 0; 26 for (int j = 1; j <= n; j++) 27 { 28 p3 = p2; p2 = p1; 29 p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j; 30 } 31 // p1 is legendre polynomial 32 33 pp = n * (z*p1-p2) / (z*z - 1); 34 z1 = z; 35 z = z1-p1/pp; 36 37 if (fabs (z - z1) < 1e-14) break; 38 } 39 40 xi[i-1] = 0.5 * (1 - z); 41 xi[n-i] = 0.5 * (1 + z); 42 wi[i-1] = wi[n-i] = 1.0 / ( (1 - z * z) * pp * pp); 43 } 44 } 45 46 47 48 // ---------------------------------------------------------------------------- 49 // PolynomialBasis 50 // ---------------------------------------------------------------------------- 51 52 CalcLegendrePolynomials(double x)53 void PolynomialBasis :: CalcLegendrePolynomials (double x) 54 { 55 double p1 = 1.0, p2 = 0.0, p3; 56 57 lp[0] = 1.0; 58 59 for (int j=1; j<=order; j++) 60 { 61 p3=p2; p2=p1; 62 p1=((2.0*j-1.0)*(2*x-1)*p2-(j-1.0)*p3)/j; 63 lp[j] = p1; 64 } 65 } 66 67 CalcScaledLegendrePolynomials(double x,double t)68 void PolynomialBasis :: CalcScaledLegendrePolynomials (double x, double t) 69 { 70 double p1 = 1.0, p2 = 0.0, p3; 71 72 lp[0] = 1.0; 73 74 double x2mt = 2*x-t; 75 double t2 = t*t; 76 77 for (int j=1; j<=order; j++) 78 { 79 p3=p2; p2=p1; 80 p1=((2.0*j-1.0)*x2mt*p2-t2*(j-1.0)*p3)/j; 81 lp[j] = p1; 82 } 83 } 84 85 CalcDLegendrePolynomials(double x)86 void PolynomialBasis :: CalcDLegendrePolynomials (double x) 87 { 88 double p1 = 0., p2 = 0., p3; 89 90 dlp[0] = 0.; 91 92 for (int j = 1; j <= order-1; j++) 93 { 94 p3=p2; p2=p1; 95 p1=((2.*j-1.)*(2*lp[j-1]+(2*x-1)*p2)-(j-1.)*p3)/j; 96 dlp[j] = p1; 97 } 98 } 99 100 CalcF(double x)101 void PolynomialBasis :: CalcF (double x) 102 { 103 CalcLegendrePolynomials (x); 104 105 for (int j = 0; j<=order-2; j++) 106 f[j] = (lp[j+2]-lp[j])/(2.0*(j+1)+1)/2.0; 107 } 108 109 CalcDf(double x)110 void PolynomialBasis :: CalcDf (double x) 111 { 112 CalcLegendrePolynomials (x); 113 114 for (int j = 0; j <= order-2; j++) 115 df[j] = lp[j+1]; 116 } 117 118 CalcFDf(double x)119 void PolynomialBasis :: CalcFDf (double x) 120 { 121 CalcLegendrePolynomials (x); 122 123 for (int j = 0; j<=order-2; j++) 124 { 125 f[j] = (lp[j+2]-lp[j])/(2.0*(j+1)+1)/2.0; 126 df[j] = lp[j+1]; 127 } 128 } 129 130 CalcDDf(double x)131 void PolynomialBasis :: CalcDDf (double x) 132 { 133 CalcLegendrePolynomials (x); 134 CalcDLegendrePolynomials (x); 135 136 for (int j = 0; j <= order-2; j++) ddf[j] = dlp[j+1]; 137 } 138 139 CalcFScaled(double x,double t)140 void PolynomialBasis :: CalcFScaled (double x, double t) 141 { 142 double tt = t*t; 143 double x2mt = 2*x-t; 144 145 146 147 148 double p1 = 0.5*x2mt, p2 = -0.5, p3 = 0.0; 149 for (int j=2; j<=order; j++) 150 { 151 p3=p2; p2=p1; 152 p1=((2.0*j-3.0) * x2mt * p2 - tt*(j-3.0)*p3)/j; 153 f[j-2] = p1; 154 } 155 156 /* 157 CalcF (x/t); 158 double fac = t*t; 159 for (int j = 0; j<=order-2; j++, fac*=t) 160 f[j] *= fac; 161 */ 162 } 163 CalcFScaled(int p,double x,double t,double * values)164 void PolynomialBasis :: CalcFScaled (int p, double x, double t, double * values) 165 { 166 double tt = t*t; 167 double x2mt = 2*x-t; 168 169 double p1 = 0.5*x2mt, p2 = -0.5, p3 = 0.0; 170 for (int j=2; j<=p; j++) 171 { 172 p3=p2; p2=p1; 173 p1=((2.0*j-3.0) * x2mt * p2 - tt*(j-3.0)*p3)/j; 174 values[j-2] = p1; 175 } 176 177 /* 178 CalcF (x/t); 179 double fac = t*t; 180 for (int j = 0; j<=order-2; j++, fac*=t) 181 f[j] *= fac; 182 */ 183 } 184 185 186 187 188 // ---------------------------------------------------------------------------- 189 // BaseFiniteElement1D 190 // ---------------------------------------------------------------------------- 191 192 CalcVertexShapes()193 void BaseFiniteElement1D :: CalcVertexShapes () 194 { 195 vshape[0] = xi(0); 196 vshape[1] = 1-xi(0); 197 198 vdshape[0] = 1; 199 vdshape[1] = -1; 200 201 /* 202 if (edgeorient == -1) 203 { 204 Swap (vshape[0], vshape[1]); 205 Swap (vdshape[0], vdshape[1]); 206 } 207 */ 208 209 } 210 211 CalcEdgeShapes()212 void BaseFiniteElement1D :: CalcEdgeShapes () 213 { 214 b.SetOrder (edgeorder); 215 if (edgeorient == 1) 216 b.CalcFDf( 1-xi(0) ); 217 else 218 b.CalcFDf( xi(0) ); 219 220 for (int k = 2; k <= edgeorder; k++) 221 { 222 eshape[k-2] = b.GetF(k); 223 edshape[k-2] = -b.GetDf(k); 224 } 225 } 226 227 CalcEdgeLaplaceShapes()228 void BaseFiniteElement1D :: CalcEdgeLaplaceShapes () 229 { 230 b.SetOrder (edgeorder); 231 if (edgeorient == 1) 232 b.CalcDDf( 1-xi(0) ); 233 else 234 b.CalcDDf( xi(0) ); 235 236 for (int k = 2; k <= edgeorder; k++) 237 eddshape[k-2] = b.GetDDf(k); 238 } 239 240 241 242 243 // ---------------------------------------------------------------------------- 244 // BaseFiniteElement2D 245 // ---------------------------------------------------------------------------- 246 247 SetElementNumber(int aelnr)248 void BaseFiniteElement2D :: SetElementNumber (int aelnr) 249 { 250 int locmaxedgeorder = -1; 251 252 BaseFiniteElement<2> :: SetElementNumber (aelnr); 253 const Element2d & elem = mesh[(SurfaceElementIndex) (elnr-1)]; 254 top.GetSurfaceElementEdges (elnr, &(edgenr[0]), &(edgeorient[0])); 255 facenr = top.GetSurfaceElementFace (elnr); 256 faceorient = top.GetSurfaceElementFaceOrientation (elnr); 257 258 for (int v = 0; v < nvertices; v++) 259 vertexnr[v] = elem[v]; 260 261 for (int e = 0; e < nedges; e++) 262 { 263 edgeorder[e] = curv.GetEdgeOrder (edgenr[e]-1); // 1-based 264 locmaxedgeorder = max2 (edgeorder[e], locmaxedgeorder); 265 } 266 267 faceorder = curv.GetFaceOrder (facenr-1); // 1-based 268 CalcNFaceShapes (); 269 270 if (locmaxedgeorder > maxedgeorder) 271 { 272 maxedgeorder = locmaxedgeorder; 273 eshape.SetSize(nedges * (maxedgeorder-1)); 274 edshape.SetSize(nedges * (maxedgeorder-1)); 275 } 276 277 if (faceorder > maxfaceorder) 278 { 279 maxfaceorder = faceorder; 280 fshape.SetSize( nfaceshapes ); // number of face shape functions 281 fdshape.SetSize( nfaceshapes ); 282 fddshape.SetSize ( nfaceshapes ); 283 } 284 }; 285 286 287 // ---------------------------------------------------------------------------- 288 // BaseFiniteElement3D 289 // ---------------------------------------------------------------------------- 290 291 SetElementNumber(int aelnr)292 void BaseFiniteElement3D :: SetElementNumber (int aelnr) 293 { 294 int locmaxedgeorder = -1; 295 int locmaxfaceorder = -1; 296 int v, f, e; 297 298 BaseFiniteElement<3> :: SetElementNumber (aelnr); 299 Element elem = mesh[(ElementIndex) (elnr-1)]; 300 top.GetElementEdges (elnr, &(edgenr[0]), &(edgeorient[0])); 301 top.GetElementFaces (elnr, &(facenr[0]), &(faceorient[0])); 302 303 for (v = 0; v < nvertices; v++) 304 vertexnr[v] = elem[v]; 305 306 for (f = 0; f < nfaces; f++) 307 { 308 surfacenr[f] = top.GetFace2SurfaceElement (facenr[f]); 309 // surfaceorient[f] = top.GetSurfaceElementFaceOrientation (surfacenr[f]); 310 } 311 312 for (e = 0; e < nedges; e++) 313 { 314 edgeorder[e] = curv.GetEdgeOrder (edgenr[e]-1); // 1-based 315 locmaxedgeorder = max (edgeorder[e], locmaxedgeorder); 316 } 317 318 for (f = 0; f < nfaces; f++) 319 { 320 faceorder[f] = curv.GetFaceOrder (facenr[f]-1); // 1-based 321 locmaxfaceorder = max (faceorder[f], locmaxfaceorder); 322 } 323 324 CalcNFaceShapes (); 325 326 if (locmaxedgeorder > maxedgeorder) 327 { 328 maxedgeorder = locmaxedgeorder; 329 eshape.SetSize(nedges * (maxedgeorder-1)); 330 edshape.SetSize(nedges * (maxedgeorder-1)); 331 } 332 333 if (locmaxfaceorder > maxfaceorder) 334 { 335 maxfaceorder = locmaxfaceorder; 336 fshape.SetSize( nfaces * (maxfaceorder-1) * (maxfaceorder-1)); // number of face shape functions 337 fdshape.SetSize( nfaces * (maxfaceorder-1) * (maxfaceorder-1)); 338 } 339 }; 340 341 342 343 344 // ---------------------------------------------------------------------------- 345 // FETrig 346 // ---------------------------------------------------------------------------- 347 348 SetReferencePoint(Point<2> axi)349 void FETrig :: SetReferencePoint (Point<2> axi) 350 { 351 BaseFiniteElement2D :: SetReferencePoint (axi); 352 lambda(0) = xi(0); 353 lambda(1) = xi(1); 354 lambda(2) = 1-xi(0)-xi(1); 355 356 dlambda(0,0) = 1; dlambda(0,1) = 0; 357 dlambda(1,0) = 0; dlambda(1,1) = 1; 358 dlambda(2,0) = -1; dlambda(2,1) = -1; 359 } 360 361 SetVertexSingularity(int v,int exponent)362 void FETrig :: SetVertexSingularity (int v, int exponent) 363 { 364 int i; 365 if (1-lambda(v) < EPSILON) return; 366 367 Point<3> lamold = lambda; 368 369 Vec<2> dfac; 370 371 double fac = pow(1-lambda(v),exponent-1); 372 373 for (i = 0; i < 2; i++) 374 { 375 dfac(i) = -(exponent-1)*pow(1-lambda(v),exponent-2)*dlambda(v,i); 376 dlambda(v,i) *= exponent * pow(1-lambda(v),exponent-1); 377 } 378 379 lambda(v) = 1-pow(1-lambda(v),exponent); 380 381 for (i = 0; i < nvertices; i++) 382 { 383 if (i == v) continue; 384 for (int j = 0; j < 2; j++) 385 dlambda(i,j) = dlambda(i,j) * fac + lamold(i) * dfac(j); 386 387 lambda(i) *= fac; 388 } 389 } 390 391 392 CalcVertexShapes()393 void FETrig :: CalcVertexShapes () 394 { 395 for (int v = 0; v < nvertices; v++) 396 { 397 vshape[v] = lambda(v); 398 vdshape[v](0) = dlambda(v,0); 399 vdshape[v](1) = dlambda(v,1); 400 } 401 } 402 403 CalcEdgeShapes()404 void FETrig :: CalcEdgeShapes () 405 { 406 int index = 0; 407 for (int e = 0; e < nedges; e++) 408 { 409 if (edgeorder[e] <= 1) continue; 410 411 int i0 = eledge[e][0]-1; 412 int i1 = eledge[e][1]-1; 413 414 double arg = lambda(i0) + lambda(i1); // = 1-lambda[i2]; 415 416 if (fabs(arg) < EPSILON) // division by 0 417 { 418 for (int k = 2; k <= edgeorder[e]; k++) 419 { 420 eshape[index] = 0; 421 edshape[index] = Vec<2>(0,0); 422 index++; 423 } 424 continue; 425 } 426 427 if (edgeorient[e] == -1) Swap (i0, i1); // reverse orientation 428 429 double xi = lambda(i1)/arg; 430 431 b1.SetOrder (edgeorder[e]); 432 b1.CalcFDf (xi); 433 434 double decay = arg; 435 double ddecay; 436 437 double l0 = lambda(i0); 438 double l0x = dlambda(i0,0); 439 double l0y = dlambda(i0,1); 440 441 double l1 = lambda(i1); 442 double l1x = dlambda(i1,0); 443 double l1y = dlambda(i1,1); 444 445 for (int k = 2; k <= edgeorder[e]; k++) 446 { 447 ddecay = k*decay; 448 decay *= arg; 449 450 eshape[index] = b1.GetF (k) * decay; 451 edshape[index] = Vec<2> 452 (b1.GetDf(k) * (l1x*arg - l1*(l0x+l1x)) * 453 decay / (arg * arg) + b1.GetF(k) * ddecay * (l0x+l1x), 454 b1.GetDf(k) * (l1y*arg - l1*(l0y+l1y)) * 455 decay / (arg * arg) + b1.GetF(k) * ddecay * (l0y+l1y)); 456 index++; 457 } 458 } 459 // (*testout) << "index = " << index << endl; 460 // (*testout) << "eshape = " << eshape << ", edshape = " << edshape << endl; 461 462 463 // eshape = 0.0; 464 // edshape = 0.0; 465 466 /* 467 int index = 0; 468 for (int e = 0; e < nedges; e++) 469 { 470 if (edgeorder[e] <= 1) continue; 471 472 int i0 = eledge[e][0]-1; 473 int i1 = eledge[e][1]-1; 474 475 if (edgeorient[e] == -1) Swap (i0, i1); // reverse orientation 476 477 double x = lambda(i1)-lambda(i0); 478 double y = 1-lambda(i0)-lambda(i1); 479 double fy = (1-y)*(1-y); 480 481 // double p3 = 0, p3x = 0, p3y = 0; 482 // double p2 = -1, p2x = 0, p2y = 0; 483 // double p1 = x, p1x = 1, p1y = 0; 484 485 double p3 = 0, p3x = 0, p3y = 0; 486 double p2 = -0.5, p2x = 0, p2y = 0; 487 double p1 = 0.5*x, p1x = 0.5, p1y = 0; 488 489 for (int j=2; j<= edgeorder[e]; j++) 490 { 491 p3=p2; p3x = p2x; p3y = p2y; 492 p2=p1; p2x = p1x; p2y = p1y; 493 double c1 = (2.0*j-3) / j; 494 double c2 = (j-3.0) / j; 495 496 p1 = c1 * x * p2 - c2 * fy * p3; 497 p1x = c1 * p2 + c1 * x * p2x - c2 * fy * p3x; 498 p1y = c1 * x * p2y - (c2 * 2 * (y-1) * p3 + c2 * fy * p3y); 499 500 eshape[index] = p1; 501 for (int k = 0; k < 2; k++) 502 edshape[index](k) = 503 p1x * ( dlambda(i1,k) - dlambda(i0,k) ) + 504 p1y * (-dlambda(i1,k) - dlambda(i0,k) ); 505 index++; 506 } 507 508 } 509 // (*testout) << "eshape = " << eshape << ", edshape = " << edshape << endl; 510 */ 511 } 512 513 CalcFaceShapes()514 void FETrig :: CalcFaceShapes () 515 { 516 int index = 0; 517 518 int i0 = elface[0][0]-1; 519 int i1 = elface[0][1]-1; 520 int i2 = elface[0][2]-1; 521 522 // sort lambda_i's by the corresponing vertex numbers 523 524 if (vertexnr[i1] < vertexnr[i0]) Swap (i0, i1); 525 if (vertexnr[i2] < vertexnr[i0]) Swap (i0, i2); 526 if (vertexnr[i2] < vertexnr[i1]) Swap (i1, i2); 527 528 double arg = lambda(i1) + lambda(i2); 529 530 if (fabs(arg) < EPSILON) // division by 0 531 { 532 for (int k = 0; k < nfaceshapes; k++) 533 { 534 fshape[index] = 0; 535 fdshape[index] = Vec<2>(0,0); 536 index++; 537 } 538 return; 539 } 540 541 b1.SetOrder (faceorder); 542 b2.SetOrder (faceorder); 543 544 b1.CalcFDf (lambda(i0)); 545 b2.CalcFDf (lambda(i2)/arg); 546 547 double decay = 1; 548 double ddecay; 549 550 double l0 = lambda(i0); 551 double l1 = lambda(i1); 552 double l2 = lambda(i2); 553 double dl0x = dlambda(i0,0); 554 double dl0y = dlambda(i0,1); 555 double dl1x = dlambda(i1,0); 556 double dl1y = dlambda(i1,1); 557 double dl2x = dlambda(i2,0); 558 double dl2y = dlambda(i2,1); 559 560 double dargx = dl1x + dl2x; 561 double dargy = dl1y + dl2y; 562 563 for (int n1 = 2; n1 < faceorder; n1++) 564 { 565 ddecay = (n1-1)*decay; 566 decay *= arg; 567 568 for (int n0 = 2; n0 < faceorder-n1+2; n0++) 569 { 570 fshape[index] = b1.GetF(n0) * b2.GetF(n1) * decay; 571 fdshape[index] = Vec<2> 572 (b1.GetDf(n0) * dl0x * b2.GetF(n1) * decay + 573 b1.GetF(n0) * b2.GetDf(n1) * (dl2x * arg - l2 * dargx)/(arg*arg) * decay + 574 b1.GetF(n0) * b2.GetF(n1) * ddecay * dargx, 575 b1.GetDf(n0) * dl0y * b2.GetF(n1) * decay + 576 b1.GetF(n0) * b2.GetDf(n1) * (dl2y * arg - l2 * dargy)/(arg*arg) * decay + 577 b1.GetF(n0) * b2.GetF(n1) * ddecay * dargy); 578 index++; 579 } 580 } 581 } 582 583 584 CalcFaceLaplaceShapes()585 void FETrig :: CalcFaceLaplaceShapes () 586 { 587 int index = 0; 588 589 int i0 = elface[0][0]-1; 590 int i1 = elface[0][1]-1; 591 int i2 = elface[0][2]-1; 592 593 if (vertexnr[i1] < vertexnr[i0]) Swap (i0, i1); 594 if (vertexnr[i2] < vertexnr[i0]) Swap (i0, i2); 595 if (vertexnr[i2] < vertexnr[i1]) Swap (i1, i2); 596 597 double arg = lambda(i1) + lambda(i2); 598 599 if (fabs(arg) < EPSILON) // division by 0 600 { 601 for (int k = 0; k < nfaceshapes; k++) 602 fddshape[k] = 0; 603 return; 604 } 605 606 b1.SetOrder (faceorder); 607 b2.SetOrder (faceorder); 608 609 b1.CalcFDf (lambda(i0)); 610 b1.CalcDDf (lambda(i0)); 611 b2.CalcFDf (lambda(i2)/arg); 612 b2.CalcDDf (lambda(i2)/arg); 613 614 double decay = 1; 615 double ddecay = 0; 616 double dddecay; 617 618 double l0 = lambda(i0); 619 double l1 = lambda(i1); 620 double l2 = lambda(i2); 621 double dl0x = dlambda(i0,0); 622 double dl0y = dlambda(i0,1); 623 double dl1x = dlambda(i1,0); 624 double dl1y = dlambda(i1,1); 625 double dl2x = dlambda(i2,0); 626 double dl2y = dlambda(i2,1); 627 628 double dargx = dl1x + dl2x; 629 double dargy = dl1y + dl2y; 630 631 for (int n1 = 2; n1 < faceorder; n1++) 632 { 633 dddecay = (n1-1)*ddecay; 634 ddecay = (n1-1)*decay; 635 decay *= arg; 636 637 for (int n0 = 2; n0 < faceorder-n1+2; n0++) 638 { 639 fddshape[index] = 640 641 // b1.GetDf(n0) * dl0x * b2.GetF(n1) * decay .... derived 642 643 b1.GetDDf(n0) * dl0x * dl0x * b2.GetF(n1) * decay + 644 b1.GetDf(n0) * dl0x * b2.GetDf(n1) * (dl2x * arg - l2 * dargx)/(arg*arg) * decay + 645 b1.GetDf(n0) * dl0x * b2.GetF(n1) * ddecay * dargx + 646 647 648 // b1.GetF(n0) * b2.GetDf(n1) * (dl2x * arg - l2 * dargx)/(arg*arg) * decay ... derived 649 650 b1.GetDf(n0) * dl0x * b2.GetDf(n1) * (dl2x * arg - l2 * dargx)/(arg*arg) * decay + 651 b1.GetF(n0) * b2.GetDDf(n1) * pow((dl2x * arg - l2 * dargx)/(arg*arg),2) * decay + 652 b1.GetF(n0) * b2.GetDf(n1) * (-2*dargx/arg) * (dl2x * arg - l2 * dargx)/(arg*arg) * decay + 653 b1.GetF(n0) * b2.GetDf(n1) * (dl2x * arg - l2 * dargx)/(arg*arg) * ddecay * dargx + 654 655 656 // b1.GetF(n0) * b2.GetF(n1) * ddecay * dargx ... derived 657 658 b1.GetDf(n0) * dl0x * b2.GetF(n1) * ddecay * dargx + 659 b1.GetF(n0) * b2.GetDf(n1) * (dl2x * arg - l2 * dargx)/(arg*arg) * ddecay * dargx + 660 b1.GetF(n0) * b2.GetF(n1) * dddecay * dargx * dargx + 661 662 663 // b1.GetDf(n0) * dl0y * b2.GetF(n1) * decay ... derived 664 665 b1.GetDDf(n0) * dl0y * dl0y * b2.GetF(n1) * decay + 666 b1.GetDf(n0) * dl0y * b2.GetDf(n1) * (dl2y * arg - l2 * dargy)/(arg*arg) * decay + 667 b1.GetDf(n0) * dl0y * b2.GetF(n1) * ddecay * dargy + 668 669 670 // b1.GetF(n0) * b2.GetDf(n1) * (dl2y * arg - l2 * dargy)/(arg*arg) * decay ... derived 671 672 b1.GetDf(n0) * dl0y * b2.GetDf(n1) * (dl2y * arg - l2 * dargy)/(arg*arg) * decay + 673 b1.GetF(n0) * b2.GetDDf(n1) * pow((dl2y * arg - l2 * dargy)/(arg*arg),2) * decay + 674 b1.GetF(n0) * b2.GetDf(n1) * (-2*dargy/arg) * (dl2y * arg - l2 * dargy)/(arg*arg) * decay + 675 b1.GetF(n0) * b2.GetDf(n1) * (dl2y * arg - l2 * dargy)/(arg*arg) * ddecay * dargy + 676 677 678 // b1.GetF(n0) * b2.GetF(n1) * ddecay * dargy ... derived 679 680 b1.GetDf(n0) * dl0y * b2.GetF(n1) * ddecay * dargy + 681 b1.GetF(n0) * b2.GetDf(n1) * (dl2y * arg - l2 * dargy)/(arg*arg) * ddecay * dargy + 682 b1.GetF(n0) * b2.GetF(n1) * dddecay * dargy * dargy; 683 684 index++; 685 } 686 } 687 } 688 689 690 691 // ---------------------------------------------------------------------------- 692 // FEQuad 693 // ---------------------------------------------------------------------------- 694 695 CalcVertexShapes()696 void FEQuad :: CalcVertexShapes () 697 { 698 vshape[0] = (1-xi(0)) * (1-xi(1)); 699 vshape[1] = ( xi(0)) * (1-xi(1)); 700 vshape[2] = ( xi(0)) * ( xi(1)); 701 vshape[3] = (1-xi(0)) * ( xi(1)); 702 703 vdshape[0] = Vec<2> ( -(1-xi(1)), -(1-xi(0)) ); 704 vdshape[1] = Vec<2> ( (1-xi(1)), -( xi(0)) ); 705 vdshape[2] = Vec<2> ( ( xi(1)), ( xi(0)) ); 706 vdshape[3] = Vec<2> ( -( xi(1)), (1-xi(0)) ); 707 } 708 709 CalcEdgeShapes()710 void FEQuad :: CalcEdgeShapes () 711 { 712 int index = 0; 713 714 double arg0[4] = { xi(0), 1-xi(0), 1-xi(1), xi(1) }; 715 double arg1[4] = { 1-xi(1), xi(1), 1-xi(0), xi(0) }; 716 double darg0[4] = { 1, -1, -1, 1 }; 717 double darg1[4] = { -1, 1, -1, 1 }; 718 719 for (int e = 0; e < nedges; e++) 720 { 721 b1.SetOrder (edgeorder[e]); 722 b1.CalcFDf (edgeorient[e] == 1 ? arg0[e] : 1-arg0[e]); 723 724 double decay = arg1[e]; 725 double ddecay; 726 727 for (int k = 2; k <= edgeorder[e]; k++, index++) 728 { 729 ddecay = k*decay; 730 decay *= arg1[e]; 731 // linear decay 732 eshape[index] = b1.GetF(k) * arg1[e]; 733 734 if (e < 2) 735 edshape[index] = Vec<2> 736 (darg0[e] * edgeorient[e] * b1.GetDf(k) * arg1[e], 737 b1.GetF(k) * darg1[e]); 738 else 739 edshape[index] = Vec<2> 740 (b1.GetF(k) * darg1[e], 741 darg0[e] * edgeorient[e] * b1.GetDf(k) * arg1[e]); 742 743 // exponential decay 744 /* eshape[index] = b1.GetF(k) * decay; 745 746 if (e < 2) 747 edshape[index] = Vec<2> 748 (darg0[e] * edgeorient[e] * b1.GetDf(k) * decay, 749 b1.GetF(k) * ddecay * darg1[e]); 750 else 751 edshape[index] = Vec<2> 752 (b1.GetF(k) * ddecay * darg1[e], 753 darg0[e] * edgeorient[e] * b1.GetDf(k) * decay); 754 */ 755 } 756 } 757 } 758 759 CalcFaceShapes()760 void FEQuad :: CalcFaceShapes () 761 { 762 int index = 0; 763 764 // find index of point with smallest number 765 766 int i0 = 0; 767 for (int i = 1; i < 4; i++) 768 if (vertexnr[elface[0][i]-1] < vertexnr[elface[0][i0]-1]) i0 = i; 769 770 double x; 771 double y; 772 double dxx; 773 double dxy; 774 double dyx; 775 double dyy; 776 777 switch (i0) 778 { 779 case 0: 780 x = xi(0); y = xi(1); 781 dxx = 1; dxy = 0; 782 dyx = 0; dyy = 1; 783 break; 784 case 1: 785 x = xi(1); y = 1-xi(0); 786 dxx = 0; dxy = 1; 787 dyx = -1; dyy = 0; 788 break; 789 case 2: 790 x = 1-xi(0); y = 1-xi(1); 791 dxx = -1; dxy = 0; 792 dyx = 0; dyy = -1; 793 break; 794 case 3: 795 x = 1-xi(1); y = xi(0); 796 dxx = 0; dxy =-1; 797 dyx = 1; dyy = 0; 798 break; 799 } 800 801 if (vertexnr[elface[0][(i0+1) % 4]-1] > vertexnr[elface[0][(i0+3) % 4]-1]) 802 { 803 Swap (x,y); 804 Swap (dxx, dyx); 805 Swap (dxy, dyy); 806 } 807 808 b1.SetOrder (faceorder); 809 b2.SetOrder (faceorder); 810 811 b1.CalcFDf (x); 812 b2.CalcFDf (y); 813 814 for (int n0 = 2; n0 <= faceorder; n0++) 815 for (int n1 = 2; n1 <= faceorder; n1++) 816 { 817 fshape[index] = b1.GetF(n0) * b2.GetF(n1); 818 fdshape[index] = Vec<2> (b1.GetDf(n0) * dxx * b2.GetF(n1) + b1.GetF(n0) * b2.GetDf(n1) * dyx, 819 b1.GetDf(n0) * dxy * b2.GetF(n1) + b1.GetF(n0) * b2.GetDf(n1) * dyy); 820 index++; 821 } 822 } 823 824 CalcFaceLaplaceShapes()825 void FEQuad :: CalcFaceLaplaceShapes () 826 { 827 int index = 0; 828 829 // find index of point with smallest number 830 831 int i0 = 0; 832 for (int i = 1; i < 4; i++) 833 if (vertexnr[elface[0][i]-1] < vertexnr[elface[0][i0]-1]) i0 = i; 834 835 double x; 836 double y; 837 double dxx; 838 double dxy; 839 double dyx; 840 double dyy; 841 842 switch (i0) 843 { 844 case 0: 845 x = xi(0); y = xi(1); 846 dxx = 1; dxy = 0; 847 dyx = 0; dyy = 1; 848 break; 849 case 1: 850 x = xi(1); y = 1-xi(0); 851 dxx = 0; dxy = 1; 852 dyx = -1; dyy = 0; 853 break; 854 case 2: 855 x = 1-xi(0); y = 1-xi(1); 856 dxx = -1; dxy = 0; 857 dyx = 0; dyy = -1; 858 break; 859 case 3: 860 x = 1-xi(1); y = xi(0); 861 dxx = 0; dxy =-1; 862 dyx = 1; dyy = 0; 863 break; 864 } 865 866 if (vertexnr[elface[0][(i0+1) % 4]-1] > vertexnr[elface[0][(i0+3) % 4]-1]) 867 { 868 Swap (x,y); 869 Swap (dxx, dyx); 870 Swap (dxy, dyy); 871 } 872 873 b1.SetOrder (faceorder); 874 b2.SetOrder (faceorder); 875 876 b1.CalcFDf (x); 877 b1.CalcDDf (x); 878 b2.CalcFDf (y); 879 b2.CalcDDf (y); 880 881 for (int n0 = 2; n0 <= faceorder; n0++) 882 for (int n1 = 2; n1 <= faceorder; n1++) 883 { 884 fddshape[index] = 885 b1.GetDDf(n0) * dxx * dxx * b2.GetF(n1) + 886 2* b1.GetDf(n0) * dxx * b2.GetDf(n1) * dyx + 887 b1.GetF(n0) * b2.GetDDf(n1) * dyx * dyx + 888 889 b1.GetDDf(n0) * dxy * dxy * b2.GetF(n1) + 890 2* b1.GetDf(n0) * dxy * b2.GetDf(n1) * dyy + 891 b1.GetF(n0) * b2.GetDDf(n1) * dyy * dyy; 892 893 index++; 894 } 895 } 896 897 898 899 // ---------------------------------------------------------------------------- 900 // FETet 901 // ---------------------------------------------------------------------------- 902 903 SetReferencePoint(Point<3> axi)904 void FETet :: SetReferencePoint (Point<3> axi) 905 { 906 BaseFiniteElement3D :: SetReferencePoint (axi); 907 908 lambda(0) = xi(0); 909 lambda(1) = xi(1); 910 lambda(2) = xi(2); 911 lambda(3) = 1-xi(0)-xi(1)-xi(2); 912 913 dlambda(0,0) = 1; dlambda(0,1) = 0; dlambda(0,2) = 0; 914 dlambda(1,0) = 0; dlambda(1,1) = 1; dlambda(1,2) = 0; 915 dlambda(2,0) = 0; dlambda(2,1) = 0; dlambda(2,2) = 1; 916 dlambda(3,0) = -1; dlambda(3,1) = -1; dlambda(3,2) = -1; 917 } 918 919 CalcVertexShapes()920 void FETet :: CalcVertexShapes () 921 { 922 for (int v = 0; v < nvertices; v++) 923 { 924 vshape[v] = lambda(v); 925 vdshape[v](0) = dlambda(v,0); 926 vdshape[v](1) = dlambda(v,1); 927 vdshape[v](2) = dlambda(v,2); 928 } 929 } 930 CalcVertexShapesOnly()931 void FETet :: CalcVertexShapesOnly () 932 { 933 for (int v = 0; v < nvertices; v++) 934 vshape[v] = lambda(v); 935 } 936 937 CalcEdgeShapes()938 void FETet :: CalcEdgeShapes () 939 { 940 int index = 0; 941 942 for (int e = 0; e < nedges; e++) 943 { 944 int i0 = eledge[e][0]-1; 945 int i1 = eledge[e][1]-1; 946 947 double arg = lambda(i0)+lambda(i1); 948 949 if (fabs(arg) < EPSILON) // division by 0 950 { 951 for (int k = 2; k <= edgeorder[e]; k++) 952 { 953 eshape[index] = 0; 954 edshape[index] = Vec<3>(0,0,0); 955 index++; 956 } 957 continue; 958 } 959 960 if (edgeorient[e] == -1) Swap (i0, i1); 961 962 double xi = lambda[i1]/arg; 963 964 b1.SetOrder (edgeorder[e]); 965 b1.CalcFDf (xi); 966 967 double decay = arg; 968 double ddecay; 969 970 double l0 = lambda(i0); 971 double dl0x = dlambda(i0,0); 972 double dl0y = dlambda(i0,1); 973 double dl0z = dlambda(i0,2); 974 975 double l1 = lambda(i1); 976 double dl1x = dlambda(i1,0); 977 double dl1y = dlambda(i1,1); 978 double dl1z = dlambda(i1,2); 979 980 double dargx = dl0x + dl1x; 981 double dargy = dl0y + dl1y; 982 double dargz = dl0z + dl1z; 983 984 for (int k = 2; k <= edgeorder[e]; k++) 985 { 986 ddecay = k*decay; 987 decay *= arg; 988 989 eshape[index] = b1.GetF (k) * decay; 990 edshape[index] = Vec<3> 991 (b1.GetDf(k) * (dl1x*arg - l1*dargx) * 992 decay / (arg * arg) + b1.GetF(k) * ddecay * dargx, 993 b1.GetDf(k) * (dl1y*arg - l1*dargy) * 994 decay / (arg * arg) + b1.GetF(k) * ddecay * dargy, 995 b1.GetDf(k) * (dl1z*arg - l1*dargz) * 996 decay / (arg * arg) + b1.GetF(k) * ddecay * dargz); 997 998 index++; 999 } 1000 } 1001 1002 1003 1004 /* 1005 int index = 0; 1006 for (int e = 0; e < nedges; e++) 1007 { 1008 if (edgeorder[e] <= 1) continue; 1009 1010 int i0 = eledge[e][0]-1; 1011 int i1 = eledge[e][1]-1; 1012 1013 if (edgeorient[e] == -1) Swap (i0, i1); // reverse orientation 1014 1015 double x = lambda(i1)-lambda(i0); 1016 double y = 1-lambda(i0)-lambda(i1); 1017 double fy = (1-y)*(1-y); 1018 1019 // double p3 = 0, p3x = 0, p3y = 0; 1020 // double p2 = -1, p2x = 0, p2y = 0; 1021 // double p1 = x, p1x = 1, p1y = 0; 1022 1023 double p3 = 0, p3x = 0, p3y = 0; 1024 double p2 = -0.5, p2x = 0, p2y = 0; 1025 double p1 = 0.5*x, p1x = 0.5, p1y = 0; 1026 1027 for (int j=2; j<= edgeorder[e]; j++) 1028 { 1029 p3=p2; p3x = p2x; p3y = p2y; 1030 p2=p1; p2x = p1x; p2y = p1y; 1031 double c1 = (2.0*j-3) / j; 1032 double c2 = (j-3.0) / j; 1033 1034 p1 = c1 * x * p2 - c2 * fy * p3; 1035 p1x = c1 * p2 + c1 * x * p2x - c2 * fy * p3x; 1036 p1y = c1 * x * p2y - (c2 * 2 * (y-1) * p3 + c2 * fy * p3y); 1037 1038 eshape[index] = p1; 1039 for (int k = 0; k < 3; k++) 1040 edshape[index](k) = 1041 p1x * ( dlambda(i1,k) - dlambda(i0,k) ) + 1042 p1y * (-dlambda(i1,k) - dlambda(i0,k) ); 1043 index++; 1044 } 1045 1046 } 1047 1048 */ 1049 1050 1051 1052 1053 } 1054 1055 1056 1057 CalcEdgeShapesOnly()1058 void FETet :: CalcEdgeShapesOnly () 1059 { 1060 int index = 0; 1061 1062 for (int e = 0; e < nedges; e++) 1063 { 1064 int i0 = eledge[e][0]-1; 1065 int i1 = eledge[e][1]-1; 1066 1067 if (edgeorient[e] == -1) Swap (i0, i1); 1068 1069 double arg = lambda(i0)+lambda(i1); 1070 double xi = lambda[i1]; 1071 1072 /* 1073 b1.SetOrder (edgeorder[e]); 1074 b1.CalcFScaled (xi, arg); 1075 1076 for (int k = 2; k <= edgeorder[e]; k++, index++) 1077 eshape[index] = b1.GetF (k); 1078 */ 1079 b1.CalcFScaled (edgeorder[e], xi, arg, &eshape[index]); 1080 index += edgeorder[e]-1; 1081 } 1082 } 1083 1084 1085 1086 1087 1088 CalcFaceShapes()1089 void FETet :: CalcFaceShapes () 1090 { 1091 int index = 0; 1092 1093 for (int f = 0; f < nfaces; f++) 1094 { 1095 if (faceorder[f] <= 2) continue; 1096 1097 int i0 = elface[f][0]-1; 1098 int i1 = elface[f][1]-1; 1099 int i2 = elface[f][2]-1; 1100 1101 if (vertexnr[i1] < vertexnr[i0]) Swap (i0, i1); 1102 if (vertexnr[i2] < vertexnr[i0]) Swap (i0, i2); 1103 if (vertexnr[i2] < vertexnr[i1]) Swap (i1, i2); 1104 1105 double arg = lambda(i1) + lambda(i2); 1106 double arg2 = lambda(i0) + lambda(i1) + lambda(i2); 1107 1108 if (fabs(arg) < EPSILON || fabs(arg2) < EPSILON) // division by 0 1109 { 1110 for (int k = 0; k < nfaceshapes[f]; k++) 1111 { 1112 fshape[index] = 0; 1113 fdshape[index] = Vec<3>(0,0,0); 1114 index++; 1115 } 1116 continue; 1117 } 1118 1119 b1.SetOrder (faceorder[f]); 1120 b2.SetOrder (faceorder[f]); 1121 1122 b1.CalcFDf (lambda(i0)/arg2); 1123 b2.CalcFDf (lambda(i2)/arg); 1124 1125 double decay = 1; 1126 double ddecay; 1127 1128 double l0 = lambda(i0); 1129 double l1 = lambda(i1); 1130 double l2 = lambda(i2); 1131 double dl0x = dlambda(i0,0); 1132 double dl0y = dlambda(i0,1); 1133 double dl0z = dlambda(i0,2); 1134 double dl1x = dlambda(i1,0); 1135 double dl1y = dlambda(i1,1); 1136 double dl1z = dlambda(i1,2); 1137 double dl2x = dlambda(i2,0); 1138 double dl2y = dlambda(i2,1); 1139 double dl2z = dlambda(i2,2); 1140 1141 double dargx = dl1x + dl2x; 1142 double dargy = dl1y + dl2y; 1143 double dargz = dl1z + dl2z; 1144 double darg2x = dl0x + dl1x + dl2x; 1145 double darg2y = dl0y + dl1y + dl2y; 1146 double darg2z = dl0z + dl1z + dl2z; 1147 1148 for (int n1 = 2; n1 < faceorder[f]; n1++) 1149 { 1150 ddecay = (n1-1)*decay; 1151 decay *= arg; 1152 1153 double decay2 = arg2; 1154 double ddecay2; 1155 1156 for (int n0 = 2; n0 < faceorder[f]-n1+2; n0++) 1157 { 1158 ddecay2 = n0*decay2; 1159 decay2 *= arg2; 1160 1161 fshape[index] = b1.GetF(n0) * b2.GetF(n1) * decay * decay2; 1162 fdshape[index] = Vec<3> 1163 (b1.GetDf(n0) * (dl0x * arg2 - l0 * darg2x)/(arg2*arg2) * b2.GetF(n1) * decay * decay2 + 1164 b1.GetF(n0) * b2.GetDf(n1) * (dl2x * arg - l2 * dargx)/(arg*arg) * decay * decay2 + 1165 b1.GetF(n0) * b2.GetF(n1) * ddecay * dargx * decay2 + 1166 b1.GetF(n0) * b2.GetF(n1) * decay * ddecay2 * darg2x, 1167 1168 b1.GetDf(n0) * (dl0y * arg2 - l0 * darg2y)/(arg2*arg2) * b2.GetF(n1) * decay * decay2 + 1169 b1.GetF(n0) * b2.GetDf(n1) * (dl2y * arg - l2 * dargy)/(arg*arg) * decay * decay2 + 1170 b1.GetF(n0) * b2.GetF(n1) * ddecay * dargy * decay2 + 1171 b1.GetF(n0) * b2.GetF(n1) * decay * ddecay2 * darg2y, 1172 1173 b1.GetDf(n0) * (dl0z * arg2 - l0 * darg2z)/(arg2*arg2) * b2.GetF(n1) * decay * decay2 + 1174 b1.GetF(n0) * b2.GetDf(n1) * (dl2z * arg - l2 * dargz)/(arg*arg) * decay * decay2 + 1175 b1.GetF(n0) * b2.GetF(n1) * ddecay * dargz * decay2 + 1176 b1.GetF(n0) * b2.GetF(n1) * decay * ddecay2 * darg2z); 1177 1178 index++; 1179 } 1180 } 1181 } 1182 } 1183 1184 1185 1186 1187 // ---------------------------------------------------------------------------- 1188 // FEPrism 1189 // ---------------------------------------------------------------------------- 1190 1191 SetReferencePoint(Point<3> axi)1192 void FEPrism :: SetReferencePoint (Point<3> axi) 1193 { 1194 BaseFiniteElement3D :: SetReferencePoint (axi); 1195 1196 lambda(0) = xi(0); 1197 lambda(1) = xi(1); 1198 lambda(2) = 1-xi(0)-xi(1); 1199 lambda(3) = xi(2); 1200 1201 dlambda(0,0) = 1; dlambda(0,1) = 0; dlambda(0,2) = 0; 1202 dlambda(1,0) = 0; dlambda(1,1) = 1; dlambda(1,2) = 0; 1203 dlambda(2,0) = -1; dlambda(2,1) = -1; dlambda(2,2) = 0; 1204 dlambda(3,0) = 0; dlambda(3,1) = 0; dlambda(3,2) = 1; 1205 } 1206 1207 CalcVertexShapes()1208 void FEPrism :: CalcVertexShapes () 1209 { 1210 double zcomp = 1-lambda(3); 1211 1212 for (int v = 0; v < nvertices; v++) 1213 { 1214 if (v == 3) zcomp = 1-zcomp; 1215 1216 vshape[v] = lambda(v % 3) * zcomp; 1217 vdshape[v](0) = dlambda(v % 3,0) * zcomp; 1218 vdshape[v](1) = dlambda(v % 3,1) * zcomp; 1219 vdshape[v](2) = lambda(v % 3) * (-dlambda(3,2)); 1220 1221 if (v >= 3) vdshape[v](2) *= -1; 1222 } 1223 } 1224 1225 CalcEdgeShapes()1226 void FEPrism :: CalcEdgeShapes () 1227 { 1228 int index = 0; 1229 int e; 1230 // triangle edge shapes 1231 1232 for (e = 0; e < 6; e++) 1233 { 1234 int i0 = (eledge[e][0]-1) % 3; 1235 int i1 = (eledge[e][1]-1) % 3; 1236 1237 double arg = lambda[i0]+lambda[i1]; 1238 1239 if (fabs(arg) < EPSILON) // division by 0 1240 { 1241 for (int k = 2; k <= edgeorder[e]; k++) 1242 { 1243 eshape[index] = 0; 1244 edshape[index] = Vec<3>(0,0,0); 1245 index++; 1246 } 1247 continue; 1248 } 1249 1250 if (edgeorient[e] == -1) Swap (i0, i1); 1251 1252 double xi = lambda[i1]/arg; 1253 1254 b1.SetOrder (edgeorder[e]); 1255 b1.CalcFDf (xi); 1256 1257 double decay = arg; 1258 double ddecay; 1259 1260 double zarg = e < 3 ? (1-lambda(3)) : lambda(3); 1261 double zcomp = zarg; 1262 double dzcomp; 1263 1264 double l0 = lambda(i0); 1265 double dl0x = dlambda(i0,0); 1266 double dl0y = dlambda(i0,1); 1267 1268 double l1 = lambda(i1); 1269 double dl1x = dlambda(i1,0); 1270 double dl1y = dlambda(i1,1); 1271 1272 double dargx = dl0x + dl1x; 1273 double dargy = dl0y + dl1y; 1274 1275 1276 /* 1277 1278 for (int k = 2; k <= edgeorder[e]; k++) 1279 { 1280 ddecay = k*decay; 1281 decay *= arg; 1282 1283 dzcomp = k*zcomp; 1284 zcomp *= zarg; 1285 1286 eshape[index] = b1.GetF (k) * decay * zcomp; 1287 edshape[index] = Vec<3> 1288 ((b1.GetDf(k) * (dl1x*arg - l1*dargx) * 1289 decay / (arg * arg) + b1.GetF(k) * ddecay * dargx) * zcomp, 1290 (b1.GetDf(k) * (dl1y*arg - l1*dargy) * 1291 decay / (arg * arg) + b1.GetF(k) * ddecay * dargy) * zcomp, 1292 b1.GetF(k) * decay * dzcomp * (e < 3 ? -1 : 1)); 1293 index++; 1294 } 1295 */ 1296 1297 // JS linear in z-direction (for thin plates !) 1298 for (int k = 2; k <= edgeorder[e]; k++) 1299 { 1300 ddecay = k*decay; 1301 decay *= arg; 1302 1303 // dzcomp = k*zcomp; 1304 // zcomp *= zarg; 1305 dzcomp = 1; 1306 zcomp = zarg; 1307 1308 eshape[index] = b1.GetF (k) * decay * zcomp; 1309 edshape[index] = Vec<3> 1310 ((b1.GetDf(k) * (dl1x*arg - l1*dargx) * 1311 decay / (arg * arg) + b1.GetF(k) * ddecay * dargx) * zcomp, 1312 (b1.GetDf(k) * (dl1y*arg - l1*dargy) * 1313 decay / (arg * arg) + b1.GetF(k) * ddecay * dargy) * zcomp, 1314 b1.GetF(k) * decay * dzcomp * (e < 3 ? -1 : 1)); 1315 index++; 1316 } 1317 1318 1319 1320 } 1321 1322 // rectangle edge shapes 1323 1324 for (e = 6; e < nedges; e++) 1325 { 1326 int i0 = eledge[e][0]-1; 1327 1328 double arg = lambda[i0]; 1329 1330 if (fabs(arg) < EPSILON) // division by 0 1331 { 1332 for (int k = 2; k <= edgeorder[e]; k++) 1333 { 1334 eshape[index] = 0.; 1335 edshape[index] = Vec<3>(0.,0.,0.); 1336 index++; 1337 } 1338 continue; 1339 } 1340 1341 double xi = lambda[3]; 1342 1343 if (edgeorient[e] == -1) xi = 1-xi; 1344 1345 b1.SetOrder (edgeorder[e]); 1346 b1.CalcFDf (xi); 1347 1348 double decay = arg; 1349 double ddecay; 1350 1351 double l0 = lambda(i0); 1352 double l0x = dlambda(i0,0); 1353 double l0y = dlambda(i0,1); 1354 1355 for (int k = 2; k <= edgeorder[e]; k++) 1356 { 1357 ddecay = k*decay; 1358 decay *= arg; 1359 1360 eshape[index] = b1.GetF (k) * decay; 1361 edshape[index] = Vec<3> 1362 (b1.GetF(k) * ddecay * l0x, 1363 b1.GetF(k) * ddecay * l0y, 1364 b1.GetDf(k) * edgeorient[e] * decay); 1365 index++; 1366 } 1367 } 1368 } 1369 1370 CalcFaceShapes()1371 void FEPrism :: CalcFaceShapes () 1372 { 1373 int index = 0; 1374 int f; 1375 1376 // triangle face parts 1377 1378 for (f = 0; f < 2; f++) 1379 { 1380 int i0 = elface[f][0]-1; 1381 int i1 = elface[f][1]-1; 1382 int i2 = elface[f][2]-1; 1383 1384 if (vertexnr[i1] < vertexnr[i0]) Swap (i0, i1); 1385 if (vertexnr[i2] < vertexnr[i0]) Swap (i0, i2); 1386 if (vertexnr[i2] < vertexnr[i1]) Swap (i1, i2); 1387 1388 i0 = i0 % 3; 1389 i1 = i1 % 3; 1390 i2 = i2 % 3; 1391 1392 double arg = lambda(i1) + lambda(i2); 1393 1394 if (fabs(arg) < EPSILON) // division by 0 1395 { 1396 for (int k = 0; k < nfaceshapes[f]; k++) 1397 { 1398 fshape[index] = 0; 1399 fdshape[index] = Vec<3>(0,0,0); 1400 index++; 1401 } 1402 continue; 1403 } 1404 1405 b1.SetOrder (faceorder[f]); 1406 b2.SetOrder (faceorder[f]); 1407 1408 b1.CalcFDf (lambda(i0)); 1409 b2.CalcFDf (lambda(i2)/arg); 1410 1411 double decay = 1; 1412 double ddecay; 1413 1414 double l0 = lambda(i0); 1415 double l1 = lambda(i1); 1416 double l2 = lambda(i2); 1417 double dl0x = dlambda(i0,0); 1418 double dl0y = dlambda(i0,1); 1419 double dl0z = dlambda(i0,2); 1420 double dl1x = dlambda(i1,0); 1421 double dl1y = dlambda(i1,1); 1422 double dl1z = dlambda(i1,2); 1423 double dl2x = dlambda(i2,0); 1424 double dl2y = dlambda(i2,1); 1425 double dl2z = dlambda(i2,2); 1426 1427 double dargx = dl1x + dl2x; 1428 double dargy = dl1y + dl2y; 1429 double dargz = dl1z + dl2z; 1430 1431 double arg2 = f == 0 ? 1-xi(2) : xi(2); 1432 double darg2z = f == 0 ? -1 : 1; 1433 1434 for (int n1 = 2; n1 < faceorder[f]; n1++) 1435 { 1436 ddecay = (n1-1)*decay; 1437 decay *= arg; 1438 1439 double decay2 = arg2; 1440 double ddecay2; 1441 1442 for (int n0 = 2; n0 < faceorder[f]-n1+2; n0++) 1443 { 1444 ddecay2 = n0*decay2; 1445 decay2 *= arg2; 1446 1447 fshape[index] = b1.GetF(n0) * b2.GetF(n1) * decay * decay2; 1448 fdshape[index] = Vec<3> 1449 (b1.GetDf(n0) * dl0x * b2.GetF(n1) * decay * decay2 + 1450 b1.GetF(n0) * b2.GetDf(n1) * (dl2x * arg - l2 * dargx)/(arg*arg) * decay * decay2 + 1451 b1.GetF(n0) * b2.GetF(n1) * ddecay * dargx * decay2, 1452 1453 b1.GetDf(n0) * dl0y * b2.GetF(n1) * decay * decay2 + 1454 b1.GetF(n0) * b2.GetDf(n1) * (dl2y * arg - l2 * dargy)/(arg*arg) * decay * decay2 + 1455 b1.GetF(n0) * b2.GetF(n1) * ddecay * dargy * decay2, 1456 1457 b1.GetDf(n0) * dl0z * b2.GetF(n1) * decay * decay2 + 1458 b1.GetF(n0) * b2.GetDf(n1) * (dl2z * arg - l2 * dargz)/(arg*arg) * decay * decay2 + 1459 b1.GetF(n0) * b2.GetF(n1) * ddecay * dargz * decay2 + 1460 b1.GetF(n0) * b2.GetF(n1) * decay * ddecay2 * darg2z); 1461 1462 index++; 1463 } 1464 } 1465 } 1466 1467 1468 // quad parts 1469 1470 for (f = 2; f < nfaces; f++) 1471 { 1472 // find index of point with smallest number 1473 1474 int i, i0 = 0; 1475 for (i = 1; i < 4; i++) 1476 if (vertexnr[elface[f][i]-1] < vertexnr[elface[f][i0]-1]) i0 = i; 1477 1478 double arg = 0; 1479 double dargx = 0; 1480 double dargy = 0; 1481 double dargz = 0; 1482 for (i = 0; i < 4; i++) 1483 { 1484 arg += lambda((elface[f][i]-1) % 3)/2.0; 1485 dargx += dlambda((elface[f][i]-1) % 3,0)/2.0; 1486 dargy += dlambda((elface[f][i]-1) % 3,1)/2.0; 1487 dargz += dlambda((elface[f][i]-1) % 3,2)/2.0; 1488 } 1489 1490 if (fabs(arg) < EPSILON) // division by 0 1491 { 1492 for (int k = 0; k < nfaceshapes[f]; k++) 1493 { 1494 fshape[index] = 0; 1495 fdshape[index] = Vec<3>(0,0,0); 1496 index++; 1497 } 1498 continue; 1499 } 1500 1501 int i1 = (i0+3) % 4; 1502 int j = (elface[f][i0]-1) % 3; 1503 1504 double lam = lambda(j)/arg; 1505 double dlamx = (dlambda(j,0)*arg-lambda(j)*dargx)/(arg*arg); 1506 double dlamy = (dlambda(j,1)*arg-lambda(j)*dargy)/(arg*arg); 1507 double dlamz = (dlambda(j,2)*arg-lambda(j)*dargz)/(arg*arg); 1508 1509 double x; 1510 double z; 1511 double dxx; 1512 double dxy; 1513 double dxz; 1514 double dzx; 1515 double dzy; 1516 double dzz; 1517 1518 int ratvar; 1519 /* 1520 switch (i0) 1521 { 1522 case 0: 1523 x = xi(2); z = lam; 1524 1525 dxx = 0; dxy = 0; dxz = 1; 1526 dzx = dlamx; dzy = dlamy; dzz = dlamz; 1527 ratvar = 1; 1528 break; 1529 case 1: 1530 x = 1-lam; z = xi(2); 1531 dxx = -dlamx; dxy = -dlamy; dxz = -dlamz; 1532 dzx = 0; dzy = 0; dzz = 1; 1533 ratvar = 0; 1534 break; 1535 case 2: 1536 x = 1-xi(2); z = 1-lam; 1537 dxx = 0; dxy = 0; dxz = -1; 1538 dzx = -dlamx; dzy = -dlamy; dzz = -dlamz; 1539 ratvar = 1; 1540 break; 1541 case 3: 1542 x = lam; z = 1-xi(2); 1543 dxx = dlamx; dxy = dlamy; dxz = dlamz; 1544 dzx = 0; dzy = 0; dzz = -1; 1545 ratvar = 0; 1546 break; 1547 } 1548 */ 1549 1550 ratvar = 0; 1551 x = 1-lam; 1552 dxx = -dlamx; dxy = -dlamy; dxz = -dlamz; 1553 if (i0 == 0 || i0 == 1) 1554 { 1555 z = xi(2); 1556 dzx = 0; dzy = 0; dzz = 1; 1557 } 1558 else 1559 { 1560 z = 1-xi(2); 1561 dzx = 0; dzy = 0; dzz = -1; 1562 } 1563 1564 int ix = i0 ^ 1; 1565 int iz = 3-i0; 1566 1567 if (vertexnr[elface[f][ix]-1] > vertexnr[elface[f][iz]-1]) 1568 { 1569 Swap (x,z); 1570 Swap (dxx, dzx); 1571 Swap (dxy, dzy); 1572 Swap (dxz, dzz); 1573 ratvar = 1-ratvar; 1574 } 1575 1576 b1.SetOrder (faceorder[f]); 1577 b2.SetOrder (faceorder[f]); 1578 1579 b1.CalcFDf (x); 1580 b2.CalcFDf (z); 1581 1582 double decay = arg; 1583 double ddecay; 1584 1585 for (int n0 = 2; n0 <= faceorder[f]; n0++) 1586 { 1587 ddecay = n0*decay; 1588 decay *= arg; 1589 1590 if (ratvar == 1) decay = arg; 1591 1592 for (int n1 = 2; n1 <= faceorder[f]; n1++) 1593 { 1594 if (ratvar == 1) 1595 { 1596 ddecay = n1*decay; 1597 decay *= arg; 1598 } 1599 1600 fshape[index] = b1.GetF(n0) * b2.GetF(n1) * decay; 1601 fdshape[index] = Vec<3> 1602 (b1.GetDf(n0) * dxx * b2.GetF(n1) * decay + 1603 b1.GetF(n0) * b2.GetDf(n1) * dzx * decay + 1604 b1.GetF(n0) * b2.GetF(n1) * ddecay * dargx, 1605 1606 b1.GetDf(n0) * dxy * b2.GetF(n1) * decay + 1607 b1.GetF(n0) * b2.GetDf(n1) * dzy * decay + 1608 b1.GetF(n0) * b2.GetF(n1) * ddecay * dargy, 1609 1610 b1.GetDf(n0) * dxz * b2.GetF(n1) * decay + 1611 b1.GetF(n0) * b2.GetDf(n1) * dzz * decay + 1612 b1.GetF(n0) * b2.GetF(n1) * ddecay * dargz); 1613 1614 index++; 1615 } 1616 } 1617 } 1618 } 1619 1620 1621 1622 // ---------------------------------------------------------------------------- 1623 // FEPyramid 1624 // ---------------------------------------------------------------------------- 1625 1626 SetReferencePoint(Point<3> axi)1627 void FEPyramid :: SetReferencePoint (Point<3> axi) 1628 { 1629 BaseFiniteElement3D :: SetReferencePoint (axi); 1630 } 1631 1632 CalcVertexShapes()1633 void FEPyramid :: CalcVertexShapes () 1634 { 1635 double x = xi(0); 1636 double y = xi(1); 1637 double z = xi(2); 1638 1639 if (z == 1.) z = 1-1e-10; 1640 vshape[0] = (1-z-x)*(1-z-y) / (1-z); 1641 vshape[1] = x*(1-z-y) / (1-z); 1642 vshape[2] = x*y / (1-z); 1643 vshape[3] = (1-z-x)*y / (1-z); 1644 vshape[4] = z; 1645 1646 double z1 = 1-z; 1647 double z2 = z1*z1; 1648 1649 vdshape[0] = Vec<3>( -(z1-y)/z1, -(z1-x)/z1, ((x+y+2*z-2)*z1+(z1-y)*(z1-x))/z2 ); 1650 vdshape[1] = Vec<3>( (z1-y)/z1, -x/z1, (-x*z1+x*(z1-y))/z2 ); 1651 vdshape[2] = Vec<3>( y/z1, x/z1, x*y/z2 ); 1652 vdshape[3] = Vec<3>( -y/z1, (z1-x)/z1, (-y*z1+y*(z1-x))/z2 ); 1653 vdshape[4] = Vec<3>( 0, 0, 1 ); 1654 } 1655 1656 CalcEdgeShapes()1657 void FEPyramid :: CalcEdgeShapes () 1658 { 1659 int index = 0; 1660 1661 for (int e = 0; e < GetNEdges(); e++) 1662 { 1663 for (int k = 2; k <= edgeorder[e]; k++) 1664 { 1665 eshape[index] = 0; 1666 edshape[index] = Vec<3>(0,0,0); 1667 index++; 1668 } 1669 } 1670 } 1671 1672 CalcFaceShapes()1673 void FEPyramid :: CalcFaceShapes () 1674 { 1675 int index = 0; 1676 1677 for (int f = 0; f < GetNFaces(); f++) 1678 { 1679 for (int k = 0; k < nfaceshapes[f]; k++) 1680 { 1681 fshape[index] = 0; 1682 fdshape[index] = Vec<3>(0,0,0); 1683 index++; 1684 } 1685 } 1686 } 1687 1688 1689 1690 1691 1692 // ---------------------------------------------------------------------------- 1693 // FEHex 1694 // ---------------------------------------------------------------------------- 1695 1696 SetReferencePoint(Point<3> axi)1697 void FEHex :: SetReferencePoint (Point<3> axi) 1698 { 1699 BaseFiniteElement3D :: SetReferencePoint (axi); 1700 } 1701 1702 CalcVertexShapes()1703 void FEHex :: CalcVertexShapes () 1704 { 1705 double x = xi(0); 1706 double y = xi(1); 1707 double z = xi(2); 1708 1709 vshape[0] = (1-x)*(1-y)*(1-z); 1710 vshape[1] = x *(1-y)*(1-z); 1711 vshape[2] = x * y *(1-z); 1712 vshape[3] = (1-x)* y *(1-z); 1713 vshape[4] = (1-x)*(1-y)* z; 1714 vshape[5] = x *(1-y)* z; 1715 vshape[6] = x * y * z; 1716 vshape[7] = (1-x)* y * z; 1717 1718 vdshape[0] = Vec<3>(-(1-y)*(1-z), -(1-x)*(1-z), -(1-x)*(1-y)); 1719 vdshape[1] = Vec<3>( (1-y)*(1-z), -x *(1-z), -x *(1-y)); 1720 vdshape[2] = Vec<3>( y *(1-z), x *(1-z), -x * y); 1721 vdshape[3] = Vec<3>( -y *(1-z), (1-x)*(1-z), -(1-x)*y); 1722 vdshape[4] = Vec<3>(-(1-y)* z, -(1-x)* z, (1-x)*(1-y)); 1723 vdshape[5] = Vec<3>( (1-y)* z, -x * z, x *(1-y)); 1724 vdshape[6] = Vec<3>( y * z, x * z, x * y); 1725 vdshape[7] = Vec<3>( -y * z, (1-x)* z, (1-x)*y); 1726 1727 } 1728 1729 CalcEdgeShapes()1730 void FEHex :: CalcEdgeShapes () 1731 { 1732 int index = 0; 1733 1734 for (int e = 0; e < GetNEdges(); e++) 1735 { 1736 for (int k = 2; k <= edgeorder[e]; k++) 1737 { 1738 eshape[index] = 0; 1739 edshape[index] = Vec<3>(0,0,0); 1740 index++; 1741 } 1742 } 1743 } 1744 1745 CalcFaceShapes()1746 void FEHex :: CalcFaceShapes () 1747 { 1748 int index = 0; 1749 1750 for (int f = 0; f < GetNFaces(); f++) 1751 { 1752 for (int k = 0; k < nfaceshapes[f]; k++) 1753 { 1754 fshape[index] = 0; 1755 fdshape[index] = Vec<3>(0,0,0); 1756 index++; 1757 } 1758 } 1759 } 1760 1761 1762 1763 1764 1765 1766 1767 1768 IsSurfaceElementCurved(int elnr) const1769 int CurvedElements :: IsSurfaceElementCurved (int elnr) const 1770 { 1771 if (mesh.coarsemesh) 1772 { 1773 const HPRefElement & hpref_el = 1774 (*mesh.hpelements) [ mesh[(SurfaceElementIndex) elnr].hp_elnr]; 1775 1776 return mesh.coarsemesh->GetCurvedElements().IsSurfaceElementCurved (hpref_el.coarse_elnr); 1777 } 1778 1779 1780 1781 1782 Element2d elem = mesh[(SurfaceElementIndex) elnr]; 1783 1784 switch (elem.GetType()) 1785 { 1786 case TRIG: 1787 { 1788 FETrig fe2d(*this); 1789 fe2d.SetElementNumber (elnr+1); 1790 return (fe2d.IsCurved()); 1791 break; 1792 } 1793 1794 case QUAD: 1795 { 1796 FEQuad fe2d(*this); 1797 fe2d.SetElementNumber (elnr+1); 1798 return (fe2d.IsCurved()); 1799 break; 1800 } 1801 1802 } 1803 return 0; 1804 } 1805 1806 1807 IsElementCurved(int elnr) const1808 int CurvedElements :: IsElementCurved (int elnr) const 1809 { 1810 if (mesh.coarsemesh) 1811 { 1812 const HPRefElement & hpref_el = 1813 (*mesh.hpelements) [ mesh[(ElementIndex) elnr].hp_elnr]; 1814 1815 return mesh.coarsemesh->GetCurvedElements().IsElementCurved (hpref_el.coarse_elnr); 1816 } 1817 1818 1819 1820 Element elem = mesh[(ElementIndex) elnr]; 1821 1822 switch (elem.GetType()) 1823 { 1824 case TET: 1825 { 1826 FETet fe3d(*this); 1827 fe3d.SetElementNumber (elnr+1); 1828 return (fe3d.IsCurved()); 1829 break; 1830 } 1831 1832 case PRISM: 1833 { 1834 FEPrism fe3d(*this); 1835 fe3d.SetElementNumber (elnr+1); 1836 return (fe3d.IsCurved()); 1837 break; 1838 } 1839 1840 case PYRAMID: 1841 { 1842 FEPyramid fe3d(*this); 1843 fe3d.SetElementNumber (elnr+1); 1844 return (fe3d.IsCurved()); 1845 break; 1846 } 1847 1848 case HEX: 1849 { 1850 FEHex fe3d(*this); 1851 fe3d.SetElementNumber (elnr+1); 1852 return (fe3d.IsCurved()); 1853 break; 1854 } 1855 1856 } 1857 1858 return 0; 1859 } 1860 1861 CalcSegmentTransformation(double xi,int segnr,Point<3> * x,Vec<3> * dxdxi)1862 void CurvedElements :: CalcSegmentTransformation (double xi, int segnr, 1863 Point<3> * x, Vec<3> * dxdxi) 1864 { 1865 if (mesh.coarsemesh) 1866 { 1867 const HPRefElement & hpref_el = 1868 (*mesh.hpelements) [ mesh[(SegmentIndex) segnr].hp_elnr]; 1869 1870 // xi umrechnen 1871 double lami[2]; 1872 lami[0] = xi; 1873 lami[1] = 1-xi; 1874 1875 double coarse_xi = 0; 1876 for (int i = 0; i < 2; i++) 1877 coarse_xi += hpref_el.param[i][0] * lami[i]; 1878 1879 mesh.coarsemesh->GetCurvedElements().CalcSegmentTransformation (coarse_xi, hpref_el.coarse_elnr, x, dxdxi); 1880 return; 1881 } 1882 1883 1884 // xi = 1-xi; // or, this way ????? JS 1885 FESegm segm (*this); 1886 segm.SetElementNumber (segnr+1); 1887 segm.SetReferencePoint (Point<1>(xi)); 1888 1889 // segm.CalcVertexShapes (x != NULL, dxdxi != NULL); 1890 segm.CalcVertexShapes (); 1891 1892 if (x) 1893 { 1894 (*x) = Point<3>(0,0,0); 1895 for (int v = 0; v < 2; v++) 1896 (*x) = (*x) + segm.GetVertexShape(v) * mesh.Point(segm.GetVertexNr(v)); 1897 } 1898 1899 if (dxdxi) 1900 { 1901 (*dxdxi) = Vec<3>(0,0,0); 1902 for (int v = 0; v < 2; v++) 1903 (*dxdxi) = (*dxdxi) + segm.GetVertexDShape(v) * mesh.Point(segm.GetVertexNr(v)); 1904 } 1905 1906 if (segm.GetEdgeOrder() > 1) 1907 { 1908 // segm.CalcEdgeShapes (x != NULL, dxdxi != NULL); 1909 segm.CalcEdgeShapes (); 1910 1911 if (x) 1912 { 1913 int gindex = edgecoeffsindex[segm.GetEdgeNr()-1]; 1914 for (int k = 2; k <= segm.GetEdgeOrder(); k++, gindex++) 1915 (*x) = (*x) + segm.GetEdgeShape(k-2) * edgecoeffs[gindex]; 1916 } 1917 1918 if (dxdxi) 1919 { 1920 int gindex = edgecoeffsindex[segm.GetEdgeNr()-1]; 1921 for (int k = 2; k <= segm.GetEdgeOrder(); k++, gindex++) 1922 (*dxdxi) = (*dxdxi) + segm.GetEdgeDShape(k-2) * edgecoeffs[gindex]; 1923 } 1924 } 1925 } 1926 1927 1928 CalcMultiPointSegmentTransformation(ARRAY<double> * xi,int segnr,ARRAY<Point<3>> * x,ARRAY<Vec<3>> * dxdxi)1929 void CurvedElements :: CalcMultiPointSegmentTransformation (ARRAY<double> * xi, int segnr, 1930 ARRAY<Point<3> > * x, 1931 ARRAY<Vec<3> > * dxdxi) 1932 { 1933 int size = xi->Size(); 1934 1935 if (mesh.coarsemesh) 1936 { 1937 const HPRefElement & hpref_el = 1938 (*mesh.hpelements) [ mesh[(SegmentIndex) segnr].hp_elnr]; 1939 1940 // xi umrechnen 1941 ARRAY< Point<2> > lami; 1942 lami.SetSize (size); 1943 for (int i = 0; i<size; i++) 1944 { 1945 lami[i](0) = (*xi)[i]; 1946 lami[i](1) = 1-(*xi)[i]; 1947 } 1948 1949 ARRAY<double> coarse_xi; 1950 coarse_xi.SetSize (size); 1951 coarse_xi = 0; 1952 1953 1954 for (int j = 0; j < 2; j++) 1955 for (int i = 0; i<size; i++) 1956 coarse_xi[i] += hpref_el.param[j][0] * lami[i](j); 1957 1958 mesh.coarsemesh->GetCurvedElements().CalcMultiPointSegmentTransformation 1959 (&coarse_xi, hpref_el.coarse_elnr, x, dxdxi); 1960 return; 1961 } 1962 1963 1964 1965 FESegm segm (*this); 1966 segm.SetElementNumber (segnr+1); 1967 x->SetSize (size); 1968 dxdxi->SetSize (size); 1969 1970 for (int i = 0; i < size; i++) 1971 { 1972 segm.SetReferencePoint (Point<1>((*xi)[i])); 1973 1974 // segm.CalcVertexShapes (x != NULL, dxdxi != NULL); 1975 segm.CalcVertexShapes (); 1976 1977 (*x)[i] = Point<3>(0,0,0); 1978 (*dxdxi)[i] = Vec<3>(0,0,0); 1979 1980 for (int v = 0; v < 2; v++) 1981 { 1982 (*x)[i] = (*x)[i] + segm.GetVertexShape(v) * mesh.Point(segm.GetVertexNr(v)); 1983 (*dxdxi)[i] = (*dxdxi)[i] + segm.GetVertexDShape(v) * mesh.Point(segm.GetVertexNr(v)); 1984 } 1985 1986 if (segm.GetEdgeOrder() > 1) 1987 { 1988 // segm.CalcEdgeShapes (x != NULL, dxdxi != NULL); 1989 segm.CalcEdgeShapes (); 1990 1991 int gindex = edgecoeffsindex[segm.GetEdgeNr()-1]; 1992 for (int k = 2; k <= segm.GetEdgeOrder(); k++, gindex++) 1993 { 1994 (*x)[i] = (*x)[i] + segm.GetEdgeShape(k-2) * edgecoeffs[gindex]; 1995 (*dxdxi)[i] = (*dxdxi)[i] + segm.GetEdgeDShape(k-2) * edgecoeffs[gindex]; 1996 } 1997 } 1998 } 1999 } 2000 2001 2002 2003 2004 CalcSurfaceTransformation(Point<2> xi,int elnr,Point<3> * x,Mat<3,2> * dxdxi)2005 void CurvedElements :: CalcSurfaceTransformation (Point<2> xi, int elnr, 2006 Point<3> * x, Mat<3,2> * dxdxi) 2007 { 2008 if (mesh.coarsemesh) 2009 { 2010 const HPRefElement & hpref_el = 2011 (*mesh.hpelements) [ mesh[(SurfaceElementIndex) elnr].hp_elnr]; 2012 2013 // xi umrechnen 2014 double lami[4]; 2015 FlatVector vlami(4, lami); 2016 vlami = 0; 2017 mesh[(SurfaceElementIndex) elnr] . GetShapeNew (xi, vlami); 2018 2019 Mat<2,2> trans; 2020 Mat<3,2> dxdxic; 2021 if (dxdxi) 2022 { 2023 MatrixFixWidth<2> dlami(4); 2024 dlami = 0; 2025 mesh[(SurfaceElementIndex) elnr] . GetDShapeNew (xi, dlami); 2026 2027 trans = 0; 2028 for (int k = 0; k < 2; k++) 2029 for (int l = 0; l < 2; l++) 2030 { 2031 double sum = 0; 2032 for (int i = 0; i < 4; i++) 2033 sum += hpref_el.param[i][l] * dlami(i, k); 2034 trans(l,k) = sum; 2035 } 2036 } 2037 2038 Point<2> coarse_xi(0,0); 2039 for (int i = 0; i < 4; i++) 2040 { 2041 coarse_xi(0) += hpref_el.param[i][0] * lami[i]; 2042 coarse_xi(1) += hpref_el.param[i][1] * lami[i]; 2043 } 2044 mesh.coarsemesh->GetCurvedElements().CalcSurfaceTransformation (coarse_xi, hpref_el.coarse_elnr, x, &dxdxic); 2045 2046 if (dxdxi) 2047 *dxdxi = dxdxic * trans; 2048 2049 return; 2050 } 2051 2052 2053 2054 2055 const Element2d & elem = mesh[(SurfaceElementIndex) elnr]; 2056 2057 BaseFiniteElement2D * fe2d; 2058 2059 // char locmem[max2(sizeof(FEQuad), sizeof(FETrig))]; 2060 char locmemtrig[sizeof(FETrig)]; 2061 char locmemquad[sizeof(FEQuad)]; 2062 switch (elem.GetType()) 2063 { 2064 case TRIG: fe2d = new (locmemtrig) FETrig (*this); break; 2065 case QUAD: fe2d = new (locmemquad) FEQuad (*this); break; 2066 } 2067 2068 // fe2d->SetElementNumber (elnr+1); 2069 fe2d->SetReferencePoint (xi); 2070 fe2d->CalcVertexShapes (); 2071 2072 if (x) 2073 { 2074 (*x) = Point<3>(0,0,0); 2075 for (int v = 0; v < fe2d->GetNVertices(); v++) 2076 (*x) = (*x) + fe2d->GetVertexShape(v) * mesh.Point(elem[v]); 2077 // (*x) = (*x) + fe2d->GetVertexShape(v) * mesh.Point(fe2d->GetVertexNr(v)); 2078 } 2079 2080 if (dxdxi) 2081 { 2082 for (int i = 0; i < 3; i++) 2083 for (int j = 0; j < 2; j++) 2084 (*dxdxi)(i,j) = 0; 2085 2086 for (int v = 0; v < elem.GetNV(); v++) 2087 for (int i = 0; i < 3; i++) 2088 for (int j = 0; j < 2; j++) 2089 (*dxdxi)(i,j) += fe2d->GetVertexDShape(v)(j) * mesh.Point(elem[v]).X(i+1); 2090 // (*dxdxi)(i,j) += fe2d->GetVertexDShape(v)(j) * mesh.Point(fe2d->GetVertexNr(v)).X(i+1); 2091 } 2092 2093 if (IsHighOrder()) 2094 { 2095 fe2d->SetElementNumber (elnr+1); 2096 2097 // fe2d->CalcEdgeShapes (x != NULL, dxdxi != NULL); 2098 fe2d->CalcEdgeShapes (); 2099 if (x) 2100 { 2101 int index = 0; 2102 for (int e = 0; e < fe2d->GetNEdges(); e++) 2103 { 2104 int gindex = edgecoeffsindex[fe2d->GetEdgeNr(e)-1]; 2105 2106 for (int k = 2; k <= fe2d->GetEdgeOrder(e); k++, index++, gindex++) 2107 (*x) = (*x) + fe2d->GetEdgeShape(index) * edgecoeffs[gindex]; 2108 } 2109 } 2110 2111 if (dxdxi) 2112 { 2113 int index = 0; 2114 for (int e = 0; e < fe2d->GetNEdges(); e++) 2115 { 2116 int gindex = edgecoeffsindex[fe2d->GetEdgeNr(e)-1]; 2117 for (int k = 2; k <= fe2d->GetEdgeOrder(e); k++, index++, gindex++) 2118 for (int i = 0; i < 3; i++) 2119 for (int j = 0; j < 2; j++) 2120 (*dxdxi)(i,j) += fe2d->GetEdgeDShape(index)(j) * edgecoeffs[gindex](i); 2121 } 2122 } 2123 2124 if (mesh.GetDimension() == 3) 2125 { 2126 // fe2d->CalcFaceShapes (x != NULL, dxdxi != NULL); 2127 fe2d->CalcFaceShapes (); 2128 2129 if (x) 2130 { 2131 int gindex = facecoeffsindex[fe2d->GetFaceNr()-1]; 2132 for (int index = 0; index < fe2d->GetNFaceShapes(); index++, gindex++) 2133 { 2134 (*x) = (*x) + fe2d->GetFaceShape(index) * facecoeffs[gindex]; 2135 } 2136 } 2137 2138 if (dxdxi) 2139 { 2140 int gindex = facecoeffsindex[fe2d->GetFaceNr()-1]; 2141 for (int index = 0; index < fe2d->GetNFaceShapes(); index++, gindex++) 2142 for (int i = 0; i < 3; i++) 2143 for (int j = 0; j < 2; j++) 2144 (*dxdxi)(i,j) += fe2d->GetFaceDShape(index)(j) * facecoeffs[gindex](i); 2145 } 2146 } 2147 } 2148 2149 fe2d -> ~BaseFiniteElement2D(); 2150 } 2151 2152 2153 2154 2155 2156 2157 2158 CalcMultiPointSurfaceTransformation(ARRAY<Point<2>> * xi,int elnr,ARRAY<Point<3>> * x,ARRAY<Mat<3,2>> * dxdxi)2159 void CurvedElements :: CalcMultiPointSurfaceTransformation (ARRAY< Point<2> > * xi, int elnr, 2160 ARRAY< Point<3> > * x, 2161 ARRAY< Mat<3,2> > * dxdxi) 2162 { 2163 2164 int size = xi->Size(); 2165 2166 if (mesh.coarsemesh) 2167 { 2168 const HPRefElement & hpref_el = 2169 (*mesh.hpelements) [ mesh[(SurfaceElementIndex) elnr].hp_elnr]; 2170 2171 // xi umrechnen 2172 ARRAY< Point<2> > coarse_xi; 2173 ARRAY< Mat<2,2> > trans; 2174 ARRAY< Mat<3,2> > dxdxic; 2175 2176 coarse_xi.SetSize (size); 2177 coarse_xi = 0; 2178 trans.SetSize (size); 2179 dxdxic.SetSize (size); 2180 2181 for (int mp = 0; mp<size; mp++) 2182 { 2183 double lami[4]; 2184 FlatVector vlami(4, lami); 2185 vlami = 0; 2186 mesh[(SurfaceElementIndex) elnr] . GetShapeNew ((*xi)[mp], vlami); 2187 2188 MatrixFixWidth<2> dlami(4); 2189 dlami = 0; 2190 mesh[(SurfaceElementIndex) elnr] . GetDShapeNew ((*xi)[mp], dlami); 2191 2192 trans[mp] = 0; 2193 for (int k = 0; k < 2; k++) 2194 for (int l = 0; l < 2; l++) 2195 { 2196 double sum = 0; 2197 for (int i = 0; i < 4; i++) 2198 sum += hpref_el.param[i][l] * dlami(i, k); 2199 trans[mp](l,k) = sum; 2200 } 2201 2202 for (int i = 0; i < 4; i++) 2203 { 2204 coarse_xi[mp](0) += hpref_el.param[i][0] * lami[i]; 2205 coarse_xi[mp](1) += hpref_el.param[i][1] * lami[i]; 2206 } 2207 } 2208 2209 mesh.coarsemesh->GetCurvedElements().CalcMultiPointSurfaceTransformation (&coarse_xi, hpref_el.coarse_elnr, x, &dxdxic); 2210 2211 for (int mp = 0; mp < size; mp++) 2212 { 2213 (*dxdxi)[mp] = dxdxic[mp] * trans[mp]; 2214 } 2215 2216 return; 2217 } 2218 2219 2220 x->SetSize (size); 2221 dxdxi->SetSize (size); 2222 2223 const Element2d & elem = mesh[(SurfaceElementIndex) elnr]; 2224 2225 BaseFiniteElement2D * fe2d; 2226 2227 // char locmem[max2(sizeof(FEQuad), sizeof(FETrig))]; 2228 char locmemtrig[sizeof(FETrig)]; 2229 char locmemquad[sizeof(FEQuad)]; 2230 switch (elem.GetType()) 2231 { 2232 case TRIG: fe2d = new (locmemtrig) FETrig (*this); break; 2233 case QUAD: fe2d = new (locmemquad) FEQuad (*this); break; 2234 } 2235 2236 fe2d->SetElementNumber (elnr+1); 2237 2238 for (int mp = 0; mp < size; mp++) 2239 { 2240 fe2d->SetReferencePoint ((*xi)[mp]); 2241 fe2d->CalcVertexShapes (); 2242 2243 (*x)[mp] = Point<3>(0,0,0); 2244 2245 for (int i = 0; i < 3; i++) 2246 for (int j = 0; j < 2; j++) 2247 (*dxdxi)[mp](i,j) = 0; 2248 2249 for (int v = 0; v < fe2d->GetNVertices(); v++) 2250 { 2251 (*x)[mp] = (*x)[mp] + fe2d->GetVertexShape(v) * mesh.Point(fe2d->GetVertexNr(v)); 2252 for (int i = 0; i < 3; i++) 2253 for (int j = 0; j < 2; j++) 2254 (*dxdxi)[mp](i,j) += fe2d->GetVertexDShape(v)(j) * mesh.Point(fe2d->GetVertexNr(v)).X(i+1); 2255 } 2256 2257 if (IsHighOrder()) 2258 { 2259 // fe2d->CalcEdgeShapes (x != NULL, dxdxi != NULL); 2260 fe2d->CalcEdgeShapes (); 2261 2262 int index = 0; 2263 for (int e = 0; e < fe2d->GetNEdges(); e++) 2264 { 2265 int gindex = edgecoeffsindex[fe2d->GetEdgeNr(e)-1]; 2266 2267 for (int k = 2; k <= fe2d->GetEdgeOrder(e); k++, index++, gindex++) 2268 { 2269 (*x)[mp] = (*x)[mp] + fe2d->GetEdgeShape(index) * edgecoeffs[gindex]; 2270 for (int i = 0; i < 3; i++) 2271 for (int j = 0; j < 2; j++) 2272 (*dxdxi)[mp](i,j) += fe2d->GetEdgeDShape(index)(j) * edgecoeffs[gindex](i); 2273 } 2274 } 2275 2276 2277 if (mesh.GetDimension() == 3) 2278 { 2279 // fe2d->CalcFaceShapes (x != NULL, dxdxi != NULL); 2280 fe2d->CalcFaceShapes (); 2281 2282 int gindex = facecoeffsindex[fe2d->GetFaceNr()-1]; 2283 for (int index = 0; index < fe2d->GetNFaceShapes(); index++, gindex++) 2284 { 2285 (*x)[mp] = (*x)[mp] + fe2d->GetFaceShape(index) * facecoeffs[gindex]; 2286 for (int i = 0; i < 3; i++) 2287 for (int j = 0; j < 2; j++) 2288 (*dxdxi)[mp](i,j) += fe2d->GetFaceDShape(index)(j) * facecoeffs[gindex](i); 2289 } 2290 } 2291 } 2292 } 2293 2294 fe2d -> ~BaseFiniteElement2D(); 2295 } 2296 2297 2298 2299 2300 2301 2302 2303 2304 2305 2306 2307 2308 2309 2310 2311 2312 2313 2314 CalcElementTransformation(Point<3> xi,int elnr,Point<3> * x,Mat<3,3> * dxdxi)2315 void CurvedElements :: CalcElementTransformation (Point<3> xi, int elnr, 2316 Point<3> * x, Mat<3,3> * dxdxi) 2317 { 2318 2319 if (mesh.coarsemesh) 2320 { 2321 const HPRefElement & hpref_el = 2322 (*mesh.hpelements) [ mesh[(ElementIndex) elnr].hp_elnr]; 2323 2324 2325 /* 2326 *testout << " Curved Element " << elnr << endl; 2327 *testout << " hpref_el.coarse_elnr " << hpref_el.coarse_elnr << endl; 2328 *testout << " hpref_el.param = " << endl; 2329 for(int j=0;j< hpref_el.np; j++) 2330 { 2331 for(int k=0;k<3;k++) 2332 *testout << hpref_el.param[j][k] << "\t"; 2333 *testout << endl; 2334 } 2335 */ 2336 2337 double lami[8]; 2338 FlatVector vlami(8, lami); 2339 vlami = 0; 2340 mesh[(ElementIndex) elnr] . GetShapeNew (xi, vlami); 2341 2342 Point<3> coarse_xi(0,0,0); 2343 for (int i = 0; i < hpref_el.np; i++) 2344 for (int l = 0; l < 3; l++) 2345 coarse_xi(l) += hpref_el.param[i][l] * lami[i]; 2346 2347 Mat<3,3> trans, dxdxic; 2348 if (dxdxi) 2349 { 2350 MatrixFixWidth<3> dlami(8); 2351 dlami = 0; 2352 mesh[(ElementIndex) elnr] . GetDShapeNew (xi, dlami); 2353 2354 trans = 0; 2355 for (int k = 0; k < 3; k++) 2356 for (int l = 0; l < 3; l++) 2357 { 2358 double sum = 0; 2359 for (int i = 0; i < hpref_el.np; i++) 2360 sum += hpref_el.param[i][l] * dlami(i, k); 2361 trans(l,k) = sum; 2362 } 2363 } 2364 /* 2365 *testout << " x " << x << endl; 2366 // << "\t" << x.X(2) << "\t" << x.X(3) << endl; 2367 *testout << " Element Trafo " << coarse_xi(0) << " \t " << coarse_xi(1) << " \t " << coarse_xi(2) << endl; 2368 */ 2369 2370 2371 2372 mesh.coarsemesh->GetCurvedElements().CalcElementTransformation (coarse_xi, hpref_el.coarse_elnr, x, &dxdxic); 2373 2374 if (dxdxi) 2375 *dxdxi = dxdxic * trans; 2376 return; 2377 } 2378 2379 2380 2381 2382 2383 2384 2385 2386 2387 Element elem = mesh[(ElementIndex) elnr]; 2388 BaseFiniteElement3D * fe3d; 2389 2390 // char locmem[max2(sizeof(FETet), sizeof(FEPrism))]; 2391 char locmemtet[sizeof(FETet)]; 2392 char locmemprism[sizeof(FEPrism)]; 2393 char locmempyramid[sizeof(FEPyramid)]; 2394 char locmemhex[sizeof(FEHex)]; 2395 switch (elem.GetType()) 2396 { 2397 case TET: fe3d = new (locmemtet) FETet (*this); break; 2398 case PRISM: fe3d = new (locmemprism) FEPrism (*this); break; 2399 case PYRAMID: fe3d = new (locmempyramid) FEPyramid (*this); break; 2400 case HEX: fe3d = new (locmemhex) FEHex (*this); break; 2401 } 2402 2403 // fe3d->SetElementNumber (elnr+1); 2404 fe3d->SetReferencePoint (xi); 2405 2406 fe3d->CalcVertexShapes (); 2407 // fe3d->CalcVertexShapes (x != NULL, dxdxi != NULL); 2408 2409 2410 if (x) 2411 { 2412 (*x) = Point<3>(0,0,0); 2413 for (int v = 0; v < fe3d->GetNVertices(); v++) 2414 (*x) += fe3d->GetVertexShape(v) * Vec<3> (mesh.Point(elem[v])); 2415 } 2416 2417 if (dxdxi) 2418 { 2419 for (int i = 0; i < 3; i++) 2420 for (int j = 0; j < 3; j++) 2421 (*dxdxi)(i,j) = 0; 2422 2423 for (int v = 0; v < fe3d->GetNVertices(); v++) 2424 for (int i = 0; i < 3; i++) 2425 for (int j = 0; j < 3; j++) 2426 (*dxdxi)(i,j) += fe3d->GetVertexDShape(v)(j) * mesh.Point(elem[v]).X(i+1); 2427 } 2428 2429 if (IsHighOrder()) 2430 { 2431 fe3d->SetElementNumber (elnr+1); 2432 // fe3d->CalcEdgeShapes (x != NULL, dxdxi != NULL); 2433 fe3d->CalcEdgeShapes (); 2434 2435 2436 if (x) 2437 { 2438 int index = 0; 2439 for (int e = 0; e < fe3d->GetNEdges(); e++) 2440 { 2441 int gindex = edgecoeffsindex[fe3d->GetEdgeNr(e)-1]; 2442 for (int k = 2; k <= fe3d->GetEdgeOrder(e); k++, index++, gindex++) 2443 (*x) += fe3d->GetEdgeShape(index) * edgecoeffs[gindex]; 2444 } 2445 } 2446 2447 if (dxdxi) 2448 { 2449 int index = 0; 2450 for (int e = 0; e < fe3d->GetNEdges(); e++) 2451 { 2452 int gindex = edgecoeffsindex[fe3d->GetEdgeNr(e)-1]; 2453 for (int k = 2; k <= fe3d->GetEdgeOrder(e); k++, index++, gindex++) 2454 for (int i = 0; i < 3; i++) 2455 for (int j = 0; j < 3; j++) 2456 (*dxdxi)(i,j) += fe3d->GetEdgeDShape(index)(j) * edgecoeffs[gindex](i); 2457 } 2458 } 2459 2460 2461 // cout << "switched off face mapping" << endl; 2462 if (mesh.GetDimension() == 3) // JS 2463 { 2464 fe3d->CalcFaceShapes (); 2465 // fe3d->CalcFaceShapes (x != NULL, dxdxi != NULL); 2466 2467 if (x) 2468 { 2469 int index = 0; 2470 for (int f = 0; f < fe3d->GetNFaces(); f++) 2471 { 2472 int gindex = facecoeffsindex[fe3d->GetFaceNr(f)-1]; 2473 for (int k = 0; k < fe3d->GetNFaceShapes(f); k++, index++, gindex++) 2474 (*x) += fe3d->GetFaceShape(index) * facecoeffs[gindex]; 2475 } 2476 } 2477 2478 if (dxdxi) 2479 { 2480 int index = 0; 2481 for (int f = 0; f < fe3d->GetNFaces(); f++) 2482 { 2483 int gindex = facecoeffsindex[fe3d->GetFaceNr(f)-1]; 2484 for (int k = 0; k < fe3d->GetNFaceShapes(f); k++, index++, gindex++) 2485 for (int i = 0; i < 3; i++) 2486 for (int j = 0; j < 3; j++) 2487 (*dxdxi)(i,j) += fe3d->GetFaceDShape(index)(j) * facecoeffs[gindex](i); 2488 } 2489 } 2490 } 2491 } 2492 2493 fe3d -> ~BaseFiniteElement3D(); 2494 } 2495 2496 2497 2498 2499 2500 2501 2502 2503 2504 #ifdef ROBERT 2505 CalcMultiPointElementTransformation(ARRAY<Point<3>> * xi,int elnr,ARRAY<Point<3>> * x,ARRAY<Mat<3,3>> * dxdxi)2506 void CurvedElements :: CalcMultiPointElementTransformation (ARRAY< Point<3> > * xi, int elnr, 2507 ARRAY< Point<3> > * x, 2508 ARRAY< Mat<3,3> > * dxdxi) 2509 { 2510 int size = (*xi).Size(); 2511 x->SetSize (size); 2512 2513 if (dxdxi) dxdxi->SetSize (size); 2514 2515 if (mesh.coarsemesh) 2516 { 2517 const HPRefElement & hpref_el = 2518 (*mesh.hpelements) [ mesh[(ElementIndex) elnr].hp_elnr]; 2519 2520 ARRAY< Mat<3,3> > trans, dxdxic; 2521 if (dxdxi) 2522 { 2523 trans.SetSize (size); 2524 dxdxic.SetSize (size); 2525 } 2526 2527 ARRAY<Point<3> > coarse_xi(size); 2528 2529 double lami[8]; 2530 FlatVector vlami(8, lami); 2531 const Element & el = mesh[(ElementIndex) elnr]; 2532 int el_np = el.GetNP(); 2533 2534 for (int mp = 0; mp < size; mp++) 2535 { 2536 el.GetShapeNew ((*xi)[mp], vlami); 2537 2538 Point<3> pc(0,0,0); 2539 for (int i = 0; i < el_np; i++) 2540 for (int l = 0; l < 3; l++) 2541 pc(l) += hpref_el.param[i][l] * lami[i]; 2542 2543 coarse_xi[mp] = pc; 2544 2545 if (dxdxi) 2546 { 2547 MatrixFixWidth<3> dlami(8); 2548 dlami = 0; 2549 mesh[(ElementIndex) elnr] . GetDShapeNew ((*xi)[mp], dlami); 2550 2551 trans[mp] = 0; 2552 for (int k = 0; k < 3; k++) 2553 for (int l = 0; l < 3; l++) 2554 { 2555 double sum = 0; 2556 for (int i = 0; i < 8; i++) 2557 sum += hpref_el.param[i][l] * dlami(i, k); 2558 trans[mp](l,k) = sum; 2559 } 2560 } 2561 } 2562 2563 if (dxdxi) 2564 mesh.coarsemesh->GetCurvedElements().CalcMultiPointElementTransformation (&coarse_xi, hpref_el.coarse_elnr, x, &dxdxic); 2565 else 2566 mesh.coarsemesh->GetCurvedElements().CalcMultiPointElementTransformation (&coarse_xi, hpref_el.coarse_elnr, x, 0); 2567 2568 if (dxdxi) 2569 for (int mp = 0; mp < size; mp++) 2570 (*dxdxi)[mp] = dxdxic[mp] * trans[mp]; 2571 2572 return; 2573 } 2574 2575 2576 2577 2578 2579 2580 2581 2582 2583 Element elem = mesh[(ElementIndex) elnr]; 2584 BaseFiniteElement3D * fe3d; 2585 2586 // char locmem[max2(sizeof(FETet), sizeof(FEPrism))]; 2587 char locmemtet[sizeof(FETet)]; 2588 char locmemprism[sizeof(FEPrism)]; 2589 char locmempyramid[sizeof(FEPyramid)]; 2590 char locmemhex[sizeof(FEHex)]; 2591 switch (elem.GetType()) 2592 { 2593 case TET: fe3d = new (locmemtet) FETet (*this); break; 2594 case PRISM: fe3d = new (locmemprism) FEPrism (*this); break; 2595 case PYRAMID: fe3d = new (locmempyramid) FEPyramid (*this); break; 2596 case HEX: fe3d = new (locmemhex) FEHex (*this); break; 2597 } 2598 2599 fe3d->SetElementNumber (elnr+1); 2600 2601 2602 for (int mp = 0; mp < size; mp++) 2603 { 2604 fe3d->SetReferencePoint ((*xi)[mp]); 2605 2606 fe3d->CalcVertexShapes (); 2607 // fe3d->CalcVertexShapes (x != NULL, dxdxi != NULL); 2608 2609 (*x)[mp] = Point<3>(0,0,0); 2610 if (dxdxi) 2611 for (int i = 0; i < 3; i++) 2612 for (int j = 0; j < 3; j++) 2613 (*dxdxi)[mp](i,j) = 0; 2614 2615 for (int v = 0; v < fe3d->GetNVertices(); v++) 2616 { 2617 (*x)[mp] += fe3d->GetVertexShape(v) * Vec<3> (mesh.Point(fe3d->GetVertexNr(v))); 2618 if (dxdxi) 2619 for (int i = 0; i < 3; i++) 2620 for (int j = 0; j < 3; j++) 2621 (*dxdxi)[mp](i,j) += fe3d->GetVertexDShape(v)(j) * mesh.Point(fe3d->GetVertexNr(v)).X(i+1); 2622 } 2623 2624 2625 if (IsHighOrder()) 2626 { 2627 // fe3d->CalcEdgeShapes (x != NULL, dxdxi != NULL); 2628 fe3d->CalcEdgeShapes (); 2629 2630 int index = 0; 2631 for (int e = 0; e < fe3d->GetNEdges(); e++) 2632 { 2633 int gindex = edgecoeffsindex[fe3d->GetEdgeNr(e)-1]; 2634 for (int k = 2; k <= fe3d->GetEdgeOrder(e); k++, index++, gindex++) 2635 { 2636 (*x)[mp] += fe3d->GetEdgeShape(index) * edgecoeffs[gindex]; 2637 if (dxdxi) 2638 for (int i = 0; i < 3; i++) 2639 for (int j = 0; j < 3; j++) 2640 (*dxdxi)[mp](i,j) += fe3d->GetEdgeDShape(index)(j) * edgecoeffs[gindex](i); 2641 } 2642 } 2643 2644 if (mesh.GetDimension() == 3) 2645 { 2646 fe3d->CalcFaceShapes (); 2647 // fe3d->CalcFaceShapes (x != NULL, dxdxi != NULL); 2648 2649 int index = 0; 2650 for (int f = 0; f < fe3d->GetNFaces(); f++) 2651 { 2652 int gindex = facecoeffsindex[fe3d->GetFaceNr(f)-1]; 2653 for (int k = 0; k < fe3d->GetNFaceShapes(f); k++, index++, gindex++) 2654 { 2655 (*x)[mp] += fe3d->GetFaceShape(index) * facecoeffs[gindex]; 2656 if (dxdxi) 2657 for (int i = 0; i < 3; i++) 2658 for (int j = 0; j < 3; j++) 2659 (*dxdxi)[mp](i,j) += fe3d->GetFaceDShape(index)(j) * facecoeffs[gindex](i); 2660 } 2661 } 2662 } 2663 } 2664 } 2665 2666 fe3d -> ~BaseFiniteElement3D(); 2667 } 2668 2669 #endif 2670 2671 2672 2673 2674 2675 CalcMultiPointElementTransformation(ARRAY<Point<3>> * xi,int elnr,ARRAY<Point<3>> * x,ARRAY<Mat<3,3>> * dxdxi)2676 void CurvedElements :: CalcMultiPointElementTransformation (ARRAY< Point<3> > * xi, int elnr, 2677 ARRAY< Point<3> > * x, 2678 ARRAY< Mat<3,3> > * dxdxi) 2679 { 2680 int size = (*xi).Size(); 2681 x->SetSize (size); 2682 2683 if (dxdxi) dxdxi->SetSize (size); 2684 2685 if (mesh.coarsemesh) 2686 { 2687 const HPRefElement & hpref_el = 2688 (*mesh.hpelements) [ mesh[(ElementIndex) elnr].hp_elnr]; 2689 2690 ARRAY< Mat<3,3> > trans, dxdxic; 2691 if (dxdxi) 2692 { 2693 trans.SetSize (size); 2694 dxdxic.SetSize (size); 2695 } 2696 2697 ARRAY<Point<3> > coarse_xi(size); 2698 2699 double lami[8]; 2700 FlatVector vlami(8, lami); 2701 const Element & el = mesh[(ElementIndex) elnr]; 2702 int el_np = el.GetNP(); 2703 2704 for (int mp = 0; mp < size; mp++) 2705 { 2706 el.GetShapeNew ((*xi)[mp], vlami); 2707 2708 Point<3> pc(0,0,0); 2709 for (int i = 0; i < el_np; i++) 2710 for (int l = 0; l < 3; l++) 2711 pc(l) += hpref_el.param[i][l] * lami[i]; 2712 2713 coarse_xi[mp] = pc; 2714 2715 if (dxdxi) 2716 { 2717 MatrixFixWidth<3> dlami(8); 2718 dlami = 0; 2719 mesh[(ElementIndex) elnr] . GetDShapeNew ((*xi)[mp], dlami); 2720 2721 trans[mp] = 0; 2722 for (int k = 0; k < 3; k++) 2723 for (int l = 0; l < 3; l++) 2724 { 2725 double sum = 0; 2726 for (int i = 0; i < 8; i++) 2727 sum += hpref_el.param[i][l] * dlami(i, k); 2728 trans[mp](l,k) = sum; 2729 } 2730 } 2731 } 2732 2733 if (dxdxi) 2734 mesh.coarsemesh->GetCurvedElements().CalcMultiPointElementTransformation (&coarse_xi, hpref_el.coarse_elnr, x, &dxdxic); 2735 else 2736 mesh.coarsemesh->GetCurvedElements().CalcMultiPointElementTransformation (&coarse_xi, hpref_el.coarse_elnr, x, 0); 2737 2738 if (dxdxi) 2739 for (int mp = 0; mp < size; mp++) 2740 (*dxdxi)[mp] = dxdxic[mp] * trans[mp]; 2741 2742 return; 2743 } 2744 2745 2746 2747 2748 2749 2750 2751 2752 2753 Element elem = mesh[(ElementIndex) elnr]; 2754 BaseFiniteElement3D * fe3d; 2755 2756 // char locmem[max2(sizeof(FETet), sizeof(FEPrism))]; 2757 char locmemtet[sizeof(FETet)]; 2758 char locmemprism[sizeof(FEPrism)]; 2759 char locmempyramid[sizeof(FEPyramid)]; 2760 char locmemhex[sizeof(FEHex)]; 2761 switch (elem.GetType()) 2762 { 2763 case TET: fe3d = new (locmemtet) FETet (*this); break; 2764 case PRISM: fe3d = new (locmemprism) FEPrism (*this); break; 2765 case PYRAMID: fe3d = new (locmempyramid) FEPyramid (*this); break; 2766 case HEX: fe3d = new (locmemhex) FEHex (*this); break; 2767 } 2768 2769 fe3d->SetElementNumber (elnr+1); 2770 2771 if (dxdxi) 2772 for (int mp = 0; mp < size; mp++) 2773 (*dxdxi)[mp] = 0.0; 2774 2775 Vec<3> vertices[8]; 2776 ArrayMem<Vec<3>, 100> loc_edgecoefs(0); 2777 ArrayMem<Vec<3>, 100> loc_facecoefs(0); 2778 2779 for (int v = 0; v < fe3d->GetNVertices(); v++) 2780 vertices[v] = Vec<3> (mesh.Point(fe3d->GetVertexNr(v))); 2781 2782 if (IsHighOrder()) 2783 { 2784 for (int e = 0; e < fe3d->GetNEdges(); e++) 2785 { 2786 int gindex = edgecoeffsindex[fe3d->GetEdgeNr(e)-1]; 2787 for (int k = 2; k <= fe3d->GetEdgeOrder(e); k++, gindex++) 2788 loc_edgecoefs.Append (edgecoeffs[gindex]); 2789 } 2790 2791 if (mesh.GetDimension() == 3) 2792 for (int f = 0; f < fe3d->GetNFaces(); f++) 2793 { 2794 int gindex = facecoeffsindex[fe3d->GetFaceNr(f)-1]; 2795 for (int k = 0; k < fe3d->GetNFaceShapes(f); k++, gindex++) 2796 loc_facecoefs.Append (facecoeffs[gindex]); 2797 } 2798 } 2799 2800 2801 2802 if (!dxdxi) 2803 2804 for (int mp = 0; mp < size; mp++) 2805 { 2806 fe3d->SetReferencePoint ((*xi)[mp]); 2807 fe3d->CalcVertexShapesOnly (); 2808 2809 if (IsHighOrder()) 2810 { 2811 fe3d->CalcEdgeShapesOnly (); 2812 if (mesh.GetDimension() == 3) 2813 fe3d->CalcFaceShapes (); 2814 } 2815 2816 Point<3> hp (0,0,0); 2817 for (int v = 0; v < fe3d->GetNVertices(); v++) 2818 hp += fe3d->GetVertexShape(v) * vertices[v]; 2819 for (int index = 0; index < loc_edgecoefs.Size(); index++) 2820 hp += fe3d->GetEdgeShape(index) * loc_edgecoefs[index]; 2821 for (int index = 0; index < loc_facecoefs.Size(); index++) 2822 hp += fe3d->GetFaceShape(index) * loc_facecoefs[index]; // edgecoeffs[gindex]; 2823 (*x)[mp] = hp; 2824 } 2825 2826 else 2827 2828 for (int mp = 0; mp < size; mp++) 2829 { 2830 fe3d->SetReferencePoint ((*xi)[mp]); 2831 fe3d->CalcVertexShapes (); 2832 2833 if (IsHighOrder()) 2834 { 2835 fe3d->CalcEdgeShapes (); 2836 if (mesh.GetDimension() == 3) 2837 fe3d->CalcFaceShapes (); 2838 } 2839 2840 Point<3> hp (0,0,0); 2841 for (int v = 0; v < fe3d->GetNVertices(); v++) 2842 hp += fe3d->GetVertexShape(v) * vertices[v]; 2843 for (int index = 0; index < loc_edgecoefs.Size(); index++) 2844 hp += fe3d->GetEdgeShape(index) * loc_edgecoefs[index]; 2845 for (int index = 0; index < loc_facecoefs.Size(); index++) 2846 hp += fe3d->GetFaceShape(index) * loc_facecoefs[index]; // edgecoeffs[gindex]; 2847 (*x)[mp] = hp; 2848 2849 2850 for (int v = 0; v < fe3d->GetNVertices(); v++) 2851 for (int i = 0; i < 3; i++) 2852 for (int j = 0; j < 3; j++) 2853 (*dxdxi)[mp](i,j) += fe3d->GetVertexDShape(v)(j) * vertices[v](i); 2854 2855 for (int index = 0; index < loc_edgecoefs.Size(); index++) 2856 for (int i = 0; i < 3; i++) 2857 for (int j = 0; j < 3; j++) 2858 (*dxdxi)[mp](i,j) += fe3d->GetEdgeDShape(index)(j) * loc_edgecoefs[index](i); 2859 2860 for (int index = 0; index < loc_facecoefs.Size(); index++) 2861 for (int i = 0; i < 3; i++) 2862 for (int j = 0; j < 3; j++) 2863 (*dxdxi)[mp](i,j) += fe3d->GetFaceDShape(index)(j) * loc_facecoefs[index](i); 2864 } 2865 2866 fe3d -> ~BaseFiniteElement3D(); 2867 } 2868 2869 2870 2871 /* 2872 class Init 2873 { 2874 PolynomialBasis b; 2875 public: 2876 Init (); 2877 }; 2878 2879 Init :: Init() 2880 { 2881 ofstream edge("edge.out"); 2882 b.SetOrder(10); 2883 for (double x = 0; x <= 1; x += 0.01) 2884 { 2885 b.CalcFScaled(x, 1); 2886 edge << x; 2887 for (int j = 2; j <= 5; j++) 2888 edge << " " << b.GetF(j); 2889 edge << endl; 2890 } 2891 } 2892 2893 Init in; 2894 */ 2895 2896 } // namespace netgen 2897 2898 2899 #endif 2900