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_HEX_I1_FEM) |\n" \
94 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 <<
"| 2) Basis values for VALUE and DIV 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 hexNodes(0,0) = -1.0; hexNodes(0,1) = -1.0; hexNodes(0,2) = -1.0;
119 hexNodes(1,0) = 1.0; hexNodes(1,1) = -1.0; hexNodes(1,2) = -1.0;
120 hexNodes(2,0) = 1.0; hexNodes(2,1) = 1.0; hexNodes(2,2) = -1.0;
121 hexNodes(3,0) = -1.0; hexNodes(3,1) = 1.0; hexNodes(3,2) = -1.0;
123 hexNodes(4,0) = -1.0; hexNodes(4,1) = -1.0; hexNodes(4,2) = 1.0;
124 hexNodes(5,0) = 1.0; hexNodes(5,1) = -1.0; hexNodes(5,2) = 1.0;
125 hexNodes(6,0) = 1.0; hexNodes(6,1) = 1.0; hexNodes(6,2) = 1.0;
126 hexNodes(7,0) = -1.0; hexNodes(7,1) = 1.0; hexNodes(7,2) = 1.0;
128 hexNodes(8,0) = 0.0; hexNodes(8,1) = 0.0; hexNodes(8,2) = 0.0;
130 hexNodes(9,0) = 1.0; hexNodes(9,1) = 0.0; hexNodes(9,2) = 0.0;
131 hexNodes(10,0)= -1.0; hexNodes(10,1)= 0.0; hexNodes(10,2)= 0.0;
133 hexNodes(11,0)= 0.0; hexNodes(11,1)= 1.0; hexNodes(11,2)= 0.0;
134 hexNodes(12,0)= 0.0; hexNodes(12,1)= -1.0; hexNodes(12,2)= 0.0;
136 hexNodes(13,0)= 0.0; hexNodes(13,1)= 0.0; hexNodes(13,2)= 1.0;
137 hexNodes(14,0)= 0.0; hexNodes(14,1)= 0.0; hexNodes(14,2)= -1.0;
147 INTREPID_TEST_COMMAND( hexBasis.
getValues(vals, hexNodes, OPERATOR_GRAD), throwCounter, nException );
150 INTREPID_TEST_COMMAND( hexBasis.
getValues(vals, hexNodes, OPERATOR_CURL), throwCounter, nException );
155 INTREPID_TEST_COMMAND( hexBasis.
getDofOrdinal(3,0,0), throwCounter, nException );
157 INTREPID_TEST_COMMAND( hexBasis.
getDofOrdinal(1,1,1), throwCounter, nException );
159 INTREPID_TEST_COMMAND( hexBasis.
getDofOrdinal(0,4,1), throwCounter, nException );
161 INTREPID_TEST_COMMAND( hexBasis.
getDofTag(12), throwCounter, nException );
163 INTREPID_TEST_COMMAND( hexBasis.
getDofTag(-1), throwCounter, nException );
165 #ifdef HAVE_INTREPID_DEBUG
169 INTREPID_TEST_COMMAND( hexBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
173 INTREPID_TEST_COMMAND( hexBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
177 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
181 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals2, hexNodes, OPERATOR_DIV), throwCounter, nException );
185 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException );
189 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals4, hexNodes, OPERATOR_DIV), throwCounter, nException );
193 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals5, hexNodes, OPERATOR_VALUE), throwCounter, nException );
197 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals6, hexNodes, OPERATOR_DIV), throwCounter, nException );
201 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals7, hexNodes, OPERATOR_VALUE), throwCounter, nException );
205 catch (
const std::logic_error & err) {
206 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
207 *outStream << err.what() <<
'\n';
208 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
213 if (throwCounter != nException) {
215 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
220 <<
"===============================================================================\n"\
221 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
222 <<
"===============================================================================\n";
225 std::vector<std::vector<int> > allTags = hexBasis.
getAllDofTags();
228 for (
unsigned i = 0; i < allTags.size(); i++) {
229 int bfOrd = hexBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
231 std::vector<int> myTag = hexBasis.
getDofTag(bfOrd);
232 if( !( (myTag[0] == allTags[i][0]) &&
233 (myTag[1] == allTags[i][1]) &&
234 (myTag[2] == allTags[i][2]) &&
235 (myTag[3] == allTags[i][3]) ) ) {
237 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
238 *outStream <<
" getDofOrdinal( {"
239 << allTags[i][0] <<
", "
240 << allTags[i][1] <<
", "
241 << allTags[i][2] <<
", "
242 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
243 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
247 << myTag[3] <<
"}\n";
253 std::vector<int> myTag = hexBasis.
getDofTag(bfOrd);
254 int myBfOrd = hexBasis.
getDofOrdinal(myTag[0], myTag[1], myTag[2]);
255 if( bfOrd != myBfOrd) {
257 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
258 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
262 << myTag[3] <<
"} but getDofOrdinal({"
266 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
270 catch (
const std::logic_error & err){
271 *outStream << err.what() <<
"\n\n";
277 <<
"===============================================================================\n"\
278 <<
"| TEST 3: correctness of basis function values |\n"\
279 <<
"===============================================================================\n";
281 outStream -> precision(20);
284 double basisValues[] = {
286 0.,-0.25,0., 0.,0.,0., 0.,0.,0., -0.25,0.,0., 0.,0.,-0.25, 0.,0.,0.,
287 0.,-0.25,0., 0.25,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,-0.25, 0.,0.,0.,
288 0.,0.,0., 0.25,0.,0., 0.,0.25,0., 0.,0.,0., 0.,0.,-0.25, 0.,0.,0.,
289 0.,0.,0., 0.,0.,0., 0.,0.25,0., -0.25,0.,0., 0.,0.,-0.25, 0.,0.,0.,
291 0.,-0.25,0., 0.,0.,0., 0.,0.,0., -0.25,0.,0., 0.,0.,0., 0.,0.,0.25,
292 0.,-0.25,0., 0.25,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.25,
293 0.,0.,0., 0.25,0.,0., 0.,0.25,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.25,
294 0.,0.,0., 0.,0.,0., 0.,0.25,0., -0.25,0.,0., 0.,0.,0., 0.,0.,0.25,
296 0.,-0.125,0., 0.125,0.,0., 0.,0.125,0., -0.125,0.,0., 0.,0.,-0.125, 0.,0.,0.125,
298 0.,-0.125,0., 0.25,0.,0., 0.,0.125,0., 0.,0.,0., 0.,0.,-0.125, 0.,0.,0.125,
299 0.,-0.125,0., 0.,0.,0., 0.,0.125,0., -0.25,0.,0., 0.,0.,-0.125, 0.,0.,0.125,
301 0.,0.,0., 0.125,0.,0., 0.,0.25,0., -0.125,0.,0., 0.,0.,-0.125, 0.,0.,0.125,
302 0.,-0.25,0., 0.125,0.,0., 0.,0.,0., -0.125,0.,0., 0.,0.,-0.125, 0.,0.,0.125,
304 0.,-0.125,0., 0.125,0.,0., 0.,0.125,0., -0.125,0.,0., 0.,0.,0., 0.,0.,0.25,
305 0.,-0.125,0., 0.125,0.,0., 0.,0.125,0., -0.125,0.,0., 0.,0.,-0.25, 0.,0.,0.
309 double basisDivs[] = {
311 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
312 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
313 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
314 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
316 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
317 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
318 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
319 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
321 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
323 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
324 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
326 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
327 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
329 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
330 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
336 int numPoints = hexNodes.dimension(0);
344 vals.
resize(numFields, numPoints, spaceDim);
345 hexBasis.
getValues(vals, hexNodes, OPERATOR_VALUE);
346 for (
int i = 0; i < numFields; i++) {
347 for (
int j = 0; j < numPoints; j++) {
348 for (
int k = 0; k < spaceDim; k++) {
351 int l = k + i * spaceDim + j * spaceDim * numFields;
352 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
354 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
357 *outStream <<
" At multi-index { ";
358 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
359 *outStream <<
"} computed value: " << vals(i,j,k)
360 <<
" but reference value: " << basisValues[l] <<
"\n";
367 vals.
resize(numFields, numPoints);
368 hexBasis.
getValues(vals, hexNodes, OPERATOR_DIV);
369 for (
int i = 0; i < numFields; i++) {
370 for (
int j = 0; j < numPoints; j++) {
371 int l = i + j * numFields;
372 if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) {
374 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
377 *outStream <<
" At multi-index { ";
378 *outStream << i <<
" ";*outStream << j <<
" ";
379 *outStream <<
"} computed divergence component: " << vals(i,j)
380 <<
" but reference divergence component: " << basisDivs[l] <<
"\n";
388 catch (
const std::logic_error & err) {
389 *outStream << err.what() <<
"\n\n";
395 <<
"===============================================================================\n"\
396 <<
"| TEST 4: correctness of DoF locations |\n"\
397 <<
"===============================================================================\n";
400 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
402 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
410 #ifdef HAVE_INTREPID_DEBUG
412 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
414 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
416 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
419 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
421 if (throwCounter != nException) {
423 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
428 normals(0,0) = 0.0; normals(0,1) = -4.0; normals(0,2) = 0.0;
429 normals(1,0) = 4.0; normals(1,1) = 0.0; normals(1,2) = 0.0;
430 normals(2,0) = 0.0; normals(2,1) = 4.0; normals(2,2) = 0.0;
431 normals(3,0) = -4.0; normals(3,1) = 0.0; normals(3,2) = 0.0;
432 normals(4,0) = 0.0; normals(4,1) = 0.0; normals(4,2) = -4.0;
433 normals(5,0) = 0.0; normals(5,1) = 0.0; normals(5,2) = 4.0;
435 basis->getValues(bvals, cvals, OPERATOR_VALUE);
437 for (
int i=0; i<bvals.dimension(0); i++) {
438 for (
int j=0; j<bvals.dimension(1); j++) {
441 for(
int d=0;d<spaceDim;d++)
442 normal += bvals(i,j,d)*normals(j,d);
444 if ((i != j) && (std::abs(normal - 0.0) > INTREPID_TOL)) {
446 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);
447 *outStream << buffer;
449 else if ((i == j) && (std::abs(normal - 1.0) > INTREPID_TOL)) {
451 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);
452 *outStream << buffer;
458 catch (
const std::logic_error & err){
459 *outStream << err.what() <<
"\n\n";
464 std::cout <<
"End Result: TEST FAILED\n";
466 std::cout <<
"End Result: TEST PASSED\n";
469 std::cout.copyfmt(oldFormatState);
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron 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...
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_HEX_I1_FEM class.
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedron cell...
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...