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_Cn_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 using namespace std;
60 using namespace Intrepid;
61 
62 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
63 { \
64  ++nException; \
65  try { \
66  S ; \
67  } \
68  catch (const std::logic_error & err) { \
69  ++throwCounter; \
70  *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
71  *outStream << err.what() << '\n'; \
72  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
73  }; \
74 }
75 
76 
77 int main(int argc, char *argv[]) {
78 
79  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
80 
81  // This little trick lets us print to std::cout only if
82  // a (dummy) command-line argument is provided.
83  int iprint = argc - 1;
84  Teuchos::RCP<std::ostream> outStream;
85  Teuchos::oblackholestream bhs; // outputs nothing
86  if (iprint > 0)
87  outStream = Teuchos::rcp(&std::cout, false);
88  else
89  outStream = Teuchos::rcp(&bhs, false);
90 
91  // Save the format state of the original std::cout.
92  Teuchos::oblackholestream oldFormatState;
93  oldFormatState.copyfmt(std::cout);
94 
95  *outStream \
96  << "===============================================================================\n" \
97  << "| |\n" \
98  << "| Unit Test (Basis_HGRAD_LINE_Cn_FEM) |\n" \
99  << "| |\n" \
100  << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
101  << "| 2) Basis values for VALUE, GRAD, CURL, and Dk operators |\n" \
102  << "| |\n" \
103  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
104  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
105  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
106  << "| |\n" \
107  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
108  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
109  << "| |\n" \
110  << "===============================================================================\n"\
111  << "| TEST 1: Basis creation, exception testing |\n"\
112  << "===============================================================================\n";
113 
114 
115  // Define basis and error flag
116  shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >()); // create cell topology
117  const int deg = 5;
118 
119  // get the points for the basis
120  FieldContainer<double> pts(PointTools::getLatticeSize(line,deg),1);
121  PointTools::getLattice<double,FieldContainer<double> >(pts,line,deg);
122 
123  Basis_HGRAD_LINE_Cn_FEM<double, FieldContainer<double> > lineBasis(deg, pts);
124  int errorFlag = 0;
125 
126  // Initialize throw counter for exception testing
127  int nException = 0;
128  int throwCounter = 0;
129 
130  // Define array containing vertices of the reference Line and a few other points
131  int numIntervals = 100;
132  FieldContainer<double> lineNodes(numIntervals+1, 1);
133  for (int i=0; i<numIntervals+1; i++) {
134  lineNodes(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
135  }
136 
137  // Generic array for the output values; needs to be properly resized depending on the operator type
138  FieldContainer<double> vals;
139 
140  try{
141  // exception #1: DIV cannot be applied to scalar functions
142  // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
143  vals.resize(lineBasis.getCardinality(), lineNodes.dimension(0) );
144  //INTREPID_TEST_COMMAND( lineBasis.getValues(vals, lineNodes, OPERATOR_DIV), throwCounter, nException );
145 
146  // Exceptions 1-5: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
147  // getDofTag() to access invalid array elements thereby causing bounds check exception
148  // exception #1
149  INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(2,0,0), throwCounter, nException );
150  // exception #2
151  INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,1,1), throwCounter, nException );
152  // exception #3
153  INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,0,7), throwCounter, nException );
154  // exception #4
155  INTREPID_TEST_COMMAND( lineBasis.getDofTag(6), throwCounter, nException );
156  // exception #5
157  INTREPID_TEST_COMMAND( lineBasis.getDofTag(-1), throwCounter, nException );
158  // not an exception
159  INTREPID_TEST_COMMAND( lineBasis.getDofTag(5), throwCounter, nException ); --nException;
160 
161 #ifdef HAVE_INTREPID_DEBUG
162  // Exceptions 6-14 test exception handling with incorrectly dimensioned input/output arrays
163  // exception #6: input points array must be of rank-2
164  FieldContainer<double> badPoints1(4, 5, 3);
165  INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
166 
167  // exception #7: dimension 1 in the input point array must equal space dimension of the cell
168  FieldContainer<double> badPoints2(4, 3);
169  INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
170 
171  // exception #8: output values must be of rank-2 for OPERATOR_VALUE
172  FieldContainer<double> badVals1(4, 3, 1);
173  INTREPID_TEST_COMMAND( lineBasis.getValues(badVals1, lineNodes, OPERATOR_VALUE), throwCounter, nException );
174 
175  // exception #9: output values must be of rank-3 for OPERATOR_GRAD
176  FieldContainer<double> badVals2(4, 3);
177  INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_GRAD), throwCounter, nException );
178 
179  // exception #10: output values must be of rank-2 for OPERATOR_D1
180  INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_D1), throwCounter, nException );
181 
182  // exception #11: incorrect 0th dimension of output array (must equal number of basis functions)
183  FieldContainer<double> badVals3(lineBasis.getCardinality() + 1, lineNodes.dimension(0));
184  INTREPID_TEST_COMMAND( lineBasis.getValues(badVals3, lineNodes, OPERATOR_VALUE), throwCounter, nException );
185 
186  // exception #12: incorrect 1st dimension of output array (must equal number of points)
187  FieldContainer<double> badVals4(lineBasis.getCardinality(), lineNodes.dimension(0) + 1);
188  INTREPID_TEST_COMMAND( lineBasis.getValues(badVals4, lineNodes, OPERATOR_VALUE), throwCounter, nException );
189 
190  // exception #13: incorrect 2nd dimension of output array (must equal spatial dimension)
191  FieldContainer<double> badVals5(lineBasis.getCardinality(), lineNodes.dimension(0), 2);
192  INTREPID_TEST_COMMAND( lineBasis.getValues(badVals5, lineNodes, OPERATOR_GRAD), throwCounter, nException );
193 
194 
195  // not an exception
196  FieldContainer<double> goodVals2(lineBasis.getCardinality(), lineNodes.dimension(0));
197  INTREPID_TEST_COMMAND( lineBasis.getValues(goodVals2, lineNodes, OPERATOR_VALUE), throwCounter, nException ); --nException;
198 #endif
199 
200  }
201  catch (const std::logic_error & err) {
202  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
203  *outStream << err.what() << '\n';
204  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
205  errorFlag = -1000;
206  };
207 
208  // Check if number of thrown exceptions matches the one we expect
209  if (throwCounter != nException) {
210  errorFlag++;
211  *outStream << std::setw(70) << "FAILURE! Incorrect number of exceptions." << "\n";
212  }
213 
214  *outStream \
215  << "\n"
216  << "===============================================================================\n" \
217  << "| TEST 2: correctness of basis function values |\n" \
218  << "===============================================================================\n";
219  outStream -> precision(20);
220 
221 
222  try {
223  // Check Kronecker property for Lagrange polynomials.
224  int maxorder = 4;
225 
226  for (int ordi=1; ordi <= maxorder; ordi++) {
227  FieldContainer<double> pts(PointTools::getLatticeSize(line,ordi),1);
228  PointTools::getLattice<double,FieldContainer<double> >(pts,line,ordi);
229  Basis_HGRAD_LINE_Cn_FEM<double, FieldContainer<double> > lineBasis(ordi, pts);
230  FieldContainer<double> vals(lineBasis.getCardinality(),pts.dimension(0));
231  lineBasis.getValues(vals,pts,OPERATOR_VALUE);
232  for (int i=0;i<lineBasis.getCardinality();i++) {
233  for (int j=0;j<pts.dimension(0);j++) {
234  if ( i==j && std::abs( vals(i,j) - 1.0 ) > INTREPID_TOL ) {
235  errorFlag++;
236  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
237  *outStream << " Basis function " << i << " does not have unit value at its node\n";
238  }
239  if ( i!=j && std::abs( vals(i,j) ) > INTREPID_TOL ) {
240  errorFlag++;
241  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
242  *outStream << " Basis function " << i << " does not vanish at node " << j << "\n";
243  }
244  }
245  }
246  }
247  }
248  // Catch unexpected errors
249  catch (const std::logic_error & err) {
250  *outStream << err.what() << "\n\n";
251  errorFlag = -1000;
252  };
253 
254  if (errorFlag != 0)
255  std::cout << "End Result: TEST FAILED\n";
256  else
257  std::cout << "End Result: TEST PASSED\n";
258 
259  // reset format state of std::cout
260  std::cout.copyfmt(oldFormatState);
261 
262  return errorFlag;
263 }
Intrepid::Basis_HGRAD_LINE_Cn_FEM
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1...
Definition: Intrepid_HGRAD_LINE_Cn_FEM.hpp:78
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