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_QUAD_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_HCURL_QUAD_I1_FEM) |\n" \
93  << "| |\n" \
94  << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95  << "| 2) Basis values for VALUE and CURL 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_HCURL_QUAD_I1_FEM<double, FieldContainer<double> > quadBasis;
110  int errorFlag = 0;
111 
112  // Initialize throw counter for exception testing
113  int nException = 0;
114  int throwCounter = 0;
115 
116  // Array with the 4 vertices of the reference Quadrilateral, its center and 4 more points
117  FieldContainer<double> quadNodes(9, 2);
118  quadNodes(0,0) = -1.0; quadNodes(0,1) = -1.0;
119  quadNodes(1,0) = 1.0; quadNodes(1,1) = -1.0;
120  quadNodes(2,0) = 1.0; quadNodes(2,1) = 1.0;
121  quadNodes(3,0) = -1.0; quadNodes(3,1) = 1.0;
122 
123  quadNodes(4,0) = 0.0; quadNodes(4,1) = 0.0;
124  quadNodes(5,0) = 0.0; quadNodes(5,1) = -0.5;
125  quadNodes(6,0) = 0.0; quadNodes(6,1) = 0.5;
126  quadNodes(7,0) = -0.5; quadNodes(7,1) = 0.0;
127  quadNodes(8,0) = 0.5; quadNodes(8,1) = 0.0;
128 
129  // Generic array for the output values; needs to be properly resized depending on the operator type
130  FieldContainer<double> vals;
131 
132  try{
133  // exception #1: GRAD cannot be applied to HCURL functions
134  // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
135  vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0), 4 );
136  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD), throwCounter, nException );
137 
138  // exception #2: DIV cannot be applied to HCURL functions
139  // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
140  vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0) );
141  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException );
142 
143  // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
144  // getDofTag() to access invalid array elements thereby causing bounds check exception
145  // exception #3
146  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException );
147  // exception #4
148  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException );
149  // exception #5
150  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,1), throwCounter, nException );
151  // exception #6
152  INTREPID_TEST_COMMAND( quadBasis.getDofTag(12), throwCounter, nException );
153  // exception #7
154  INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException );
155 
156 #ifdef HAVE_INTREPID_DEBUG
157  // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays
158  // exception #8: input points array must be of rank-2
159  FieldContainer<double> badPoints1(4, 5, 3);
160  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
161 
162  // exception #9 dimension 1 in the input point array must equal space dimension of the cell
163  FieldContainer<double> badPoints2(4, quadBasis.getBaseCellTopology().getDimension() + 1);
164  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
165 
166  // exception #10 output values must be of rank-3 for OPERATOR_VALUE in 2D
167  FieldContainer<double> badVals1(4, 3);
168  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
169 
170  FieldContainer<double> badCurls1(4,3,2);
171  // exception #11 output values must be of rank-2 for OPERATOR_CURL
172  INTREPID_TEST_COMMAND( quadBasis.getValues(badCurls1, quadNodes, OPERATOR_CURL), throwCounter, nException );
173 
174  // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
175  FieldContainer<double> badVals2(quadBasis.getCardinality() + 1, quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension());
176  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
177 
178  // exception #13 incorrect 1st dimension of output array (must equal number of points)
179  FieldContainer<double> badVals3(quadBasis.getCardinality(), quadNodes.dimension(0) + 1, quadBasis.getBaseCellTopology().getDimension() );
180  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
181 
182  // exception #14: incorrect 2nd dimension of output array for VALUE (must equal the space dimension)
183  FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() - 1);
184  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
185 
186  // exception #15: D2 cannot be applied to HCURL functions
187  // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
188  vals.resize(quadBasis.getCardinality(),
189  quadNodes.dimension(0),
190  Intrepid::getDkCardinality(OPERATOR_D2, quadBasis.getBaseCellTopology().getDimension()));
191  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_D2), throwCounter, nException );
192 #endif
193 
194  }
195  catch (const std::logic_error & err) {
196  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
197  *outStream << err.what() << '\n';
198  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
199  errorFlag = -1000;
200  };
201 
202  // Check if number of thrown exceptions matches the one we expect
203  // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match.
204  if (throwCounter != nException) {
205  errorFlag++;
206  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
207  }
208 //#endif
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 = quadBasis.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 = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
222 
223  std::vector<int> myTag = quadBasis.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 < quadBasis.getCardinality(); bfOrd++) {
245  std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
246  int myBfOrd = quadBasis.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 4x2 correct basis set values at an evaluation point: (P,F,D) layout
276  double basisValues[] = {
277  0.500000, 0, 0, 0, 0, 0, 0, -0.500000, 0.500000, 0, 0, 0.500000, 0, \
278  0, 0, 0, 0, 0, 0, 0.500000, -0.500000, 0, 0, 0, 0, 0, 0, 0, \
279  -0.500000, 0, 0, -0.500000, 0.250000, 0, 0, 0.250000, -0.250000, 0, \
280  0, -0.250000, 0.375000, 0, 0, 0.250000, -0.125000, 0, 0, -0.250000, \
281  0.125000, 0, 0, 0.250000, -0.375000, 0, 0, -0.250000, 0.250000, 0, 0, \
282  0.125000, -0.250000, 0, 0, -0.375000, 0.250000, 0, 0, 0.375000, \
283  -0.250000, 0, 0, -0.125000
284  };
285 
286  // CURL: correct values in (F,P) format
287  double basisCurls[] = {
288  0.25, 0.25, 0.25, 0.25,
289  0.25, 0.25, 0.25, 0.25,
290  0.25, 0.25, 0.25, 0.25,
291  0.25, 0.25, 0.25, 0.25,
292  0.25, 0.25, 0.25, 0.25,
293  0.25, 0.25, 0.25, 0.25,
294  0.25, 0.25, 0.25, 0.25,
295  0.25, 0.25, 0.25, 0.25,
296  0.25, 0.25, 0.25, 0.25
297  };
298 
299 
300  try{
301 
302  // Dimensions for the output arrays:
303  int numFields = quadBasis.getCardinality();
304  int numPoints = quadNodes.dimension(0);
305  int spaceDim = quadBasis.getBaseCellTopology().getDimension();
306 
307  // Generic array for values and curls that will be properly sized before each call
308  FieldContainer<double> vals;
309 
310  // Check VALUE of basis functions: resize vals to rank-3 container:
311  vals.resize(numFields, numPoints, spaceDim);
312  quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE);
313  for (int i = 0; i < numFields; i++) {
314  for (int j = 0; j < numPoints; j++) {
315  for (int k = 0; k < spaceDim; k++) {
316 
317  // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k
318  int l = k + i * spaceDim + j * spaceDim * numFields;
319  if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
320  errorFlag++;
321  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
322 
323  // Output the multi-index of the value where the error is:
324  *outStream << " At multi-index { ";
325  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
326  *outStream << "} computed value: " << vals(i,j,k)
327  << " but reference value: " << basisValues[l] << "\n";
328  }
329  }
330  }
331  }
332 
333  // Check CURL of basis function: resize vals to rank-2 container
334  vals.resize(numFields, numPoints);
335  quadBasis.getValues(vals, quadNodes, OPERATOR_CURL);
336  for (int i = 0; i < numFields; i++) {
337  for (int j = 0; j < numPoints; j++) {
338  int l = i + j * numFields;
339  if (std::abs(vals(i,j) - basisCurls[l]) > INTREPID_TOL) {
340  errorFlag++;
341  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
342 
343  // Output the multi-index of the value where the error is:
344  *outStream << " At multi-index { ";
345  *outStream << i << " ";*outStream << j << " ";
346  *outStream << "} computed curl component: " << vals(i,j)
347  << " but reference curl component: " << basisCurls[l] << "\n";
348  }
349  }
350  }
351  }
352 
353  // Catch unexpected errors
354  catch (const std::logic_error & err) {
355  *outStream << err.what() << "\n\n";
356  errorFlag = -1000;
357  };
358 
359  *outStream \
360  << "\n"
361  << "===============================================================================\n"\
362  << "| TEST 4: correctness of DoF locations |\n"\
363  << "===============================================================================\n";
364 
365  try{
366  Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
367  Teuchos::rcp(new Basis_HCURL_QUAD_I1_FEM<double, FieldContainer<double> >);
368  Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
369  Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
370 
371  int spaceDim = 2;
372  FieldContainer<double> cvals;
373  FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality(),2); // last dimension is spatial dim
374 
375  // Check exceptions.
376 #ifdef HAVE_INTREPID_DEBUG
377  cvals.resize(1,2,3);
378  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
379  cvals.resize(3,2);
380  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
381  cvals.resize(4,3);
382  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
383 #endif
384  cvals.resize(4,spaceDim);
385  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
386  // Check if number of thrown exceptions matches the one we expect
387  if (throwCounter != nException) {
388  errorFlag++;
389  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
390  }
391 
392  // Check mathematical correctness
393  FieldContainer<double> tangents(basis->getCardinality(),spaceDim); // tangents at each point basis point
394  tangents(0,0) = 2.0; tangents(0,1) = 0.0;
395  tangents(1,0) = 0.0; tangents(1,1) = 2.0;
396  tangents(2,0) = -2.0; tangents(2,1) = 0.0;
397  tangents(3,0) = 0.0; tangents(3,1) = -2.0;
398 
399  basis->getValues(bvals, cvals, OPERATOR_VALUE);
400  char buffer[120];
401  for (int i=0; i<bvals.dimension(0); i++) {
402  for (int j=0; j<bvals.dimension(1); j++) {
403 
404  double tangent = 0.0;
405  for(int d=0;d<spaceDim;d++)
406  tangent += bvals(i,j,d)*tangents(j,d);
407 
408  if ((i != j) && (std::abs(tangent - 0.0) > INTREPID_TOL)) {
409  errorFlag++;
410  sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), tangent, 0.0);
411  *outStream << buffer;
412  }
413  else if ((i == j) && (std::abs(tangent - 1.0) > INTREPID_TOL)) {
414  errorFlag++;
415  sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), tangent, 1.0);
416  *outStream << buffer;
417  }
418  }
419  }
420 
421  }
422  catch (const std::logic_error & err){
423  *outStream << err.what() << "\n\n";
424  errorFlag = -1000;
425  };
426 
427 
428  if (errorFlag != 0)
429  std::cout << "End Result: TEST FAILED\n";
430  else
431  std::cout << "End Result: TEST PASSED\n";
432 
433  // reset format state of std::cout
434  std::cout.copyfmt(oldFormatState);
435 
436  return errorFlag;
437 }
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::DofCoordsInterface
This is an interface class for bases whose degrees of freedom can be associated with spatial location...
Definition: Intrepid_Basis.hpp:353
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_HCURL_QUAD_I1_FEM.hpp
Header file for the Intrepid::HCURL_QUAD_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_HCURL_QUAD_I1_FEM
Implementation of the default H(curl)-compatible FEM basis of degree 1 on Quadrilateral cell...
Definition: Intrepid_HCURL_QUAD_I1_FEM.hpp:99
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_HCURL_QUAD_I1_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Quadrilateral cell.
Definition: Intrepid_HCURL_QUAD_I1_FEMDef.hpp:92
Generated by   doxygen 1.8.5