52 #include "Teuchos_oblackholestream.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_GlobalMPISession.hpp"
55 #include "Teuchos_ScalarTraits.hpp"
56 #include <Kokkos_Core.hpp>
59 using namespace Intrepid;
81 const shards::CellTopology& parentCell,
82 const Kokkos::View<double**>& subcParamVert_A,
83 const Kokkos::View<double**>& subcParamVert_B,
85 const Teuchos::RCP<std::ostream>& outStream);
88 int main(
int argc,
char *argv[]) {
90 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
93 typedef shards::CellTopology CellTopology;
96 int iprint = argc - 1;
98 Teuchos::RCP<std::ostream> outStream;
99 Teuchos::oblackholestream bhs;
102 outStream = Teuchos::rcp(&std::cout,
false);
104 outStream = Teuchos::rcp(&bhs,
false);
107 Teuchos::oblackholestream oldFormatState;
108 oldFormatState.copyfmt(std::cout);
111 <<
"===============================================================================\n" \
113 <<
"| Unit Test CellTools |\n" \
115 <<
"| 1) Edge parametrizations |\n" \
116 <<
"| 2) Face parametrizations |\n" \
117 <<
"| 3) Edge tangents |\n" \
118 <<
"| 4) Face tangents and normals |\n" \
120 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \
121 <<
"| Denis Ridzal (dridzal@sandia.gov), or |\n" \
122 <<
"| Kara Peterson (kjpeter@sandia.gov) |\n" \
124 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
125 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
127 <<
"===============================================================================\n";
133 Kokkos::View<double**> cube_1(
"cube_1",2, 1);
139 Kokkos::View<double**> simplex_2(
"simplex_2",3, 2);
140 simplex_2(0, 0) = 0.0; simplex_2(0, 1) = 0.0;
141 simplex_2(1, 0) = 1.0; simplex_2(1, 1) = 0.0;
142 simplex_2(2, 0) = 0.0; simplex_2(2, 1) = 1.0;
146 Kokkos::View<double**> cube_2(
"cube_2",4, 2);
147 cube_2(0, 0) = -1.0; cube_2(0, 1) = -1.0;
148 cube_2(1, 0) = 1.0; cube_2(1, 1) = -1.0;
149 cube_2(2, 0) = 1.0; cube_2(2, 1) = 1.0;
150 cube_2(3, 0) = -1.0; cube_2(3, 1) = 1.0;
154 std::vector<shards::CellTopology> allTopologies;
155 shards::getTopologies(allTopologies);
165 <<
"===============================================================================\n"\
166 <<
"| Test 1: edge parametrizations: |\n"\
167 <<
"===============================================================================\n\n";
172 for(
int topoOrd = 0; topoOrd < (int)allTopologies.size(); topoOrd++){
175 if(allTopologies[topoOrd].getDimension() > 1 && CellTools::hasReferenceCell(allTopologies[topoOrd]) ){
176 *outStream <<
" Testing edge parametrization for " << allTopologies[topoOrd].getName() <<
"\n";
178 allTopologies[topoOrd],
189 <<
"===============================================================================\n"\
190 <<
"| Test 2: face parametrizations: |\n"\
191 <<
"===============================================================================\n\n";
196 for(
int topoOrd = 0; topoOrd < (int)allTopologies.size(); topoOrd++){
199 if(allTopologies[topoOrd].getDimension() > 2 && CellTools::hasReferenceCell(allTopologies[topoOrd]) ){
200 *outStream <<
" Testing face parametrization for cell topology " << allTopologies[topoOrd].getName() <<
"\n";
202 allTopologies[topoOrd],
218 std::vector<shards::CellTopology> standardBaseTopologies;
219 shards::getTopologies(standardBaseTopologies, 4, shards::STANDARD_CELL, shards::BASE_TOPOLOGY);
222 CellTopology paramEdge (shards::getCellTopologyData<shards::Line<2> >() );
223 CellTopology paramTriFace (shards::getCellTopologyData<shards::Triangle<3> >() );
224 CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
231 <<
"===============================================================================\n"\
232 <<
"| Test 3: edge tangents/normals for stand. cells with base topologies: |\n"\
233 <<
"===============================================================================\n\n";
235 std::vector<shards::CellTopology>::iterator cti;
238 Teuchos::RCP<Cubature<double> > edgeCubature = cubFactory.
create(paramEdge, 6);
239 int cubDim = edgeCubature -> getDimension();
240 int numCubPoints = edgeCubature -> getNumPoints();
245 edgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights);
249 for(cti = standardBaseTopologies.begin(); cti !=standardBaseTopologies.end(); ++cti){
252 if( ( (*cti).getDimension() >= 2) && ( (*cti).getKey() != shards::Pyramid<5>::key) ){
254 int cellDim = (*cti).getDimension();
255 int vCount = (*cti).getVertexCount();
256 Kokkos::View<double**> refCellVertices(
"refCellVertices",vCount, cellDim);
257 CellTools::getReferenceSubcellVertices(refCellVertices, cellDim, 0, (*cti) );
259 *outStream <<
" Testing edge tangents";
260 if(cellDim == 2) { *outStream <<
" and normals"; }
261 *outStream <<
" for cell topology " << (*cti).getName() <<
"\n";
265 Kokkos::View<double***> physCellVertices(
"physCellVertices",1, vCount, cellDim);
269 for(
int v = 0; v < vCount; v++){
270 for(
int d = 0; d < cellDim; d++){
271 double delta = Teuchos::ScalarTraits<double>::random()/8.0;
272 physCellVertices(0, v, d) = refCellVertices(v, d) + delta;
277 Kokkos::View<double**> refEdgePoints(
"refEdgePoints",numCubPoints, cellDim);
278 Kokkos::View<double****> edgePointsJacobians(
"edgePointsJacobians",1, numCubPoints, cellDim, cellDim);
279 Kokkos::View<double***> edgePointTangents(
"edgePointTangents",1, numCubPoints, cellDim);
280 Kokkos::View<double***> edgePointNormals(
"edgePointNormals",1, numCubPoints, cellDim);
283 for(
int edgeOrd = 0; edgeOrd < (int)(*cti).getEdgeCount(); edgeOrd++){
290 CellTools::mapToReferenceSubcell(refEdgePoints, paramEdgePoints, 1, edgeOrd, (*cti) );
291 CellTools::setJacobian(edgePointsJacobians, refEdgePoints, physCellVertices, (*cti) );
292 CellTools::getPhysicalEdgeTangents(edgePointTangents, edgePointsJacobians, edgeOrd, (*cti));
300 int v0ord = (*cti).getNodeMap(1, edgeOrd, 0);
301 int v1ord = (*cti).getNodeMap(1, edgeOrd, 1);
303 for(
int pt = 0; pt < numCubPoints; pt++){
306 Kokkos::View<double*> edgeBenchmarkTangents(
"edgeBenchmarkTangents",3);
308 for(
int d = 0; d < cellDim; d++){
309 edgeBenchmarkTangents(d) = (physCellVertices(0, v1ord, d) - physCellVertices(0, v0ord, d))/2.0;
312 if( abs(edgeBenchmarkTangents(d) - edgePointTangents(0, pt, d)) > INTREPID_THRESHOLD ){
315 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
316 <<
" Edge tangent computation by CellTools failed for: \n"
317 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
318 <<
" Edge ordinal = " << edgeOrd <<
"\n"
319 <<
" Edge point number = " << pt <<
"\n"
320 <<
" Tangent coordinate = " << d <<
"\n"
321 <<
" CellTools value = " << edgePointTangents(0, pt, d) <<
"\n"
322 <<
" Benchmark value = " << edgeBenchmarkTangents(d) <<
"\n\n";
328 CellTools::getPhysicalSideNormals(edgePointNormals, edgePointsJacobians, edgeOrd, (*cti));
329 if( abs(edgeBenchmarkTangents(1) - edgePointNormals(0, pt, 0)) > INTREPID_THRESHOLD ){
332 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
333 <<
" Edge Normal computation by CellTools failed for: \n"
334 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
335 <<
" Edge ordinal = " << edgeOrd <<
"\n"
336 <<
" Edge point number = " << pt <<
"\n"
337 <<
" Normal coordinate = " << 0 <<
"\n"
338 <<
" CellTools value = " << edgePointNormals(0, pt, 0) <<
"\n"
339 <<
" Benchmark value = " << edgeBenchmarkTangents(1) <<
"\n\n";
341 if( abs(edgeBenchmarkTangents(0) + edgePointNormals(0, pt, 1)) > INTREPID_THRESHOLD ){
344 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
345 <<
" Edge Normal computation by CellTools failed for: \n"
346 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
347 <<
" Edge ordinal = " << edgeOrd <<
"\n"
348 <<
" Edge point number = " << pt <<
"\n"
349 <<
" Normal coordinate = " << 1 <<
"\n"
350 <<
" CellTools value = " << edgePointNormals(0, pt, 1) <<
"\n"
351 <<
" Benchmark value = " << -edgeBenchmarkTangents(0) <<
"\n\n";
363 <<
"===============================================================================\n"\
364 <<
"| Test 4: face/side normals for stand. 3D cells with base topologies: | |\n"\
365 <<
"===============================================================================\n\n";
369 Teuchos::RCP<Cubature<double> > triFaceCubature = cubFactory.
create(paramTriFace, 6);
370 Teuchos::RCP<Cubature<double> > quadFaceCubature = cubFactory.
create(paramQuadFace, 6);
372 int faceCubDim = triFaceCubature -> getDimension();
373 int numTriFaceCubPoints = triFaceCubature -> getNumPoints();
374 int numQuadFaceCubPoints = quadFaceCubature -> getNumPoints();
382 triFaceCubature -> getCubature(paramTriFacePointsFC, paramTriFaceWeightsFC);
383 quadFaceCubature -> getCubature(paramQuadFacePointsFC, paramQuadFaceWeightsFC);
384 Kokkos::View<double**> paramTriFacePoints(¶mTriFacePointsFC[0],paramTriFacePointsFC.dimension(0),paramTriFacePointsFC.dimension(1));
385 Kokkos::View<double*> paramTriFaceWeights(¶mTriFaceWeightsFC[0],paramTriFaceWeightsFC.dimension(0));
386 Kokkos::View<double**> paramQuadFacePoints(¶mQuadFacePointsFC[0],paramQuadFacePointsFC.dimension(0),paramQuadFacePointsFC.dimension(1));
387 Kokkos::View<double*> paramQuadFaceWeights(¶mQuadFaceWeightsFC[0],paramQuadFaceWeightsFC.dimension(0));
389 for(cti = standardBaseTopologies.begin(); cti !=standardBaseTopologies.end(); ++cti){
392 if( ( (*cti).getDimension() == 3) && ( (*cti).getKey() != shards::Pyramid<5>::key) ){
394 int cellDim = (*cti).getDimension();
395 int vCount = (*cti).getVertexCount();
396 Kokkos::View<double**> refCellVertices(
"refCellVertices",vCount, cellDim);
397 CellTools::getReferenceSubcellVertices(refCellVertices, cellDim, 0, (*cti) );
399 *outStream <<
" Testing face/side normals for cell topology " << (*cti).getName() <<
"\n";
402 Kokkos::View<double***> physCellVertices(
"physCellVertices",1, vCount, cellDim);
406 for(
int v = 0; v < vCount; v++){
407 for(
int d = 0; d < cellDim; d++){
408 double delta = Teuchos::ScalarTraits<double>::random()/8.0;
409 physCellVertices(0, v, d) = refCellVertices(v, d) + delta;
415 Kokkos::View<double**> refTriFacePoints(
"refTriFacePoints",numTriFaceCubPoints, cellDim);
416 Kokkos::View<double**> refQuadFacePoints(
"refQuadFacePoints",numQuadFaceCubPoints, cellDim);
417 Kokkos::View<double****> triFacePointsJacobians(
"triFacePointsJacobians",1, numTriFaceCubPoints, cellDim, cellDim);
418 Kokkos::View<double****> quadFacePointsJacobians(
"quadFacePointsJacobians",1, numQuadFaceCubPoints, cellDim, cellDim);
419 Kokkos::View<double***> triFacePointNormals(
"triFacePointNormals",1, numTriFaceCubPoints, cellDim);
420 Kokkos::View<double***> triSidePointNormals(
"triSidePointNormals",1, numTriFaceCubPoints, cellDim);
421 Kokkos::View<double***> quadFacePointNormals(
"quadFacePointNormals",1, numQuadFaceCubPoints, cellDim);
422 Kokkos::View<double***> quadSidePointNormals(
"quadSidePointNormals",1, numQuadFaceCubPoints, cellDim);
426 for(
int faceOrd = 0; faceOrd < (int)(*cti).getSideCount(); faceOrd++){
430 switch( (*cti).getCellTopologyData(2, faceOrd) -> key ) {
432 case shards::Triangle<3>::key:
435 CellTools::mapToReferenceSubcell(refTriFacePoints, paramTriFacePoints, 2, faceOrd, (*cti) );
436 CellTools::setJacobian(triFacePointsJacobians, refTriFacePoints, physCellVertices, (*cti) );
437 CellTools::getPhysicalFaceNormals(triFacePointNormals, triFacePointsJacobians, faceOrd, (*cti));
438 CellTools::getPhysicalSideNormals(triSidePointNormals, triFacePointsJacobians, faceOrd, (*cti));
445 int v0ord = (*cti).getNodeMap(2, faceOrd, 0);
446 int v1ord = (*cti).getNodeMap(2, faceOrd, 1);
447 int v2ord = (*cti).getNodeMap(2, faceOrd, 2);
451 for(
int pt = 0; pt < numTriFaceCubPoints; pt++){
452 Kokkos::View<double*> tanX(
"tanX",3), tanY(
"tanY",3), faceNormal(
"faceNormal",3);
453 for(
int d = 0; d < cellDim; d++){
454 tanX(d) = (physCellVertices(0, v1ord, d) - physCellVertices(0, v0ord, d));
455 tanY(d) = (physCellVertices(0, v2ord, d) - physCellVertices(0, v0ord, d));
461 for(
int d = 0; d < cellDim; d++){
464 if( abs(faceNormal(d) - triFacePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){
467 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
468 <<
" Face normal computation by CellTools failed for: \n"
469 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
470 <<
" Face Topology = " << (*cti).getCellTopologyData(2, faceOrd) -> name <<
"\n"
471 <<
" Face ordinal = " << faceOrd <<
"\n"
472 <<
" Face point number = " << pt <<
"\n"
473 <<
" Normal coordinate = " << d <<
"\n"
474 <<
" CellTools value = " << triFacePointNormals(0, pt, d)
475 <<
" Benchmark value = " << faceNormal(d) <<
"\n\n";
478 if( abs(faceNormal(d) - triSidePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){
481 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
482 <<
" Side normal computation by CellTools failed for: \n"
483 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
484 <<
" Side Topology = " << (*cti).getCellTopologyData(2, faceOrd) -> name <<
"\n"
485 <<
" Side ordinal = " << faceOrd <<
"\n"
486 <<
" Side point number = " << pt <<
"\n"
487 <<
" Normal coordinate = " << d <<
"\n"
488 <<
" CellTools value = " << triSidePointNormals(0, pt, d)
489 <<
" Benchmark value = " << faceNormal(d) <<
"\n\n";
496 case shards::Quadrilateral<4>::key:
499 CellTools::mapToReferenceSubcell(refQuadFacePoints, paramQuadFacePoints, 2, faceOrd, (*cti) );
500 CellTools::setJacobian(quadFacePointsJacobians, refQuadFacePoints, physCellVertices, (*cti) );
501 CellTools::getPhysicalFaceNormals(quadFacePointNormals, quadFacePointsJacobians, faceOrd, (*cti));
502 CellTools::getPhysicalSideNormals(quadSidePointNormals, quadFacePointsJacobians, faceOrd, (*cti));
511 int v0ord = (*cti).getNodeMap(2, faceOrd, 0);
512 int v1ord = (*cti).getNodeMap(2, faceOrd, 1);
513 int v2ord = (*cti).getNodeMap(2, faceOrd, 2);
514 int v3ord = (*cti).getNodeMap(2, faceOrd, 3);
517 for(
int pt = 0; pt < numTriFaceCubPoints; pt++){
518 Kokkos::View<double*> tanX(
"tanX",3), tanY(
"tanY",3), faceNormal(
"faceNormal",3);
519 for(
int d = 0; d < cellDim; d++){
520 tanX(d) = (physCellVertices(0, v0ord, d)*(-1.0 + paramQuadFacePoints(pt,1) ) +
521 physCellVertices(0, v1ord, d)*( 1.0 - paramQuadFacePoints(pt,1) ) +
522 physCellVertices(0, v2ord, d)*( 1.0 + paramQuadFacePoints(pt,1) ) +
523 physCellVertices(0, v3ord, d)*(-1.0 - paramQuadFacePoints(pt,1) ) )/4.0;
525 tanY(d) = (physCellVertices(0, v0ord, d)*(-1.0 + paramQuadFacePoints(pt,0) ) +
526 physCellVertices(0, v1ord, d)*(-1.0 - paramQuadFacePoints(pt,0) ) +
527 physCellVertices(0, v2ord, d)*( 1.0 + paramQuadFacePoints(pt,0) ) +
528 physCellVertices(0, v3ord, d)*( 1.0 - paramQuadFacePoints(pt,0) ) )/4.0;
533 for(
int d = 0; d < cellDim; d++){
536 if( abs(faceNormal(d) - quadFacePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){
539 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
540 <<
" Face normal computation by CellTools failed for: \n"
541 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
542 <<
" Face Topology = " << (*cti).getCellTopologyData(2, faceOrd) -> name <<
"\n"
543 <<
" Face ordinal = " << faceOrd <<
"\n"
544 <<
" Face point number = " << pt <<
"\n"
545 <<
" Normal coordinate = " << d <<
"\n"
546 <<
" CellTools value = " << quadFacePointNormals(0, pt, d)
547 <<
" Benchmark value = " << faceNormal(d) <<
"\n\n";
550 if( abs(faceNormal(d) - quadSidePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){
553 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
554 <<
" Side normal computation by CellTools failed for: \n"
555 <<
" Cell Topology = " << (*cti).getName() <<
"\n"
556 <<
" Side Topology = " << (*cti).getCellTopologyData(2, faceOrd) -> name <<
"\n"
557 <<
" Side ordinal = " << faceOrd <<
"\n"
558 <<
" Side point number = " << pt <<
"\n"
559 <<
" Normal coordinate = " << d <<
"\n"
560 <<
" CellTools value = " << quadSidePointNormals(0, pt, d)
561 <<
" Benchmark value = " << faceNormal(d) <<
"\n\n";
569 *outStream <<
" Face normals test failure: face topology not supported \n\n";
579 catch (std::logic_error err) {
580 *outStream << err.what() <<
"\n";
586 std::cout <<
"End Result: TEST FAILED\n";
588 std::cout <<
"End Result: TEST PASSED\n";
591 std::cout.copyfmt(oldFormatState);
599 const shards::CellTopology& parentCell,
600 const Kokkos::View<double**>& subcParamVert_A,
601 const Kokkos::View<double**>& subcParamVert_B,
603 const Teuchos::RCP<std::ostream>& outStream){
606 int cellDim = parentCell.getDimension();
607 int subcCount = parentCell.getSubcellCount(subcDim);
611 for(
int subcOrd = 0; subcOrd < subcCount; subcOrd++){
612 int subcVertexCount = parentCell.getVertexCount(subcDim, subcOrd);
616 Kokkos::View<double**> refSubcellVertices(
"refSubcellVertices",subcVertexCount, cellDim);
617 Kokkos::View<double**> mappedParamVertices(
"mappedParamVertices",subcVertexCount, cellDim);
634 else if(subcDim == 2) {
637 if(subcVertexCount == 3){
645 else if(subcVertexCount == 4){
655 for(
int subcVertOrd = 0; subcVertOrd < subcVertexCount; subcVertOrd++){
656 for(
int dim = 0; dim < cellDim; dim++){
658 if(mappedParamVertices(subcVertOrd, dim) != refSubcellVertices(subcVertOrd, dim) ) {
661 << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n"
662 <<
" Cell Topology = " << parentCell.getName() <<
"\n"
663 <<
" Parametrization of subcell " << subcOrd <<
" which is "
664 << parentCell.getName(subcDim,subcOrd) <<
" failed for vertex " << subcVertOrd <<
":\n"
665 <<
" 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.