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_TET_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) or |\n" \
98 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
100 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
101 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
103 <<
"===============================================================================\n"\
104 <<
"| TEST 1: Basis creation, exception testing |\n"\
105 <<
"===============================================================================\n";
113 int throwCounter = 0;
117 tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0;
118 tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0;
119 tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0;
120 tetNodes(3,0) = 0.0; tetNodes(3,1) = 0.0; tetNodes(3,2) = 1.0;
121 tetNodes(4,0) = 0.5; tetNodes(4,1) = 0.0; tetNodes(4,2) = 0.0;
122 tetNodes(5,0) = 0.5; tetNodes(5,1) = 0.5; tetNodes(5,2) = 0.0;
123 tetNodes(6,0) = 0.0; tetNodes(6,1) = 0.5; tetNodes(6,2) = 0.0;
124 tetNodes(7,0) = 0.0; tetNodes(7,1) = 0.0; tetNodes(7,2) = 0.5;
125 tetNodes(8,0) = 0.5; tetNodes(8,1) = 0.0; tetNodes(8,2) = 0.5;
126 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 );
142 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, tetNodes, OPERATOR_DIV), throwCounter, nException );
147 INTREPID_TEST_COMMAND( tetBasis.
getDofOrdinal(3,0,0), throwCounter, nException );
149 INTREPID_TEST_COMMAND( tetBasis.
getDofOrdinal(1,1,1), throwCounter, nException );
151 INTREPID_TEST_COMMAND( tetBasis.
getDofOrdinal(0,4,1), throwCounter, nException );
153 INTREPID_TEST_COMMAND( tetBasis.
getDofTag(7), throwCounter, nException );
155 INTREPID_TEST_COMMAND( tetBasis.
getDofTag(-1), throwCounter, nException );
157 #ifdef HAVE_INTREPID_DEBUG
161 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals,badPoints1,OPERATOR_VALUE), throwCounter, nException );
165 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals,badPoints2,OPERATOR_VALUE), throwCounter, nException );
169 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals1,tetNodes,OPERATOR_VALUE), throwCounter, nException );
172 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals1,tetNodes,OPERATOR_CURL), throwCounter, nException );
176 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals2,tetNodes,OPERATOR_VALUE), throwCounter, nException );
180 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals3,tetNodes,OPERATOR_VALUE), throwCounter, nException );
184 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals4,tetNodes,OPERATOR_VALUE), throwCounter, nException );
187 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals4,tetNodes,OPERATOR_CURL), throwCounter, nException );
191 catch (
const std::logic_error & err) {
192 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
193 *outStream << err.what() <<
'\n';
194 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
200 if (throwCounter != nException) {
202 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
207 <<
"===============================================================================\n"\
208 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
209 <<
"===============================================================================\n";
212 std::vector<std::vector<int> > allTags = tetBasis.
getAllDofTags();
215 for (
unsigned i = 0; i < allTags.size(); i++) {
216 int bfOrd = tetBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
218 std::vector<int> myTag = tetBasis.
getDofTag(bfOrd);
219 if( !( (myTag[0] == allTags[i][0]) &&
220 (myTag[1] == allTags[i][1]) &&
221 (myTag[2] == allTags[i][2]) &&
222 (myTag[3] == allTags[i][3]) ) ) {
224 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
225 *outStream <<
" getDofOrdinal( {"
226 << allTags[i][0] <<
", "
227 << allTags[i][1] <<
", "
228 << allTags[i][2] <<
", "
229 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
230 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
234 << myTag[3] <<
"}\n";
240 std::vector<int> myTag = tetBasis.
getDofTag(bfOrd);
241 int myBfOrd = tetBasis.
getDofOrdinal(myTag[0], myTag[1], myTag[2]);
242 if( bfOrd != myBfOrd) {
244 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
245 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
249 << myTag[3] <<
"} but getDofOrdinal({"
253 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
257 catch (
const std::logic_error & err){
258 *outStream << err.what() <<
"\n\n";
264 <<
"===============================================================================\n"\
265 <<
"| TEST 3: correctness of basis function values |\n"\
266 <<
"===============================================================================\n";
268 outStream -> precision(20);
271 double basisValues[] = {
273 1.0,0.,0., 0.,0.,0., 0.,-1.0,0., 0.,0.,1.0,
276 1.0,1.0,1.0, 0.,1.,0., 0.,0.,0., 0.,0.,0.,
279 0.,0.,0., -1.,0.,0., -1.0,-1.0,-1.0,
280 0.,0.,0., 0.,0.,0., 0.,0.,1.,
282 0.,0.,0., 0.,0.,0., 0.,0.,0., 1.0,1.0,1.0,
283 -1.,0.,0., 0.,-1.,0.,
286 1.0,0.5,0.5, 0.,0.5,0., 0.,-0.5,0.,
287 0.,0.,0.5, 0.,0.,0.5, 0.,0.,0.,
289 0.5,0.5,0.5, -0.5,0.5,0.,
290 -0.5,-0.5,-0.5, 0.,0.,0., 0.,0.,0.5, 0.,0.,0.5,
292 0.5,0.,0., -0.5,0.,0., -0.5,-1.0,-0.5,
293 0.,0.,0.5, 0.,0.,0., 0.,0.,0.5,
295 0.5,0.,0., 0.,0.,0., 0.,-0.5,0., 0.5,0.5,1.0,
296 -0.5,0.,0., 0.,-0.5,0.,
298 0.5,0.5,0.5, 0.,0.5,0., 0.,0.,0., 0.5,0.5,0.5,
299 -0.5,0.,0.5, 0.,-0.5,0.,
301 0.,0.,0., -0.5,0.,0., -0.5,-0.5,-0.5, 0.5,0.5,0.5,
302 -0.5,0.,0., 0.,-0.5,0.5
306 double basisCurls[] = {
308 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, -2.0,2.0,0.,
309 0.,-2.0,0., 2.0,0.,0.,
311 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, -2.0,2.0,0.,
312 0.,-2.0,0., 2.0,0.,0.,
314 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, -2.0,2.0,0.,
315 0.,-2.0,0., 2.0,0.,0.,
317 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, -2.0,2.0,0.,
318 0.,-2.0,0., 2.0,0.,0.,
321 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, -2.0,2.0,0.,
322 0.,-2.0,0., 2.0,0.,0.,
324 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, -2.0,2.0,0.,
325 0.,-2.0,0., 2.0,0.,0.,
327 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, -2.0,2.0,0.,
328 0.,-2.0,0., 2.0,0.,0.,
330 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, -2.0,2.0,0.,
331 0.,-2.0,0., 2.0,0.,0.,
333 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, -2.0,2.0,0.,
334 0.,-2.0,0., 2.0,0.,0.,
336 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, -2.0,2.0,0.,
337 0.,-2.0,0., 2.0,0.,0.,
344 int numPoints = tetNodes.dimension(0);
351 vals.
resize(numFields, numPoints, spaceDim);
352 tetBasis.
getValues(vals, tetNodes, OPERATOR_VALUE);
353 for (
int i = 0; i < numFields; i++) {
354 for (
int j = 0; j < numPoints; j++) {
355 for (
int k = 0; k < spaceDim; k++) {
358 int l = k + i * spaceDim + j * spaceDim * numFields;
359 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
361 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
364 *outStream <<
" At multi-index { ";
365 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
366 *outStream <<
"} computed value: " << vals(i,j,k)
367 <<
" but reference value: " << basisValues[l] <<
"\n";
374 vals.
resize(numFields, numPoints, spaceDim);
375 tetBasis.
getValues(vals, tetNodes, OPERATOR_CURL);
376 for (
int i = 0; i < numFields; i++) {
377 for (
int j = 0; j < numPoints; j++) {
378 for (
int k = 0; k < spaceDim; k++) {
381 int l = k + i * spaceDim + j * spaceDim * numFields;
382 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
384 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
387 *outStream <<
" At multi-index { ";
388 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
389 *outStream <<
"} computed curl component: " << vals(i,j,k)
390 <<
" but reference curl component: " << basisCurls[l] <<
"\n";
399 catch (
const std::logic_error & err) {
400 *outStream << err.what() <<
"\n\n";
406 <<
"===============================================================================\n"\
407 <<
"| TEST 4: correctness of DoF locations |\n"\
408 <<
"===============================================================================\n";
411 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
413 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
421 #ifdef HAVE_INTREPID_DEBUG
423 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
425 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
427 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
430 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
432 if (throwCounter != nException) {
434 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
439 tangents(0,0) = 1.0; tangents(0,1) = 0.0; tangents(0,2) = 0.0;
440 tangents(1,0) = -1.0; tangents(1,1) = 1.0; tangents(1,2) = 0.0;
441 tangents(2,0) = 0.0; tangents(2,1) = -1.0; tangents(2,2) = 0.0;
442 tangents(3,0) = 0.0; tangents(3,1) = 0.0; tangents(3,2) = 1.0;
443 tangents(4,0) = -1.0; tangents(4,1) = 0.0; tangents(4,2) = 1.0;
444 tangents(5,0) = 0.0; tangents(5,1) = -1.0; tangents(5,2) = 1.0;
446 basis->getValues(bvals, cvals, OPERATOR_VALUE);
448 for (
int i=0; i<bvals.dimension(0); i++) {
449 for (
int j=0; j<bvals.dimension(1); j++) {
451 double tangent = 0.0;
452 for(
int d=0;d<spaceDim;d++)
453 tangent += bvals(i,j,d)*tangents(j,d);
455 if ((i != j) && (std::abs(tangent - 0.0) > INTREPID_TOL)) {
457 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);
458 *outStream << buffer;
460 else if ((i == j) && (std::abs(tangent - 1.0) > INTREPID_TOL)) {
462 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);
463 *outStream << buffer;
469 catch (
const std::logic_error & err){
470 *outStream << err.what() <<
"\n\n";
475 std::cout <<
"End Result: TEST FAILED\n";
477 std::cout <<
"End Result: TEST PASSED\n";
480 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 Tetrahedron cell.
Implementation of the default H(curl)-compatible FEM basis of degree 1 on Tetrahedron 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_TET_I1_FEM class.