52 #include "Teuchos_oblackholestream.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_GlobalMPISession.hpp"
55 #include "Teuchos_ScalarTraits.hpp"
58 using namespace Intrepid;
80 const shards::CellTopology& parentCell,
84 const Teuchos::RCP<std::ostream>& outStream);
87 int main(
int argc,
char *argv[]) {
89 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
92 typedef shards::CellTopology CellTopology;
95 int iprint = argc - 1;
97 Teuchos::RCP<std::ostream> outStream;
98 Teuchos::oblackholestream bhs;
101 outStream = Teuchos::rcp(&std::cout,
false);
103 outStream = Teuchos::rcp(&bhs,
false);
106 Teuchos::oblackholestream oldFormatState;
107 oldFormatState.copyfmt(std::cout);
110 <<
"===============================================================================\n" \
112 <<
"| Unit Test CellTools |\n" \
114 <<
"| 1) Edge parametrizations |\n" \
115 <<
"| 2) Face parametrizations |\n" \
116 <<
"| 3) Edge tangents |\n" \
117 <<
"| 4) Face tangents and normals |\n" \
119 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \
120 <<
"| Denis Ridzal (dridzal@sandia.gov), or |\n" \
121 <<
"| Kara Peterson (kjpeter@sandia.gov) |\n" \
123 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
124 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
126 <<
"===============================================================================\n";
139 simplex_2(0, 0) = 0.0; simplex_2(0, 1) = 0.0;
140 simplex_2(1, 0) = 1.0; simplex_2(1, 1) = 0.0;
141 simplex_2(2, 0) = 0.0; simplex_2(2, 1) = 1.0;
146 cube_2(0, 0) = -1.0; cube_2(0, 1) = -1.0;
147 cube_2(1, 0) = 1.0; cube_2(1, 1) = -1.0;
148 cube_2(2, 0) = 1.0; cube_2(2, 1) = 1.0;
149 cube_2(3, 0) = -1.0; cube_2(3, 1) = 1.0;
153 std::vector<shards::CellTopology> allTopologies;
154 shards::getTopologies(allTopologies);
164 <<
"===============================================================================\n"\
165 <<
"| Test 1: edge parametrizations: |\n"\
166 <<
"===============================================================================\n\n";
171 for(
int topoOrd = 0; topoOrd < (int)allTopologies.size(); topoOrd++){
174 if(allTopologies[topoOrd].getDimension() > 1 && CellTools::hasReferenceCell(allTopologies[topoOrd]) ){
175 *outStream <<
" Testing edge parametrization for " << allTopologies[topoOrd].getName() <<
"\n";
177 allTopologies[topoOrd],
188 <<
"===============================================================================\n"\
189 <<
"| Test 2: face parametrizations: |\n"\
190 <<
"===============================================================================\n\n";
195 for(
int topoOrd = 0; topoOrd < (int)allTopologies.size(); topoOrd++){
198 if(allTopologies[topoOrd].getDimension() > 2 && CellTools::hasReferenceCell(allTopologies[topoOrd]) ){
199 *outStream <<
" Testing face parametrization for cell topology " << allTopologies[topoOrd].getName() <<
"\n";
201 allTopologies[topoOrd],
217 std::vector<shards::CellTopology> standardBaseTopologies;
218 shards::getTopologies(standardBaseTopologies, 4, shards::STANDARD_CELL, shards::BASE_TOPOLOGY);
221 CellTopology paramEdge (shards::getCellTopologyData<shards::Line<2> >() );
222 CellTopology paramTriFace (shards::getCellTopologyData<shards::Triangle<3> >() );
223 CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
230 <<
"===============================================================================\n"\
231 <<
"| Test 3: edge tangents/normals for stand. cells with base topologies: |\n"\
232 <<
"===============================================================================\n\n";
234 std::vector<shards::CellTopology>::iterator cti;
237 Teuchos::RCP<Cubature<double> > edgeCubature = cubFactory.
create(paramEdge, 6);
238 int cubDim = edgeCubature -> getDimension();
239 int numCubPoints = edgeCubature -> getNumPoints();
244 edgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights);
248 for(cti = standardBaseTopologies.begin(); cti !=standardBaseTopologies.end(); ++cti){
251 if( ( (*cti).getDimension() >= 2) && ( (*cti).getKey() != shards::Pyramid<5>::key) ){
253 int cellDim = (*cti).getDimension();
254 int vCount = (*cti).getVertexCount();
256 CellTools::getReferenceSubcellVertices(refCellVertices, cellDim, 0, (*cti) );
258 *outStream <<
" Testing edge tangents";
259 if(cellDim == 2) { *outStream <<
" and normals"; }
260 *outStream <<
" for cell topology " << (*cti).getName() <<
"\n";
268 for(
int v = 0; v < vCount; v++){
269 for(
int d = 0; d < cellDim; d++){
270 double delta = Teuchos::ScalarTraits<double>::random()/8.0;
271 physCellVertices(0, v, d) = refCellVertices(v, d) + delta;
282 for(
int edgeOrd = 0; edgeOrd < (int)(*cti).getEdgeCount(); edgeOrd++){
289 CellTools::mapToReferenceSubcell(refEdgePoints, paramEdgePoints, 1, edgeOrd, (*cti) );
290 CellTools::setJacobian(edgePointsJacobians, refEdgePoints, physCellVertices, (*cti) );
291 CellTools::getPhysicalEdgeTangents(edgePointTangents, edgePointsJacobians, edgeOrd, (*cti));
299 int v0ord = (*cti).getNodeMap(1, edgeOrd, 0);
300 int v1ord = (*cti).getNodeMap(1, edgeOrd, 1);
302 for(
int pt = 0; pt < numCubPoints; pt++){
307 for(
int d = 0; d < cellDim; d++){
308 edgeBenchmarkTangents(d) = (physCellVertices(0, v1ord, d) - physCellVertices(0, v0ord, d))/2.0;
311 if( abs(edgeBenchmarkTangents(d) - edgePointTangents(0, pt, d)) > INTREPID_THRESHOLD ){
314 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
315 <<
" Edge tangent computation by CellTools failed for: \n"
316 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
317 <<
" Edge ordinal = " << edgeOrd <<
"\n"
318 <<
" Edge point number = " << pt <<
"\n"
319 <<
" Tangent coordinate = " << d <<
"\n"
320 <<
" CellTools value = " << edgePointTangents(0, pt, d) <<
"\n"
321 <<
" Benchmark value = " << edgeBenchmarkTangents(d) <<
"\n\n";
327 CellTools::getPhysicalSideNormals(edgePointNormals, edgePointsJacobians, edgeOrd, (*cti));
328 if( abs(edgeBenchmarkTangents(1) - edgePointNormals(0, pt, 0)) > INTREPID_THRESHOLD ){
331 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
332 <<
" Edge Normal computation by CellTools failed for: \n"
333 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
334 <<
" Edge ordinal = " << edgeOrd <<
"\n"
335 <<
" Edge point number = " << pt <<
"\n"
336 <<
" Normal coordinate = " << 0 <<
"\n"
337 <<
" CellTools value = " << edgePointNormals(0, pt, 0) <<
"\n"
338 <<
" Benchmark value = " << edgeBenchmarkTangents(1) <<
"\n\n";
340 if( abs(edgeBenchmarkTangents(0) + edgePointNormals(0, pt, 1)) > INTREPID_THRESHOLD ){
343 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
344 <<
" Edge Normal computation by CellTools failed for: \n"
345 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
346 <<
" Edge ordinal = " << edgeOrd <<
"\n"
347 <<
" Edge point number = " << pt <<
"\n"
348 <<
" Normal coordinate = " << 1 <<
"\n"
349 <<
" CellTools value = " << edgePointNormals(0, pt, 1) <<
"\n"
350 <<
" Benchmark value = " << -edgeBenchmarkTangents(0) <<
"\n\n";
362 <<
"===============================================================================\n"\
363 <<
"| Test 4: face/side normals for stand. 3D cells with base topologies: | |\n"\
364 <<
"===============================================================================\n\n";
368 Teuchos::RCP<Cubature<double> > triFaceCubature = cubFactory.
create(paramTriFace, 6);
369 Teuchos::RCP<Cubature<double> > quadFaceCubature = cubFactory.
create(paramQuadFace, 6);
371 int faceCubDim = triFaceCubature -> getDimension();
372 int numTriFaceCubPoints = triFaceCubature -> getNumPoints();
373 int numQuadFaceCubPoints = quadFaceCubature -> getNumPoints();
381 triFaceCubature -> getCubature(paramTriFacePoints, paramTriFaceWeights);
382 quadFaceCubature -> getCubature(paramQuadFacePoints, paramQuadFaceWeights);
386 for(cti = standardBaseTopologies.begin(); cti !=standardBaseTopologies.end(); ++cti){
389 if( ( (*cti).getDimension() == 3) && ( (*cti).getKey() != shards::Pyramid<5>::key) ){
391 int cellDim = (*cti).getDimension();
392 int vCount = (*cti).getVertexCount();
394 CellTools::getReferenceSubcellVertices(refCellVertices, cellDim, 0, (*cti) );
396 *outStream <<
" Testing face/side normals for cell topology " << (*cti).getName() <<
"\n";
403 for(
int v = 0; v < vCount; v++){
404 for(
int d = 0; d < cellDim; d++){
405 double delta = Teuchos::ScalarTraits<double>::random()/8.0;
406 physCellVertices(0, v, d) = refCellVertices(v, d) + delta;
423 for(
int faceOrd = 0; faceOrd < (int)(*cti).getSideCount(); faceOrd++){
427 switch( (*cti).getCellTopologyData(2, faceOrd) -> key ) {
429 case shards::Triangle<3>::key:
432 CellTools::mapToReferenceSubcell(refTriFacePoints, paramTriFacePoints, 2, faceOrd, (*cti) );
433 CellTools::setJacobian(triFacePointsJacobians, refTriFacePoints, physCellVertices, (*cti) );
434 CellTools::getPhysicalFaceNormals(triFacePointNormals, triFacePointsJacobians, faceOrd, (*cti));
435 CellTools::getPhysicalSideNormals(triSidePointNormals, triFacePointsJacobians, faceOrd, (*cti));
442 int v0ord = (*cti).getNodeMap(2, faceOrd, 0);
443 int v1ord = (*cti).getNodeMap(2, faceOrd, 1);
444 int v2ord = (*cti).getNodeMap(2, faceOrd, 2);
448 for(
int pt = 0; pt < numTriFaceCubPoints; pt++){
450 for(
int d = 0; d < cellDim; d++){
451 tanX(d) = (physCellVertices(0, v1ord, d) - physCellVertices(0, v0ord, d));
452 tanY(d) = (physCellVertices(0, v2ord, d) - physCellVertices(0, v0ord, d));
458 for(
int d = 0; d < cellDim; d++){
461 if( abs(faceNormal(d) - triFacePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){
464 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
465 <<
" Face normal computation by CellTools failed for: \n"
466 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
467 <<
" Face Topology = " << (*cti).getCellTopologyData(2, faceOrd) -> name <<
"\n"
468 <<
" Face ordinal = " << faceOrd <<
"\n"
469 <<
" Face point number = " << pt <<
"\n"
470 <<
" Normal coordinate = " << d <<
"\n"
471 <<
" CellTools value = " << triFacePointNormals(0, pt, d)
472 <<
" Benchmark value = " << faceNormal(d) <<
"\n\n";
475 if( abs(faceNormal(d) - triSidePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){
478 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
479 <<
" Side normal computation by CellTools failed for: \n"
480 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
481 <<
" Side Topology = " << (*cti).getCellTopologyData(2, faceOrd) -> name <<
"\n"
482 <<
" Side ordinal = " << faceOrd <<
"\n"
483 <<
" Side point number = " << pt <<
"\n"
484 <<
" Normal coordinate = " << d <<
"\n"
485 <<
" CellTools value = " << triSidePointNormals(0, pt, d)
486 <<
" Benchmark value = " << faceNormal(d) <<
"\n\n";
493 case shards::Quadrilateral<4>::key:
496 CellTools::mapToReferenceSubcell(refQuadFacePoints, paramQuadFacePoints, 2, faceOrd, (*cti) );
497 CellTools::setJacobian(quadFacePointsJacobians, refQuadFacePoints, physCellVertices, (*cti) );
498 CellTools::getPhysicalFaceNormals(quadFacePointNormals, quadFacePointsJacobians, faceOrd, (*cti));
499 CellTools::getPhysicalSideNormals(quadSidePointNormals, quadFacePointsJacobians, faceOrd, (*cti));
508 int v0ord = (*cti).getNodeMap(2, faceOrd, 0);
509 int v1ord = (*cti).getNodeMap(2, faceOrd, 1);
510 int v2ord = (*cti).getNodeMap(2, faceOrd, 2);
511 int v3ord = (*cti).getNodeMap(2, faceOrd, 3);
514 for(
int pt = 0; pt < numTriFaceCubPoints; pt++){
516 for(
int d = 0; d < cellDim; d++){
517 tanX(d) = (physCellVertices(0, v0ord, d)*(-1.0 + paramQuadFacePoints(pt,1) ) +
518 physCellVertices(0, v1ord, d)*( 1.0 - paramQuadFacePoints(pt,1) ) +
519 physCellVertices(0, v2ord, d)*( 1.0 + paramQuadFacePoints(pt,1) ) +
520 physCellVertices(0, v3ord, d)*(-1.0 - paramQuadFacePoints(pt,1) ) )/4.0;
522 tanY(d) = (physCellVertices(0, v0ord, d)*(-1.0 + paramQuadFacePoints(pt,0) ) +
523 physCellVertices(0, v1ord, d)*(-1.0 - paramQuadFacePoints(pt,0) ) +
524 physCellVertices(0, v2ord, d)*( 1.0 + paramQuadFacePoints(pt,0) ) +
525 physCellVertices(0, v3ord, d)*( 1.0 - paramQuadFacePoints(pt,0) ) )/4.0;
530 for(
int d = 0; d < cellDim; d++){
533 if( abs(faceNormal(d) - quadFacePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){
536 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
537 <<
" Face normal computation by CellTools failed for: \n"
538 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
539 <<
" Face Topology = " << (*cti).getCellTopologyData(2, faceOrd) -> name <<
"\n"
540 <<
" Face ordinal = " << faceOrd <<
"\n"
541 <<
" Face point number = " << pt <<
"\n"
542 <<
" Normal coordinate = " << d <<
"\n"
543 <<
" CellTools value = " << quadFacePointNormals(0, pt, d)
544 <<
" Benchmark value = " << faceNormal(d) <<
"\n\n";
547 if( abs(faceNormal(d) - quadSidePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){
550 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
551 <<
" Side normal computation by CellTools failed for: \n"
552 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
553 <<
" Side Topology = " << (*cti).getCellTopologyData(2, faceOrd) -> name <<
"\n"
554 <<
" Side ordinal = " << faceOrd <<
"\n"
555 <<
" Side point number = " << pt <<
"\n"
556 <<
" Normal coordinate = " << d <<
"\n"
557 <<
" CellTools value = " << quadSidePointNormals(0, pt, d)
558 <<
" Benchmark value = " << faceNormal(d) <<
"\n\n";
566 *outStream <<
" Face normals test failure: face topology not supported \n\n";
576 catch (std::logic_error err) {
577 *outStream << err.what() <<
"\n";
583 std::cout <<
"End Result: TEST FAILED\n";
585 std::cout <<
"End Result: TEST PASSED\n";
588 std::cout.copyfmt(oldFormatState);
596 const shards::CellTopology& parentCell,
600 const Teuchos::RCP<std::ostream>& outStream){
603 int cellDim = parentCell.getDimension();
604 int subcCount = parentCell.getSubcellCount(subcDim);
608 for(
int subcOrd = 0; subcOrd < subcCount; subcOrd++){
609 int subcVertexCount = parentCell.getVertexCount(subcDim, subcOrd);
631 else if(subcDim == 2) {
634 if(subcVertexCount == 3){
642 else if(subcVertexCount == 4){
652 for(
int subcVertOrd = 0; subcVertOrd < subcVertexCount; subcVertOrd++){
653 for(
int dim = 0; dim < cellDim; dim++){
655 if(mappedParamVertices(subcVertOrd, dim) != refSubcellVertices(subcVertOrd, dim) ) {
658 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
659 <<
" Cell Topology = " << parentCell.getName() <<
"\n"
660 <<
" Parametrization of subcell " << subcOrd <<
" which is "
661 << parentCell.getName(subcDim,subcOrd) <<
" failed for vertex " << subcVertOrd <<
":\n"
662 <<
" parametrization map fails to map correctly coordinate " << dim <<
" of that vertex\n\n";
Header file for utility class to provide multidimensional containers.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
A factory class that generates specific instances of cubatures.
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > °ree)
Factory method.
void testSubcellParametrizations(int &errorFlag, const shards::CellTopology &parentCell, const FieldContainer< double > &subcParamVert_A, const FieldContainer< double > &subcParamVert_B, const int subcDim, const Teuchos::RCP< std::ostream > &outStream)
Maps the vertices of the subcell parametrization domain to that subcell.