Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/example/CellTools/example_03.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 
00055 using namespace std;
00056 using namespace Intrepid;
00057 
00058 int main(int argc, char *argv[]) {
00059 
00060   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00061 
00062   typedef CellTools<double>       CellTools;
00063   typedef shards::CellTopology    CellTopology;
00064 
00065  std::cout \
00066   << "===============================================================================\n" \
00067   << "|                                                                             |\n" \
00068   << "|                   Example use of the CellTools class                        |\n" \
00069   << "|                                                                             |\n" \
00070   << "|  1) Definition of integration points on edge and face worksets              |\n" \
00071   << "|  2) Computation of face normals and edge tangents on face and edge worksets |\n" \
00072   << "|  2) Computation of side normals on face and edge worksets                   |\n" \
00073   << "|                                                                             |\n" \
00074   << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov)                      |\n" \
00075   << "|                      Denis Ridzal (dridzal@sandia.gov), or                  |\n" \
00076   << "|                      Kara Peterson (kjpeter@sandia.gov)                     |\n" \
00077   << "|                                                                             |\n" \
00078   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00079   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00080   << "|                                                                             |\n" \
00081   << "===============================================================================\n"\
00082   << "|  EXAMPLE 1: Definition of integration points on a Triangle edge workset     |\n"\
00083   << "===============================================================================\n";
00084  
00085  /*
00086   *  1. Common tasks for getting integration points on edge worksets
00087   */
00088   
00089   // Step 1.a: Define CubatureFactory
00090   DefaultCubatureFactory<double>  cubFactory;   
00091   
00092   // Step 1.b: Define topology of the edge parametrization domain [-1,1]
00093   CellTopology paramEdge(shards::getCellTopologyData<shards::Line<2> >() );
00094   
00095   // Step 1.c: Define storage for cubature points and weights on [-1,1]
00096   FieldContainer<double> paramEdgeWeights;
00097   FieldContainer<double> paramEdgePoints;
00098   
00099   // Step 1.d: Define storage for cubature points on a reference edge
00100   FieldContainer<double> refEdgePoints;
00101   
00102   // Step 1.f: Define storage for cubature points on workset edges
00103   FieldContainer<double> worksetEdgePoints;
00104 
00105 
00106   
00107   /*
00108    *  2. Define edge workset comprising of 4 edges corresponding to reference edge 2 on Triangle<3>
00109    */
00110   
00111   // Step 2.a: Specify cell topology of the parent cell
00112   CellTopology triangle_3( shards::getCellTopologyData<shards::Triangle<3> >() );
00113   
00114   // Step 2.b: Specify the vertices of 4 parent Triangle<3> cells
00115   int worksetSize    = 4;
00116   int pCellNodeCount = triangle_3.getVertexCount();
00117   int pCellDim       = triangle_3.getDimension();
00118   
00119   FieldContainer<double> triNodes(worksetSize, pCellNodeCount, pCellDim);
00120   // Triangle 0
00121   triNodes(0, 0, 0) = 0.5;   triNodes(0, 0, 1) = 2.0;
00122   triNodes(0, 1, 0) = 0.0;   triNodes(0, 1, 1) = 1.0;
00123   triNodes(0, 2, 0) = 1.0;   triNodes(0, 2, 1) = 1.0;
00124   // Triangle 1
00125   triNodes(1, 0, 0) = 1.0;   triNodes(1, 0, 1) = 1.0;               
00126   triNodes(1, 1, 0) = 1.0;   triNodes(1, 1, 1) = 0.0;
00127   triNodes(1, 2, 0) = 0.0;   triNodes(1, 2, 1) = 0.0;
00128   // Triangle 2
00129   triNodes(2, 0, 0) = 0.0;   triNodes(2, 0, 1) = 0.0;
00130   triNodes(2, 1, 0) = 1.0;   triNodes(2, 1, 1) = 1.0;
00131   triNodes(2, 2, 0) = 0.0;   triNodes(2, 2, 1) = 1.0;
00132   // Triangle 3
00133   triNodes(3, 0, 0) = 1.0;   triNodes(3, 0, 1) = 1.0;
00134   triNodes(3, 1, 0) = 2.5;   triNodes(3, 1, 1) = 1.5;
00135   triNodes(3, 2, 0) = 0.5;   triNodes(3, 2, 1) = 2.0;
00136   
00137   // Step 2.c: Specify dimension and ordinal (relative to ref. cell) of the subcells in the workset
00138   int subcellDim = 1;
00139   int subcellOrd = 2;  
00140   
00141   
00142   
00143   /*
00144    *  3. Obtain the desired Gauss rule on the edge parametrization domain [-1,1]
00145    */
00146   
00147   // Step 3.a: selects Gauss rule of order 6 on [-1,1]
00148   Teuchos::RCP<Cubature<double> > triEdgeCubature = cubFactory.create(paramEdge, 6); 
00149   
00150   // Step 3.b allocate storage for cubature points on [-1,1]
00151   int cubDim       = triEdgeCubature -> getDimension();
00152   int numCubPoints = triEdgeCubature -> getNumPoints();
00153   
00154   // Arrays must be properly sized for the specified set of Gauss points
00155   paramEdgePoints.resize(numCubPoints, cubDim);
00156   paramEdgeWeights.resize(numCubPoints);
00157   triEdgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights);
00158 
00159   
00160   
00161   /*
00162    *  4. Map Gauss points from [-1,1] to every edge in the workset
00163    */
00164 
00165   // Step 4.a: Allocate storage for integration points on the reference edge
00166   refEdgePoints.resize(numCubPoints, pCellDim);
00167  
00168   // Step 4.b: Allocate storage for integration points on the edge workset
00169   worksetEdgePoints.resize(worksetSize, numCubPoints, pCellDim);
00170   
00171   // Step 4.c: Map Gauss points to reference edge: paramEdgePoints -> refEdgePoints
00172   CellTools::mapToReferenceSubcell(refEdgePoints,
00173                                    paramEdgePoints,
00174                                    subcellDim,                      
00175                                    subcellOrd,
00176                                    triangle_3);
00177   
00178   // Step 4.d: Map Gauss points from ref. edge to edge workset: refEdgePoints -> worksetEdgePoints
00179   CellTools::mapToPhysicalFrame(worksetEdgePoints,
00180                                 refEdgePoints,
00181                                 triNodes,
00182                                 triangle_3);
00183   
00184   
00185   
00186   /*
00187    *  5. Print Gauss points on the edges of the workset
00188    */
00189   for(int pCell = 0; pCell < worksetSize; pCell++){
00190       
00191     CellTools::printWorksetSubcell(triNodes, triangle_3, pCell, subcellDim, subcellOrd);
00192      
00193       for(int pt = 0; pt < numCubPoints; pt++){
00194         std::cout << "\t 1D Gauss point " 
00195         << std::setw(12) << std::right << paramEdgePoints(pt, 0) << std::setw(10) << "  -->  " << "(" 
00196         << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 0) << ", " 
00197         << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 1) << ")\n";
00198       }    
00199       std::cout << "\n";      
00200   }//pCell
00201   
00202   
00203   
00204   std::cout << "\n" \
00205     << "===============================================================================\n"\
00206     << "|  EXAMPLE 2: Definition of integration points on a Quadrilateral edge workset|\n"\
00207     << "===============================================================================\n";
00208   
00209   /*
00210    *  2. Define edge workset comprising of 2 edges corresponding to reference edge 2 on Quadrilateral<4>
00211    *     (Can reuse Step 1, Example 1 because edge parametrization domain is always [-1,1])
00212    */
00213   
00214   // Step 2.a: Specify cell topology of the parent cell
00215   CellTopology quadrilateral_4( shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00216   
00217   // Step 2.b: Specify the vertices of 2 parent Quadrilateral<4> cells
00218   worksetSize    = 2;
00219   pCellNodeCount = quadrilateral_4.getVertexCount();
00220   pCellDim       = quadrilateral_4.getDimension();
00221   
00222   FieldContainer<double> quadNodes(worksetSize, pCellNodeCount, pCellDim);
00223   // Quadrilateral 0
00224   quadNodes(0, 0, 0) = 1.00;   quadNodes(0, 0, 1) = 1.00;               
00225   quadNodes(0, 1, 0) = 2.00;   quadNodes(0, 1, 1) = 0.75;
00226   quadNodes(0, 2, 0) = 1.75;   quadNodes(0, 2, 1) = 2.00;  
00227   quadNodes(0, 3, 0) = 1.25;   quadNodes(0, 3, 1) = 2.00; 
00228   // Quadrilateral 1
00229   quadNodes(1, 0, 0) = 2.00;   quadNodes(1, 0, 1) = 0.75;               
00230   quadNodes(1, 1, 0) = 3.00;   quadNodes(1, 1, 1) = 1.25;
00231   quadNodes(1, 2, 0) = 2.75;   quadNodes(1, 2, 1) = 2.25;
00232   quadNodes(1, 3, 0) = 1.75;   quadNodes(1, 3, 1) = 2.00;
00233   
00234   // Step 2.c: Specify dimension and ordinal (relative to ref. cell) of the subcells in the workset
00235   subcellDim = 1;
00236   subcellOrd = 2;  
00237   
00238   /*
00239    *  3. Obtain the desired Gauss rule on the edge parametrization domain [-1,1]
00240    */
00241   
00242   // Step 3.a: selects Gauss rule of order 4 on [-1,1]
00243   Teuchos::RCP<Cubature<double> > quadEdgeCubature = cubFactory.create(paramEdge, 6); 
00244   
00245   // Step 3.b: allocate storage for cubature points 
00246   cubDim       = quadEdgeCubature -> getDimension();
00247   numCubPoints = quadEdgeCubature -> getNumPoints();
00248   
00249   // Arrays must be properly sized for the specified set of Gauss points
00250   paramEdgePoints.resize(numCubPoints, cubDim);
00251   paramEdgeWeights.resize(numCubPoints);
00252   quadEdgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights);
00253   
00254   /*
00255    *  4. Map Gauss points from [-1,1] to every edge in the workset
00256    */
00257   
00258   // Step 4.a: Allocate storage for integration points on the reference edge
00259   refEdgePoints.resize(numCubPoints, pCellDim);
00260   
00261   // Step 4.b: Allocate storage for integration points on the edge workset
00262   worksetEdgePoints.resize(worksetSize, numCubPoints, pCellDim);
00263   
00264   // Step 4.c: Map Gauss points to reference edge: paramEdgePoints -> refEdgePoints
00265   CellTools::mapToReferenceSubcell(refEdgePoints,
00266                                    paramEdgePoints,
00267                                    subcellDim,                      
00268                                    subcellOrd,
00269                                    quadrilateral_4);
00270   
00271   // Step 4.d: Map Gauss points from ref. edge to edge workset: refEdgePoints -> worksetEdgePoints
00272   CellTools::mapToPhysicalFrame(worksetEdgePoints,
00273                                 refEdgePoints,
00274                                 quadNodes,
00275                                 quadrilateral_4);
00276   
00277   /*
00278    *  5. Print Gauss points on the edges of the workset
00279    */
00280   for(int pCell = 0; pCell < worksetSize; pCell++){
00281     
00282     CellTools::printWorksetSubcell(quadNodes, quadrilateral_4, pCell, subcellDim, subcellOrd);
00283     
00284     for(int pt = 0; pt < numCubPoints; pt++){
00285       std::cout << "\t 1D Gauss point " 
00286       << std::setw(12) << std::right << paramEdgePoints(pt, 0) << std::setw(10) << "  -->  " << "(" 
00287       << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 0) << ", " 
00288       << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 1) << ")\n";
00289     }    
00290     std::cout << "\n";      
00291   }//pCell
00292   
00293   
00294   
00295   std::cout << "\n" \
00296     << "===============================================================================\n"\
00297     << "| EXAMPLE 3: Edge tangents at Gauss points on a Quadrilateral edge workset    |\n"\
00298     << "===============================================================================\n";
00299   
00300   /*  This task requires Gauss points on edge parametrization domain [-1,1] and on the 
00301    *  reference edge whose ordinal matches the edge workset ordinal. This repeats the first few
00302    *  steps from Example 2:
00303    *  
00304    *  1. Define cubature factory and topology for edge parametrization domain [-1,1];
00305    *     (Can reuse Step 1, Example 1 because edge parametrization domain is always [-1,1])
00306    *  2. Define an edge workset;
00307    *  3. Obtain the desired Gauss rule on the edge parametrization domain [-1,1];
00308    *  4. Map Gauss points from [-1,1] to reference edge 
00309    *     NOTE: this example only demonstrates computation of the edge tangents and so, Gauss 
00310    *     points on the edge workset are not needed. Thus we skip mapping Gauss points from 
00311    *     reference edge to edge workset
00312    *
00313    *  5. Compute Jacobians at Gauss points on reference edge for all cells in the workset:
00314    */
00315   
00316   // Step 5.a: Define and allocate storage for workset Jacobians
00317   FieldContainer<double> worksetJacobians(worksetSize, numCubPoints, pCellDim, pCellDim);
00318   
00319   // Step 5.b: Compute Jacobians at Gauss pts. on reference edge for all parent cells:
00320   CellTools::setJacobian(worksetJacobians,
00321                          refEdgePoints,
00322                          quadNodes,
00323                          quadrilateral_4);
00324   /*
00325    * 6. Get the (non-normalized) edge tangents for the edge workset:
00326    */
00327   // Step 6.a: Allocate storage for edge tangents
00328   FieldContainer<double> edgeTangents(worksetSize, numCubPoints, pCellDim);
00329   
00330   // Step 6.b: Compute the edge tangents:
00331   CellTools::getPhysicalEdgeTangents(edgeTangents,
00332                                      worksetJacobians,
00333                                      subcellOrd,
00334                                      quadrilateral_4); 
00335   
00336   // Step 6.c: Print edge tangents at Gauss points on workset edges (these Gauss points were computed in Example 2)
00337   std::cout 
00338     << "Edge tangents computed by CellTools::getPhysicalEdgeTangents.\n"
00339     << "Local edge ordinal = " << subcellOrd <<"\n";
00340   
00341   for(int pCell = 0; pCell < worksetSize; pCell++){
00342     
00343     CellTools::printWorksetSubcell(quadNodes, quadrilateral_4, pCell, subcellDim, subcellOrd);
00344     
00345     for(int pt = 0; pt < numCubPoints; pt++){
00346       std::cout << "\t 2D Gauss point: (" 
00347       << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 0) << ", " 
00348       << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 1) << ")  " 
00349       << std::setw(8) << " edge tangent:  " << "(" 
00350       << std::setw(8) << std::right << edgeTangents(pCell, pt, 0) << ", " 
00351       << std::setw(8) << std::right << edgeTangents(pCell, pt, 1) << ")\n";
00352     }    
00353     std::cout << "\n";      
00354   }//pCell
00355 
00356   std::cout << "\n" \
00357     << "===============================================================================\n"\
00358     << "| EXAMPLE 4: Side normals at Gauss points on a Quadrilateral side workset     |\n"\
00359     << "===============================================================================\n";
00360 
00361   /* For this task we reuse the edge workset from Example 3 as a side workset and compute
00362    * normals to the sides in that set. This task repeats the first 5 steps from Example 3
00363    * The only difference is that in Step 6 we call CellTools::getPhysicalSideNormals
00364    */
00365   /*
00366    * 6. Get the (non-normalized) side normals for the side (edge) workset:
00367    */
00368   // Step 6.a: Allocate storage for side normals
00369   FieldContainer<double> sideNormals(worksetSize, numCubPoints, pCellDim);
00370   
00371   // Step 6.b: Compute the side normals:
00372   CellTools::getPhysicalSideNormals(sideNormals,
00373                                     worksetJacobians,
00374                                     subcellOrd,
00375                                     quadrilateral_4); 
00376   
00377   // Step 6.c: Print side normals at Gauss points on workset sides (these Gauss points were computed in Example 2)
00378   std::cout 
00379     << "Side normals computed by CellTools::getPhysicalSideNormals.\n"
00380     << "Local side (edge) ordinal = " << subcellOrd <<"\n";
00381   
00382   for(int pCell = 0; pCell < worksetSize; pCell++){
00383     
00384     CellTools::printWorksetSubcell(quadNodes, quadrilateral_4, pCell, subcellDim, subcellOrd);
00385     
00386     for(int pt = 0; pt < numCubPoints; pt++){
00387       std::cout << "\t 2D Gauss point: (" 
00388       << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 0) << ", " 
00389       << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 1) << ")  " 
00390       << std::setw(8) << " side normal:  " << "(" 
00391       << std::setw(8) << std::right << sideNormals(pCell, pt, 0) << ", " 
00392       << std::setw(8) << std::right << sideNormals(pCell, pt, 1) << ")\n";
00393     }    
00394     std::cout << "\n";      
00395   }//pCell
00396   
00397     
00398   std::cout << "\n" \
00399     << "===============================================================================\n"\
00400     << "| EXAMPLE 5: Definition of integration points on a Hexahedron edge workset    |\n"\
00401     << "===============================================================================\n";
00402   
00403   /*
00404    *  2. Define edge workset comprising of 1 edge corresponding to reference edge 10 on Hexahedron<8>
00405    *     (Can reuse Step 1, Example 1 because edge parametrization domain is always [-1,1])
00406    */
00407   
00408   // Step 2.a: Specify cell topology of the parent cell
00409   CellTopology hexahedron_8( shards::getCellTopologyData<shards::Hexahedron<8> >() );
00410   
00411   // Step 2.b: Specify the vertices of the parent Hexahedron<8> cell
00412   worksetSize    = 1;
00413   pCellNodeCount = hexahedron_8.getVertexCount();
00414   pCellDim       = hexahedron_8.getDimension();
00415 
00416   FieldContainer<double> hexNodes(worksetSize, pCellNodeCount, pCellDim);
00417   // bottom face vertices
00418   hexNodes(0, 0, 0) = 0.00;   hexNodes(0, 0, 1) = 0.00,   hexNodes(0, 0, 2) = 0.00;          
00419   hexNodes(0, 1, 0) = 1.00;   hexNodes(0, 1, 1) = 0.00,   hexNodes(0, 1, 2) = 0.00;
00420   hexNodes(0, 2, 0) = 1.00;   hexNodes(0, 2, 1) = 1.00,   hexNodes(0, 2, 2) = 0.00;
00421   hexNodes(0, 3, 0) = 0.00;   hexNodes(0, 3, 1) = 1.00,   hexNodes(0, 3, 2) = 0.00;
00422   // top face vertices
00423   hexNodes(0, 4, 0) = 0.00;   hexNodes(0, 4, 1) = 0.00,   hexNodes(0, 4, 2) = 1.00;          
00424   hexNodes(0, 5, 0) = 1.00;   hexNodes(0, 5, 1) = 0.00,   hexNodes(0, 5, 2) = 1.00;
00425   hexNodes(0, 6, 0) = 1.00;   hexNodes(0, 6, 1) = 1.00,   hexNodes(0, 6, 2) = 0.75;
00426   hexNodes(0, 7, 0) = 0.00;   hexNodes(0, 7, 1) = 1.00,   hexNodes(0, 7, 2) = 1.00;
00427   
00428   // An alternative hex obtained by intersection of the unit cube [0,1]^3 and the plane
00429   // z = 1 - 1/4x - 1/4y. The top face (local ordinal 5) of the resulting hex lies in this plane
00430   // and has the same normal vector parallel to (1/4,1/4,1). This workset allows to test face normals
00431   FieldContainer<double> hexNodesAlt(worksetSize, pCellNodeCount, pCellDim);
00432   // bottom face vertices
00433   hexNodesAlt(0, 0, 0) = 0.00;   hexNodesAlt(0, 0, 1) = 0.00,   hexNodesAlt(0, 0, 2) = 0.00;          
00434   hexNodesAlt(0, 1, 0) = 1.00;   hexNodesAlt(0, 1, 1) = 0.00,   hexNodesAlt(0, 1, 2) = 0.00;
00435   hexNodesAlt(0, 2, 0) = 1.00;   hexNodesAlt(0, 2, 1) = 1.00,   hexNodesAlt(0, 2, 2) = 0.00;
00436   hexNodesAlt(0, 3, 0) = 0.00;   hexNodesAlt(0, 3, 1) = 1.00,   hexNodesAlt(0, 3, 2) = 0.00;
00437   // top face vertices
00438   hexNodesAlt(0, 4, 0) = 0.00;   hexNodesAlt(0, 4, 1) = 0.00,   hexNodesAlt(0, 4, 2) = 1.00;          
00439   hexNodesAlt(0, 5, 0) = 1.00;   hexNodesAlt(0, 5, 1) = 0.00,   hexNodesAlt(0, 5, 2) = 0.75;
00440   hexNodesAlt(0, 6, 0) = 1.00;   hexNodesAlt(0, 6, 1) = 1.00,   hexNodesAlt(0, 6, 2) = 0.50;
00441   hexNodesAlt(0, 7, 0) = 0.00;   hexNodesAlt(0, 7, 1) = 1.00,   hexNodesAlt(0, 7, 2) = 0.75;  
00442   
00443   // Step 2.c: Specify the edge ordinal, relative to the reference cell, of the edge workset
00444   subcellDim = 1;
00445   subcellOrd = 5;  
00446   
00447   /*
00448    *  3. Obtain the desired Gauss rule on the edge parametrization domain [-1,1]
00449    */
00450   
00451   // Step 3.a: selects Gauss rule of order 4 on [-1,1]
00452   Teuchos::RCP<Cubature<double> > hexEdgeCubature = cubFactory.create(paramEdge, 4); 
00453   
00454   // Step 3.b allocate storage for cubature points
00455   cubDim       = hexEdgeCubature -> getDimension();
00456   numCubPoints = hexEdgeCubature -> getNumPoints();
00457   
00458   // Arrays must be properly sized for the specified set of Gauss points
00459   paramEdgePoints.resize(numCubPoints, cubDim);
00460   paramEdgeWeights.resize(numCubPoints);
00461   hexEdgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights);
00462   
00463   /*
00464    *  4. Map Gauss points from [-1,1] to every edge in the workset
00465    */
00466   
00467   // Step 4.a: Allocate storage for integration points on the reference edge
00468   refEdgePoints.resize(numCubPoints, pCellDim);
00469   
00470   // Step 4.b: Allocate storage for integration points on the edge workset
00471   worksetEdgePoints.resize(worksetSize, numCubPoints, pCellDim);
00472   
00473   // Step 4.c: Map Gauss points to reference edge: paramEdgePoints -> refEdgePoints
00474   CellTools::mapToReferenceSubcell(refEdgePoints,
00475                                    paramEdgePoints,
00476                                    subcellDim,                      
00477                                    subcellOrd,
00478                                    hexahedron_8);
00479   
00480   // Step 4.d: Map Gauss points from ref. edge to edge workset: refEdgePoints -> worksetEdgePoints
00481   CellTools::mapToPhysicalFrame(worksetEdgePoints,
00482                                 refEdgePoints,
00483                                 hexNodes,
00484                                 hexahedron_8);
00485   
00486   /*
00487    *  5. Print Gauss points on the edges of the workset
00488    */
00489   for(int pCell = 0; pCell < worksetSize; pCell++){
00490     
00491     CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd);
00492     
00493     for(int pt = 0; pt < numCubPoints; pt++){
00494       std::cout << "\t 1D Gauss point " 
00495       << std::setw(12) << std::right << paramEdgePoints(pt, 0) << std::setw(10) << "  -->  " << "(" 
00496       << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 0) << ", " 
00497       << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 1) << ", " 
00498       << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 2) << ")\n";
00499     }    
00500     std::cout << "\n";      
00501   }//pCell
00502   
00503   
00504   
00505   std::cout << "\n" \
00506     << "===============================================================================\n"\
00507     << "| EXAMPLE 6: Edge tangents at Gauss points on a Hexahedron edge workset       |\n"\
00508     << "===============================================================================\n";
00509   
00510   /*  This task requires Gauss points on edge parametrization domain [-1,1] and on the 
00511    *  reference edge whose ordinal matches the edge workset ordinal. This repeats the first few
00512    *  steps from Example 5:
00513    *  
00514    *  1. Define cubature factory and topology for edge parametrization domain [-1,1];
00515    *     (Can reuse Step 1, Example 1 because edge parametrization domain is always [-1,1])
00516    *  2. Define an edge workset;
00517    *  3. Obtain the desired Gauss rule on the edge parametrization domain [-1,1];
00518    *  4. Map Gauss points from [-1,1] to reference edge 
00519    *     NOTE: this example only demonstrates computation of the edge tangents and so, Gauss 
00520    *     points on the edge workset are not needed. Thus we skip mapping Gauss points from 
00521    *     reference edge to edge workset
00522    *
00523    *  5. Compute Jacobians at Gauss points on reference edge for all cells in the workset:
00524    */
00525   
00526   // Step 5.a: Define and allocate storage for workset Jacobians
00527   worksetJacobians.resize(worksetSize, numCubPoints, pCellDim, pCellDim);
00528   
00529   // Step 5.b: Compute Jacobians at Gauss pts. on reference edge for all parent cells:
00530   CellTools::setJacobian(worksetJacobians,
00531                          refEdgePoints,
00532                          hexNodes,
00533                          hexahedron_8);
00534   
00535   /*
00536    * 6. Get the (non-normalized) edge tangents for the edge workset:
00537    */
00538   // Step 6.a: Allocate storage for edge tangents
00539   edgeTangents.resize(worksetSize, numCubPoints, pCellDim);
00540   
00541   // Step 6.b: Compute the edge tangents:
00542   CellTools::getPhysicalEdgeTangents(edgeTangents,
00543                                      worksetJacobians,
00544                                      subcellOrd,
00545                                      hexahedron_8); 
00546   
00547   // Step 6.c: Print edge tangents at Gauss points on workset edges (these Gauss points were computed in Example 5)
00548   std::cout 
00549     << "Edge tangents computed by CellTools::getPhysicalEdgeTangents.\n"
00550     << "Local edge ordinal = " << subcellOrd <<"\n";
00551   
00552   for(int pCell = 0; pCell < worksetSize; pCell++){
00553     
00554     CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd);
00555     
00556     for(int pt = 0; pt < numCubPoints; pt++){
00557       std::cout << "\t 3D Gauss point: (" 
00558       << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 0) << ", " 
00559       << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 1) << ", " 
00560       << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 2) << ")  " 
00561       << std::setw(8) << " edge tangent:  " << "(" 
00562       << std::setw(8) << std::right << edgeTangents(pCell, pt, 0) << ", " 
00563       << std::setw(8) << std::right << edgeTangents(pCell, pt, 1) << ", " 
00564       << std::setw(8) << std::right << edgeTangents(pCell, pt, 2) << ")\n";
00565     }    
00566     std::cout << "\n";      
00567   }//pCell
00568   
00569    
00570  std::cout << "\n" \
00571     << "===============================================================================\n"\
00572     << "| EXAMPLE 7: Definition of integration points on a Hexahedron face workset    |\n"\
00573     << "===============================================================================\n";
00574   /*
00575    *  1. Common tasks for getting integration points on face worksets: 
00576    *  1.a: can reuse the cubature factory, but need parametrization domain for the faces
00577    */
00578   
00579   // Step 1.b: Define topology of the face parametrization domain as [-1,1]x[-1,1]
00580   CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00581   
00582   // Step 1.c: Define storage for cubature points and weights on [-1,1]x[-1,1]
00583   FieldContainer<double> paramFaceWeights;
00584   FieldContainer<double> paramFacePoints;
00585   
00586   // Step 1.d: Define storage for cubature points on a reference face
00587   FieldContainer<double> refFacePoints;
00588   
00589   // Step 1.f: Define storage for cubature points on workset faces
00590   FieldContainer<double> worksetFacePoints;
00591     
00592   /*
00593    *  2. Define face workset comprising of 1 face corresponding to reference face 5 on Hexahedron<8>
00594    *  2.a: Reuse the parent cell topology from Example 3
00595    *  2.b: Reuse the vertices from Example 3
00596    */
00597   
00598   // Step 2.c: Specify dimension and ordinal (relative to ref. cell) of the subcells in the workset
00599   subcellDim = 2;
00600   subcellOrd = 5;  
00601   
00602   /*
00603    *  3. Obtain the desired Gauss rule on the face parametrization domain [-1,1]x[-1,1]
00604    */
00605   
00606   // Step 3.a: selects Gauss rule of order 3 on [-1,1]x[-1,1]
00607   Teuchos::RCP<Cubature<double> > hexFaceCubature = cubFactory.create(paramQuadFace, 3); 
00608   
00609   // Step 3.b allocate storage for cubature points on [-1,1]x[-1,1]
00610   cubDim       = hexFaceCubature -> getDimension();
00611   numCubPoints = hexFaceCubature -> getNumPoints();
00612   
00613   // Arrays must be properly sized for the specified set of Gauss points
00614   paramFacePoints.resize(numCubPoints, cubDim);
00615   paramFaceWeights.resize(numCubPoints);
00616   hexFaceCubature -> getCubature(paramFacePoints, paramFaceWeights);
00617   
00618   /*
00619    *  4. Map Gauss points from [-1,1]x[-1,1] to every face in the workset
00620    */
00621     
00622   // Step 4.a: Allocate storage for integration points on the reference face
00623   refFacePoints.resize(numCubPoints, pCellDim);
00624 
00625   // Step 4.b: Allocate storage for integration points on the face in the workset
00626   worksetFacePoints.resize(worksetSize, numCubPoints, pCellDim);
00627 
00628   // Step 4.c: Map Gauss points to reference face: paramFacePoints -> refFacePoints
00629   CellTools::mapToReferenceSubcell(refFacePoints,
00630                                    paramFacePoints,
00631                                    subcellDim,                      
00632                                    subcellOrd,
00633                                    hexahedron_8);
00634 
00635   // Step 4.d: Map Gauss points from ref. face to face workset: refFacePoints -> worksetFacePoints
00636   CellTools::mapToPhysicalFrame(worksetFacePoints,
00637                                 refFacePoints,
00638                                 hexNodesAlt,
00639                                 hexahedron_8);
00640   
00641   
00642   /*
00643    *  5. Print Gauss points on the faces of the workset
00644    */
00645   for(int pCell = 0; pCell < worksetSize; pCell++){
00646     
00647     CellTools::printWorksetSubcell(hexNodesAlt, hexahedron_8, pCell, subcellDim, subcellOrd);
00648     
00649     for(int pt = 0; pt < numCubPoints; pt++){
00650       std::cout << "\t 2D Gauss point (" 
00651       << std::setw(8) << std::right <<  paramFacePoints(pt, 0) << ", "
00652       << std::setw(8) << std::right <<  paramFacePoints(pt, 1) << ")  " 
00653       << std::setw(8) << " -->  " << "(" 
00654       << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 0) << ", " 
00655       << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 1) << ", " 
00656       << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 2) << ")\n";
00657     }    
00658     std::cout << "\n";      
00659   }//pCell
00660   
00661   
00662   std::cout << "\n" \
00663     << "===============================================================================\n"\
00664     << "| EXAMPLE 8: Face normals at Gauss points on a Hexahedron face workset        |\n"\
00665     << "===============================================================================\n";
00666   
00667   /*  This task requires Gauss points on face parametrization domain [-1,1]x[-1,1] and on the 
00668    *  reference face whose ordinal matches the face workset ordinal. This repeats the first few
00669    *  steps from Example 5:
00670    *  
00671    *  1. Define cubature factory and topology for face parametrization domain [-1,1]x[-1,1];
00672    *  2. Define a face workset;
00673    *  3. Obtain the desired Gauss rule on the face parametrization domain [-1,1]x[-1,1];
00674    *  4. Map Gauss points from [-1,1]x[-1,1] to reference face
00675    *     NOTE: this example only demonstrates computation of the face normals and so, Gaus 
00676    *     points on the face workset are not needed. Thus we skip mapping Gauss points from 
00677    *     reference face to face workset
00678    *
00679    *  5. Compute Jacobians at Gauss points on reference face for all cells in the workset:
00680    */
00681   
00682   // Step 5.a: Define and allocate storage for workset Jacobians (reuse FC from example 5)
00683   worksetJacobians.resize(worksetSize, numCubPoints, pCellDim, pCellDim);
00684   
00685   // Step 5.b: Compute Jacobians at Gauss pts. on reference face for all parent cells:
00686   CellTools::setJacobian(worksetJacobians,
00687                          refFacePoints,
00688                          hexNodesAlt,
00689                          hexahedron_8);
00690   
00691   
00692   /*
00693    * 6. Get the (non-normalized) face normals for the face workset directly
00694    */
00695   // Step 6.a: Allocate storage for face normals
00696   FieldContainer<double> faceNormals(worksetSize, numCubPoints, pCellDim);
00697   
00698   // Step 6.b: Compute the face normals
00699   CellTools::getPhysicalFaceNormals(faceNormals,
00700                                     worksetJacobians,
00701                                     subcellOrd,
00702                                     hexahedron_8);
00703   
00704   // Step 6.c: Print face normals at Gauss points on workset faces (these Gauss points were computed in Example 5)
00705   std::cout 
00706     << "Face normals computed by CellTools::getPhysicalFaceNormals\n"
00707     << "Local face ordinal = " << subcellOrd <<"\n";
00708 
00709   for(int pCell = 0; pCell < worksetSize; pCell++){
00710     
00711     CellTools::printWorksetSubcell(hexNodesAlt, hexahedron_8, pCell, subcellDim, subcellOrd);
00712     
00713     for(int pt = 0; pt < numCubPoints; pt++){
00714       std::cout << "\t 3D Gauss point: (" 
00715       << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 0) << ", " 
00716       << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 1) << ", " 
00717       << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 2) << ")  " 
00718       << std::setw(8) << " outer normal:  " << "(" 
00719       << std::setw(8) << std::right << faceNormals(pCell, pt, 0) << ", " 
00720       << std::setw(8) << std::right << faceNormals(pCell, pt, 1) << ", " 
00721       << std::setw(8) << std::right << faceNormals(pCell, pt, 2) << ")\n";
00722     }    
00723     std::cout << "\n";      
00724   }//pCell
00725     
00726   /*
00727    * 7. Get the (non-normalized) face normals for the face workset using the face tangents. This may
00728    *    be useful if, for whatever reason,  face tangents are needed independently 
00729    */
00730   // Step 7.a: Allocate storage for face tangents
00731   FieldContainer<double> uFaceTan(worksetSize, numCubPoints, pCellDim);
00732   FieldContainer<double> vFaceTan(worksetSize, numCubPoints, pCellDim);
00733   
00734   // Step 7.b: Compute face tangents
00735   CellTools::getPhysicalFaceTangents(uFaceTan,
00736                                      vFaceTan,
00737                                      worksetJacobians,
00738                                      subcellOrd,
00739                                      hexahedron_8);
00740   
00741   // Step 7.c: Face outer normals (relative to parent cell) are uTan x vTan:
00742   RealSpaceTools<double>::vecprod(faceNormals, uFaceTan, vFaceTan);
00743   
00744   // Step 7.d: Print face normals at Gauss points on workset faces (these points were computed in Example 7)
00745   std::cout << "Face normals computed by CellTools::getPhysicalFaceTangents followed by vecprod\n";
00746   for(int pCell = 0; pCell < worksetSize; pCell++){
00747     
00748     CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd);
00749     
00750     for(int pt = 0; pt < numCubPoints; pt++){
00751       std::cout << "\t 3D Gauss point: (" 
00752       << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 0) << ", " 
00753       << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 1) << ", " 
00754       << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 2) << ")  " 
00755       << std::setw(8) << " outer normal:  " << "(" 
00756       << std::setw(8) << std::right << faceNormals(pCell, pt, 0) << ", " 
00757       << std::setw(8) << std::right << faceNormals(pCell, pt, 1) << ", " 
00758       << std::setw(8) << std::right << faceNormals(pCell, pt, 2) << ")\n";
00759     }    
00760     std::cout << "\n";      
00761   }//pCell
00762   
00763   
00764   std::cout << "\n" \
00765     << "===============================================================================\n"\
00766     << "| EXAMPLE 8: Side normals at Gauss points on a Hexahedron side workset        |\n"\
00767     << "===============================================================================\n";
00768   
00769   /* For this task we reuse the edge workset from Example 7 as a side workset and compute
00770    * normals to the sides in that set. This task repeats the first 5 steps from Example 3
00771    * The only difference is that in Step 6 we call CellTools::getPhysicalSideNormals
00772    */
00773   
00774   /*
00775    * 6. Get the (non-normalized) side normals for the side (face) workset:
00776    */
00777   // Step 6.a: Allocate storage for side normals
00778   sideNormals.resize(worksetSize, numCubPoints, pCellDim);
00779   
00780   // Step 6.b: Compute the side normals:
00781   CellTools::getPhysicalSideNormals(sideNormals,
00782                                     worksetJacobians,
00783                                     subcellOrd,
00784                                     hexahedron_8); 
00785   
00786   // Step 7.d: Print side normals at Gauss points on workset sides (these points were computed in Example 7)
00787   std::cout << "Side normals computed by CellTools::getPhysicalSideNormals\n";
00788   for(int pCell = 0; pCell < worksetSize; pCell++){
00789     
00790     CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd);
00791     
00792     for(int pt = 0; pt < numCubPoints; pt++){
00793       std::cout << "\t 3D Gauss point: (" 
00794       << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 0) << ", " 
00795       << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 1) << ", " 
00796       << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 2) << ")  " 
00797       << std::setw(8) << " side normal:  " << "(" 
00798       << std::setw(8) << std::right << sideNormals(pCell, pt, 0) << ", " 
00799       << std::setw(8) << std::right << sideNormals(pCell, pt, 1) << ", " 
00800       << std::setw(8) << std::right << sideNormals(pCell, pt, 2) << ")\n";
00801     }    
00802     std::cout << "\n";      
00803   }//pCell
00804   
00805   return 0;
00806 }