|
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 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or 00025 // Denis Ridzal (dridzal@sandia.gov). 00026 // 00027 // ************************************************************************ 00028 // @HEADER 00029 00030 00035 #include "Intrepid_CellTools.hpp" 00036 #include "Intrepid_FieldContainer.hpp" 00037 #include "Intrepid_DefaultCubatureFactory.hpp" 00038 #include "Teuchos_oblackholestream.hpp" 00039 #include "Teuchos_RCP.hpp" 00040 #include "Teuchos_GlobalMPISession.hpp" 00041 #include "Teuchos_ScalarTraits.hpp" 00042 00043 using namespace std; 00044 using namespace Intrepid; 00045 00046 00065 void testSubcellParametrizations(int& errorFlag, 00066 const shards::CellTopology& parentCell, 00067 const FieldContainer<double>& subcParamVert_A, 00068 const FieldContainer<double>& subcParamVert_B, 00069 const int subcDim, 00070 const Teuchos::RCP<std::ostream>& outStream); 00071 00072 00073 int main(int argc, char *argv[]) { 00074 00075 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00076 00077 typedef CellTools<double> CellTools; 00078 typedef shards::CellTopology CellTopology; 00079 00080 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided. 00081 int iprint = argc - 1; 00082 00083 Teuchos::RCP<std::ostream> outStream; 00084 Teuchos::oblackholestream bhs; // outputs nothing 00085 00086 if (iprint > 0) 00087 outStream = Teuchos::rcp(&std::cout, false); 00088 else 00089 outStream = Teuchos::rcp(&bhs, false); 00090 00091 // Save the format state of the original std::cout. 00092 Teuchos::oblackholestream oldFormatState; 00093 oldFormatState.copyfmt(std::cout); 00094 00095 *outStream \ 00096 << "===============================================================================\n" \ 00097 << "| |\n" \ 00098 << "| Unit Test CellTools |\n" \ 00099 << "| |\n" \ 00100 << "| 1) Edge parametrizations |\n" \ 00101 << "| 2) Face parametrizations |\n" \ 00102 << "| 3) Edge tangents |\n" \ 00103 << "| 4) Face tangents and normals |\n" \ 00104 << "| |\n" \ 00105 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \ 00106 << "| Denis Ridzal (dridzal@sandia.gov), or |\n" \ 00107 << "| Kara Peterson (kjpeter@sandia.gov) |\n" \ 00108 << "| |\n" \ 00109 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00110 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00111 << "| |\n" \ 00112 << "===============================================================================\n"; 00113 00114 int errorFlag = 0; 00115 00116 00117 // Vertices of the parametrization domain for 1-subcells: standard 1-cube [-1,1] 00118 FieldContainer<double> cube_1(2, 1); 00119 cube_1(0,0) = -1.0; 00120 cube_1(1,0) = 1.0; 00121 00122 00123 // Vertices of the parametrization domain for triangular faces: the standard 2-simplex 00124 FieldContainer<double> simplex_2(3, 2); 00125 simplex_2(0, 0) = 0.0; simplex_2(0, 1) = 0.0; 00126 simplex_2(1, 0) = 1.0; simplex_2(1, 1) = 0.0; 00127 simplex_2(2, 0) = 0.0; simplex_2(2, 1) = 1.0; 00128 00129 00130 // Vertices of the parametrization domain for quadrilateral faces: the standard 2-cube 00131 FieldContainer<double> cube_2(4, 2); 00132 cube_2(0, 0) = -1.0; cube_2(0, 1) = -1.0; 00133 cube_2(1, 0) = 1.0; cube_2(1, 1) = -1.0; 00134 cube_2(2, 0) = 1.0; cube_2(2, 1) = 1.0; 00135 cube_2(3, 0) = -1.0; cube_2(3, 1) = 1.0; 00136 00137 00138 // Pull all available topologies from Shards 00139 std::vector<shards::CellTopology> allTopologies; 00140 shards::getTopologies(allTopologies); 00141 00142 00143 // Set to 1 for edge and 2 for face tests 00144 int subcDim; 00145 00146 try{ 00147 00148 *outStream \ 00149 << "\n" 00150 << "===============================================================================\n"\ 00151 << "| Test 1: edge parametrizations: |\n"\ 00152 << "===============================================================================\n\n"; 00153 00154 subcDim = 1; 00155 00156 // Loop over the cell topologies 00157 for(int topoOrd = 0; topoOrd < (int)allTopologies.size(); topoOrd++){ 00158 00159 // Test only 2D and 3D topologies that have reference cells, e.g., exclude Line, Pentagon, etc. 00160 if(allTopologies[topoOrd].getDimension() > 1 && CellTools::hasReferenceCell(allTopologies[topoOrd]) ){ 00161 *outStream << " Testing edge parametrization for " << allTopologies[topoOrd].getName() <<"\n"; 00162 testSubcellParametrizations(errorFlag, 00163 allTopologies[topoOrd], 00164 cube_1, 00165 cube_1, 00166 subcDim, 00167 outStream); 00168 } 00169 } 00170 00171 00172 *outStream \ 00173 << "\n" 00174 << "===============================================================================\n"\ 00175 << "| Test 2: face parametrizations: |\n"\ 00176 << "===============================================================================\n\n"; 00177 00178 subcDim = 2; 00179 00180 // Loop over the cell topologies 00181 for(int topoOrd = 0; topoOrd < (int)allTopologies.size(); topoOrd++){ 00182 00183 // Test only 3D topologies that have reference cells 00184 if(allTopologies[topoOrd].getDimension() > 2 && CellTools::hasReferenceCell(allTopologies[topoOrd]) ){ 00185 *outStream << " Testing face parametrization for cell topology " << allTopologies[topoOrd].getName() <<"\n"; 00186 testSubcellParametrizations(errorFlag, 00187 allTopologies[topoOrd], 00188 simplex_2, 00189 cube_2, 00190 subcDim, 00191 outStream); 00192 } 00193 } 00194 00195 00196 /*********************************************************************************************** 00197 * 00198 * Common for test 3 and 4: edge tangents and face normals for standard cells with base topo 00199 * 00200 **********************************************************************************************/ 00201 00202 // Allocate storage and extract all standard cells with base topologies 00203 std::vector<shards::CellTopology> standardBaseTopologies; 00204 shards::getTopologies(standardBaseTopologies, 4, shards::STANDARD_CELL, shards::BASE_TOPOLOGY); 00205 00206 // Define topologies for the edge and face parametrization domains. (faces are Tri or Quad) 00207 CellTopology paramEdge (shards::getCellTopologyData<shards::Line<2> >() ); 00208 CellTopology paramTriFace (shards::getCellTopologyData<shards::Triangle<3> >() ); 00209 CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() ); 00210 00211 // Define CubatureFactory: 00212 DefaultCubatureFactory<double> cubFactory; 00213 00214 *outStream \ 00215 << "\n" 00216 << "===============================================================================\n"\ 00217 << "| Test 3: edge tangents/normals for stand. cells with base topologies: |\n"\ 00218 << "===============================================================================\n\n"; 00219 // This test loops over standard cells with base topologies, creates a set of nodes and tests tangents/normals 00220 std::vector<shards::CellTopology>::iterator cti; 00221 00222 // Define cubature on the edge parametrization domain: 00223 Teuchos::RCP<Cubature<double> > edgeCubature = cubFactory.create(paramEdge, 6); 00224 int cubDim = edgeCubature -> getDimension(); 00225 int numCubPoints = edgeCubature -> getNumPoints(); 00226 00227 // Allocate storage for cubature points and weights on edge parameter domain and fill with points: 00228 FieldContainer<double> paramEdgePoints(numCubPoints, cubDim); 00229 FieldContainer<double> paramEdgeWeights(numCubPoints); 00230 edgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights); 00231 00232 00233 // Loop over admissible topologies 00234 for(cti = standardBaseTopologies.begin(); cti !=standardBaseTopologies.end(); ++cti){ 00235 00236 // Exclude 0D (node), 1D (Line) and Pyramid<5> cells 00237 if( ( (*cti).getDimension() >= 2) && ( (*cti).getKey() != shards::Pyramid<5>::key) ){ 00238 00239 int cellDim = (*cti).getDimension(); 00240 int vCount = (*cti).getVertexCount(); 00241 FieldContainer<double> refCellVertices(vCount, cellDim); 00242 CellTools::getReferenceSubcellVertices(refCellVertices, cellDim, 0, (*cti) ); 00243 00244 *outStream << " Testing edge tangents"; 00245 if(cellDim == 2) { *outStream << " and normals"; } 00246 *outStream <<" for cell topology " << (*cti).getName() <<"\n"; 00247 00248 00249 // Array for physical cell vertices ( must have rank 3 for setJacobians) 00250 FieldContainer<double> physCellVertices(1, vCount, cellDim); 00251 00252 // Randomize reference cell vertices by moving them up to +/- (1/8) units along their 00253 // coordinate axis. Guaranteed to be non-degenerate for standard cells with base topology 00254 for(int v = 0; v < vCount; v++){ 00255 for(int d = 0; d < cellDim; d++){ 00256 double delta = Teuchos::ScalarTraits<double>::random()/8.0; 00257 physCellVertices(0, v, d) = refCellVertices(v, d) + delta; 00258 } //for d 00259 }// for v 00260 00261 // Allocate storage for cub. points on a ref. edge; Jacobians, phys. edge tangents/normals 00262 FieldContainer<double> refEdgePoints(numCubPoints, cellDim); 00263 FieldContainer<double> edgePointsJacobians(1, numCubPoints, cellDim, cellDim); 00264 FieldContainer<double> edgePointTangents(1, numCubPoints, cellDim); 00265 FieldContainer<double> edgePointNormals(1, numCubPoints, cellDim); 00266 00267 // Loop over edges: 00268 for(int edgeOrd = 0; edgeOrd < (int)(*cti).getEdgeCount(); edgeOrd++){ 00269 /* 00270 * Compute tangents on the specified physical edge using CellTools: 00271 * 1. Map points from edge parametrization domain to ref. edge with specified ordinal 00272 * 2. Compute parent cell Jacobians at ref. edge points 00273 * 3. Compute physical edge tangents 00274 */ 00275 CellTools::mapToReferenceSubcell(refEdgePoints, paramEdgePoints, 1, edgeOrd, (*cti) ); 00276 CellTools::setJacobian(edgePointsJacobians, refEdgePoints, physCellVertices, (*cti) ); 00277 CellTools::getPhysicalEdgeTangents(edgePointTangents, edgePointsJacobians, edgeOrd, (*cti)); 00278 /* 00279 * Compute tangents directly using parametrization of phys. edge and compare with CellTools tangents. 00280 * 1. Get edge vertices 00281 * 2. For affine edges tangent coordinates are given by F'(t) = (V1-V0)/2 00282 * (for now we only test affine edges, but later we will test edges for cells 00283 * with extended topologies.) 00284 */ 00285 int v0ord = (*cti).getNodeMap(1, edgeOrd, 0); 00286 int v1ord = (*cti).getNodeMap(1, edgeOrd, 1); 00287 00288 for(int pt = 0; pt < numCubPoints; pt++){ 00289 00290 // Temp storage for directly computed edge tangents 00291 FieldContainer<double> edgeBenchmarkTangents(3); 00292 00293 for(int d = 0; d < cellDim; d++){ 00294 edgeBenchmarkTangents(d) = (physCellVertices(0, v1ord, d) - physCellVertices(0, v0ord, d))/2.0; 00295 00296 // Compare with d-component of edge tangent by CellTools 00297 if( abs(edgeBenchmarkTangents(d) - edgePointTangents(0, pt, d)) > INTREPID_THRESHOLD ){ 00298 errorFlag++; 00299 *outStream 00300 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00301 << " Edge tangent computation by CellTools failed for: \n" 00302 << " Cell Topology = " << (*cti).getName() << "\n" 00303 << " Edge ordinal = " << edgeOrd << "\n" 00304 << " Edge point number = " << pt << "\n" 00305 << " Tangent coordinate = " << d << "\n" 00306 << " CellTools value = " << edgePointTangents(0, pt, d) << "\n" 00307 << " Benchmark value = " << edgeBenchmarkTangents(d) << "\n\n"; 00308 } 00309 } // for d 00310 00311 // Test side normals for 2D cells only: edge normal has coordinates (t1, -t0) 00312 if(cellDim == 2) { 00313 CellTools::getPhysicalSideNormals(edgePointNormals, edgePointsJacobians, edgeOrd, (*cti)); 00314 if( abs(edgeBenchmarkTangents(1) - edgePointNormals(0, pt, 0)) > INTREPID_THRESHOLD ){ 00315 errorFlag++; 00316 *outStream 00317 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00318 << " Edge Normal computation by CellTools failed for: \n" 00319 << " Cell Topology = " << (*cti).getName() << "\n" 00320 << " Edge ordinal = " << edgeOrd << "\n" 00321 << " Edge point number = " << pt << "\n" 00322 << " Normal coordinate = " << 0 << "\n" 00323 << " CellTools value = " << edgePointNormals(0, pt, 0) << "\n" 00324 << " Benchmark value = " << edgeBenchmarkTangents(1) << "\n\n"; 00325 } 00326 if( abs(edgeBenchmarkTangents(0) + edgePointNormals(0, pt, 1)) > INTREPID_THRESHOLD ){ 00327 errorFlag++; 00328 *outStream 00329 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00330 << " Edge Normal computation by CellTools failed for: \n" 00331 << " Cell Topology = " << (*cti).getName() << "\n" 00332 << " Edge ordinal = " << edgeOrd << "\n" 00333 << " Edge point number = " << pt << "\n" 00334 << " Normal coordinate = " << 1 << "\n" 00335 << " CellTools value = " << edgePointNormals(0, pt, 1) << "\n" 00336 << " Benchmark value = " << -edgeBenchmarkTangents(0) << "\n\n"; 00337 } 00338 } // edge normals 00339 } // for pt 00340 }// for edgeOrd 00341 }// if admissible cell 00342 }// for cti 00343 00344 00345 00346 *outStream \ 00347 << "\n" 00348 << "===============================================================================\n"\ 00349 << "| Test 4: face/side normals for stand. 3D cells with base topologies: | |\n"\ 00350 << "===============================================================================\n\n"; 00351 // This test loops over standard 3D cells with base topologies, creates a set of nodes and tests normals 00352 00353 // Define cubature on the edge parametrization domain: 00354 Teuchos::RCP<Cubature<double> > triFaceCubature = cubFactory.create(paramTriFace, 6); 00355 Teuchos::RCP<Cubature<double> > quadFaceCubature = cubFactory.create(paramQuadFace, 6); 00356 00357 int faceCubDim = triFaceCubature -> getDimension(); 00358 int numTriFaceCubPoints = triFaceCubature -> getNumPoints(); 00359 int numQuadFaceCubPoints = quadFaceCubature -> getNumPoints(); 00360 00361 // Allocate storage for cubature points and weights on face parameter domain and fill with points: 00362 FieldContainer<double> paramTriFacePoints(numTriFaceCubPoints, faceCubDim); 00363 FieldContainer<double> paramTriFaceWeights(numTriFaceCubPoints); 00364 FieldContainer<double> paramQuadFacePoints(numQuadFaceCubPoints, faceCubDim); 00365 FieldContainer<double> paramQuadFaceWeights(numQuadFaceCubPoints); 00366 00367 triFaceCubature -> getCubature(paramTriFacePoints, paramTriFaceWeights); 00368 quadFaceCubature -> getCubature(paramQuadFacePoints, paramQuadFaceWeights); 00369 00370 00371 // Loop over admissible topologies 00372 for(cti = standardBaseTopologies.begin(); cti !=standardBaseTopologies.end(); ++cti){ 00373 00374 // Exclude 2D and Pyramid<5> cells 00375 if( ( (*cti).getDimension() == 3) && ( (*cti).getKey() != shards::Pyramid<5>::key) ){ 00376 00377 int cellDim = (*cti).getDimension(); 00378 int vCount = (*cti).getVertexCount(); 00379 FieldContainer<double> refCellVertices(vCount, cellDim); 00380 CellTools::getReferenceSubcellVertices(refCellVertices, cellDim, 0, (*cti) ); 00381 00382 *outStream << " Testing face/side normals for cell topology " << (*cti).getName() <<"\n"; 00383 00384 // Array for physical cell vertices ( must have rank 3 for setJacobians) 00385 FieldContainer<double> physCellVertices(1, vCount, cellDim); 00386 00387 // Randomize reference cell vertices by moving them up to +/- (1/8) units along their 00388 // coordinate axis. Guaranteed to be non-degenerate for standard cells with base topology 00389 for(int v = 0; v < vCount; v++){ 00390 for(int d = 0; d < cellDim; d++){ 00391 double delta = Teuchos::ScalarTraits<double>::random()/8.0; 00392 physCellVertices(0, v, d) = refCellVertices(v, d) + delta; 00393 } //for d 00394 }// for v 00395 00396 // Allocate storage for cub. points on a ref. face; Jacobians, phys. face normals and 00397 // benchmark normals. 00398 FieldContainer<double> refTriFacePoints(numTriFaceCubPoints, cellDim); 00399 FieldContainer<double> refQuadFacePoints(numQuadFaceCubPoints, cellDim); 00400 FieldContainer<double> triFacePointsJacobians(1, numTriFaceCubPoints, cellDim, cellDim); 00401 FieldContainer<double> quadFacePointsJacobians(1, numQuadFaceCubPoints, cellDim, cellDim); 00402 FieldContainer<double> triFacePointNormals(1, numTriFaceCubPoints, cellDim); 00403 FieldContainer<double> triSidePointNormals(1, numTriFaceCubPoints, cellDim); 00404 FieldContainer<double> quadFacePointNormals(1, numQuadFaceCubPoints, cellDim); 00405 FieldContainer<double> quadSidePointNormals(1, numQuadFaceCubPoints, cellDim); 00406 00407 00408 // Loop over faces: 00409 for(int faceOrd = 0; faceOrd < (int)(*cti).getSideCount(); faceOrd++){ 00410 00411 // This test presently includes only Triangle<3> and Quadrilateral<4> faces. Once we support 00412 // cells with extended topologies we will add their faces as well. 00413 switch( (*cti).getTopology(2, faceOrd) -> key ) { 00414 00415 case shards::Triangle<3>::key: 00416 { 00417 // Compute face normals using CellTools 00418 CellTools::mapToReferenceSubcell(refTriFacePoints, paramTriFacePoints, 2, faceOrd, (*cti) ); 00419 CellTools::setJacobian(triFacePointsJacobians, refTriFacePoints, physCellVertices, (*cti) ); 00420 CellTools::getPhysicalFaceNormals(triFacePointNormals, triFacePointsJacobians, faceOrd, (*cti)); 00421 CellTools::getPhysicalSideNormals(triSidePointNormals, triFacePointsJacobians, faceOrd, (*cti)); 00422 /* 00423 * Compute face normals using direct linear parametrization of the face: the map from 00424 * standard 2-simplex to physical Triangle<3> face in 3D is 00425 * F(x,y) = V0 + (V1-V0)x + (V2-V0)*y 00426 * Face normal is vector product Tx X Ty where Tx = (V1-V0); Ty = (V2-V0) 00427 */ 00428 int v0ord = (*cti).getNodeMap(2, faceOrd, 0); 00429 int v1ord = (*cti).getNodeMap(2, faceOrd, 1); 00430 int v2ord = (*cti).getNodeMap(2, faceOrd, 2); 00431 00432 // Loop over face points: redundant for affine faces, but CellTools gives one vector 00433 // per point so need to check all points anyways. 00434 for(int pt = 0; pt < numTriFaceCubPoints; pt++){ 00435 FieldContainer<double> tanX(3), tanY(3), faceNormal(3); 00436 for(int d = 0; d < cellDim; d++){ 00437 tanX(d) = (physCellVertices(0, v1ord, d) - physCellVertices(0, v0ord, d)); 00438 tanY(d) = (physCellVertices(0, v2ord, d) - physCellVertices(0, v0ord, d)); 00439 }// for d 00440 00441 RealSpaceTools<double>::vecprod(faceNormal, tanX, tanY); 00442 00443 // Compare direct normal with d-component of the face/side normal by CellTools 00444 for(int d = 0; d < cellDim; d++){ 00445 00446 // face normal method 00447 if( abs(faceNormal(d) - triFacePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){ 00448 errorFlag++; 00449 *outStream 00450 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00451 << " Face normal computation by CellTools failed for: \n" 00452 << " Cell Topology = " << (*cti).getName() << "\n" 00453 << " Face Topology = " << (*cti).getTopology(2, faceOrd) -> name << "\n" 00454 << " Face ordinal = " << faceOrd << "\n" 00455 << " Face point number = " << pt << "\n" 00456 << " Normal coordinate = " << d << "\n" 00457 << " CellTools value = " << triFacePointNormals(0, pt, d) 00458 << " Benchmark value = " << faceNormal(d) << "\n\n"; 00459 } 00460 //side normal method 00461 if( abs(faceNormal(d) - triSidePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){ 00462 errorFlag++; 00463 *outStream 00464 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00465 << " Side normal computation by CellTools failed for: \n" 00466 << " Cell Topology = " << (*cti).getName() << "\n" 00467 << " Side Topology = " << (*cti).getTopology(2, faceOrd) -> name << "\n" 00468 << " Side ordinal = " << faceOrd << "\n" 00469 << " Side point number = " << pt << "\n" 00470 << " Normal coordinate = " << d << "\n" 00471 << " CellTools value = " << triSidePointNormals(0, pt, d) 00472 << " Benchmark value = " << faceNormal(d) << "\n\n"; 00473 } 00474 } // for d 00475 } // for pt 00476 } 00477 break; 00478 00479 case shards::Quadrilateral<4>::key: 00480 { 00481 // Compute face normals using CellTools 00482 CellTools::mapToReferenceSubcell(refQuadFacePoints, paramQuadFacePoints, 2, faceOrd, (*cti) ); 00483 CellTools::setJacobian(quadFacePointsJacobians, refQuadFacePoints, physCellVertices, (*cti) ); 00484 CellTools::getPhysicalFaceNormals(quadFacePointNormals, quadFacePointsJacobians, faceOrd, (*cti)); 00485 CellTools::getPhysicalSideNormals(quadSidePointNormals, quadFacePointsJacobians, faceOrd, (*cti)); 00486 /* 00487 * Compute face normals using direct bilinear parametrization of the face: the map from 00488 * [-1,1]^2 to physical Quadrilateral<4> face in 3D is 00489 * F(x,y) = ((V0+V1+V2+V3) + (-V0+V1+V2-V3)*X + (-V0-V1+V2+V3)*Y + (V0-V1+V2-V3)*X*Y)/4 00490 * Face normal is vector product Tx X Ty where 00491 * Tx = ((-V0+V1+V2-V3) + (V0-V1+V2-V3)*Y)/4 00492 * Ty = ((-V0-V1+V2+V3) + (V0-V1+V2-V3)*X)/4 00493 */ 00494 int v0ord = (*cti).getNodeMap(2, faceOrd, 0); 00495 int v1ord = (*cti).getNodeMap(2, faceOrd, 1); 00496 int v2ord = (*cti).getNodeMap(2, faceOrd, 2); 00497 int v3ord = (*cti).getNodeMap(2, faceOrd, 3); 00498 00499 // Loop over face points (redundant for affine faces, but needed for later when we handle non-affine ones) 00500 for(int pt = 0; pt < numTriFaceCubPoints; pt++){ 00501 FieldContainer<double> tanX(3), tanY(3), faceNormal(3); 00502 for(int d = 0; d < cellDim; d++){ 00503 tanX(d) = (physCellVertices(0, v0ord, d)*(-1.0 + paramQuadFacePoints(pt,1) ) + 00504 physCellVertices(0, v1ord, d)*( 1.0 - paramQuadFacePoints(pt,1) ) + 00505 physCellVertices(0, v2ord, d)*( 1.0 + paramQuadFacePoints(pt,1) ) + 00506 physCellVertices(0, v3ord, d)*(-1.0 - paramQuadFacePoints(pt,1) ) )/4.0; 00507 00508 tanY(d) = (physCellVertices(0, v0ord, d)*(-1.0 + paramQuadFacePoints(pt,0) ) + 00509 physCellVertices(0, v1ord, d)*(-1.0 - paramQuadFacePoints(pt,0) ) + 00510 physCellVertices(0, v2ord, d)*( 1.0 + paramQuadFacePoints(pt,0) ) + 00511 physCellVertices(0, v3ord, d)*( 1.0 - paramQuadFacePoints(pt,0) ) )/4.0; 00512 }// for d 00513 00514 RealSpaceTools<double>::vecprod(faceNormal, tanX, tanY); 00515 // Compare direct normal with d-component of the face/side normal by CellTools 00516 for(int d = 0; d < cellDim; d++){ 00517 00518 // face normal method 00519 if( abs(faceNormal(d) - quadFacePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){ 00520 errorFlag++; 00521 *outStream 00522 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00523 << " Face normal computation by CellTools failed for: \n" 00524 << " Cell Topology = " << (*cti).getName() << "\n" 00525 << " Face Topology = " << (*cti).getTopology(2, faceOrd) -> name << "\n" 00526 << " Face ordinal = " << faceOrd << "\n" 00527 << " Face point number = " << pt << "\n" 00528 << " Normal coordinate = " << d << "\n" 00529 << " CellTools value = " << quadFacePointNormals(0, pt, d) 00530 << " Benchmark value = " << faceNormal(d) << "\n\n"; 00531 } 00532 //side normal method 00533 if( abs(faceNormal(d) - quadSidePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){ 00534 errorFlag++; 00535 *outStream 00536 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00537 << " Side normal computation by CellTools failed for: \n" 00538 << " Cell Topology = " << (*cti).getName() << "\n" 00539 << " Side Topology = " << (*cti).getTopology(2, faceOrd) -> name << "\n" 00540 << " Side ordinal = " << faceOrd << "\n" 00541 << " Side point number = " << pt << "\n" 00542 << " Normal coordinate = " << d << "\n" 00543 << " CellTools value = " << quadSidePointNormals(0, pt, d) 00544 << " Benchmark value = " << faceNormal(d) << "\n\n"; 00545 } 00546 } // for d 00547 }// for pt 00548 }// case Quad 00549 break; 00550 default: 00551 errorFlag++; 00552 *outStream << " Face normals test failure: face topology not supported \n\n"; 00553 } // switch 00554 }// for faceOrd 00555 }// if admissible 00556 }// for cti 00557 }// try 00558 00559 //============================================================================================// 00560 // Wrap up test: check if the test broke down unexpectedly due to an exception // 00561 //============================================================================================// 00562 catch (std::logic_error err) { 00563 *outStream << err.what() << "\n"; 00564 errorFlag = -1000; 00565 }; 00566 00567 00568 if (errorFlag != 0) 00569 std::cout << "End Result: TEST FAILED\n"; 00570 else 00571 std::cout << "End Result: TEST PASSED\n"; 00572 00573 // reset format state of std::cout 00574 std::cout.copyfmt(oldFormatState); 00575 00576 return errorFlag; 00577 } 00578 00579 00580 00581 void testSubcellParametrizations(int& errorFlag, 00582 const shards::CellTopology& parentCell, 00583 const FieldContainer<double>& subcParamVert_A, 00584 const FieldContainer<double>& subcParamVert_B, 00585 const int subcDim, 00586 const Teuchos::RCP<std::ostream>& outStream){ 00587 00588 // Get cell dimension and subcell count 00589 int cellDim = parentCell.getDimension(); 00590 int subcCount = parentCell.getSubcellCount(subcDim); 00591 00592 00593 // Loop over subcells of the specified dimension 00594 for(int subcOrd = 0; subcOrd < subcCount; subcOrd++){ 00595 int subcVertexCount = parentCell.getVertexCount(subcDim, subcOrd); 00596 00597 00598 // Storage for correct reference subcell vertices and for the images of the parametrization domain points 00599 FieldContainer<double> refSubcellVertices(subcVertexCount, cellDim); 00600 FieldContainer<double> mappedParamVertices(subcVertexCount, cellDim); 00601 00602 00603 // Retrieve correct reference subcell vertices 00604 CellTools<double>::getReferenceSubcellVertices(refSubcellVertices, subcDim, subcOrd, parentCell); 00605 00606 00607 // Map vertices of the parametrization domain to 1 or 2-subcell with ordinal subcOrd 00608 // For edges parametrization domain is always 1-cube passed as "subcParamVert_A" 00609 if(subcDim == 1) { 00610 CellTools<double>::mapToReferenceSubcell(mappedParamVertices, 00611 subcParamVert_A, 00612 subcDim, 00613 subcOrd, 00614 parentCell); 00615 } 00616 // For faces need to treat Triangle and Quadrilateral faces separately 00617 else if(subcDim == 2) { 00618 00619 // domain "subcParamVert_A" is the standard 2-simplex 00620 if(subcVertexCount == 3){ 00621 CellTools<double>::mapToReferenceSubcell(mappedParamVertices, 00622 subcParamVert_A, 00623 subcDim, 00624 subcOrd, 00625 parentCell); 00626 } 00627 // Domain "subcParamVert_B" is the standard 2-cube 00628 else if(subcVertexCount == 4){ 00629 CellTools<double>::mapToReferenceSubcell(mappedParamVertices, 00630 subcParamVert_B, 00631 subcDim, 00632 subcOrd, 00633 parentCell); 00634 } 00635 } 00636 00637 // Compare the images of the parametrization domain vertices with the true vertices. 00638 for(int subcVertOrd = 0; subcVertOrd < subcVertexCount; subcVertOrd++){ 00639 for(int dim = 0; dim < cellDim; dim++){ 00640 00641 if(mappedParamVertices(subcVertOrd, dim) != refSubcellVertices(subcVertOrd, dim) ) { 00642 errorFlag++; 00643 *outStream 00644 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00645 << " Cell Topology = " << parentCell.getName() << "\n" 00646 << " Parametrization of subcell " << subcOrd << " which is " 00647 << parentCell.getName(subcDim,subcOrd) << " failed for vertex " << subcVertOrd << ":\n" 00648 << " parametrization map fails to map correctly coordinate " << dim << " of that vertex\n\n"; 00649 00650 }//if 00651 }// for dim 00652 }// for subcVertOrd 00653 }// for subcOrd 00654 00655 } 00656 00657 00658 00659 00660 00661
1.7.4