Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/test/Cell/test

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
http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/test/Cell/test_01.cpp
