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 
131  try{
132  // Generic array for the output values; needs to be properly resized depending on the operator type
134 
135  // exception #1: GRAD cannot be applied to HDIV functions
136  // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
137  vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0), 3 );
138  INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_GRAD), throwCounter, nException );
139 
140  // exception #2: CURL cannot be applied to HDIV functions
141  INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_CURL), 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( tetBasis.getDofOrdinal(3,0,0), throwCounter, nException );
147  // exception #4
148  INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(1,1,1), throwCounter, nException );
149  // exception #5
150  INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(0,4,1), throwCounter, nException );
151  // exception #6
152  INTREPID_TEST_COMMAND( tetBasis.getDofTag(7), throwCounter, nException );
153  // exception #7
154  INTREPID_TEST_COMMAND( tetBasis.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( tetBasis.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, 2);
164  INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
165 
166  // exception #10 output values must be of rank-3 for OPERATOR_VALUE
167  FieldContainer<double> badVals1(4, 3);
168  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
169 
170  // exception #11 output values must be of rank-2 for OPERATOR_DIV
171  FieldContainer<double> badVals2(4, 3, 1);
172  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_VALUE), throwCounter, nException );
173 
174  // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
175  FieldContainer<double> badVals3(tetBasis.getCardinality() + 1, tetNodes.dimension(0), 3);
176  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
177 
178  // exception #13 incorrect 0th dimension of output array (must equal number of basis functions)
179  FieldContainer<double> badVals4(tetBasis.getCardinality() + 1, tetNodes.dimension(0));
180  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals4, tetNodes, OPERATOR_DIV), throwCounter, nException );
181 
182  // exception #14 incorrect 1st dimension of output array (must equal number of points)
183  FieldContainer<double> badVals5(tetBasis.getCardinality(), tetNodes.dimension(0) + 1, 3);
184  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals5, tetNodes, OPERATOR_VALUE), throwCounter, nException );
185 
186  // exception #15 incorrect 1st dimension of output array (must equal number of points)
187  FieldContainer<double> badVals6(tetBasis.getCardinality(), tetNodes.dimension(0) + 1);
188  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_DIV), throwCounter, nException );
189 
190  // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
191  FieldContainer<double> badVals7(tetBasis.getCardinality(), tetNodes.dimension(0), 4);
192  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals7, tetNodes, OPERATOR_VALUE), throwCounter, nException );
193 #endif
194 
195  }
196  catch (const std::logic_error & err) {
197  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
198  *outStream << err.what() << '\n';
199  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
200  errorFlag = -1000;
201  };
202 
203  // Check if number of thrown exceptions matches the one we expect
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 = tetBasis.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 = tetBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
221 
222  std::vector<int> myTag = tetBasis.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 < tetBasis.getCardinality(); bfOrd++) {
244  std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
245  int myBfOrd = tetBasis.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 basis values in (P,F,D) format: each row gives the 4x3 correct basis set values
275  // at an evaluation point. Note that getValues returns results as an (F,P,D) array.
276  double basisValues[] = {
277  // 4 vertices
278  0.,-2.0,0., 0.,0.,0., -2.0,0.,0., 0.,0.,-2.0,
279  2.0,-2.0,0., 2.0,0.,0., 0.,0.,0., 2.0,0.,-2.0,
280  0.,0.,0., 0.,2.0,0., -2.0,2.0,0., 0,2.0,-2.0,
281  0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, 0.,0.,0.,
282  // 6 edge midpoints
283  1.0,-2.0,0., 1.0,0.,0., -1.0,0.,0., 1.0,0.,-2.0,
284  1.0,-1.0,0., 1.0,1.0,0., -1.0,1.0,0., 1.0,1.0,-2.0,
285  0.,-1.0,0., 0.,1.0,0., -2.0,1.0,0., 0.,1.0,-2.0,
286  0.,-2.0,1.0, 0.,0.,1.0, -2.0,0.,1.0, 0.,0.,-1.0,
287  1.0,-2.0,1.0, 1.0,0.,1.0, -1.0,0.,1.0, 1.0,0.,-1.0,
288  0.,-1.0,1.0, 0.,1.0,1.0, -2.0,1.0,1.0, 0.,1.0,-1.0
289  // bf0 bf1 bf2 bf3
290  };
291 
292  // DIV: each row gives the 4 correct values of the divergence of the 4 basis functions
293  double basisDivs[] = {
294  // 4 vertices
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.0, 6.0, 6.0, 6.0,
299  // 6 edge midpoints
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  6.0, 6.0, 6.0, 6.0
306  };
307 
308  try{
309 
310  // Dimensions for the output arrays:
311  int numFields = tetBasis.getCardinality();
312  int numPoints = tetNodes.dimension(0);
313  int spaceDim = tetBasis.getBaseCellTopology().getDimension();
314 
315  // Generic array for values and curls that will be properly sized before each call
317 
318  // Check VALUE of basis functions: resize vals to rank-3 container:
319  vals.resize(numFields, numPoints, spaceDim);
320  tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE);
321  for (int i = 0; i < numFields; i++) {
322  for (int j = 0; j < numPoints; j++) {
323  for (int k = 0; k < spaceDim; k++) {
324  // basisValues is (P,F,D) array so its multiindex is (j,i,k) and not (i,j,k)!
325  int l = k + i * spaceDim + j * spaceDim * numFields;
326  if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
327  errorFlag++;
328  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
329 
330  // Output the multi-index of the value where the error is:
331  *outStream << " At (Field,Point,Dim) multi-index { ";
332  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
333  *outStream << "} computed value: " << vals(i,j,k)
334  << " but reference value: " << basisValues[l] << "\n";
335  }
336  }
337  }
338  }
339 
340 
341  // Check DIV of basis function: resize vals to rank-2 container
342  vals.resize(numFields, numPoints);
343  tetBasis.getValues(vals, tetNodes, OPERATOR_DIV);
344  for (int i = 0; i < numFields; i++) {
345  for (int j = 0; j < numPoints; j++) {
346  int l = i + j * numFields;
347  if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) {
348  errorFlag++;
349  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
350 
351  // Output the multi-index of the value where the error is:
352  *outStream << " At multi-index { ";
353  *outStream << i << " ";*outStream << j << " ";
354  *outStream << "} computed divergence component: " << vals(i,j)
355  << " but reference divergence component: " << basisDivs[l] << "\n";
356  }
357  }
358  }
359 
360  }
361 
362  // Catch unexpected errors
363  catch (const std::logic_error & err) {
364  *outStream << err.what() << "\n\n";
365  errorFlag = -1000;
366  };
367 
368  *outStream \
369  << "\n"
370  << "===============================================================================\n"\
371  << "| TEST 4: correctness of DoF locations |\n"\
372  << "===============================================================================\n";
373 
374  try{
375  Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
376  Teuchos::rcp(new Basis_HDIV_TET_I1_FEM<double, FieldContainer<double> >);
377  Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
378  Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
379 
380  int spaceDim = 3;
382  FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality(),spaceDim); // last dimension is spatial dim
383 
384  // Check exceptions.
385 #ifdef HAVE_INTREPID_DEBUG
386  cvals.resize(1,2,3);
387  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
388  cvals.resize(3,2);
389  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
390  cvals.resize(4,2);
391  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
392 #endif
393  cvals.resize(4,spaceDim);
394  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
395  // Check if number of thrown exceptions matches the one we expect
396  if (throwCounter != nException) {
397  errorFlag++;
398  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
399  }
400 
401  // Check mathematical correctness
402  FieldContainer<double> normals(basis->getCardinality(),spaceDim); // normals at each point basis point
403  normals(0,0) = 0.0; normals(0,1) = -0.5; normals(0,2) = 0.0;
404  normals(1,0) = 0.5; normals(1,1) = 0.5; normals(1,2) = 0.5;
405  normals(2,0) = -0.5; normals(2,1) = 0.0; normals(2,2) = 0.0;
406  normals(3,0) = 0.0; normals(3,1) = 0.0; normals(3,2) = -0.5;
407 
408  basis->getValues(bvals, cvals, OPERATOR_VALUE);
409  char buffer[120];
410  for (int i=0; i<bvals.dimension(0); i++) {
411  for (int j=0; j<bvals.dimension(1); j++) {
412 
413  double normal = 0.0;
414  for(int d=0;d<spaceDim;d++)
415  normal += bvals(i,j,d)*normals(j,d);
416 
417  if ((i != j) && (std::abs(normal - 0.0) > INTREPID_TOL)) {
418  errorFlag++;
419  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);
420  *outStream << buffer;
421  }
422  else if ((i == j) && (std::abs(normal - 1.0) > INTREPID_TOL)) {
423  errorFlag++;
424  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);
425  *outStream << buffer;
426  }
427  }
428  }
429 
430  }
431  catch (const std::logic_error & err){
432  *outStream << err.what() << "\n\n";
433  errorFlag = -1000;
434  };
435 
436  if (errorFlag != 0)
437  std::cout << "End Result: TEST FAILED\n";
438  else
439  std::cout << "End Result: TEST PASSED\n";
440 
441  // reset format state of std::cout
442  std::cout.copyfmt(oldFormatState);
443 
444  return errorFlag;
445 }
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...