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