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_WEDGE_C2_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;
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;
132 wedgeNodes(12,0)= 0.5; wedgeNodes(12,1)= 0.0; wedgeNodes(12,2)= 1.0;
133 wedgeNodes(13,0)= 0.5; wedgeNodes(13,1)= 0.5; wedgeNodes(13,2)= 1.0;
134 wedgeNodes(14,0)= 0.0; wedgeNodes(14,1)= 0.5; wedgeNodes(14,2)= 1.0;
135 wedgeNodes(15,0)= 0.5; wedgeNodes(15,1)= 0.0; wedgeNodes(15,2)= 0.0;
136 wedgeNodes(16,0)= 0.5; wedgeNodes(16,1)= 0.5; wedgeNodes(16,2)= 0.0;
137 wedgeNodes(17,0)= 0.0; wedgeNodes(17,1)= 0.5; wedgeNodes(17,2)= 0.0;
147 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
152 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
157 INTREPID_TEST_COMMAND( wedgeBasis.
getDofOrdinal(3,0,0), throwCounter, nException );
159 INTREPID_TEST_COMMAND( wedgeBasis.
getDofOrdinal(1,1,1), throwCounter, nException );
161 INTREPID_TEST_COMMAND( wedgeBasis.
getDofOrdinal(0,9,0), throwCounter, nException );
163 INTREPID_TEST_COMMAND( wedgeBasis.
getDofTag(18), throwCounter, nException );
165 INTREPID_TEST_COMMAND( wedgeBasis.
getDofTag(-1), throwCounter, nException );
167 #ifdef HAVE_INTREPID_DEBUG 171 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
175 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
179 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
183 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals2, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
186 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals2, wedgeNodes, OPERATOR_D1), throwCounter, nException );
189 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals2, wedgeNodes, OPERATOR_D2), throwCounter, nException );
193 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
197 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals4, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
201 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals5, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
205 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals6, wedgeNodes, OPERATOR_D2), throwCounter, nException );
208 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals6, wedgeNodes, OPERATOR_D3), throwCounter, nException );
212 catch (
const std::logic_error & err) {
213 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
214 *outStream << err.what() <<
'\n';
215 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
220 if (throwCounter != nException) {
222 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
227 <<
"===============================================================================\n"\
228 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
229 <<
"===============================================================================\n";
232 std::vector<std::vector<int> > allTags = wedgeBasis.
getAllDofTags();
235 for (
unsigned i = 0; i < allTags.size(); i++) {
236 int bfOrd = wedgeBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
238 std::vector<int> myTag = wedgeBasis.
getDofTag(bfOrd);
239 if( !( (myTag[0] == allTags[i][0]) &&
240 (myTag[1] == allTags[i][1]) &&
241 (myTag[2] == allTags[i][2]) &&
242 (myTag[3] == allTags[i][3]) ) ) {
244 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
245 *outStream <<
" getDofOrdinal( {" 246 << allTags[i][0] <<
", " 247 << allTags[i][1] <<
", " 248 << allTags[i][2] <<
", " 249 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
250 *outStream <<
" getDofTag(" << bfOrd <<
") = { " 254 << myTag[3] <<
"}\n";
259 for(
int bfOrd = 0; bfOrd < wedgeBasis.
getCardinality(); bfOrd++) {
260 std::vector<int> myTag = wedgeBasis.
getDofTag(bfOrd);
261 int myBfOrd = wedgeBasis.
getDofOrdinal(myTag[0], myTag[1], myTag[2]);
262 if( bfOrd != myBfOrd) {
264 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
265 *outStream <<
" getDofTag(" << bfOrd <<
") = { " 269 << myTag[3] <<
"} but getDofOrdinal({" 273 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
277 catch (
const std::logic_error & err){
278 *outStream << err.what() <<
"\n\n";
284 <<
"===============================================================================\n"\
285 <<
"| TEST 3: correctness of basis function values |\n"\
286 <<
"===============================================================================\n";
288 outStream -> precision(20);
291 double basisValues[] = {
292 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, \
293 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, \
294 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, \
295 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
296 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
297 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
298 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
299 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, \
300 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, \
301 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, \
302 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
303 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
304 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
305 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
306 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00 };
309 std::string fileName;
310 std::ifstream dataFile;
313 std::vector<double> basisGrads;
315 fileName =
"./testdata/WEDGE_C2_GradVals.dat";
316 dataFile.open(fileName.c_str());
317 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
318 ">>> ERROR (HGRAD_WEDGE_C2/test01): could not open GRAD values data file, test aborted.");
319 while (!dataFile.eof() ){
322 std::getline(dataFile, line);
323 stringstream data_line(line);
324 while(data_line >> temp){
325 basisGrads.push_back(temp);
336 std::vector<double> basisD2;
338 fileName =
"./testdata/WEDGE_C2_D2Vals.dat";
339 dataFile.open(fileName.c_str());
340 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
341 ">>> ERROR (HGRAD_WEDGE_C2/test01): could not open D2 values data file, test aborted.");
343 while (!dataFile.eof() ){
346 std::getline(dataFile, line);
347 stringstream data_line(line);
348 while(data_line >> temp){
349 basisD2.push_back(temp);
357 std::vector<double> basisD3;
359 fileName =
"./testdata/WEDGE_C2_D3Vals.dat";
360 dataFile.open(fileName.c_str());
361 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
362 ">>> ERROR (HGRAD_WEDGE_C2/test01): could not open D3 values data file, test aborted.");
364 while (!dataFile.eof() ){
367 std::getline(dataFile, line);
368 stringstream data_line(line);
369 while(data_line >> temp){
370 basisD3.push_back(temp);
378 std::vector<double> basisD4;
380 fileName =
"./testdata/WEDGE_C2_D4Vals.dat";
381 dataFile.open(fileName.c_str());
382 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
383 ">>> ERROR (HGRAD_WEDGE_C2/test01): could not open D4 values data file, test aborted.");
385 while (!dataFile.eof() ){
388 std::getline(dataFile, line);
389 stringstream data_line(line);
390 while(data_line >> temp){
391 basisD4.push_back(temp);
402 int numPoints = wedgeNodes.dimension(0);
409 vals.
resize(numFields, numPoints);
410 wedgeBasis.
getValues(vals, wedgeNodes, OPERATOR_VALUE);
411 for (
int i = 0; i < numFields; i++) {
412 for (
int j = 0; j < numPoints; j++) {
413 int l = i + j * numFields;
414 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
416 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
419 *outStream <<
" At multi-index { ";
420 *outStream << i <<
" ";*outStream << j <<
" ";
421 *outStream <<
"} computed value: " << vals(i,j)
422 <<
" but reference value: " << basisValues[l] <<
"\n";
430 vals.
resize(numFields, numPoints, spaceDim);
431 wedgeBasis.
getValues(vals, wedgeNodes, OPERATOR_GRAD);
432 for (
int i = 0; i < numFields; i++) {
433 for (
int j = 0; j < numPoints; j++) {
434 for (
int k = 0; k < spaceDim; k++) {
437 int l = k + j * spaceDim + i * spaceDim * numPoints;
438 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
440 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
443 *outStream <<
" At multi-index { ";
444 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
445 *outStream <<
"} computed grad component: " << vals(i,j,k)
446 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
453 wedgeBasis.
getValues(vals, wedgeNodes, OPERATOR_D1);
454 for (
int i = 0; i < numFields; i++) {
455 for (
int j = 0; j < numPoints; j++) {
456 for (
int k = 0; k < spaceDim; k++) {
459 int l = k + j * spaceDim + i * spaceDim * numPoints;
460 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
462 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
465 *outStream <<
" At multi-index { ";
466 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
467 *outStream <<
"} computed D1 component: " << vals(i,j,k)
468 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
477 vals.
resize(numFields, numPoints, D2cardinality);
478 wedgeBasis.
getValues(vals, wedgeNodes, OPERATOR_D2);
479 for (
int i = 0; i < numFields; i++) {
480 for (
int j = 0; j < numPoints; j++) {
481 for (
int k = 0; k < D2cardinality; k++) {
484 int l = k + j * D2cardinality + i * D2cardinality * numPoints;
485 if (std::abs(vals(i,j,k) - basisD2[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 D2 component: " << vals(i,j,k)
493 <<
" but reference D2 component: " << basisD2[l] <<
"\n";
502 vals.
resize(numFields, numPoints, D3cardinality);
503 wedgeBasis.
getValues(vals, wedgeNodes, OPERATOR_D3);
505 for (
int i = 0; i < numFields; i++) {
506 for (
int j = 0; j < numPoints; j++) {
507 for (
int k = 0; k < D3cardinality; k++) {
510 int l = k + j * D3cardinality + i * D3cardinality * numPoints;
511 if (std::abs(vals(i,j,k) - basisD3[l]) > INTREPID_TOL) {
513 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
516 *outStream <<
" At multi-index { ";
517 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
518 *outStream <<
"} computed D3 component: " << vals(i,j,k)
519 <<
" but reference D3 component: " << basisD3[l] <<
"\n";
528 vals.
resize(numFields, numPoints, D4cardinality);
529 wedgeBasis.
getValues(vals, wedgeNodes, OPERATOR_D4);
530 for (
int i = 0; i < numFields; i++) {
531 for (
int j = 0; j < numPoints; j++) {
532 for (
int k = 0; k < D4cardinality; k++) {
535 int l = k + j * D4cardinality + i * D4cardinality * numPoints;
536 if (std::abs(vals(i,j,k) - basisD4[l]) > INTREPID_TOL) {
538 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
541 *outStream <<
" At multi-index { ";
542 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
543 *outStream <<
"} computed D4 component: " << vals(i,j,k)
544 <<
" but reference D4 component: " << basisD2[l] <<
"\n";
552 for(
EOperator op = OPERATOR_D5; op < OPERATOR_MAX; op++) {
556 vals.
resize(numFields, numPoints, DkCardin);
558 wedgeBasis.
getValues(vals, wedgeNodes, op);
559 for (
int i = 0; i < vals.
size(); i++) {
560 if (std::abs(vals[i]) > INTREPID_TOL) {
562 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
565 std::vector<int> myIndex;
568 *outStream <<
" At multi-index { ";
569 for(
int j = 0; j < vals.
rank(); j++) {
570 *outStream << myIndex[j] <<
" ";
572 *outStream <<
"} computed D"<< ord <<
" component: " << vals[i]
573 <<
" but reference D" << ord <<
" component: 0 \n";
580 catch (
const std::logic_error & err) {
581 *outStream << err.what() <<
"\n\n";
586 std::cout <<
"End Result: TEST FAILED\n";
588 std::cout <<
"End Result: TEST PASSED\n";
591 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) ...
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Wedge cell.
void getMultiIndex(int &i0, const int valueEnum) const
Returns the multi-index of a value, based on its enumeration, as a list, for rank-1 containers...
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Header file for utility class to provide multidimensional containers.
int rank() const
Return rank of the FieldContainer = number of indices used to tag the multi-indexed value...
int getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Wedge cell.
Header file for the Intrepid::G_WEDGE_C2_FEM class.
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...
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...