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;
137 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, tetNodes, OPERATOR_GRAD), throwCounter, nException );
140 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, tetNodes, OPERATOR_CURL), 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,1), throwCounter, nException );
151 INTREPID_TEST_COMMAND( tetBasis.
getDofTag(7), 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_VALUE), throwCounter, nException );
175 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
179 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals4, tetNodes, OPERATOR_DIV), throwCounter, nException );
183 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals5, tetNodes, OPERATOR_VALUE), throwCounter, nException );
187 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals6, tetNodes, OPERATOR_DIV), throwCounter, nException );
191 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals7, tetNodes, OPERATOR_VALUE), throwCounter, nException );
195 catch (
const std::logic_error & err) {
196 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
197 *outStream << err.what() <<
'\n';
198 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
203 if (throwCounter != nException) {
205 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
210 <<
"===============================================================================\n"\
211 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
212 <<
"===============================================================================\n";
215 std::vector<std::vector<int> > allTags = tetBasis.
getAllDofTags();
218 for (
unsigned i = 0; i < allTags.size(); i++) {
219 int bfOrd = tetBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
221 std::vector<int> myTag = tetBasis.
getDofTag(bfOrd);
222 if( !( (myTag[0] == allTags[i][0]) &&
223 (myTag[1] == allTags[i][1]) &&
224 (myTag[2] == allTags[i][2]) &&
225 (myTag[3] == allTags[i][3]) ) ) {
227 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
228 *outStream <<
" getDofOrdinal( {"
229 << allTags[i][0] <<
", "
230 << allTags[i][1] <<
", "
231 << allTags[i][2] <<
", "
232 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
233 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
237 << myTag[3] <<
"}\n";
243 std::vector<int> myTag = tetBasis.
getDofTag(bfOrd);
244 int myBfOrd = tetBasis.
getDofOrdinal(myTag[0], myTag[1], myTag[2]);
245 if( bfOrd != myBfOrd) {
247 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
248 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
252 << myTag[3] <<
"} but getDofOrdinal({"
256 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
260 catch (
const std::logic_error & err){
261 *outStream << err.what() <<
"\n\n";
267 <<
"===============================================================================\n"\
268 <<
"| TEST 3: correctness of basis function values |\n"\
269 <<
"===============================================================================\n";
271 outStream -> precision(20);
275 double basisValues[] = {
277 0.,-2.0,0., 0.,0.,0., -2.0,0.,0., 0.,0.,-2.0,
278 2.0,-2.0,0., 2.0,0.,0., 0.,0.,0., 2.0,0.,-2.0,
279 0.,0.,0., 0.,2.0,0., -2.0,2.0,0., 0,2.0,-2.0,
280 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, 0.,0.,0.,
282 1.0,-2.0,0., 1.0,0.,0., -1.0,0.,0., 1.0,0.,-2.0,
283 1.0,-1.0,0., 1.0,1.0,0., -1.0,1.0,0., 1.0,1.0,-2.0,
284 0.,-1.0,0., 0.,1.0,0., -2.0,1.0,0., 0.,1.0,-2.0,
285 0.,-2.0,1.0, 0.,0.,1.0, -2.0,0.,1.0, 0.,0.,-1.0,
286 1.0,-2.0,1.0, 1.0,0.,1.0, -1.0,0.,1.0, 1.0,0.,-1.0,
287 0.,-1.0,1.0, 0.,1.0,1.0, -2.0,1.0,1.0, 0.,1.0,-1.0
292 double basisDivs[] = {
311 int numPoints = tetNodes.dimension(0);
318 vals.
resize(numFields, numPoints, spaceDim);
319 tetBasis.
getValues(vals, tetNodes, OPERATOR_VALUE);
320 for (
int i = 0; i < numFields; i++) {
321 for (
int j = 0; j < numPoints; j++) {
322 for (
int k = 0; k < spaceDim; k++) {
324 int l = k + i * spaceDim + j * spaceDim * numFields;
325 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
327 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
330 *outStream <<
" At (Field,Point,Dim) multi-index { ";
331 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
332 *outStream <<
"} computed value: " << vals(i,j,k)
333 <<
" but reference value: " << basisValues[l] <<
"\n";
341 vals.
resize(numFields, numPoints);
342 tetBasis.
getValues(vals, tetNodes, OPERATOR_DIV);
343 for (
int i = 0; i < numFields; i++) {
344 for (
int j = 0; j < numPoints; j++) {
345 int l = i + j * numFields;
346 if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) {
348 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
351 *outStream <<
" At multi-index { ";
352 *outStream << i <<
" ";*outStream << j <<
" ";
353 *outStream <<
"} computed divergence component: " << vals(i,j)
354 <<
" but reference divergence component: " << basisDivs[l] <<
"\n";
362 catch (
const std::logic_error & err) {
363 *outStream << err.what() <<
"\n\n";
369 <<
"===============================================================================\n"\
370 <<
"| TEST 4: correctness of DoF locations |\n"\
371 <<
"===============================================================================\n";
374 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
376 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
384 #ifdef HAVE_INTREPID_DEBUG
386 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
388 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
390 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
393 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
395 if (throwCounter != nException) {
397 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
402 normals(0,0) = 0.0; normals(0,1) = -0.5; normals(0,2) = 0.0;
403 normals(1,0) = 0.5; normals(1,1) = 0.5; normals(1,2) = 0.5;
404 normals(2,0) = -0.5; normals(2,1) = 0.0; normals(2,2) = 0.0;
405 normals(3,0) = 0.0; normals(3,1) = 0.0; normals(3,2) = -0.5;
407 basis->getValues(bvals, cvals, OPERATOR_VALUE);
409 for (
int i=0; i<bvals.dimension(0); i++) {
410 for (
int j=0; j<bvals.dimension(1); j++) {
413 for(
int d=0;d<spaceDim;d++)
414 normal += bvals(i,j,d)*normals(j,d);
416 if ((i != j) && (std::abs(normal - 0.0) > INTREPID_TOL)) {
418 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);
419 *outStream << buffer;
421 else if ((i == j) && (std::abs(normal - 1.0) > INTREPID_TOL)) {
423 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);
424 *outStream << buffer;
430 catch (
const std::logic_error & err){
431 *outStream << err.what() <<
"\n\n";
436 std::cout <<
"End Result: TEST FAILED\n";
438 std::cout <<
"End Result: TEST PASSED\n";
441 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...