52#include "Shards_CellTopology.hpp"
53#include "Teuchos_GlobalMPISession.hpp"
56using namespace Intrepid;
58int main(
int argc,
char *argv[]) {
60 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
63 typedef shards::CellTopology CellTopology;
66 <<
"===============================================================================\n" \
68 <<
"| Example use of the CellTools class |\n" \
70 <<
"| 1) Definition of integration points on edge and face worksets |\n" \
71 <<
"| 2) Computation of face normals and edge tangents on face and edge worksets |\n" \
72 <<
"| 2) Computation of side normals on face and edge worksets |\n" \
74 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \
75 <<
"| Denis Ridzal (dridzal@sandia.gov), or |\n" \
76 <<
"| Kara Peterson (kjpeter@sandia.gov) |\n" \
78 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
79 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
81 <<
"===============================================================================\n"\
82 <<
"| EXAMPLE 1: Definition of integration points on a Triangle edge workset |\n"\
83 <<
"===============================================================================\n";
90 DefaultCubatureFactory<double> cubFactory;
93 CellTopology paramEdge(shards::getCellTopologyData<shards::Line<2> >() );
96 FieldContainer<double> paramEdgeWeights;
97 FieldContainer<double> paramEdgePoints;
100 FieldContainer<double> refEdgePoints;
103 FieldContainer<double> worksetEdgePoints;
112 CellTopology triangle_3( shards::getCellTopologyData<shards::Triangle<3> >() );
116 int pCellNodeCount = triangle_3.getVertexCount();
117 int pCellDim = triangle_3.getDimension();
119 FieldContainer<double> triNodes(worksetSize, pCellNodeCount, pCellDim);
121 triNodes(0, 0, 0) = 0.5; triNodes(0, 0, 1) = 2.0;
122 triNodes(0, 1, 0) = 0.0; triNodes(0, 1, 1) = 1.0;
123 triNodes(0, 2, 0) = 1.0; triNodes(0, 2, 1) = 1.0;
125 triNodes(1, 0, 0) = 1.0; triNodes(1, 0, 1) = 1.0;
126 triNodes(1, 1, 0) = 1.0; triNodes(1, 1, 1) = 0.0;
127 triNodes(1, 2, 0) = 0.0; triNodes(1, 2, 1) = 0.0;
129 triNodes(2, 0, 0) = 0.0; triNodes(2, 0, 1) = 0.0;
130 triNodes(2, 1, 0) = 1.0; triNodes(2, 1, 1) = 1.0;
131 triNodes(2, 2, 0) = 0.0; triNodes(2, 2, 1) = 1.0;
133 triNodes(3, 0, 0) = 1.0; triNodes(3, 0, 1) = 1.0;
134 triNodes(3, 1, 0) = 2.5; triNodes(3, 1, 1) = 1.5;
135 triNodes(3, 2, 0) = 0.5; triNodes(3, 2, 1) = 2.0;
148 Teuchos::RCP<Cubature<double> > triEdgeCubature = cubFactory.create(paramEdge, 6);
151 int cubDim = triEdgeCubature -> getDimension();
152 int numCubPoints = triEdgeCubature -> getNumPoints();
155 paramEdgePoints.resize(numCubPoints, cubDim);
156 paramEdgeWeights.resize(numCubPoints);
157 triEdgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights);
166 refEdgePoints.resize(numCubPoints, pCellDim);
169 worksetEdgePoints.resize(worksetSize, numCubPoints, pCellDim);
189 for(
int pCell = 0; pCell < worksetSize; pCell++){
193 for(
int pt = 0; pt < numCubPoints; pt++){
194 std::cout <<
"\t 1D Gauss point "
195 << std::setw(12) << std::right << paramEdgePoints(pt, 0) << std::setw(10) <<
" --> " <<
"("
196 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 0) <<
", "
197 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 1) <<
")\n";
205 <<
"===============================================================================\n"\
206 <<
"| EXAMPLE 2: Definition of integration points on a Quadrilateral edge workset|\n"\
207 <<
"===============================================================================\n";
215 CellTopology quadrilateral_4( shards::getCellTopologyData<shards::Quadrilateral<4> >() );
219 pCellNodeCount = quadrilateral_4.getVertexCount();
220 pCellDim = quadrilateral_4.getDimension();
222 FieldContainer<double> quadNodes(worksetSize, pCellNodeCount, pCellDim);
224 quadNodes(0, 0, 0) = 1.00; quadNodes(0, 0, 1) = 1.00;
225 quadNodes(0, 1, 0) = 2.00; quadNodes(0, 1, 1) = 0.75;
226 quadNodes(0, 2, 0) = 1.75; quadNodes(0, 2, 1) = 2.00;
227 quadNodes(0, 3, 0) = 1.25; quadNodes(0, 3, 1) = 2.00;
229 quadNodes(1, 0, 0) = 2.00; quadNodes(1, 0, 1) = 0.75;
230 quadNodes(1, 1, 0) = 3.00; quadNodes(1, 1, 1) = 1.25;
231 quadNodes(1, 2, 0) = 2.75; quadNodes(1, 2, 1) = 2.25;
232 quadNodes(1, 3, 0) = 1.75; quadNodes(1, 3, 1) = 2.00;
243 Teuchos::RCP<Cubature<double> > quadEdgeCubature = cubFactory.create(paramEdge, 6);
246 cubDim = quadEdgeCubature -> getDimension();
247 numCubPoints = quadEdgeCubature -> getNumPoints();
250 paramEdgePoints.resize(numCubPoints, cubDim);
251 paramEdgeWeights.resize(numCubPoints);
252 quadEdgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights);
259 refEdgePoints.resize(numCubPoints, pCellDim);
262 worksetEdgePoints.resize(worksetSize, numCubPoints, pCellDim);
280 for(
int pCell = 0; pCell < worksetSize; pCell++){
284 for(
int pt = 0; pt < numCubPoints; pt++){
285 std::cout <<
"\t 1D Gauss point "
286 << std::setw(12) << std::right << paramEdgePoints(pt, 0) << std::setw(10) <<
" --> " <<
"("
287 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 0) <<
", "
288 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 1) <<
")\n";
296 <<
"===============================================================================\n"\
297 <<
"| EXAMPLE 3: Edge tangents at Gauss points on a Quadrilateral edge workset |\n"\
298 <<
"===============================================================================\n";
317 FieldContainer<double> worksetJacobians(worksetSize, numCubPoints, pCellDim, pCellDim);
320 CellTools::setJacobian(worksetJacobians,
328 FieldContainer<double> edgeTangents(worksetSize, numCubPoints, pCellDim);
338 <<
"Edge tangents computed by CellTools::getPhysicalEdgeTangents.\n"
339 <<
"Local edge ordinal = " << subcellOrd <<
"\n";
341 for(
int pCell = 0; pCell < worksetSize; pCell++){
345 for(
int pt = 0; pt < numCubPoints; pt++){
346 std::cout <<
"\t 2D Gauss point: ("
347 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 0) <<
", "
348 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 1) <<
") "
349 << std::setw(8) <<
" edge tangent: " <<
"("
350 << std::setw(8) << std::right << edgeTangents(pCell, pt, 0) <<
", "
351 << std::setw(8) << std::right << edgeTangents(pCell, pt, 1) <<
")\n";
357 <<
"===============================================================================\n"\
358 <<
"| EXAMPLE 4: Side normals at Gauss points on a Quadrilateral side workset |\n"\
359 <<
"===============================================================================\n";
369 FieldContainer<double> sideNormals(worksetSize, numCubPoints, pCellDim);
379 <<
"Side normals computed by CellTools::getPhysicalSideNormals.\n"
380 <<
"Local side (edge) ordinal = " << subcellOrd <<
"\n";
382 for(
int pCell = 0; pCell < worksetSize; pCell++){
386 for(
int pt = 0; pt < numCubPoints; pt++){
387 std::cout <<
"\t 2D Gauss point: ("
388 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 0) <<
", "
389 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 1) <<
") "
390 << std::setw(8) <<
" side normal: " <<
"("
391 << std::setw(8) << std::right << sideNormals(pCell, pt, 0) <<
", "
392 << std::setw(8) << std::right << sideNormals(pCell, pt, 1) <<
")\n";
399 <<
"===============================================================================\n"\
400 <<
"| EXAMPLE 5: Definition of integration points on a Hexahedron edge workset |\n"\
401 <<
"===============================================================================\n";
409 CellTopology hexahedron_8( shards::getCellTopologyData<shards::Hexahedron<8> >() );
413 pCellNodeCount = hexahedron_8.getVertexCount();
414 pCellDim = hexahedron_8.getDimension();
416 FieldContainer<double> hexNodes(worksetSize, pCellNodeCount, pCellDim);
418 hexNodes(0, 0, 0) = 0.00; hexNodes(0, 0, 1) = 0.00, hexNodes(0, 0, 2) = 0.00;
419 hexNodes(0, 1, 0) = 1.00; hexNodes(0, 1, 1) = 0.00, hexNodes(0, 1, 2) = 0.00;
420 hexNodes(0, 2, 0) = 1.00; hexNodes(0, 2, 1) = 1.00, hexNodes(0, 2, 2) = 0.00;
421 hexNodes(0, 3, 0) = 0.00; hexNodes(0, 3, 1) = 1.00, hexNodes(0, 3, 2) = 0.00;
423 hexNodes(0, 4, 0) = 0.00; hexNodes(0, 4, 1) = 0.00, hexNodes(0, 4, 2) = 1.00;
424 hexNodes(0, 5, 0) = 1.00; hexNodes(0, 5, 1) = 0.00, hexNodes(0, 5, 2) = 1.00;
425 hexNodes(0, 6, 0) = 1.00; hexNodes(0, 6, 1) = 1.00, hexNodes(0, 6, 2) = 0.75;
426 hexNodes(0, 7, 0) = 0.00; hexNodes(0, 7, 1) = 1.00, hexNodes(0, 7, 2) = 1.00;
431 FieldContainer<double> hexNodesAlt(worksetSize, pCellNodeCount, pCellDim);
433 hexNodesAlt(0, 0, 0) = 0.00; hexNodesAlt(0, 0, 1) = 0.00, hexNodesAlt(0, 0, 2) = 0.00;
434 hexNodesAlt(0, 1, 0) = 1.00; hexNodesAlt(0, 1, 1) = 0.00, hexNodesAlt(0, 1, 2) = 0.00;
435 hexNodesAlt(0, 2, 0) = 1.00; hexNodesAlt(0, 2, 1) = 1.00, hexNodesAlt(0, 2, 2) = 0.00;
436 hexNodesAlt(0, 3, 0) = 0.00; hexNodesAlt(0, 3, 1) = 1.00, hexNodesAlt(0, 3, 2) = 0.00;
438 hexNodesAlt(0, 4, 0) = 0.00; hexNodesAlt(0, 4, 1) = 0.00, hexNodesAlt(0, 4, 2) = 1.00;
439 hexNodesAlt(0, 5, 0) = 1.00; hexNodesAlt(0, 5, 1) = 0.00, hexNodesAlt(0, 5, 2) = 0.75;
440 hexNodesAlt(0, 6, 0) = 1.00; hexNodesAlt(0, 6, 1) = 1.00, hexNodesAlt(0, 6, 2) = 0.50;
441 hexNodesAlt(0, 7, 0) = 0.00; hexNodesAlt(0, 7, 1) = 1.00, hexNodesAlt(0, 7, 2) = 0.75;
452 Teuchos::RCP<Cubature<double> > hexEdgeCubature = cubFactory.create(paramEdge, 4);
455 cubDim = hexEdgeCubature -> getDimension();
456 numCubPoints = hexEdgeCubature -> getNumPoints();
459 paramEdgePoints.resize(numCubPoints, cubDim);
460 paramEdgeWeights.resize(numCubPoints);
461 hexEdgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights);
468 refEdgePoints.resize(numCubPoints, pCellDim);
471 worksetEdgePoints.resize(worksetSize, numCubPoints, pCellDim);
489 for(
int pCell = 0; pCell < worksetSize; pCell++){
493 for(
int pt = 0; pt < numCubPoints; pt++){
494 std::cout <<
"\t 1D Gauss point "
495 << std::setw(12) << std::right << paramEdgePoints(pt, 0) << std::setw(10) <<
" --> " <<
"("
496 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 0) <<
", "
497 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 1) <<
", "
498 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 2) <<
")\n";
506 <<
"===============================================================================\n"\
507 <<
"| EXAMPLE 6: Edge tangents at Gauss points on a Hexahedron edge workset |\n"\
508 <<
"===============================================================================\n";
527 worksetJacobians.resize(worksetSize, numCubPoints, pCellDim, pCellDim);
530 CellTools::setJacobian(worksetJacobians,
539 edgeTangents.resize(worksetSize, numCubPoints, pCellDim);
549 <<
"Edge tangents computed by CellTools::getPhysicalEdgeTangents.\n"
550 <<
"Local edge ordinal = " << subcellOrd <<
"\n";
552 for(
int pCell = 0; pCell < worksetSize; pCell++){
556 for(
int pt = 0; pt < numCubPoints; pt++){
557 std::cout <<
"\t 3D Gauss point: ("
558 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 0) <<
", "
559 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 1) <<
", "
560 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 2) <<
") "
561 << std::setw(8) <<
" edge tangent: " <<
"("
562 << std::setw(8) << std::right << edgeTangents(pCell, pt, 0) <<
", "
563 << std::setw(8) << std::right << edgeTangents(pCell, pt, 1) <<
", "
564 << std::setw(8) << std::right << edgeTangents(pCell, pt, 2) <<
")\n";
571 <<
"===============================================================================\n"\
572 <<
"| EXAMPLE 7: Definition of integration points on a Hexahedron face workset |\n"\
573 <<
"===============================================================================\n";
580 CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
583 FieldContainer<double> paramFaceWeights;
584 FieldContainer<double> paramFacePoints;
587 FieldContainer<double> refFacePoints;
590 FieldContainer<double> worksetFacePoints;
607 Teuchos::RCP<Cubature<double> > hexFaceCubature = cubFactory.create(paramQuadFace, 3);
610 cubDim = hexFaceCubature -> getDimension();
611 numCubPoints = hexFaceCubature -> getNumPoints();
614 paramFacePoints.resize(numCubPoints, cubDim);
615 paramFaceWeights.resize(numCubPoints);
616 hexFaceCubature -> getCubature(paramFacePoints, paramFaceWeights);
623 refFacePoints.resize(numCubPoints, pCellDim);
626 worksetFacePoints.resize(worksetSize, numCubPoints, pCellDim);
645 for(
int pCell = 0; pCell < worksetSize; pCell++){
649 for(
int pt = 0; pt < numCubPoints; pt++){
650 std::cout <<
"\t 2D Gauss point ("
651 << std::setw(8) << std::right << paramFacePoints(pt, 0) <<
", "
652 << std::setw(8) << std::right << paramFacePoints(pt, 1) <<
") "
653 << std::setw(8) <<
" --> " <<
"("
654 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 0) <<
", "
655 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 1) <<
", "
656 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 2) <<
")\n";
663 <<
"===============================================================================\n"\
664 <<
"| EXAMPLE 8: Face normals at Gauss points on a Hexahedron face workset |\n"\
665 <<
"===============================================================================\n";
683 worksetJacobians.resize(worksetSize, numCubPoints, pCellDim, pCellDim);
686 CellTools::setJacobian(worksetJacobians,
696 FieldContainer<double> faceNormals(worksetSize, numCubPoints, pCellDim);
706 <<
"Face normals computed by CellTools::getPhysicalFaceNormals\n"
707 <<
"Local face ordinal = " << subcellOrd <<
"\n";
709 for(
int pCell = 0; pCell < worksetSize; pCell++){
713 for(
int pt = 0; pt < numCubPoints; pt++){
714 std::cout <<
"\t 3D Gauss point: ("
715 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 0) <<
", "
716 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 1) <<
", "
717 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 2) <<
") "
718 << std::setw(8) <<
" outer normal: " <<
"("
719 << std::setw(8) << std::right << faceNormals(pCell, pt, 0) <<
", "
720 << std::setw(8) << std::right << faceNormals(pCell, pt, 1) <<
", "
721 << std::setw(8) << std::right << faceNormals(pCell, pt, 2) <<
")\n";
731 FieldContainer<double> uFaceTan(worksetSize, numCubPoints, pCellDim);
732 FieldContainer<double> vFaceTan(worksetSize, numCubPoints, pCellDim);
745 std::cout <<
"Face normals computed by CellTools::getPhysicalFaceTangents followed by vecprod\n";
746 for(
int pCell = 0; pCell < worksetSize; pCell++){
750 for(
int pt = 0; pt < numCubPoints; pt++){
751 std::cout <<
"\t 3D Gauss point: ("
752 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 0) <<
", "
753 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 1) <<
", "
754 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 2) <<
") "
755 << std::setw(8) <<
" outer normal: " <<
"("
756 << std::setw(8) << std::right << faceNormals(pCell, pt, 0) <<
", "
757 << std::setw(8) << std::right << faceNormals(pCell, pt, 1) <<
", "
758 << std::setw(8) << std::right << faceNormals(pCell, pt, 2) <<
")\n";
765 <<
"===============================================================================\n"\
766 <<
"| EXAMPLE 8: Side normals at Gauss points on a Hexahedron side workset |\n"\
767 <<
"===============================================================================\n";
778 sideNormals.resize(worksetSize, numCubPoints, pCellDim);
787 std::cout <<
"Side normals computed by CellTools::getPhysicalSideNormals\n";
788 for(
int pCell = 0; pCell < worksetSize; pCell++){
792 for(
int pt = 0; pt < numCubPoints; pt++){
793 std::cout <<
"\t 3D Gauss point: ("
794 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 0) <<
", "
795 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 1) <<
", "
796 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 2) <<
") "
797 << std::setw(8) <<
" side normal: " <<
"("
798 << std::setw(8) << std::right << sideNormals(pCell, pt, 0) <<
", "
799 << std::setw(8) << std::right << sideNormals(pCell, pt, 1) <<
", "
800 << std::setw(8) << std::right << sideNormals(pCell, pt, 2) <<
")\n";
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.