Intrepid
http://trilinos.sandia.gov/packages/docs/r11.2/packages/intrepid/example/CellTools/example_04.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 "Shards_CellTopology.hpp"
00053 #include "Teuchos_GlobalMPISession.hpp"
00054 
00059 void vField(double& v1, double& v2, double& v3, 
00060             const double& x, const double& y, const double& z);
00061 
00062 using namespace std;
00063 using namespace Intrepid;
00064 
00065 int main(int argc, char *argv[]) {
00066 
00067   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00068 
00069   typedef CellTools<double>       CellTools;
00070   typedef RealSpaceTools<double>  RealSpaceTools;
00071   typedef shards::CellTopology    CellTopology;
00072 
00073  std::cout \
00074   << "===============================================================================\n" \
00075   << "|                                                                             |\n" \
00076   << "|                   Example use of the CellTools class                        |\n" \
00077   << "|                                                                             |\n" \
00078   << "|  1) Computation of face flux, for a given vector field, on a face workset   |\n" \
00079   << "|  2) Computation of edge circulation, for a given vector field, on a face    |\n" \
00080   << "|     workset.                                                                |\n" \
00081   << "|                                                                             |\n" \
00082   << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov)                      |\n" \
00083   << "|                      Denis Ridzal (dridzal@sandia.gov), or                  |\n" \
00084   << "|                      Kara Peterson (kjpeter@sandia.gov)                     |\n" \
00085   << "|                                                                             |\n" \
00086   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00087   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00088   << "|                                                                             |\n" \
00089   << "===============================================================================\n"\
00090   << "|  EXAMPLE 1: Computation of face flux on a face workset                      |\n"\
00091   << "===============================================================================\n";
00092 
00093   
00108   /*************************************************************************************************
00109     *
00110     *  Step 0. Face workset comprising of 1 face of a Hexahedron<8> cell
00111     *
00112     ************************************************************************************************/
00113   
00114   //   Step 0.a: Specify cell topology of the parent cell
00115   CellTopology hexahedron_8( shards::getCellTopologyData<shards::Hexahedron<8> >() );
00116   
00117   //   Step 0.b: Specify the vertices of the parent Hexahedron<8> cell
00118   int worksetSize    = 2;
00119   int pCellNodeCount = hexahedron_8.getVertexCount();
00120   int pCellDim       = hexahedron_8.getDimension();
00121   
00122   FieldContainer<double> hexNodes(worksetSize, pCellNodeCount, pCellDim);
00123   // cell 0 bottom face vertices:
00124   hexNodes(0, 0, 0) = 0.00;   hexNodes(0, 0, 1) = 0.00,   hexNodes(0, 0, 2) = 0.00;          
00125   hexNodes(0, 1, 0) = 1.00;   hexNodes(0, 1, 1) = 0.00,   hexNodes(0, 1, 2) = 0.00;
00126   hexNodes(0, 2, 0) = 1.00;   hexNodes(0, 2, 1) = 1.00,   hexNodes(0, 2, 2) = 0.00;
00127   hexNodes(0, 3, 0) = 0.00;   hexNodes(0, 3, 1) = 1.00,   hexNodes(0, 3, 2) = 0.00;
00128   // cell 0 top face vertices
00129   hexNodes(0, 4, 0) = 0.00;   hexNodes(0, 4, 1) = 0.00,   hexNodes(0, 4, 2) = 1.00;          
00130   hexNodes(0, 5, 0) = 1.00;   hexNodes(0, 5, 1) = 0.00,   hexNodes(0, 5, 2) = 1.00;
00131   hexNodes(0, 6, 0) = 1.00;   hexNodes(0, 6, 1) = 1.00,   hexNodes(0, 6, 2) = 1.00;
00132   hexNodes(0, 7, 0) = 0.00;   hexNodes(0, 7, 1) = 1.00,   hexNodes(0, 7, 2) = 1.00;
00133   
00134   // cell 1 bottom face vertices:
00135   hexNodes(1, 0, 0) = 0.00;   hexNodes(1, 0, 1) = 0.00,   hexNodes(1, 0, 2) = 0.00;          
00136   hexNodes(1, 1, 0) = 1.00;   hexNodes(1, 1, 1) = 0.00,   hexNodes(1, 1, 2) = 0.00;
00137   hexNodes(1, 2, 0) = 1.00;   hexNodes(1, 2, 1) = 1.00,   hexNodes(1, 2, 2) = 0.00;
00138   hexNodes(1, 3, 0) = 0.00;   hexNodes(1, 3, 1) = 1.00,   hexNodes(1, 3, 2) = 0.00;
00139   // cell 1 top face vertices
00140   hexNodes(1, 4, 0) = 0.00;   hexNodes(1, 4, 1) = 0.00,   hexNodes(1, 4, 2) = 1.00;          
00141   hexNodes(1, 5, 0) = 1.00;   hexNodes(1, 5, 1) = 0.00,   hexNodes(1, 5, 2) = 1.00;
00142   hexNodes(1, 6, 0) = 1.00;   hexNodes(1, 6, 1) = 1.00,   hexNodes(1, 6, 2) = 0.75;
00143   hexNodes(1, 7, 0) = 0.00;   hexNodes(1, 7, 1) = 1.00,   hexNodes(1, 7, 2) = 1.00;
00144   
00145   
00146   //   Step 0.c: Specify the face ordinal, relative to the reference cell, of the face workset
00147   int subcellDim = 2;
00148   int subcellOrd = 1;  
00149   
00150   
00151   
00152   /*************************************************************************************************
00153     *
00154     *  Step 1:    Obtain Gauss points on workset faces for Hexahedron<8> topology
00155     *       1.1   Define cubature factory, face parametrization domain and arrays for cubature points
00156     *       1.2   Select Gauss rule on D = [-1,1]x[-1,1] 
00157     *       1.3   Map Gauss points from D to reference face and workset faces
00158     *
00159     ************************************************************************************************/
00160   
00161   //   Step 1.1.a: Define CubatureFactory
00162   DefaultCubatureFactory<double>  cubFactory;   
00163   
00164   //   Step 1.1.b: Define topology of the face parametrization domain as [-1,1]x[-1,1]
00165   CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00166   
00167   //   Step 1.1.c: Define storage for cubature points and weights on [-1,1]x[-1,1]
00168   FieldContainer<double> paramGaussWeights;
00169   FieldContainer<double> paramGaussPoints;
00170   
00171   //   Step 1.1.d: Define storage for cubature points on a reference face
00172   FieldContainer<double> refGaussPoints;
00173   
00174   //   Step 1.1.f: Define storage for cubature points on workset faces
00175   FieldContainer<double> worksetGaussPoints;
00176 
00177   //----------------
00178   
00179   //   Step 1.2.a: selects Gauss rule of order 3 on [-1,1]x[-1,1]
00180   Teuchos::RCP<Cubature<double> > hexFaceCubature = cubFactory.create(paramQuadFace, 3); 
00181   
00182   //   Step 1.2.b allocate storage for cubature points on [-1,1]x[-1,1]
00183   int cubDim       = hexFaceCubature -> getDimension();
00184   int numCubPoints = hexFaceCubature -> getNumPoints();
00185   
00186   // Arrays must be properly sized for the specified set of Gauss points
00187   paramGaussPoints.resize(numCubPoints, cubDim);
00188   paramGaussWeights.resize(numCubPoints);
00189   hexFaceCubature -> getCubature(paramGaussPoints, paramGaussWeights);
00190   
00191   //----------------
00192   
00193   //   Step 1.3.a: Allocate storage for Gauss points on the reference face
00194   refGaussPoints.resize(numCubPoints, pCellDim);
00195 
00196   //   Step 1.3.b: Allocate storage for Gauss points on the face in the workset
00197   worksetGaussPoints.resize(worksetSize, numCubPoints, pCellDim);
00198 
00199   //   Step 1.3.c: Map Gauss points to reference face: paramGaussPoints -> refGaussPoints
00200   CellTools::mapToReferenceSubcell(refGaussPoints,
00201                                    paramGaussPoints,
00202                                    subcellDim,                      
00203                                    subcellOrd,
00204                                    hexahedron_8);
00205 
00206   //   Step 1.3.d: Map Gauss points from ref. face to face workset: refGaussPoints -> worksetGaussPoints
00207   CellTools::mapToPhysicalFrame(worksetGaussPoints,
00208                                 refGaussPoints,
00209                                 hexNodes,
00210                                 hexahedron_8);
00211   
00212   
00213   
00214   /*************************************************************************************************
00215     *
00216     *  Step 2.   Obtain (non-normalized) face normals at cubature points on workset faces
00217     *       2.1  Compute parent cell Jacobians at Gauss points on workset faces
00218     *       2.2  Compute face tangents on workset faces and their vector product
00219     *
00220     ************************************************************************************************/
00221   
00222   //   Step 2.1.a: Define and allocate storage for workset Jacobians
00223   FieldContainer<double> worksetJacobians(worksetSize, numCubPoints, pCellDim, pCellDim);
00224   
00225   //   Step 2.1.b: Compute Jacobians at Gauss pts. on reference face for all parent cells:
00226   CellTools::setJacobian(worksetJacobians,
00227                          refGaussPoints,
00228                          hexNodes,
00229                          hexahedron_8);
00230   
00231   //----------------
00232   
00233   //   Step 2.2.a: Allocate storage for face tangents and face normals
00234   FieldContainer<double> worksetFaceTu(worksetSize, numCubPoints, pCellDim);
00235   FieldContainer<double> worksetFaceTv(worksetSize, numCubPoints, pCellDim);
00236   FieldContainer<double> worksetFaceN(worksetSize, numCubPoints, pCellDim);
00237   
00238   //   Step 2.2.b: Compute face tangents
00239   CellTools::getPhysicalFaceTangents(worksetFaceTu,
00240                                      worksetFaceTv,
00241                                      worksetJacobians,
00242                                      subcellOrd,
00243                                      hexahedron_8);
00244   
00245   //   Step 2.2.c: Face outer normals (relative to parent cell) are uTan x vTan:
00246   RealSpaceTools::vecprod(worksetFaceN, worksetFaceTu, worksetFaceTv);
00247   
00248   
00249   
00250   /*************************************************************************************************
00251     *
00252     *  Step 3.   Evaluate the vector field u(x,y,z) at cubature points on workset faces
00253     *
00254     ************************************************************************************************/
00255   
00256   //   Step 3.a:  Allocate storage for vector field values at Gauss points on workset faces
00257   FieldContainer<double> worksetVFieldVals(worksetSize, numCubPoints, pCellDim);
00258   
00259   //   Step 3.b:  Compute vector field at Gauss points: here we take u(x,y,z) = (x,y,z)
00260   for(int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){
00261     for(int ptOrd = 0; ptOrd < numCubPoints; ptOrd++){
00262       
00263       double x = worksetGaussPoints(pCellOrd, ptOrd, 0);
00264       double y = worksetGaussPoints(pCellOrd, ptOrd, 1);
00265       double z = worksetGaussPoints(pCellOrd, ptOrd, 2);
00266 
00267       vField(worksetVFieldVals(pCellOrd, ptOrd, 0), 
00268              worksetVFieldVals(pCellOrd, ptOrd, 1), 
00269              worksetVFieldVals(pCellOrd, ptOrd, 2),  x, y, z);
00270       
00271     }// pt
00272   }//cell
00273   
00274 
00275   /*************************************************************************************************
00276     *
00277     *  Step 4.   Compute dot product of u(x,y,z) and the face normals, times the cubature weights
00278     *
00279     ************************************************************************************************/
00280   
00281   // Allocate storage for dot product of vector field and face normals at Gauss points
00282   FieldContainer<double> worksetFieldDotNormal(worksetSize, numCubPoints);
00283   
00284   // Compute the dot product
00285   RealSpaceTools::dot(worksetFieldDotNormal, worksetVFieldVals, worksetFaceN);
00286   
00287   // Allocate storage for face fluxes on the workset
00288   FieldContainer<double> worksetFluxes(worksetSize);
00289   
00290   //----------------
00291   
00292   // Integration loop (temporary)
00293   for(int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){
00294     worksetFluxes(pCellOrd) = 0.0;
00295     for(int pt = 0; pt < numCubPoints; pt++ ){
00296       worksetFluxes(pCellOrd) += worksetFieldDotNormal(pCellOrd, pt)*paramGaussWeights(pt);
00297     }// pt
00298   }//cell
00299   
00300   std::cout << "Face fluxes on workset faces : \n\n";
00301   for(int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){
00302     
00303     CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCellOrd, subcellDim, subcellOrd);
00304     std::cout << " Flux = " << worksetFluxes(pCellOrd) << "\n\n";
00305     
00306   }
00307   
00308   
00309   
00310   /*************************************************************************************************
00311     *
00312     *  Optional: print Gauss points and face normals at Gauss points
00313     *
00314     ************************************************************************************************/
00315 
00316   // Print Gauss points on [-1,1]x[-1,1] and their images on workset faces
00317   std::cout \
00318     << "===============================================================================\n" \
00319     << "|                        Gauss points on workset faces:                       |\n" \
00320     << "===============================================================================\n";
00321   for(int pCell = 0; pCell < worksetSize; pCell++){
00322     
00323     CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd);
00324     
00325     for(int pt = 0; pt < numCubPoints; pt++){
00326       std::cout << "\t 2D Gauss point (" 
00327       << std::setw(8) << std::right <<  paramGaussPoints(pt, 0) << ", "
00328       << std::setw(8) << std::right <<  paramGaussPoints(pt, 1) << ")  " 
00329       << std::setw(8) << " -->  " << "(" 
00330       << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 0) << ", " 
00331       << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 1) << ", " 
00332       << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 2) << ")\n";
00333     }    
00334     std::cout << "\n\n";      
00335   }//pCell
00336   
00337   
00338   // Print face normals at Gauss points on workset faces
00339   std::cout \
00340     << "===============================================================================\n" \
00341     << "|          Face normals (non-unit) at Gauss points on workset faces:          |\n" \
00342     << "===============================================================================\n";
00343   for(int pCell = 0; pCell < worksetSize; pCell++){
00344     
00345     CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd);
00346     
00347     for(int pt = 0; pt < numCubPoints; pt++){
00348       std::cout << "\t 3D Gauss point: (" 
00349       << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 0) << ", " 
00350       << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 1) << ", " 
00351       << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 2) << ")" 
00352       << std::setw(8) << " out. normal:  " << "(" 
00353       << std::setw(8) << std::right << worksetFaceN(pCell, pt, 0) << ", " 
00354       << std::setw(8) << std::right << worksetFaceN(pCell, pt, 1) << ", " 
00355       << std::setw(8) << std::right << worksetFaceN(pCell, pt, 2) << ")\n";
00356     }    
00357     std::cout << "\n";      
00358   }//pCell
00359   
00360   return 0;
00361 }
00362 
00363 /*************************************************************************************************
00364  *
00365  *  Definition of the vector field function
00366  *
00367  ************************************************************************************************/
00368 
00369 
00370 void vField(double& v1, double& v2, double& v3, const double& x, const double& y, const double& z)
00371 {
00372   v1 = x;
00373   v2 = y;
00374   v3 = z;
00375 }
00376 
00377