50#include "Teuchos_oblackholestream.hpp"
51#include "Teuchos_RCP.hpp"
52#include "Teuchos_GlobalMPISession.hpp"
55using 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"; \
71int 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_HGRAD_TET_C2_FEM) |\n" \
94 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 <<
"| 2) Basis values for VALUE, GRAD, and Dk 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 tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0;
119 tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0;
120 tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0;
121 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_CURL), throwCounter, nException );
143 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, tetNodes, OPERATOR_DIV), throwCounter, nException );
148 INTREPID_TEST_COMMAND( tetBasis.
getDofOrdinal(3,0,0), throwCounter, nException );
150 INTREPID_TEST_COMMAND( tetBasis.
getDofOrdinal(1,1,1), throwCounter, nException );
152 INTREPID_TEST_COMMAND( tetBasis.
getDofOrdinal(0,4,0), throwCounter, nException );
154 INTREPID_TEST_COMMAND( tetBasis.
getDofTag(10), throwCounter, nException );
156 INTREPID_TEST_COMMAND( tetBasis.
getDofTag(-1), throwCounter, nException );
158#ifdef HAVE_INTREPID_DEBUG
162 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
166 INTREPID_TEST_COMMAND( tetBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
170 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
174 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals2, tetNodes, OPERATOR_GRAD), throwCounter, nException );
177 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals2, tetNodes, OPERATOR_D1), throwCounter, nException );
180 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals2, tetNodes, OPERATOR_D2), throwCounter, nException );
184 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
188 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals4, tetNodes, OPERATOR_VALUE), throwCounter, nException );
192 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals5, tetNodes, OPERATOR_GRAD), throwCounter, nException );
196 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals6, tetNodes, OPERATOR_D1), throwCounter, nException );
199 INTREPID_TEST_COMMAND( tetBasis.
getValues(badVals6, tetNodes, OPERATOR_D2), throwCounter, nException );
203 catch (
const std::logic_error & err) {
204 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
205 *outStream << err.what() <<
'\n';
206 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
211 if (throwCounter != nException) {
213 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
218 <<
"===============================================================================\n"\
219 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
220 <<
"===============================================================================\n";
223 std::vector<std::vector<int> > allTags = tetBasis.
getAllDofTags();
226 for (
unsigned i = 0; i < allTags.size(); i++) {
227 int bfOrd = tetBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
229 std::vector<int> myTag = tetBasis.
getDofTag(bfOrd);
230 if( !( (myTag[0] == allTags[i][0]) &&
231 (myTag[1] == allTags[i][1]) &&
232 (myTag[2] == allTags[i][2]) &&
233 (myTag[3] == allTags[i][3]) ) ) {
235 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
236 *outStream <<
" getDofOrdinal( {"
237 << allTags[i][0] <<
", "
238 << allTags[i][1] <<
", "
239 << allTags[i][2] <<
", "
240 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
241 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
245 << myTag[3] <<
"}\n";
251 std::vector<int> myTag = tetBasis.
getDofTag(bfOrd);
252 int myBfOrd = tetBasis.
getDofOrdinal(myTag[0], myTag[1], myTag[2]);
253 if( bfOrd != myBfOrd) {
255 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
256 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
260 << myTag[3] <<
"} but getDofOrdinal({"
264 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
268 catch (
const std::logic_error & err){
269 *outStream << err.what() <<
"\n\n";
275 <<
"===============================================================================\n"\
276 <<
"| TEST 3: correctness of basis function values |\n"\
277 <<
"===============================================================================\n";
279 outStream -> precision(20);
282 double basisValues[] = {
283 1.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00000, 0, 0, 0, 0, 0, 0, 0, \
284 0, 0, 0, 1.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00000, 0, 0, 0, 0, \
285 0, 0, 0, 0, 0, 0, 1.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00000, 0, \
286 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
287 1.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00000, 0, 0, 0, 0, 0, 0, 0, \
291 double basisGrads[] = {
292 -3.00000, -3.00000, -3.00000, 1.00000, 1.00000, 1.00000, 1.00000, \
293 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, -1.00000, -1.00000, \
294 -1.00000, 1.00000, 1.00000, 1.00000, -1.00000, -1.00000, -1.00000, \
295 -1.00000, -1.00000, -1.00000, 1.00000, 1.00000, 1.00000, 1.00000, \
296 1.00000, 1.00000, -1.00000, 0, 0, 3.00000, 0, 0, -1.00000, 0, 0, \
297 -1.00000, 0, 0, 1.00000, 0, 0, 1.00000, 0, 0, -1.00000, 0, 0, \
298 -1.00000, 0, 0, 1.00000, 0, 0, -1.00000, 0, 0, 0, -1.00000, 0, 0, \
299 -1.00000, 0, 0, 3.00000, 0, 0, -1.00000, 0, 0, -1.00000, 0, 0, \
300 1.00000, 0, 0, 1.00000, 0, 0, -1.00000, 0, 0, -1.00000, 0, 0, \
301 1.00000, 0, 0, 0, -1.00000, 0, 0, -1.00000, 0, 0, -1.00000, 0, 0, \
302 3.00000, 0, 0, -1.00000, 0, 0, -1.00000, 0, 0, -1.00000, 0, 0, \
303 1.00000, 0, 0, 1.00000, 0, 0, 1.00000, 4.00000, 0, 0, -4.00000, \
304 -4.00000, -4.00000, 0, 0, 0, 0, 0, 0, 0, -2.00000, -2.00000, \
305 -2.00000, -2.00000, -2.00000, 2.00000, 0, 0, 2.00000, 0, 0, -2.00000, \
306 -2.00000, -2.00000, 0, 0, 0, 0, 0, 0, 0, 4.00000, 0, 4.00000, 0, 0, \
307 0, 0, 0, 0, 2.00000, 0, 2.00000, 2.00000, 0, 2.00000, 0, 0, 0, 0, 0, \
308 0, 2.00000, 0, 2.00000, 0, 0, 0, 4.00000, 0, 0, 0, 0, -4.00000, \
309 -4.00000, -4.00000, 0, 0, 0, 0, 2.00000, 0, -2.00000, -2.00000, \
310 -2.00000, -2.00000, 0, -2.00000, 0, 2.00000, 0, 0, 0, 0, -2.00000, \
311 -2.00000, -2.00000, 0, 0, 4.00000, 0, 0, 0, 0, 0, 0, -4.00000, \
312 -4.00000, -4.00000, 0, 0, 2.00000, 0, 0, 0, 0, 0, 2.00000, -2.00000, \
313 -2.00000, 0, -2.00000, -2.00000, -2.00000, -2.00000, -2.00000, \
314 -2.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 4.00000, 0, 0, 0, 0, \
315 2.00000, 0, 0, 2.00000, 0, 0, 0, 2.00000, 0, 0, 2.00000, 0, 2.00000, \
316 2.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.00000, 0, 4.00000, 0, 0, 0, \
317 0, 0, 0, 2.00000, 0, 0, 2.00000, 0, 2.00000, 0, 0, 2.00000, 0, 0, \
322 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
323 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
324 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
325 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
326 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
327 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
328 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
329 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
330 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 0, 0, 0, 0, 0, 4.00000, \
331 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, \
332 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, \
333 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, \
334 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, \
335 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, \
336 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, \
337 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 0, 0, 4.00000, \
338 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, \
339 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, \
340 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, \
341 0, 0, 4.00000, -8.00000, -4.00000, -4.00000, 0, 0, 0, -8.00000, \
342 -4.00000, -4.00000, 0, 0, 0, -8.00000, -4.00000, -4.00000, 0, 0, 0, \
343 -8.00000, -4.00000, -4.00000, 0, 0, 0, -8.00000, -4.00000, -4.00000, \
344 0, 0, 0, -8.00000, -4.00000, -4.00000, 0, 0, 0, -8.00000, -4.00000, \
345 -4.00000, 0, 0, 0, -8.00000, -4.00000, -4.00000, 0, 0, 0, -8.00000, \
346 -4.00000, -4.00000, 0, 0, 0, -8.00000, -4.00000, -4.00000, 0, 0, 0, \
347 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, \
348 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, \
349 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, \
350 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, -4.00000, 0, -8.00000, -4.00000, \
351 0, 0, -4.00000, 0, -8.00000, -4.00000, 0, 0, -4.00000, 0, -8.00000, \
352 -4.00000, 0, 0, -4.00000, 0, -8.00000, -4.00000, 0, 0, -4.00000, 0, \
353 -8.00000, -4.00000, 0, 0, -4.00000, 0, -8.00000, -4.00000, 0, 0, \
354 -4.00000, 0, -8.00000, -4.00000, 0, 0, -4.00000, 0, -8.00000, \
355 -4.00000, 0, 0, -4.00000, 0, -8.00000, -4.00000, 0, 0, -4.00000, 0, \
356 -8.00000, -4.00000, 0, 0, 0, -4.00000, 0, -4.00000, -8.00000, 0, 0, \
357 -4.00000, 0, -4.00000, -8.00000, 0, 0, -4.00000, 0, -4.00000, \
358 -8.00000, 0, 0, -4.00000, 0, -4.00000, -8.00000, 0, 0, -4.00000, 0, \
359 -4.00000, -8.00000, 0, 0, -4.00000, 0, -4.00000, -8.00000, 0, 0, \
360 -4.00000, 0, -4.00000, -8.00000, 0, 0, -4.00000, 0, -4.00000, \
361 -8.00000, 0, 0, -4.00000, 0, -4.00000, -8.00000, 0, 0, -4.00000, 0, \
362 -4.00000, -8.00000, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, \
363 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, \
364 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, \
365 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 0, 0, \
366 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, \
367 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, \
368 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, \
377 int numPoints = tetNodes.dimension(0);
384 vals.
resize(numFields, numPoints);
385 tetBasis.
getValues(vals, tetNodes, OPERATOR_VALUE);
386 for (
int i = 0; i < numFields; i++) {
387 for (
int j = 0; j < numPoints; j++) {
388 int l = i + j * numFields;
389 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
391 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
394 *outStream <<
" At multi-index { ";
395 *outStream << i <<
" ";*outStream << j <<
" ";
396 *outStream <<
"} computed value: " << vals(i,j)
397 <<
" but reference value: " << basisValues[l] <<
"\n";
403 vals.
resize(numFields, numPoints, spaceDim);
404 tetBasis.
getValues(vals, tetNodes, OPERATOR_GRAD);
405 for (
int i = 0; i < numFields; i++) {
406 for (
int j = 0; j < numPoints; j++) {
407 for (
int k = 0; k < spaceDim; k++) {
410 int l = k + j * spaceDim + i * spaceDim * numPoints;
411 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
413 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
416 *outStream <<
" At multi-index { ";
417 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
418 *outStream <<
"} computed grad component: " << vals(i,j,k)
419 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
426 tetBasis.
getValues(vals, tetNodes, OPERATOR_D1);
427 for (
int i = 0; i < numFields; i++) {
428 for (
int j = 0; j < numPoints; j++) {
429 for (
int k = 0; k < spaceDim; k++) {
432 int l = k + j * spaceDim + i * spaceDim * numPoints;
433 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
435 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
438 *outStream <<
" At multi-index { ";
439 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
440 *outStream <<
"} computed D1 component: " << vals(i,j,k)
441 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
449 vals.
resize(numFields, numPoints, D2cardinality);
450 tetBasis.
getValues(vals, tetNodes, OPERATOR_D2);
451 for (
int i = 0; i < numFields; i++) {
452 for (
int j = 0; j < numPoints; j++) {
453 for (
int k = 0; k < D2cardinality; k++) {
456 int l = k + j * D2cardinality + i * D2cardinality * numPoints;
457 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
459 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
462 *outStream <<
" At multi-index { ";
463 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
464 *outStream <<
"} computed D2 component: " << vals(i,j,k)
465 <<
" but reference D2 component: " << basisD2[l] <<
"\n";
473 for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
477 vals.
resize(numFields, numPoints, DkCardin);
480 for (
int i = 0; i < vals.
size(); i++) {
481 if (std::abs(vals[i]) > INTREPID_TOL) {
483 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
486 std::vector<int> myIndex;
489 *outStream <<
" At multi-index { ";
490 for(
int j = 0; j < vals.
rank(); j++) {
491 *outStream << myIndex[j] <<
" ";
493 *outStream <<
"} computed D"<< ord <<
" component: " << vals[i]
494 <<
" but reference D" << ord <<
" component: 0 \n";
501 catch (
const std::logic_error & err) {
502 *outStream << err.what() <<
"\n\n";
507 std::cout <<
"End Result: TEST FAILED\n";
509 std::cout <<
"End Result: TEST PASSED\n";
512 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_TET_C2_FEM class.
int getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Tetrahedron cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Tetrahedron cell.
virtual int getDofOrdinal(const int subcDim, const int subcOrd, const int subcDofOrd)
DoF tag to ordinal lookup.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
virtual int getCardinality() const
Returns cardinality of the basis.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
void getMultiIndex(int &i0, const int valueEnum) const
Returns the multi-index of a value, based on its enumeration, as a list, for rank-1 containers.
int rank() const
Return rank of the FieldContainer = number of indices used to tag the multi-indexed value.