Intrepid
test_01.cpp
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 
49 #include "Intrepid_FieldContainer.hpp"
50 #include "Intrepid_HGRAD_LINE_Hermite_FEM.hpp"
51 #include "Intrepid_DefaultCubatureFactory.hpp"
52 #include "Intrepid_ArrayTools.hpp"
53 #include "Intrepid_PointTools.hpp"
54 #include "Intrepid_FunctionSpaceTools.hpp"
55 #include "Teuchos_oblackholestream.hpp"
56 #include "Teuchos_RCP.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
58 
59 #include <iomanip>
60 
61 using namespace std;
62 using 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 
79 int 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 
135  FieldContainer<double> pts(PointTools::getLatticeSize(line,deg),1);
136  PointTools::getLattice<double,FieldContainer<double> >(pts,line,deg);
137  Basis_HGRAD_LINE_Hermite_FEM<double, FieldContainer<double> > lineBasis(pts);
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
144  FieldContainer<double> vals;
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;
253  FieldContainer<double> pts(PointTools::getLatticeSize(line,npts),1);
254  PointTools::getLattice<double,FieldContainer<double> >(pts,line,npts);
255  Basis_HGRAD_LINE_Hermite_FEM<double, FieldContainer<double> > lineBasis(pts);
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
389  FieldContainer<double> pts(PointTools::getLatticeSize(line,npts),1);
390  PointTools::getGaussPoints<double,FieldContainer<double> >(pts,npts);
391  Basis_HGRAD_LINE_Hermite_FEM<double, FieldContainer<double> > lineBasis(pts);
392 
393  int C = lineBasis.getCardinality();
394  int n = C/2;
395 
396  FieldContainer<double> f0(n);
397  FieldContainer<double> f1(n);
398  FieldContainer<double> f2(n);
399  FieldContainer<double> f3(n);
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 }
Intrepid::Basis_HGRAD_LINE_Hermite_FEM
Implements Hermite interpolant basis of degree n on the reference Line cell. The basis has cardinalit...
Definition: Intrepid_HGRAD_LINE_Hermite_FEM.hpp:93
Intrepid_HGRAD_LINE_Hermite_FEM.hpp
Header file for the Intrepid::HGRAD_LINE_Hermite_FEM class.
Intrepid_PointTools.hpp
Header file for utility class to provide point tools, such as barycentric coordinates, equispaced lattices, and warp-blend point distrubtions.
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid_ArrayTools.hpp
Header file for utility class to provide array tools, such as tensor contractions, etc.
Intrepid_DefaultCubatureFactory.hpp
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Intrepid::FieldContainer::resize
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
Definition: Intrepid_FieldContainerDef.hpp:1137
Intrepid_FunctionSpaceTools.hpp
Header file for the Intrepid::FunctionSpaceTools class.
Generated by   doxygen 1.8.5