55 #include "Teuchos_oblackholestream.hpp"
56 #include "Teuchos_RCP.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
60 int main(
int argc,
char *argv[]) {
62 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
66 int iprint = argc - 1;
67 Teuchos::RCP<std::ostream> outStream;
68 Teuchos::oblackholestream bhs;
70 outStream = Teuchos::rcp(&std::cout,
false);
72 outStream = Teuchos::rcp(&bhs,
false);
75 Teuchos::oblackholestream oldFormatState;
76 oldFormatState.copyfmt(std::cout);
79 <<
"===============================================================================\n" \
81 <<
"| Unit Test (CubatureControlVolume) |\n" \
82 <<
"| (CubatureControlVolumeSide) |\n" \
83 <<
"| (CubatureControlVolumeBoundary) |\n" \
85 <<
"| 1) Correctness of cubature points and weights |\n"\
86 <<
"| 2) Comparison of sub-control volume weights and primary cell volume |\n"\
87 <<
"| 3) Control volume integration |\n"\
89 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
90 <<
"| Denis Ridzal (dridzal@sandia.gov) or |\n" \
91 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
93 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
94 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
96 <<
"===============================================================================\n"\
97 <<
"| TEST 1: correctness of cubature points and weights |\n"\
98 <<
"===============================================================================\n";
102 Teuchos::RCP<shards::CellTopology> cellTopo;
103 Teuchos::RCP<Intrepid::Cubature<double,Intrepid::FieldContainer<double> > > controlVolCub;
104 Teuchos::RCP<Intrepid::Cubature<double,Intrepid::FieldContainer<double> > > controlVolSideCub;
105 Teuchos::RCP<Intrepid::Cubature<double,Intrepid::FieldContainer<double> > > controlVolBCCub;
112 cellTopo = Teuchos::rcp(
new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<> >()));
116 int numPoints = controlVolCub->getNumPoints();
120 int numSidePoints = controlVolSideCub->getNumPoints();
125 int numBCPoints = controlVolBCCub->getNumPoints();
129 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0;
130 cellCoords(0,1,0) = 0.5; cellCoords(0,1,1) = 0.0;
131 cellCoords(0,2,0) = 0.5; cellCoords(0,2,1) = 1.0;
132 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 1.0;
133 cellCoords(1,0,0) = 0.5; cellCoords(1,0,1) = 0.0;
134 cellCoords(1,1,0) = 1.0; cellCoords(1,1,1) = 0.0;
135 cellCoords(1,2,0) = 1.0; cellCoords(1,2,1) = 1.0;
136 cellCoords(1,3,0) = 0.5; cellCoords(1,3,1) = 1.0;
139 double exactCubPoints[] = {
140 0.125, 0.25, 0.375, 0.25,
141 0.375, 0.75, 0.125, 0.75,
142 0.625, 0.25, 0.875, 0.25,
143 0.875, 0.75, 0.625, 0.75
147 double exactCubWeights[] = {
148 0.125, 0.125, 0.125, 0.125,
149 0.125, 0.125, 0.125, 0.125
153 double exactSideCubPoints[] = {
154 0.25, 0.25, 0.375, 0.5,
155 0.25, 0.75, 0.125, 0.5,
156 0.75, 0.25, 0.875, 0.5,
157 0.75, 0.75, 0.625, 0.5
161 double exactSideCubWeights[] = {
163 -0.5, 0.0, 0.0,-0.25,
169 double exactBCCubPoints[] = {
170 0.375, 1.0, 0.125, 1.0,
171 0.875, 1.0, 0.625, 1.0
175 double exactBCCubWeights[] = {
176 0.25, 0.25, 0.25, 0.25
180 if (numPoints != exactPoints) {
182 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
183 *outStream <<
"Number of cubature points: " << numPoints;
184 *outStream <<
" does not equal correct number: " << exactPoints <<
"\n";
187 int exactBCPoints = 2;
188 if (numBCPoints != exactBCPoints) {
190 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
191 *outStream <<
"Number of boundary cubature points: " << numBCPoints;
192 *outStream <<
" does not equal correct number: " << exactBCPoints <<
"\n";
196 int numCells = 2;
int spaceDim = 2;
199 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
205 controlVolSideCub->getCubature(sideCubPoints,sideCubWeights,cellCoords);
211 controlVolBCCub->getCubature(bcCubPoints,bcCubWeights,cellCoords);
217 for (
int i=0; i<numCells; i++) {
218 for (
int j=0; j<numPoints; j++) {
219 for (
int k=0; k<spaceDim; k++) {
222 if (std::abs(cubPoints(i,j,k) - exactCubPoints[countp]) > Intrepid::INTREPID_TOL) {
224 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
225 *outStream <<
" At multi-index { ";
226 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
227 *outStream <<
"} computed point: " << cubPoints(i,j,k)
228 <<
" but reference value: " << exactCubPoints[countp] <<
"\n";
232 if (std::abs(sideCubPoints(i,j,k) - exactSideCubPoints[countp]) > Intrepid::INTREPID_TOL) {
234 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
237 *outStream <<
" At multi-index { ";
238 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
239 *outStream <<
"} computed point: " << sideCubPoints(i,j,k)
240 <<
" but reference value: " << exactSideCubPoints[countp] <<
"\n";
244 if (std::abs(sideCubWeights(i,j,k) - exactSideCubWeights[countp]) > Intrepid::INTREPID_TOL) {
246 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
249 *outStream <<
" At multi-index { ";
250 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
251 *outStream <<
"} computed weight: " << sideCubWeights(i,j,k)
252 <<
" but reference value: " << exactSideCubWeights[countp] <<
"\n";
260 if (std::abs(cubWeights(i,j) - exactCubWeights[countw]) > Intrepid::INTREPID_TOL) {
262 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
263 *outStream <<
" At multi-index { ";
264 *outStream << i <<
" ";*outStream << j <<
" ";
265 *outStream <<
"} computed weight: " << cubWeights(i,j)
266 <<
" but reference value: " << exactCubWeights[countw] <<
"\n";
271 for (
int j=0; j<numBCPoints; j++) {
272 for (
int k=0; k<spaceDim; k++) {
275 if (std::abs(bcCubPoints(i,j,k) - exactBCCubPoints[countbcp]) > Intrepid::INTREPID_TOL) {
277 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
278 *outStream <<
" At multi-index { ";
279 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
280 *outStream <<
"} computed point: " << bcCubPoints(i,j,k)
281 <<
" but reference value: " << exactBCCubPoints[countbcp] <<
"\n";
287 if (std::abs(bcCubWeights(i,j) - exactBCCubWeights[countbcw]) > Intrepid::INTREPID_TOL) {
289 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
290 *outStream <<
" At multi-index { ";
291 *outStream << i <<
" ";*outStream << j <<
" ";
292 *outStream <<
"} computed weight: " << bcCubWeights(i,j)
293 <<
" but reference value: " << exactBCCubWeights[countbcw] <<
"\n";
299 catch (std::exception err) {
300 *outStream << err.what() <<
"\n";
308 cellTopo = Teuchos::rcp(
new shards::CellTopology(shards::getCellTopologyData< shards::Triangle<> >()));
312 int numPoints = controlVolCub->getNumPoints();
316 int numSidePoints = controlVolSideCub->getNumPoints();
321 int numBCPoints = controlVolBCCub->getNumPoints();
325 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0;
326 cellCoords(0,1,0) = 0.5; cellCoords(0,1,1) = 0.0;
327 cellCoords(0,2,0) = 0.0; cellCoords(0,2,1) = 0.5;
328 cellCoords(1,0,0) = 0.5; cellCoords(1,0,1) = 0.0;
329 cellCoords(1,1,0) = 0.5; cellCoords(1,1,1) = 0.5;
330 cellCoords(1,2,0) = 0.0; cellCoords(1,2,1) = 0.5;
333 double exactCubPoints[] = {
334 0.1041666666666667, 0.1041666666666667, 0.2916666666666667, 0.1041666666666667,
335 0.1041666666666667, 0.2916666666666667, 0.3958333333333333, 0.2083333333333333,
336 0.3958333333333333, 0.3958333333333333, 0.2083333333333333, 0.3958333333333333
340 double exactCubWeights[] = {
341 0.0416666666666667, 0.0416666666666667, 0.0416666666666667,
342 0.0416666666666667, 0.0416666666666667, 0.0416666666666667
346 double exactSideCubPoints[] = {
347 0.2083333333333333, 0.0833333333333333, 0.2083333333333333, 0.2083333333333333,
348 0.0833333333333333, 0.2083333333333333, 0.4166666666666667, 0.2916666666666667,
349 0.2916666666666667, 0.4166666666666667, 0.2916666666666667, 0.2916666666666667
353 double exactSideCubWeights[] = {
354 0.1666666666666667, 0.0833333333333333,-0.0833333333333333, 0.0833333333333333,
355 -0.0833333333333333,-0.1666666666666667, 0.0833333333333333, 0.1666666666666667,
356 -0.1666666666666667,-0.0833333333333333, 0.0833333333333333,-0.0833333333333333
360 double exactBCCubPoints[] = {
361 0.375, 0.125, 0.125, 0.375,
362 0.375, 0.5, 0.125, 0.5
366 double exactBCCubWeights[] = {
367 0.353553390593274, 0.353553390593274, 0.25, 0.25
372 if (numPoints != exactPoints) {
374 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
375 *outStream <<
"Number of cubature points: " << numPoints;
376 *outStream <<
" does not equal correct number: " << exactPoints <<
"\n";
378 if (numSidePoints != exactPoints) {
380 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
381 *outStream <<
"Number of side cubature points: " << numSidePoints;
382 *outStream <<
" does not equal correct number: " << exactPoints <<
"\n";
385 int exactBCPoints = 2;
386 if (numBCPoints != exactBCPoints) {
388 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
389 *outStream <<
"Number of boundary cubature points: " << numBCPoints;
390 *outStream <<
" does not equal correct number: " << exactBCPoints <<
"\n";
394 int numCells = 2;
int spaceDim = 2;
397 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
403 controlVolSideCub->getCubature(sideCubPoints,sideCubWeights,cellCoords);
409 controlVolBCCub->getCubature(bcCubPoints,bcCubWeights,cellCoords);
415 for (
int i=0; i<numCells; i++) {
416 for (
int j=0; j<numPoints; j++) {
417 for (
int k=0; k<spaceDim; k++) {
420 if (std::abs(cubPoints(i,j,k) - exactCubPoints[countp]) > Intrepid::INTREPID_TOL) {
422 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
423 *outStream <<
" At multi-index { ";
424 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
425 *outStream <<
"} computed point: " << cubPoints(i,j,k)
426 <<
" but reference value: " << exactCubPoints[countp] <<
"\n";
430 if (std::abs(sideCubPoints(i,j,k) - exactSideCubPoints[countp]) > Intrepid::INTREPID_TOL) {
432 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
435 *outStream <<
" At multi-index { ";
436 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
437 *outStream <<
"} computed point: " << sideCubPoints(i,j,k)
438 <<
" but reference value: " << exactSideCubPoints[countp] <<
"\n";
442 if (std::abs(sideCubWeights(i,j,k) - exactSideCubWeights[countp]) > Intrepid::INTREPID_TOL) {
444 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
447 *outStream <<
" At multi-index { ";
448 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
449 *outStream <<
"} computed weight: " << sideCubWeights(i,j,k)
450 <<
" but reference value: " << exactSideCubWeights[countp] <<
"\n";
458 if (std::abs(cubWeights(i,j) - exactCubWeights[countw]) > Intrepid::INTREPID_TOL) {
460 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
461 *outStream <<
" At multi-index { ";
462 *outStream << i <<
" ";*outStream << j <<
" ";
463 *outStream <<
"} computed weight: " << cubWeights(i,j)
464 <<
" but reference value: " << exactCubWeights[countw] <<
"\n";
469 for (
int j=0; j<numBCPoints; j++) {
470 for (
int k=0; k<spaceDim; k++) {
473 if (std::abs(bcCubPoints(i,j,k) - exactBCCubPoints[countbcp]) > Intrepid::INTREPID_TOL) {
475 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
476 *outStream <<
" At multi-index { ";
477 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
478 *outStream <<
"} computed point: " << bcCubPoints(i,j,k)
479 <<
" but reference value: " << exactBCCubPoints[countbcp] <<
"\n";
485 if (std::abs(bcCubWeights(i,j) - exactBCCubWeights[countbcw]) > Intrepid::INTREPID_TOL) {
487 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
488 *outStream <<
" At multi-index { ";
489 *outStream << i <<
" ";*outStream << j <<
" ";
490 *outStream <<
"} computed weight: " << bcCubWeights(i,j)
491 <<
" but reference value: " << bcCubWeights(i,j) - exactBCCubWeights[countbcw] <<
"\n";
498 catch (std::exception err) {
499 *outStream << err.what() <<
"\n";
507 cellTopo = Teuchos::rcp(
new shards::CellTopology(shards::getCellTopologyData< shards::Tetrahedron<> >()));
511 int numPoints = controlVolCub->getNumPoints();
515 int numSidePoints = controlVolSideCub->getNumPoints();
520 int numBCPoints = controlVolBCCub->getNumPoints();
524 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0; cellCoords(0,0,2) = 0.0;
525 cellCoords(0,1,0) = 1.0; cellCoords(0,1,1) = 0.0; cellCoords(0,1,2) = 0.0;
526 cellCoords(0,2,0) = 0.0; cellCoords(0,2,1) = 1.0; cellCoords(0,2,2) = 0.0;
527 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 0.0; cellCoords(0,3,2) = 1.0;
530 double exactCubPoints[] = {
531 0.17708333333333333, 0.17708333333333333, 0.17708333333333333,
532 0.46875000000000000, 0.17708333333333333, 0.17708333333333333,
533 0.17708333333333333, 0.46875000000000000, 0.17708333333333333,
534 0.17708333333333333, 0.17708333333333333, 0.46875000000000000
538 double exactCubWeights[] = {
539 0.0416666666666667, 0.0416666666666667,
540 0.0416666666666667, 0.0416666666666667
544 double exactSideCubPoints[] = {
545 0.3541666666666667, 0.1458333333333333, 0.1458333333333333,
546 0.3541666666666667, 0.3541666666666667, 0.1458333333333333,
547 0.1458333333333333, 0.3541666666666667, 0.1458333333333333,
548 0.1458333333333333, 0.1458333333333333, 0.3541666666666667,
549 0.3541666666666667, 0.1458333333333333, 0.3541666666666667,
550 0.1458333333333333, 0.3541666666666667, 0.3541666666666667
554 double exactSideCubWeights[] = {
555 0.0833333333333333, 0.0416666666666667, 0.041666666666667,
556 -0.0416666666666667, 0.0416666666666667, 0.000000000000000,
557 -0.0416666666666667,-0.0833333333333333,-0.041666666666667,
558 0.0416666666666667, 0.0416666666666667, 0.083333333333333,
559 -0.0416666666666667, 0.0000000000000000, 0.041666666666667,
560 0.0000000000000000,-0.0416666666666667, 0.041666666666667
564 double exactBCCubPoints[] = {
565 0.208333333333333, 0.00, 0.208333333333333,
566 0.583333333333333, 0.00, 0.208333333333333,
567 0.208333333333333, 0.00, 0.583333333333333,
571 double exactBCCubWeights[] = {
572 0.166666666666667, 0.166666666666667, 0.166666666666667
577 if (numPoints != exactPoints) {
579 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
580 *outStream <<
"Number of cubature points: " << numPoints;
581 *outStream <<
" does not equal correct number: " << exactPoints <<
"\n";
585 if (numSidePoints != exactPoints) {
587 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
588 *outStream <<
"Number of side cubature points: " << numSidePoints;
589 *outStream <<
" does not equal correct number: " << exactPoints <<
"\n";
592 int exactBCPoints = 3;
593 if (numBCPoints != exactBCPoints) {
595 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
596 *outStream <<
"Number of boundary cubature points: " << numBCPoints;
597 *outStream <<
" does not equal correct number: " << exactBCPoints <<
"\n";
601 int numCells = 1;
int spaceDim = 3;
604 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
610 controlVolSideCub->getCubature(sideCubPoints,sideCubWeights,cellCoords);
616 controlVolBCCub->getCubature(bcCubPoints,bcCubWeights,cellCoords);
622 for (
int i=0; i<numCells; i++) {
623 for (
int j=0; j<numPoints; j++) {
624 for (
int k=0; k<spaceDim; k++) {
627 if (std::abs(cubPoints(i,j,k) - exactCubPoints[countp]) > Intrepid::INTREPID_TOL) {
629 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
630 *outStream <<
" At multi-index { ";
631 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
632 *outStream <<
"} computed point: " << cubPoints(i,j,k)
633 <<
" but reference value: " << exactCubPoints[countp] <<
"\n";
637 if (std::abs(sideCubPoints(i,j,k) - exactSideCubPoints[countp]) > Intrepid::INTREPID_TOL) {
639 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
642 *outStream <<
" At multi-index { ";
643 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
644 *outStream <<
"} computed point: " << sideCubPoints(i,j,k)
645 <<
" but reference value: " << exactSideCubPoints[countp] <<
"\n";
649 if (std::abs(sideCubWeights(i,j,k) - exactSideCubWeights[countp]) > Intrepid::INTREPID_TOL) {
651 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
654 *outStream <<
" At multi-index { ";
655 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
656 *outStream <<
"} computed weight: " << sideCubWeights(i,j,k)
657 <<
" but reference value: " << exactSideCubWeights[countp] <<
"\n";
665 if (std::abs(cubWeights(i,j) - exactCubWeights[countw]) > Intrepid::INTREPID_TOL) {
667 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
668 *outStream <<
" At multi-index { ";
669 *outStream << i <<
" ";*outStream << j <<
" ";
670 *outStream <<
"} computed weight: " << cubWeights(i,j)
671 <<
" but reference value: " << exactCubWeights[countw] <<
"\n";
672 *outStream << cubWeights(i,j)-exactCubWeights[countp] <<
"\n";
677 for (
int j=0; j<numBCPoints; j++) {
678 for (
int k=0; k<spaceDim; k++) {
681 if (std::abs(bcCubPoints(i,j,k) - exactBCCubPoints[countbcp]) > Intrepid::INTREPID_TOL) {
683 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
684 *outStream <<
" At multi-index { ";
685 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
686 *outStream <<
"} computed point: " << bcCubPoints(i,j,k)
687 <<
" but reference value: " << exactBCCubPoints[countbcp] <<
"\n";
693 if (std::abs(bcCubWeights(i,j) - exactBCCubWeights[countbcw]) > Intrepid::INTREPID_TOL) {
695 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
696 *outStream <<
" At multi-index { ";
697 *outStream << i <<
" ";*outStream << j <<
" ";
698 *outStream <<
"} computed weight: " << bcCubWeights(i,j)
699 <<
" but reference value: " << exactBCCubWeights[countbcw] <<
"\n";
707 catch (std::exception err) {
708 *outStream << err.what() <<
"\n";
716 cellTopo = Teuchos::rcp(
new shards::CellTopology(shards::getCellTopologyData< shards::Hexahedron<> >()));
720 int numPoints = controlVolCub->getNumPoints();
724 int numSidePoints = controlVolSideCub->getNumPoints();
729 int numBCPoints = controlVolBCCub->getNumPoints();
733 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0; cellCoords(0,0,2) = 0.0;
734 cellCoords(0,1,0) = 2.0; cellCoords(0,1,1) = 0.0; cellCoords(0,1,2) = 0.0;
735 cellCoords(0,2,0) = 2.0; cellCoords(0,2,1) = 1.5; cellCoords(0,2,2) = 0.0;
736 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 1.5; cellCoords(0,3,2) = 0.0;
737 cellCoords(0,4,0) = 0.0; cellCoords(0,4,1) = 0.0; cellCoords(0,4,2) = 1.0;
738 cellCoords(0,5,0) = 2.0; cellCoords(0,5,1) = 0.0; cellCoords(0,5,2) = 1.0;
739 cellCoords(0,6,0) = 2.0; cellCoords(0,6,1) = 1.5; cellCoords(0,6,2) = 1.0;
740 cellCoords(0,7,0) = 0.0; cellCoords(0,7,1) = 1.5; cellCoords(0,7,2) = 1.0;
743 double exactCubPoints[] = {
744 0.5, 0.375, 0.25, 1.5, 0.375, 0.25,
745 1.5, 1.125, 0.25, 0.5, 1.125, 0.25,
746 0.5, 0.375, 0.75, 1.5, 0.375, 0.75,
747 1.5, 1.125, 0.75, 0.5, 1.125, 0.75
751 double exactCubWeights[] = {
752 0.375, 0.375, 0.375, 0.375,
753 0.375, 0.375, 0.375, 0.375
757 double exactSideCubPoints[] = {
758 1.0, 0.375, 0.25, 1.5, 0.750, 0.25,
759 1.0, 1.125, 0.25, 0.5, 0.750, 0.25,
760 1.0, 0.375, 0.75, 1.5, 0.750, 0.75,
761 1.0, 1.125, 0.75, 0.5, 0.750, 0.75,
762 0.5, 0.375, 0.50, 1.5, 0.375, 0.50,
763 1.5, 1.125, 0.50, 0.5, 1.125, 0.50
767 double exactSideCubWeights[] = {
768 0.375, 0.00, 0.00, 0.00, 0.50, 0.00,
769 -0.375, 0.00, 0.00, 0.00,-0.50, 0.00,
770 0.375, 0.00, 0.00, 0.00, 0.50, 0.00,
771 -0.375, 0.00, 0.00, 0.00,-0.50, 0.00,
772 0.000, 0.00, 0.75, 0.00, 0.00, 0.75,
773 0.000, 0.00, 0.75, 0.00, 0.00, 0.75,
778 if (numPoints != exactPoints) {
780 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
781 *outStream <<
"Number of cubature points: " << numPoints;
782 *outStream <<
" does not equal correct number: " << exactPoints <<
"\n";
786 if (numSidePoints != exactPoints) {
788 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
789 *outStream <<
"Number of side cubature points: " << numSidePoints;
790 *outStream <<
" does not equal correct number: " << exactPoints <<
"\n";
794 double exactBCCubPoints[] = {
802 double exactBCCubWeights[] = {
807 int numCells = 1;
int spaceDim = 3;
810 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
816 controlVolSideCub->getCubature(sideCubPoints,sideCubWeights,cellCoords);
822 controlVolBCCub->getCubature(bcCubPoints,bcCubWeights,cellCoords);
828 for (
int i=0; i<numCells; i++) {
829 for (
int j=0; j<numPoints; j++) {
830 for (
int k=0; k<spaceDim; k++) {
833 if (std::abs(cubPoints(i,j,k) - exactCubPoints[countp]) > Intrepid::INTREPID_TOL) {
835 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
836 *outStream <<
" At multi-index { ";
837 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
838 *outStream <<
"} computed point: " << cubPoints(i,j,k)
839 <<
" but reference value: " << exactCubPoints[countp] <<
"\n";
843 if (std::abs(sideCubPoints(i,j,k) - exactSideCubPoints[countp]) > Intrepid::INTREPID_TOL) {
845 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
848 *outStream <<
" At multi-index { ";
849 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
850 *outStream <<
"} computed point: " << sideCubPoints(i,j,k)
851 <<
" but reference value: " << exactSideCubPoints[countp] <<
"\n";
855 if (std::abs(sideCubWeights(i,j,k) - exactSideCubWeights[countp]) > Intrepid::INTREPID_TOL) {
857 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
860 *outStream <<
" At multi-index { ";
861 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
862 *outStream <<
"} computed weight: " << sideCubWeights(i,j,k)
863 <<
" but reference value: " << exactSideCubWeights[countp] <<
"\n";
871 if (std::abs(cubWeights(i,j) - exactCubWeights[countw]) > Intrepid::INTREPID_TOL) {
873 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
874 *outStream <<
" At multi-index { ";
875 *outStream << i <<
" ";*outStream << j <<
" ";
876 *outStream <<
"} computed weight: " << cubWeights(i,j)
877 <<
" but reference value: " << exactCubWeights[countw] <<
"\n";
878 *outStream << cubWeights(i,j)-exactCubWeights[countp] <<
"\n";
883 for (
int j=0; j<numBCPoints; j++) {
884 for (
int k=0; k<spaceDim; k++) {
887 if (std::abs(bcCubPoints(i,j,k) - exactBCCubPoints[countbcp]) > Intrepid::INTREPID_TOL) {
889 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
890 *outStream <<
" At multi-index { ";
891 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
892 *outStream <<
"} computed point: " << bcCubPoints(i,j,k)
893 <<
" but reference value: " << exactBCCubPoints[countbcp] <<
"\n";
899 if (std::abs(bcCubWeights(i,j) - exactBCCubWeights[countbcw]) > Intrepid::INTREPID_TOL) {
901 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
902 *outStream <<
" At multi-index { ";
903 *outStream << i <<
" ";*outStream << j <<
" ";
904 *outStream <<
"} computed weight: " << bcCubWeights(i,j)
905 <<
" but reference value: " << exactBCCubWeights[countbcw] <<
"\n";
912 catch (std::exception err) {
913 *outStream << err.what() <<
"\n";
918 <<
"===============================================================================\n"\
919 <<
"| TEST 2: comparison of sub-control volume weights and primary cell volume |\n"\
920 <<
"===============================================================================\n";
926 cellTopo = Teuchos::rcp(
new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<> >()));
930 int numPoints = controlVolCub->getNumPoints();
934 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0;
935 cellCoords(0,1,0) = 2.4; cellCoords(0,1,1) = 0.0;
936 cellCoords(0,2,0) = 2.4; cellCoords(0,2,1) = 3.1;
937 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 3.1;
939 double exact_area = 2.4*3.1;
942 int numCells = 1;
int spaceDim = 2;
945 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
948 double total_area = 0.0;
949 for (
int i=0; i<numPoints; i++) {
950 total_area += cubWeights(0,i);
952 if (std::abs(total_area - exact_area) > Intrepid::INTREPID_TOL) {
954 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
955 *outStream <<
" Sum of sub-control volume areas: ";
956 *outStream << total_area;
957 *outStream <<
" does not equal quad primary cell area: " << exact_area <<
"\n";
961 catch (std::exception err) {
962 *outStream << err.what() <<
"\n";
970 cellTopo = Teuchos::rcp(
new shards::CellTopology(shards::getCellTopologyData< shards::Triangle<> >()));
974 int numPoints = controlVolCub->getNumPoints();
978 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0;
979 cellCoords(0,1,0) = 3.6; cellCoords(0,1,1) = 0.0;
980 cellCoords(0,2,0) = 0.0; cellCoords(0,2,1) = 2.8;
982 double exact_area = 0.5*3.6*2.8;
985 int numCells = 1;
int spaceDim = 2;
988 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
991 double total_area = 0.0;
992 for (
int i=0; i<numPoints; i++) {
993 total_area += cubWeights(0,i);
995 if (std::abs(total_area - exact_area) > Intrepid::INTREPID_TOL) {
997 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
998 *outStream <<
" Sum of sub-control volume areas: ";
999 *outStream << total_area;
1000 *outStream <<
" does not equal triangle primary cell area: " << exact_area <<
"\n";
1004 catch (std::exception err) {
1005 *outStream << err.what() <<
"\n";
1013 cellTopo = Teuchos::rcp(
new shards::CellTopology(shards::getCellTopologyData< shards::Hexahedron<> >()));
1017 int numPoints = controlVolCub->getNumPoints();
1021 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0; cellCoords(0,0,2) = 0.0;
1022 cellCoords(0,1,0) = 2.4; cellCoords(0,1,1) = 0.0; cellCoords(0,1,2) = 0.0;
1023 cellCoords(0,2,0) = 2.4; cellCoords(0,2,1) = 3.1; cellCoords(0,2,2) = 0.0;
1024 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 3.1; cellCoords(0,3,2) = 0.0;
1025 cellCoords(0,4,0) = 0.0; cellCoords(0,4,1) = 0.0; cellCoords(0,4,2) = 1.7;
1026 cellCoords(0,5,0) = 2.4; cellCoords(0,5,1) = 0.0; cellCoords(0,5,2) = 1.7;
1027 cellCoords(0,6,0) = 2.4; cellCoords(0,6,1) = 3.1; cellCoords(0,6,2) = 1.7;
1028 cellCoords(0,7,0) = 0.0; cellCoords(0,7,1) = 3.1; cellCoords(0,7,2) = 1.7;
1030 double exact_vol = 2.4*3.1*1.7;
1033 int numCells = 1;
int spaceDim = 3;
1036 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
1039 double total_vol = 0.0;
1040 for (
int i=0; i<numPoints; i++) {
1041 total_vol += cubWeights(0,i);
1043 if (std::abs(total_vol - exact_vol) > Intrepid::INTREPID_TOL) {
1045 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
1046 *outStream <<
" Sum of sub-control volumes: ";
1047 *outStream << total_vol;
1048 *outStream <<
" does not equal hexahedron primary cell volume: " << exact_vol <<
"\n";
1052 catch (std::exception err) {
1053 *outStream << err.what() <<
"\n";
1061 cellTopo = Teuchos::rcp(
new shards::CellTopology(shards::getCellTopologyData< shards::Tetrahedron<> >()));
1065 int numPoints = controlVolCub->getNumPoints();
1069 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0; cellCoords(0,0,2) = 0.0;
1070 cellCoords(0,1,0) = 3.6; cellCoords(0,1,1) = 0.0; cellCoords(0,1,2) = 0.0;
1071 cellCoords(0,2,0) = 0.0; cellCoords(0,2,1) = 2.8; cellCoords(0,2,2) = 0.0;
1072 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 2.8; cellCoords(0,3,2) = 1.7;
1074 double exact_vol = (0.5*3.6*2.8)*1.7/3.0;
1077 int numCells = 1;
int spaceDim = 3;
1080 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
1083 double total_vol = 0.0;
1084 for (
int i=0; i<numPoints; i++) {
1085 total_vol += cubWeights(0,i);
1087 if (std::abs(total_vol - exact_vol) > Intrepid::INTREPID_TOL) {
1089 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
1090 *outStream <<
" Sum of sub-control volumes: ";
1091 *outStream << total_vol;
1092 *outStream <<
" does not equal tetrahedron primary cell volume: " << exact_vol <<
"\n";
1096 catch (std::exception err) {
1097 *outStream << err.what() <<
"\n";
1102 <<
"===============================================================================\n"\
1103 <<
"| TEST 3: control volume integration |\n"\
1104 <<
"===============================================================================\n";
1110 cellTopo = Teuchos::rcp(
new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<> >()));
1114 int numPoints = controlVolCub->getNumPoints();
1118 int numSidePoints = controlVolSideCub->getNumPoints();
1122 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0;
1123 cellCoords(0,1,0) = 0.5; cellCoords(0,1,1) = 0.0;
1124 cellCoords(0,2,0) = 0.41; cellCoords(0,2,1) = 0.58;
1125 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 0.5;
1126 cellCoords(1,0,0) = 0.5; cellCoords(1,0,1) = 0.0;
1127 cellCoords(1,1,0) = 1.0; cellCoords(1,1,1) = 0.0;
1128 cellCoords(1,2,0) = 1.0; cellCoords(1,2,1) = 0.5;
1129 cellCoords(1,3,0) = 0.41; cellCoords(1,3,1) = 0.58;
1130 cellCoords(2,0,0) = 0.0; cellCoords(2,0,1) = 0.5;
1131 cellCoords(2,1,0) = 0.41; cellCoords(2,1,1) = 0.58;
1132 cellCoords(2,2,0) = 0.5; cellCoords(2,2,1) = 1.0;
1133 cellCoords(2,3,0) = 0.0; cellCoords(2,3,1) = 1.0;
1134 cellCoords(3,0,0) = 0.41; cellCoords(3,0,1) = 0.58;
1135 cellCoords(3,1,0) = 1.0; cellCoords(3,1,1) = 0.5;
1136 cellCoords(3,2,0) = 1.0; cellCoords(3,2,1) = 1.0;
1137 cellCoords(3,3,0) = 0.5; cellCoords(3,3,1) = 1.0;
1140 int numCells = 4;
int spaceDim = 2;
1143 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
1147 controlVolSideCub->getCubature(sideCubPoints,sideCubWeights,cellCoords);
1153 double a = 2.1;
double b = 1.4;
1155 for (
int i = 0; i < numCells; i++) {
1156 for (
int j = 0; j < numSidePoints; j++) {
1157 F(i,j,0) = a*sideCubPoints(i,j,0);
1158 F(i,j,1) = b*sideCubPoints(i,j,1);
1163 double surfaceInt = 0.0;
1166 surfaceInt += - F(0,1,0)*sideCubWeights(0,1,0) - F(0,1,1)*sideCubWeights(0,1,1)
1167 + F(0,2,0)*sideCubWeights(0,2,0) + F(0,2,1)*sideCubWeights(0,2,1);
1170 surfaceInt += - F(1,2,0)*sideCubWeights(1,2,0) - F(1,2,1)*sideCubWeights(1,2,1)
1171 + F(1,3,0)*sideCubWeights(1,3,0) + F(1,3,1)*sideCubWeights(1,3,1);
1174 surfaceInt += - F(2,0,0)*sideCubWeights(2,0,0) - F(2,0,1)*sideCubWeights(2,0,1)
1175 + F(2,1,0)*sideCubWeights(2,1,0) + F(2,1,1)*sideCubWeights(2,1,1);
1178 surfaceInt += - F(3,3,0)*sideCubWeights(3,3,0) - F(3,3,1)*sideCubWeights(3,3,1)
1179 + F(3,0,0)*sideCubWeights(3,0,0) + F(3,0,1)*sideCubWeights(3,0,1);
1182 double volumeInt = 0.0;
1185 volumeInt += (a+b)*cubWeights(0,2);
1188 volumeInt += (a+b)*cubWeights(1,3);
1191 volumeInt += (a+b)*cubWeights(2,1);
1194 volumeInt += (a+b)*cubWeights(3,0);
1196 if (std::abs(surfaceInt - volumeInt) > Intrepid::INTREPID_TOL) {
1198 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
1199 *outStream <<
" Integral of (F cdot n) over surface : ";
1200 *outStream << surfaceInt;
1201 *outStream <<
" does not equal integral of div F over volume: " << volumeInt <<
"\n";
1205 catch (std::exception err) {
1206 *outStream << err.what() <<
"\n";
1214 cellTopo = Teuchos::rcp(
new shards::CellTopology(shards::getCellTopologyData< shards::Triangle<> >()));
1218 int numPoints = controlVolCub->getNumPoints();
1222 int numSidePoints = controlVolSideCub->getNumPoints();
1226 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0;
1227 cellCoords(0,1,0) = 1.0; cellCoords(0,1,1) = 0.0;
1228 cellCoords(0,2,0) = 0.41; cellCoords(0,2,1) = 0.58;
1229 cellCoords(1,0,0) = 1.0; cellCoords(1,0,1) = 0.0;
1230 cellCoords(1,1,0) = 1.0; cellCoords(1,1,1) = 1.0;
1231 cellCoords(1,2,0) = 0.41; cellCoords(1,2,1) = 0.58;
1232 cellCoords(2,0,0) = 1.0; cellCoords(2,0,1) = 1.0;
1233 cellCoords(2,1,0) = 0.0; cellCoords(2,1,1) = 1.0;
1234 cellCoords(2,2,0) = 0.41; cellCoords(2,2,1) = 0.58;
1235 cellCoords(3,0,0) = 0.0; cellCoords(3,0,1) = 1.0;
1236 cellCoords(3,1,0) = 0.0; cellCoords(3,1,1) = 0.0;
1237 cellCoords(3,2,0) = 0.41; cellCoords(3,2,1) = 0.58;
1240 int numCells = 4;
int spaceDim = 2;
1243 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
1247 controlVolSideCub->getCubature(sideCubPoints,sideCubWeights,cellCoords);
1253 double a = 2.1;
double b = 1.4;
1255 for (
int i = 0; i < numCells; i++) {
1256 for (
int j = 0; j < numSidePoints; j++) {
1257 F(i,j,0) = a*sideCubPoints(i,j,0);
1258 F(i,j,1) = b*sideCubPoints(i,j,1);
1263 double surfaceInt = 0.0;
1266 surfaceInt += - F(0,1,0)*sideCubWeights(0,1,0) - F(0,1,1)*sideCubWeights(0,1,1)
1267 + F(0,2,0)*sideCubWeights(0,2,0) + F(0,2,1)*sideCubWeights(0,2,1);
1270 surfaceInt += - F(1,1,0)*sideCubWeights(1,1,0) - F(1,1,1)*sideCubWeights(1,1,1)
1271 + F(1,2,0)*sideCubWeights(1,2,0) + F(1,2,1)*sideCubWeights(1,2,1);
1274 surfaceInt += - F(2,1,0)*sideCubWeights(2,1,0) - F(2,1,1)*sideCubWeights(2,1,1)
1275 + F(2,2,0)*sideCubWeights(2,2,0) + F(2,2,1)*sideCubWeights(2,2,1);
1278 surfaceInt += - F(3,1,0)*sideCubWeights(3,1,0) - F(3,1,1)*sideCubWeights(3,1,1)
1279 + F(3,2,0)*sideCubWeights(3,2,0) + F(3,2,1)*sideCubWeights(3,2,1);
1282 double volumeInt = 0.0;
1285 volumeInt += (a+b)*cubWeights(0,2);
1288 volumeInt += (a+b)*cubWeights(1,2);
1291 volumeInt += (a+b)*cubWeights(2,2);
1294 volumeInt += (a+b)*cubWeights(3,2);
1296 if (std::abs(surfaceInt - volumeInt) > Intrepid::INTREPID_TOL) {
1298 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
1299 *outStream <<
" Integral of (F cdot n) over surface : ";
1300 *outStream << surfaceInt;
1301 *outStream <<
" does not equal integral of div F over volume: " << volumeInt <<
"\n";
1305 catch (std::exception err) {
1306 *outStream << err.what() <<
"\n";
1312 std::cout <<
"End Result: TEST FAILED\n";
1314 std::cout <<
"End Result: TEST PASSED\n";
1317 std::cout.copyfmt(oldFormatState);
Defines cubature (integration) rules over Neumann boundaries for control volume method.
Header file for the Intrepid::CubatureControlVolume class.
Header file for the Intrepid::CubatureControlVolumeSide class.
Header file for the Intrepid::CubatureControlVolumeBoundary class.
Defines cubature (integration) rules over control volumes.
Defines cubature (integration) rules over control volumes.