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_HCURL_HEX_In_FEM.hpp"
50 #include "Intrepid_PointTools.hpp"
51 #include "Teuchos_oblackholestream.hpp"
52 #include "Teuchos_RCP.hpp"
53 #include "Teuchos_GlobalMPISession.hpp"
54 
55 using namespace std;
56 using namespace Intrepid;
57 
58 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
59 { \
60  ++nException; \
61  try { \
62  S ; \
63  } \
64  catch (const std::logic_error & err) { \
65  ++throwCounter; \
66  *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
67  *outStream << err.what() << '\n'; \
68  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
69  }; \
70 }
71 
72 int main(int argc, char *argv[]) {
73 
74  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
75 
76  // This little trick lets us print to std::cout only if
77  // a (dummy) command-line argument is provided.
78  int iprint = argc - 1;
79  Teuchos::RCP<std::ostream> outStream;
80  Teuchos::oblackholestream bhs; // outputs nothing
81  if (iprint > 0)
82  outStream = Teuchos::rcp(&std::cout, false);
83  else
84  outStream = Teuchos::rcp(&bhs, false);
85 
86  // Save the format state of the original std::cout.
87  Teuchos::oblackholestream oldFormatState;
88  oldFormatState.copyfmt(std::cout);
89 
90  *outStream \
91  << "===============================================================================\n" \
92  << "| |\n" \
93  << "| Unit Test (Basis_HCURL_HEX_In_FEM) |\n" \
94  << "| |\n" \
95  << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
96  << "| 2) Basis values for VALUE and CURL 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
110  const int deg = 1;
111  shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
112  FieldContainer<double> closedPts(PointTools::getLatticeSize(line,deg),1);
113  FieldContainer<double> openPts(PointTools::getLatticeSize(line,deg+1,1),1);
114  PointTools::getLattice<double,FieldContainer<double> >(closedPts,line,deg);
115  PointTools::getLattice<double,FieldContainer<double> >(openPts,line,deg+1,1);
116 
117  Basis_HCURL_HEX_In_FEM<double, FieldContainer<double> > hexBasis(deg,closedPts,openPts);
118  int errorFlag = 0;
119 
120  // Initialize throw counter for exception testing
121  int nException = 0;
122  int throwCounter = 0;
123 
124  // Define array containing the 8 vertices of the reference HEX, its center and 6 face centers
125  FieldContainer<double> hexNodes(8, 3);
126  hexNodes(0,0) = -1.0; hexNodes(0,1) = -1.0; hexNodes(0,2) = -1.0;
127  hexNodes(1,0) = 1.0; hexNodes(1,1) = -1.0; hexNodes(1,2) = -1.0;
128  hexNodes(2,0) = -1.0; hexNodes(2,1) = 1.0; hexNodes(2,2) = -1.0;
129  hexNodes(3,0) = 1.0; hexNodes(3,1) = 1.0; hexNodes(3,2) = -1.0;
130  hexNodes(4,0) = -1.0; hexNodes(4,1) = -1.0; hexNodes(4,2) = 1.0;
131  hexNodes(5,0) = 1.0; hexNodes(5,1) = -1.0; hexNodes(5,2) = 1.0;
132  hexNodes(6,0) = -1.0; hexNodes(6,1) = 1.0; hexNodes(6,2) = 1.0;
133  hexNodes(7,0) = 1.0; hexNodes(7,1) = 1.0; hexNodes(7,2) = 1.0;
134 
135 
136 
137  // Generic array for the output values; needs to be properly resized depending on the operator type
138  FieldContainer<double> vals;
139 
140  try{
141  // exception #1: GRAD cannot be applied to HCURL functions
142  // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
143  vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0), 4 );
144  INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD), throwCounter, nException );
145 
146  // exception #2: DIV cannot be applied to HCURL functions
147  // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
148  vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0) );
149  INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_DIV), throwCounter, nException );
150 
151  // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
152  // getDofTag() to access invalid array elements thereby causing bounds check exception
153  // exception #3
154  INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(3,0,0), throwCounter, nException );
155  // exception #4
156  INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(1,1,1), throwCounter, nException );
157  // exception #5
158  INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(0,4,1), throwCounter, nException );
159  // exception #6
160  INTREPID_TEST_COMMAND( hexBasis.getDofTag(12), throwCounter, nException );
161  // exception #7
162  INTREPID_TEST_COMMAND( hexBasis.getDofTag(-1), throwCounter, nException );
163 
164 #ifdef HAVE_INTREPID_DEBUG
165  // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays
166  // exception #8: input points array must be of rank-2
167  FieldContainer<double> badPoints1(4, 5, 3);
168  INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
169 
170  // exception #9 dimension 1 in the input point array must equal space dimension of the cell
171  FieldContainer<double> badPoints2(4, 2);
172  INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
173 
174  // exception #10 output values must be of rank-3 for OPERATOR_VALUE
175  FieldContainer<double> badVals1(4, 3);
176  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
177 
178  // exception #11 output values must be of rank-3 for OPERATOR_CURL
179  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_CURL), throwCounter, nException );
180 
181  // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
182  FieldContainer<double> badVals2(hexBasis.getCardinality() + 1, hexNodes.dimension(0), 3);
183  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
184 
185  // exception #13 incorrect 1st dimension of output array (must equal number of points)
186  FieldContainer<double> badVals3(hexBasis.getCardinality(), hexNodes.dimension(0) + 1, 3);
187  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
188 
189  // exception #14: incorrect 2nd dimension of output array (must equal the space dimension)
190  FieldContainer<double> badVals4(hexBasis.getCardinality(), hexNodes.dimension(0), 4);
191  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
192 
193  // exception #15: incorrect 2nd dimension of output array (must equal the space dimension)
194  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_CURL), throwCounter, nException ) ;
195 
196  // exception #16: D2 cannot be applied to HCURL functions
197  // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
198  vals.resize(hexBasis.getCardinality(),
199  hexNodes.dimension(0),
200  Intrepid::getDkCardinality(OPERATOR_D2, hexBasis.getBaseCellTopology().getDimension()));
201  INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_D2), throwCounter, nException );
202 
203 #endif
204 
205  }
206  catch (const std::logic_error & err) {
207  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
208  *outStream << err.what() << '\n';
209  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
210  errorFlag = -1000;
211  };
212 
213  // Check if number of thrown exceptions matches the one we expect
214  // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match.
215  if (throwCounter != nException) {
216  errorFlag++;
217  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
218  }
219 //#endif
220 
221  *outStream \
222  << "\n"
223  << "===============================================================================\n"\
224  << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
225  << "===============================================================================\n";
226 
227  try{
228  std::vector<std::vector<int> > allTags = hexBasis.getAllDofTags();
229 
230  // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
231  for (unsigned i = 0; i < allTags.size(); i++) {
232  int bfOrd = hexBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
233 
234 // for (unsigned j=0;j<4;j++) std::cout << allTags[i][j] << " "; std::cout << std::endl;
235 
236  std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
237  if( !( (myTag[0] == allTags[i][0]) &&
238  (myTag[1] == allTags[i][1]) &&
239  (myTag[2] == allTags[i][2]) &&
240  (myTag[3] == allTags[i][3]) ) ) {
241  errorFlag++;
242  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
243  *outStream << " getDofOrdinal( {"
244  << allTags[i][0] << ", "
245  << allTags[i][1] << ", "
246  << allTags[i][2] << ", "
247  << allTags[i][3] << "}) = " << bfOrd <<" but \n";
248  *outStream << " getDofTag(" << bfOrd << ") = { "
249  << myTag[0] << ", "
250  << myTag[1] << ", "
251  << myTag[2] << ", "
252  << myTag[3] << "}\n";
253  }
254  }
255 
256  // Now do the same but loop over basis functions
257  for( int bfOrd = 0; bfOrd < hexBasis.getCardinality(); bfOrd++) {
258  std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
259  int myBfOrd = hexBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
260  if( bfOrd != myBfOrd) {
261  errorFlag++;
262  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
263  *outStream << " getDofTag(" << bfOrd << ") = { "
264  << myTag[0] << ", "
265  << myTag[1] << ", "
266  << myTag[2] << ", "
267  << myTag[3] << "} but getDofOrdinal({"
268  << myTag[0] << ", "
269  << myTag[1] << ", "
270  << myTag[2] << ", "
271  << myTag[3] << "} ) = " << myBfOrd << "\n";
272  }
273  }
274  }
275  catch (const std::logic_error & err){
276  *outStream << err.what() << "\n\n";
277  errorFlag = -1000;
278  };
279 
280  *outStream \
281  << "\n"
282  << "===============================================================================\n"\
283  << "| TEST 3: correctness of basis function values |\n"\
284  << "===============================================================================\n";
285 
286  outStream -> precision(20);
287 
288  // VALUE: Each row pair gives the 12x3 correct basis set values at an evaluation point: (P,F,D) layout
289  double basisValues[] = {
290  1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
291  0,0,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
292  0,0,0, 0,0,0, 0,0,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0,
293  0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 1,0,0, 1,0,0,
294  0,1,0, 0,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
295  0,0,0, 0,1,0, 0,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
296  0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,1,0, 0,0,0, 0,1,0, 0,0,0,
297  0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,1,0, 0,0,0, 0,1,0,
298  0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0,
299  0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0,
300  0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0,
301  0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1
302  };
303 
304  // CURL
305 
306  double basisCurls[] = {
307  0,-0.5,0.5, 0,-0.5,0.5, 0,0,0.5, 0,0,0.5, 0,-0.5,0, 0,-0.5,0, 0,0,0, 0,0,0,
308  0,0,-0.5, 0,0,-0.5, 0,-0.5,-0.5, 0,-0.5,-0.5, 0,0,0, 0,0,0, 0,-0.5,0, 0,-0.5,0,
309  0,0.5,0, 0,0.5,0, 0,0,0, 0,0,0, 0,0.5,0.5, 0,0.5,0.5, 0,0,0.5, 0,0,0.5,
310  0,0,0, 0,0,0, 0,0.5,0, 0,0.5,0, 0,0,-0.5, 0,0,-0.5, 0,0.5,-0.5, 0,0.5,-0.5,
311  // y-component basis functions
312  // first y-component bf is (0,1/4(1-x)(1-z),0)
313  // curl is (1/4(1-x),0,-1/4(1-z))
314  0.5,0,-0.5, 0,0,-0.5, 0.5,0,-0.5, 0,0,-0.5, 0.5,0,0, 0,0,0, 0.5,0,0, 0,0,0,
315  // second y-component bf is (0,1/4(1+x)(1-z),0)
316  // curl is (1/4(1+x),0,1/4(1-z))
317  0,0,0.5, 0.5,0,0.5, 0,0,0.5, 0.5,0,0.5, 0,0,0, 0.5,0,0, 0,0,0, 0.5,0,0,
318  // third y-component bf is (0,1/4(1-x)(1+z),0)
319  // curl is (-1/4(1-x),0,-1/4(1+z))
320  -0.5,0,0, 0,0,0, -0.5,0,0, 0,0,0, -0.5,0,-0.5, 0,0,-0.5, -0.5,0,-0.5, 0,0,-0.5,
321  // fourth y-component bf is (0,1/4(1+x)(1+z),0)
322  // curl is (-1/4(1+x),0,1/4(1+z))
323  0.0,0,0, -0.5,0,0, 0.0,0,0, -0.5,0,0, 0.0,0,0.5, -0.5,0,0.5, 0.0,0,0.5, -0.5,0,0.5,
324  // first z-component bf is (0,0,1/4(1-x)(1-y))
325  // curl is (-1/4(1-x),1/4(1-y),0)
326  -0.5,0.5,0, 0,0.5,0, -0.5,0,0, 0,0,0, -0.5,0.5,0, 0,0.5,0, -0.5,0,0, 0,0,0,
327  // second z-component bf is (0,0,1/4(1+x)(1-y))
328  // curl is (-1/4(1+x),1/4(1-y),0)
329  0.0,-0.5,0, -0.5,-0.5,0, 0,0,0, -0.5,0,0, 0,-0.5,0, -0.5,-0.5,0, 0,0,0, -0.5,0,0,
330  // third z-component bf is (0,0,1/4(1-x)(1+y))
331  // curl is (1/4(1-x),1/4(1+y),0)
332  0.5,0,0, 0,0,0, 0.5,0.5,0, 0,0.5,0, 0.5,0,0, 0,0,0, 0.5,0.5,0, 0,0.5,0,
333  // fourth z-component bf is (0,0,1/4(1+x)(1+y))
334  // curl is (1/4(1+x),-1/4(1+y),0)
335  0,0,0, 0.5,0,0, 0,-0.5,0, 0.5,-0.5,0, 0,0,0, 0.5,0,0, 0,-0.5,0, 0.5,-0.5,0
336  };
337 
338 
339  try{
340 
341  // Dimensions for the output arrays:
342  int numFields = hexBasis.getCardinality();
343  int numPoints = hexNodes.dimension(0);
344  int spaceDim = hexBasis.getBaseCellTopology().getDimension();
345 
346  // Generic array for values and curls that will be properly sized before each call
347  FieldContainer<double> vals;
348 
349  // Check VALUE of basis functions: resize vals to rank-3 container:
350  vals.resize(numFields, numPoints, spaceDim);
351  hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE);
352  for (int i = 0; i < numFields; i++) {
353  for (int j = 0; j < numPoints; j++) {
354  for (int k = 0; k < spaceDim; k++) {
355 
356  // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k
357  int l = k + i * spaceDim * numPoints + j * spaceDim;
358  if (std::abs(vals(i,j,k) - basisValues[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 value: " << vals(i,j,k)
366  << " but reference value: " << basisValues[l] << "\n";
367  }
368  }
369  }
370  }
371 
372  // Check CURL of basis function: resize vals to rank-3 container
373  vals.resize(numFields, numPoints, spaceDim);
374  hexBasis.getValues(vals, hexNodes, OPERATOR_CURL);
375  for (int i = 0; i < numFields; i++) {
376  for (int j = 0; j < numPoints; j++) {
377  for (int k = 0; k < spaceDim; k++) {
378  // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k
379  int l = k + i * spaceDim * numPoints + j * spaceDim;
380  if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
381  errorFlag++;
382  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
383 
384  // Output the multi-index of the value where the error is:
385  *outStream << " At multi-index { ";
386  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
387  *outStream << "} computed curl component: " << vals(i,j,k)
388  << " but reference curl component: " << basisCurls[l] << "\n";
389  }
390  }
391  }
392  }
393 
394  }
395 
396  // Catch unexpected errors
397  catch (const std::logic_error & err) {
398  *outStream << err.what() << "\n\n";
399  errorFlag = -1000;
400  };
401 
402  if (errorFlag != 0)
403  std::cout << "End Result: TEST FAILED\n";
404  else
405  std::cout << "End Result: TEST PASSED\n";
406 
407  // reset format state of std::cout
408  std::cout.copyfmt(oldFormatState);
409 
410  return errorFlag;
411 }
Intrepid_PointTools.hpp
Header file for utility class to provide point tools, such as barycentric coordinates, equispaced lattices, and warp-blend point distrubtions.
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
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_HCURL_HEX_In_FEM
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedral cell...
Definition: Intrepid_HCURL_HEX_In_FEM.hpp:69
Intrepid_HCURL_HEX_In_FEM.hpp
Header file for the Intrepid::HCURL_HEX_In_FEM class.
Generated by   doxygen 1.8.5