49#include "Intrepid_HGRAD_QUAD_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_QUAD_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_QUAD_C1_FEM<double, FieldContainer<double> > quadBasis;
114 int throwCounter = 0;
117 FieldContainer<double> quadNodes(5, 2);
118 quadNodes(0,0) = -1.0; quadNodes(0,1) = -1.0;
119 quadNodes(1,0) = 1.0; quadNodes(1,1) = -1.0;
120 quadNodes(2,0) = 1.0; quadNodes(2,1) = 1.0;
121 quadNodes(3,0) = -1.0; quadNodes(3,1) = 1.0;
122 quadNodes(4,0) = 0.0; quadNodes(4,1) = 0.0;
127 FieldContainer<double> vals;
131 vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0));
132 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException );
137 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException );
139 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException );
141 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,0), throwCounter, nException );
143 INTREPID_TEST_COMMAND( quadBasis.getDofTag(5), throwCounter, nException );
145 INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException );
147#ifdef HAVE_INTREPID_DEBUG
150 FieldContainer<double> badPoints1(4, 5, 3);
151 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
154 FieldContainer<double> badPoints2(4, 3);
155 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
158 FieldContainer<double> badVals1(4, 3, 1);
159 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
162 FieldContainer<double> badVals2(4, 3);
163 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_GRAD), throwCounter, nException );
166 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_CURL), throwCounter, nException );
169 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_D2), throwCounter, nException );
172 FieldContainer<double> badVals3(quadBasis.getCardinality() + 1, quadNodes.dimension(0));
173 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException );
176 FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0) + 1);
177 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException );
180 FieldContainer<double> badVals5(quadBasis.getCardinality(), quadNodes.dimension(0), 4);
181 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals5, quadNodes, OPERATOR_GRAD), throwCounter, nException );
184 FieldContainer<double> badVals6(quadBasis.getCardinality(), quadNodes.dimension(0), 40);
185 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D2), throwCounter, nException );
188 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D3), throwCounter, nException );
192 catch (
const std::logic_error & err) {
193 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
194 *outStream << err.what() <<
'\n';
195 *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 = quadBasis.getAllDofTags();
215 for (
unsigned i = 0; i < allTags.size(); i++) {
216 int bfOrd = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
218 std::vector<int> myTag = quadBasis.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";
239 for(
int bfOrd = 0; bfOrd < quadBasis.getCardinality(); bfOrd++) {
240 std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
241 int myBfOrd = quadBasis.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[] = {
280 double basisGrads[] = {
281 -0.5, -0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5,
282 -0.5, 0.0, 0.5, -0.5, 0.0, 0.5, 0.0, 0.0,
283 0.0, 0.0, 0.0, -0.5, 0.5, 0.5, -0.5, 0.0,
284 0.0, -0.5, 0.0, 0.0, 0.5, 0.0, -0.5, 0.5,
285 -0.25,-0.25, 0.25,-0.25, 0.25, 0.25, -0.25, 0.25
289 double basisCurls[] = {
290 -0.5, 0.5, 0.0, -0.5, 0.0, 0.0, 0.5, 0.0,
291 0.0, 0.5, -0.5, -0.5, 0.5, 0.0, 0.0, 0.0,
292 0.0, 0.0, -0.5, 0.0, 0.5, -0.5, 0.0, 0.5,
293 -0.5, 0.0, 0.0, 0.0, 0.0, -0.5, 0.5, 0.5,
294 -0.25, 0.25, -0.25,-0.25, 0.25,-0.25, 0.25, 0.25
299 0.0, 0.25, 0.0, 0.0,-0.25, 0.0, 0.0, 0.25, 0.0, 0.0,-0.25, 0.0,
300 0.0, 0.25, 0.0, 0.0,-0.25, 0.0, 0.0, 0.25, 0.0, 0.0,-0.25, 0.0,
301 0.0, 0.25, 0.0, 0.0,-0.25, 0.0, 0.0, 0.25, 0.0, 0.0,-0.25, 0.0,
302 0.0, 0.25, 0.0, 0.0,-0.25, 0.0, 0.0, 0.25, 0.0, 0.0,-0.25, 0.0,
303 0.0, 0.25, 0.0, 0.0,-0.25, 0.0, 0.0, 0.25, 0.0, 0.0,-0.25, 0.0,
309 int numFields = quadBasis.getCardinality();
310 int numPoints = quadNodes.dimension(0);
311 int spaceDim = quadBasis.getBaseCellTopology().getDimension();
315 FieldContainer<double> vals;
318 vals.resize(numFields, numPoints);
319 quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE);
320 for (
int i = 0; i < numFields; i++) {
321 for (
int j = 0; j < numPoints; j++) {
322 int l = i + j * numFields;
323 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
325 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
328 *outStream <<
" At multi-index { ";
329 *outStream << i <<
" ";*outStream << j <<
" ";
330 *outStream <<
"} computed value: " << vals(i,j)
331 <<
" but reference value: " << basisValues[l] <<
"\n";
337 vals.resize(numFields, numPoints, spaceDim);
338 quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD);
339 for (
int i = 0; i < numFields; i++) {
340 for (
int j = 0; j < numPoints; j++) {
341 for (
int k = 0; k < spaceDim; k++) {
342 int l = k + i * spaceDim + j * spaceDim * numFields;
343 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
345 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
348 *outStream <<
" At multi-index { ";
349 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
350 *outStream <<
"} computed grad component: " << vals(i,j,k)
351 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
359 quadBasis.getValues(vals, quadNodes, OPERATOR_D1);
360 for (
int i = 0; i < numFields; i++) {
361 for (
int j = 0; j < numPoints; j++) {
362 for (
int k = 0; k < spaceDim; k++) {
363 int l = k + i * spaceDim + j * spaceDim * numFields;
364 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
366 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
369 *outStream <<
" At multi-index { ";
370 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
371 *outStream <<
"} computed D1 component: " << vals(i,j,k)
372 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
380 vals.resize(numFields, numPoints, spaceDim);
381 quadBasis.getValues(vals, quadNodes, OPERATOR_CURL);
382 for (
int i = 0; i < numFields; i++) {
383 for (
int j = 0; j < numPoints; j++) {
384 for (
int k = 0; k < spaceDim; k++) {
385 int l = k + i * spaceDim + j * spaceDim * numFields;
386 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
388 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
391 *outStream <<
" At multi-index { ";
392 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
393 *outStream <<
"} computed curl component: " << vals(i,j,k)
394 <<
" but reference curl component: " << basisCurls[l] <<
"\n";
401 vals.resize(numFields, numPoints, D2Cardin);
402 quadBasis.getValues(vals, quadNodes, OPERATOR_D2);
403 for (
int i = 0; i < numFields; i++) {
404 for (
int j = 0; j < numPoints; j++) {
405 for (
int k = 0; k < D2Cardin; k++) {
406 int l = k + i * D2Cardin + j * D2Cardin * numFields;
407 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
409 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
412 *outStream <<
" At multi-index { ";
413 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
414 *outStream <<
"} computed D2 component: " << vals(i,j,k)
415 <<
" but reference D2 component: " << basisD2[l] <<
"\n";
423 for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
427 vals.resize(numFields, numPoints, DkCardin);
429 quadBasis.getValues(vals, quadNodes, op);
430 for (
int i = 0; i < vals.size(); i++) {
431 if (std::abs(vals[i]) > INTREPID_TOL) {
433 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
436 std::vector<int> myIndex;
437 vals.getMultiIndex(myIndex,i);
439 *outStream <<
" At multi-index { ";
440 for(
int j = 0; j < vals.rank(); j++) {
441 *outStream << myIndex[j] <<
" ";
443 *outStream <<
"} computed D"<< ord <<
" component: " << vals[i]
444 <<
" but reference D" << ord <<
" component: 0 \n";
450 catch (
const std::logic_error & err) {
451 *outStream << err.what() <<
"\n\n";
458 <<
"===============================================================================\n"\
459 <<
"| TEST 4: correctness of DoF locations |\n"\
460 <<
"===============================================================================\n";
463 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
465 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
466 Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
468 FieldContainer<double> cvals;
469 FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality());
472#ifdef HAVE_INTREPID_DEBUG
474 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
476 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
478 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
481 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
483 if (throwCounter != nException) {
485 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
489 basis->getValues(bvals, cvals, OPERATOR_VALUE);
491 for (
int i=0; i<bvals.dimension(0); i++) {
492 for (
int j=0; j<bvals.dimension(1); j++) {
493 if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) {
495 sprintf(buffer,
"\nValue of basis function %d at (%6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), bvals(i,j), 0.0);
496 *outStream << buffer;
498 else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) {
500 sprintf(buffer,
"\nValue of basis function %d at (%6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), bvals(i,j), 1.0);
501 *outStream << buffer;
507 catch (
const std::logic_error & err){
508 *outStream << err.what() <<
"\n\n";
513 std::cout <<
"End Result: TEST FAILED\n";
515 std::cout <<
"End Result: TEST PASSED\n";
518 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.
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 1 on Quadrilateral cell.