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 
50 #include "Intrepid_PointTools.hpp"
51 #include "Teuchos_oblackholestream.hpp"
52 #include "Teuchos_RCP.hpp"
53 #include "Teuchos_GlobalMPISession.hpp"
54 
55 using namespace std;
56 using namespace Intrepid;
57 
58 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
59 { \
60  ++nException; \
61  try { \
62  S ; \
63  } \
64  catch (const std::logic_error & err) { \
65  ++throwCounter; \
66  *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
67  *outStream << err.what() << '\n'; \
68  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
69  }; \
70 }
71 
72 int main(int argc, char *argv[]) {
73 
74  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
75 
76  // This little trick lets us print to std::cout only if
77  // a (dummy) command-line argument is provided.
78  int iprint = argc - 1;
79  Teuchos::RCP<std::ostream> outStream;
80  Teuchos::oblackholestream bhs; // outputs nothing
81  if (iprint > 0)
82  outStream = Teuchos::rcp(&std::cout, false);
83  else
84  outStream = Teuchos::rcp(&bhs, false);
85 
86  // Save the format state of the original std::cout.
87  Teuchos::oblackholestream oldFormatState;
88  oldFormatState.copyfmt(std::cout);
89 
90  *outStream \
91  << "===============================================================================\n" \
92  << "| |\n" \
93  << "| Unit Test (Basis_HCURL_QUAD_In_FEM) |\n" \
94  << "| |\n" \
95  << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
96  << "| 2) Basis values for VALUE and CURL operators |\n" \
97  << "| |\n" \
98  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
99  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
100  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
101  << "| |\n" \
102  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
103  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
104  << "| |\n" \
105  << "===============================================================================\n"\
106  << "| TEST 1: Basis creation, exception testing |\n"\
107  << "===============================================================================\n";
108 
109  // Define basis and error flag
110  const int deg = 1;
111  shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
112 
113  Basis_HCURL_QUAD_In_FEM<double, FieldContainer<double> > quadBasis(deg,POINTTYPE_EQUISPACED);
114  int errorFlag = 0;
115 
116  // Initialize throw counter for exception testing
117  int nException = 0;
118  int throwCounter = 0;
119 
120  // Array with the 4 vertices of the reference Quadrilateral, its center and 4 more points
121  FieldContainer<double> quadNodes(9, 2);
122  quadNodes(0,0) = -1.0; quadNodes(0,1) = -1.0;
123  quadNodes(1,0) = 0.0; quadNodes(1,1) = -1.0;
124  quadNodes(2,0) = 1.0; quadNodes(2,1) = -1.0;
125  quadNodes(3,0) = -1.0; quadNodes(3,1) = 0.0;
126  quadNodes(4,0) = 0.0; quadNodes(4,1) = 0.0;
127  quadNodes(5,0) = 1.0; quadNodes(5,1) = 0.0;
128  quadNodes(6,0) = -1.0; quadNodes(6,1) = 1.0;
129  quadNodes(7,0) = 0.0; quadNodes(7,1) = 1.0;
130  quadNodes(8,0) = 1.0; quadNodes(8,1) = 1.0;
131 
132 
133  // Generic array for the output values; needs to be properly resized depending on the operator type
135 
136  try{
137  // exception #1: GRAD cannot be applied to HCURL functions
138  // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
139  vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0), 4 );
140  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD), throwCounter, nException );
141 
142  // exception #2: DIV cannot be applied to HCURL functions
143  // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
144  vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0) );
145  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException );
146 
147  // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
148  // getDofTag() to access invalid array elements thereby causing bounds check exception
149  // exception #3
150  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException );
151  // exception #4
152  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException );
153  // exception #5
154  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,1), throwCounter, nException );
155  // exception #6
156  INTREPID_TEST_COMMAND( quadBasis.getDofTag(12), throwCounter, nException );
157  // exception #7
158  INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException );
159 
160 #ifdef HAVE_INTREPID_DEBUG
161  // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays
162  // exception #8: input points array must be of rank-2
163  FieldContainer<double> badPoints1(4, 5, 3);
164  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
165 
166  // exception #9 dimension 1 in the input point array must equal space dimension of the cell
167  FieldContainer<double> badPoints2(4, quadBasis.getBaseCellTopology().getDimension() + 1);
168  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
169 
170  // exception #10 output values must be of rank-3 for OPERATOR_VALUE in 2D
171  FieldContainer<double> badVals1(4, 3);
172  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
173 
174  FieldContainer<double> badCurls1(4,3,2);
175  // exception #11 output values must be of rank-2 for OPERATOR_CURL
176  INTREPID_TEST_COMMAND( quadBasis.getValues(badCurls1, quadNodes, OPERATOR_CURL), throwCounter, nException );
177 
178  // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
179  FieldContainer<double> badVals2(quadBasis.getCardinality() + 1, quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension());
180  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
181 
182  // exception #13 incorrect 1st dimension of output array (must equal number of points)
183  FieldContainer<double> badVals3(quadBasis.getCardinality(), quadNodes.dimension(0) + 1, quadBasis.getBaseCellTopology().getDimension() );
184  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
185 
186  // exception #14: incorrect 2nd dimension of output array for VALUE (must equal the space dimension)
187  FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() - 1);
188  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
189 
190  // exception #15: D2 cannot be applied to HCURL functions
191  // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
192  vals.resize(quadBasis.getCardinality(),
193  quadNodes.dimension(0),
194  Intrepid::getDkCardinality(OPERATOR_D2, quadBasis.getBaseCellTopology().getDimension()));
195  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_D2), throwCounter, nException );
196 #endif
197 
198  }
199  catch (const std::logic_error & err) {
200  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
201  *outStream << err.what() << '\n';
202  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
203  errorFlag = -1000;
204  };
205 
206  // Check if number of thrown exceptions matches the one we expect
207  // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match.
208  if (throwCounter != nException) {
209  errorFlag++;
210  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
211  }
212 //#endif
213 
214  *outStream \
215  << "\n"
216  << "===============================================================================\n"\
217  << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
218  << "===============================================================================\n";
219 
220  try{
221  std::vector<std::vector<int> > allTags = quadBasis.getAllDofTags();
222 
223  // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
224  for (unsigned i = 0; i < allTags.size(); i++) {
225  int bfOrd = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
226 
227  std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
228  if( !( (myTag[0] == allTags[i][0]) &&
229  (myTag[1] == allTags[i][1]) &&
230  (myTag[2] == allTags[i][2]) &&
231  (myTag[3] == allTags[i][3]) ) ) {
232  errorFlag++;
233  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
234  *outStream << " getDofOrdinal( {"
235  << allTags[i][0] << ", "
236  << allTags[i][1] << ", "
237  << allTags[i][2] << ", "
238  << allTags[i][3] << "}) = " << bfOrd <<" but \n";
239  *outStream << " getDofTag(" << bfOrd << ") = { "
240  << myTag[0] << ", "
241  << myTag[1] << ", "
242  << myTag[2] << ", "
243  << myTag[3] << "}\n";
244  }
245  }
246 
247  // Now do the same but loop over basis functions
248  for( int bfOrd = 0; bfOrd < quadBasis.getCardinality(); bfOrd++) {
249  std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
250  int myBfOrd = quadBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
251  if( bfOrd != myBfOrd) {
252  errorFlag++;
253  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
254  *outStream << " getDofTag(" << bfOrd << ") = { "
255  << myTag[0] << ", "
256  << myTag[1] << ", "
257  << myTag[2] << ", "
258  << myTag[3] << "} but getDofOrdinal({"
259  << myTag[0] << ", "
260  << myTag[1] << ", "
261  << myTag[2] << ", "
262  << myTag[3] << "} ) = " << myBfOrd << "\n";
263  }
264  }
265  }
266  catch (const std::logic_error & err){
267  *outStream << err.what() << "\n\n";
268  errorFlag = -1000;
269  };
270 
271  *outStream \
272  << "\n"
273  << "===============================================================================\n"\
274  << "| TEST 3: correctness of basis function values |\n"\
275  << "===============================================================================\n";
276 
277  outStream -> precision(20);
278 
279  // VALUE: Each row pair gives the 4x2 correct basis set values at an evaluation point: (P,F,D) layout
280  double basisValues[] = {
281  // first bf, first row of points
282  1,0,1,0,1,0,
283  // first bf, second row of points
284  0.5,0,0.5,0,0.5,0,
285  // first bf, third row of points
286  0,0,0,0,0,0,
287  // second bf, first row of points
288  0,0,0,0,0,0,
289  // second bf, second row of points
290  0.5,0,0.5,0,0.5,0,
291  // second bf, third row of points
292  1,0,1,0,1,0,
293  // third bf, first row of points
294  0,1,0,0.5,0,0,
295  // third bf, second row of points
296  0,1,0,0.5,0,0,
297  // third bf, third row of points
298  0,1,0,0.5,0,0,
299  // fourth bf, first row of points
300  0,0,0,0.5,0,1,
301  // fourth bf, second row of points
302  0,0,0,0.5,0,1,
303  // fourth bf, third row of points
304  0,0,0,0.5,0,1,
305  };
306 
307  // CURL: correct values in (F,P) format
308  double basisCurls[] = {
309  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
310  -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
311  -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
312  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5
313  };
314 
315 
316  try{
317  // Dimensions for the output arrays:
318  int numFields = quadBasis.getCardinality();
319  int numPoints = quadNodes.dimension(0);
320  int spaceDim = quadBasis.getBaseCellTopology().getDimension();
321 
322  // Generic array for values and curls that will be properly sized before each call
324 
325  // Check VALUE of basis functions: resize vals to rank-3 container:
326  vals.resize(numFields, numPoints, spaceDim);
327  quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE);
328  for (int i = 0; i < numFields; i++) {
329  for (int j = 0; j < numPoints; j++) {
330  for (int k = 0; k < spaceDim; k++) {
331 
332  // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k
333  int l = k + i * spaceDim * numPoints + j * spaceDim;
334  if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
335  errorFlag++;
336  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
337 
338  // Output the multi-index of the value where the error is:
339  *outStream << " At multi-index { ";
340  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
341  *outStream << "} computed value: " << vals(i,j,k)
342  << " but reference value: " << basisValues[l] << "\n";
343  }
344  }
345  }
346  }
347 
348  // Check CURL of basis function: resize vals to rank-2 container
349  vals.resize(numFields, numPoints);
350  quadBasis.getValues(vals, quadNodes, OPERATOR_CURL);
351  for (int i = 0; i < numFields; i++) {
352  for (int j = 0; j < numPoints; j++) {
353  int l = j + i * numPoints;
354  if (std::abs(vals(i,j) - basisCurls[l]) > INTREPID_TOL) {
355  errorFlag++;
356  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
357 
358  // Output the multi-index of the value where the error is:
359  *outStream << " At multi-index { ";
360  *outStream << i << " ";*outStream << j << " ";
361  *outStream << "} computed curl component: " << vals(i,j)
362  << " but reference curl component: " << basisCurls[l] << "\n";
363  }
364  }
365  }
366 
367  }
368 
369  // Catch unexpected errors
370  catch (const std::logic_error & err) {
371  *outStream << err.what() << "\n\n";
372  errorFlag = -1000;
373  };
374 
375  if (errorFlag != 0)
376  std::cout << "End Result: TEST FAILED\n";
377  else
378  std::cout << "End Result: TEST PASSED\n";
379 
380  // reset format state of std::cout
381  std::cout.copyfmt(oldFormatState);
382 
383  return errorFlag;
384 }
int main(int argc, char *argv[])
Performs a code-code comparison to FIAT for Nedelec bases on tets (values and curls) ...
Definition: test_01.cpp:65
Header file for utility class to provide point tools, such as barycentric coordinates, equispaced lattices, and warp-blend point distrubtions.
Header file for utility class to provide multidimensional containers.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
Header file for the Intrepid::HCURL_QUAD_In_FEM class.
Implementation of the default H(div)-compatible FEM basis of degree 1 on Quadrilateral cell...