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_HGRAD_WEDGE_I2_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_WEDGE_I2_FEM) |\n" \
93  << "| |\n" \
94  << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95  << "| 2) Basis values for VALUE, GRAD, 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_I2_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  // Nodes of Wedge<15>: vertices, edge midpoints
117  FieldContainer<double> wedgeNodes(15, 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.5; wedgeNodes(6,1) = 0.0; wedgeNodes(6,2) = -1.0;
126  wedgeNodes(7,0) = 0.5; wedgeNodes(7,1) = 0.5; wedgeNodes(7,2) = -1.0;
127  wedgeNodes(8,0) = 0.0; wedgeNodes(8,1) = 0.5; wedgeNodes(8,2) = -1.0;
128  wedgeNodes(9,0) = 0.0; wedgeNodes(9,1) = 0.0; wedgeNodes(9,2) = 0.0;
129  wedgeNodes(10,0)= 1.0; wedgeNodes(10,1)= 0.0; wedgeNodes(10,2)= 0.0;
130  wedgeNodes(11,0)= 0.0; wedgeNodes(11,1)= 1.0; wedgeNodes(11,2)= 0.0;
131  wedgeNodes(12,0)= 0.5; wedgeNodes(12,1)= 0.0; wedgeNodes(12,2)= 1.0;
132  wedgeNodes(13,0)= 0.5; wedgeNodes(13,1)= 0.5; wedgeNodes(13,2)= 1.0;
133  wedgeNodes(14,0)= 0.0; wedgeNodes(14,1)= 0.5; wedgeNodes(14,2)= 1.0;
134 
135 
136  // Generic array for the output values; needs to be properly resized depending on the operator type
137  FieldContainer<double> vals;
138 
139  try{
140  // exception #1: CURL cannot be applied to scalar functions
141  // resize vals to rank-3 container with dimensions (num. points, num. basis functions)
142  vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 );
143  INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
144 
145  // exception #2: DIV cannot be applied to scalar functions
146  // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
147  vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) );
148  INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
149 
150  // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
151  // getDofTag() to access invalid array elements thereby causing bounds check exception
152  // exception #3
153  INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException );
154  // exception #4
155  INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException );
156  // exception #5
157  INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,9,0), throwCounter, nException );
158  // exception #6
159  INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(15), throwCounter, nException );
160  // exception #7
161  INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException );
162 
163 #ifdef HAVE_INTREPID_DEBUG
164  // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays
165  // exception #8: input points array must be of rank-2
166  FieldContainer<double> badPoints1(4, 5, 3);
167  INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
168 
169  // exception #9 dimension 1 in the input point array must equal space dimension of the cell
170  FieldContainer<double> badPoints2(4, wedgeBasis.getBaseCellTopology().getDimension() + 1);
171  INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
172 
173  // exception #10 output values must be of rank-2 for OPERATOR_VALUE
174  FieldContainer<double> badVals1(4, 3, 1);
175  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
176 
177  // exception #11 output values must be of rank-3 for OPERATOR_GRAD
178  FieldContainer<double> badVals2(4, 3);
179  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
180 
181  // exception #12 output values must be of rank-3 for OPERATOR_D1
182  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D1), throwCounter, nException );
183 
184  // exception #13 output values must be of rank-3 for OPERATOR_D2
185  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D2), throwCounter, nException );
186 
187  // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
188  FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0));
189  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
190 
191  // exception #15 incorrect 1st dimension of output array (must equal number of points)
192  FieldContainer<double> badVals4(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1);
193  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
194 
195  // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
196  FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), wedgeBasis.getBaseCellTopology().getDimension() - 1);
197  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
198 
199  // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D)
200  FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 40);
201  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D2), throwCounter, nException );
202 
203  // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D)
204  INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D3), throwCounter, nException );
205 #endif
206 
207  }
208  catch (const std::logic_error & err) {
209  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
210  *outStream << err.what() << '\n';
211  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
212  errorFlag = -1000;
213  };
214 
215  // Check if number of thrown exceptions matches the one we expect - 18
216  if (throwCounter != nException) {
217  errorFlag++;
218  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
219  }
220 
221  *outStream \
222  << "\n"
223  << "===============================================================================\n"\
224  << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
225  << "===============================================================================\n";
226 
227  try{
228  std::vector<std::vector<int> > allTags = wedgeBasis.getAllDofTags();
229 
230  // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
231  for (unsigned i = 0; i < allTags.size(); i++) {
232  int bfOrd = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
233 
234  std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
235  if( !( (myTag[0] == allTags[i][0]) &&
236  (myTag[1] == allTags[i][1]) &&
237  (myTag[2] == allTags[i][2]) &&
238  (myTag[3] == allTags[i][3]) ) ) {
239  errorFlag++;
240  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
241  *outStream << " getDofOrdinal( {"
242  << allTags[i][0] << ", "
243  << allTags[i][1] << ", "
244  << allTags[i][2] << ", "
245  << allTags[i][3] << "}) = " << bfOrd <<" but \n";
246  *outStream << " getDofTag(" << bfOrd << ") = { "
247  << myTag[0] << ", "
248  << myTag[1] << ", "
249  << myTag[2] << ", "
250  << myTag[3] << "}\n";
251  }
252  }
253 
254  // Now do the same but loop over basis functions
255  for( int bfOrd = 0; bfOrd < wedgeBasis.getCardinality(); bfOrd++) {
256  std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
257  int myBfOrd = wedgeBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
258  if( bfOrd != myBfOrd) {
259  errorFlag++;
260  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
261  *outStream << " getDofTag(" << bfOrd << ") = { "
262  << myTag[0] << ", "
263  << myTag[1] << ", "
264  << myTag[2] << ", "
265  << myTag[3] << "} but getDofOrdinal({"
266  << myTag[0] << ", "
267  << myTag[1] << ", "
268  << myTag[2] << ", "
269  << myTag[3] << "} ) = " << myBfOrd << "\n";
270  }
271  }
272  }
273  catch (const std::logic_error & err){
274  *outStream << err.what() << "\n\n";
275  errorFlag = -1000;
276  };
277 
278  *outStream \
279  << "\n"
280  << "===============================================================================\n"\
281  << "| TEST 3: correctness of basis function values |\n"\
282  << "===============================================================================\n";
283 
284  outStream -> precision(20);
285 
286  // VALUE: correct basis function values in (F,P) format
287  double basisValues[] = {
288  1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
289  0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
290  0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
291  0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
292  0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
293  0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
294  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
295  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
296  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
297  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
298  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,\
299  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,\
300  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,\
301  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,\
302  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0};
303 
304  // GRAD, D1, D2, and D3 test values are stored in files due to their large size
305  std::string fileName;
306  std::ifstream dataFile;
307 
308  // GRAD and D1 values are stored in (F,P,D) format in a data file. Read file and do the test
309  std::vector<double> basisGrads; // Flat array for the gradient values.
310 
311  fileName = "./testdata/WEDGE_I2_GradVals.dat";
312  dataFile.open(fileName.c_str());
313  TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
314  ">>> ERROR (HGRAD_WEDGE_I2/test01): could not open GRAD values data file, test aborted.");
315  while (!dataFile.eof() ){
316  double temp;
317  string line; // string for one line of input file
318  std::getline(dataFile, line); // get next line from file
319  stringstream data_line(line); // convert to stringstream
320  while(data_line >> temp){ // extract value from line
321  basisGrads.push_back(temp); // push into vector
322  }
323  }
324  // It turns out that just closing and then opening the ifstream variable does not reset it
325  // and subsequent open() command fails. One fix is to explicitely clear the ifstream, or
326  // scope the variables.
327  dataFile.close();
328  dataFile.clear();
329 
330 
331  //D2: flat array with the values of D2 applied to basis functions. Multi-index is (F,P,D2cardinality)
332  std::vector<double> basisD2;
333 
334  fileName = "./testdata/WEDGE_I2_D2Vals.dat";
335  dataFile.open(fileName.c_str());
336  TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
337  ">>> ERROR (HGRAD_WEDGE_I2/test01): could not open D2 values data file, test aborted.");
338 
339  while (!dataFile.eof() ){
340  double temp;
341  string line; // string for one line of input file
342  std::getline(dataFile, line); // get next line from file
343  stringstream data_line(line); // convert to stringstream
344  while(data_line >> temp){ // extract value from line
345  basisD2.push_back(temp); // push into vector
346  }
347  }
348  dataFile.close();
349  dataFile.clear();
350 
351 
352  //D3: flat array with the values of D3 applied to basis functions. Multi-index is (F,P,D3cardinality)
353  std::vector<double> basisD3;
354 
355  fileName = "./testdata/WEDGE_I2_D3Vals.dat";
356  dataFile.open(fileName.c_str());
357  TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
358  ">>> ERROR (HGRAD_WEDGE_I2/test01): could not open D3 values data file, test aborted.");
359 
360  while (!dataFile.eof() ){
361  double temp;
362  string line; // string for one line of input file
363  std::getline(dataFile, line); // get next line from file
364  stringstream data_line(line); // convert to stringstream
365  while(data_line >> temp){ // extract value from line
366  basisD3.push_back(temp); // push into vector
367  }
368  }
369  dataFile.close();
370  dataFile.clear();
371 
372  try{
373 
374  // Dimensions for the output arrays:
375  int numFields = wedgeBasis.getCardinality();
376  int numPoints = wedgeNodes.dimension(0);
377  int spaceDim = wedgeBasis.getBaseCellTopology().getDimension();
378 
379  // Generic array for values, grads, curls, etc. that will be properly sized before each call
380  FieldContainer<double> vals;
381 
382  // Check VALUE of basis functions: resize vals to rank-2 container:
383  vals.resize(numFields, numPoints);
384  wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE);
385  for (int i = 0; i < numFields; i++) {
386  for (int j = 0; j < numPoints; j++) {
387  int l = i + j * numFields;
388  if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
389  errorFlag++;
390  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
391 
392  // Output the multi-index of the value where the error is:
393  *outStream << " At multi-index { ";
394  *outStream << i << " ";*outStream << j << " ";
395  *outStream << "} computed value: " << vals(i,j)
396  << " but reference value: " << basisValues[l] << "\n";
397  }
398  }
399  }
400 
401 
402 
403  // Check GRAD of basis function: resize vals to rank-3 container
404  vals.resize(numFields, numPoints, spaceDim);
405  wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD);
406  for (int i = 0; i < numFields; i++) {
407  for (int j = 0; j < numPoints; j++) {
408  for (int k = 0; k < spaceDim; k++) {
409 
410  // basisGrads is (F,P,D), compute offset:
411  int l = k + j * spaceDim + i * spaceDim * numPoints;
412  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
413  errorFlag++;
414  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
415 
416  // Output the multi-index of the value where the error is:
417  *outStream << " At multi-index { ";
418  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
419  *outStream << "} computed grad component: " << vals(i,j,k)
420  << " but reference grad component: " << basisGrads[l] << "\n";
421  }
422  }
423  }
424  }
425 
426  // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
427  wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D1);
428  for (int i = 0; i < numFields; i++) {
429  for (int j = 0; j < numPoints; j++) {
430  for (int k = 0; k < spaceDim; k++) {
431 
432  // basisGrads is (F,P,D), compute offset:
433  int l = k + j * spaceDim + i * spaceDim * numPoints;
434  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
435  errorFlag++;
436  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
437 
438  // Output the multi-index of the value where the error is:
439  *outStream << " At multi-index { ";
440  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
441  *outStream << "} computed D1 component: " << vals(i,j,k)
442  << " but reference D1 component: " << basisGrads[l] << "\n";
443  }
444  }
445  }
446  }
447 
448 
449  // Check D2 of basis function
450  int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
451  vals.resize(numFields, numPoints, D2cardinality);
452  wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D2);
453  for (int i = 0; i < numFields; i++) {
454  for (int j = 0; j < numPoints; j++) {
455  for (int k = 0; k < D2cardinality; k++) {
456 
457  // basisD2 is (F,P,Dk), compute offset:
458  int l = k + j * D2cardinality + i * D2cardinality * numPoints;
459  if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
460  errorFlag++;
461  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
462 
463  // Output the multi-index of the value where the error is:
464  *outStream << " At multi-index { ";
465  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
466  *outStream << "} computed D2 component: " << vals(i,j,k)
467  << " but reference D2 component: " << basisD2[l] << "\n";
468  }
469  }
470  }
471  }
472 
473 
474  // Check D3 of basis function
475  int D3cardinality = Intrepid::getDkCardinality(OPERATOR_D3, spaceDim);
476  vals.resize(numFields, numPoints, D3cardinality);
477  wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D3);
478 
479  for (int i = 0; i < numFields; i++) {
480  for (int j = 0; j < numPoints; j++) {
481  for (int k = 0; k < D3cardinality; k++) {
482 
483  // basisD3 is (F,P,Dk), compute offset:
484  int l = k + j * D3cardinality + i * D3cardinality * numPoints;
485  if (std::abs(vals(i,j,k) - basisD3[l]) > INTREPID_TOL) {
486  errorFlag++;
487  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
488 
489  // Output the multi-index of the value where the error is:
490  *outStream << " At multi-index { ";
491  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
492  *outStream << "} computed D3 component: " << vals(i,j,k)
493  << " but reference D3 component: " << basisD3[l] << "\n";
494  }
495  }
496  }
497  }
498 
499 
500  // Check all higher derivatives - must be zero.
501  for(EOperator op = OPERATOR_D4; op < OPERATOR_MAX; op++) {
502 
503  // The last dimension is the number of kth derivatives and needs to be resized for every Dk
504  int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
505  vals.resize(numFields, numPoints, DkCardin);
506 
507  wedgeBasis.getValues(vals, wedgeNodes, op);
508  for (int i = 0; i < vals.size(); i++) {
509  if (std::abs(vals[i]) > INTREPID_TOL) {
510  errorFlag++;
511  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
512 
513  // Get the multi-index of the value where the error is and the operator order
514  std::vector<int> myIndex;
515  vals.getMultiIndex(myIndex,i);
516  int ord = Intrepid::getOperatorOrder(op);
517  *outStream << " At multi-index { ";
518  for(int j = 0; j < vals.rank(); j++) {
519  *outStream << myIndex[j] << " ";
520  }
521  *outStream << "} computed D"<< ord <<" component: " << vals[i]
522  << " but reference D" << ord << " component: 0 \n";
523  }
524  }
525  }
526  }
527 
528  // Catch unexpected errors
529  catch (const std::logic_error & err) {
530  *outStream << err.what() << "\n\n";
531  errorFlag = -1000;
532  };
533 
534  if (errorFlag != 0)
535  std::cout << "End Result: TEST FAILED\n";
536  else
537  std::cout << "End Result: TEST PASSED\n";
538 
539  // reset format state of std::cout
540  std::cout.copyfmt(oldFormatState);
541 
542  return errorFlag;
543 }
Intrepid::FieldContainer::size
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
Definition: Intrepid_FieldContainerDef.hpp:601
Intrepid_HGRAD_WEDGE_I2_FEM.hpp
Header file for the Intrepid::HGRAD_WEDGE_I2_FEM class.
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::Basis_HGRAD_WEDGE_I2_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_I2_FEMDef.hpp:106
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::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::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_I2_FEM
Implementation of an H(grad)-compatible FEM basis of degree 2 on Wedge cell.
Definition: Intrepid_HGRAD_WEDGE_I2_FEM.hpp:110
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::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
Generated by   doxygen 1.8.5