49 #include "Intrepid_HGRAD_QUAD_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_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";
114 int throwCounter = 0;
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;
131 INTREPID_TEST_COMMAND( quadBasis.
getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException );
136 INTREPID_TEST_COMMAND( quadBasis.
getDofOrdinal(3,0,0), throwCounter, nException );
138 INTREPID_TEST_COMMAND( quadBasis.
getDofOrdinal(1,1,1), throwCounter, nException );
140 INTREPID_TEST_COMMAND( quadBasis.
getDofOrdinal(0,4,0), throwCounter, nException );
142 INTREPID_TEST_COMMAND( quadBasis.
getDofTag(5), throwCounter, nException );
144 INTREPID_TEST_COMMAND( quadBasis.
getDofTag(-1), throwCounter, nException );
146 #ifdef HAVE_INTREPID_DEBUG
150 INTREPID_TEST_COMMAND( quadBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
154 INTREPID_TEST_COMMAND( quadBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
158 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
162 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals2, quadNodes, OPERATOR_GRAD), throwCounter, nException );
165 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals2, quadNodes, OPERATOR_CURL), throwCounter, nException );
168 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals2, quadNodes, OPERATOR_D2), throwCounter, nException );
172 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException );
176 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException );
180 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals5, quadNodes, OPERATOR_GRAD), throwCounter, nException );
184 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals6, quadNodes, OPERATOR_D2), throwCounter, nException );
187 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals6, quadNodes, OPERATOR_D3), throwCounter, nException );
191 catch (
const std::logic_error & err) {
192 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
193 *outStream << err.what() <<
'\n';
194 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
199 if (throwCounter != nException) {
201 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
206 <<
"===============================================================================\n"\
207 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
208 <<
"===============================================================================\n";
211 std::vector<std::vector<int> > allTags = quadBasis.
getAllDofTags();
214 for (
unsigned i = 0; i < allTags.size(); i++) {
215 int bfOrd = quadBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
217 std::vector<int> myTag = quadBasis.
getDofTag(bfOrd);
218 if( !( (myTag[0] == allTags[i][0]) &&
219 (myTag[1] == allTags[i][1]) &&
220 (myTag[2] == allTags[i][2]) &&
221 (myTag[3] == allTags[i][3]) ) ) {
223 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
224 *outStream <<
" getDofOrdinal( {"
225 << allTags[i][0] <<
", "
226 << allTags[i][1] <<
", "
227 << allTags[i][2] <<
", "
228 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
229 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
233 << myTag[3] <<
"}\n";
238 for(
int bfOrd = 0; bfOrd < quadBasis.
getCardinality(); bfOrd++) {
239 std::vector<int> myTag = quadBasis.
getDofTag(bfOrd);
240 int myBfOrd = quadBasis.
getDofOrdinal(myTag[0], myTag[1], myTag[2]);
241 if( bfOrd != myBfOrd) {
243 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
244 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
248 << myTag[3] <<
"} but getDofOrdinal({"
252 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
256 catch (
const std::logic_error & err){
257 *outStream << err.what() <<
"\n\n";
263 <<
"===============================================================================\n"\
264 <<
"| TEST 3: correctness of basis function values |\n"\
265 <<
"===============================================================================\n";
267 outStream -> precision(20);
270 double basisValues[] = {
279 double basisGrads[] = {
280 -0.5, -0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5,
281 -0.5, 0.0, 0.5, -0.5, 0.0, 0.5, 0.0, 0.0,
282 0.0, 0.0, 0.0, -0.5, 0.5, 0.5, -0.5, 0.0,
283 0.0, -0.5, 0.0, 0.0, 0.5, 0.0, -0.5, 0.5,
284 -0.25,-0.25, 0.25,-0.25, 0.25, 0.25, -0.25, 0.25
288 double basisCurls[] = {
289 -0.5, 0.5, 0.0, -0.5, 0.0, 0.0, 0.5, 0.0,
290 0.0, 0.5, -0.5, -0.5, 0.5, 0.0, 0.0, 0.0,
291 0.0, 0.0, -0.5, 0.0, 0.5, -0.5, 0.0, 0.5,
292 -0.5, 0.0, 0.0, 0.0, 0.0, -0.5, 0.5, 0.5,
293 -0.25, 0.25, -0.25,-0.25, 0.25,-0.25, 0.25, 0.25
298 0.0, 0.25, 0.0, 0.0,-0.25, 0.0, 0.0, 0.25, 0.0, 0.0,-0.25, 0.0,
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,
309 int numPoints = quadNodes.dimension(0);
311 int D2Cardin = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
317 vals.
resize(numFields, numPoints);
318 quadBasis.
getValues(vals, quadNodes, OPERATOR_VALUE);
319 for (
int i = 0; i < numFields; i++) {
320 for (
int j = 0; j < numPoints; j++) {
321 int l = i + j * numFields;
322 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
324 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
327 *outStream <<
" At multi-index { ";
328 *outStream << i <<
" ";*outStream << j <<
" ";
329 *outStream <<
"} computed value: " << vals(i,j)
330 <<
" but reference value: " << basisValues[l] <<
"\n";
336 vals.
resize(numFields, numPoints, spaceDim);
337 quadBasis.
getValues(vals, quadNodes, OPERATOR_GRAD);
338 for (
int i = 0; i < numFields; i++) {
339 for (
int j = 0; j < numPoints; j++) {
340 for (
int k = 0; k < spaceDim; k++) {
341 int l = k + i * spaceDim + j * spaceDim * numFields;
342 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
344 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
347 *outStream <<
" At multi-index { ";
348 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
349 *outStream <<
"} computed grad component: " << vals(i,j,k)
350 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
358 quadBasis.
getValues(vals, quadNodes, OPERATOR_D1);
359 for (
int i = 0; i < numFields; i++) {
360 for (
int j = 0; j < numPoints; j++) {
361 for (
int k = 0; k < spaceDim; k++) {
362 int l = k + i * spaceDim + j * spaceDim * numFields;
363 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
365 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
368 *outStream <<
" At multi-index { ";
369 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
370 *outStream <<
"} computed D1 component: " << vals(i,j,k)
371 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
379 vals.
resize(numFields, numPoints, spaceDim);
380 quadBasis.
getValues(vals, quadNodes, OPERATOR_CURL);
381 for (
int i = 0; i < numFields; i++) {
382 for (
int j = 0; j < numPoints; j++) {
383 for (
int k = 0; k < spaceDim; k++) {
384 int l = k + i * spaceDim + j * spaceDim * numFields;
385 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
387 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
390 *outStream <<
" At multi-index { ";
391 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
392 *outStream <<
"} computed curl component: " << vals(i,j,k)
393 <<
" but reference curl component: " << basisCurls[l] <<
"\n";
400 vals.
resize(numFields, numPoints, D2Cardin);
401 quadBasis.
getValues(vals, quadNodes, OPERATOR_D2);
402 for (
int i = 0; i < numFields; i++) {
403 for (
int j = 0; j < numPoints; j++) {
404 for (
int k = 0; k < D2Cardin; k++) {
405 int l = k + i * D2Cardin + j * D2Cardin * numFields;
406 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
408 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
411 *outStream <<
" At multi-index { ";
412 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
413 *outStream <<
"} computed D2 component: " << vals(i,j,k)
414 <<
" but reference D2 component: " << basisD2[l] <<
"\n";
422 for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
425 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
426 vals.
resize(numFields, numPoints, DkCardin);
428 quadBasis.
getValues(vals, quadNodes, op);
429 for (
int i = 0; i < vals.
size(); i++) {
430 if (std::abs(vals[i]) > INTREPID_TOL) {
432 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
435 std::vector<int> myIndex;
437 int ord = Intrepid::getOperatorOrder(op);
438 *outStream <<
" At multi-index { ";
439 for(
int j = 0; j < vals.
rank(); j++) {
440 *outStream << myIndex[j] <<
" ";
442 *outStream <<
"} computed D"<< ord <<
" component: " << vals[i]
443 <<
" but reference D" << ord <<
" component: 0 \n";
449 catch (
const std::logic_error & err) {
450 *outStream << err.what() <<
"\n\n";
457 <<
"===============================================================================\n"\
458 <<
"| TEST 4: correctness of DoF locations |\n"\
459 <<
"===============================================================================\n";
462 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
464 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
471 #ifdef HAVE_INTREPID_DEBUG
473 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
475 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
477 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
480 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
482 if (throwCounter != nException) {
484 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
488 basis->getValues(bvals, cvals, OPERATOR_VALUE);
490 for (
int i=0; i<bvals.dimension(0); i++) {
491 for (
int j=0; j<bvals.dimension(1); j++) {
492 if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) {
494 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);
495 *outStream << buffer;
497 else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) {
499 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);
500 *outStream << buffer;
506 catch (
const std::logic_error & err){
507 *outStream << err.what() <<
"\n\n";
512 std::cout <<
"End Result: TEST FAILED\n";
514 std::cout <<
"End Result: TEST PASSED\n";
517 std::cout.copyfmt(oldFormatState);
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
virtual int getCardinality() const
Returns cardinality of the basis.
This is an interface class for bases whose degrees of freedom can be associated with spatial location...
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Header file for utility class to provide multidimensional containers.
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...
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Quadrilateral cell...
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Quadrilateral 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...