54 #include "Teuchos_oblackholestream.hpp"
55 #include "Teuchos_RCP.hpp"
56 #include "Teuchos_GlobalMPISession.hpp"
59 using namespace Intrepid;
61 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
67 catch (const std::logic_error & err) { \
69 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
70 *outStream << err.what() << '\n'; \
71 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
76 int main(
int argc,
char *argv[]) {
78 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
82 int iprint = argc - 1;
83 Teuchos::RCP<std::ostream> outStream;
84 Teuchos::oblackholestream bhs;
86 outStream = Teuchos::rcp(&std::cout,
false);
88 outStream = Teuchos::rcp(&bhs,
false);
91 Teuchos::oblackholestream oldFormatState;
92 oldFormatState.copyfmt(std::cout);
95 <<
"===============================================================================\n" \
97 <<
"| Unit Test (Basis_HGRAD_LINE_Cn_FEM_JACOBI) |\n" \
99 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
100 <<
"| 2) Basis values for VALUE, GRAD, CURL, and Dk operators |\n" \
102 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
103 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
104 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
106 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
107 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
109 <<
"===============================================================================\n"\
110 <<
"| TEST 1: Basis creation, exception testing |\n"\
111 <<
"===============================================================================\n";
115 double alpha = 0.0, beta = 0.0;
121 int throwCounter = 0;
124 int numIntervals = 100;
126 for (
int i=0; i<numIntervals+1; i++) {
127 lineNodes(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
137 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(2,0,0), throwCounter, nException );
139 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,1,1), throwCounter, nException );
141 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,0,7), throwCounter, nException );
143 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,0,5), throwCounter, nException ); --nException;
145 INTREPID_TEST_COMMAND( lineBasis.getDofTag(6), throwCounter, nException );
147 INTREPID_TEST_COMMAND( lineBasis.getDofTag(-1), throwCounter, nException );
149 INTREPID_TEST_COMMAND( lineBasis.getDofTag(5), throwCounter, nException ); --nException;
150 #ifdef HAVE_INTREPID_DEBUG
154 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
158 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
162 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals1, lineNodes, OPERATOR_VALUE), throwCounter, nException );
166 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_GRAD), throwCounter, nException );
169 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_CURL), throwCounter, nException );
172 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_DIV), throwCounter, nException );
175 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_D1), throwCounter, nException );
179 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals3, lineNodes, OPERATOR_VALUE), throwCounter, nException );
183 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals4, lineNodes, OPERATOR_VALUE), throwCounter, nException );
187 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals5, lineNodes, OPERATOR_GRAD), throwCounter, nException );
191 INTREPID_TEST_COMMAND( lineBasis.getValues(goodVals2, lineNodes, OPERATOR_VALUE), throwCounter, nException ); --nException;
195 catch (
const std::logic_error & err) {
196 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
197 *outStream << err.what() <<
'\n';
198 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
203 if (throwCounter != nException) {
205 *outStream << std::setw(70) <<
"FAILURE! Incorrect number of exceptions." <<
"\n";
211 <<
"===============================================================================\n"\
212 <<
"| TEST 3: orthogonality of basis functions |\n"\
213 <<
"===============================================================================\n";
215 outStream -> precision(20);
223 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
224 for (
int ordi=0; ordi < maxorder; ordi++) {
226 Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasisLeft =
229 for (
int ordj=0; ordj < maxorder; ordj++) {
232 Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasisRight =
236 Teuchos::RCP<Cubature<double> > lineCub = cubFactory.
create(line, ordi+ordj);
237 int numPoints = lineCub->getNumPoints();
241 lineCub->getCubature(cubPoints, cubWeights);
243 for (
int i=0; i<numPoints; i++) { cubWeightsC(0,i) = cubWeights(i); }
247 int numFieldsLeft = lineBasisLeft ->getCardinality();
248 int numFieldsRight = lineBasisRight->getCardinality();
250 valsRight(numFieldsRight,numPoints);
251 lineBasisLeft ->getValues(valsLeft, cubPoints, OPERATOR_VALUE);
252 lineBasisRight->getValues(valsRight, cubPoints, OPERATOR_VALUE);
256 valsRightC(1, numFieldsRight,numPoints),
257 massMatrix(1, numFieldsLeft, numFieldsRight);
258 ArrayTools::cloneFields<double>(valsLeftC, valsLeft);
259 ArrayTools::cloneFields<double>(valsRightC, valsRight);
260 ArrayTools::scalarMultiplyDataField<double>(valsRightC, cubWeightsC, valsRightC);
261 FunctionSpaceTools::integrate<double>(massMatrix, valsLeftC, valsRightC, COMP_CPP);
264 for (
int i=0; i<numFieldsLeft; i++) {
265 for (
int j=0; j<numFieldsRight; j++) {
268 if ( std::abs(massMatrix(0,i,j)-(double)(2.0/(2.0*j+1.0))) > INTREPID_TOL ) {
269 *outStream <<
"Incorrect ii (\"diagonal\") value for i=" << i <<
", j=" << j <<
": "
270 << massMatrix(0,i,j) <<
" != " <<
"2/(2*" << j <<
"+1)\n\n";
275 if ( std::abs(massMatrix(0,i,j)) > INTREPID_TOL ) {
276 *outStream <<
"Incorrect ij (\"off-diagonal\") value for i=" << i <<
", j=" << j <<
": "
277 << massMatrix(0,i,j) <<
" != " <<
"0\n\n";
289 catch (
const std::logic_error & err) {
290 *outStream << err.what() <<
"\n\n";
296 <<
"===============================================================================\n"\
297 <<
"| TEST 4: correctness of basis function derivatives |\n"\
298 <<
"===============================================================================\n";
300 outStream -> precision(20);
303 double basisValues[] = {
304 1.000000000000000, 1.000000000000000, 1.000000000000000, \
305 1.000000000000000, -1.000000000000000, -0.3333333333333333, \
306 0.3333333333333333, 1.000000000000000, 1.000000000000000, \
307 -0.3333333333333333, -0.3333333333333333, 1.000000000000000, \
308 -1.000000000000000, 0.4074074074074074, -0.4074074074074074, \
311 double basisD1Values[] =
312 {0, 0, 0, 0, 1.000000000000000, 1.000000000000000, 1.000000000000000, \
313 1.000000000000000, -3.000000000000000, -1.000000000000000, \
314 1.000000000000000, 3.000000000000000, 6.000000000000000, \
315 -0.6666666666666667, -0.6666666666666667, 6.000000000000000};
317 double basisD2Values[] =
318 {0, 0, 0, 0, 0, 0, 0, 0, 3.000000000000000, 3.000000000000000, \
319 3.000000000000000, 3.000000000000000, -15.00000000000000, \
320 -5.000000000000000, 5.000000000000000, 15.00000000000000};
322 double basisD3Values[] =
323 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.00000000000000, \
324 15.00000000000000, 15.00000000000000, 15.00000000000000};
330 int numIntervals = 3;
333 for (
int i=0; i<numIntervals+1; i++) {
334 lineNodes3(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
336 int numFields = lineBasis3.getCardinality();
340 vals.
resize(numFields, numPoints);
341 lineBasis3.getValues(vals,lineNodes3,OPERATOR_VALUE);
342 for (
int i = 0; i < numFields; i++) {
343 for (
int j = 0; j < numPoints; j++) {
346 int l = j + i * numPoints;
347 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
349 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
352 *outStream <<
" At multi-index { ";
353 *outStream << i <<
" ";*outStream << j <<
" ";
354 *outStream <<
"} computed value: " << vals(i,j)
355 <<
" but reference value: " << basisValues[l] <<
"\n";
361 vals.
resize(numFields, numPoints,1);
362 lineBasis3.getValues(vals,lineNodes3,OPERATOR_D1);
363 for (
int i = 0; i < numFields; i++) {
364 for (
int j = 0; j < numPoints; j++) {
367 int l = j + i * numPoints;
368 if (std::abs(vals(i,j,0) - basisD1Values[l]) > INTREPID_TOL) {
370 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
373 *outStream <<
" At multi-index { ";
374 *outStream << i <<
" ";*outStream << j <<
" ";
375 *outStream <<
"} computed value: " << vals(i,j,0)
376 <<
" but reference value: " << basisD1Values[l] <<
"\n";
381 vals.
resize(numFields, numPoints,1);
382 lineBasis3.getValues(vals,lineNodes3,OPERATOR_D2);
383 for (
int i = 0; i < numFields; i++) {
384 for (
int j = 0; j < numPoints; j++) {
387 int l = j + i * numPoints;
388 if (std::abs(vals(i,j,0) - basisD2Values[l]) > INTREPID_TOL) {
390 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
393 *outStream <<
" At multi-index { ";
394 *outStream << i <<
" ";*outStream << j <<
" ";
395 *outStream <<
"} computed value: " << vals(i,j,0)
396 <<
" but reference value: " << basisD2Values[l] <<
"\n";
401 vals.
resize(numFields, numPoints,1);
402 lineBasis3.getValues(vals,lineNodes3,OPERATOR_D3);
403 for (
int i = 0; i < numFields; i++) {
404 for (
int j = 0; j < numPoints; j++) {
407 int l = j + i * numPoints;
408 if (std::abs(vals(i,j,0) - basisD3Values[l]) > INTREPID_TOL) {
410 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
413 *outStream <<
" At multi-index { ";
414 *outStream << i <<
" ";*outStream << j <<
" ";
415 *outStream <<
"} computed value: " << vals(i,j,0)
416 <<
" but reference value: " << basisD3Values[l] <<
"\n";
422 catch (
const std::logic_error & err) {
423 *outStream << err.what() <<
"\n\n";
429 std::cout <<
"End Result: TEST FAILED\n";
431 std::cout <<
"End Result: TEST PASSED\n";
434 std::cout.copyfmt(oldFormatState);
Header file for the Intrepid::HGRAD_LINE_Cn_FEM class.
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1...
int dimension(const int whichDim) const
Returns the specified dimension.
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...
A factory class that generates specific instances of cubatures.
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > °ree)
Factory method.