50 #include "Teuchos_oblackholestream.hpp" 51 #include "Teuchos_RCP.hpp" 52 #include "Teuchos_GlobalMPISession.hpp" 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_PYR_I2_FEM) |\n" \
94 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 <<
"| 2) Basis values for VALUE, GRAD, CURL, 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;
118 pyrNodes(0,0) = -1.0; pyrNodes(0,1) = -1.0; pyrNodes(0,2) = 0;
119 pyrNodes(1,0) = 1.0; pyrNodes(1,1) = -1.0; pyrNodes(1,2) = 0;
120 pyrNodes(2,0) = 1.0; pyrNodes(2,1) = 1.0; pyrNodes(2,2) = 0;
121 pyrNodes(3,0) = -1.0; pyrNodes(3,1) = 1.0; pyrNodes(3,2) = 0;
122 pyrNodes(4,0) = 0.0; pyrNodes(4,1) = 0.0; pyrNodes(4,2) = 1.0;
123 pyrNodes(5,0) = 0.0; pyrNodes(5,1) = -1.0; pyrNodes(5,2) = 0.0;
124 pyrNodes(6,0) = 1.0; pyrNodes(6,1) = 0.0; pyrNodes(6,2) = 0.0;
125 pyrNodes(7,0) = 0.0; pyrNodes(7,1) = 1.0; pyrNodes(7,2) = 0.0;
126 pyrNodes(8,0) = -1.0; pyrNodes(8,1) = 0.0; pyrNodes(8,2) = 0.0;
127 pyrNodes(9,0) = -0.5; pyrNodes(9,1) = -0.5; pyrNodes(9,2) = 0.5;
128 pyrNodes(10,0) = 0.5; pyrNodes(10,1) = -0.5; pyrNodes(10,2) = 0.5;
129 pyrNodes(11,0) = 0.5; pyrNodes(11,1) = 0.5; pyrNodes(11,2) = 0.5;
130 pyrNodes(12,0) = -0.5; pyrNodes(12,1) = 0.5; pyrNodes(12,2) = 0.5;
132 pyrNodes(13,0) = 0.25; pyrNodes(13,1) = 0.5; pyrNodes(13,2) = 0.2;
133 pyrNodes(14,0) = -0.7 ; pyrNodes(14,1) = 0.3; pyrNodes(14,2) = 0.3;
134 pyrNodes(15,0) = 0.; pyrNodes(15,1) = -0.05; pyrNodes(15,2) = 0.95;
135 pyrNodes(16,0) = -0.15; pyrNodes(16,1) = -0.2; pyrNodes(16,2) = 0.75;
136 pyrNodes(17,0) = -0.4; pyrNodes(17,1) = 0.9; pyrNodes(17,2) = 0.0;
146 INTREPID_TEST_COMMAND( pyrBasis.
getValues(vals, pyrNodes, OPERATOR_CURL), throwCounter, nException );
151 INTREPID_TEST_COMMAND( pyrBasis.
getValues(vals, pyrNodes, OPERATOR_DIV), throwCounter, nException );
156 INTREPID_TEST_COMMAND( pyrBasis.
getDofOrdinal(3,0,0), throwCounter, nException );
158 INTREPID_TEST_COMMAND( pyrBasis.
getDofOrdinal(1,1,1), throwCounter, nException );
160 INTREPID_TEST_COMMAND( pyrBasis.
getDofOrdinal(0,6,0), throwCounter, nException );
162 INTREPID_TEST_COMMAND( pyrBasis.
getDofTag(14), throwCounter, nException );
164 INTREPID_TEST_COMMAND( pyrBasis.
getDofTag(-1), throwCounter, nException );
166 #ifdef HAVE_INTREPID_DEBUG 170 INTREPID_TEST_COMMAND( pyrBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
174 INTREPID_TEST_COMMAND( pyrBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
178 INTREPID_TEST_COMMAND( pyrBasis.
getValues(badVals1, pyrNodes, OPERATOR_VALUE), throwCounter, nException );
182 INTREPID_TEST_COMMAND( pyrBasis.
getValues(badVals2, pyrNodes, OPERATOR_GRAD), throwCounter, nException );
185 INTREPID_TEST_COMMAND( pyrBasis.
getValues(badVals2, pyrNodes, OPERATOR_D1), throwCounter, nException );
192 INTREPID_TEST_COMMAND( pyrBasis.
getValues(badVals3, pyrNodes, OPERATOR_VALUE), throwCounter, nException );
196 INTREPID_TEST_COMMAND( pyrBasis.
getValues(badVals4, pyrNodes, OPERATOR_VALUE), throwCounter, nException );
200 INTREPID_TEST_COMMAND( pyrBasis.
getValues(badVals5, pyrNodes, OPERATOR_GRAD), throwCounter, nException );
211 catch (
const std::logic_error & err) {
212 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
213 *outStream << err.what() <<
'\n';
214 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
219 if (throwCounter != nException) {
221 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
226 <<
"===============================================================================\n"\
227 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
228 <<
"===============================================================================\n";
231 std::vector<std::vector<int> > allTags = pyrBasis.
getAllDofTags();
234 for (
unsigned i = 0; i < allTags.size(); i++) {
235 int bfOrd = pyrBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
237 std::vector<int> myTag = pyrBasis.
getDofTag(bfOrd);
238 if( !( (myTag[0] == allTags[i][0]) &&
239 (myTag[1] == allTags[i][1]) &&
240 (myTag[2] == allTags[i][2]) &&
241 (myTag[3] == allTags[i][3]) ) ) {
243 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
244 *outStream <<
" getDofOrdinal( {" 245 << allTags[i][0] <<
", " 246 << allTags[i][1] <<
", " 247 << allTags[i][2] <<
", " 248 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
249 *outStream <<
" getDofTag(" << bfOrd <<
") = { " 253 << myTag[3] <<
"}\n";
259 std::vector<int> myTag = pyrBasis.
getDofTag(bfOrd);
260 int myBfOrd = pyrBasis.
getDofOrdinal(myTag[0], myTag[1], myTag[2]);
261 if( bfOrd != myBfOrd) {
263 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
264 *outStream <<
" getDofTag(" << bfOrd <<
") = { " 268 << myTag[3] <<
"} but getDofOrdinal({" 272 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
276 catch (
const std::logic_error & err){
277 *outStream << err.what() <<
"\n\n";
283 <<
"===============================================================================\n"\
284 <<
"| TEST 3: correctness of basis function values |\n"\
285 <<
"===============================================================================\n";
287 outStream -> precision(20);
290 std::string fileName;
291 std::ifstream dataFile;
294 std::vector<double> basisValues;
295 fileName =
"./testdata/PYR_I2_Vals.dat";
296 dataFile.open(fileName.c_str());
297 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
298 ">>> ERROR (HGRAD_PYR_I2/test01): could not open basis values data file, test aborted.");
299 while (!dataFile.eof() ){
302 std::getline(dataFile, line);
303 stringstream data_line(line);
304 while(data_line >> temp){
305 basisValues.push_back(temp);
316 std::vector<double> basisGrads;
317 fileName =
"./testdata/PYR_I2_GradVals.dat";
318 dataFile.open(fileName.c_str());
319 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
320 ">>> ERROR (HGRAD_PYR_I2/test01): could not open GRAD values data file, test aborted.");
321 while (!dataFile.eof() ){
324 std::getline(dataFile, line);
325 stringstream data_line(line);
326 while(data_line >> temp){
327 basisGrads.push_back(temp);
337 std::vector<double> basisD2;
339 fileName =
"./testdata/PYR_I2_D2Vals.dat";
340 dataFile.open(fileName.c_str());
341 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
342 ">>> ERROR (HGRAD_PYR_I2/test01): could not open D2 values data file, test aborted.");
344 while (!dataFile.eof() ){
347 std::getline(dataFile, line);
348 stringstream data_line(line);
349 while(data_line >> temp){
350 basisD2.push_back(temp);
360 int numPoints = pyrNodes.dimension(0);
368 vals.
resize(numFields, numPoints);
369 pyrBasis.
getValues(vals, pyrNodes, OPERATOR_VALUE);
370 for (
int i = 0; i < numFields; i++) {
371 for (
int j = 0; j < numPoints; j++) {
372 int l = j + i * numPoints;
374 if (std::abs(vals(i,j) - basisValues[l]) > 100.*INTREPID_TOL) {
376 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
379 *outStream <<
" At multi-index { ";
380 *outStream << i <<
" ";*outStream << j <<
" ";
381 *outStream <<
"} computed value: " << vals(i,j)
382 <<
" but reference value: " << basisValues[l] <<
"\n";
390 vals.
resize(numFields, numPoints, spaceDim);
391 pyrBasis.
getValues(vals, pyrNodes, OPERATOR_GRAD);
392 for (
int i = 0; i < numFields; i++) {
393 for (
int j = 0; j < numPoints; j++) {
394 for (
int k = 0; k < spaceDim; k++) {
396 int l = i + j * numFields + k * numFields * numPoints;
398 if (std::abs(vals(i,j,k) - basisGrads[l]) > 100.*INTREPID_TOL) {
400 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
403 *outStream <<
" At multi-index { ";
404 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
405 *outStream <<
"} computed grad component: " << vals(i,j,k)
406 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
415 pyrBasis.
getValues(vals, pyrNodes, OPERATOR_D1);
416 for (
int i = 0; i < numFields; i++) {
417 for (
int j = 0; j < numPoints; j++) {
418 for (
int k = 0; k < spaceDim; k++) {
419 int l = i + j * numFields + k * numFields * numPoints;
421 if (std::abs(vals(i,j,k) - basisGrads[l]) > 100.*INTREPID_TOL) {
423 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
426 *outStream <<
" At multi-index { ";
427 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
428 *outStream <<
"} computed D1 component: " << vals(i,j,k)
429 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
438 vals.
resize(numFields, numPoints, D2Cardin);
439 pyrBasis.
getValues(vals, pyrNodes, OPERATOR_D2);
441 for (
int i = 0; i < numFields; i++) {
442 for (
int j = 0; j < numPoints; j++) {
444 for (
int k = 0; k < D2Cardin; k++) {
446 int l = i + j * numFields + k * numFields * numPoints;
450 if (std::abs(vals(i,j,k) - basisD2[l]) > 100.*INTREPID_TOL) {
452 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
455 *outStream <<
" At multi-index { ";
456 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
457 *outStream <<
"} computed D2 component: " << vals(i,j,k)
458 <<
" but reference D2 component: " << basisD2[l] <<
"\n";
473 catch (
const std::logic_error & err) {
474 *outStream << err.what() <<
"\n\n";
479 std::cout <<
"End Result: TEST FAILED\n";
481 std::cout <<
"End Result: TEST PASSED\n";
484 std::cout.copyfmt(oldFormatState);
virtual int getCardinality() const
Returns cardinality of the basis.
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
int main(int argc, char *argv[])
Performs a code-code comparison to FIAT for Nedelec bases on tets (values and curls) ...
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Header file for utility class to provide multidimensional containers.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
Implementation of an H(grad)-compatible FEM basis of degree 2 on a Pyramid cell.
virtual int getDofOrdinal(const int subcDim, const int subcOrd, const int subcDofOrd)
DoF tag to ordinal lookup.
Header file for the Intrepid::HGRAD_PYR_I2_FEM class.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Pyramid cell.