49#include "Intrepid_HGRAD_TET_C1_FEM.hpp"
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_TET_C1_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";
109 Basis_HGRAD_TET_C1_FEM<double, FieldContainer<double> > tetBasis;
114 int throwCounter = 0;
117 FieldContainer<double> tetNodes(8, 3);
118 tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0;
119 tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0;
120 tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0;
121 tetNodes(3,0) = 0.0; tetNodes(3,1) = 0.0; tetNodes(3,2) = 1.0;
122 tetNodes(4,0) = 0.25; tetNodes(4,1) = 0.5; tetNodes(4,2) = 0.1;
123 tetNodes(5,0) = 0.5; tetNodes(5,1) = 0.5; tetNodes(5,2) = 0.0;
124 tetNodes(6,0) = 0.5; tetNodes(6,1) = 0.0; tetNodes(6,2) = 0.5;
125 tetNodes(7,0) = 0.0; tetNodes(7,1) = 0.5; tetNodes(7,2) = 0.5;
129 FieldContainer<double> vals;
134 vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0), 4 );
135 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException );
139 vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0) );
140 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_DIV), throwCounter, nException );
145 INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(3,0,0), throwCounter, nException );
147 INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(1,1,1), throwCounter, nException );
149 INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(0,4,0), throwCounter, nException );
151 INTREPID_TEST_COMMAND( tetBasis.getDofTag(5), throwCounter, nException );
153 INTREPID_TEST_COMMAND( tetBasis.getDofTag(-1), throwCounter, nException );
155#ifdef HAVE_INTREPID_DEBUG
158 FieldContainer<double> badPoints1(4, 5, 3);
159 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
162 FieldContainer<double> badPoints2(4, tetBasis.getBaseCellTopology().getDimension() - 1);
163 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
166 FieldContainer<double> badVals1(4, 3, 1);
167 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
170 FieldContainer<double> badVals2(4, 3);
171 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_GRAD), throwCounter, nException );
174 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_D1), throwCounter, nException );
177 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_D2), throwCounter, nException );
180 FieldContainer<double> badVals3(tetBasis.getCardinality() + 1, tetNodes.dimension(0));
181 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
184 FieldContainer<double> badVals4(tetBasis.getCardinality(), tetNodes.dimension(0) + 1);
185 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals4, tetNodes, OPERATOR_VALUE), throwCounter, nException );
188 FieldContainer<double> badVals5(tetBasis.getCardinality(), tetNodes.dimension(0), tetBasis.getBaseCellTopology().getDimension() + 1);
189 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals5, tetNodes, OPERATOR_GRAD), throwCounter, nException );
192 FieldContainer<double> badVals6(tetBasis.getCardinality(), tetNodes.dimension(0), 40);
193 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_D1), throwCounter, nException );
196 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_D2), throwCounter, nException );
200 catch (
const std::logic_error & err) {
201 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
202 *outStream << err.what() <<
'\n';
203 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
208 if (throwCounter != nException) {
210 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
215 <<
"===============================================================================\n"\
216 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
217 <<
"===============================================================================\n";
220 std::vector<std::vector<int> > allTags = tetBasis.getAllDofTags();
223 for (
unsigned i = 0; i < allTags.size(); i++) {
224 int bfOrd = tetBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
226 std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
227 if( !( (myTag[0] == allTags[i][0]) &&
228 (myTag[1] == allTags[i][1]) &&
229 (myTag[2] == allTags[i][2]) &&
230 (myTag[3] == allTags[i][3]) ) ) {
232 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
233 *outStream <<
" getDofOrdinal( {"
234 << allTags[i][0] <<
", "
235 << allTags[i][1] <<
", "
236 << allTags[i][2] <<
", "
237 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
238 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
242 << myTag[3] <<
"}\n";
247 for(
int bfOrd = 0; bfOrd < tetBasis.getCardinality(); bfOrd++) {
248 std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
249 int myBfOrd = tetBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
250 if( bfOrd != myBfOrd) {
252 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
253 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
257 << myTag[3] <<
"} but getDofOrdinal({"
261 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
265 catch (
const std::logic_error & err){
266 *outStream << err.what() <<
"\n\n";
272 <<
"===============================================================================\n"\
273 <<
"| TEST 3: correctness of basis function values |\n"\
274 <<
"===============================================================================\n";
276 outStream -> precision(20);
279 double basisValues[] = {
291 double basisGrads[] = {
292 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
293 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
294 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
295 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
296 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
297 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
298 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
299 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
305 int numFields = tetBasis.getCardinality();
306 int numPoints = tetNodes.dimension(0);
307 int spaceDim = tetBasis.getBaseCellTopology().getDimension();
310 FieldContainer<double> vals;
313 vals.resize(numFields, numPoints);
314 tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE);
315 for (
int i = 0; i < numFields; i++) {
316 for (
int j = 0; j < numPoints; j++) {
317 int l = i + j * numFields;
318 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
320 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
323 *outStream <<
" At multi-index { ";
324 *outStream << i <<
" ";*outStream << j <<
" ";
325 *outStream <<
"} computed value: " << vals(i,j)
326 <<
" but reference value: " << basisValues[l] <<
"\n";
332 vals.resize(numFields, numPoints, spaceDim);
333 tetBasis.getValues(vals, tetNodes, OPERATOR_GRAD);
334 for (
int i = 0; i < numFields; i++) {
335 for (
int j = 0; j < numPoints; j++) {
336 for (
int k = 0; k < spaceDim; k++) {
337 int l = k + i * spaceDim + j * spaceDim * numFields;
338 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
340 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
343 *outStream <<
" At multi-index { ";
344 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
345 *outStream <<
"} computed grad component: " << vals(i,j,k)
346 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
353 tetBasis.getValues(vals, tetNodes, OPERATOR_D1);
354 for (
int i = 0; i < numFields; i++) {
355 for (
int j = 0; j < numPoints; j++) {
356 for (
int k = 0; k < spaceDim; k++) {
357 int l = k + i * spaceDim + j * spaceDim * numFields;
358 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
360 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
363 *outStream <<
" At multi-index { ";
364 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
365 *outStream <<
"} computed D1 component: " << vals(i,j,k)
366 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
373 for(EOperator op = OPERATOR_D2; op < OPERATOR_MAX; op++) {
377 vals.resize(numFields, numPoints, DkCardin);
379 tetBasis.getValues(vals, tetNodes, op);
380 for (
int i = 0; i < vals.size(); i++) {
381 if (std::abs(vals[i]) > INTREPID_TOL) {
383 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
386 std::vector<int> myIndex;
387 vals.getMultiIndex(myIndex,i);
389 *outStream <<
" At multi-index { ";
390 for(
int j = 0; j < vals.rank(); j++) {
391 *outStream << myIndex[j] <<
" ";
393 *outStream <<
"} computed D"<< ord <<
" component: " << vals[i]
394 <<
" but reference D" << ord <<
" component: 0 \n";
401 catch (
const std::logic_error & err) {
402 *outStream << err.what() <<
"\n\n";
407 std::cout <<
"End Result: TEST FAILED\n";
409 std::cout <<
"End Result: TEST PASSED\n";
412 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.
int getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.