Intrepid
test_01.cpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
55#include "Teuchos_oblackholestream.hpp"
56#include "Teuchos_RCP.hpp"
57#include "Teuchos_GlobalMPISession.hpp"
58
59#include <iomanip>
60
61using namespace std;
62using namespace Intrepid;
63
64#define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
65{ \
66 ++nException; \
67 try { \
68 S ; \
69 } \
70 catch (const std::logic_error & err) { \
71 ++throwCounter; \
72 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
73 *outStream << err.what() << '\n'; \
74 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
75 }; \
76}
77
78
79int main(int argc, char *argv[]) {
80
81 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
82
83 // This little trick lets us print to std::cout only if
84 // a (dummy) command-line argument is provided.
85 int iprint = argc - 1;
86 Teuchos::RCP<std::ostream> outStream;
87 Teuchos::oblackholestream bhs; // outputs nothing
88 if (iprint > 0)
89 outStream = Teuchos::rcp(&std::cout, false);
90 else
91 outStream = Teuchos::rcp(&bhs, false);
92 // Save the format state of the original std::cout.
93 Teuchos::oblackholestream oldFormatState;
94 oldFormatState.copyfmt(std::cout);
95
96 *outStream \
97 << "===============================================================================\n" \
98 << "| |\n" \
99 << "| Unit Test (Basis_HGRAD_LINE_Hermite_FEM) |\n" \
100 << "| |\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" \
107 << "| |\n" \
108 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
109 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
110 << "| |\n" \
111 << "===============================================================================\n"\
112 << "| TEST 1: Basis creation, exception testing |\n"\
113 << "===============================================================================\n";
114
115
116 // Define basis and error flag
117 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >()); // create cell topology
118
119 int errorFlag = 0;
120
121 // Initialize throw counter for exception testing
122 int nException = 0;
123 int throwCounter = 0;
124
125 // Define array containing vertices of the reference Line and a few other points
126 int numIntervals = 100;
127 FieldContainer<double> lineNodes(numIntervals+1, 1);
128 for (int i=0; i<numIntervals+1; i++) {
129 lineNodes(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
130 }
131
132 try {
133 const int deg = 5;
134
136 PointTools::getLattice<double,FieldContainer<double> >(pts,line,deg);
138
139 *outStream << std::endl;
140
141 lineBasis.printTags( *outStream );
142
143 // Generic array for the output values; needs to be properly resized depending on the operator type
145
146 *outStream << "lineBasis.getCardinality() = " << lineBasis.getCardinality() << std::endl;
147 *outStream << "lineBasis.getDegree() = " << lineBasis.getDegree() << std::endl;
148 *outStream << "lineBasis.getBaseCellTopology() = " << lineBasis.getBaseCellTopology() << std::endl;
149 *outStream << "lineBasis.getBasisType() = " << lineBasis.getBasisType() << std::endl;
150 *outStream << "lineBasis.getCoordinateSystem() = " << lineBasis.getCoordinateSystem() << std::endl;
151 *outStream << std::endl;
152
153
154#ifdef HAVE_INTREPID_DEBUG
155 // exception #1: DIV cannot be applied to scalar functions
156 // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
157 vals.resize(lineBasis.getCardinality(), lineNodes.dimension(0) );
158 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, lineNodes, OPERATOR_DIV), throwCounter, nException );
159#endif // HAVE_INTREPID_DEBUG
160
161 // Exceptions #2-8: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
162 // getDofTag() to access invalid array elements thereby causing bounds check exception
163
164 // exception #2 - There are no two-dimensional subcells
165 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(2,0,0), throwCounter, nException );
166
167 // exception #3 - There are at most two zero-dimensional subcells
168 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(0,2,0), throwCounter, nException );
169
170 // exception #4 - Zero-dimensional subcells have at most 2 DoF
171 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(0,0,2), throwCounter, nException );
172
173 // exception #5 - There is at most one one-dimensional subcell
174 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,1,0), throwCounter, nException );
175
176 // exception #6 - One-dimensional subcell cannot have DoF ordinal larger than
177 // cardinality-1
178 int C = lineBasis.getCardinality();
179 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,0,C), throwCounter, nException );
180
181 // exception #7 - DoF cannot exceed cardinality
182 INTREPID_TEST_COMMAND( lineBasis.getDofTag(C), throwCounter, nException );
183
184 // exception #8 - No negative indices
185 INTREPID_TEST_COMMAND( lineBasis.getDofTag(-1), throwCounter, nException );
186
187 // not an exception
188 INTREPID_TEST_COMMAND( lineBasis.getDofTag(5), throwCounter, nException ); --nException;
189
190#ifdef HAVE_INTREPID_DEBUG
191 // Exceptions 9-16 test exception handling with incorrectly dimensioned input/output arrays
192 // exception #9: input points array must be of rank-2
193 FieldContainer<double> badPoints1(4, 5, 3);
194 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
195
196 // exception #10: dimension 1 in the input point array must equal space dimension of the cell
197 FieldContainer<double> badPoints2(4, 3);
198 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
199
200 // exception #11: output values must be of rank-2 for OPERATOR_VALUE
201 FieldContainer<double> badVals1(4, 3, 1);
202 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals1, lineNodes, OPERATOR_VALUE), throwCounter, nException );
203
204 // exception #12: output values must be of rank-3 for OPERATOR_GRAD
205 FieldContainer<double> badVals2(4, 3);
206 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_GRAD), throwCounter, nException );
207
208 // exception #13: output values must be of rank-2 for OPERATOR_D1
209 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_D1), throwCounter, nException );
210
211 // exception #14: incorrect 0th dimension of output array (must equal number of basis functions)
212 FieldContainer<double> badVals3(lineBasis.getCardinality() + 1, lineNodes.dimension(0));
213 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals3, lineNodes, OPERATOR_VALUE), throwCounter, nException );
214
215 // exception #15: incorrect 1st dimension of output array (must equal number of points)
216 FieldContainer<double> badVals4(lineBasis.getCardinality(), lineNodes.dimension(0) + 1);
217 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals4, lineNodes, OPERATOR_VALUE), throwCounter, nException );
218
219 // exception #16: incorrect 2nd dimension of output array (must equal spatial dimension)
220 FieldContainer<double> badVals5(lineBasis.getCardinality(), lineNodes.dimension(0), 2);
221 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals5, lineNodes, OPERATOR_GRAD), throwCounter, nException );
222
223 // not an exception
224 FieldContainer<double> goodVals2(lineBasis.getCardinality(), lineNodes.dimension(0));
225 INTREPID_TEST_COMMAND( lineBasis.getValues(goodVals2, lineNodes, OPERATOR_VALUE), throwCounter, nException ); --nException;
226#endif // HAVE_INTREPID_DEBUG
227
228 }
229
230
231 catch (const std::logic_error & err) {
232 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
233 *outStream << err.what() << '\n';
234 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
235 errorFlag = -1000;
236 };
237 // Check if number of thrown exceptions matches the one we expect
238 if (throwCounter != nException) {
239 errorFlag++;
240 *outStream << std::setw(70) << "FAILURE! Incorrect number of exceptions." << "\n";
241 }
242
243 *outStream \
244 << "\n"
245 << "===============================================================================\n" \
246 << "| TEST 2: Correctness of basis function values |\n" \
247 << "===============================================================================\n";
248 outStream -> precision(20);
249
250 try {
251
252 int npts=5;
254 PointTools::getLattice<double,FieldContainer<double> >(pts,line,npts);
256
257 FieldContainer<double> vals(lineBasis.getCardinality(),npts+1);
258 FieldContainer<double> ders(lineBasis.getCardinality(),npts+1,1);
259
260 lineBasis.getValues(vals,pts,OPERATOR_VALUE);
261 lineBasis.getValues(ders,pts,OPERATOR_D1);
262
263 int C = lineBasis.getCardinality();
264 int n = C/2;
265
266 // Loop over basis functions
267 for (int i=0; i<n; i++) {
268
269 // Loop over interpolation points
270 for (int j=0; j<n; j++) {
271
272 if ( i == j ) {
273
274 if( std::abs( vals(2*i,j) - 1.0 ) > INTREPID_TOL ) {
275 errorFlag++;
276 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
277 *outStream << " Basis function " << 2*i << " does not have unit value at its node\n";
278 }
279
280 if( std::abs( ders(2*i,j,0) ) > INTREPID_TOL ) {
281 errorFlag++;
282 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
283 *outStream << " Basis function " << 2*i+1 << " does not have zero first derivative at its node\n";
284 }
285
286 if( std::abs( vals(2*i+1,j) ) > INTREPID_TOL ) {
287 errorFlag++;
288 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
289 *outStream << " Basis function " << 2*i+1 << " does not have zero value at its node\n";
290 }
291
292 if( std::abs( ders(2*i+1,j,0) - 1.0 ) > INTREPID_TOL ) {
293 errorFlag++;
294 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
295 *outStream << " Basis function " << 2*i+1 << " does not have unit first derivative at its node\n";
296 }
297 }
298 else { // i != j
299
300 if( std::abs( vals(2*i,j) ) > INTREPID_TOL ) {
301 errorFlag++;
302 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
303 *outStream << " Basis function " << 2*i << " does not vanish at node " << j << "\n";
304 }
305
306 if( std::abs( ders(2*i,j,0) ) > INTREPID_TOL ) {
307 errorFlag++;
308 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
309 *outStream << " Basis function " << 2*i << " does not have zero first derivative at node " << j << "\n";
310 }
311
312 if( std::abs( vals(2*i+1,j) ) > INTREPID_TOL ) {
313 errorFlag++;
314 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
315 *outStream << " Basis function " << 2*i+1 << " does not vanish at node " << j << "\n";
316 }
317
318 if( std::abs( ders(2*i+1,j,0) ) > INTREPID_TOL ) {
319 errorFlag++;
320 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
321 *outStream << " Basis function " << 2*i+1 << " does not have zero first derivative at node " << j << "\n";
322 }
323
324 } // end else i != j
325
326 } // end loop over interpolation points
327
328 } // end loop over basis functions
329
330
331 *outStream << std::setprecision(4);
332
333 *outStream << "\n\nBasis function values on nodes\n" << std::endl;
334
335 // Loop over basis functions
336 for (int i=0; i<C; i++) {
337
338 // Loop over interpolation points
339 for (int j=0; j<n; j++) {
340
341 *outStream << std::setw(12) << vals(i,j);
342 }
343
344 *outStream << std::endl;
345 }
346
347
348 *outStream << "\n\nBasis function derivatives on nodes\n" << std::endl;
349
350 // Loop over basis functions
351 for (int i=0; i<C; i++) {
352
353 // Loop over interpolation points
354 for (int j=0; j<n; j++) {
355
356 *outStream << std::setw(12) << ders(i,j,0);
357 }
358
359 *outStream << std::endl;
360 }
361
362 } // end try
363
364 // Catch unexpected errors
365 catch (const std::logic_error & err) {
366 *outStream << err.what() << "\n\n";
367 errorFlag = -1000;
368 };
369
370 *outStream \
371 << "\n"
372 << "===============================================================================\n" \
373 << "| TEST 3: Correctness of basis function higher derivatives via numerical |\n" \
374 << "| differentiation. |\n" \
375 << "| |\n" \
376 << "| Let f(x_i) = sin(x_i), f'(x_i) = cos(x_i) |\n" \
377 << "| |\n" \
378 << "| and compare the second and third derivatives obtained by differentiating |\n" \
379 << "| the Hermite interpolating polynomial with analytical values |\n" \
380 << "| |\n" \
381 << "===============================================================================\n";
382 outStream -> precision(20);
383
384
385 try {
386
387 int npts = 6;
388 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >()); // create cell topology
390 PointTools::getGaussPoints<double,FieldContainer<double> >(pts,npts);
392
393 int C = lineBasis.getCardinality();
394 int n = C/2;
395
400
401 FieldContainer<double> der2(C,n,1);
402 lineBasis.getValues(der2,pts,OPERATOR_D2);
403
404 FieldContainer<double> der3(C,n,1);
405 lineBasis.getValues(der3,pts,OPERATOR_D3);
406
407 // Loop over interpolation points
408 for( int j=0; j<n; ++j ) {
409 f0(j) = std::sin(pts(j,0));
410 f1(j) = std::cos(pts(j,0));
411 }
412
413 double error2 = 0;
414 double error3 = 0;
415
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);
420
421 f3(j) += f0(i)*der3(2*i,j,0);
422 f3(j) += f1(i)*der3(2*i+1,j,0);
423 }
424
425 error2 += std::pow(f0(j)+f2(j),2);
426 error3 += std::pow(f1(j)+f3(j),2);
427 }
428
429 error2 = std::sqrt(error2);
430 error3 = std::sqrt(error3);
431
432 *outStream << std::setprecision(16);
433
434 int width = 24;
435 std::string bar(20,'-');
436
437
438 *outStream << "\n\n"
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;
444
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;
450
451
452 // Loop over interpolation points
453 for (int j=0; j<n; j++) {
454
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;
460 }
461
462 double errtol = 1e-9;
463
464
465 *outStream << std::endl;
466 *outStream << "|f+f\"| = " << error2 << std::endl;
467 *outStream << "|f'+f'\"| = " << error3 << std::endl;
468
469 if( error2 > errtol ) {
470 errorFlag++;
471 *outStream << std::setw(70) << "FAILURE! Second derivative not within tolerance " << errtol << "\n";
472 }
473 if( error3 > errtol ) {
474 errorFlag++;
475 *outStream << std::setw(70) << "FAILURE! Third derivative not within tolerance " << errtol << "\n";
476 }
477
478 }
479
480 // Catch unexpected errors
481 catch (const std::logic_error & err) {
482 *outStream << err.what() << "\n\n";
483 errorFlag = -1000;
484 };
485
486 if (errorFlag != 0)
487 std::cout << "End Result: TEST FAILED\n";
488 else
489 std::cout << "End Result: TEST PASSED\n";
490
491 // reset format state of std::cout
492 std::cout.copyfmt(oldFormatState);
493
494 return errorFlag;
495}
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for the Intrepid::HGRAD_LINE_Hermite_FEM class.
Header file for utility class to provide point tools, such as barycentric coordinates,...
Implements Hermite interpolant basis of degree n on the reference Line cell. The basis has cardinalit...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
static int getLatticeSize(const shards::CellTopology &cellType, const int order, const int offset=0)
Computes the number of points in a lattice of a given order on a simplex (currently disabled for othe...