|
Intrepid
|
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 }
1.7.6.1