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_HEX_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";
114 int throwCounter = 0;
119 hexNodes(0, 0) = -1.0; hexNodes(0, 1) = -1.0; hexNodes(0, 2) = -1.0;
120 hexNodes(1, 0) = 1.0; hexNodes(1, 1) = -1.0; hexNodes(1, 2) = -1.0;
121 hexNodes(2, 0) = 1.0; hexNodes(2, 1) = 1.0; hexNodes(2, 2) = -1.0;
122 hexNodes(3, 0) = -1.0; hexNodes(3, 1) = 1.0; hexNodes(3, 2) = -1.0;
124 hexNodes(4, 0) = -1.0; hexNodes(4, 1) = -1.0; hexNodes(4, 2) = 1.0;
125 hexNodes(5, 0) = 1.0; hexNodes(5, 1) = -1.0; hexNodes(5, 2) = 1.0;
126 hexNodes(6, 0) = 1.0; hexNodes(6, 1) = 1.0; hexNodes(6, 2) = 1.0;
127 hexNodes(7, 0) = -1.0; hexNodes(7, 1) = 1.0; hexNodes(7, 2) = 1.0;
130 hexNodes(8, 0) = 0.0; hexNodes(8, 1) = -1.0; hexNodes(8, 2) = -1.0;
131 hexNodes(9, 0) = 1.0; hexNodes(9, 1) = 0.0; hexNodes(9, 2) = -1.0;
132 hexNodes(10,0) = 0.0; hexNodes(10,1) = 1.0; hexNodes(10,2) = -1.0;
133 hexNodes(11,0) = -1.0; hexNodes(11,1) = 0.0; hexNodes(11,2) = -1.0;
134 hexNodes(12,0) = -1.0; hexNodes(12,1) = -1.0; hexNodes(12,2) = 0.0;
135 hexNodes(13,0) = 1.0; hexNodes(13,1) = -1.0; hexNodes(13,2) = 0.0;
136 hexNodes(14,0) = 1.0; hexNodes(14,1) = 1.0; hexNodes(14,2) = 0.0;
137 hexNodes(15,0) = -1.0; hexNodes(15,1) = 1.0; hexNodes(15,2) = 0.0;
138 hexNodes(16,0) = 0.0; hexNodes(16,1) = -1.0; hexNodes(16,2) = 1.0;
139 hexNodes(17,0) = 1.0; hexNodes(17,1) = 0.0; hexNodes(17,2) = 1.0;
140 hexNodes(18,0) = 0.0; hexNodes(18,1) = 1.0; hexNodes(18,2) = 1.0;
141 hexNodes(19,0) = -1.0; hexNodes(19,1) = 0.0; hexNodes(19,2) = 1.0;
150 INTREPID_TEST_COMMAND( hexBasis.
getValues(vals, hexNodes, OPERATOR_CURL), throwCounter, nException );
155 INTREPID_TEST_COMMAND( hexBasis.
getValues(vals, hexNodes, OPERATOR_DIV), throwCounter, nException );
160 INTREPID_TEST_COMMAND( hexBasis.
getDofOrdinal(3,10,0), throwCounter, nException );
162 INTREPID_TEST_COMMAND( hexBasis.
getDofOrdinal(1,2,1), throwCounter, nException );
164 INTREPID_TEST_COMMAND( hexBasis.
getDofOrdinal(0,4,1), throwCounter, nException );
166 INTREPID_TEST_COMMAND( hexBasis.
getDofTag(21), throwCounter, nException );
168 INTREPID_TEST_COMMAND( hexBasis.
getDofTag(-1), throwCounter, nException );
170#ifdef HAVE_INTREPID_DEBUG
174 INTREPID_TEST_COMMAND( hexBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
178 INTREPID_TEST_COMMAND( hexBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
182 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
186 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals2, hexNodes, OPERATOR_GRAD), throwCounter, nException );
189 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals2, hexNodes, OPERATOR_D1), throwCounter, nException );
192 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals2, hexNodes, OPERATOR_D2), throwCounter, nException );
196 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException );
200 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException );
204 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals5, hexNodes, OPERATOR_GRAD), throwCounter, nException );
208 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals6, hexNodes, OPERATOR_D2), throwCounter, nException );
212 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals7, hexNodes, OPERATOR_D3), throwCounter, nException );
216 catch (
const std::logic_error & err) {
217 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
218 *outStream << err.what() <<
'\n';
219 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
225 if (throwCounter != nException) {
227 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
232 <<
"===============================================================================\n"\
233 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
234 <<
"===============================================================================\n";
237 std::vector<std::vector<int> > allTags = hexBasis.
getAllDofTags();
240 for (
unsigned i = 0; i < allTags.size(); i++) {
241 int bfOrd = hexBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
243 std::vector<int> myTag = hexBasis.
getDofTag(bfOrd);
244 if( !( (myTag[0] == allTags[i][0]) &&
245 (myTag[1] == allTags[i][1]) &&
246 (myTag[2] == allTags[i][2]) &&
247 (myTag[3] == allTags[i][3]) ) ) {
249 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
250 *outStream <<
" getDofOrdinal( {"
251 << allTags[i][0] <<
", "
252 << allTags[i][1] <<
", "
253 << allTags[i][2] <<
", "
254 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
255 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
259 << myTag[3] <<
"}\n";
265 std::vector<int> myTag = hexBasis.
getDofTag(bfOrd);
266 int myBfOrd = hexBasis.
getDofOrdinal(myTag[0], myTag[1], myTag[2]);
267 if( bfOrd != myBfOrd) {
269 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
270 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
274 << myTag[3] <<
"} but getDofOrdinal({"
278 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
282 catch (
const std::logic_error & err){
283 *outStream << err.what() <<
"\n\n";
290 <<
"===============================================================================\n"\
291 <<
"| TEST 3: correctness of basis function values |\n"\
292 <<
"===============================================================================\n";
294 outStream -> precision(20);
297 double basisValues[] = {
298 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
299 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
300 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
301 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
302 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
303 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
304 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
305 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
306 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
307 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
308 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
309 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, \
310 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, \
311 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, \
312 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, \
313 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, \
314 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, \
315 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, \
316 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, \
317 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0 };
321 std::string fileName;
322 std::ifstream dataFile;
326 std::vector<double> basisGrads;
328 fileName =
"./testdata/HEX_I2_GradVals.dat";
329 dataFile.open(fileName.c_str());
330 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
331 ">>> ERROR (HGRAD_HEX_I2/test01): could not open GRAD values data file, test aborted.");
332 while (!dataFile.eof() ){
335 std::getline(dataFile, line);
336 stringstream data_line(line);
337 while(data_line >> temp){
338 basisGrads.push_back(temp);
349 std::vector<double> basisD2;
350 fileName =
"./testdata/HEX_I2_D2Vals.dat";
351 dataFile.open(fileName.c_str());
352 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
353 ">>> ERROR (HGRAD_HEX_I2/test01): could not open D2 values data file, test aborted.");
354 while (!dataFile.eof() ){
357 std::getline(dataFile, line);
358 stringstream data_line(line);
359 while(data_line >> temp){
360 basisD2.push_back(temp);
368 std::vector<double> basisD3;
370 fileName =
"./testdata/HEX_I2_D3Vals.dat";
371 dataFile.open(fileName.c_str());
372 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
373 ">>> ERROR (HGRAD_HEX_I2/test01): could not open D3 values data file, test aborted.");
375 while (!dataFile.eof() ){
378 std::getline(dataFile, line);
379 stringstream data_line(line);
380 while(data_line >> temp){
381 basisD3.push_back(temp);
393 int numPoints = hexNodes.dimension(0);
400 vals.
resize(numFields, numPoints);
401 hexBasis.
getValues(vals, hexNodes, OPERATOR_VALUE);
402 for (
int i = 0; i < numFields; i++) {
403 for (
int j = 0; j < numPoints; j++) {
404 int l = i + j * numFields;
405 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
407 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
410 *outStream <<
" At multi-index { ";
411 *outStream << i <<
" ";*outStream << j <<
" ";
412 *outStream <<
"} computed value: " << vals(i,j)
413 <<
" but reference value: " << basisValues[l] <<
"\n";
420 vals.
resize(numFields, numPoints, spaceDim);
421 hexBasis.
getValues(vals, hexNodes, OPERATOR_GRAD);
422 for (
int i = 0; i < numFields; i++) {
423 for (
int j = 0; j < numPoints; j++) {
424 for (
int k = 0; k < spaceDim; k++) {
427 int l = k + j * spaceDim + i * spaceDim * numPoints;
429 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
431 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
434 *outStream <<
" At multi-index { ";
435 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
436 *outStream <<
"} computed grad component: " << vals(i,j,k)
437 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
445 hexBasis.
getValues(vals, hexNodes, OPERATOR_D1);
446 for (
int i = 0; i < numFields; i++) {
447 for (
int j = 0; j < numPoints; j++) {
448 for (
int k = 0; k < spaceDim; k++) {
451 int l = k + j * spaceDim + i * spaceDim * numPoints;
452 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
454 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
457 *outStream <<
" At multi-index { ";
458 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
459 *outStream <<
"} computed D1 component: " << vals(i,j,k)
460 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
470 vals.
resize(numFields, numPoints, D2cardinality);
471 hexBasis.
getValues(vals, hexNodes, OPERATOR_D2);
472 for (
int i = 0; i < numFields; i++) {
473 for (
int j = 0; j < numPoints; j++) {
474 for (
int k = 0; k < D2cardinality; k++) {
477 int l = k + j * D2cardinality + i * D2cardinality * numPoints;
478 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
480 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
483 *outStream <<
" At multi-index { ";
484 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
485 *outStream <<
"} computed D2 component: " << vals(i,j,k)
486 <<
" but reference D2 component: " << basisD2[l] <<
"\n";
495 vals.
resize(numFields, numPoints, D3cardinality);
496 hexBasis.
getValues(vals, hexNodes, OPERATOR_D3);
498 for (
int i = 0; i < numFields; i++) {
499 for (
int j = 0; j < numPoints; j++) {
500 for (
int k = 0; k < D3cardinality; k++) {
503 int l = k + j * D3cardinality + i * D3cardinality * numPoints;
504 if (std::abs(vals(i,j,k) - basisD3[l]) > INTREPID_TOL) {
506 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
509 *outStream <<
" At multi-index { ";
510 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
511 *outStream <<
"} computed D3 component: " << vals(i,j,k)
512 <<
" but reference D3 component: " << basisD3[l] <<
"\n";
521 catch (
const std::logic_error & err) {
522 *outStream << err.what() <<
"\n\n";
528 <<
"===============================================================================\n"\
529 <<
"| TEST 4: correctness of DoF locations |\n"\
530 <<
"===============================================================================\n";
533 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
535 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
536 Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
542#ifdef HAVE_INTREPID_DEBUG
544 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
546 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
548 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
551 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
553 if (throwCounter != nException) {
555 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
559 basis->getValues(bvals, cvals, OPERATOR_VALUE);
561 for (
int i=0; i<bvals.dimension(0); i++) {
562 for (
int j=0; j<bvals.dimension(1); j++) {
563 if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) {
565 sprintf(buffer,
"\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), bvals(i,j), 0.0);
566 *outStream << buffer;
568 else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) {
570 sprintf(buffer,
"\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), bvals(i,j), 1.0);
571 *outStream << buffer;
577 catch (
const std::logic_error & err){
578 *outStream << err.what() <<
"\n\n";
583 std::cout <<
"End Result: TEST FAILED\n";
585 std::cout <<
"End Result: TEST PASSED\n";
588 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_HEX_I2_FEM class.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
Implementation of the serendipity-family H(grad)-compatible FEM basis of degree 2 on a Hexahedron cel...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
virtual int getDofOrdinal(const int subcDim, const int subcOrd, const int subcDofOrd)
DoF tag to ordinal lookup.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
virtual int getCardinality() const
Returns cardinality of the basis.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.