65int main(
int argc,
char *argv[]) {
67 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
71 typedef shards::CellTopology CellTopology;
74 <<
"===============================================================================\n" \
76 <<
"| Example use of the CellTools class |\n" \
78 <<
"| 1) Computation of face flux, for a given vector field, on a face workset |\n" \
79 <<
"| 2) Computation of edge circulation, for a given vector field, on a face |\n" \
82 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \
83 <<
"| Denis Ridzal (dridzal@sandia.gov), or |\n" \
84 <<
"| Kara Peterson (kjpeter@sandia.gov) |\n" \
86 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
87 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
89 <<
"===============================================================================\n"\
90 <<
"| EXAMPLE 1: Computation of face flux on a face workset |\n"\
91 <<
"===============================================================================\n";
115 CellTopology hexahedron_8( shards::getCellTopologyData<shards::Hexahedron<8> >() );
119 int pCellNodeCount = hexahedron_8.getVertexCount();
120 int pCellDim = hexahedron_8.getDimension();
122 FieldContainer<double> hexNodes(worksetSize, pCellNodeCount, pCellDim);
124 hexNodes(0, 0, 0) = 0.00; hexNodes(0, 0, 1) = 0.00, hexNodes(0, 0, 2) = 0.00;
125 hexNodes(0, 1, 0) = 1.00; hexNodes(0, 1, 1) = 0.00, hexNodes(0, 1, 2) = 0.00;
126 hexNodes(0, 2, 0) = 1.00; hexNodes(0, 2, 1) = 1.00, hexNodes(0, 2, 2) = 0.00;
127 hexNodes(0, 3, 0) = 0.00; hexNodes(0, 3, 1) = 1.00, hexNodes(0, 3, 2) = 0.00;
129 hexNodes(0, 4, 0) = 0.00; hexNodes(0, 4, 1) = 0.00, hexNodes(0, 4, 2) = 1.00;
130 hexNodes(0, 5, 0) = 1.00; hexNodes(0, 5, 1) = 0.00, hexNodes(0, 5, 2) = 1.00;
131 hexNodes(0, 6, 0) = 1.00; hexNodes(0, 6, 1) = 1.00, hexNodes(0, 6, 2) = 1.00;
132 hexNodes(0, 7, 0) = 0.00; hexNodes(0, 7, 1) = 1.00, hexNodes(0, 7, 2) = 1.00;
135 hexNodes(1, 0, 0) = 0.00; hexNodes(1, 0, 1) = 0.00, hexNodes(1, 0, 2) = 0.00;
136 hexNodes(1, 1, 0) = 1.00; hexNodes(1, 1, 1) = 0.00, hexNodes(1, 1, 2) = 0.00;
137 hexNodes(1, 2, 0) = 1.00; hexNodes(1, 2, 1) = 1.00, hexNodes(1, 2, 2) = 0.00;
138 hexNodes(1, 3, 0) = 0.00; hexNodes(1, 3, 1) = 1.00, hexNodes(1, 3, 2) = 0.00;
140 hexNodes(1, 4, 0) = 0.00; hexNodes(1, 4, 1) = 0.00, hexNodes(1, 4, 2) = 1.00;
141 hexNodes(1, 5, 0) = 1.00; hexNodes(1, 5, 1) = 0.00, hexNodes(1, 5, 2) = 1.00;
142 hexNodes(1, 6, 0) = 1.00; hexNodes(1, 6, 1) = 1.00, hexNodes(1, 6, 2) = 0.75;
143 hexNodes(1, 7, 0) = 0.00; hexNodes(1, 7, 1) = 1.00, hexNodes(1, 7, 2) = 1.00;
162 DefaultCubatureFactory<double> cubFactory;
165 CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
168 FieldContainer<double> paramGaussWeights;
169 FieldContainer<double> paramGaussPoints;
172 FieldContainer<double> refGaussPoints;
175 FieldContainer<double> worksetGaussPoints;
180 Teuchos::RCP<Cubature<double> > hexFaceCubature = cubFactory.create(paramQuadFace, 3);
183 int cubDim = hexFaceCubature -> getDimension();
184 int numCubPoints = hexFaceCubature -> getNumPoints();
187 paramGaussPoints.resize(numCubPoints, cubDim);
188 paramGaussWeights.resize(numCubPoints);
189 hexFaceCubature -> getCubature(paramGaussPoints, paramGaussWeights);
194 refGaussPoints.resize(numCubPoints, pCellDim);
197 worksetGaussPoints.resize(worksetSize, numCubPoints, pCellDim);
224 FieldContainer<double> worksetJacobians(worksetSize, numCubPoints, pCellDim, pCellDim);
227 CellTools::setJacobian(worksetJacobians,
235 FieldContainer<double> worksetFaceTu(worksetSize, numCubPoints, pCellDim);
236 FieldContainer<double> worksetFaceTv(worksetSize, numCubPoints, pCellDim);
237 FieldContainer<double> worksetFaceN(worksetSize, numCubPoints, pCellDim);
258 FieldContainer<double> worksetVFieldVals(worksetSize, numCubPoints, pCellDim);
261 for(
int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){
262 for(
int ptOrd = 0; ptOrd < numCubPoints; ptOrd++){
264 double x = worksetGaussPoints(pCellOrd, ptOrd, 0);
265 double y = worksetGaussPoints(pCellOrd, ptOrd, 1);
266 double z = worksetGaussPoints(pCellOrd, ptOrd, 2);
268 vField(worksetVFieldVals(pCellOrd, ptOrd, 0),
269 worksetVFieldVals(pCellOrd, ptOrd, 1),
270 worksetVFieldVals(pCellOrd, ptOrd, 2), x, y, z);
283 FieldContainer<double> worksetFieldDotNormal(worksetSize, numCubPoints);
289 FieldContainer<double> worksetFluxes(worksetSize);
294 for(
int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){
295 worksetFluxes(pCellOrd) = 0.0;
296 for(
int pt = 0; pt < numCubPoints; pt++ ){
297 worksetFluxes(pCellOrd) += worksetFieldDotNormal(pCellOrd, pt)*paramGaussWeights(pt);
301 std::cout <<
"Face fluxes on workset faces : \n\n";
302 for(
int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){
305 std::cout <<
" Flux = " << worksetFluxes(pCellOrd) <<
"\n\n";
319 <<
"===============================================================================\n" \
320 <<
"| Gauss points on workset faces: |\n" \
321 <<
"===============================================================================\n";
322 for(
int pCell = 0; pCell < worksetSize; pCell++){
326 for(
int pt = 0; pt < numCubPoints; pt++){
327 std::cout <<
"\t 2D Gauss point ("
328 << std::setw(8) << std::right << paramGaussPoints(pt, 0) <<
", "
329 << std::setw(8) << std::right << paramGaussPoints(pt, 1) <<
") "
330 << std::setw(8) <<
" --> " <<
"("
331 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 0) <<
", "
332 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 1) <<
", "
333 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 2) <<
")\n";
341 <<
"===============================================================================\n" \
342 <<
"| Face normals (non-unit) at Gauss points on workset faces: |\n" \
343 <<
"===============================================================================\n";
344 for(
int pCell = 0; pCell < worksetSize; pCell++){
348 for(
int pt = 0; pt < numCubPoints; pt++){
349 std::cout <<
"\t 3D Gauss point: ("
350 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 0) <<
", "
351 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 1) <<
", "
352 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 2) <<
")"
353 << std::setw(8) <<
" out. normal: " <<
"("
354 << std::setw(8) << std::right << worksetFaceN(pCell, pt, 0) <<
", "
355 << std::setw(8) << std::right << worksetFaceN(pCell, pt, 1) <<
", "
356 << std::setw(8) << std::right << worksetFaceN(pCell, pt, 2) <<
")\n";