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 
48 #include "Intrepid_FieldContainer.hpp"
49 #include "Intrepid_HGRAD_TET_C1_FEM.hpp"
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_HGRAD_TET_C1_FEM) |\n" \
93  << "| |\n" \
94  << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95  << "| 2) Basis values for VALUE, GRAD, CURL, and Dk 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
109  Basis_HGRAD_TET_C1_FEM<double, FieldContainer<double> > tetBasis;
110  int errorFlag = 0;
111 
112  // Initialize throw counter for exception testing
113  int nException = 0;
114  int throwCounter = 0;
115 
116  // Define array containing the 4 vertices of the reference TET and its center.
117  FieldContainer<double> tetNodes(8, 3);
118  tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0;
119  tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0;
120  tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0;
121  tetNodes(3,0) = 0.0; tetNodes(3,1) = 0.0; tetNodes(3,2) = 1.0;
122  tetNodes(4,0) = 0.25; tetNodes(4,1) = 0.5; tetNodes(4,2) = 0.1;
123  tetNodes(5,0) = 0.5; tetNodes(5,1) = 0.5; tetNodes(5,2) = 0.0;
124  tetNodes(6,0) = 0.5; tetNodes(6,1) = 0.0; tetNodes(6,2) = 0.5;
125  tetNodes(7,0) = 0.0; tetNodes(7,1) = 0.5; tetNodes(7,2) = 0.5;
126 
127 
128  // Generic array for the output values; needs to be properly resized depending on the operator type
129  FieldContainer<double> vals;
130 
131  try{
132  // exception #1: CURL cannot be applied to scalar functions
133  // resize vals to rank-3 container with dimensions (num. points, num. basis functions, arbitrary)
134  vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0), 4 );
135  INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException );
136 
137  // exception #2: DIV cannot be applied to scalar functions
138  // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
139  vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0) );
140  INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_DIV), 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,0), throwCounter, nException );
150  // exception #6
151  INTREPID_TEST_COMMAND( tetBasis.getDofTag(5), throwCounter, nException );
152  // exception #7
153  INTREPID_TEST_COMMAND( tetBasis.getDofTag(-1), throwCounter, nException );
154 
155 #ifdef HAVE_INTREPID_DEBUG
156  // Exceptions 8-18 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, tetBasis.getBaseCellTopology().getDimension() - 1);
163  INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
164 
165  // exception #10 output values must be of rank-2 for OPERATOR_VALUE
166  FieldContainer<double> badVals1(4, 3, 1);
167  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
168 
169  // exception #11 output values must be of rank-3 for OPERATOR_GRAD
170  FieldContainer<double> badVals2(4, 3);
171  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_GRAD), throwCounter, nException );
172 
173  // exception #12 output values must be of rank-3 for OPERATOR_D1
174  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_D1), throwCounter, nException );
175 
176  // exception #13 output values must be of rank-3 for OPERATOR_D2
177  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_D2), throwCounter, nException );
178 
179  // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
180  FieldContainer<double> badVals3(tetBasis.getCardinality() + 1, tetNodes.dimension(0));
181  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
182 
183  // exception #15 incorrect 1st dimension of output array (must equal number of points)
184  FieldContainer<double> badVals4(tetBasis.getCardinality(), tetNodes.dimension(0) + 1);
185  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals4, tetNodes, OPERATOR_VALUE), throwCounter, nException );
186 
187  // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
188  FieldContainer<double> badVals5(tetBasis.getCardinality(), tetNodes.dimension(0), tetBasis.getBaseCellTopology().getDimension() + 1);
189  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals5, tetNodes, OPERATOR_GRAD), throwCounter, nException );
190 
191  // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D)
192  FieldContainer<double> badVals6(tetBasis.getCardinality(), tetNodes.dimension(0), 40);
193  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_D1), throwCounter, nException );
194 
195  // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D)
196  INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_D2), throwCounter, nException );
197 #endif
198 
199  }
200  catch (const std::logic_error & err) {
201  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
202  *outStream << err.what() << '\n';
203  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
204  errorFlag = -1000;
205  };
206 
207  // Check if number of thrown exceptions matches the one we expect
208  if (throwCounter != nException) {
209  errorFlag++;
210  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
211  }
212 
213  *outStream \
214  << "\n"
215  << "===============================================================================\n"\
216  << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
217  << "===============================================================================\n";
218 
219  try{
220  std::vector<std::vector<int> > allTags = tetBasis.getAllDofTags();
221 
222  // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
223  for (unsigned i = 0; i < allTags.size(); i++) {
224  int bfOrd = tetBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
225 
226  std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
227  if( !( (myTag[0] == allTags[i][0]) &&
228  (myTag[1] == allTags[i][1]) &&
229  (myTag[2] == allTags[i][2]) &&
230  (myTag[3] == allTags[i][3]) ) ) {
231  errorFlag++;
232  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
233  *outStream << " getDofOrdinal( {"
234  << allTags[i][0] << ", "
235  << allTags[i][1] << ", "
236  << allTags[i][2] << ", "
237  << allTags[i][3] << "}) = " << bfOrd <<" but \n";
238  *outStream << " getDofTag(" << bfOrd << ") = { "
239  << myTag[0] << ", "
240  << myTag[1] << ", "
241  << myTag[2] << ", "
242  << myTag[3] << "}\n";
243  }
244  }
245 
246  // Now do the same but loop over basis functions
247  for( int bfOrd = 0; bfOrd < tetBasis.getCardinality(); bfOrd++) {
248  std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
249  int myBfOrd = tetBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
250  if( bfOrd != myBfOrd) {
251  errorFlag++;
252  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
253  *outStream << " getDofTag(" << bfOrd << ") = { "
254  << myTag[0] << ", "
255  << myTag[1] << ", "
256  << myTag[2] << ", "
257  << myTag[3] << "} but getDofOrdinal({"
258  << myTag[0] << ", "
259  << myTag[1] << ", "
260  << myTag[2] << ", "
261  << myTag[3] << "} ) = " << myBfOrd << "\n";
262  }
263  }
264  }
265  catch (const std::logic_error & err){
266  *outStream << err.what() << "\n\n";
267  errorFlag = -1000;
268  };
269 
270  *outStream \
271  << "\n"
272  << "===============================================================================\n"\
273  << "| TEST 3: correctness of basis function values |\n"\
274  << "===============================================================================\n";
275 
276  outStream -> precision(20);
277 
278  // VALUE: Each row gives the 4 correct basis set values at an evaluation point
279  double basisValues[] = {
280  1.0, 0.0, 0.0, 0.0,
281  0.0, 1.0, 0.0, 0.0,
282  0.0, 0.0, 1.0, 0.0,
283  0.0, 0.0, 0.0, 1.0,
284  0.15,0.25,0.5, 0.1,
285  0.0, 0.5, 0.5, 0.0,
286  0.0, 0.5, 0.0, 0.5,
287  0.0, 0.0, 0.5, 0.5
288  };
289 
290  // GRAD and D1: each row gives the 3 x 4 correct values of the gradients of the 4 basis functions
291  double basisGrads[] = {
292  -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
293  -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
294  -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
295  -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
296  -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
297  -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
298  -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
299  -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
300  };
301 
302  try{
303 
304  // Dimensions for the output arrays:
305  int numFields = tetBasis.getCardinality();
306  int numPoints = tetNodes.dimension(0);
307  int spaceDim = tetBasis.getBaseCellTopology().getDimension();
308 
309  // Generic array for values, grads, curls, etc. that will be properly sized before each call
310  FieldContainer<double> vals;
311 
312  // Check VALUE of basis functions: resize vals to rank-2 container:
313  vals.resize(numFields, numPoints);
314  tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE);
315  for (int i = 0; i < numFields; i++) {
316  for (int j = 0; j < numPoints; j++) {
317  int l = i + j * numFields;
318  if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
319  errorFlag++;
320  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
321 
322  // Output the multi-index of the value where the error is:
323  *outStream << " At multi-index { ";
324  *outStream << i << " ";*outStream << j << " ";
325  *outStream << "} computed value: " << vals(i,j)
326  << " but reference value: " << basisValues[l] << "\n";
327  }
328  }
329  }
330 
331  // Check GRAD of basis function: resize vals to rank-3 container
332  vals.resize(numFields, numPoints, spaceDim);
333  tetBasis.getValues(vals, tetNodes, OPERATOR_GRAD);
334  for (int i = 0; i < numFields; i++) {
335  for (int j = 0; j < numPoints; j++) {
336  for (int k = 0; k < spaceDim; k++) {
337  int l = k + i * spaceDim + j * spaceDim * numFields;
338  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
339  errorFlag++;
340  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
341 
342  // Output the multi-index of the value where the error is:
343  *outStream << " At multi-index { ";
344  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
345  *outStream << "} computed grad component: " << vals(i,j,k)
346  << " but reference grad component: " << basisGrads[l] << "\n";
347  }
348  }
349  }
350  }
351 
352  // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
353  tetBasis.getValues(vals, tetNodes, OPERATOR_D1);
354  for (int i = 0; i < numFields; i++) {
355  for (int j = 0; j < numPoints; j++) {
356  for (int k = 0; k < spaceDim; k++) {
357  int l = k + i * spaceDim + j * spaceDim * numFields;
358  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
359  errorFlag++;
360  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
361 
362  // Output the multi-index of the value where the error is:
363  *outStream << " At multi-index { ";
364  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
365  *outStream << "} computed D1 component: " << vals(i,j,k)
366  << " but reference D1 component: " << basisGrads[l] << "\n";
367  }
368  }
369  }
370  }
371 
372  // Check all higher derivatives - must be zero.
373  for(EOperator op = OPERATOR_D2; op < OPERATOR_MAX; op++) {
374 
375  // The last dimension is the number of kth derivatives and needs to be resized for every Dk
376  int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
377  vals.resize(numFields, numPoints, DkCardin);
378 
379  tetBasis.getValues(vals, tetNodes, op);
380  for (int i = 0; i < vals.size(); i++) {
381  if (std::abs(vals[i]) > INTREPID_TOL) {
382  errorFlag++;
383  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
384 
385  // Get the multi-index of the value where the error is and the operator order
386  std::vector<int> myIndex;
387  vals.getMultiIndex(myIndex,i);
388  int ord = Intrepid::getOperatorOrder(op);
389  *outStream << " At multi-index { ";
390  for(int j = 0; j < vals.rank(); j++) {
391  *outStream << myIndex[j] << " ";
392  }
393  *outStream << "} computed D"<< ord <<" component: " << vals[i]
394  << " but reference D" << ord << " component: 0 \n";
395  }
396  }
397  }
398  }
399 
400  // Catch unexpected errors
401  catch (const std::logic_error & err) {
402  *outStream << err.what() << "\n\n";
403  errorFlag = -1000;
404  };
405 
406  if (errorFlag != 0)
407  std::cout << "End Result: TEST FAILED\n";
408  else
409  std::cout << "End Result: TEST PASSED\n";
410 
411  // reset format state of std::cout
412  std::cout.copyfmt(oldFormatState);
413 
414  return errorFlag;
415 }
Intrepid::FieldContainer::size
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
Definition: Intrepid_FieldContainerDef.hpp:601
Intrepid::Basis::getDofTag
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
Definition: Intrepid_BasisDef.hpp:78
Intrepid::Basis::getCardinality
virtual int getCardinality() const
Returns cardinality of the basis.
Definition: Intrepid_BasisDef.hpp:100
Intrepid::Basis::getAllDofTags
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Definition: Intrepid_BasisDef.hpp:89
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid::FieldContainer::getMultiIndex
void getMultiIndex(int &i0, const int valueEnum) const
Returns the multi-index of a value, based on its enumeration, as a list, for rank-1 containers...
Definition: Intrepid_FieldContainerDef.hpp:863
Intrepid::FieldContainer::rank
int rank() const
Return rank of the FieldContainer = number of indices used to tag the multi-indexed value...
Definition: Intrepid_FieldContainerDef.hpp:594
Intrepid::Basis_HGRAD_TET_C1_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Tetrahedron cell.
Definition: Intrepid_HGRAD_TET_C1_FEMDef.hpp:94
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::Basis::getDofOrdinal
virtual int getDofOrdinal(const int subcDim, const int subcOrd, const int subcDofOrd)
DoF tag to ordinal lookup.
Definition: Intrepid_BasisDef.hpp:51
Intrepid::Basis::getBaseCellTopology
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
Definition: Intrepid_BasisDef.hpp:112
Intrepid::Basis_HGRAD_TET_C1_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Tetrahedron cell...
Definition: Intrepid_HGRAD_TET_C1_FEM.hpp:81
Generated by   doxygen 1.8.5