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;
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
151 INTREPID_TEST_COMMAND( quadBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
155 INTREPID_TEST_COMMAND( quadBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
159 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
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 );
173 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException );
177 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException );
181 INTREPID_TEST_COMMAND( quadBasis.
getValues(badVals5, quadNodes, OPERATOR_GRAD), throwCounter, nException );
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,
310 int numPoints = quadNodes.dimension(0);
312 int D2Cardin = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
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++) {
426 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
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;
438 int ord = Intrepid::getOperatorOrder(op);
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 =
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);
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...