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_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";
476 int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
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";
501 int D3cardinality = Intrepid::getDkCardinality(OPERATOR_D3, spaceDim);
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";
527 int D4cardinality = Intrepid::getDkCardinality(OPERATOR_D4, spaceDim);
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++) {
555 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
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;
567 int ord = Intrepid::getOperatorOrder(op);
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);
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
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 the default H(grad)-compatible FEM basis of degree 2 on Wedge cell.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Header file for utility class to provide multidimensional containers.
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...
int rank() const
Return rank of the FieldContainer = number of indices used to tag the multi-indexed value...
Header file for the Intrepid::G_WEDGE_C2_FEM class.
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.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Wedge cell.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...