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_HDIV_WEDGE_I1_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_HDIV_WEDGE_I1_FEM) |\n" \
93  << "| |\n" \
94  << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95  << "| 2) Basis values for VALUE and DIV 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_HDIV_WEDGE_I1_FEM<double, FieldContainer<double> > wedgeBasis;
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 6 vertices of the reference WEDGE and 6 other points.
117  FieldContainer<double> wedgeNodes(12, 3);
118  wedgeNodes(0,0) = 0.0; wedgeNodes(0,1) = 0.0; wedgeNodes(0,2) = -1.0;
119  wedgeNodes(1,0) = 1.0; wedgeNodes(1,1) = 0.0; wedgeNodes(1,2) = -1.0;
120  wedgeNodes(2,0) = 0.0; wedgeNodes(2,1) = 1.0; wedgeNodes(2,2) = -1.0;
121  wedgeNodes(3,0) = 0.0; wedgeNodes(3,1) = 0.0; wedgeNodes(3,2) = 1.0;
122  wedgeNodes(4,0) = 1.0; wedgeNodes(4,1) = 0.0; wedgeNodes(4,2) = 1.0;
123  wedgeNodes(5,0) = 0.0; wedgeNodes(5,1) = 1.0; wedgeNodes(5,2) = 1.0;
124 
125  wedgeNodes(6,0) = 0.25; wedgeNodes(6,1) = 0.5; wedgeNodes(6,2) = -1.0;
126  wedgeNodes(7,0) = 0.5; wedgeNodes(7,1) = 0.25; wedgeNodes(7,2) = 0.0;
127  wedgeNodes(8,0) = 0.25; wedgeNodes(8,1) = 0.25; wedgeNodes(8,2) = 1.0;
128  wedgeNodes(9,0) = 0.25; wedgeNodes(9,1) = 0.0; wedgeNodes(9,2) = 0.75;
129  wedgeNodes(10,0)= 0.0; wedgeNodes(10,1)= 0.5; wedgeNodes(10,2)= -0.25;
130  wedgeNodes(11,0)= 0.5; wedgeNodes(11,1)= 0.5; wedgeNodes(11,2)= 0.0;
131 
132  // Generic array for the output values; needs to be properly resized depending on the operator type
133  FieldContainer<double> vals;
134 
135  try{
136  // exception #1: GRAD cannot be applied to HDIV functions
137  // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
138  vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 );
139  INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
140 
141  // exception #2: CURL cannot be applied to HDIV functions
142  INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_CURL), throwCounter, nException );
143 
144  // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
145  // getDofTag() to access invalid array elements thereby causing bounds check exception
146  // exception #3
147  INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException );
148  // exception #4
149  INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException );
150  // exception #5
151  INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,4,1), throwCounter, nException );
152  // exception #6
153  INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(11), throwCounter, nException );
154  // exception #7
155  INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException );
156 
157 #ifdef HAVE_INTREPID_DEBUG
158  // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays
159  // exception #8: input points array must be of rank-2
160  FieldContainer<double> badPoints1(4, 5, 3);
161  INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
162 
163  // exception #9 dimension 1 in the input point array must equal space dimension of the cell
164  FieldContainer<double> badPoints2(4, 2);
165  INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
166 
167  // exception #10 output values must be of rank-3 for OPERATOR_VALUE
168  FieldContainer<double> badVals1(4, 3);
169  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
170 
171  // exception #11 output values must be of rank-2 for OPERATOR_DIV
172  FieldContainer<double> badVals2(4, 3, 1);
173  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
174 
175  // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
176  FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0), 3);
177  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
178 
179  // exception #13 incorrect 0th dimension of output array (must equal number of basis functions)
180  FieldContainer<double> badVals4(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0));
181  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
182 
183  // exception #14 incorrect 1st dimension of output array (must equal number of points)
184  FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1, 3);
185  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
186 
187  // exception #15 incorrect 1st dimension of output array (must equal number of points)
188  FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1);
189  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
190 
191  // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
192  FieldContainer<double> badVals7(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 4);
193  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals7, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
194 #endif
195 
196  }
197  catch (const std::logic_error & err) {
198  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
199  *outStream << err.what() << '\n';
200  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
201  errorFlag = -1000;
202  };
203 
204  // Check if number of thrown exceptions matches the one we expect
205  if (throwCounter != nException) {
206  errorFlag++;
207  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
208  }
209 
210  *outStream \
211  << "\n"
212  << "===============================================================================\n"\
213  << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
214  << "===============================================================================\n";
215 
216  try{
217  std::vector<std::vector<int> > allTags = wedgeBasis.getAllDofTags();
218 
219  // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
220  for (unsigned i = 0; i < allTags.size(); i++) {
221  int bfOrd = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
222 
223  std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
224  if( !( (myTag[0] == allTags[i][0]) &&
225  (myTag[1] == allTags[i][1]) &&
226  (myTag[2] == allTags[i][2]) &&
227  (myTag[3] == allTags[i][3]) ) ) {
228  errorFlag++;
229  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
230  *outStream << " getDofOrdinal( {"
231  << allTags[i][0] << ", "
232  << allTags[i][1] << ", "
233  << allTags[i][2] << ", "
234  << allTags[i][3] << "}) = " << bfOrd <<" but \n";
235  *outStream << " getDofTag(" << bfOrd << ") = { "
236  << myTag[0] << ", "
237  << myTag[1] << ", "
238  << myTag[2] << ", "
239  << myTag[3] << "}\n";
240  }
241  }
242 
243  // Now do the same but loop over basis functions
244  for( int bfOrd = 0; bfOrd < wedgeBasis.getCardinality(); bfOrd++) {
245  std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
246  int myBfOrd = wedgeBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
247  if( bfOrd != myBfOrd) {
248  errorFlag++;
249  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
250  *outStream << " getDofTag(" << bfOrd << ") = { "
251  << myTag[0] << ", "
252  << myTag[1] << ", "
253  << myTag[2] << ", "
254  << myTag[3] << "} but getDofOrdinal({"
255  << myTag[0] << ", "
256  << myTag[1] << ", "
257  << myTag[2] << ", "
258  << myTag[3] << "} ) = " << myBfOrd << "\n";
259  }
260  }
261  }
262  catch (const std::logic_error & err){
263  *outStream << err.what() << "\n\n";
264  errorFlag = -1000;
265  };
266 
267  *outStream \
268  << "\n"
269  << "===============================================================================\n"\
270  << "| TEST 3: correctness of basis function values |\n"\
271  << "===============================================================================\n";
272 
273  outStream -> precision(20);
274 
275  // VALUE: Each row pair gives the 5x3 correct basis set values at an evaluation point
276  double basisValues[] = {
277  0, -0.500000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, -2.00000, 0, 0, 0, \
278  0.500000, -0.500000, 0, 0.500000, 0, 0, 0, 0, 0, 0, 0, -2.00000, 0, \
279  0, 0, 0, 0, 0, 0, 0.500000, 0, -0.500000, 0.500000, 0, 0, 0, \
280  -2.00000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, \
281  0, 0, 0, 2.00000, 0.500000, -0.500000, 0, 0.500000, 0, 0, 0, 0, 0, 0, \
282  0, 0, 0, 0, 2.00000, 0, 0, 0, 0, 0.500000, 0, -0.500000, 0.500000, 0, \
283  0, 0, 0, 0, 0, 2.00000, 0.125000, -0.250000, 0, 0.125000, 0.250000, \
284  0, -0.375000, 0.250000, 0, 0, 0, -2.00000, 0, 0, 0, 0.250000, \
285  -0.375000, 0, 0.250000, 0.125000, 0, -0.250000, 0.125000, 0, 0, 0, \
286  -1.00000, 0, 0, 1.00000, 0.125000, -0.375000, 0, 0.125000, 0.125000, \
287  0, -0.375000, 0.125000, 0, 0, 0, 0, 0, 0, 2.00000, 0.125000, \
288  -0.500000, 0, 0.125000, 0, 0, -0.375000, 0, 0, 0, 0, -0.250000, 0, 0, \
289  1.75000, 0, -0.250000, 0, 0, 0.250000, 0, -0.500000, 0.250000, 0, 0, \
290  0, -1.25000, 0, 0, 0.750000, 0.250000, -0.250000, 0, 0.250000, \
291  0.250000, 0, -0.250000, 0.250000, 0, 0, 0, -1.00000, 0, 0, 1.00000};
292 
293  // DIV: each row pair gives the 5 correct values of the divergence of the 5 basis functions
294  double basisDivs[] = {
295  // 6 vertices
296  1.0, 1.0, 1.0, 1.0, 1.0,
297  1.0, 1.0, 1.0, 1.0, 1.0,
298  1.0, 1.0, 1.0, 1.0, 1.0,
299  1.0, 1.0, 1.0, 1.0, 1.0,
300  1.0, 1.0, 1.0, 1.0, 1.0,
301  1.0, 1.0, 1.0, 1.0, 1.0,
302  // 6 other points
303  1.0, 1.0, 1.0, 1.0, 1.0,
304  1.0, 1.0, 1.0, 1.0, 1.0,
305  1.0, 1.0, 1.0, 1.0, 1.0,
306  1.0, 1.0, 1.0, 1.0, 1.0,
307  1.0, 1.0, 1.0, 1.0, 1.0,
308  1.0, 1.0, 1.0, 1.0, 1.0
309  };
310 
311  try{
312 
313  // Dimensions for the output arrays:
314  int numFields = wedgeBasis.getCardinality();
315  int numPoints = wedgeNodes.dimension(0);
316  int spaceDim = wedgeBasis.getBaseCellTopology().getDimension();
317 
318  // Generic array for values and curls that will be properly sized before each call
319  FieldContainer<double> vals;
320 
321  // Check VALUE of basis functions: resize vals to rank-3 container:
322  vals.resize(numFields, numPoints, spaceDim);
323  wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE);
324  for (int i = 0; i < numFields; i++) {
325  for (int j = 0; j < numPoints; j++) {
326  for (int k = 0; k < spaceDim; k++) {
327  int l = k + i * spaceDim + j * spaceDim * numFields;
328  if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
329  errorFlag++;
330  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
331 
332  // Output the multi-index of the value where the error is:
333  *outStream << " At multi-index { ";
334  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
335  *outStream << "} computed value: " << vals(i,j,k)
336  << " but reference value: " << basisValues[l] << "\n";
337  }
338  }
339  }
340  }
341 
342 
343  // Check DIV of basis function: resize vals to rank-2 container
344  vals.resize(numFields, numPoints);
345  wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV);
346  for (int i = 0; i < numFields; i++) {
347  for (int j = 0; j < numPoints; j++) {
348  int l = i + j * numFields;
349  if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) {
350  errorFlag++;
351  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
352 
353  // Output the multi-index of the value where the error is:
354  *outStream << " At multi-index { ";
355  *outStream << i << " ";*outStream << j << " ";
356  *outStream << "} computed divergence component: " << vals(i,j)
357  << " but reference divergence component: " << basisDivs[l] << "\n";
358  }
359  }
360  }
361 
362  }
363 
364  // Catch unexpected errors
365  catch (const std::logic_error & err) {
366  *outStream << err.what() << "\n\n";
367  errorFlag = -1000;
368  };
369 
370  if (errorFlag != 0)
371  std::cout << "End Result: TEST FAILED\n";
372  else
373  std::cout << "End Result: TEST PASSED\n";
374 
375  // reset format state of std::cout
376  std::cout.copyfmt(oldFormatState);
377 
378  return errorFlag;
379 }
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::Basis_HDIV_WEDGE_I1_FEM
Implementation of the default H(div)-compatible FEM basis of degree 1 on Wedge cell.
Definition: Intrepid_HDIV_WEDGE_I1_FEM.hpp:103
Intrepid::Basis_HDIV_WEDGE_I1_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Wedge cell.
Definition: Intrepid_HDIV_WEDGE_I1_FEMDef.hpp:94
Intrepid_HDIV_WEDGE_I1_FEM.hpp
Header file for the Intrepid::HDIV_WEDGE_I1_FEM class.
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
Generated by   doxygen 1.8.5