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 
49 #include "Intrepid_HGRAD_PYR_C1_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_HGRAD_PYR_C1_FEM) |\n" \
93  << "| |\n" \
94  << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95  << "| 2) Basis values for VALUE, GRAD, CURL, and Dk 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
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 PYRMID and some other points.
117  FieldContainer<double> pyrNodes(10, 3);
118  pyrNodes(0,0) = -1.0; pyrNodes(0,1) = -1.0; pyrNodes(0,2) = 0;
119  pyrNodes(1,0) = 1.0; pyrNodes(1,1) = -1.0; pyrNodes(1,2) = 0;
120  pyrNodes(2,0) = 1.0; pyrNodes(2,1) = 1.0; pyrNodes(2,2) = 0;
121  pyrNodes(3,0) = -1.0; pyrNodes(3,1) = 1.0; pyrNodes(3,2) = 0;
122  pyrNodes(4,0) = 0.0; pyrNodes(4,1) = 0.0; pyrNodes(4,2) = 1.0;
123 
124  pyrNodes(5,0) = 0.25; pyrNodes(5,1) = 0.5; pyrNodes(5,2) = 0.2;
125  pyrNodes(6,0) = -0.7 ; pyrNodes(6,1) = 0.3; pyrNodes(6,2) = 0.3;
126  pyrNodes(7,0) = 0.; pyrNodes(7,1) = -0.05; pyrNodes(7,2) = 0.95;
127  pyrNodes(8,0) = -0.15; pyrNodes(8,1) = -0.2; pyrNodes(8,2) = 0.75;
128  pyrNodes(9,0) = -0.4; pyrNodes(9,1) = 0.9; pyrNodes(9,2) = 0.0;
129 
130 
131  // Generic array for the output values; needs to be properly resized depending on the operator type
133 
134  try{
135  // exception #1: CURL cannot be applied to scalar functions
136  // resize vals to rank-3 container with dimensions (num. points, num. basis functions)
137  vals.resize(pyrBasis.getCardinality(), pyrNodes.dimension(0), 3 );
138  INTREPID_TEST_COMMAND( pyrBasis.getValues(vals, pyrNodes, OPERATOR_CURL), throwCounter, nException );
139 
140  // exception #2: DIV cannot be applied to scalar functions
141  // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
142  vals.resize(pyrBasis.getCardinality(), pyrNodes.dimension(0) );
143  INTREPID_TEST_COMMAND( pyrBasis.getValues(vals, pyrNodes, OPERATOR_DIV), 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( pyrBasis.getDofOrdinal(3,0,0), throwCounter, nException );
149  // exception #4
150  INTREPID_TEST_COMMAND( pyrBasis.getDofOrdinal(1,1,1), throwCounter, nException );
151  // exception #5
152  INTREPID_TEST_COMMAND( pyrBasis.getDofOrdinal(0,6,0), throwCounter, nException );
153  // exception #6
154  INTREPID_TEST_COMMAND( pyrBasis.getDofTag(7), throwCounter, nException );
155  // exception #7
156  INTREPID_TEST_COMMAND( pyrBasis.getDofTag(-1), throwCounter, nException );
157 
158 #ifdef HAVE_INTREPID_DEBUG
159  // Exceptions 8-18 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( pyrBasis.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( pyrBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
167 
168  // exception #10 output values must be of rank-2 for OPERATOR_VALUE
169  FieldContainer<double> badVals1(4, 3, 1);
170  INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals1, pyrNodes, OPERATOR_VALUE), throwCounter, nException );
171 
172  // exception #11 output values must be of rank-3 for OPERATOR_GRAD
173  FieldContainer<double> badVals2(4, 3);
174  INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals2, pyrNodes, OPERATOR_GRAD), throwCounter, nException );
175 
176  // exception #12 output values must be of rank-3 for OPERATOR_D1
177  INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals2, pyrNodes, OPERATOR_D1), throwCounter, nException );
178 
179  // exception #13 output values must be of rank-3 for OPERATOR_D2
180  //INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals2, pyrNodes, OPERATOR_D2), throwCounter, nException );
181 
182  // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
183  FieldContainer<double> badVals3(pyrBasis.getCardinality() + 1, pyrNodes.dimension(0));
184  INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals3, pyrNodes, OPERATOR_VALUE), throwCounter, nException );
185 
186  // exception #15 incorrect 1st dimension of output array (must equal number of points)
187  FieldContainer<double> badVals4(pyrBasis.getCardinality(), pyrNodes.dimension(0) + 1);
188  INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals4, pyrNodes, OPERATOR_VALUE), throwCounter, nException );
189 
190  // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
191  FieldContainer<double> badVals5(pyrBasis.getCardinality(), pyrNodes.dimension(0), 4);
192  INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals5, pyrNodes, OPERATOR_GRAD), throwCounter, nException );
193 
194  // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D)
195  //FieldContainer<double> badVals6(pyrBasis.getCardinality(), pyrNodes.dimension(0), 40);
196  //INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals6, pyrNodes, OPERATOR_D2), throwCounter, nException );
197 
198  // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D)
199  //INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals6, pyrNodes, OPERATOR_D3), throwCounter, nException );
200 #endif
201 
202  }
203  catch (const std::logic_error & err) {
204  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
205  *outStream << err.what() << '\n';
206  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
207  errorFlag = -1000;
208  };
209 
210  // Check if number of thrown exceptions matches the one we expect - 18
211  if (throwCounter != nException) {
212  errorFlag++;
213  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
214  }
215 
216  *outStream \
217  << "\n"
218  << "===============================================================================\n"\
219  << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
220  << "===============================================================================\n";
221 
222  try{
223  std::vector<std::vector<int> > allTags = pyrBasis.getAllDofTags();
224 
225  // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
226  for (unsigned i = 0; i < allTags.size(); i++) {
227  int bfOrd = pyrBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
228 
229  std::vector<int> myTag = pyrBasis.getDofTag(bfOrd);
230  if( !( (myTag[0] == allTags[i][0]) &&
231  (myTag[1] == allTags[i][1]) &&
232  (myTag[2] == allTags[i][2]) &&
233  (myTag[3] == allTags[i][3]) ) ) {
234  errorFlag++;
235  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
236  *outStream << " getDofOrdinal( {"
237  << allTags[i][0] << ", "
238  << allTags[i][1] << ", "
239  << allTags[i][2] << ", "
240  << allTags[i][3] << "}) = " << bfOrd <<" but \n";
241  *outStream << " getDofTag(" << bfOrd << ") = { "
242  << myTag[0] << ", "
243  << myTag[1] << ", "
244  << myTag[2] << ", "
245  << myTag[3] << "}\n";
246  }
247  }
248 
249  // Now do the same but loop over basis functions
250  for( int bfOrd = 0; bfOrd < pyrBasis.getCardinality(); bfOrd++) {
251  std::vector<int> myTag = pyrBasis.getDofTag(bfOrd);
252  int myBfOrd = pyrBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
253  if( bfOrd != myBfOrd) {
254  errorFlag++;
255  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
256  *outStream << " getDofTag(" << bfOrd << ") = { "
257  << myTag[0] << ", "
258  << myTag[1] << ", "
259  << myTag[2] << ", "
260  << myTag[3] << "} but getDofOrdinal({"
261  << myTag[0] << ", "
262  << myTag[1] << ", "
263  << myTag[2] << ", "
264  << myTag[3] << "} ) = " << myBfOrd << "\n";
265  }
266  }
267  }
268  catch (const std::logic_error & err){
269  *outStream << err.what() << "\n\n";
270  errorFlag = -1000;
271  };
272 
273  *outStream \
274  << "\n"
275  << "===============================================================================\n"\
276  << "| TEST 3: correctness of basis function values |\n"\
277  << "===============================================================================\n";
278 
279  outStream -> precision(20);
280 
281  // VALUE: Each row gives the 4 correct basis set values at an evaluation point
282  double basisValues[] = {
283  1.0, 0.0, 0.0, 0.0, 0.0,
284  0.0, 1.0, 0.0, 0.0, 0.0,
285  0.0, 0.0, 1.0, 0.0, 0.0,
286  0.0, 0.0, 0.0, 1.0, 0.0,
287  0.0, 0.0, 0.0, 0.0, 1.0,
288  //
289  0.0515625, 0.0984375, 0.4265625, 0.2234375, 0.2,
290  0.2, 0, 0., 0.5, 0.3,
291  0.025, 0.025, 0., 0, 0.95,
292  0.18, 0.045, 0.005, 0.02, 0.75,
293  0.035, 0.015, 0.285, 0.665, 0.,
294  };
295 
296  // GRAD and D1: each row gives the 3 x 5 correct values of the gradients of the 5 basis functions
297  double basisGrads[] = {
298  -0.5, -0.5, 0.0, 0.5, 0.0, -0.5, 0.0, 0.0, 0.0, 0.0, 0.5, -0.5, 0.0, 0.0, 1.0, \
299  -0.5, 0.0, -0.5, 0.5, -0.5, 0.0, 0.0, 0.5, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, \
300  0.0, 0.0, 0.0, 0.0, -0.5, -0.5, 0.5, 0.5, 0.0, -0.5, 0.0, -0.5, 0.0, 0.0, 1.0, \
301  0.0, -0.5, -0.5, 0.0, 0.0, 0.0, 0.5, 0.0, -0.5, -0.5, 0.5, 0.0, 0.0, 0.0, 1.0, \
302  -0.25,-0.25,-0.25, 0.25,-0.25,-0.25, 0.25, 0.25,-0.25,-0.25, 0.25,-0.25, 0.0, 0.0, 1.0, \
303  -0.09375, -0.171875, -0.201171875, 0.09375, -0.328125, -0.298828125, 0.40625, 0.328125, -0.201171875, -0.40625, 0.171875, -0.298828125, 0.0, 0.0, 1.0, \
304  -0.1428571428571429, -0.5, -0.3571428571428571, 0.1428571428571429, 0.0, -0.1428571428571429, 0.3571428571428571, 0.0, -0.3571428571428571, -0.3571428571428571, 0.5, -0.1428571428571429, 0.0, 0.0, 1.0, \
305  -0.5, -0.25, -0.25, 0.5, -0.25, -0.25, 0.0, 0.25, -0.25, 0., 0.25, -0.25, 0.0, 0.0, 1.0, \
306  -0.45, -0.4, -0.13, 0.45, -0.1, -0.37, 0.05, 0.1, -0.13, -0.05, 0.4, -0.37, 0.0, 0.0, 1.0, \
307  -0.025, -0.35, -0.34, 0.025, -0.15, -0.16, 0.475, 0.15, -0.34, -0.475, 0.35, -0.16, 0.0, 0.0, 1.0
308  };
309 
310 
311  //D2: flat array with the values of D2 applied to basis functions. Multi-index is (P,F,K)
312  const double eps = std::numeric_limits<double>::epsilon( );
313  double basisD2[] = {
314  0, 0.25,-0.25, 0,-0.25, 0.5, 0,-0.25, 0.25, 0, 0.25,-0.5, 0, 0.25,-0.25, 0,-0.25, 0.5, 0,-0.25, 0.25, 0, 0.25,-0.5, 0, 0, 0, 0, 0, 0, \
315  0, 0.25,-0.25, 0, 0.25,-0.5, 0,-0.25, 0.25, 0,-0.25, 0.5, 0, 0.25,-0.25, 0, 0.25,-0.5, 0,-0.25, 0.25, 0,-0.25, 0.5, 0, 0, 0, 0, 0, 0, \
316  0, 0.25, 0.25, 0, 0.25, 0.5, 0,-0.25,-0.25, 0,-0.25,-0.5, 0, 0.25, 0.25, 0, 0.25, 0.5, 0,-0.25,-0.25, 0,-0.25,-0.5, 0, 0, 0, 0, 0, 0, \
317  0, 0.25, 0.25, 0,-0.25,-0.5, 0,-0.25,-0.25, 0, 0.25, 0.5, 0, 0.25, 0.25, 0,-0.25,-0.5, 0,-0.25,-0.25, 0, 0.25, 0.5, 0, 0, 0, 0, 0, 0, \
318  0, 0.25/eps, 0, 0, 0, 0, 0,-0.25/eps, 0, 0, 0, 0, 0, 0.25/eps, 0, 0, 0, 0, 0,-0.25/eps, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, \
319  0, 0.3125, 0.1953125, 0, 0.09765625, 0.1220703125, 0,-0.3125,-0.1953125, 0,-0.09765625,-0.1220703125, 0, 0.3125, 0.1953125, 0, 0.09765625, 0.1220703125, 0,-0.3125,-0.1953125, 0,-0.09765625,-0.1220703125, 0, 0, 0, 0, 0, 0, \
320  0, 0.3571428571428571, 0.1530612244897959, 0,-0.3571428571428571,-0.306122448979592, 0,-0.3571428571428572,-0.1530612244897959, 0, 0.3571428571428571, 0.306122448979592, 0, 0.3571428571428571, 0.1530612244897959, 0,-0.3571428571428571,-0.306122448979592, 0,-0.3571428571428571,-0.1530612244897959, 0, 0.3571428571428571, 0.306122448979592, 0, 0, 0, 0, 0, 0, \
321  0, 5,-5, 0, 0, 0, 0,-5, 5, 0, 0, 0, 0, 5,-5, 0, 0, 0, 0, -5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
322  0, 1,-0.8, 0,-0.6, 0.96, 0, -1, 0.8, 0, 0.6,-0.96, 0,1,-0.8, 0, -0.6, 0.96, 0, -1, 0.8, 0, 0.6, -0.96, 0, 0, 0, 0, 0, 0, \
323  0, 0.25, 0.225, 0,-0.1,-0.18, 0,-0.25,-0.225, 0, 0.1, 0.18, 0, 0.25, 0.225, 0,-0.1,-0.18, 0,-0.25,-0.225,0,0.1,0.18, 0, 0, 0, 0, 0, 0
324  };
325  try{
326 
327  // Dimensions for the output arrays:
328  int numFields = pyrBasis.getCardinality();
329  int numPoints = pyrNodes.dimension(0);
330  int spaceDim = pyrBasis.getBaseCellTopology().getDimension();
331  int D2Cardin = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
332 
333  // Generic array for values, grads, curls, etc. that will be properly sized before each call
335 
336  // Check VALUE of basis functions: resize vals to rank-2 container:
337  vals.resize(numFields, numPoints);
338  pyrBasis.getValues(vals, pyrNodes, OPERATOR_VALUE);
339  for (int i = 0; i < numFields; i++) {
340  for (int j = 0; j < numPoints; j++) {
341  int l = i + j * numFields;
342  if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
343  errorFlag++;
344  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
345 
346  // Output the multi-index of the value where the error is:
347  *outStream << " At multi-index { ";
348  *outStream << i << " ";*outStream << j << " ";
349  *outStream << "} computed value: " << vals(i,j)
350  << " but reference value: " << basisValues[l] << "\n";
351  }
352  }
353  }
354 
355  // Check GRAD of basis function: resize vals to rank-3 container
356  vals.resize(numFields, numPoints, spaceDim);
357  pyrBasis.getValues(vals, pyrNodes, OPERATOR_GRAD);
358  for (int i = 0; i < numFields; i++) {
359  for (int j = 0; j < numPoints; j++) {
360  for (int k = 0; k < spaceDim; k++) {
361  int l = k + i * spaceDim + j * spaceDim * numFields;
362  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
363  errorFlag++;
364  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
365 
366  // Output the multi-index of the value where the error is:
367  *outStream << " At multi-index { ";
368  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
369  *outStream << "} computed grad component: " << vals(i,j,k)
370  << " but reference grad component: " << basisGrads[l] << "\n";
371  }
372  }
373  }
374  }
375 
376  // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
377  pyrBasis.getValues(vals, pyrNodes, OPERATOR_D1);
378  for (int i = 0; i < numFields; i++) {
379  for (int j = 0; j < numPoints; j++) {
380  for (int k = 0; k < spaceDim; k++) {
381  int l = k + i * spaceDim + j * spaceDim * numFields;
382  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
383  errorFlag++;
384  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
385 
386  // Output the multi-index of the value where the error is:
387  *outStream << " At multi-index { ";
388  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
389  *outStream << "} computed D1 component: " << vals(i,j,k)
390  << " but reference D1 component: " << basisGrads[l] << "\n";
391  }
392  }
393  }
394  }
395 
396  // Check D2 of basis function
397  vals.resize(numFields, numPoints, D2Cardin);
398  pyrBasis.getValues(vals, pyrNodes, OPERATOR_D2);
399  for (int i = 0; i < numFields; i++) {
400  for (int j = 0; j < numPoints; j++) {
401  if(j == 4) continue; //derivatives are singular when z = 1.
402  for (int k = 0; k < D2Cardin; k++) {
403  int l = k + i * D2Cardin + j * D2Cardin * numFields;
404  if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
405  errorFlag++;
406  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
407 
408  // Output the multi-index of the value where the error is:
409  *outStream << " At multi-index { ";
410  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
411  *outStream << "} computed D2 component: " << vals(i,j,k)
412  << " but reference D2 component: " << basisD2[l] << "\n";
413  }
414  }
415  }
416  }
417  }
418 
419 
420  // Catch unexpected errors
421  catch (const std::logic_error & err) {
422  *outStream << err.what() << "\n\n";
423  errorFlag = -1000;
424  };
425 
426  if (errorFlag != 0)
427  std::cout << "End Result: TEST FAILED\n";
428  else
429  std::cout << "End Result: TEST PASSED\n";
430 
431  // reset format state of std::cout
432  std::cout.copyfmt(oldFormatState);
433 
434  return errorFlag;
435 }
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
virtual int getCardinality() const
Returns cardinality of the basis.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Header file for utility class to provide multidimensional containers.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Pyramid cell.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Pyramid cell...
virtual int getDofOrdinal(const int subcDim, const int subcOrd, const int subcDofOrd)
DoF tag to ordinal lookup.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...