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";
110 Basis_HGRAD_LINE_C1_FEM<double, FieldContainer<double> > lineBasis;
115 int throwCounter = 0;
118 FieldContainer<double> lineNodes(4, 1);
119 lineNodes(0,0) = -1.0;
120 lineNodes(1,0) = 1.0;
121 lineNodes(2,0) = 0.0;
122 lineNodes(3,0) = 0.5;
125 FieldContainer<double> vals;
129 vals.resize(lineBasis.getCardinality(), lineNodes.dimension(0) );
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
147 FieldContainer<double> badPoints1(4, 5, 3);
148 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
151 FieldContainer<double> badPoints2(4, 2);
152 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
155 FieldContainer<double> badVals1(4, 3, 1);
156 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals1, lineNodes, OPERATOR_VALUE), throwCounter, nException );
159 FieldContainer<double> badVals2(4, 3);
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 );
172 FieldContainer<double> badVals3(lineBasis.getCardinality() + 1, lineNodes.dimension(0));
173 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals3, lineNodes, OPERATOR_VALUE), throwCounter, nException );
176 FieldContainer<double> badVals4(lineBasis.getCardinality(), lineNodes.dimension(0) + 1);
177 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals4, lineNodes, OPERATOR_VALUE), throwCounter, nException );
180 FieldContainer<double> badVals5(lineBasis.getCardinality(), lineNodes.dimension(0), 2);
181 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals5, lineNodes, OPERATOR_GRAD), throwCounter, nException );
184 FieldContainer<double> badVals6(lineBasis.getCardinality(), lineNodes.dimension(0), 40);
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[] = {
289 int numFields = lineBasis.getCardinality();
290 int numPoints = lineNodes.dimension(0);
291 int spaceDim = lineBasis.getBaseCellTopology().getDimension();
294 FieldContainer<double> vals;
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;
415 vals.getMultiIndex(myIndex,i);
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.