Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Cell/test_01.cpp
Go to the documentation of this file.
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