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_HDIV_TET_I1_FEM) |\n" \
93  << "| |\n" \
94  << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95  << "| 3) Exception tests on input arguments and invalid operators |\n" \
96  << "| 2) Basis values for VALUE and DIV 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
111  int errorFlag = 0;
112 
113  // Initialize throw counter for exception testing
114  int nException = 0;
115  int throwCounter = 0;
116 
117  // Define array containing the 4 vertices of the reference TET and its 6 edge midpoints.
118  FieldContainer<double> tetNodes(10, 3);
119  tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0;
120  tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0;
121  tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0;
122  tetNodes(3,0) = 0.0; tetNodes(3,1) = 0.0; tetNodes(3,2) = 1.0;
123  tetNodes(4,0) = 0.5; tetNodes(4,1) = 0.0; tetNodes(4,2) = 0.0;
124  tetNodes(5,0) = 0.5; tetNodes(5,1) = 0.5; tetNodes(5,2) = 0.0;
125  tetNodes(6,0) = 0.0; tetNodes(6,1) = 0.5; tetNodes(6,2) = 0.0;
126  tetNodes(7,0) = 0.0; tetNodes(7,1) = 0.0; tetNodes(7,2) = 0.5;
127  tetNodes(8,0) = 0.5; tetNodes(8,1) = 0.0; tetNodes(8,2) = 0.5;
128  tetNodes(9,0) = 0.0; tetNodes(9,1) = 0.5; tetNodes(9,2) = 0.5;
129 
130  // Generic array for the output values; needs to be properly resized depending on the operator type
132 
133  try{
134  // exception #1: GRAD cannot be applied to HDIV functions
135  // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
136  vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0), 3 );
137  INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_GRAD), throwCounter, nException );
138 
139  // exception #2: CURL cannot be applied to HDIV functions
140  INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException );
141 
142  // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
143  // getDofTag() to access invalid array elements thereby causing bounds check exception
144  // exception #3
145  INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(3,0,0), throwCounter, nException );
146  // exception #4
147  INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(1,1,1), throwCounter, nException );
148  // exception #5
149  INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(0,4,1), throwCounter, nException );
150  // exception #6
151  INTREPID_TEST_COMMAND( tetBasis.getDofTag(7), throwCounter, nException );
152  // exception #7
153  INTREPID_TEST_COMMAND( tetBasis.getDofTag(-1), throwCounter, nException );
154 
155 #ifdef HAVE_INTREPID_DEBUG
156  // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays
157  // exception #8: input points array must be of rank-2
158  FieldContainer<double> badPoints1(4, 5, 3);
159  INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
160 
161  // exception #9 dimension 1 in the input point array must equal space dimension of the cell
162  FieldContainer<double> badPoints2(4, 2);
163  INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
164 
165  // exception #10 output values must be of rank-3 for OPERATOR_VALUE
166  FieldContainer<double> badVals1(4, 3);
167  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
168 
169  // exception #11 output values must be of rank-2 for OPERATOR_DIV
170  FieldContainer<double> badVals2(4, 3, 1);
171  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_VALUE), throwCounter, nException );
172 
173  // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
174  FieldContainer<double> badVals3(tetBasis.getCardinality() + 1, tetNodes.dimension(0), 3);
175  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
176 
177  // exception #13 incorrect 0th dimension of output array (must equal number of basis functions)
178  FieldContainer<double> badVals4(tetBasis.getCardinality() + 1, tetNodes.dimension(0));
179  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals4, tetNodes, OPERATOR_DIV), throwCounter, nException );
180 
181  // exception #14 incorrect 1st dimension of output array (must equal number of points)
182  FieldContainer<double> badVals5(tetBasis.getCardinality(), tetNodes.dimension(0) + 1, 3);
183  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals5, tetNodes, OPERATOR_VALUE), throwCounter, nException );
184 
185  // exception #15 incorrect 1st dimension of output array (must equal number of points)
186  FieldContainer<double> badVals6(tetBasis.getCardinality(), tetNodes.dimension(0) + 1);
187  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_DIV), throwCounter, nException );
188 
189  // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
190  FieldContainer<double> badVals7(tetBasis.getCardinality(), tetNodes.dimension(0), 4);
191  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals7, tetNodes, OPERATOR_VALUE), 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  if (throwCounter != nException) {
204  errorFlag++;
205  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
206  }
207 
208  *outStream \
209  << "\n"
210  << "===============================================================================\n"\
211  << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
212  << "===============================================================================\n";
213 
214  try{
215  std::vector<std::vector<int> > allTags = tetBasis.getAllDofTags();
216 
217  // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
218  for (unsigned i = 0; i < allTags.size(); i++) {
219  int bfOrd = tetBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
220 
221  std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
222  if( !( (myTag[0] == allTags[i][0]) &&
223  (myTag[1] == allTags[i][1]) &&
224  (myTag[2] == allTags[i][2]) &&
225  (myTag[3] == allTags[i][3]) ) ) {
226  errorFlag++;
227  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
228  *outStream << " getDofOrdinal( {"
229  << allTags[i][0] << ", "
230  << allTags[i][1] << ", "
231  << allTags[i][2] << ", "
232  << allTags[i][3] << "}) = " << bfOrd <<" but \n";
233  *outStream << " getDofTag(" << bfOrd << ") = { "
234  << myTag[0] << ", "
235  << myTag[1] << ", "
236  << myTag[2] << ", "
237  << myTag[3] << "}\n";
238  }
239  }
240 
241  // Now do the same but loop over basis functions
242  for( int bfOrd = 0; bfOrd < tetBasis.getCardinality(); bfOrd++) {
243  std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
244  int myBfOrd = tetBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
245  if( bfOrd != myBfOrd) {
246  errorFlag++;
247  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
248  *outStream << " getDofTag(" << bfOrd << ") = { "
249  << myTag[0] << ", "
250  << myTag[1] << ", "
251  << myTag[2] << ", "
252  << myTag[3] << "} but getDofOrdinal({"
253  << myTag[0] << ", "
254  << myTag[1] << ", "
255  << myTag[2] << ", "
256  << myTag[3] << "} ) = " << myBfOrd << "\n";
257  }
258  }
259  }
260  catch (const std::logic_error & err){
261  *outStream << err.what() << "\n\n";
262  errorFlag = -1000;
263  };
264 
265  *outStream \
266  << "\n"
267  << "===============================================================================\n"\
268  << "| TEST 3: correctness of basis function values |\n"\
269  << "===============================================================================\n";
270 
271  outStream -> precision(20);
272 
273  // VALUE: Correct basis values in (P,F,D) format: each row gives the 4x3 correct basis set values
274  // at an evaluation point. Note that getValues returns results as an (F,P,D) array.
275  double basisValues[] = {
276  // 4 vertices
277  0.,-2.0,0., 0.,0.,0., -2.0,0.,0., 0.,0.,-2.0,
278  2.0,-2.0,0., 2.0,0.,0., 0.,0.,0., 2.0,0.,-2.0,
279  0.,0.,0., 0.,2.0,0., -2.0,2.0,0., 0,2.0,-2.0,
280  0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, 0.,0.,0.,
281  // 6 edge midpoints
282  1.0,-2.0,0., 1.0,0.,0., -1.0,0.,0., 1.0,0.,-2.0,
283  1.0,-1.0,0., 1.0,1.0,0., -1.0,1.0,0., 1.0,1.0,-2.0,
284  0.,-1.0,0., 0.,1.0,0., -2.0,1.0,0., 0.,1.0,-2.0,
285  0.,-2.0,1.0, 0.,0.,1.0, -2.0,0.,1.0, 0.,0.,-1.0,
286  1.0,-2.0,1.0, 1.0,0.,1.0, -1.0,0.,1.0, 1.0,0.,-1.0,
287  0.,-1.0,1.0, 0.,1.0,1.0, -2.0,1.0,1.0, 0.,1.0,-1.0
288  // bf0 bf1 bf2 bf3
289  };
290 
291  // DIV: each row gives the 4 correct values of the divergence of the 4 basis functions
292  double basisDivs[] = {
293  // 4 vertices
294  6.0, 6.0, 6.0, 6.0,
295  6.0, 6.0, 6.0, 6.0,
296  6.0, 6.0, 6.0, 6.0,
297  6.0, 6.0, 6.0, 6.0,
298  // 6 edge midpoints
299  6.0, 6.0, 6.0, 6.0,
300  6.0, 6.0, 6.0, 6.0,
301  6.0, 6.0, 6.0, 6.0,
302  6.0, 6.0, 6.0, 6.0,
303  6.0, 6.0, 6.0, 6.0,
304  6.0, 6.0, 6.0, 6.0
305  };
306 
307  try{
308 
309  // Dimensions for the output arrays:
310  int numFields = tetBasis.getCardinality();
311  int numPoints = tetNodes.dimension(0);
312  int spaceDim = tetBasis.getBaseCellTopology().getDimension();
313 
314  // Generic array for values and curls that will be properly sized before each call
316 
317  // Check VALUE of basis functions: resize vals to rank-3 container:
318  vals.resize(numFields, numPoints, spaceDim);
319  tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE);
320  for (int i = 0; i < numFields; i++) {
321  for (int j = 0; j < numPoints; j++) {
322  for (int k = 0; k < spaceDim; k++) {
323  // basisValues is (P,F,D) array so its multiindex is (j,i,k) and not (i,j,k)!
324  int l = k + i * spaceDim + j * spaceDim * numFields;
325  if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
326  errorFlag++;
327  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
328 
329  // Output the multi-index of the value where the error is:
330  *outStream << " At (Field,Point,Dim) multi-index { ";
331  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
332  *outStream << "} computed value: " << vals(i,j,k)
333  << " but reference value: " << basisValues[l] << "\n";
334  }
335  }
336  }
337  }
338 
339 
340  // Check DIV of basis function: resize vals to rank-2 container
341  vals.resize(numFields, numPoints);
342  tetBasis.getValues(vals, tetNodes, OPERATOR_DIV);
343  for (int i = 0; i < numFields; i++) {
344  for (int j = 0; j < numPoints; j++) {
345  int l = i + j * numFields;
346  if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) {
347  errorFlag++;
348  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
349 
350  // Output the multi-index of the value where the error is:
351  *outStream << " At multi-index { ";
352  *outStream << i << " ";*outStream << j << " ";
353  *outStream << "} computed divergence component: " << vals(i,j)
354  << " but reference divergence component: " << basisDivs[l] << "\n";
355  }
356  }
357  }
358 
359  }
360 
361  // Catch unexpected errors
362  catch (const std::logic_error & err) {
363  *outStream << err.what() << "\n\n";
364  errorFlag = -1000;
365  };
366 
367  *outStream \
368  << "\n"
369  << "===============================================================================\n"\
370  << "| TEST 4: correctness of DoF locations |\n"\
371  << "===============================================================================\n";
372 
373  try{
374  Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
375  Teuchos::rcp(new Basis_HDIV_TET_I1_FEM<double, FieldContainer<double> >);
376  Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
377  Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
378 
379  int spaceDim = 3;
381  FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality(),spaceDim); // last dimension is spatial dim
382 
383  // Check exceptions.
384 #ifdef HAVE_INTREPID_DEBUG
385  cvals.resize(1,2,3);
386  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
387  cvals.resize(3,2);
388  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
389  cvals.resize(4,2);
390  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
391 #endif
392  cvals.resize(4,spaceDim);
393  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
394  // Check if number of thrown exceptions matches the one we expect
395  if (throwCounter != nException) {
396  errorFlag++;
397  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
398  }
399 
400  // Check mathematical correctness
401  FieldContainer<double> normals(basis->getCardinality(),spaceDim); // normals at each point basis point
402  normals(0,0) = 0.0; normals(0,1) = -0.5; normals(0,2) = 0.0;
403  normals(1,0) = 0.5; normals(1,1) = 0.5; normals(1,2) = 0.5;
404  normals(2,0) = -0.5; normals(2,1) = 0.0; normals(2,2) = 0.0;
405  normals(3,0) = 0.0; normals(3,1) = 0.0; normals(3,2) = -0.5;
406 
407  basis->getValues(bvals, cvals, OPERATOR_VALUE);
408  char buffer[120];
409  for (int i=0; i<bvals.dimension(0); i++) {
410  for (int j=0; j<bvals.dimension(1); j++) {
411 
412  double normal = 0.0;
413  for(int d=0;d<spaceDim;d++)
414  normal += bvals(i,j,d)*normals(j,d);
415 
416  if ((i != j) && (std::abs(normal - 0.0) > INTREPID_TOL)) {
417  errorFlag++;
418  sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), normal, 0.0);
419  *outStream << buffer;
420  }
421  else if ((i == j) && (std::abs(normal - 1.0) > INTREPID_TOL)) {
422  errorFlag++;
423  sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), normal, 1.0);
424  *outStream << buffer;
425  }
426  }
427  }
428 
429  }
430  catch (const std::logic_error & err){
431  *outStream << err.what() << "\n\n";
432  errorFlag = -1000;
433  };
434 
435  if (errorFlag != 0)
436  std::cout << "End Result: TEST FAILED\n";
437  else
438  std::cout << "End Result: TEST PASSED\n";
439 
440  // reset format state of std::cout
441  std::cout.copyfmt(oldFormatState);
442 
443  return errorFlag;
444 }
Implementation of the default H(div)-compatible FEM basis of degree 1 on Tetrahedron 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...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Tetrahedron cell.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HDIV_TET_I1_FEM class.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
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...