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";
109 Basis_HGRAD_TET_C2_FEM<double, FieldContainer<double> > tetBasis;
114 int throwCounter = 0;
117 FieldContainer<double> tetNodes(10, 3);
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;
132 FieldContainer<double> vals;
137 vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0), 4 );
138 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException );
142 vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0) );
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
161 FieldContainer<double> badPoints1(4, 5, 3);
162 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
165 FieldContainer<double> badPoints2(4, tetBasis.getBaseCellTopology().getDimension() - 1);
166 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
169 FieldContainer<double> badVals1(4, 3, 1);
170 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
173 FieldContainer<double> badVals2(4, 3);
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 );
183 FieldContainer<double> badVals3(tetBasis.getCardinality() + 1, tetNodes.dimension(0));
184 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
187 FieldContainer<double> badVals4(tetBasis.getCardinality(), tetNodes.dimension(0) + 1);
188 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals4, tetNodes, OPERATOR_VALUE), throwCounter, nException );
191 FieldContainer<double> badVals5(tetBasis.getCardinality(), tetNodes.dimension(0), tetBasis.getBaseCellTopology().getDimension() + 1);
192 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals5, tetNodes, OPERATOR_GRAD), throwCounter, nException );
195 FieldContainer<double> badVals6(tetBasis.getCardinality(), tetNodes.dimension(0), 40);
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";
250 for(
int bfOrd = 0; bfOrd < tetBasis.getCardinality(); bfOrd++) {
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, \
376 int numFields = tetBasis.getCardinality();
377 int numPoints = tetNodes.dimension(0);
378 int spaceDim = tetBasis.getBaseCellTopology().getDimension();
381 FieldContainer<double> vals;
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);
479 tetBasis.getValues(vals, tetNodes, op);
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;
487 vals.getMultiIndex(myIndex,i);
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.