Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/example/CellTools/example_02.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 
00050 #include "Intrepid_CellTools.hpp"
00051 #include "Intrepid_FieldContainer.hpp"
00052 #include "Intrepid_DefaultCubatureFactory.hpp"
00053 #include "Teuchos_GlobalMPISession.hpp"
00054 
00055 #include "Shards_CellTopology.hpp"
00056 
00057 #include "Teuchos_RCP.hpp"
00058 
00059 using namespace std;
00060 using namespace Intrepid;
00061 using namespace shards;
00062 
00063 int main(int argc, char *argv[]) {
00064 
00065   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00066 
00067   typedef CellTools<double>       CellTools;
00068   typedef shards::CellTopology    CellTopology;
00069   
00070   cout \
00071   << "===============================================================================\n" \
00072   << "|                                                                             |\n" \
00073   << "|                   Example use of the CellTools class                        |\n" \
00074   << "|                                                                             |\n" \
00075   << "|  1) Reference edge parametrizations                                         |\n" \
00076   << "|  2) Reference face parametrizations                                         |\n" \
00077   << "|                                                                             |\n" \
00078   << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov)                      |\n" \
00079   << "|                      Denis Ridzal (dridzal@sandia.gov), or                  |\n" \
00080   << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00081   << "|                                                                             |\n" \
00082   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00083   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00084   << "|                                                                             |\n" \
00085   << "===============================================================================\n"\
00086   << "| Summary:                                                                    |\n"\
00087   << "| Reference edge parametrizations map [-1,1] to the edges of reference cells. |\n"\
00088   << "| They are used to define, e.g., integration points on the edges of 2D and 3D |\n"\
00089   << "| reference cells. Edge parametrizations for special 2D  cells such as Beam   |\n"\
00090   << "| and ShellLine, are also supported.                                          |\n"\
00091   << "===============================================================================\n";
00092  
00093   /* 
00094     Specification of integration points on 1-subcells (edges) of reference cells. Edges are 
00095     parametrized by [-1,1] and integration points on an edge are defined by mapping integration
00096     points from the parametrization domain [-1,1] to a specific edge on the reference cell.
00097    
00098    1. Common tasks: definition of integration points in the edge parametrization domain [-1,1]
00099       These steps are independent of parent cell topology:
00100    
00101       a. Instantiate a CubatureFactory object to create cubatures (needed for face maps too)
00102       b. Define parametrization domain for the edges as having Line<2> cell topology. This is 
00103          required by the CubatureFactory in order to select cubature points and weights from 
00104          the reference line [-1,1]
00105       c. Use CubatureFactory to select cubature of the desired degree for the Line<2> topology
00106       d. Allocate containers for the cubature points and weights.
00107    
00108    2. Parent cell topology specific tasks
00109    
00110       a. Select the parent cell topology
00111       b. Allocate containers for the images of the integration points on [-1,1] on the edges
00112       c. Apply the edge parametrization map to the pointss in [-1,1]
00113    */
00114   
00115   // Step 1.a (Define CubatureFactory)
00116   DefaultCubatureFactory<double>  cubFactory;   
00117   
00118   
00119   // Step 1.b (Define the topology of the parametrization domain)
00120   CellTopology edgeParam(shards::getCellTopologyData<shards::Line<2> >() );
00121   
00122   
00123   // Step 1.c (selects Gauss rule of order 6 on [-1,1]) 
00124   Teuchos::RCP<Cubature<double> > edgeParamCubature = cubFactory.create(edgeParam, 6); 
00125   
00126   
00127   // Step 1.d (allocate storage for cubature points)
00128   int cubDim       = edgeParamCubature -> getDimension();
00129   int numCubPoints = edgeParamCubature -> getNumPoints();
00130   
00131   FieldContainer<double> edgeParamCubPts(numCubPoints, cubDim);
00132   FieldContainer<double> edgeParamCubWts(numCubPoints);
00133   edgeParamCubature -> getCubature(edgeParamCubPts, edgeParamCubWts);
00134     
00135   
00136   
00137   std::cout \
00138   << "===============================================================================\n"\
00139   << "| EXAMPLE 1.1                                                                 |\n"
00140   << "| Edge parametrizations for standard 2D cells: Triangle                       |\n"\
00141   << "===============================================================================\n";
00142     
00143   // Step 2.a (select reference cell topology)
00144   CellTopology triangle_3(getCellTopologyData<Triangle<3> >() );
00145   
00146   
00147   // Step 2.b (allocate storage for points on an edge of the reference cell)
00148   FieldContainer<double> triEdgePoints(numCubPoints, triangle_3.getDimension() );
00149     
00150   
00151   // Step 2.c (same points are mapped to all edges, can also map different set to each edge) 
00152   for(int edgeOrd = 0; edgeOrd < (int)triangle_3.getEdgeCount(); edgeOrd++){
00153     
00154     CellTools::mapToReferenceSubcell(triEdgePoints, 
00155                                      edgeParamCubPts,
00156                                      1,
00157                                      edgeOrd,
00158                                      triangle_3);
00159     
00160     // Optional: print the vertices of the reference edge
00161     CellTools::printSubcellVertices(1, edgeOrd, triangle_3);
00162         
00163     for(int pt = 0; pt < numCubPoints; pt++){
00164       std::cout << "\t Parameter point " 
00165       << std::setw(12) << std::right << edgeParamCubPts(pt, 0) << std::setw(10) << "  -->  " << "(" 
00166       << std::setw(10) << std::right << triEdgePoints(pt, 0) << " , " 
00167       << std::setw(10) << std::right << triEdgePoints(pt, 1) << ")\n";
00168     }    
00169     std::cout << "\n";
00170   }
00171   
00172   
00173   
00174   std::cout \
00175     << "===============================================================================\n"\
00176     << "| EXAMPLE 1.2                                                                 |\n"
00177     << "| Edge parametrizations for standard 2D cells: Quadrilateral                  |\n"\
00178     << "===============================================================================\n";
00179   
00180   // Step 2.a (select reference cell topology)
00181   CellTopology quad_4(getCellTopologyData<Quadrilateral<4> >() );
00182   
00183   
00184   // Step 2.b (allocate storage for points on an edge of the reference cell)
00185   FieldContainer<double> quadEdgePoints(numCubPoints, quad_4.getDimension() );
00186   
00187   
00188   // Step 2.c (same points are mapped to all edges, can also map different set to each edge) 
00189   for(int edgeOrd = 0; edgeOrd < (int)quad_4.getEdgeCount(); edgeOrd++){
00190     
00191     CellTools::mapToReferenceSubcell(quadEdgePoints, 
00192                                      edgeParamCubPts,
00193                                      1,
00194                                      edgeOrd,
00195                                      quad_4);
00196     
00197     // Optional: print the vertices of the reference edge
00198     CellTools::printSubcellVertices(1, edgeOrd, quad_4);
00199     
00200     for(int pt = 0; pt < numCubPoints; pt++){
00201       std::cout << "\t Parameter point " 
00202       << std::setw(12) << std::right << edgeParamCubPts(pt, 0) << std::setw(10) << "  -->  " << "(" 
00203       << std::setw(10) << std::right << quadEdgePoints(pt, 0) << " , " 
00204       << std::setw(10) << std::right << quadEdgePoints(pt, 1) << ")\n";
00205     }    
00206     std::cout << "\n";
00207   }
00208   
00209   
00210   /* 
00211     Specification of integration points on 2-subcells (faces) of reference cells. Reference cells
00212     can have triangular, quadrilateral or a mixture of triangular and quadrilateral faces. Thus,
00213     parametrization domain of a face depends on that face's topology and is either the standard 
00214     2-simplex {(0,0), (1,0), (0,1)} for triangular faces or the standard 2-cube [-1,1]^2 for 
00215     quadrilateral faces. 
00216       
00217    1. Common tasks: definition of integration points in the standard 2-simplex and the standard
00218       2-cube. These steps are independent of parent cell topology:
00219    
00220       a. Instantiate a CubatureFactory object to create cubatures (already done!)
00221       b. Define parametrization domain for traingular faces as having Triangle<3> topology and for
00222          quadrilateral faces - as having Quadrilateral<4> topology. This is required by the 
00223          CubatureFactory in order to select cubature points and weights from the appropriate
00224          face parametrization domain. 
00225       c. Use CubatureFactory to select cubature of the desired degree for Triangle<3> and 
00226          Quadrilateral<4> topologies
00227       d. Allocate containers for the cubature points and weights on the parametrization domains.
00228    
00229    2. Parent cell topology specific tasks
00230    
00231       a. Select the parent cell topology
00232       b. Allocate containers for the images of the integration points from the parametrization
00233          domains on the reference faces
00234       c. Apply the face parametrization map to the points in the parametrization domain
00235    */
00236   
00237   // Step 1.b (Define the topology of the parametrization domain)
00238   CellTopology triFaceParam(shards::getCellTopologyData<shards::Triangle<3> >() );
00239   CellTopology quadFaceParam(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00240   
00241   
00242   // Step 1.c (selects Gauss rule of order 3 on [-1,1]^2 and a rule of order 3 on Triangle) 
00243   Teuchos::RCP<Cubature<double> > triFaceParamCubature = cubFactory.create(triFaceParam, 3); 
00244   Teuchos::RCP<Cubature<double> > quadFaceParamCubature = cubFactory.create(quadFaceParam, 3); 
00245   
00246   
00247   // Step 1.d - Triangle faces (allocate storage for cubature points)
00248   int triFaceCubDim    = triFaceParamCubature -> getDimension();
00249   int triFaceNumCubPts = triFaceParamCubature -> getNumPoints();
00250   
00251   FieldContainer<double> triFaceParamCubPts(triFaceNumCubPts, triFaceCubDim);
00252   FieldContainer<double> triFaceParamCubWts(triFaceNumCubPts);
00253   triFaceParamCubature -> getCubature(triFaceParamCubPts, triFaceParamCubWts);
00254   
00255   
00256   // Step 1.d - Quadrilateral faces (allocate storage for cubature points)
00257   int quadFaceCubDim    = quadFaceParamCubature -> getDimension();
00258   int quadFaceNumCubPts = quadFaceParamCubature -> getNumPoints();
00259   
00260   FieldContainer<double> quadFaceParamCubPts(quadFaceNumCubPts, quadFaceCubDim);
00261   FieldContainer<double> quadFaceParamCubWts(quadFaceNumCubPts);
00262   quadFaceParamCubature -> getCubature(quadFaceParamCubPts, quadFaceParamCubWts);
00263   
00264   
00265   
00266   std::cout \
00267     << "===============================================================================\n"\
00268     << "| EXAMPLE 2.1                                                                 |\n"
00269     << "| Face parametrizations for standard 3D cells: Tetrahedron                    |\n"\
00270     << "===============================================================================\n";
00271   
00272   // Step 2.a (select reference cell topology)
00273   CellTopology tet_4(getCellTopologyData<Tetrahedron<4> >() );
00274   
00275   
00276   // Step 2.b (allocate storage for points on a face of the reference cell)
00277   FieldContainer<double> tetFacePoints(triFaceNumCubPts, tet_4.getDimension() );
00278   
00279   
00280   // Step 2.c (same points are mapped to all faces, can also map different set to each face) 
00281   for(int faceOrd = 0; faceOrd < (int)tet_4.getSideCount(); faceOrd++){
00282     
00283     CellTools::mapToReferenceSubcell(tetFacePoints, 
00284                                      triFaceParamCubPts,
00285                                      2,
00286                                      faceOrd,
00287                                      tet_4);
00288     
00289     // Optional: print the vertices of the reference face 
00290     CellTools::printSubcellVertices(2, faceOrd, tet_4);
00291     
00292     for(int pt = 0; pt < triFaceNumCubPts; pt++){
00293       std::cout << "\t Parameter point (" 
00294       << std::setw(10) << std::right <<  triFaceParamCubPts(pt, 0) << " , "
00295       << std::setw(10) << std::right <<  triFaceParamCubPts(pt, 1) << ")  " 
00296       << std::setw(10) << " -->  " << "(" 
00297       << std::setw(10) << std::right << tetFacePoints(pt, 0) << " , " 
00298       << std::setw(10) << std::right << tetFacePoints(pt, 1) << " , " 
00299       << std::setw(10) << std::right << tetFacePoints(pt, 2) << ")\n";
00300     }    
00301     std::cout << "\n";
00302   }
00303   
00304   
00305   
00306   std::cout \
00307     << "===============================================================================\n"\
00308     << "| EXAMPLE 2.2                                                                 |\n"
00309     << "| Face parametrizations for standard 3D cells: Wedge                          |\n"\
00310     << "| Example of a reference cell that has two different kinds of faces           |\n"\
00311     << "===============================================================================\n";
00312   
00313   
00314   // Step 2.a (select reference cell topology)
00315   CellTopology wedge_6(getCellTopologyData<Wedge<6> >() );
00316   
00317   
00318   // Step 2.b (allocate storage for points on a face of the reference cell)
00319   //          Wedge<6> has Triangle<3> and Quadrilateral<4> faces. Two different arrays are needed
00320   //          to store the points on each face because different types integration rules are used
00321   //          on their respective parametrization domains and numbers of points defined by these
00322   //          rules do not necessarily match.
00323   FieldContainer<double> wedgeTriFacePoints(triFaceNumCubPts, wedge_6.getDimension() );
00324   FieldContainer<double> wedgeQuadFacePoints(quadFaceNumCubPts, wedge_6.getDimension() );
00325 
00326   
00327   // Step 2.c (for Wedge<6> need to distinguish Triangle<3> and Quadrilateral<4> faces) 
00328   for(int faceOrd = 0; faceOrd < (int)wedge_6.getSideCount(); faceOrd++){
00329     
00330     // Optional: print the vertices of the reference face 
00331     CellTools::printSubcellVertices(2, faceOrd, wedge_6);
00332     
00333     if( wedge_6.getKey(2, faceOrd) == shards::Triangle<3>::key ){
00334       CellTools::mapToReferenceSubcell(wedgeTriFacePoints, 
00335                                        triFaceParamCubPts,
00336                                        2,
00337                                        faceOrd,
00338                                        wedge_6);
00339       
00340       for(int pt = 0; pt < triFaceNumCubPts; pt++){
00341         std::cout << "\t Parameter point (" 
00342         << std::setw(10) << std::right <<  triFaceParamCubPts(pt, 0) << " , "
00343         << std::setw(10) << std::right <<  triFaceParamCubPts(pt, 1) << ")  "
00344         << std::setw(10) << "   -->    " << "(" 
00345         << std::setw(10) << std::right << wedgeTriFacePoints(pt, 0) << " , " 
00346         << std::setw(10) << std::right << wedgeTriFacePoints(pt, 1) << " , " 
00347         << std::setw(10) << std::right << wedgeTriFacePoints(pt, 2) << ")\n";
00348       }    
00349       std::cout << "\n";
00350     }
00351     else if(wedge_6.getKey(2, faceOrd) == shards::Quadrilateral<4>::key) {
00352       CellTools::mapToReferenceSubcell(wedgeQuadFacePoints, 
00353                                        quadFaceParamCubPts,
00354                                        2,
00355                                        faceOrd,
00356                                        wedge_6);
00357       
00358       for(int pt = 0; pt < quadFaceNumCubPts; pt++){
00359         std::cout << "\t Parameter point (" 
00360         << std::setw(10) << std::right <<  quadFaceParamCubPts(pt, 0) << " , "
00361         << std::setw(10) << std::right <<  quadFaceParamCubPts(pt, 1) << ")  "
00362         << std::setw(10) << "   -->    " << "(" 
00363         << std::setw(10) << std::right << wedgeQuadFacePoints(pt, 0) << " , " 
00364         << std::setw(10) << std::right << wedgeQuadFacePoints(pt, 1) << " , " 
00365         << std::setw(10) << std::right << wedgeQuadFacePoints(pt, 2) << ")\n";
00366       }    
00367       std::cout << "\n";
00368     }
00369     else {
00370       std::cout << " Invalid face encountered \n"; 
00371     }
00372   }
00373   
00374   return 0;
00375 }