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