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