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 
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 (std::logic_error err) { \
65  ++throwCounter; \
66  *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
67  *outStream << err.what() << '\n'; \
68  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
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_TRI_C2_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 (dridzal@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 
109  // Define basis and error flag
111  int errorFlag = 0;
112 
113  // Initialize throw counter for exception testing
114  int nException = 0;
115  int throwCounter = 0;
116 
117  // The unisolvent set for the default Lagrangian basis of degree 2 uses the standard equispaced
118  // lattice on Triangle consisting of the 3 vertices and the 3 edge midpoints. To check basis
119  // correctness it suffices to evaluate each basis function at the lattice points, defined below,
120  // but we throw in an additional random point.
121  FieldContainer<double> triNodes(7, 2);
122  triNodes(0,0) = 0.0; triNodes(0,1) = 0.0;
123  triNodes(1,0) = 1.0; triNodes(1,1) = 0.0;
124  triNodes(2,0) = 0.0; triNodes(2,1) = 1.0;
125  triNodes(3,0) = 0.5; triNodes(3,1) = 0.0;
126  triNodes(4,0) = 0.5; triNodes(4,1) = 0.5;
127  triNodes(5,0) = 0.0; triNodes(5,1) = 0.5;
128  triNodes(6,0) = 0.1732; triNodes(6,1) = 0.6714;
129 
130  // Generic array for the output values; needs to be properly resized depending on the operator type
132 
133  try{
134  // exception #1: DIV cannot be applied to scalar functions
135  // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
136  vals.resize(triBasis.getCardinality(), triNodes.dimension(0) );
137  INTREPID_TEST_COMMAND( triBasis.getValues(vals, triNodes, OPERATOR_DIV), throwCounter, nException );
138 
139  // Exceptions 2-6: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
140  // getDofTag() to access invalid array elements thereby causing bounds check exception
141  // exception #2
142  INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(2,0,0), throwCounter, nException );
143  // exception #3
144  INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(1,1,1), throwCounter, nException );
145  // exception #4
146  INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(0,4,0), throwCounter, nException );
147  // exception #5
148  INTREPID_TEST_COMMAND( triBasis.getDofTag(6), throwCounter, nException );
149  // exception #6
150  INTREPID_TEST_COMMAND( triBasis.getDofTag(-1), throwCounter, nException );
151 
152 #ifdef HAVE_INTREPID_DEBUG
153  // Exceptions 7-17 test exception handling with incorrectly dimensioned input/output arrays
154  // exception #7: input points array must be of rank-2
155  FieldContainer<double> badPoints1(4, 5, 3);
156  INTREPID_TEST_COMMAND( triBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
157 
158  // exception #8 dimension 1 in the input point array must equal space dimension of the cell
159  FieldContainer<double> badPoints2(4, triBasis.getBaseCellTopology().getDimension() + 1);
160  INTREPID_TEST_COMMAND( triBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
161 
162  // exception #9 output values must be of rank-2 for OPERATOR_VALUE
163  FieldContainer<double> badVals1(4, 3, 1);
164  INTREPID_TEST_COMMAND( triBasis.getValues(badVals1, triNodes, OPERATOR_VALUE), throwCounter, nException );
165 
166  // exception #10 output values must be of rank-3 for OPERATOR_GRAD
167  FieldContainer<double> badVals2(4, 3);
168  INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_GRAD), throwCounter, nException );
169 
170  // exception #11 output values must be of rank-3 for OPERATOR_CURL
171  INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_CURL), throwCounter, nException );
172 
173  // exception #12 output values must be of rank-3 for OPERATOR_D2
174  INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_D2), throwCounter, nException );
175 
176  // exception #13 incorrect 1st dimension of output array (must equal number of basis functions)
177  FieldContainer<double> badVals3(triBasis.getCardinality() + 1, triNodes.dimension(0));
178  INTREPID_TEST_COMMAND( triBasis.getValues(badVals3, triNodes, OPERATOR_VALUE), throwCounter, nException );
179 
180  // exception #14 incorrect 0th dimension of output array (must equal number of points)
181  FieldContainer<double> badVals4(triBasis.getCardinality(), triNodes.dimension(0) + 1);
182  INTREPID_TEST_COMMAND( triBasis.getValues(badVals4, triNodes, OPERATOR_VALUE), throwCounter, nException );
183 
184  // exception #15: incorrect 2nd dimension of output array (must equal the space dimension)
185  FieldContainer<double> badVals5(triBasis.getCardinality(), triNodes.dimension(0), triBasis.getBaseCellTopology().getDimension() + 1);
186  INTREPID_TEST_COMMAND( triBasis.getValues(badVals5, triNodes, OPERATOR_GRAD), throwCounter, nException );
187 
188  // exception #16: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D)
189  FieldContainer<double> badVals6(triBasis.getCardinality(), triNodes.dimension(0), 40);
190  INTREPID_TEST_COMMAND( triBasis.getValues(badVals6, triNodes, OPERATOR_D2), throwCounter, nException );
191 
192  // exception #17: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D)
193  INTREPID_TEST_COMMAND( triBasis.getValues(badVals6, triNodes, OPERATOR_D3), throwCounter, nException );
194 #endif
195 
196  }
197  catch (std::logic_error err) {
198  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
199  *outStream << err.what() << '\n';
200  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
201  errorFlag = -1000;
202  };
203 
204  // Check if number of thrown exceptions matches the one we expect
205  if (throwCounter != nException) {
206  errorFlag++;
207  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
208  }
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 = triBasis.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 = triBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
222 
223  std::vector<int> myTag = triBasis.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 < triBasis.getCardinality(); bfOrd++) {
245  std::vector<int> myTag = triBasis.getDofTag(bfOrd);
246  int myBfOrd = triBasis.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 (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: Correct basis values in (F,P) format: each row gives the 6 correct basis values at
276  // the lattice points (3 vertices and 3 edge midpoints)
277  double basisValues[] = {
278  1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.10710168,
279  0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -0.11320352,
280  0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.23015592,
281  0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.10766112,
282  0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.46514592,
283  0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.41734224
284  };
285 
286  // GRAD and D1: Correct gradients and D1 in (F,P,D) format: each row contains 6x2 values of
287  // gradients of basis functions. Same array is used to check correctness of CURL.
288  double basisGrads[] = {
289  // V0 V1 V2 M0 M1 M2 RP
290  -3.0, -3.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0,-1.0, 0.3784, 0.3784,
291  -1.0, 0.0, 3.0, 0.0, -1.0, 0.0, 1.0, 0.0, 1.0, 0.0, -1.0, 0.0, -0.3072, 0.0,
292  0.0, -1.0, 0.0,-1.0, 0.0, 3.0, 0.0, -1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.6856,
293  4.0, 0.0, -4.0,-4.0, 0.0, 0.0, 0.0, -2.0, -2.0,-2.0, 2.0, 0.0, -0.0712, -0.6928,
294  0.0, 0.0, 0.0, 4.0, 4.0, 0.0, 0.0, 2.0, 2.0, 2.0, 2.0, 0.0, 2.6856, 0.6928,
295  0.0, 4.0, 0.0, 0.0, -4.0,-4.0, 0.0, 2.0, -2.0,-2.0, -2.0, 0.0, -2.6856, -2.064
296  };
297 
298  // D2: Correct multiset of second order partials in (F,P,Dk) format. D2 cardinality = 3 for 2D
299  double basisD2[] = {
300  // V0 V1 V2 M0 M1 M2 RP
301  4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0,
302  4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0,
303  0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0,
304  -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0,
305  0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0,
306  0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0
307  };
308  try{
309 
310  // Dimensions for the output arrays:
311  int numFields = triBasis.getCardinality();
312  int numPoints = triNodes.dimension(0);
313  int spaceDim = triBasis.getBaseCellTopology().getDimension();
314 
315  // Generic array for values, grads, curls, etc. that will be properly sized before each call
317 
318  // Check VALUE of basis functions: resize vals to rank-2 container:
319  vals.resize(numFields, numPoints);
320  triBasis.getValues(vals, triNodes, OPERATOR_VALUE);
321  for (int i = 0; i < numFields; i++) {
322  for (int j = 0; j < numPoints; j++) {
323 
324  // Compute offset for (F,P) container
325  int l = j + i * numPoints;
326  if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
327  errorFlag++;
328  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
329 
330  // Output the multi-index of the value where the error is:
331  *outStream << " At multi-index { ";
332  *outStream << i << " ";*outStream << j << " ";
333  *outStream << "} computed value: " << vals(i,j)
334  << " but reference value: " << basisValues[l] << "\n";
335  }
336  }
337  }
338 
339  // Check GRAD of basis function: resize vals to rank-3 container
340  vals.resize(numFields, numPoints, spaceDim);
341  triBasis.getValues(vals, triNodes, OPERATOR_GRAD);
342  for (int i = 0; i < numFields; i++) {
343  for (int j = 0; j < numPoints; j++) {
344  for (int k = 0; k < spaceDim; k++) {
345  // basisGrads is (F,P,D), compute offset:
346  int l = k + j * spaceDim + i * spaceDim * numPoints;
347  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
348  errorFlag++;
349  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
350 
351  // Output the multi-index of the value where the error is:
352  *outStream << " At multi-index { ";
353  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
354  *outStream << "} computed grad component: " << vals(i,j,k)
355  << " but reference grad component: " << basisGrads[l] << "\n";
356  }
357  }
358  }
359  }
360 
361  // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
362  triBasis.getValues(vals, triNodes, OPERATOR_D1);
363  for (int i = 0; i < numFields; i++) {
364  for (int j = 0; j < numPoints; j++) {
365  for (int k = 0; k < spaceDim; k++) {
366  // basisGrads is (F,P,D), compute offset:
367  int l = k + j * spaceDim + i * spaceDim * numPoints;
368  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
369  errorFlag++;
370  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
371 
372  // Output the multi-index of the value where the error is:
373  *outStream << " At multi-index { ";
374  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
375  *outStream << "} computed D1 component: " << vals(i,j,k)
376  << " but reference D1 component: " << basisGrads[l] << "\n";
377  }
378  }
379  }
380  }
381 
382  // Check CURL of basis function: resize vals just for illustration!
383  vals.resize(numFields, numPoints, spaceDim);
384  triBasis.getValues(vals, triNodes, OPERATOR_CURL);
385  for (int i = 0; i < numFields; i++) {
386  for (int j = 0; j < numPoints; j++) {
387  // We will use "rotated" basisGrads to check CURL: get offsets to extract (u_y, -u_x)
388  int curl_0 = 1 + j * spaceDim + i * spaceDim * numPoints; // position of y-derivative
389  int curl_1 = 0 + j * spaceDim + i * spaceDim * numPoints; // position of x-derivative
390 
391  double curl_value_0 = basisGrads[curl_0];
392  double curl_value_1 =-basisGrads[curl_1];
393  if (std::abs(vals(i,j,0) - curl_value_0) > INTREPID_TOL) {
394  errorFlag++;
395  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
396  // Output the multi-index of the value where the error is:
397  *outStream << " At multi-index { ";
398  *outStream << i << " ";*outStream << j << " ";*outStream << 0 << " ";
399  *outStream << "} computed curl component: " << vals(i,j,0)
400  << " but reference curl component: " << curl_value_0 << "\n";
401  }
402  if (std::abs(vals(i,j,1) - curl_value_1) > INTREPID_TOL) {
403  errorFlag++;
404  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
405  // Output the multi-index of the value where the error is:
406  *outStream << " At multi-index { ";
407  *outStream << i << " ";*outStream << j << " ";*outStream << 1 << " ";
408  *outStream << "} computed curl component: " << vals(i,j,1)
409  << " but reference curl component: " << curl_value_1 << "\n";
410  }
411  }
412  }
413 
414  // Check D2 of basis function
415  int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
416  vals.resize(numFields, numPoints, D2cardinality);
417  triBasis.getValues(vals, triNodes, OPERATOR_D2);
418  for (int i = 0; i < numFields; i++) {
419  for (int j = 0; j < numPoints; j++) {
420  for (int k = 0; k < D2cardinality; k++) {
421 
422  // basisD2 is (F,P,Dk), compute offset:
423  int l = k + j * D2cardinality + i * D2cardinality * numPoints;
424  if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
425  errorFlag++;
426  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
427 
428  // Output the multi-index of the value where the error is:
429  *outStream << " At multi-index { ";
430  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
431  *outStream << "} computed D2 component: " << vals(i,j,k)
432  << " but reference D2 component: " << basisD2[l] << "\n";
433  }
434  }
435  }
436  }
437 
438  // Check all higher derivatives - must be zero.
439  for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
440 
441  // The last dimension is the number of kth derivatives and needs to be resized for every Dk
442  int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
443  vals.resize(numFields, numPoints, DkCardin);
444 
445  triBasis.getValues(vals, triNodes, op);
446  for (int i = 0; i < vals.size(); i++) {
447  if (std::abs(vals[i]) > INTREPID_TOL) {
448  errorFlag++;
449  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
450 
451  // Get the multi-index of the value where the error is and the operator order
452  std::vector<int> myIndex;
453  vals.getMultiIndex(myIndex,i);
454  int ord = Intrepid::getOperatorOrder(op);
455  *outStream << " At multi-index { ";
456  for(int j = 0; j < vals.rank(); j++) {
457  *outStream << myIndex[j] << " ";
458  }
459  *outStream << "} computed D"<< ord <<" component: " << vals[i]
460  << " but reference D" << ord << " component: 0 \n";
461  }
462  }
463  }
464  }
465 
466  // Catch unexpected errors
467  catch (std::logic_error err) {
468  *outStream << err.what() << "\n\n";
469  errorFlag = -1000;
470  };
471 
472  if (errorFlag != 0)
473  std::cout << "End Result: TEST FAILED\n";
474  else
475  std::cout << "End Result: TEST PASSED\n";
476 
477  // reset format state of std::cout
478  std::cout.copyfmt(oldFormatState);
479 
480  return errorFlag;
481 }
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
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 getMultiIndex(int &i0, const int valueEnum) const
Returns the multi-index of a value, based on its enumeration, as a list, for rank-1 containers...
int rank() const
Return rank of the FieldContainer = number of indices used to tag the multi-indexed value...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
Header file for the Intrepid::HGRAD_TRI_C2_FEM class.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
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...
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Triangle cell...