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