|
Intrepid
|
00001 // @HEADER 00002 // ************************************************************************ 00003 // 00004 // Intrepid Package 00005 // Copyright (2007) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // Redistribution and use in source and binary forms, with or without 00011 // modification, are permitted provided that the following conditions are 00012 // met: 00013 // 00014 // 1. Redistributions of source code must retain the above copyright 00015 // notice, this list of conditions and the following disclaimer. 00016 // 00017 // 2. Redistributions in binary form must reproduce the above copyright 00018 // notice, this list of conditions and the following disclaimer in the 00019 // documentation and/or other materials provided with the distribution. 00020 // 00021 // 3. Neither the name of the Corporation nor the names of the 00022 // contributors may be used to endorse or promote products derived from 00023 // this software without specific prior written permission. 00024 // 00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00036 // 00037 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) 00038 // Denis Ridzal (dridzal@sandia.gov), or 00039 // Kara Peterson (kjpeter@sandia.gov) 00040 // 00041 // ************************************************************************ 00042 // @HEADER 00043 00048 #ifdef _MSC_VER 00049 #include "winmath.h" 00050 #endif 00051 00052 00053 namespace Intrepid { 00054 00055 template<class Scalar, class ArrayType> 00056 void PointTools::getLattice(ArrayType &pts , 00057 const shards::CellTopology& cellType , 00058 const int order , 00059 const int offset , 00060 const EPointType pointType ) 00061 { 00062 switch (pointType) { 00063 case POINTTYPE_EQUISPACED: 00064 getEquispacedLattice<Scalar,ArrayType>( cellType , pts , order , offset ); 00065 break; 00066 case POINTTYPE_WARPBLEND: 00067 getWarpBlendLattice<Scalar,ArrayType>( cellType , pts , order , offset ); 00068 break; 00069 default: 00070 TEUCHOS_TEST_FOR_EXCEPTION( true , 00071 std::invalid_argument , 00072 "PointTools::getLattice: invalid EPointType" ); 00073 } 00074 return; 00075 } 00076 00077 template<class Scalar, class ArrayType> 00078 void PointTools::getGaussPoints( ArrayType &pts , 00079 const int order ) 00080 { 00081 00082 Scalar *z = new Scalar[order+1]; 00083 Scalar *w = new Scalar[order+1]; 00084 00085 IntrepidPolylib::zwgj( z , w , order + 1 , 0.0 , 0.0 ); 00086 for (int i=0;i<order+1;i++) { 00087 pts(i,0) = z[i]; 00088 } 00089 00090 delete []z; 00091 delete []w; 00092 } 00093 00094 00095 template<class Scalar, class ArrayType> 00096 void PointTools::getEquispacedLattice(const shards::CellTopology& cellType , 00097 ArrayType &points , 00098 const int order , 00099 const int offset ) 00100 00101 { 00102 switch (cellType.getKey()) { 00103 case shards::Tetrahedron<4>::key: 00104 case shards::Tetrahedron<8>::key: 00105 case shards::Tetrahedron<10>::key: 00106 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) ) 00107 || points.dimension(1) != 3 , 00108 std::invalid_argument , 00109 ">>> ERROR(PointTools::getEquispacedLattice): points argument is ill-sized." ); 00110 getEquispacedLatticeTetrahedron<Scalar,ArrayType>( points , order , offset ); 00111 break; 00112 case shards::Triangle<3>::key: 00113 case shards::Triangle<4>::key: 00114 case shards::Triangle<6>::key: 00115 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) ) 00116 || points.dimension(1) != 2 , 00117 std::invalid_argument , 00118 ">>> ERROR(PointTools::getEquispacedLattice): points argument is ill-sized." ); 00119 getEquispacedLatticeTriangle<Scalar,ArrayType>( points , order , offset ); 00120 break; 00121 case shards::Line<2>::key: 00122 case shards::Line<3>::key: 00123 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) ) 00124 || points.dimension(1) != 1 , 00125 std::invalid_argument , 00126 ">>> ERROR(PointTools::getEquispacedLattice): points argument is ill-sized." ); 00127 getEquispacedLatticeLine<Scalar,ArrayType>( points , order , offset ); 00128 break; 00129 default: 00130 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , 00131 ">>> ERROR (Intrepid::PointTools::getEquispacedLattice): Illegal cell type" ); 00132 } 00133 00134 } 00135 00136 template<class Scalar, class ArrayType> 00137 void PointTools::getWarpBlendLattice( const shards::CellTopology& cellType , 00138 ArrayType &points , 00139 const int order , 00140 const int offset ) 00141 00142 { 00143 switch (cellType.getKey()) { 00144 case shards::Tetrahedron<4>::key: 00145 case shards::Tetrahedron<8>::key: 00146 case shards::Tetrahedron<10>::key: 00147 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) ) 00148 || points.dimension(1) != 3 , 00149 std::invalid_argument , 00150 ">>> ERROR(PointTools::getWarpBlendLattice): points argument is ill-sized." ); 00151 getWarpBlendLatticeTetrahedron<Scalar,ArrayType>( points , order , offset ); 00152 break; 00153 case shards::Triangle<3>::key: 00154 case shards::Triangle<4>::key: 00155 case shards::Triangle<6>::key: 00156 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) ) 00157 || points.dimension(1) != 2 , 00158 std::invalid_argument , 00159 ">>> ERROR(PointTools::getWarpBlendLattice): points argument is ill-sized." ); 00160 getWarpBlendLatticeTriangle<Scalar,ArrayType>( points , order , offset ); 00161 break; 00162 case shards::Line<2>::key: 00163 case shards::Line<3>::key: 00164 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) ) 00165 || points.dimension(1) != 1 , 00166 std::invalid_argument , 00167 ">>> ERROR(PointTools::getWarpBlendLattice): points argument is ill-sized." ); 00168 getWarpBlendLatticeLine<Scalar,ArrayType>( points , order , offset ); 00169 break; 00170 default: 00171 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , 00172 ">>> ERROR (Intrepid::PointTools::getWarpBlendLattice): Illegal cell type" ); 00173 } 00174 00175 } 00176 00177 template<class Scalar, class ArrayType> 00178 void PointTools::getEquispacedLatticeLine( ArrayType &points , 00179 const int order , 00180 const int offset ) 00181 { 00182 TEUCHOS_TEST_FOR_EXCEPTION( order < 0 , 00183 std::invalid_argument , 00184 ">>> ERROR (Intrepid::PointTools::getEquispacedLatticeLine): order must be positive" ); 00185 if (order == 0) { 00186 points(0,0) = 0.0; 00187 } 00188 else { 00189 const Scalar h = 2.0 / order; 00190 00191 for (int i=offset;i<=order-offset;i++) { 00192 points(i-offset,0) = -1.0 + h * (Scalar) i; 00193 } 00194 } 00195 00196 return; 00197 } 00198 00199 template<class Scalar, class ArrayType> 00200 void PointTools::getEquispacedLatticeTriangle( ArrayType &points , 00201 const int order , 00202 const int offset ) 00203 { 00204 TEUCHOS_TEST_FOR_EXCEPTION( order <= 0 , 00205 std::invalid_argument , 00206 ">>> ERROR (Intrepid::PointTools::getEquispacedLatticeLine): order must be positive" ); 00207 00208 const Scalar h = 1.0 / order; 00209 int cur = 0; 00210 00211 for (int i=offset;i<=order-offset;i++) { 00212 for (int j=offset;j<=order-i-offset;j++) { 00213 points(cur,0) = (Scalar)0.0 + (Scalar) j * h ; 00214 points(cur,1) = (Scalar)0.0 + (Scalar) i * h; 00215 cur++; 00216 } 00217 } 00218 00219 return; 00220 } 00221 00222 template<class Scalar, class ArrayType> 00223 void PointTools::getEquispacedLatticeTetrahedron( ArrayType &points , 00224 const int order , 00225 const int offset ) 00226 { 00227 TEUCHOS_TEST_FOR_EXCEPTION( (order <= 0) , 00228 std::invalid_argument , 00229 ">>> ERROR (Intrepid::PointTools::getEquispacedLatticeTetrahedron): order must be positive" ); 00230 00231 const Scalar h = 1.0 / order; 00232 int cur = 0; 00233 00234 for (int i=offset;i<=order-offset;i++) { 00235 for (int j=offset;j<=order-i-offset;j++) { 00236 for (int k=offset;k<=order-i-j-offset;k++) { 00237 points(cur,0) = (Scalar) k * h; 00238 points(cur,1) = (Scalar) j * h; 00239 points(cur,2) = (Scalar) i * h; 00240 cur++; 00241 } 00242 } 00243 } 00244 00245 return; 00246 } 00247 00248 template<class Scalar, class ArrayType> 00249 void PointTools::getWarpBlendLatticeLine( ArrayType &points , 00250 const int order , 00251 const int offset ) 00252 { 00253 Scalar *z = new Scalar[order+1]; 00254 Scalar *w = new Scalar[order+1]; 00255 00256 // order is order of polynomial degree. The Gauss-Lobatto points are accurate 00257 // up to degree 2*i-1 00258 IntrepidPolylib::zwglj( z , w , order+1 , 0.0 , 0.0 ); 00259 00260 for (int i=offset;i<=order-offset;i++) { 00261 points(i-offset,0) = z[i]; 00262 } 00263 00264 delete []z; 00265 delete []w; 00266 00267 return; 00268 } 00269 00270 template<class Scalar, class ArrayType> 00271 void PointTools::warpFactor( const int order , 00272 const ArrayType &xnodes , 00273 const ArrayType &xout , 00274 ArrayType &warp) 00275 { 00276 TEUCHOS_TEST_FOR_EXCEPTION( ( warp.dimension(0) != xout.dimension(0) ) , 00277 std::invalid_argument , 00278 ">>> ERROR (PointTools::warpFactor): xout and warp must be same size." ); 00279 00280 warp.initialize(); 00281 00282 FieldContainer<Scalar> d( xout.dimension(0) ); 00283 d.initialize(); 00284 00285 FieldContainer<Scalar> xeq( order + 1 ,1); 00286 PointTools::getEquispacedLatticeLine<Scalar,ArrayType>( xeq , order , 0 ); 00287 xeq.resize( order + 1 ); 00288 00289 TEUCHOS_TEST_FOR_EXCEPTION( ( xeq.dimension(0) != xnodes.dimension(0) ) , 00290 std::invalid_argument , 00291 ">>> ERROR (PointTools::warpFactor): xeq and xnodes must be same size." ); 00292 00293 for (int i=0;i<=order;i++) { 00294 00295 for (int k=0;k<d.dimension(0);k++) { 00296 d(k) = xnodes(i) - xeq(i); 00297 } 00298 00299 for (int j=1;j<order;j++) { 00300 if (i != j) { 00301 for (int k=0;k<d.dimension(0);k++) { 00302 d(k) = d(k) * ( (xout(k)-xeq(j)) / (xeq(i)-xeq(j)) ); 00303 } 00304 } 00305 } 00306 00307 // deflate end roots 00308 if ( i != 0 ) { 00309 for (int k=0;k<d.dimension(0);k++) { 00310 d(k) = -d(k) / (xeq(i) - xeq(0)); 00311 } 00312 } 00313 00314 if (i != order ) { 00315 for (int k=0;k<d.dimension(0);k++) { 00316 d(k) = d(k) / (xeq(i) - xeq(order)); 00317 } 00318 } 00319 00320 for (int k=0;k<d.dimension(0);k++) { 00321 warp(k) += d(k); 00322 } 00323 00324 } 00325 00326 00327 return; 00328 } 00329 00330 template<class Scalar, class ArrayType> 00331 void PointTools::getWarpBlendLatticeTriangle( ArrayType &points , 00332 const int order , 00333 const int offset ) 00334 { 00335 /* get Gauss-Lobatto points */ 00336 Intrepid::FieldContainer<Scalar> gaussX( order + 1 , 1 ); 00337 00338 PointTools::getWarpBlendLatticeLine<Scalar,FieldContainer<Scalar> >( gaussX , order , 0 ); 00339 00340 gaussX.resize(gaussX.dimension(0)); 00341 00342 Scalar alpopt[] = {0.0000,0.0000,1.4152,0.1001,0.2751,0.9800,1.0999, 00343 1.2832,1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258}; 00344 00345 Scalar alpha; 00346 00347 if (order >= 1 && order < 16) { 00348 alpha = alpopt[order-1]; 00349 } 00350 else { 00351 alpha = 5.0 / 3.0; 00352 } 00353 00354 const int p = order; /* switch to Warburton's notation */ 00355 int N = (p+1)*(p+2)/2; 00356 00357 /* equidistributed nodes on equilateral triangle */ 00358 Intrepid::FieldContainer<Scalar> L1( N ); 00359 Intrepid::FieldContainer<Scalar> L2( N ); 00360 Intrepid::FieldContainer<Scalar> L3( N ); 00361 Intrepid::FieldContainer<Scalar> X(N); 00362 Intrepid::FieldContainer<Scalar> Y(N); 00363 00364 int sk = 0; 00365 for (int n=1;n<=p+1;n++) { 00366 for (int m=1;m<=p+2-n;m++) { 00367 L1(sk) = (n-1) / (Scalar)p; 00368 L3(sk) = (m-1) / (Scalar)p; 00369 L2(sk) = 1.0 - L1(sk) - L3(sk); 00370 sk++; 00371 } 00372 } 00373 00374 for (int n=0;n<N;n++) { 00375 X(n) = -L2(n) + L3(n); 00376 Y(n) = (-L2(n) - L3(n) + 2*L1(n))/1.7320508075688772; 00377 } 00378 00379 /* get blending function for each node at each edge */ 00380 Intrepid::FieldContainer<Scalar> blend1(N); 00381 Intrepid::FieldContainer<Scalar> blend2(N); 00382 Intrepid::FieldContainer<Scalar> blend3(N); 00383 00384 for (int n=0;n<N;n++) { 00385 blend1(n) = 4.0 * L2(n) * L3(n); 00386 blend2(n) = 4.0 * L1(n) * L3(n); 00387 blend3(n) = 4.0 * L1(n) * L2(n); 00388 } 00389 00390 /* get difference of each barycentric coordinate */ 00391 Intrepid::FieldContainer<Scalar> L3mL2(N); 00392 Intrepid::FieldContainer<Scalar> L1mL3(N); 00393 Intrepid::FieldContainer<Scalar> L2mL1(N); 00394 00395 for (int k=0;k<N;k++) { 00396 L3mL2(k) = L3(k)-L2(k); 00397 L1mL3(k) = L1(k)-L3(k); 00398 L2mL1(k) = L2(k)-L1(k); 00399 } 00400 00401 FieldContainer<Scalar> warpfactor1(N); 00402 FieldContainer<Scalar> warpfactor2(N); 00403 FieldContainer<Scalar> warpfactor3(N); 00404 00405 warpFactor<Scalar,FieldContainer<Scalar> >( order , gaussX , L3mL2 , warpfactor1 ); 00406 warpFactor<Scalar,FieldContainer<Scalar> >( order , gaussX , L1mL3 , warpfactor2 ); 00407 warpFactor<Scalar,FieldContainer<Scalar> >( order , gaussX , L2mL1 , warpfactor3 ); 00408 00409 FieldContainer<Scalar> warp1(N); 00410 FieldContainer<Scalar> warp2(N); 00411 FieldContainer<Scalar> warp3(N); 00412 00413 for (int k=0;k<N;k++) { 00414 warp1(k) = blend1(k) * warpfactor1(k) * 00415 ( 1.0 + alpha * alpha * L1(k) * L1(k) ); 00416 warp2(k) = blend2(k) * warpfactor2(k) * 00417 ( 1.0 + alpha * alpha * L2(k) * L2(k) ); 00418 warp3(k) = blend3(k) * warpfactor3(k) * 00419 ( 1.0 + alpha * alpha * L3(k) * L3(k) ); 00420 } 00421 00422 for (int k=0;k<N;k++) { 00423 X(k) += 1.0 * warp1(k) + cos( 2.0 * M_PI / 3.0 ) * warp2(k) + cos(4*M_PI/3.0) * warp3(k); 00424 Y(k) += 0.0 * warp1(k) + sin( 2.0 * M_PI / 3.0 ) * warp2(k) + sin( 4*M_PI/3.0) * warp3(k); 00425 } 00426 00427 FieldContainer<Scalar> warXY(N,2); 00428 00429 for (int k=0;k<N;k++) { 00430 warXY(k,0) = X(k); 00431 warXY(k,1) = Y(k); 00432 } 00433 00434 00435 // finally, convert the warp-blend points to the correct triangle 00436 FieldContainer<Scalar> warburtonVerts(1,3,2); 00437 warburtonVerts(0,0,0) = -1.0; 00438 warburtonVerts(0,0,1) = -1.0/sqrt(3.0); 00439 warburtonVerts(0,1,0) = 1.0; 00440 warburtonVerts(0,1,1) = -1.0/sqrt(3.0); 00441 warburtonVerts(0,2,0) = 0.0; 00442 warburtonVerts(0,2,1) = 2.0/sqrt(3.0); 00443 00444 FieldContainer<Scalar> refPts(N,2); 00445 00446 Intrepid::CellTools<Scalar>::mapToReferenceFrame( refPts , 00447 warXY , 00448 warburtonVerts , 00449 shards::getCellTopologyData< shards::Triangle<3> >(), 00450 0 ); 00451 00452 // now write from refPts into points, taking care of offset 00453 int noffcur = 0; // index into refPts 00454 int offcur = 0; // index int points 00455 for (int i=0;i<=order;i++) { 00456 for (int j=0;j<=order-i;j++) { 00457 if ( (i >= offset) && (i <= order-offset) && 00458 (j >= offset) && (j <= order-i-offset) ) { 00459 points(offcur,0) = refPts(noffcur,0); 00460 points(offcur,1) = refPts(noffcur,1); 00461 offcur++; 00462 } 00463 noffcur++; 00464 } 00465 } 00466 00467 return; 00468 } 00469 00470 00471 template<class Scalar, class ArrayType> 00472 void PointTools::warpShiftFace3D( const int order , 00473 const Scalar pval , 00474 const ArrayType &L1, 00475 const ArrayType &L2, 00476 const ArrayType &L3, 00477 const ArrayType &L4, 00478 ArrayType &dxy) 00479 { 00480 evalshift<Scalar,ArrayType>(order,pval,L2,L3,L4,dxy); 00481 return; 00482 } 00483 00484 template<class Scalar, class ArrayType> 00485 void PointTools::evalshift( const int order , 00486 const Scalar pval , 00487 const ArrayType &L1 , 00488 const ArrayType &L2 , 00489 const ArrayType &L3 , 00490 ArrayType &dxy ) 00491 { 00492 // get Gauss-Lobatto-nodes 00493 FieldContainer<Scalar> gaussX(order+1,1); 00494 PointTools::getWarpBlendLatticeLine<Scalar,FieldContainer<Scalar> >( gaussX , order , 0 ); 00495 gaussX.resize(order+1); 00496 const int N = L1.dimension(0); 00497 00498 // Warburton code reverses them 00499 for (int k=0;k<=order;k++) { 00500 gaussX(k) *= -1.0; 00501 } 00502 00503 // blending function at each node for each edge 00504 FieldContainer<Scalar> blend1(N); 00505 FieldContainer<Scalar> blend2(N); 00506 FieldContainer<Scalar> blend3(N); 00507 00508 for (int i=0;i<N;i++) { 00509 blend1(i) = L2(i) * L3(i); 00510 blend2(i) = L1(i) * L3(i); 00511 blend3(i) = L1(i) * L2(i); 00512 } 00513 00514 // amount of warp for each node for each edge 00515 FieldContainer<Scalar> warpfactor1(N); 00516 FieldContainer<Scalar> warpfactor2(N); 00517 FieldContainer<Scalar> warpfactor3(N); 00518 00519 // differences of barycentric coordinates 00520 FieldContainer<Scalar> L3mL2(N); 00521 FieldContainer<Scalar> L1mL3(N); 00522 FieldContainer<Scalar> L2mL1(N); 00523 00524 for (int k=0;k<N;k++) { 00525 L3mL2(k) = L3(k)-L2(k); 00526 L1mL3(k) = L1(k)-L3(k); 00527 L2mL1(k) = L2(k)-L1(k); 00528 } 00529 00530 evalwarp<Scalar,FieldContainer<Scalar> >( warpfactor1 , order , gaussX , L3mL2 ); 00531 evalwarp<Scalar,FieldContainer<Scalar> >( warpfactor2 , order , gaussX , L1mL3 ); 00532 evalwarp<Scalar,FieldContainer<Scalar> >( warpfactor3 , order , gaussX , L2mL1 ); 00533 00534 for (int k=0;k<N;k++) { 00535 warpfactor1(k) *= 4.0; 00536 warpfactor2(k) *= 4.0; 00537 warpfactor3(k) *= 4.0; 00538 } 00539 00540 FieldContainer<Scalar> warp1(N); 00541 FieldContainer<Scalar> warp2(N); 00542 FieldContainer<Scalar> warp3(N); 00543 00544 for (int k=0;k<N;k++) { 00545 warp1(k) = blend1(k) * warpfactor1(k) * 00546 ( 1.0 + pval * pval * L1(k) * L1(k) ); 00547 warp2(k) = blend2(k) * warpfactor2(k) * 00548 ( 1.0 + pval * pval * L2(k) * L2(k) ); 00549 warp3(k) = blend3(k) * warpfactor3(k) * 00550 ( 1.0 + pval * pval * L3(k) * L3(k) ); 00551 } 00552 00553 for (int k=0;k<N;k++) { 00554 dxy(k,0) = 1.0 * warp1(k) + cos( 2.0 * M_PI / 3.0 ) * warp2(k) + cos( 4.0*M_PI/3.0 ) * warp3(k); 00555 dxy(k,1) = 0.0 * warp1(k) + sin( 2.0 * M_PI / 3.0 ) * warp2(k) + sin( 4.0*M_PI/3.0 ) * warp3(k); 00556 } 00557 00558 return; 00559 00560 } 00561 00562 /* one-d edge warping function */ 00563 template<class Scalar, class ArrayType> 00564 void PointTools::evalwarp( ArrayType &warp , 00565 const int order , 00566 const ArrayType &xnodes , 00567 const ArrayType &xout ) 00568 { 00569 FieldContainer<Scalar> xeq(order+1); 00570 FieldContainer<Scalar> d(xout.dimension(0)); 00571 00572 d.initialize(); 00573 00574 for (int i=0;i<=order;i++) { 00575 xeq(i) = -1.0 + 2.0 * ( order - i ) / order; 00576 } 00577 00578 00579 00580 for (int i=0;i<=order;i++) { 00581 d.initialize( xnodes(i) - xeq(i) ); 00582 for (int j=1;j<order;j++) { 00583 if (i!=j) { 00584 for (int k=0;k<d.dimension(0);k++) { 00585 d(k) = d(k) * (xout(k)-xeq(j))/(xeq(i)-xeq(j)); 00586 } 00587 } 00588 } 00589 if (i!=0) { 00590 for (int k=0;k<d.dimension(0);k++) { 00591 d(k) = -d(k)/(xeq(i)-xeq(0)); 00592 } 00593 } 00594 if (i!=order) { 00595 for (int k=0;k<d.dimension(0);k++) { 00596 d(k) = d(k)/(xeq(i)-xeq(order)); 00597 } 00598 } 00599 00600 for (int k=0;k<d.dimension(0);k++) { 00601 warp(k) += d(k); 00602 } 00603 } 00604 00605 return; 00606 } 00607 00608 00609 template<class Scalar, class ArrayType> 00610 void PointTools::getWarpBlendLatticeTetrahedron(ArrayType &points , 00611 const int order , 00612 const int offset ) 00613 { 00614 Scalar alphastore[] = { 0,0,0,0.1002, 1.1332,1.5608,1.3413,1.2577,1.1603, 00615 1.10153,0.6080,0.4523,0.8856,0.8717,0.9655}; 00616 Scalar alpha; 00617 00618 if (order <= 15) { 00619 alpha = alphastore[order-1]; 00620 } 00621 else { 00622 alpha = 1.0; 00623 } 00624 00625 const int N = (order+1)*(order+2)*(order+3)/6; 00626 Scalar tol = 1.e-10; 00627 00628 FieldContainer<Scalar> shift(N,3); 00629 shift.initialize(); 00630 00631 /* create 3d equidistributed nodes on Warburton tet */ 00632 FieldContainer<Scalar> equipoints(N,3); 00633 int sk = 0; 00634 for (int n=0;n<=order;n++) { 00635 for (int m=0;m<=order-n;m++) { 00636 for (int q=0;q<=order-n-m;q++) { 00637 equipoints(sk,0) = -1.0 + (q * 2.0 ) / order; 00638 equipoints(sk,1) = -1.0 + (m * 2.0 ) / order; 00639 equipoints(sk,2) = -1.0 + (n * 2.0 ) / order; 00640 sk++; 00641 } 00642 } 00643 } 00644 00645 00646 /* construct barycentric coordinates for equispaced lattice */ 00647 FieldContainer<Scalar> L1(N); 00648 FieldContainer<Scalar> L2(N); 00649 FieldContainer<Scalar> L3(N); 00650 FieldContainer<Scalar> L4(N); 00651 for (int i=0;i<N;i++) { 00652 L1(i) = (1.0 + equipoints(i,2)) / 2.0; 00653 L2(i) = (1.0 + equipoints(i,1)) / 2.0; 00654 L3(i) = -(1.0 + equipoints(i,0) + equipoints(i,1) + equipoints(i,2)) / 2.0; 00655 L4(i) = (1.0 + equipoints(i,0)) / 2.0; 00656 } 00657 00658 00659 /* vertices of Warburton tet */ 00660 FieldContainer<Scalar> warVerts(4,3); 00661 warVerts(0,0) = -1.0; 00662 warVerts(0,1) = -1.0/sqrt(3.0); 00663 warVerts(0,2) = -1.0/sqrt(6.0); 00664 warVerts(1,0) = 1.0; 00665 warVerts(1,1) = -1.0/sqrt(3.0); 00666 warVerts(1,2) = -1.0/sqrt(6.0); 00667 warVerts(2,0) = 0.0; 00668 warVerts(2,1) = 2.0 / sqrt(3.0); 00669 warVerts(2,2) = -1.0/sqrt(6.0); 00670 warVerts(3,0) = 0.0; 00671 warVerts(3,1) = 0.0; 00672 warVerts(3,2) = 3.0 / sqrt(6.0); 00673 00674 00675 /* tangents to faces */ 00676 FieldContainer<Scalar> t1(4,3); 00677 FieldContainer<Scalar> t2(4,3); 00678 for (int i=0;i<3;i++) { 00679 t1(0,i) = warVerts(1,i) - warVerts(0,i); 00680 t1(1,i) = warVerts(1,i) - warVerts(0,i); 00681 t1(2,i) = warVerts(2,i) - warVerts(1,i); 00682 t1(3,i) = warVerts(2,i) - warVerts(0,i); 00683 t2(0,i) = warVerts(2,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) ); 00684 t2(1,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) ); 00685 t2(2,i) = warVerts(3,i) - 0.5 * ( warVerts(1,i) + warVerts(2,i) ); 00686 t2(3,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(2,i) ); 00687 } 00688 00689 /* normalize tangents */ 00690 for (int n=0;n<4;n++) { 00691 /* Compute norm of t1(n) and t2(n) */ 00692 Scalar normt1n = 0.0; 00693 Scalar normt2n = 0.0; 00694 for (int i=0;i<3;i++) { 00695 normt1n += (t1(n,i) * t1(n,i)); 00696 normt2n += (t2(n,i) * t2(n,i)); 00697 } 00698 normt1n = sqrt(normt1n); 00699 normt2n = sqrt(normt2n); 00700 /* normalize each tangent now */ 00701 for (int i=0;i<3;i++) { 00702 t1(n,i) /= normt1n; 00703 t2(n,i) /= normt2n; 00704 } 00705 } 00706 00707 /* undeformed coordinates */ 00708 FieldContainer<Scalar> XYZ(N,3); 00709 for (int i=0;i<N;i++) { 00710 for (int j=0;j<3;j++) { 00711 XYZ(i,j) = L3(i)*warVerts(0,j) + L4(i)*warVerts(1,j) + L2(i)*warVerts(2,j) + L1(i)*warVerts(3,j); 00712 } 00713 } 00714 00715 for (int face=1;face<=4;face++) { 00716 FieldContainer<Scalar> La, Lb, Lc, Ld; 00717 FieldContainer<Scalar> warp(N,2); 00718 FieldContainer<Scalar> blend(N); 00719 FieldContainer<Scalar> denom(N); 00720 switch (face) { 00721 case 1: 00722 La = L1; Lb = L2; Lc = L3; Ld = L4; break; 00723 case 2: 00724 La = L2; Lb = L1; Lc = L3; Ld = L4; break; 00725 case 3: 00726 La = L3; Lb = L1; Lc = L4; Ld = L2; break; 00727 case 4: 00728 La = L4; Lb = L1; Lc = L3; Ld = L2; break; 00729 } 00730 00731 /* get warp tangential to face */ 00732 warpShiftFace3D<Scalar,ArrayType>(order,alpha,La,Lb,Lc,Ld,warp); 00733 00734 for (int k=0;k<N;k++) { 00735 blend(k) = Lb(k) * Lc(k) * Ld(k); 00736 } 00737 00738 for (int k=0;k<N;k++) { 00739 denom(k) = (Lb(k) + 0.5 * La(k)) * (Lc(k) + 0.5*La(k)) * (Ld(k) + 0.5 * La(k)); 00740 } 00741 00742 for (int k=0;k<N;k++) { 00743 if (denom(k) > tol) { 00744 blend(k) *= ( 1.0 + alpha * alpha * La(k) * La(k) ) / denom(k); 00745 } 00746 } 00747 00748 00749 // compute warp and blend 00750 for (int k=0;k<N;k++) { 00751 for (int j=0;j<3;j++) { 00752 shift(k,j) = shift(k,j) + blend(k) * warp(k,0) * t1(face-1,j) 00753 + blend(k) * warp(k,1) * t2(face-1,j); 00754 } 00755 } 00756 00757 for (int k=0;k<N;k++) { 00758 if (La(k) < tol && ( Lb(k) < tol || Lc(k) < tol || Ld(k) < tol )) { 00759 for (int j=0;j<3;j++) { 00760 shift(k,j) = warp(k,0) * t1(face-1,j) + warp(k,1) * t2(face-1,j); 00761 } 00762 } 00763 } 00764 00765 } 00766 00767 FieldContainer<Scalar> updatedPoints(N,3); 00768 for (int k=0;k<N;k++) { 00769 for (int j=0;j<3;j++) { 00770 updatedPoints(k,j) = XYZ(k,j) + shift(k,j); 00771 } 00772 } 00773 00774 warVerts.resize( 1 , 4 , 3 ); 00775 00776 // now we convert to Pavel's reference triangle! 00777 FieldContainer<Scalar> refPts(N,3); 00778 CellTools<Scalar>::mapToReferenceFrame( refPts ,updatedPoints , 00779 warVerts , 00780 shards::getCellTopologyData<shards::Tetrahedron<4> >() , 00781 0 ); 00782 00783 // now write from refPts into points, taking offset into account 00784 int noffcur = 0; 00785 int offcur = 0; 00786 for (int i=0;i<=order;i++) { 00787 for (int j=0;j<=order-i;j++) { 00788 for (int k=0;k<=order-i-j;k++) { 00789 if ( (i >= offset) && (i <= order-offset) && 00790 (j >= offset) && (j <= order-i-offset) && 00791 (k >= offset) && (k <= order-i-j-offset) ) { 00792 points(offcur,0) = refPts(noffcur,0); 00793 points(offcur,1) = refPts(noffcur,1); 00794 points(offcur,2) = refPts(noffcur,2); 00795 offcur++; 00796 } 00797 noffcur++; 00798 } 00799 } 00800 } 00801 00802 00803 00804 } 00805 00806 00807 } // face Intrepid
1.7.6.1