Intrepid
test_01.cpp
Go to the documentation of this file.
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_HGRAD_WEDGE_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 (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_WEDGE_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
109  Basis_HGRAD_WEDGE_C1_FEM<double, FieldContainer<double> > wedgeBasis;
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 TRIPRISM and some other points.
117  FieldContainer<double> wedgeNodes(12, 3);
118  wedgeNodes(0,0) = 0.0; wedgeNodes(0,1) = 0.0; wedgeNodes(0,2) = -1.0;
119  wedgeNodes(1,0) = 1.0; wedgeNodes(1,1) = 0.0; wedgeNodes(1,2) = -1.0;
120  wedgeNodes(2,0) = 0.0; wedgeNodes(2,1) = 1.0; wedgeNodes(2,2) = -1.0;
121  wedgeNodes(3,0) = 0.0; wedgeNodes(3,1) = 0.0; wedgeNodes(3,2) = 1.0;
122  wedgeNodes(4,0) = 1.0; wedgeNodes(4,1) = 0.0; wedgeNodes(4,2) = 1.0;
123  wedgeNodes(5,0) = 0.0; wedgeNodes(5,1) = 1.0; wedgeNodes(5,2) = 1.0;
124 
125  wedgeNodes(6,0) = 0.25; wedgeNodes(6,1) = 0.5; wedgeNodes(6,2) = -1.0;
126  wedgeNodes(7,0) = 0.5; wedgeNodes(7,1) = 0.25; wedgeNodes(7,2) = 0.0;
127  wedgeNodes(8,0) = 0.25; wedgeNodes(8,1) = 0.30; wedgeNodes(8,2) = 1.0;
128  wedgeNodes(9,0) = 0.3; wedgeNodes(9,1) = 0.0; wedgeNodes(9,2) = 0.75;
129  wedgeNodes(10,0)= 0.0; wedgeNodes(10,1)= 0.44; wedgeNodes(10,2)= -0.23;
130  wedgeNodes(11,0)= 0.4; wedgeNodes(11,1)= 0.6; wedgeNodes(11,2)= 0.0;
131 
132 
133  // Generic array for the output values; needs to be properly resized depending on the operator type
134  FieldContainer<double> vals;
135 
136  try{
137  // exception #1: CURL cannot be applied to scalar functions
138  // resize vals to rank-3 container with dimensions (num. points, num. basis functions)
139  vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 );
140  INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
141 
142  // exception #2: DIV cannot be applied to scalar functions
143  // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
144  vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) );
145  INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
146 
147  // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
148  // getDofTag() to access invalid array elements thereby causing bounds check exception
149  // exception #3
150  INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException );
151  // exception #4
152  INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException );
153  // exception #5
154  INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,6,0), throwCounter, nException );
155  // exception #6
156  INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(7), throwCounter, nException );
157  // exception #7
158  INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException );
159 
160 #ifdef HAVE_INTREPID_DEBUG
161  // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays
162  // exception #8: input points array must be of rank-2
163  FieldContainer<double> badPoints1(4, 5, 3);
164  INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
165 
166  // exception #9 dimension 1 in the input point array must equal space dimension of the cell
167  FieldContainer<double> badPoints2(4, 2);
168  INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
169 
170  // exception #10 output values must be of rank-2 for OPERATOR_VALUE
171  FieldContainer<double> badVals1(4, 3, 1);
172  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
173 
174  // exception #11 output values must be of rank-3 for OPERATOR_GRAD
175  FieldContainer<double> badVals2(4, 3);
176  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
177 
178  // exception #12 output values must be of rank-3 for OPERATOR_D1
179  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D1), throwCounter, nException );
180 
181  // exception #13 output values must be of rank-3 for OPERATOR_D2
182  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D2), throwCounter, nException );
183 
184  // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
185  FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0));
186  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
187 
188  // exception #15 incorrect 1st dimension of output array (must equal number of points)
189  FieldContainer<double> badVals4(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1);
190  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
191 
192  // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
193  FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 4);
194  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
195 
196  // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D)
197  FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 40);
198  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D2), throwCounter, nException );
199 
200  // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D)
201  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D3), throwCounter, nException );
202 #endif
203 
204  }
205  catch (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 - 18
213  if (throwCounter != nException) {
214  errorFlag++;
215  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
216  }
217 
218  *outStream \
219  << "\n"
220  << "===============================================================================\n"\
221  << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
222  << "===============================================================================\n";
223 
224  try{
225  std::vector<std::vector<int> > allTags = wedgeBasis.getAllDofTags();
226 
227  // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
228  for (unsigned i = 0; i < allTags.size(); i++) {
229  int bfOrd = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
230 
231  std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
232  if( !( (myTag[0] == allTags[i][0]) &&
233  (myTag[1] == allTags[i][1]) &&
234  (myTag[2] == allTags[i][2]) &&
235  (myTag[3] == allTags[i][3]) ) ) {
236  errorFlag++;
237  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
238  *outStream << " getDofOrdinal( {"
239  << allTags[i][0] << ", "
240  << allTags[i][1] << ", "
241  << allTags[i][2] << ", "
242  << allTags[i][3] << "}) = " << bfOrd <<" but \n";
243  *outStream << " getDofTag(" << bfOrd << ") = { "
244  << myTag[0] << ", "
245  << myTag[1] << ", "
246  << myTag[2] << ", "
247  << myTag[3] << "}\n";
248  }
249  }
250 
251  // Now do the same but loop over basis functions
252  for( int bfOrd = 0; bfOrd < wedgeBasis.getCardinality(); bfOrd++) {
253  std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
254  int myBfOrd = wedgeBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
255  if( bfOrd != myBfOrd) {
256  errorFlag++;
257  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
258  *outStream << " getDofTag(" << bfOrd << ") = { "
259  << myTag[0] << ", "
260  << myTag[1] << ", "
261  << myTag[2] << ", "
262  << myTag[3] << "} but getDofOrdinal({"
263  << myTag[0] << ", "
264  << myTag[1] << ", "
265  << myTag[2] << ", "
266  << myTag[3] << "} ) = " << myBfOrd << "\n";
267  }
268  }
269  }
270  catch (std::logic_error err){
271  *outStream << err.what() << "\n\n";
272  errorFlag = -1000;
273  };
274 
275  *outStream \
276  << "\n"
277  << "===============================================================================\n"\
278  << "| TEST 3: correctness of basis function values |\n"\
279  << "===============================================================================\n";
280 
281  outStream -> precision(20);
282 
283  // VALUE: Each row gives the 4 correct basis set values at an evaluation point
284  double basisValues[] = {
285  1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
286  0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
287  0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
288  0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
289  0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
290  0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
291  //
292  0.25, 0.25, 0.5, 0., 0., 0.,
293  0.125, 0.25, 0.125, 0.125, 0.25, 0.125,
294  0., 0., 0., 0.45, 0.25, 0.3,
295  0.0875, 0.0375, 0, 0.6125, 0.2625, 0,
296  0.3444, 0, 0.2706, 0.2156, 0, 0.1694,
297  0., 0.2, 0.3, 0., 0.2, 0.3
298  };
299 
300  // GRAD and D1: each row gives the 3 x 4 correct values of the gradients of the 4 basis functions
301  double basisGrads[] = {
302  -1., -1., -0.5, 1., 0, 0, 0, 1., 0, 0., 0., 0.5, 0., 0, 0, 0, 0., 0, \
303  -1., -1., 0., 1., 0, -0.5, 0, 1., 0, 0., 0., 0., 0., 0, 0.5, 0, 0., \
304  0, -1., -1., 0., 1., 0, 0, 0, 1., -0.5, 0., 0., 0., 0., 0, 0, 0, 0., \
305  0.5, 0., 0., -0.5, 0., 0, 0, 0, 0., 0, -1., -1., 0.5, 1., 0, 0, 0, \
306  1., 0, 0., 0., 0., 0., 0, -0.5, 0, 0., 0, -1., -1., 0., 1., 0, 0.5, \
307  0, 1., 0, 0., 0., 0., 0., 0, 0, 0, 0., -0.5, -1., -1., 0., 1., 0, 0, \
308  0, 1., 0.5, -1., -1., -0.125, 1., 0, -0.125, 0, 1., -0.25, 0., 0., \
309  0.125, 0., 0, 0.125, 0, 0., 0.25, -0.5, -0.5, -0.125, 0.5, 0, -0.25, \
310  0, 0.5, -0.125, -0.5, -0.5, 0.125, 0.5, 0, 0.25, 0, 0.5, 0.125, 0., \
311  0., -0.225, 0., 0, -0.125, 0, 0., -0.15, -1., -1., 0.225, 1., 0, \
312  0.125, 0, 1., 0.15, -0.125, -0.125, -0.35, 0.125, 0, -0.15, 0, 0.125, \
313  0, -0.875, -0.875, 0.35, 0.875, 0, 0.15, 0, 0.875, 0, -0.615, -0.615, \
314  -0.28, 0.615, 0, 0, 0, 0.615, -0.22, -0.385, -0.385, 0.28, 0.385, 0, \
315  0, 0, 0.385, 0.22, -0.5, -0.5, 0., 0.5, 0, -0.2, 0, 0.5, -0.3, -0.5, \
316  -0.5, 0., 0.5, 0, 0.2, 0, 0.5, 0.3
317  };
318  //D2: flat array with the values of D2 applied to basis functions. Multi-index is (P,F,K)
319  double basisD2[] = {
320  0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, \
321  -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, \
322  0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, \
323  -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, \
324  0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, \
325  0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, \
326  -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, \
327  0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, \
328  0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, \
329  0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, \
330  0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, \
331  0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, \
332  0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, \
333  0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, \
334  0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, \
335  0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, \
336  -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, \
337  0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, \
338  -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, \
339  0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, \
340  0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, \
341  -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, \
342  0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0
343  };
344  try{
345 
346  // Dimensions for the output arrays:
347  int numFields = wedgeBasis.getCardinality();
348  int numPoints = wedgeNodes.dimension(0);
349  int spaceDim = wedgeBasis.getBaseCellTopology().getDimension();
350  int D2Cardin = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
351 
352  // Generic array for values, grads, curls, etc. that will be properly sized before each call
353  FieldContainer<double> vals;
354 
355  // Check VALUE of basis functions: resize vals to rank-2 container:
356  vals.resize(numFields, numPoints);
357  wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE);
358  for (int i = 0; i < numFields; i++) {
359  for (int j = 0; j < numPoints; j++) {
360  int l = i + j * numFields;
361  if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
362  errorFlag++;
363  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
364 
365  // Output the multi-index of the value where the error is:
366  *outStream << " At multi-index { ";
367  *outStream << i << " ";*outStream << j << " ";
368  *outStream << "} computed value: " << vals(i,j)
369  << " but reference value: " << basisValues[l] << "\n";
370  }
371  }
372  }
373 
374  // Check GRAD of basis function: resize vals to rank-3 container
375  vals.resize(numFields, numPoints, spaceDim);
376  wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD);
377  for (int i = 0; i < numFields; i++) {
378  for (int j = 0; j < numPoints; j++) {
379  for (int k = 0; k < spaceDim; k++) {
380  int l = k + i * spaceDim + j * spaceDim * numFields;
381  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
382  errorFlag++;
383  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
384 
385  // Output the multi-index of the value where the error is:
386  *outStream << " At multi-index { ";
387  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
388  *outStream << "} computed grad component: " << vals(i,j,k)
389  << " but reference grad component: " << basisGrads[l] << "\n";
390  }
391  }
392  }
393  }
394 
395  // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
396  wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D1);
397  for (int i = 0; i < numFields; i++) {
398  for (int j = 0; j < numPoints; j++) {
399  for (int k = 0; k < spaceDim; k++) {
400  int l = k + i * spaceDim + j * spaceDim * numFields;
401  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
402  errorFlag++;
403  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
404 
405  // Output the multi-index of the value where the error is:
406  *outStream << " At multi-index { ";
407  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
408  *outStream << "} computed D1 component: " << vals(i,j,k)
409  << " but reference D1 component: " << basisGrads[l] << "\n";
410  }
411  }
412  }
413  }
414 
415  // Check D2 of basis function
416  vals.resize(numFields, numPoints, D2Cardin);
417  wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D2);
418  for (int i = 0; i < numFields; i++) {
419  for (int j = 0; j < numPoints; j++) {
420  for (int k = 0; k < D2Cardin; k++) {
421  int l = k + i * D2Cardin + j * D2Cardin * numFields;
422  if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
423  errorFlag++;
424  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
425 
426  // Output the multi-index of the value where the error is:
427  *outStream << " At multi-index { ";
428  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
429  *outStream << "} computed D2 component: " << vals(i,j,k)
430  << " but reference D2 component: " << basisD2[l] << "\n";
431  }
432  }
433  }
434  }
435 
436  // Check all higher derivatives - must be zero.
437  for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
438 
439  // The last dimension is the number of kth derivatives and needs to be resized for every Dk
440  int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
441  vals.resize(numFields, numPoints, DkCardin);
442 
443  wedgeBasis.getValues(vals, wedgeNodes, op);
444  for (int i = 0; i < vals.size(); i++) {
445  if (std::abs(vals[i]) > INTREPID_TOL) {
446  errorFlag++;
447  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
448 
449  // Get the multi-index of the value where the error is and the operator order
450  std::vector<int> myIndex;
451  vals.getMultiIndex(myIndex,i);
452  int ord = Intrepid::getOperatorOrder(op);
453  *outStream << " At multi-index { ";
454  for(int j = 0; j < vals.rank(); j++) {
455  *outStream << myIndex[j] << " ";
456  }
457  *outStream << "} computed D"<< ord <<" component: " << vals[i]
458  << " but reference D" << ord << " component: 0 \n";
459  }
460  }
461  }
462  }
463 
464  // Catch unexpected errors
465  catch (std::logic_error err) {
466  *outStream << err.what() << "\n\n";
467  errorFlag = -1000;
468  };
469 
470  if (errorFlag != 0)
471  std::cout << "End Result: TEST FAILED\n";
472  else
473  std::cout << "End Result: TEST PASSED\n";
474 
475  // reset format state of std::cout
476  std::cout.copyfmt(oldFormatState);
477 
478  return errorFlag;
479 }
Intrepid::FieldContainer::getMultiIndex
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.
Definition: Intrepid_FieldContainerDef.hpp:863
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::getCardinality
virtual int getCardinality() const
Returns cardinality of the basis.
Definition: Intrepid_BasisDef.hpp:100
Intrepid::Basis::getAllDofTags
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Definition: Intrepid_BasisDef.hpp:89
Intrepid::Basis_HGRAD_WEDGE_C1_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Wedge cell.
Definition: Intrepid_HGRAD_WEDGE_C1_FEMDef.hpp:96
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid::FieldContainer::rank
int rank() const
Return rank of the FieldContainer = number of indices used to tag the multi-indexed value.
Definition: Intrepid_FieldContainerDef.hpp:594
Intrepid::Basis_HGRAD_WEDGE_C1_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Wedge cell.
Definition: Intrepid_HGRAD_WEDGE_C1_FEM.hpp:85
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::getDkCardinality
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
Definition: Intrepid_Utils.cpp:400
Intrepid::FieldContainer< double >
Intrepid::getOperatorOrder
int getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
Definition: Intrepid_Utils.cpp:196
main
int main(int argc, char *argv[])
Performs a code-code comparison to FIAT for Nedelec bases on tets (values and curls)
Definition: test_01.cpp:65
Intrepid::EOperator
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
Definition: Intrepid_Types.hpp:206
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::FieldContainer::size
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
Definition: Intrepid_FieldContainerDef.hpp:601
Intrepid::Basis::getDofTag
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
Definition: Intrepid_BasisDef.hpp:78
Generated by   doxygen 1.8.16