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_HDIV_TET_I1_FEM) |\n" \
94 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 <<
"| 3) Exception tests on input arguments and invalid operators |\n" \
96 <<
"| 2) Basis values for VALUE and DIV operators |\n" \
98 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
99 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
100 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
102 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
103 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
105 <<
"===============================================================================\n"\
106 <<
"| TEST 1: Basis creation, exception testing |\n"\
107 <<
"===============================================================================\n";
115 int throwCounter = 0;
119 tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0;
120 tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0;
121 tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0;
122 tetNodes(3,0) = 0.0; tetNodes(3,1) = 0.0; tetNodes(3,2) = 1.0;
123 tetNodes(4,0) = 0.5; tetNodes(4,1) = 0.0; tetNodes(4,2) = 0.0;
124 tetNodes(5,0) = 0.5; tetNodes(5,1) = 0.5; tetNodes(5,2) = 0.0;
125 tetNodes(6,0) = 0.0; tetNodes(6,1) = 0.5; tetNodes(6,2) = 0.0;
126 tetNodes(7,0) = 0.0; tetNodes(7,1) = 0.0; tetNodes(7,2) = 0.5;
127 tetNodes(8,0) = 0.5; tetNodes(8,1) = 0.0; tetNodes(8,2) = 0.5;
128 tetNodes(9,0) = 0.0; tetNodes(9,1) = 0.5; tetNodes(9,2) = 0.5;
138 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, tetNodes, OPERATOR_GRAD), throwCounter, nException );
141 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException );
146 INTREPID_TEST_COMMAND( tetBasis.
getDofOrdinal(3,0,0), throwCounter, nException );
148 INTREPID_TEST_COMMAND( tetBasis.
getDofOrdinal(1,1,1), throwCounter, nException );
150 INTREPID_TEST_COMMAND( tetBasis.
getDofOrdinal(0,4,1), throwCounter, nException );
152 INTREPID_TEST_COMMAND( tetBasis.
getDofTag(7), throwCounter, nException );
154 INTREPID_TEST_COMMAND( tetBasis.
getDofTag(-1), throwCounter, nException );
156 #ifdef HAVE_INTREPID_DEBUG
160 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
164 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
168 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
172 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals2, tetNodes, OPERATOR_VALUE), throwCounter, nException );
176 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
180 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals4, tetNodes, OPERATOR_DIV), throwCounter, nException );
184 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals5, tetNodes, OPERATOR_VALUE), throwCounter, nException );
188 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals6, tetNodes, OPERATOR_DIV), throwCounter, nException );
192 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals7, tetNodes, OPERATOR_VALUE), throwCounter, nException );
196 catch (
const std::logic_error & err) {
197 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
198 *outStream << err.what() <<
'\n';
199 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
204 if (throwCounter != nException) {
206 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
211 <<
"===============================================================================\n"\
212 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
213 <<
"===============================================================================\n";
216 std::vector<std::vector<int> > allTags = tetBasis.
getAllDofTags();
219 for (
unsigned i = 0; i < allTags.size(); i++) {
220 int bfOrd = tetBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
222 std::vector<int> myTag = tetBasis.
getDofTag(bfOrd);
223 if( !( (myTag[0] == allTags[i][0]) &&
224 (myTag[1] == allTags[i][1]) &&
225 (myTag[2] == allTags[i][2]) &&
226 (myTag[3] == allTags[i][3]) ) ) {
228 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
229 *outStream <<
" getDofOrdinal( {"
230 << allTags[i][0] <<
", "
231 << allTags[i][1] <<
", "
232 << allTags[i][2] <<
", "
233 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
234 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
238 << myTag[3] <<
"}\n";
244 std::vector<int> myTag = tetBasis.
getDofTag(bfOrd);
245 int myBfOrd = tetBasis.
getDofOrdinal(myTag[0], myTag[1], myTag[2]);
246 if( bfOrd != myBfOrd) {
248 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
249 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
253 << myTag[3] <<
"} but getDofOrdinal({"
257 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
261 catch (
const std::logic_error & err){
262 *outStream << err.what() <<
"\n\n";
268 <<
"===============================================================================\n"\
269 <<
"| TEST 3: correctness of basis function values |\n"\
270 <<
"===============================================================================\n";
272 outStream -> precision(20);
276 double basisValues[] = {
278 0.,-2.0,0., 0.,0.,0., -2.0,0.,0., 0.,0.,-2.0,
279 2.0,-2.0,0., 2.0,0.,0., 0.,0.,0., 2.0,0.,-2.0,
280 0.,0.,0., 0.,2.0,0., -2.0,2.0,0., 0,2.0,-2.0,
281 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, 0.,0.,0.,
283 1.0,-2.0,0., 1.0,0.,0., -1.0,0.,0., 1.0,0.,-2.0,
284 1.0,-1.0,0., 1.0,1.0,0., -1.0,1.0,0., 1.0,1.0,-2.0,
285 0.,-1.0,0., 0.,1.0,0., -2.0,1.0,0., 0.,1.0,-2.0,
286 0.,-2.0,1.0, 0.,0.,1.0, -2.0,0.,1.0, 0.,0.,-1.0,
287 1.0,-2.0,1.0, 1.0,0.,1.0, -1.0,0.,1.0, 1.0,0.,-1.0,
288 0.,-1.0,1.0, 0.,1.0,1.0, -2.0,1.0,1.0, 0.,1.0,-1.0
293 double basisDivs[] = {
312 int numPoints = tetNodes.dimension(0);
319 vals.
resize(numFields, numPoints, spaceDim);
320 tetBasis.
getValues(vals, tetNodes, OPERATOR_VALUE);
321 for (
int i = 0; i < numFields; i++) {
322 for (
int j = 0; j < numPoints; j++) {
323 for (
int k = 0; k < spaceDim; k++) {
325 int l = k + i * spaceDim + j * spaceDim * numFields;
326 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
328 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
331 *outStream <<
" At (Field,Point,Dim) multi-index { ";
332 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
333 *outStream <<
"} computed value: " << vals(i,j,k)
334 <<
" but reference value: " << basisValues[l] <<
"\n";
342 vals.
resize(numFields, numPoints);
343 tetBasis.
getValues(vals, tetNodes, OPERATOR_DIV);
344 for (
int i = 0; i < numFields; i++) {
345 for (
int j = 0; j < numPoints; j++) {
346 int l = i + j * numFields;
347 if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) {
349 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
352 *outStream <<
" At multi-index { ";
353 *outStream << i <<
" ";*outStream << j <<
" ";
354 *outStream <<
"} computed divergence component: " << vals(i,j)
355 <<
" but reference divergence component: " << basisDivs[l] <<
"\n";
363 catch (
const std::logic_error & err) {
364 *outStream << err.what() <<
"\n\n";
370 <<
"===============================================================================\n"\
371 <<
"| TEST 4: correctness of DoF locations |\n"\
372 <<
"===============================================================================\n";
375 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
377 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
385 #ifdef HAVE_INTREPID_DEBUG
387 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
389 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
391 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
394 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
396 if (throwCounter != nException) {
398 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
403 normals(0,0) = 0.0; normals(0,1) = -0.5; normals(0,2) = 0.0;
404 normals(1,0) = 0.5; normals(1,1) = 0.5; normals(1,2) = 0.5;
405 normals(2,0) = -0.5; normals(2,1) = 0.0; normals(2,2) = 0.0;
406 normals(3,0) = 0.0; normals(3,1) = 0.0; normals(3,2) = -0.5;
408 basis->getValues(bvals, cvals, OPERATOR_VALUE);
410 for (
int i=0; i<bvals.dimension(0); i++) {
411 for (
int j=0; j<bvals.dimension(1); j++) {
414 for(
int d=0;d<spaceDim;d++)
415 normal += bvals(i,j,d)*normals(j,d);
417 if ((i != j) && (std::abs(normal - 0.0) > INTREPID_TOL)) {
419 sprintf(buffer,
"\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), normal, 0.0);
420 *outStream << buffer;
422 else if ((i == j) && (std::abs(normal - 1.0) > INTREPID_TOL)) {
424 sprintf(buffer,
"\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), normal, 1.0);
425 *outStream << buffer;
431 catch (
const std::logic_error & err){
432 *outStream << err.what() <<
"\n\n";
437 std::cout <<
"End Result: TEST FAILED\n";
439 std::cout <<
"End Result: TEST PASSED\n";
442 std::cout.copyfmt(oldFormatState);
Implementation of the default H(div)-compatible FEM basis of degree 1 on Tetrahedron cell...
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
virtual int getCardinality() const
Returns cardinality of the basis.
This is an interface class for bases whose degrees of freedom can be associated with spatial location...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Tetrahedron cell.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HDIV_TET_I1_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.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...