55 #include "Teuchos_oblackholestream.hpp"
56 #include "Teuchos_RCP.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
62 using namespace Intrepid;
64 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
70 catch (const std::logic_error & err) { \
72 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
73 *outStream << err.what() << '\n'; \
74 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
79 int main(
int argc,
char *argv[]) {
81 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
85 int iprint = argc - 1;
86 Teuchos::RCP<std::ostream> outStream;
87 Teuchos::oblackholestream bhs;
89 outStream = Teuchos::rcp(&std::cout,
false);
91 outStream = Teuchos::rcp(&bhs,
false);
93 Teuchos::oblackholestream oldFormatState;
94 oldFormatState.copyfmt(std::cout);
97 <<
"===============================================================================\n" \
99 <<
"| Unit Test (Basis_HGRAD_LINE_Hermite_FEM) |\n" \
101 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
102 <<
"| 2) Basis values for VALUE, GRAD, CURL, and Dk operators |\n" \
103 <<
"| 3) Numerical differentiation with Herite Interpolants |\n" \
104 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
105 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
106 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
108 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
109 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
111 <<
"===============================================================================\n"\
112 <<
"| TEST 1: Basis creation, exception testing |\n"\
113 <<
"===============================================================================\n";
117 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
121 PointTools::getLattice<double,FieldContainer<double> >(pts,line,deg);
124 *outStream << std::endl;
126 lineBasis.printTags( *outStream );
132 int throwCounter = 0;
135 int numIntervals = 100;
137 for (
int i=0; i<numIntervals+1; i++) {
138 lineNodes(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
144 *outStream <<
"lineBasis.getCardinality() = " << lineBasis.getCardinality() << std::endl;
145 *outStream <<
"lineBasis.getDegree() = " << lineBasis.getDegree() << std::endl;
146 *outStream <<
"lineBasis.getBaseCellTopology() = " << lineBasis.getBaseCellTopology() << std::endl;
147 *outStream <<
"lineBasis.getBasisType() = " << lineBasis.getBasisType() << std::endl;
148 *outStream <<
"lineBasis.getCoordinateSystem() = " << lineBasis.getCoordinateSystem() << std::endl;
149 *outStream << std::endl;
153 #ifdef HAVE_INTREPID_DEBUG
156 vals.
resize(lineBasis.getCardinality(), lineNodes.dimension(0) );
157 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, lineNodes, OPERATOR_DIV), throwCounter, nException );
158 #endif // HAVE_INTREPID_DEBUG
164 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(2,0,0), throwCounter, nException );
167 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(0,2,0), throwCounter, nException );
170 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(0,0,2), throwCounter, nException );
173 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,1,0), throwCounter, nException );
177 int C = lineBasis.getCardinality();
178 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,0,C), throwCounter, nException );
181 INTREPID_TEST_COMMAND( lineBasis.getDofTag(C), throwCounter, nException );
184 INTREPID_TEST_COMMAND( lineBasis.getDofTag(-1), throwCounter, nException );
187 INTREPID_TEST_COMMAND( lineBasis.getDofTag(5), throwCounter, nException ); --nException;
189 #ifdef HAVE_INTREPID_DEBUG
193 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
197 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
201 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals1, lineNodes, OPERATOR_VALUE), throwCounter, nException );
205 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_GRAD), throwCounter, nException );
208 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_D1), throwCounter, nException );
212 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals3, lineNodes, OPERATOR_VALUE), throwCounter, nException );
216 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals4, lineNodes, OPERATOR_VALUE), throwCounter, nException );
220 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals5, lineNodes, OPERATOR_GRAD), throwCounter, nException );
224 INTREPID_TEST_COMMAND( lineBasis.getValues(goodVals2, lineNodes, OPERATOR_VALUE), throwCounter, nException ); --nException;
225 #endif // HAVE_INTREPID_DEBUG
230 catch (
const std::logic_error & err) {
231 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
232 *outStream << err.what() <<
'\n';
233 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
237 if (throwCounter != nException) {
239 *outStream << std::setw(70) <<
"FAILURE! Incorrect number of exceptions." <<
"\n";
244 <<
"===============================================================================\n" \
245 <<
"| TEST 2: Correctness of basis function values |\n" \
246 <<
"===============================================================================\n";
247 outStream -> precision(20);
252 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
254 PointTools::getLattice<double,FieldContainer<double> >(pts,line,npts);
260 lineBasis.getValues(vals,pts,OPERATOR_VALUE);
261 lineBasis.getValues(ders,pts,OPERATOR_D1);
263 int C = lineBasis.getCardinality();
267 for (
int i=0; i<n; i++) {
270 for (
int j=0; j<n; j++) {
274 if( std::abs( vals(2*i,j) - 1.0 ) > INTREPID_TOL ) {
276 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
277 *outStream <<
" Basis function " << 2*i <<
" does not have unit value at its node\n";
280 if( std::abs( ders(2*i,j,0) ) > INTREPID_TOL ) {
282 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
283 *outStream <<
" Basis function " << 2*i+1 <<
" does not have zero first derivative at its node\n";
286 if( std::abs( vals(2*i+1,j) ) > INTREPID_TOL ) {
288 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
289 *outStream <<
" Basis function " << 2*i+1 <<
" does not have zero value at its node\n";
292 if( std::abs( ders(2*i+1,j,0) - 1.0 ) > INTREPID_TOL ) {
294 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
295 *outStream <<
" Basis function " << 2*i+1 <<
" does not have unit first derivative at its node\n";
300 if( std::abs( vals(2*i,j) ) > INTREPID_TOL ) {
302 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
303 *outStream <<
" Basis function " << 2*i <<
" does not vanish at node " << j <<
"\n";
306 if( std::abs( ders(2*i,j,0) ) > INTREPID_TOL ) {
308 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
309 *outStream <<
" Basis function " << 2*i <<
" does not have zero first derivative at node " << j <<
"\n";
312 if( std::abs( vals(2*i+1,j) ) > INTREPID_TOL ) {
314 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
315 *outStream <<
" Basis function " << 2*i+1 <<
" does not vanish at node " << j <<
"\n";
318 if( std::abs( ders(2*i+1,j,0) ) > INTREPID_TOL ) {
320 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
321 *outStream <<
" Basis function " << 2*i+1 <<
" does not have zero first derivative at node " << j <<
"\n";
331 *outStream << std::setprecision(4);
333 *outStream <<
"\n\nBasis function values on nodes\n" << std::endl;
336 for (
int i=0; i<C; i++) {
339 for (
int j=0; j<n; j++) {
341 *outStream << std::setw(12) << vals(i,j);
344 *outStream << std::endl;
348 *outStream <<
"\n\nBasis function derivatives on nodes\n" << std::endl;
351 for (
int i=0; i<C; i++) {
354 for (
int j=0; j<n; j++) {
356 *outStream << std::setw(12) << ders(i,j,0);
359 *outStream << std::endl;
365 catch (
const std::logic_error & err) {
366 *outStream << err.what() <<
"\n\n";
372 <<
"===============================================================================\n" \
373 <<
"| TEST 3: Correctness of basis function higher derivatives via numerical |\n" \
374 <<
"| differentiation. |\n" \
376 <<
"| Let f(x_i) = sin(x_i), f'(x_i) = cos(x_i) |\n" \
378 <<
"| and compare the second and third derivatives obtained by differentiating |\n" \
379 <<
"| the Hermite interpolating polynomial with analytical values |\n" \
381 <<
"===============================================================================\n";
382 outStream -> precision(20);
388 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
390 PointTools::getGaussPoints<double,FieldContainer<double> >(pts,npts);
393 int C = lineBasis.getCardinality();
402 lineBasis.getValues(der2,pts,OPERATOR_D2);
405 lineBasis.getValues(der3,pts,OPERATOR_D3);
408 for(
int j=0; j<n; ++j ) {
409 f0(j) = std::sin(pts(j,0));
410 f1(j) = std::cos(pts(j,0));
416 for(
int j=0; j<n; ++j ) {
417 for(
int i=0; i<n; ++i ) {
418 f2(j) += f0(i)*der2(2*i,j,0);
419 f2(j) += f1(i)*der2(2*i+1,j,0);
421 f3(j) += f0(i)*der3(2*i,j,0);
422 f3(j) += f1(i)*der3(2*i+1,j,0);
425 error2 += std::pow(f0(j)+f2(j),2);
426 error3 += std::pow(f1(j)+f3(j),2);
429 error2 = std::sqrt(error2);
430 error3 = std::sqrt(error3);
432 *outStream << std::setprecision(16);
435 std::string bar(20,
'-');
439 << std::setw(width) <<
"x_i"
440 << std::setw(width) <<
"exact f(x_i)"
441 << std::setw(width) <<
"exact f'(x_i)"
442 << std::setw(width) <<
"computed f\"(x_i)"
443 << std::setw(width) <<
"computed f'\"(x_i)" << std::endl;
445 *outStream << std::setw(width) << bar
446 << std::setw(width) << bar
447 << std::setw(width) << bar
448 << std::setw(width) << bar
449 << std::setw(width) << bar << std::endl;
453 for (
int j=0; j<n; j++) {
455 *outStream << std::setw(width) << pts(j,0)
456 << std::setw(width) << f0(j)
457 << std::setw(width) << f1(j)
458 << std::setw(width) << f2(j)
459 << std::setw(width) << f3(j) << std::endl;
462 double errtol = 1e-9;
465 *outStream << std::endl;
466 *outStream <<
"|f+f\"| = " << error2 << std::endl;
467 *outStream <<
"|f'+f'\"| = " << error3 << std::endl;
469 if( error2 > errtol ) {
471 *outStream << std::setw(70) <<
"FAILURE! Second derivative not within tolerance " << errtol <<
"\n";
473 if( error3 > errtol ) {
475 *outStream << std::setw(70) <<
"FAILURE! Third derivative not within tolerance " << errtol <<
"\n";
481 catch (
const std::logic_error & err) {
482 *outStream << err.what() <<
"\n\n";
487 std::cout <<
"End Result: TEST FAILED\n";
489 std::cout <<
"End Result: TEST PASSED\n";
492 std::cout.copyfmt(oldFormatState);
Implements Hermite interpolant basis of degree n on the reference Line cell. The basis has cardinalit...
Header file for the Intrepid::HGRAD_LINE_Hermite_FEM class.
Header file for utility class to provide multidimensional containers.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...