51 #include "Teuchos_oblackholestream.hpp"
52 #include "Teuchos_RCP.hpp"
53 #include "Teuchos_GlobalMPISession.hpp"
56 using namespace Intrepid;
58 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
64 catch (const std::logic_error & err) { \
66 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
67 *outStream << err.what() << '\n'; \
68 *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_TRI_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 (dridzal@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";
115 int throwCounter = 0;
119 triNodes(0,0) = 0.0; triNodes(0,1) = 0.0;
120 triNodes(1,0) = 1.0; triNodes(1,1) = 0.0;
121 triNodes(2,0) = 0.0; triNodes(2,1) = 1.0;
122 triNodes(3,0) = 0.5; triNodes(3,1) = 0.5;
123 triNodes(4,0) = 0.0; triNodes(4,1) = 0.75;
132 INTREPID_TEST_COMMAND( triBasis.
getValues(vals, triNodes, OPERATOR_DIV), throwCounter, nException );
137 INTREPID_TEST_COMMAND( triBasis.
getDofOrdinal(2,0,0), throwCounter, nException );
139 INTREPID_TEST_COMMAND( triBasis.
getDofOrdinal(1,1,1), throwCounter, nException );
141 INTREPID_TEST_COMMAND( triBasis.
getDofOrdinal(0,4,0), throwCounter, nException );
143 INTREPID_TEST_COMMAND( triBasis.
getDofTag(5), throwCounter, nException );
145 INTREPID_TEST_COMMAND( triBasis.
getDofTag(-1), throwCounter, nException );
147 #ifdef HAVE_INTREPID_DEBUG
151 INTREPID_TEST_COMMAND( triBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
155 INTREPID_TEST_COMMAND( triBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
159 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals1, triNodes, OPERATOR_VALUE), throwCounter, nException );
163 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals2, triNodes, OPERATOR_GRAD), throwCounter, nException );
166 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals2, triNodes, OPERATOR_CURL), throwCounter, nException );
169 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals2, triNodes, OPERATOR_D2), throwCounter, nException );
173 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals3, triNodes, OPERATOR_VALUE), throwCounter, nException );
177 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals4, triNodes, OPERATOR_VALUE), throwCounter, nException );
181 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals5, triNodes, OPERATOR_GRAD), throwCounter, nException );
185 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals6, triNodes, OPERATOR_D2), throwCounter, nException );
188 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals6, triNodes, 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 = triBasis.
getAllDofTags();
215 for (
unsigned i = 0; i < allTags.size(); i++) {
216 int bfOrd = triBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
218 std::vector<int> myTag = triBasis.
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";
240 std::vector<int> myTag = triBasis.
getDofTag(bfOrd);
241 int myBfOrd = triBasis.
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 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
282 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
283 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
284 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
285 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
289 double basisCurls[] = {
290 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0,
291 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0,
292 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0,
293 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0,
294 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0
301 int numPoints = triNodes.dimension(0);
308 vals.
resize(numFields, numPoints);
309 triBasis.
getValues(vals, triNodes, OPERATOR_VALUE);
310 for (
int i = 0; i < numFields; i++) {
311 for (
int j = 0; j < numPoints; j++) {
312 int l = i + j * numFields;
313 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
315 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
318 *outStream <<
" At multi-index { ";
319 *outStream << i <<
" ";*outStream << j <<
" ";
320 *outStream <<
"} computed value: " << vals(i,j)
321 <<
" but reference value: " << basisValues[l] <<
"\n";
327 vals.
resize(numFields, numPoints, spaceDim);
328 triBasis.
getValues(vals, triNodes, OPERATOR_GRAD);
329 for (
int i = 0; i < numFields; i++) {
330 for (
int j = 0; j < numPoints; j++) {
331 for (
int k = 0; k < spaceDim; k++) {
332 int l = k + i * spaceDim + j * spaceDim * numFields;
333 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
335 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
338 *outStream <<
" At multi-index { ";
339 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
340 *outStream <<
"} computed grad component: " << vals(i,j,k)
341 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
348 triBasis.
getValues(vals, triNodes, OPERATOR_D1);
349 for (
int i = 0; i < numFields; i++) {
350 for (
int j = 0; j < numPoints; j++) {
351 for (
int k = 0; k < spaceDim; k++) {
352 int l = k + i * spaceDim + j * spaceDim * numFields;
353 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
355 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
358 *outStream <<
" At multi-index { ";
359 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
360 *outStream <<
"} computed D1 component: " << vals(i,j,k)
361 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
369 vals.
resize(numFields, numPoints, spaceDim);
370 triBasis.
getValues(vals, triNodes, OPERATOR_CURL);
371 for (
int i = 0; i < numFields; i++) {
372 for (
int j = 0; j < numPoints; j++) {
373 for (
int k = 0; k < spaceDim; k++) {
374 int l = k + i * spaceDim + j * spaceDim * numFields;
375 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
377 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
380 *outStream <<
" At multi-index { ";
381 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
382 *outStream <<
"} computed curl component: " << vals(i,j,k)
383 <<
" but reference curl component: " << basisCurls[l] <<
"\n";
390 for(EOperator op = OPERATOR_D2; op < OPERATOR_MAX; op++) {
393 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
394 vals.
resize(numFields, numPoints, DkCardin);
397 for (
int i = 0; i < vals.
size(); i++) {
398 if (std::abs(vals[i]) > INTREPID_TOL) {
400 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
403 std::vector<int> myIndex;
405 int ord = Intrepid::getOperatorOrder(op);
406 *outStream <<
" At multi-index { ";
407 for(
int j = 0; j < vals.
rank(); j++) {
408 *outStream << myIndex[j] <<
" ";
410 *outStream <<
"} computed D"<< ord <<
" component: " << vals[i]
411 <<
" but reference D" << ord <<
" component: 0 \n";
418 catch (
const std::logic_error & err) {
419 *outStream << err.what() <<
"\n\n";
425 <<
"===============================================================================\n"\
426 <<
"| TEST 4: correctness of DoF locations |\n"\
427 <<
"===============================================================================\n";
430 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
432 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
439 #ifdef HAVE_INTREPID_DEBUG
441 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
443 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
445 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
448 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
450 if (throwCounter != nException) {
452 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
456 basis->getValues(bvals, cvals, OPERATOR_VALUE);
458 for (
int i=0; i<bvals.dimension(0); i++) {
459 for (
int j=0; j<bvals.dimension(1); j++) {
460 if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) {
462 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);
463 *outStream << buffer;
465 else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) {
467 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);
468 *outStream << buffer;
474 catch (
const std::logic_error & err){
475 *outStream << err.what() <<
"\n\n";
480 std::cout <<
"End Result: TEST FAILED\n";
482 std::cout <<
"End Result: TEST PASSED\n";
485 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...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
Header file for the Intrepid::HGRAD_TRI_C1_FEM class.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Header file for utility class to provide multidimensional containers.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Triangle cell...
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...
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
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...