50#include "Teuchos_oblackholestream.hpp"
51#include "Teuchos_RCP.hpp"
52#include "Teuchos_GlobalMPISession.hpp"
55using namespace Intrepid;
57#define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
63 catch (const std::logic_error & err) { \
65 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
66 *outStream << err.what() << '\n'; \
67 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
71int main(
int argc,
char *argv[]) {
73 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
77 int iprint = argc - 1;
78 Teuchos::RCP<std::ostream> outStream;
79 Teuchos::oblackholestream bhs;
81 outStream = Teuchos::rcp(&std::cout,
false);
83 outStream = Teuchos::rcp(&bhs,
false);
86 Teuchos::oblackholestream oldFormatState;
87 oldFormatState.copyfmt(std::cout);
90 <<
"===============================================================================\n" \
92 <<
"| Unit Test (Basis_HGRAD_WEDGE_I2_FEM) |\n" \
94 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 <<
"| 2) Basis values for VALUE, GRAD, and Dk operators |\n" \
97 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
98 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
99 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
101 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
102 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
104 <<
"===============================================================================\n"\
105 <<
"| TEST 1: Basis creation, exception testing |\n"\
106 <<
"===============================================================================\n";
109 Basis_HGRAD_WEDGE_I2_FEM<double, FieldContainer<double> > wedgeBasis;
114 int throwCounter = 0;
117 FieldContainer<double> wedgeNodes(15, 3);
118 wedgeNodes(0,0) = 0.0; wedgeNodes(0,1) = 0.0; wedgeNodes(0,2) = -1.0;
119 wedgeNodes(1,0) = 1.0; wedgeNodes(1,1) = 0.0; wedgeNodes(1,2) = -1.0;
120 wedgeNodes(2,0) = 0.0; wedgeNodes(2,1) = 1.0; wedgeNodes(2,2) = -1.0;
121 wedgeNodes(3,0) = 0.0; wedgeNodes(3,1) = 0.0; wedgeNodes(3,2) = 1.0;
122 wedgeNodes(4,0) = 1.0; wedgeNodes(4,1) = 0.0; wedgeNodes(4,2) = 1.0;
123 wedgeNodes(5,0) = 0.0; wedgeNodes(5,1) = 1.0; wedgeNodes(5,2) = 1.0;
125 wedgeNodes(6,0) = 0.5; wedgeNodes(6,1) = 0.0; wedgeNodes(6,2) = -1.0;
126 wedgeNodes(7,0) = 0.5; wedgeNodes(7,1) = 0.5; wedgeNodes(7,2) = -1.0;
127 wedgeNodes(8,0) = 0.0; wedgeNodes(8,1) = 0.5; wedgeNodes(8,2) = -1.0;
128 wedgeNodes(9,0) = 0.0; wedgeNodes(9,1) = 0.0; wedgeNodes(9,2) = 0.0;
129 wedgeNodes(10,0)= 1.0; wedgeNodes(10,1)= 0.0; wedgeNodes(10,2)= 0.0;
130 wedgeNodes(11,0)= 0.0; wedgeNodes(11,1)= 1.0; wedgeNodes(11,2)= 0.0;
131 wedgeNodes(12,0)= 0.5; wedgeNodes(12,1)= 0.0; wedgeNodes(12,2)= 1.0;
132 wedgeNodes(13,0)= 0.5; wedgeNodes(13,1)= 0.5; wedgeNodes(13,2)= 1.0;
133 wedgeNodes(14,0)= 0.0; wedgeNodes(14,1)= 0.5; wedgeNodes(14,2)= 1.0;
137 FieldContainer<double> vals;
142 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 );
143 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
147 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) );
148 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
153 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException );
155 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException );
157 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,9,0), throwCounter, nException );
159 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(15), throwCounter, nException );
161 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException );
163#ifdef HAVE_INTREPID_DEBUG
166 FieldContainer<double> badPoints1(4, 5, 3);
167 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
170 FieldContainer<double> badPoints2(4, wedgeBasis.getBaseCellTopology().getDimension() + 1);
171 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
174 FieldContainer<double> badVals1(4, 3, 1);
175 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
178 FieldContainer<double> badVals2(4, 3);
179 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
182 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D1), throwCounter, nException );
185 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D2), throwCounter, nException );
188 FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0));
189 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
192 FieldContainer<double> badVals4(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1);
193 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
196 FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), wedgeBasis.getBaseCellTopology().getDimension() - 1);
197 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
200 FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 40);
201 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D2), throwCounter, nException );
204 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D3), throwCounter, nException );
208 catch (
const std::logic_error & err) {
209 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
210 *outStream << err.what() <<
'\n';
211 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
216 if (throwCounter != nException) {
218 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
223 <<
"===============================================================================\n"\
224 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
225 <<
"===============================================================================\n";
228 std::vector<std::vector<int> > allTags = wedgeBasis.getAllDofTags();
231 for (
unsigned i = 0; i < allTags.size(); i++) {
232 int bfOrd = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
234 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
235 if( !( (myTag[0] == allTags[i][0]) &&
236 (myTag[1] == allTags[i][1]) &&
237 (myTag[2] == allTags[i][2]) &&
238 (myTag[3] == allTags[i][3]) ) ) {
240 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
241 *outStream <<
" getDofOrdinal( {"
242 << allTags[i][0] <<
", "
243 << allTags[i][1] <<
", "
244 << allTags[i][2] <<
", "
245 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
246 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
250 << myTag[3] <<
"}\n";
255 for(
int bfOrd = 0; bfOrd < wedgeBasis.getCardinality(); bfOrd++) {
256 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
257 int myBfOrd = wedgeBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
258 if( bfOrd != myBfOrd) {
260 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
261 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
265 << myTag[3] <<
"} but getDofOrdinal({"
269 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
273 catch (
const std::logic_error & err){
274 *outStream << err.what() <<
"\n\n";
280 <<
"===============================================================================\n"\
281 <<
"| TEST 3: correctness of basis function values |\n"\
282 <<
"===============================================================================\n";
284 outStream -> precision(20);
287 double basisValues[] = {
288 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
289 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
290 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
291 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
292 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
293 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
294 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
295 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
296 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
297 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
298 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,\
299 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,\
300 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,\
301 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,\
302 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0};
305 std::string fileName;
306 std::ifstream dataFile;
309 std::vector<double> basisGrads;
311 fileName =
"./testdata/WEDGE_I2_GradVals.dat";
312 dataFile.open(fileName.c_str());
313 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
314 ">>> ERROR (HGRAD_WEDGE_I2/test01): could not open GRAD values data file, test aborted.");
315 while (!dataFile.eof() ){
318 std::getline(dataFile, line);
319 stringstream data_line(line);
320 while(data_line >> temp){
321 basisGrads.push_back(temp);
332 std::vector<double> basisD2;
334 fileName =
"./testdata/WEDGE_I2_D2Vals.dat";
335 dataFile.open(fileName.c_str());
336 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
337 ">>> ERROR (HGRAD_WEDGE_I2/test01): could not open D2 values data file, test aborted.");
339 while (!dataFile.eof() ){
342 std::getline(dataFile, line);
343 stringstream data_line(line);
344 while(data_line >> temp){
345 basisD2.push_back(temp);
353 std::vector<double> basisD3;
355 fileName =
"./testdata/WEDGE_I2_D3Vals.dat";
356 dataFile.open(fileName.c_str());
357 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
358 ">>> ERROR (HGRAD_WEDGE_I2/test01): could not open D3 values data file, test aborted.");
360 while (!dataFile.eof() ){
363 std::getline(dataFile, line);
364 stringstream data_line(line);
365 while(data_line >> temp){
366 basisD3.push_back(temp);
375 int numFields = wedgeBasis.getCardinality();
376 int numPoints = wedgeNodes.dimension(0);
377 int spaceDim = wedgeBasis.getBaseCellTopology().getDimension();
380 FieldContainer<double> vals;
383 vals.resize(numFields, numPoints);
384 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE);
385 for (
int i = 0; i < numFields; i++) {
386 for (
int j = 0; j < numPoints; j++) {
387 int l = i + j * numFields;
388 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
390 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
393 *outStream <<
" At multi-index { ";
394 *outStream << i <<
" ";*outStream << j <<
" ";
395 *outStream <<
"} computed value: " << vals(i,j)
396 <<
" but reference value: " << basisValues[l] <<
"\n";
404 vals.resize(numFields, numPoints, spaceDim);
405 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD);
406 for (
int i = 0; i < numFields; i++) {
407 for (
int j = 0; j < numPoints; j++) {
408 for (
int k = 0; k < spaceDim; k++) {
411 int l = k + j * spaceDim + i * spaceDim * numPoints;
412 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
414 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
417 *outStream <<
" At multi-index { ";
418 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
419 *outStream <<
"} computed grad component: " << vals(i,j,k)
420 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
427 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D1);
428 for (
int i = 0; i < numFields; i++) {
429 for (
int j = 0; j < numPoints; j++) {
430 for (
int k = 0; k < spaceDim; k++) {
433 int l = k + j * spaceDim + i * spaceDim * numPoints;
434 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
436 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
439 *outStream <<
" At multi-index { ";
440 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
441 *outStream <<
"} computed D1 component: " << vals(i,j,k)
442 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
451 vals.resize(numFields, numPoints, D2cardinality);
452 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D2);
453 for (
int i = 0; i < numFields; i++) {
454 for (
int j = 0; j < numPoints; j++) {
455 for (
int k = 0; k < D2cardinality; k++) {
458 int l = k + j * D2cardinality + i * D2cardinality * numPoints;
459 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
461 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
464 *outStream <<
" At multi-index { ";
465 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
466 *outStream <<
"} computed D2 component: " << vals(i,j,k)
467 <<
" but reference D2 component: " << basisD2[l] <<
"\n";
476 vals.resize(numFields, numPoints, D3cardinality);
477 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D3);
479 for (
int i = 0; i < numFields; i++) {
480 for (
int j = 0; j < numPoints; j++) {
481 for (
int k = 0; k < D3cardinality; k++) {
484 int l = k + j * D3cardinality + i * D3cardinality * numPoints;
485 if (std::abs(vals(i,j,k) - basisD3[l]) > INTREPID_TOL) {
487 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
490 *outStream <<
" At multi-index { ";
491 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
492 *outStream <<
"} computed D3 component: " << vals(i,j,k)
493 <<
" but reference D3 component: " << basisD3[l] <<
"\n";
501 for(EOperator op = OPERATOR_D4; op < OPERATOR_MAX; op++) {
505 vals.resize(numFields, numPoints, DkCardin);
507 wedgeBasis.getValues(vals, wedgeNodes, op);
508 for (
int i = 0; i < vals.size(); i++) {
509 if (std::abs(vals[i]) > INTREPID_TOL) {
511 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
514 std::vector<int> myIndex;
515 vals.getMultiIndex(myIndex,i);
517 *outStream <<
" At multi-index { ";
518 for(
int j = 0; j < vals.rank(); j++) {
519 *outStream << myIndex[j] <<
" ";
521 *outStream <<
"} computed D"<< ord <<
" component: " << vals[i]
522 <<
" but reference D" << ord <<
" component: 0 \n";
529 catch (
const std::logic_error & err) {
530 *outStream << err.what() <<
"\n\n";
535 std::cout <<
"End Result: TEST FAILED\n";
537 std::cout <<
"End Result: TEST PASSED\n";
540 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_WEDGE_I2_FEM class.
int getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.