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