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";
114 int throwCounter = 0;
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;
135 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException );
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
159 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
163 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
167 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
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 );
181 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
185 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals4, tetNodes, OPERATOR_VALUE), throwCounter, nException );
189 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals5, tetNodes, OPERATOR_GRAD), throwCounter, nException );
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";
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,
306 int numPoints = tetNodes.dimension(0);
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);
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;
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.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Tetrahedron cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Tetrahedron cell.
virtual int getDofOrdinal(const int subcDim, const int subcOrd, const int subcDofOrd)
DoF tag to ordinal lookup.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
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 a templated lexicographical container for a multi-indexed scalar quantity....
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
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.