49 #include "Intrepid_HGRAD_PYR_C1_FEM.hpp"
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_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";
114 int throwCounter = 0;
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;
138 INTREPID_TEST_COMMAND( pyrBasis.
getValues(vals, pyrNodes, OPERATOR_CURL), throwCounter, nException );
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
162 INTREPID_TEST_COMMAND( pyrBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
166 INTREPID_TEST_COMMAND( pyrBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
170 INTREPID_TEST_COMMAND( pyrBasis.
getValues(badVals1, pyrNodes, OPERATOR_VALUE), throwCounter, nException );
174 INTREPID_TEST_COMMAND( pyrBasis.
getValues(badVals2, pyrNodes, OPERATOR_GRAD), throwCounter, nException );
177 INTREPID_TEST_COMMAND( pyrBasis.
getValues(badVals2, pyrNodes, OPERATOR_D1), throwCounter, nException );
184 INTREPID_TEST_COMMAND( pyrBasis.
getValues(badVals3, pyrNodes, OPERATOR_VALUE), throwCounter, nException );
188 INTREPID_TEST_COMMAND( pyrBasis.
getValues(badVals4, pyrNodes, OPERATOR_VALUE), throwCounter, nException );
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";
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
329 int numPoints = pyrNodes.dimension(0);
331 int D2Cardin = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
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);
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
virtual int getCardinality() const
Returns cardinality of the basis.
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
FEM basis evaluation on a reference Pyramid cell.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Pyramid cell...
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...