50 #include "Teuchos_oblackholestream.hpp"
51 #include "Teuchos_RCP.hpp"
52 #include "Teuchos_GlobalMPISession.hpp"
55 using 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"; \
71 int 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";
469 int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
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";
494 int D3cardinality = Intrepid::getDkCardinality(OPERATOR_D3, spaceDim);
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 =
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);
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
Implementation of the serendipity-family H(grad)-compatible FEM basis of degree 2 on a Hexahedron cel...
virtual int getCardinality() const
Returns cardinality of the basis.
This is an interface class for bases whose degrees of freedom can be associated with spatial location...
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_HEX_I2_FEM class.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
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 shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...