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_HCURL_HEX_I1_FEM) |\n" \
94 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 <<
"| 2) Basis values for VALUE and CURL 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 );
152 INTREPID_TEST_COMMAND( hexBasis.
getValues(vals, hexNodes, OPERATOR_DIV), throwCounter, nException );
157 INTREPID_TEST_COMMAND( hexBasis.
getDofOrdinal(3,0,0), throwCounter, nException );
159 INTREPID_TEST_COMMAND( hexBasis.
getDofOrdinal(1,1,1), throwCounter, nException );
161 INTREPID_TEST_COMMAND( hexBasis.
getDofOrdinal(0,4,1), throwCounter, nException );
163 INTREPID_TEST_COMMAND( hexBasis.
getDofTag(12), throwCounter, nException );
165 INTREPID_TEST_COMMAND( hexBasis.
getDofTag(-1), throwCounter, nException );
167 #ifdef HAVE_INTREPID_DEBUG
171 INTREPID_TEST_COMMAND( hexBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
175 INTREPID_TEST_COMMAND( hexBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
179 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
182 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals1, hexNodes, OPERATOR_CURL), throwCounter, nException );
186 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals2, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
190 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
194 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
197 INTREPID_TEST_COMMAND( hexBasis.
getValues(badVals4, hexNodes, OPERATOR_CURL), throwCounter, nException ) ;
202 hexNodes.dimension(0),
204 INTREPID_TEST_COMMAND( hexBasis.
getValues(vals, hexNodes, OPERATOR_D2), throwCounter, nException );
209 catch (
const std::logic_error & err) {
210 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
211 *outStream << err.what() <<
'\n';
212 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
218 if (throwCounter != nException) {
220 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
226 <<
"===============================================================================\n"\
227 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
228 <<
"===============================================================================\n";
231 std::vector<std::vector<int> > allTags = hexBasis.
getAllDofTags();
234 for (
unsigned i = 0; i < allTags.size(); i++) {
235 int bfOrd = hexBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
237 std::vector<int> myTag = hexBasis.
getDofTag(bfOrd);
238 if( !( (myTag[0] == allTags[i][0]) &&
239 (myTag[1] == allTags[i][1]) &&
240 (myTag[2] == allTags[i][2]) &&
241 (myTag[3] == allTags[i][3]) ) ) {
243 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
244 *outStream <<
" getDofOrdinal( {"
245 << allTags[i][0] <<
", "
246 << allTags[i][1] <<
", "
247 << allTags[i][2] <<
", "
248 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
249 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
253 << myTag[3] <<
"}\n";
259 std::vector<int> myTag = hexBasis.
getDofTag(bfOrd);
260 int myBfOrd = hexBasis.
getDofOrdinal(myTag[0], myTag[1], myTag[2]);
261 if( bfOrd != myBfOrd) {
263 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
264 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
268 << myTag[3] <<
"} but getDofOrdinal({"
272 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
276 catch (
const std::logic_error & err){
277 *outStream << err.what() <<
"\n\n";
283 <<
"===============================================================================\n"\
284 <<
"| TEST 3: correctness of basis function values |\n"\
285 <<
"===============================================================================\n";
287 outStream -> precision(20);
290 double basisValues[] = {
292 0.5,0.,0., 0.,0.,0., 0.,0.,0., 0.,-0.5,0., 0.,0.,0., 0.,0.,0.,
293 0.,0.,0., 0.,0.,0., 0.,0.,0.5, 0.,0.,0., 0.,0.,0., 0.,0.,0.,
295 0.5,0.,0., 0.,0.5,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.,
296 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.5, 0.,0.,0., 0.,0.,0.,
298 0.,0.,0., 0.,0.5,0., -0.5,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.,
299 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.5, 0.,0.,0.,
301 0.,0.,0., 0.,0.,0., -0.5,0.,0., 0.,-0.5,0., 0.,0.,0., 0.,0.,0.,
302 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.5,
305 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.5,0.,0., 0.,0.,0.,
306 0.,0.,0., 0.,-0.5,0., 0.,0.,0.5, 0.,0.,0., 0.,0.,0., 0.,0.,0.,
308 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.5,0.,0., 0.,0.5,0.,
309 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.5, 0.,0.,0., 0.,0.,0.,
311 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.5,0.,
312 -0.5,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.5, 0.,0.,0.,
314 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.,
315 -0.5,0.,0., 0.,-0.5,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.5,
318 0.125,0.,0., 0.,0.125,0., -0.125,0.,0., 0.,-0.125,0., 0.125,0.,0., 0.,0.125,0.,
319 -0.125,0.,0., 0.,-0.125,0., 0.,0.,0.125, 0.,0.,0.125, 0.,0.,0.125, 0.,0.,0.125,
322 0.125,0.,0., 0.,0.25,0., -0.125,0.,0., 0.,0.,0., 0.125,0.,0., 0.,0.25,0.,
323 -0.125,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.25, 0.,0.,0.25, 0.,0.,0.,
325 0.125,0.,0., 0.,0.,0., -0.125,0.,0., 0.,-0.25,0., 0.125,0.,0., 0.,0.,0.,
326 -0.125,0.,0., 0.,-0.25,0., 0.,0.,0.25, 0.,0.,0., 0.,0.,0., 0.,0.,0.25,
329 0.,0.,0., 0.,0.125,0., -0.25,0.,0., 0.,-0.125,0., 0.,0.,0., 0.,0.125,0.,
330 -0.25,0.,0., 0.,-0.125,0., 0.,0.,0., 0.,0.,0., 0.,0.,0.25, 0.,0.,0.25,
332 0.25,0.,0., 0.,0.125,0., 0.,0.,0., 0.,-0.125,0., 0.25,0.,0., 0.,0.125,0.,
333 0.,0.,0., 0.,-0.125,0., 0.,0.,0.25, 0.,0.,0.25, 0.,0.,0., 0.,0.,0.,
336 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.,0., 0.25,0.,0., 0.,0.25,0.,
337 -0.25,0.,0., 0.,-0.25,0., 0.,0.,0.125, 0.,0.,0.125, 0.,0.,0.125, 0.,0.,0.125,
339 0.25,0.,0., 0.,0.25,0., -0.25,0.,0., 0.,-0.25,0., 0.,0.,0., 0.,0.,0.,
340 0.0,0.,0., 0.,0.,0.0, 0.,0.,0.125, 0.,0.,0.125, 0.,0.,0.125, 0.,0.,0.125
344 double basisCurls[] = {
346 0.,-0.25,0.25, 0.,0.,0.25, 0.,0.,0.25, -0.25,0.,0.25, 0.,0.25,0., 0.,0.,0.,
347 0.,0.,0., 0.25,0.,0., -0.25,0.25,0., 0.,-0.25,0., 0.,0.,0., 0.25,0.,0.,
349 0.,-0.25,0.25, 0.25,0.,0.25, 0.,0.,0.25, 0.,0.,0.25, 0.,0.25,0., -0.25,0.,0.,
350 0.,0.,0., 0.,0.,0., 0.,0.25,0., -0.25,-0.25,0., 0.25,0.,0., 0.,0.,0.,
352 0.,0.,0.25, 0.25,0.,0.25, 0.,0.25,0.25, 0.,0.,0.25, 0.,0.,0., -0.25,0.,0.,
353 0.,-0.25,0., 0.,0.,0., 0.,0.,0., -0.25,0.,0., 0.25,-0.25,0., 0.,0.25,0.,
355 0.,0.,0.25, 0.,0.,0.25, 0.,0.25,0.25, -0.25,0.,0.25, 0.,0.,0., 0.,0.,0.,
356 0.,-0.25,0., 0.25,0.,0., -0.25,0.,0., 0.,0.,0., 0.,-0.25,0., 0.25,0.25,0.,
359 0.,-0.25,0., 0.,0.,0., 0.,0.,0., -0.25,0.,0., 0.,0.25,0.25, 0.,0.,0.25,
360 0.,0.,0.25, 0.25,0.,0.25, -0.25,0.25,0., 0.,-0.25,0., 0.,0.,0., 0.25,0.,0.,
362 0.,-0.25,0., 0.25,0.,0., 0.,0.,0., 0.,0.,0., 0.,0.25,0.25, -0.25,0.,0.25,
363 0.,0.,0.25, 0.,0.,0.25, 0.,0.25,0., -0.25,-0.25,0., 0.25,0.,0., 0.,0.,0.,
365 0.,0.,0., 0.25,0.,0., 0.,0.25,0., 0.,0.,0., 0.,0.,0.25, -0.25,0.,0.25,
366 0.,-0.25,0.25, 0.,0.,0.25, 0.,0.,0., -0.25,0.,0., 0.25,-0.25,0., 0.,0.25,0.,
368 0.,0.,0., 0.,0.,0., 0.,0.25,0., -0.25,0.,0., 0.,0.,0.25, 0.,0.,0.25,
369 0.,-0.25,0.25, 0.25,0.,0.25, -0.25,0.,0., 0.,0.,0., 0.,-0.25,0., 0.25,0.25,0.,
372 0.,-0.125,0.125, 0.125,0.,0.125, 0.,0.125,0.125, -0.125,0.,0.125, 0.,0.125,0.125, -0.125,0.,0.125,
373 0.,-0.125,0.125, 0.125,0.,0.125, -0.125,0.125,0., -0.125,-0.125,0., 0.125,-0.125,0., 0.125,0.125,0.,
376 0.,-0.125,0.125, 0.25,0.,0.125, 0.,0.125,0.125, 0.,0.,0.125, 0.,0.125,0.125, -0.25,0.,0.125,
377 0.,-0.125,0.125, 0.,0.,0.125, 0.,0.125,0., -0.25,-0.125,0., 0.25,-0.125,0., 0.,0.125,0.,
379 0.,-0.125,0.125, 0.,0.,0.125, 0.,0.125,0.125, -0.25,0.,0.125, 0.,0.125,0.125, 0.,0.,0.125,
380 0.,-0.125,0.125, 0.25,0.,0.125, -0.25,0.125,0., 0.,-0.125,0., 0.,-0.125,0., 0.25,0.125,0.,
383 0.,0.,0.125, 0.125,0.,0.125, 0.,0.25,0.125, -0.125,0.,0.125, 0.,0.,0.125, -0.125,0.,0.125,
384 0.,-0.25,0.125, 0.125,0.,0.125, -0.125,0.,0., -0.125,0.,0., 0.125,-0.25,0., 0.125,0.25,0.,
386 0.,-0.25,0.125, 0.125,0.,0.125, 0.,0.,0.125, -0.125,0.,0.125, 0.,0.25,0.125, -0.125,0.,0.125,
387 0.,0.,0.125, 0.125,0.,0.125, -0.125,0.25,0., -0.125,-0.25,0., 0.125,0.,0., 0.125,0.,0.,
390 0.,-0.125,0., 0.125,0.,0., 0.,0.125,0., -0.125,0.,0., 0.,0.125,0.25, -0.125,0.,0.25,
391 0.,-0.125,0.25, 0.125,0.,0.25, -0.125,0.125,0., -0.125,-0.125,0., 0.125,-0.125,0., 0.125,0.125,0.,
393 0.,-0.125,0.25, 0.125,0.,0.25, 0.,0.125,0.25, -0.125,0.,0.25, 0.,0.125,0., -0.125,0.,0.,
394 0.,-0.125,0., 0.125,0.,0., -0.125,0.125,0., -0.125,-0.125,0., 0.125,-0.125,0., 0.125,0.125,0.
402 int numPoints = hexNodes.dimension(0);
409 vals.
resize(numFields, numPoints, spaceDim);
410 hexBasis.
getValues(vals, hexNodes, OPERATOR_VALUE);
411 for (
int i = 0; i < numFields; i++) {
412 for (
int j = 0; j < numPoints; j++) {
413 for (
int k = 0; k < spaceDim; k++) {
416 int l = k + i * spaceDim + j * spaceDim * numFields;
417 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
419 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
422 *outStream <<
" At multi-index { ";
423 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
424 *outStream <<
"} computed value: " << vals(i,j,k)
425 <<
" but reference value: " << basisValues[l] <<
"\n";
432 vals.
resize(numFields, numPoints, spaceDim);
433 hexBasis.
getValues(vals, hexNodes, OPERATOR_CURL);
434 for (
int i = 0; i < numFields; i++) {
435 for (
int j = 0; j < numPoints; j++) {
436 for (
int k = 0; k < spaceDim; k++) {
439 int l = k + i * spaceDim + j * spaceDim * numFields;
440 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
442 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
445 *outStream <<
" At multi-index { ";
446 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
447 *outStream <<
"} computed curl component: " << vals(i,j,k)
448 <<
" but reference curl component: " << basisCurls[l] <<
"\n";
457 catch (
const std::logic_error & err) {
458 *outStream << err.what() <<
"\n\n";
464 <<
"===============================================================================\n"\
465 <<
"| TEST 4: correctness of DoF locations |\n"\
466 <<
"===============================================================================\n";
469 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
471 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
479 #ifdef HAVE_INTREPID_DEBUG
481 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
483 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
485 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
487 cvals.
resize(12,spaceDim);
488 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
490 if (throwCounter != nException) {
492 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
497 tangents(0,0) = 2.0; tangents(0,1) = 0.0; tangents(0,2) = 0.0;
498 tangents(1,0) = 0.0; tangents(1,1) = 2.0; tangents(1,2) = 0.0;
499 tangents(2,0) = -2.0; tangents(2,1) = 0.0; tangents(2,2) = 0.0;
500 tangents(3,0) = 0.0; tangents(3,1) = -2.0; tangents(3,2) = 0.0;
501 tangents(4,0) = 2.0; tangents(4,1) = 0.0; tangents(4,2) = 0.0;
502 tangents(5,0) = 0.0; tangents(5,1) = 2.0; tangents(5,2) = 0.0;
503 tangents(6,0) = -2.0; tangents(6,1) = 0.0; tangents(6,2) = 0.0;
504 tangents(7,0) = 0.0; tangents(7,1) = -2.0; tangents(7,2) = 0.0;
505 tangents(8,0) = 0.0; tangents(8,1) = 0.0; tangents(8,2) = 2.0;
506 tangents(9,0) = 0.0; tangents(9,1) = 0.0; tangents(9,2) = 2.0;
507 tangents(10,0) = 0.0; tangents(10,1) = 0.0; tangents(10,2) = 2.0;
508 tangents(11,0) = 0.0; tangents(11,1) = 0.0; tangents(11,2) = 2.0;
510 basis->getValues(bvals, cvals, OPERATOR_VALUE);
512 for (
int i=0; i<bvals.dimension(0); i++) {
513 for (
int j=0; j<bvals.dimension(1); j++) {
515 double tangent = 0.0;
516 for(
int d=0;d<spaceDim;d++)
517 tangent += bvals(i,j,d)*tangents(j,d);
519 if ((i != j) && (std::abs(tangent - 0.0) > INTREPID_TOL)) {
521 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), tangent, 0.0);
522 *outStream << buffer;
524 else if ((i == j) && (std::abs(tangent - 1.0) > INTREPID_TOL)) {
526 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), tangent, 1.0);
527 *outStream << buffer;
533 catch (
const std::logic_error & err){
534 *outStream << err.what() <<
"\n\n";
539 std::cout <<
"End Result: TEST FAILED\n";
541 std::cout <<
"End Result: TEST PASSED\n";
544 std::cout.copyfmt(oldFormatState);
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.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
Implementation of the default H(curl)-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...
Header file for the Intrepid::HCURL_HEX_I1_FEM class.