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 
130  try{
131  // Generic array for the output values; needs to be properly resized depending on the operator type
132  FieldContainer<double> vals;
133 
134  // exception #1: GRAD cannot be applied to HCURL functions
135  // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
136  vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0), 4 );
137  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD), throwCounter, nException );
138 
139  // exception #2: DIV cannot be applied to HCURL functions
140  // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
141  vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0) );
142  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), 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( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException );
148  // exception #4
149  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException );
150  // exception #5
151  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,1), throwCounter, nException );
152  // exception #6
153  INTREPID_TEST_COMMAND( quadBasis.getDofTag(12), throwCounter, nException );
154  // exception #7
155  INTREPID_TEST_COMMAND( quadBasis.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( quadBasis.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, quadBasis.getBaseCellTopology().getDimension() + 1);
165  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
166 
167  // exception #10 output values must be of rank-3 for OPERATOR_VALUE in 2D
168  FieldContainer<double> badVals1(4, 3);
169  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
170 
171  FieldContainer<double> badCurls1(4,3,2);
172  // exception #11 output values must be of rank-2 for OPERATOR_CURL
173  INTREPID_TEST_COMMAND( quadBasis.getValues(badCurls1, quadNodes, OPERATOR_CURL), throwCounter, nException );
174 
175  // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
176  FieldContainer<double> badVals2(quadBasis.getCardinality() + 1, quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension());
177  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
178 
179  // exception #13 incorrect 1st dimension of output array (must equal number of points)
180  FieldContainer<double> badVals3(quadBasis.getCardinality(), quadNodes.dimension(0) + 1, quadBasis.getBaseCellTopology().getDimension() );
181  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
182 
183  // exception #14: incorrect 2nd dimension of output array for VALUE (must equal the space dimension)
184  FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() - 1);
185  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
186 
187  // exception #15: D2 cannot be applied to HCURL functions
188  // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
189  vals.resize(quadBasis.getCardinality(),
190  quadNodes.dimension(0),
191  Intrepid::getDkCardinality(OPERATOR_D2, quadBasis.getBaseCellTopology().getDimension()));
192  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_D2), throwCounter, nException );
193 #endif
194 
195  }
196  catch (const std::logic_error & err) {
197  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
198  *outStream << err.what() << '\n';
199  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
200  errorFlag = -1000;
201  };
202 
203  // Check if number of thrown exceptions matches the one we expect
204  // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match.
205  if (throwCounter != nException) {
206  errorFlag++;
207  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
208  }
209 //#endif
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 = quadBasis.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 = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
223 
224  std::vector<int> myTag = quadBasis.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 < quadBasis.getCardinality(); bfOrd++) {
246  std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
247  int myBfOrd = quadBasis.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 4x2 correct basis set values at an evaluation point: (P,F,D) layout
277  double basisValues[] = {
278  0.500000, 0, 0, 0, 0, 0, 0, -0.500000, 0.500000, 0, 0, 0.500000, 0, \
279  0, 0, 0, 0, 0, 0, 0.500000, -0.500000, 0, 0, 0, 0, 0, 0, 0, \
280  -0.500000, 0, 0, -0.500000, 0.250000, 0, 0, 0.250000, -0.250000, 0, \
281  0, -0.250000, 0.375000, 0, 0, 0.250000, -0.125000, 0, 0, -0.250000, \
282  0.125000, 0, 0, 0.250000, -0.375000, 0, 0, -0.250000, 0.250000, 0, 0, \
283  0.125000, -0.250000, 0, 0, -0.375000, 0.250000, 0, 0, 0.375000, \
284  -0.250000, 0, 0, -0.125000
285  };
286 
287  // CURL: correct values in (F,P) format
288  double basisCurls[] = {
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  0.25, 0.25, 0.25, 0.25
298  };
299 
300 
301  try{
302 
303  // Dimensions for the output arrays:
304  int numFields = quadBasis.getCardinality();
305  int numPoints = quadNodes.dimension(0);
306  int spaceDim = quadBasis.getBaseCellTopology().getDimension();
307 
308  // Generic array for values and curls that will be properly sized before each call
309  FieldContainer<double> vals;
310 
311  // Check VALUE of basis functions: resize vals to rank-3 container:
312  vals.resize(numFields, numPoints, spaceDim);
313  quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE);
314  for (int i = 0; i < numFields; i++) {
315  for (int j = 0; j < numPoints; j++) {
316  for (int k = 0; k < spaceDim; k++) {
317 
318  // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k
319  int l = k + i * spaceDim + j * spaceDim * numFields;
320  if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
321  errorFlag++;
322  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
323 
324  // Output the multi-index of the value where the error is:
325  *outStream << " At multi-index { ";
326  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
327  *outStream << "} computed value: " << vals(i,j,k)
328  << " but reference value: " << basisValues[l] << "\n";
329  }
330  }
331  }
332  }
333 
334  // Check CURL of basis function: resize vals to rank-2 container
335  vals.resize(numFields, numPoints);
336  quadBasis.getValues(vals, quadNodes, OPERATOR_CURL);
337  for (int i = 0; i < numFields; i++) {
338  for (int j = 0; j < numPoints; j++) {
339  int l = i + j * numFields;
340  if (std::abs(vals(i,j) - basisCurls[l]) > INTREPID_TOL) {
341  errorFlag++;
342  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
343 
344  // Output the multi-index of the value where the error is:
345  *outStream << " At multi-index { ";
346  *outStream << i << " ";*outStream << j << " ";
347  *outStream << "} computed curl component: " << vals(i,j)
348  << " but reference curl component: " << basisCurls[l] << "\n";
349  }
350  }
351  }
352  }
353 
354  // Catch unexpected errors
355  catch (const std::logic_error & err) {
356  *outStream << err.what() << "\n\n";
357  errorFlag = -1000;
358  };
359 
360  *outStream \
361  << "\n"
362  << "===============================================================================\n"\
363  << "| TEST 4: correctness of DoF locations |\n"\
364  << "===============================================================================\n";
365 
366  try{
367  Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
368  Teuchos::rcp(new Basis_HCURL_QUAD_I1_FEM<double, FieldContainer<double> >);
369  Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
370  Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
371 
372  int spaceDim = 2;
373  FieldContainer<double> cvals;
374  FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality(),2); // last dimension is spatial dim
375 
376  // Check exceptions.
377 #ifdef HAVE_INTREPID_DEBUG
378  cvals.resize(1,2,3);
379  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
380  cvals.resize(3,2);
381  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
382  cvals.resize(4,3);
383  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
384 #endif
385  cvals.resize(4,spaceDim);
386  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
387  // Check if number of thrown exceptions matches the one we expect
388  if (throwCounter != nException) {
389  errorFlag++;
390  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
391  }
392 
393  // Check mathematical correctness
394  FieldContainer<double> tangents(basis->getCardinality(),spaceDim); // tangents at each point basis point
395  tangents(0,0) = 2.0; tangents(0,1) = 0.0;
396  tangents(1,0) = 0.0; tangents(1,1) = 2.0;
397  tangents(2,0) = -2.0; tangents(2,1) = 0.0;
398  tangents(3,0) = 0.0; tangents(3,1) = -2.0;
399 
400  basis->getValues(bvals, cvals, OPERATOR_VALUE);
401  char buffer[120];
402  for (int i=0; i<bvals.dimension(0); i++) {
403  for (int j=0; j<bvals.dimension(1); j++) {
404 
405  double tangent = 0.0;
406  for(int d=0;d<spaceDim;d++)
407  tangent += bvals(i,j,d)*tangents(j,d);
408 
409  if ((i != j) && (std::abs(tangent - 0.0) > INTREPID_TOL)) {
410  errorFlag++;
411  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);
412  *outStream << buffer;
413  }
414  else if ((i == j) && (std::abs(tangent - 1.0) > INTREPID_TOL)) {
415  errorFlag++;
416  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);
417  *outStream << buffer;
418  }
419  }
420  }
421 
422  }
423  catch (const std::logic_error & err){
424  *outStream << err.what() << "\n\n";
425  errorFlag = -1000;
426  };
427 
428 
429  if (errorFlag != 0)
430  std::cout << "End Result: TEST FAILED\n";
431  else
432  std::cout << "End Result: TEST PASSED\n";
433 
434  // reset format state of std::cout
435  std::cout.copyfmt(oldFormatState);
436 
437  return errorFlag;
438 }
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