51#include "Teuchos_oblackholestream.hpp"
52#include "Teuchos_RCP.hpp"
53#include "Teuchos_GlobalMPISession.hpp"
56using 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"; \
72int main(
int argc,
char *argv[]) {
73 Teuchos::oblackholestream oldFormatState;
74 oldFormatState.copyfmt(std::cout);
76 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
80 int iprint = argc - 1;
81 Teuchos::RCP<std::ostream> outStream;
82 Teuchos::oblackholestream bhs;
84 outStream = Teuchos::rcp(&std::cout,
false);
86 outStream = Teuchos::rcp(&bhs,
false);
91 <<
"===============================================================================\n" \
93 <<
"| Unit Test (Basis_HGRAD_LINE_C1_FEM) |\n" \
95 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
96 <<
"| 2) Basis values for VALUE, GRAD, CURL, DIV, and Dk operators |\n" \
98 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
99 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
100 <<
"| Kara Peterson (dridzal@sandia.gov). |\n" \
102 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
103 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
105 <<
"===============================================================================\n"\
106 <<
"| TEST 1: Basis creation, exception testing |\n"\
107 <<
"===============================================================================\n";
115 int throwCounter = 0;
119 lineNodes(0,0) = -1.0;
120 lineNodes(1,0) = 1.0;
121 lineNodes(2,0) = 0.0;
122 lineNodes(3,0) = 0.5;
134 INTREPID_TEST_COMMAND( lineBasis.
getDofOrdinal(2,0,0), throwCounter, nException );
136 INTREPID_TEST_COMMAND( lineBasis.
getDofOrdinal(1,1,1), throwCounter, nException );
138 INTREPID_TEST_COMMAND( lineBasis.
getDofOrdinal(0,4,0), throwCounter, nException );
140 INTREPID_TEST_COMMAND( lineBasis.
getDofTag(5), throwCounter, nException );
142 INTREPID_TEST_COMMAND( lineBasis.
getDofTag(-1), throwCounter, nException );
144#ifdef HAVE_INTREPID_DEBUG
148 INTREPID_TEST_COMMAND( lineBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
152 INTREPID_TEST_COMMAND( lineBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
156 INTREPID_TEST_COMMAND( lineBasis.
getValues(badVals1, lineNodes, OPERATOR_VALUE), throwCounter, nException );
160 INTREPID_TEST_COMMAND( lineBasis.
getValues(badVals2, lineNodes, OPERATOR_GRAD), throwCounter, nException );
163 INTREPID_TEST_COMMAND( lineBasis.
getValues(badVals2, lineNodes, OPERATOR_DIV), throwCounter, nException );
166 INTREPID_TEST_COMMAND( lineBasis.
getValues(badVals2, lineNodes, OPERATOR_CURL), throwCounter, nException );
169 INTREPID_TEST_COMMAND( lineBasis.
getValues(badVals2, lineNodes, OPERATOR_D2), throwCounter, nException );
173 INTREPID_TEST_COMMAND( lineBasis.
getValues(badVals3, lineNodes, OPERATOR_VALUE), throwCounter, nException );
177 INTREPID_TEST_COMMAND( lineBasis.
getValues(badVals4, lineNodes, OPERATOR_VALUE), throwCounter, nException );
181 INTREPID_TEST_COMMAND( lineBasis.
getValues(badVals5, lineNodes, OPERATOR_GRAD), throwCounter, nException );
185 INTREPID_TEST_COMMAND( lineBasis.
getValues(badVals6, lineNodes, OPERATOR_D2), throwCounter, nException );
188 INTREPID_TEST_COMMAND( lineBasis.
getValues(badVals6, lineNodes, 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 = lineBasis.
getAllDofTags();
215 for (
unsigned i = 0; i < allTags.size(); i++) {
216 int bfOrd = lineBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
218 std::vector<int> myTag = lineBasis.
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 < lineBasis.
getCardinality(); bfOrd++) {
240 std::vector<int> myTag = lineBasis.
getDofTag(bfOrd);
241 int myBfOrd = lineBasis.
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[] = {
279 double basisDerivs[] = {
290 int numPoints = lineNodes.dimension(0);
297 vals.
resize(numFields, numPoints);
298 lineBasis.
getValues(vals, lineNodes, OPERATOR_VALUE);
299 for (
int i = 0; i < numFields; i++) {
300 for (
int j = 0; j < numPoints; j++) {
301 int l = i + j * numFields;
302 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
304 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
307 *outStream <<
" At multi-index { ";
308 *outStream << i <<
" ";*outStream << j <<
" ";
309 *outStream <<
"} computed value: " << vals(i,j)
310 <<
" but reference value: " << basisValues[l] <<
"\n";
316 vals.
resize(numFields, numPoints, spaceDim);
317 lineBasis.
getValues(vals, lineNodes, OPERATOR_GRAD);
318 for (
int i = 0; i < numFields; i++) {
319 for (
int j = 0; j < numPoints; j++) {
320 for (
int k = 0; k < spaceDim; k++) {
321 int l = k + i * spaceDim + j * spaceDim * numFields;
322 if (std::abs(vals(i,j,k) - basisDerivs[l]) > INTREPID_TOL) {
324 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
327 *outStream <<
" At multi-index { ";
328 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
329 *outStream <<
"} computed grad component: " << vals(i,j,k)
330 <<
" but reference grad component: " << basisDerivs[l] <<
"\n";
337 lineBasis.
getValues(vals, lineNodes, OPERATOR_D1);
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) - basisDerivs[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 D1 component: " << vals(i,j,k)
350 <<
" but reference D1 component: " << basisDerivs[l] <<
"\n";
358 vals.
resize(numFields, numPoints, spaceDim);
359 lineBasis.
getValues(vals, lineNodes, OPERATOR_CURL);
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) - basisDerivs[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 curl component: " << vals(i,j,k)
372 <<
" but reference curl component: " << basisDerivs[l] <<
"\n";
380 lineBasis.
getValues(vals, lineNodes, OPERATOR_DIV);
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) - basisDerivs[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 D1 component: " << vals(i,j,k)
393 <<
" but reference D1 component: " << basisDerivs[l] <<
"\n";
401 for(EOperator op = OPERATOR_D2; op < OPERATOR_MAX; op++) {
405 vals.
resize(numFields, numPoints, DkCardin);
407 lineBasis.
getValues(vals, lineNodes, op);
408 for (
int i = 0; i < vals.
size(); i++) {
409 if (std::abs(vals[i]) > INTREPID_TOL) {
411 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
414 std::vector<int> myIndex;
417 *outStream <<
" At multi-index { ";
418 for(
int j = 0; j < vals.
rank(); j++) {
419 *outStream << myIndex[j] <<
" ";
421 *outStream <<
"} computed D"<< ord <<
" component: " << vals[i]
422 <<
" but reference D" << ord <<
" component: 0 \n";
429 catch (
const std::logic_error & err) {
430 *outStream << err.what() <<
"\n\n";
435 std::cout <<
"End Result: TEST FAILED\n";
437 std::cout <<
"End Result: TEST PASSED\n";
440 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_LINE_C1_FEM class.
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 Line cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Line cell.
virtual int getDofOrdinal(const int subcDim, const int subcOrd, const int subcDofOrd)
DoF tag to ordinal lookup.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
virtual int getCardinality() const
Returns cardinality of the basis.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
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.