49#include "Intrepid_HGRAD_PYR_C1_FEM.hpp"
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_PYR_C1_FEM) |\n" \
94 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 <<
"| 2) Basis values for VALUE, GRAD, CURL, 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_PYR_C1_FEM<double, FieldContainer<double> > pyrBasis;
114 int throwCounter = 0;
117 FieldContainer<double> pyrNodes(10, 3);
118 pyrNodes(0,0) = -1.0; pyrNodes(0,1) = -1.0; pyrNodes(0,2) = 0;
119 pyrNodes(1,0) = 1.0; pyrNodes(1,1) = -1.0; pyrNodes(1,2) = 0;
120 pyrNodes(2,0) = 1.0; pyrNodes(2,1) = 1.0; pyrNodes(2,2) = 0;
121 pyrNodes(3,0) = -1.0; pyrNodes(3,1) = 1.0; pyrNodes(3,2) = 0;
122 pyrNodes(4,0) = 0.0; pyrNodes(4,1) = 0.0; pyrNodes(4,2) = 1.0;
124 pyrNodes(5,0) = 0.25; pyrNodes(5,1) = 0.5; pyrNodes(5,2) = 0.2;
125 pyrNodes(6,0) = -0.7 ; pyrNodes(6,1) = 0.3; pyrNodes(6,2) = 0.3;
126 pyrNodes(7,0) = 0.; pyrNodes(7,1) = -0.05; pyrNodes(7,2) = 0.95;
127 pyrNodes(8,0) = -0.15; pyrNodes(8,1) = -0.2; pyrNodes(8,2) = 0.75;
128 pyrNodes(9,0) = -0.4; pyrNodes(9,1) = 0.9; pyrNodes(9,2) = 0.0;
132 FieldContainer<double> vals;
137 vals.resize(pyrBasis.getCardinality(), pyrNodes.dimension(0), 3 );
138 INTREPID_TEST_COMMAND( pyrBasis.getValues(vals, pyrNodes, OPERATOR_CURL), throwCounter, nException );
142 vals.resize(pyrBasis.getCardinality(), pyrNodes.dimension(0) );
143 INTREPID_TEST_COMMAND( pyrBasis.getValues(vals, pyrNodes, OPERATOR_DIV), throwCounter, nException );
148 INTREPID_TEST_COMMAND( pyrBasis.getDofOrdinal(3,0,0), throwCounter, nException );
150 INTREPID_TEST_COMMAND( pyrBasis.getDofOrdinal(1,1,1), throwCounter, nException );
152 INTREPID_TEST_COMMAND( pyrBasis.getDofOrdinal(0,6,0), throwCounter, nException );
154 INTREPID_TEST_COMMAND( pyrBasis.getDofTag(7), throwCounter, nException );
156 INTREPID_TEST_COMMAND( pyrBasis.getDofTag(-1), throwCounter, nException );
158#ifdef HAVE_INTREPID_DEBUG
161 FieldContainer<double> badPoints1(4, 5, 3);
162 INTREPID_TEST_COMMAND( pyrBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
165 FieldContainer<double> badPoints2(4, 2);
166 INTREPID_TEST_COMMAND( pyrBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
169 FieldContainer<double> badVals1(4, 3, 1);
170 INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals1, pyrNodes, OPERATOR_VALUE), throwCounter, nException );
173 FieldContainer<double> badVals2(4, 3);
174 INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals2, pyrNodes, OPERATOR_GRAD), throwCounter, nException );
177 INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals2, pyrNodes, OPERATOR_D1), throwCounter, nException );
183 FieldContainer<double> badVals3(pyrBasis.getCardinality() + 1, pyrNodes.dimension(0));
184 INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals3, pyrNodes, OPERATOR_VALUE), throwCounter, nException );
187 FieldContainer<double> badVals4(pyrBasis.getCardinality(), pyrNodes.dimension(0) + 1);
188 INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals4, pyrNodes, OPERATOR_VALUE), throwCounter, nException );
191 FieldContainer<double> badVals5(pyrBasis.getCardinality(), pyrNodes.dimension(0), 4);
192 INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals5, pyrNodes, OPERATOR_GRAD), 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 = pyrBasis.getAllDofTags();
226 for (
unsigned i = 0; i < allTags.size(); i++) {
227 int bfOrd = pyrBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
229 std::vector<int> myTag = pyrBasis.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 < pyrBasis.getCardinality(); bfOrd++) {
251 std::vector<int> myTag = pyrBasis.getDofTag(bfOrd);
252 int myBfOrd = pyrBasis.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.0, 0.0, 0.0, 0.0, 0.0,
284 0.0, 1.0, 0.0, 0.0, 0.0,
285 0.0, 0.0, 1.0, 0.0, 0.0,
286 0.0, 0.0, 0.0, 1.0, 0.0,
287 0.0, 0.0, 0.0, 0.0, 1.0,
289 0.0515625, 0.0984375, 0.4265625, 0.2234375, 0.2,
290 0.2, 0, 0., 0.5, 0.3,
291 0.025, 0.025, 0., 0, 0.95,
292 0.18, 0.045, 0.005, 0.02, 0.75,
293 0.035, 0.015, 0.285, 0.665, 0.,
297 double basisGrads[] = {
298 -0.5, -0.5, 0.0, 0.5, 0.0, -0.5, 0.0, 0.0, 0.0, 0.0, 0.5, -0.5, 0.0, 0.0, 1.0, \
299 -0.5, 0.0, -0.5, 0.5, -0.5, 0.0, 0.0, 0.5, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, \
300 0.0, 0.0, 0.0, 0.0, -0.5, -0.5, 0.5, 0.5, 0.0, -0.5, 0.0, -0.5, 0.0, 0.0, 1.0, \
301 0.0, -0.5, -0.5, 0.0, 0.0, 0.0, 0.5, 0.0, -0.5, -0.5, 0.5, 0.0, 0.0, 0.0, 1.0, \
302 -0.25,-0.25,-0.25, 0.25,-0.25,-0.25, 0.25, 0.25,-0.25,-0.25, 0.25,-0.25, 0.0, 0.0, 1.0, \
303 -0.09375, -0.171875, -0.201171875, 0.09375, -0.328125, -0.298828125, 0.40625, 0.328125, -0.201171875, -0.40625, 0.171875, -0.298828125, 0.0, 0.0, 1.0, \
304 -0.1428571428571429, -0.5, -0.3571428571428571, 0.1428571428571429, 0.0, -0.1428571428571429, 0.3571428571428571, 0.0, -0.3571428571428571, -0.3571428571428571, 0.5, -0.1428571428571429, 0.0, 0.0, 1.0, \
305 -0.5, -0.25, -0.25, 0.5, -0.25, -0.25, 0.0, 0.25, -0.25, 0., 0.25, -0.25, 0.0, 0.0, 1.0, \
306 -0.45, -0.4, -0.13, 0.45, -0.1, -0.37, 0.05, 0.1, -0.13, -0.05, 0.4, -0.37, 0.0, 0.0, 1.0, \
307 -0.025, -0.35, -0.34, 0.025, -0.15, -0.16, 0.475, 0.15, -0.34, -0.475, 0.35, -0.16, 0.0, 0.0, 1.0
312 const double eps = std::numeric_limits<double>::epsilon( );
314 0, 0.25,-0.25, 0,-0.25, 0.5, 0,-0.25, 0.25, 0, 0.25,-0.5, 0, 0.25,-0.25, 0,-0.25, 0.5, 0,-0.25, 0.25, 0, 0.25,-0.5, 0, 0, 0, 0, 0, 0, \
315 0, 0.25,-0.25, 0, 0.25,-0.5, 0,-0.25, 0.25, 0,-0.25, 0.5, 0, 0.25,-0.25, 0, 0.25,-0.5, 0,-0.25, 0.25, 0,-0.25, 0.5, 0, 0, 0, 0, 0, 0, \
316 0, 0.25, 0.25, 0, 0.25, 0.5, 0,-0.25,-0.25, 0,-0.25,-0.5, 0, 0.25, 0.25, 0, 0.25, 0.5, 0,-0.25,-0.25, 0,-0.25,-0.5, 0, 0, 0, 0, 0, 0, \
317 0, 0.25, 0.25, 0,-0.25,-0.5, 0,-0.25,-0.25, 0, 0.25, 0.5, 0, 0.25, 0.25, 0,-0.25,-0.5, 0,-0.25,-0.25, 0, 0.25, 0.5, 0, 0, 0, 0, 0, 0, \
318 0, 0.25/eps, 0, 0, 0, 0, 0,-0.25/eps, 0, 0, 0, 0, 0, 0.25/eps, 0, 0, 0, 0, 0,-0.25/eps, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, \
319 0, 0.3125, 0.1953125, 0, 0.09765625, 0.1220703125, 0,-0.3125,-0.1953125, 0,-0.09765625,-0.1220703125, 0, 0.3125, 0.1953125, 0, 0.09765625, 0.1220703125, 0,-0.3125,-0.1953125, 0,-0.09765625,-0.1220703125, 0, 0, 0, 0, 0, 0, \
320 0, 0.3571428571428571, 0.1530612244897959, 0,-0.3571428571428571,-0.306122448979592, 0,-0.3571428571428572,-0.1530612244897959, 0, 0.3571428571428571, 0.306122448979592, 0, 0.3571428571428571, 0.1530612244897959, 0,-0.3571428571428571,-0.306122448979592, 0,-0.3571428571428571,-0.1530612244897959, 0, 0.3571428571428571, 0.306122448979592, 0, 0, 0, 0, 0, 0, \
321 0, 5,-5, 0, 0, 0, 0,-5, 5, 0, 0, 0, 0, 5,-5, 0, 0, 0, 0, -5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
322 0, 1,-0.8, 0,-0.6, 0.96, 0, -1, 0.8, 0, 0.6,-0.96, 0,1,-0.8, 0, -0.6, 0.96, 0, -1, 0.8, 0, 0.6, -0.96, 0, 0, 0, 0, 0, 0, \
323 0, 0.25, 0.225, 0,-0.1,-0.18, 0,-0.25,-0.225, 0, 0.1, 0.18, 0, 0.25, 0.225, 0,-0.1,-0.18, 0,-0.25,-0.225,0,0.1,0.18, 0, 0, 0, 0, 0, 0
328 int numFields = pyrBasis.getCardinality();
329 int numPoints = pyrNodes.dimension(0);
330 int spaceDim = pyrBasis.getBaseCellTopology().getDimension();
334 FieldContainer<double> vals;
337 vals.resize(numFields, numPoints);
338 pyrBasis.getValues(vals, pyrNodes, OPERATOR_VALUE);
339 for (
int i = 0; i < numFields; i++) {
340 for (
int j = 0; j < numPoints; j++) {
341 int l = i + j * numFields;
342 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
344 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
347 *outStream <<
" At multi-index { ";
348 *outStream << i <<
" ";*outStream << j <<
" ";
349 *outStream <<
"} computed value: " << vals(i,j)
350 <<
" but reference value: " << basisValues[l] <<
"\n";
356 vals.resize(numFields, numPoints, spaceDim);
357 pyrBasis.getValues(vals, pyrNodes, OPERATOR_GRAD);
358 for (
int i = 0; i < numFields; i++) {
359 for (
int j = 0; j < numPoints; j++) {
360 for (
int k = 0; k < spaceDim; k++) {
361 int l = k + i * spaceDim + j * spaceDim * numFields;
362 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
364 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
367 *outStream <<
" At multi-index { ";
368 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
369 *outStream <<
"} computed grad component: " << vals(i,j,k)
370 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
377 pyrBasis.getValues(vals, pyrNodes, OPERATOR_D1);
378 for (
int i = 0; i < numFields; i++) {
379 for (
int j = 0; j < numPoints; j++) {
380 for (
int k = 0; k < spaceDim; k++) {
381 int l = k + i * spaceDim + j * spaceDim * numFields;
382 if (std::abs(vals(i,j,k) - basisGrads[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 D1 component: " << vals(i,j,k)
390 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
397 vals.resize(numFields, numPoints, D2Cardin);
398 pyrBasis.getValues(vals, pyrNodes, OPERATOR_D2);
399 for (
int i = 0; i < numFields; i++) {
400 for (
int j = 0; j < numPoints; j++) {
402 for (
int k = 0; k < D2Cardin; k++) {
403 int l = k + i * D2Cardin + j * D2Cardin * numFields;
404 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
406 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
409 *outStream <<
" At multi-index { ";
410 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
411 *outStream <<
"} computed D2 component: " << vals(i,j,k)
412 <<
" but reference D2 component: " << basisD2[l] <<
"\n";
421 catch (
const std::logic_error & err) {
422 *outStream << err.what() <<
"\n\n";
427 std::cout <<
"End Result: TEST FAILED\n";
429 std::cout <<
"End Result: TEST PASSED\n";
432 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.