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 
50 #include "Teuchos_oblackholestream.hpp"
51 #include "Teuchos_RCP.hpp"
52 #include "Teuchos_GlobalMPISession.hpp"
53 #include "Intrepid_PointTools.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 (const std::logic_error & err) { \
65  ++throwCounter; \
66  *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
67  *outStream << err.what() << '\n'; \
68  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
69  }; \
70 }
71 
72 int main(int argc, char *argv[]) {
73 
74  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
75 
76  // This little trick lets us print to std::cout only if
77  // a (dummy) command-line argument is provided.
78  int iprint = argc - 1;
79  Teuchos::RCP<std::ostream> outStream;
80  Teuchos::oblackholestream bhs; // outputs nothing
81  if (iprint > 0)
82  outStream = Teuchos::rcp(&std::cout, false);
83  else
84  outStream = Teuchos::rcp(&bhs, false);
85 
86  // Save the format state of the original std::cout.
87  Teuchos::oblackholestream oldFormatState;
88  oldFormatState.copyfmt(std::cout);
89 
90  *outStream \
91  << "===============================================================================\n" \
92  << "| |\n" \
93  << "| Unit Test (Basis_HGRAD_QUAD_Cn_FEM) |\n" \
94  << "| |\n" \
95  << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
96  << "| 2) Basis values for VALUE, GRAD, CURL, and Dk operators |\n" \
97  << "| |\n" \
98  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
99  << "| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
100  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
101  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
102  << "| |\n" \
103  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
104  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
105  << "| |\n" \
106  << "===============================================================================\n"\
107  << "| TEST 1: Basis creation, exception testing |\n"\
108  << "===============================================================================\n";
109 
110  // Define basis and error flag
111  // get points for basis
112  const int deg=2;
113  shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
114  FieldContainer<double> pts(PointTools::getLatticeSize(line,deg),1);
115  PointTools::getLattice<double,FieldContainer<double> >(pts,line,deg);
116 
117  Basis_HGRAD_QUAD_Cn_FEM<double, FieldContainer<double> > quadBasis(deg,deg,pts,pts);
118  int errorFlag = 0;
119 
120  // Initialize throw counter for exception testing
121  int nException = 0;
122  int throwCounter = 0;
123 
124  // Array with the 4 vertices, 4 edge midpoints, center of the reference QUAD and a random point.
125  FieldContainer<double> quadNodes(10, 2);
126  quadNodes(0,0) = -1.0; quadNodes(0,1) = -1.0;
127  quadNodes(1,0) = 1.0; quadNodes(1,1) = -1.0;
128  quadNodes(2,0) = 1.0; quadNodes(2,1) = 1.0;
129  quadNodes(3,0) = -1.0; quadNodes(3,1) = 1.0;
130  // edge midpoints
131  quadNodes(4,0) = 0.0; quadNodes(4,1) = -1.0;
132  quadNodes(5,0) = 1.0; quadNodes(5,1) = 0.0;
133  quadNodes(6,0) = 0.0; quadNodes(6,1) = 1.0;
134  quadNodes(7,0) = -1.0; quadNodes(7,1) = 0.0;
135  // center & random point
136  quadNodes(8,0) = 0.0; quadNodes(8,1) = 0.0;
137  quadNodes(9,0) =1./3.; quadNodes(9,1) =-3./5.;
138 
139  // Generic array for the output values; needs to be properly resized depending on the operator type
141 
142  try{
143  // exception #1: DIV cannot be applied to scalar functions
144  // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
145  vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0));
146  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException );
147 
148  // Exceptions 2-6: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
149  // getDofTag() to access invalid array elements thereby causing bounds check exception
150  // exception #2
151  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException );
152  // exception #3
153  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException );
154  // exception #4
155  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,0), throwCounter, nException );
156  // exception #5
157  INTREPID_TEST_COMMAND( quadBasis.getDofTag(10), throwCounter, nException );
158  // exception #6
159  INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException );
160 
161 #ifdef HAVE_INTREPID_DEBUG
162  // Exceptions 7- test exception handling with incorrectly dimensioned input/output arrays
163  // exception #7: input points array must be of rank-2
164  FieldContainer<double> badPoints1(4, 5, 3);
165  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
166 
167  // exception #8 dimension 1 in the input point array must equal space dimension of the cell
168  FieldContainer<double> badPoints2(4, quadBasis.getBaseCellTopology().getDimension() + 1);
169  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
170 
171  // exception #9 output values must be of rank-2 for OPERATOR_VALUE
172  FieldContainer<double> badVals1(4, 3, 1);
173  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
174 
175  // exception #10 output values must be of rank-3 for OPERATOR_GRAD
176  FieldContainer<double> badVals2(4, 3);
177  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_GRAD), throwCounter, nException );
178 
179  // exception #11 output values must be of rank-3 for OPERATOR_CURL
180  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_CURL), throwCounter, nException );
181 
182  // exception #12 output values must be of rank-3 for OPERATOR_D2
183  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_D2), throwCounter, nException );
184 
185  // exception #13 incorrect 0th dimension of output array (must equal number of basis functions)
186  FieldContainer<double> badVals3(quadBasis.getCardinality() + 1, quadNodes.dimension(0));
187  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException );
188 
189  // exception #14 incorrect 1st dimension of output array (must equal number of points in quadNodes)
190  FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0) + 1);
191  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException );
192 
193  // exception #15: incorrect 2nd dimension of output array (must equal the space dimension)
194  FieldContainer<double> badVals5(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() + 1);
195  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals5, quadNodes, OPERATOR_GRAD), throwCounter, nException );
196 
197  // exception #16: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D)
198  FieldContainer<double> badVals6(quadBasis.getCardinality(), quadNodes.dimension(0), 40);
199  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D2), throwCounter, nException );
200 
201  // exception #17: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D)
202  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D3), throwCounter, nException );
203 #endif
204 
205  }
206  catch (const std::logic_error & err) {
207  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
208  *outStream << err.what() << '\n';
209  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
210  errorFlag = -1000;
211  };
212 
213  // Check if number of thrown exceptions matches the one we expect
214  if (throwCounter != nException) {
215  errorFlag++;
216  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
217  }
218 
219  *outStream \
220  << "\n"
221  << "===============================================================================\n"\
222  << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
223  << "===============================================================================\n";
224 
225  try{
226  std::vector<std::vector<int> > allTags = quadBasis.getAllDofTags();
227 
228  // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
229  for (unsigned i = 0; i < allTags.size(); i++) {
230  int bfOrd = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
231 
232  std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
233  if( !( (myTag[0] == allTags[i][0]) &&
234  (myTag[1] == allTags[i][1]) &&
235  (myTag[2] == allTags[i][2]) &&
236  (myTag[3] == allTags[i][3]) ) ) {
237  errorFlag++;
238  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
239  *outStream << " getDofOrdinal( {"
240  << allTags[i][0] << ", "
241  << allTags[i][1] << ", "
242  << allTags[i][2] << ", "
243  << allTags[i][3] << "}) = " << bfOrd <<" but \n";
244  *outStream << " getDofTag(" << bfOrd << ") = { "
245  << myTag[0] << ", "
246  << myTag[1] << ", "
247  << myTag[2] << ", "
248  << myTag[3] << "}\n";
249  }
250  }
251 
252  // Now do the same but loop over basis functions
253  for( int bfOrd = 0; bfOrd < quadBasis.getCardinality(); bfOrd++) {
254  std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
255  int myBfOrd = quadBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
256  if( bfOrd != myBfOrd) {
257  errorFlag++;
258  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
259  *outStream << " getDofTag(" << bfOrd << ") = { "
260  << myTag[0] << ", "
261  << myTag[1] << ", "
262  << myTag[2] << ", "
263  << myTag[3] << "} but getDofOrdinal({"
264  << myTag[0] << ", "
265  << myTag[1] << ", "
266  << myTag[2] << ", "
267  << myTag[3] << "} ) = " << myBfOrd << "\n";
268  }
269  }
270  }
271  catch (const std::logic_error & err){
272  *outStream << err.what() << "\n\n";
273  errorFlag = -1000;
274  };
275 
276  *outStream \
277  << "\n"
278  << "===============================================================================\n"\
279  << "| TEST 3: correctness of basis function values |\n"\
280  << "===============================================================================\n";
281 
282  outStream -> precision(20);
283 
284  // VALUE: Correct basis values in (F,P) format:
285  double basisValues[] = {
286  1, 0, 0, 0, 0, 0, 0, 0, 0, -0.05333333333333334, \
287  0, 0, 0, 0, 1, 0, 0, 0, 0, 0.4266666666666667, \
288  0, 1, 0, 0, 0, 0, 0, 0, 0, 0.1066666666666667, \
289  0, 0, 0, 0, 0, 0, 0, 1, 0, -0.07111111111111112 , \
290  0, 0, 0, 0, 0, 0, 0, 0, 1, 0.5688888888888890, \
291  0, 0, 0, 0, 0, 1, 0, 0, 0, 0.1422222222222222 ,\
292  0, 0, 0, 1, 0, 0, 0, 0, 0, 0.01333333333333333, \
293  0, 0, 0, 0, 0, 0, 1, 0, 0, -0.1066666666666667, \
294  0, 0, 1, 0, 0, 0, 0, 0, 0, -0.02666666666666666 };
295 
296  // FIXME HERE: needs to be reordered.
297 
298  // GRAD and D1: Correct gradients and D1 in (F,P,D) format
299  // 9 basis functions, each evaluated at 10 points, with two
300  // components at each point.
301  // that looks like 10 per to me.
302  double basisGrads[] = {
303  //
304  -1.500000000000000, -1.500000000000000, 0.5000000000000000, 0, 0, 0, 0, 0.5000000000000000, -0.5000000000000000, 0, \
305  0, 0, 0, 0, 0, -0.5000000000000000, 0, 0, -0.08000000000000002, 0.1222222222222222, \
306  //
307  2.000000000000000, 0, -2.000000000000000, 0, 0, 0, 0, 0, 0, -1.500000000000000, \
308  0, 0, 0, 0.5000000000000000, 0, 0, 0, -0.5000000000000000, -0.3199999999999999, -0.9777777777777779, \
309  //
310  -0.5000000000000000, 0, 1.500000000000000, -1.500000000000000, 0, 0.5000000000000000, 0, 0, 0.5000000000000000, 0, \
311  0, -0.5000000000000000, 0, 0, 0, 0, 0, 0, 0.3999999999999999, -0.2444444444444444, \
312  //
313  0, 2.0, 0, 0, 0, 0, 0, -2.000000000000000, 0, 0, \
314  0.5000000000000000, 0, 0, 0, -1.50, 0, -0.50, 0, -0.1066666666666667, -0.1333333333333333, \
315  //
316  0, 0, 0, 0, 0, 0, 0, 0, 0, 2.0,\
317  -2.00, 0, 0, -2.0, 2.0, 0, 0, 0, -0.4266666666666667, 1.066666666666667 , \
318  //
319  0, 0, 0, 2.000000000000000, 0, -2.000000000000000, 0, 0, 0, 0, \
320  1.5, 0, 0, 0, -0.5, 0, 0.5000000000000000, 0, 0.5333333333333334, 0.2666666666666666 , \
321  //
322  0, -0.5000000000000000, 0, 0, 0.5000000000000000, 0, -1.500000000000000, 1.500000000000000, 0, 0, \
323  0, 0, -0.5000000000000000, 0, 0, 0.5000000000000000, 0, 0, 0.02000000000000000, 0.01111111111111112 , \
324  //
325  0, 0, 0, 0, -2.0, 0, 2.0, 0, 0, -0.50, \
326  0, 0, 0, 1.5, 0, 0, 0, 0.5000000000000000, 0.07999999999999997, -0.08888888888888888, \
327  //
328  0, 0, 0, -0.5000000000000000, 1.500000000000000, 1.500000000000000, -0.5000000000000000, 0, 0, 0, \
329  0, 0.5000000000000000, 0.5000000000000000, 0, 0, 0, 0, 0, -0.09999999999999998, -0.02222222222222221 \
330  //
331  };
332 
333  // D2: Correct multiset of second order partials in (F,P,Dk) format. D2 cardinality = 3 for 2D
334  // 10 quad points and 3 values per point, so
335  // each bf consists of 30 values.
336  double basisD2[] = {
337  1.0, 2.25, 1.0, 1.0, -0.75, 0, 0, 0.25, 0, 0, -0.75, 1.0, 1.0, 0.75, 0, 0, -0.25, 0, 0, -0.25, 0, 0, 0.75, 1.0, 0, 0.25, 0, 0.48, 0.1833333333333334, -0.1111111111111111,
338  //
339  -2.0, -3.0, 0, -2.0, 3.0, 0, 0, -1.0, 0, \
340  0, 1.0, 0, -2.0, 0, 1.0, 0, \
341  1.0, 0, 0, 0, 1.0, 0, -1.0, \
342  0, 0, 0, 1.0, -0.96, 0.7333333333333332, \
343  0.8888888888888890, \
344  //
345  1.0, 0.75, 0, 1.0, -2.25, 1.0, 0, 0.75, 1.0, 0, -0.25, 0, \
346  1.0, -0.75, 0, 0, -0.75, 1.0, 0, 0.25, 0, 0, 0.25, \
347  0, 0, -0.25, 0, 0.48, -0.9166666666666666, 0.2222222222222222,
348  //
349  0, -3.0, -2.0, 0, 1.0, 0, 0, -1.0, 0, 0, 3.0, \
350  -2.0, 0, -1.0, 0, 1.0, 0, 0, 0, 1.0, 0, 1.0, 0, -2.0, \
351  1.0, 0, 0, 0.6400000000000001, -0.2000000000000001, 0.2222222222222222, \
352  //
353  0, 4.0, 0, 0, -4.0, 0, 0, 4.0, 0, 0, -4.0, 0, 0, 0, \
354  -2.0, -2.0, 0, 0, 0, 0, -2.0, -2.0, 0, 0, -2.0, 0, \
355  -2.0, -1.280000000000000, -0.7999999999999998, -1.777777777777778 ,
356  //
357  0, -1.0, 0, 0, 3.0, -2.0, 0, -3.0, -2.0, 0, \
358  1.0, 0, 0, 1.0, 0, 1.0, 0, -2.0, 0, -1.0, 0, 1.0, 0, \
359  0, 1.0, 0, 0, 0.6400000000000001, 1.0, -0.4444444444444444, \
360  //
361  0, 0.75, 1.0, 0, -0.25, 0, 1.0, 0.75, 0, 1.0, -2.25, 1.0, 0, \
362  0.25, 0, 0, 0.25, 0, 1.0, -0.75, 0, 0, -0.75, 1.0, 0, \
363  -0.25, 0, -0.12, 0.01666666666666666, -0.1111111111111111, \
364  //
365  0, -1.0, 0, 0, 1.0, 0, -2.0, -3.0, 0, -2.0, 3.0, 0, 0, 0, 1.0, 0, -1.0, \
366  0, -2.0, 0, 1.0, 0, 1.0, 0, \
367  0, 0, 1.0, 0.24, 0.06666666666666665, 0.8888888888888890, \
368  //
369  0, 0.25, 0, 0, -0.75, 1.0, 1.0, 2.25, 1.0, 1.0, \
370  -0.75, 0, 0, -0.25, 0, 0, 0.75, 1.0, 1.0, \
371  0.75, 0, 0, -0.25, 0, 0, 0.25, 0, -0.12, -0.08333333333333331, 0.2222222222222222 \
372  };
373 
374  //D3: Correct multiset of second order partials in (F,P,Dk) format. D3 cardinality = 4 for 2D
375  double basisD3[] = {
376  0, -1.5, -1.5, 0, 0, -1.5, 0.5, 0, 0, 0.5,
377  0.5, 0, 0, 0.5, -1.5, 0, 0, -1.5, -0.5, 0,
378  0, -0.5, 0.5, 0, 0, 0.5, -0.5, 0, 0, -0.5,
379  -1.5, 0, 0, -0.5, -0.5, 0, 0, -1.1, -0.1666666666666667, 0,
380  //
381  0, 3.0, 2.0, 0, 0, 3.0, -2.0, 0, 0, -1.0,
382  -2.0, 0, 0, -1.0, 2.0, 0, 0, 3.0, 0, 0,
383  0, 1.0, -2.0, 0, 0, -1.0, 0, 0, 0, 1.0,
384  2.0, 0, 0, 1.0, 0, 0, 0, 2.2, -0.6666666666666665, 0,
385  //
386  0, -1.5, -0.5, 0, 0, -1.5, 1.5, 0, 0, 0.5,
387  1.5, 0, 0, 0.5, -0.5, 0, 0, -1.5, 0.5, 0,
388  0, -0.5, 1.5, 0, 0, 0.5, 0.5, 0, 0, -0.5,
389  -0.5, 0, 0, -0.5, 0.5, 0, 0, -1.1, 0.8333333333333333, 0,
390  //
391  0, 2.0, 3.0, 0, 0, 2.0, -1.0, 0, 0, -2.0,
392  -1.0, 0, 0, -2.0, 3.0, 0, 0, 2.0, 1.0, 0,
393  0, 0, -1.0, 0, 0, -2.0, 1.0, 0, 0, 0,
394  3.0, 0, 0, 0, 1.0, 0, 0, 1.2, 0.3333333333333334, 0,
395  //
396  0, -4.0, -4.0, 0, 0, -4.0, 4.0, 0, 0, 4.0,
397  4.0, 0, 0, 4.0, -4.0, 0, 0, -4.0, 0, 0,
398  0, 0, 4.0, 0, 0, 4.0, 0, 0, 0, 0,
399  -4.0, 0, 0, 0, 0, 0, 0, -2.40, 1.333333333333333, 0,
400  //
401  0, 2.0, 1.0, 0, 0, 2.0, -3.0, 0, 0, -2.0,
402  -3.0, 0, 0, -2.0, 1.0, 0, 0, 2.0, -1.0, 0,
403  0, 0, -3.0, 0, 0, -2.0, -1.0, 0, 0, 0,
404  1.0, 0, 0, 0, -1.0, 0, 0, 1.2, -1.666666666666667, 0 ,
405  //
406  0, -0.5, -1.5, 0, 0, -0.5, 0.5, 0, 0, 1.5,
407  0.5, 0, 0, 1.5, -1.5, 0, 0, -0.5, -0.5, 0,
408  0, 0.5, 0.5, 0, 0, 1.5, -0.5, 0, 0, 0.5,
409  -1.5, 0, 0, 0.5, -0.5, 0, 0, -0.09999999999999998, -0.1666666666666667, 0,
410  //
411  0, 1.0, 2.0, 0, 0, 1.0, -2.0, 0, 0, -3.0,
412  -2.0, 0, 0, -3.0, 2.0, 0, 0, 1.0, 0, 0,
413  0, -1.0, -2.0, 0, 0, -3.0, 0, 0, 0, -1.0,
414  2.0, 0, 0, -1.0, 0, 0, 0, 0.2, -0.6666666666666665, 0,
415  //
416  0, -0.5, -0.5, 0, 0, -0.5, 1.5, 0, 0, 1.5,
417  1.5, 0, 0, 1.5, -0.5, 0, 0, -0.5, 0.5, 0,
418  0, 0.5, 1.5, 0, 0, 1.5, 0.5, 0, 0, 0.5,
419  -0.5, 0, 0, 0.5, 0.5, 0, 0, -0.09999999999999998, 0.8333333333333333, 0
420  };
421  //D4: Correct multiset of second order partials in (F,P,Dk) format. D4 cardinality = 5 for 2D
422  double basisD4[] = {
423  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
424  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
425  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
426  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
427  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
428  //
429  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
430  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
431  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
432  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
433  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
434  //
435  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
436  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
437  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
438  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
439  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
440  //
441  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
442  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
443  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
444  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
445  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
446  //
447  0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0,
448  0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0,
449  0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0,
450  0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0,
451  0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0,
452  //
453  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
454  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
455  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
456  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
457  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
458  //
459  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
460  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
461  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
462  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
463  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
464  //
465  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
466  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
467  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
468  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
469  0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
470  //
471  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
472  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
473  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
474  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
475  0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0
476 };
477 
478  try{
479 
480  // Dimensions for the output arrays:
481  int numFields = quadBasis.getCardinality();
482  int numPoints = quadNodes.dimension(0);
483  int spaceDim = quadBasis.getBaseCellTopology().getDimension();
484 
485  // Generic array for values, grads, curls, etc. that will be properly sized before each call
487 
488  // Check VALUE of basis functions: resize vals to rank-2 container:
489  vals.resize(numFields, numPoints);
490  quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE);
491  for (int i = 0; i < numFields; i++) {
492  for (int j = 0; j < numPoints; j++) {
493 
494  // Compute offset for (F,P) container
495  int l = j + i * numPoints;
496  if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
497  errorFlag++;
498  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
499 
500  // Output the multi-index of the value where the error is:
501  *outStream << " At multi-index { ";
502  *outStream << i << " ";*outStream << j << " ";
503  *outStream << "} computed value: " << vals(i,j)
504  << " but reference value: " << basisValues[l] << "\n";
505  }
506  }
507  }
508 
509  // Check GRAD of basis function: resize vals to rank-3 container
510  vals.resize(numFields, numPoints, spaceDim);
511  quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD);
512  for (int i = 0; i < numFields; i++) {
513  for (int j = 0; j < numPoints; j++) {
514  for (int k = 0; k < spaceDim; k++) {
515 
516  // basisGrads is (F,P,D), compute offset:
517  int l = k + j * spaceDim + i * spaceDim * numPoints;
518  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
519  errorFlag++;
520  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
521 
522  // Output the multi-index of the value where the error is:
523  *outStream << " At multi-index { ";
524  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
525  *outStream << "} computed grad component: " << vals(i,j,k)
526  << " but reference grad component: " << basisGrads[l] << "\n";
527  }
528  }
529  }
530  }
531 
532 
533  // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
534  quadBasis.getValues(vals, quadNodes, OPERATOR_D1);
535  for (int i = 0; i < numFields; i++) {
536  for (int j = 0; j < numPoints; j++) {
537  for (int k = 0; k < spaceDim; k++) {
538 
539  // basisGrads is (F,P,D), compute offset:
540  int l = k + j * spaceDim + i * spaceDim * numPoints;
541  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
542  errorFlag++;
543  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
544 
545  // Output the multi-index of the value where the error is:
546  *outStream << " At multi-index { ";
547  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
548  *outStream << "} computed D1 component: " << vals(i,j,k)
549  << " but reference D1 component: " << basisGrads[l] << "\n";
550  }
551  }
552  }
553  }
554 
555 
556  // Check CURL of basis function: resize vals just for illustration!
557  vals.resize(numFields, numPoints, spaceDim);
558  quadBasis.getValues(vals, quadNodes, OPERATOR_CURL);
559  for (int i = 0; i < numFields; i++) {
560  for (int j = 0; j < numPoints; j++) {
561  // We will use "rotated" basisGrads to check CURL: get offsets to extract (u_y, -u_x)
562  int curl_0 = 1 + j * spaceDim + i * spaceDim * numPoints; // position of y-derivative
563  int curl_1 = 0 + j * spaceDim + i * spaceDim * numPoints; // position of x-derivative
564 
565  double curl_value_0 = basisGrads[curl_0];
566  double curl_value_1 =-basisGrads[curl_1];
567  if (std::abs(vals(i,j,0) - curl_value_0) > INTREPID_TOL) {
568  errorFlag++;
569  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
570  // Output the multi-index of the value where the error is:
571  *outStream << " At multi-index { ";
572  *outStream << i << " ";*outStream << j << " ";*outStream << 0 << " ";
573  *outStream << "} computed curl component: " << vals(i,j,0)
574  << " but reference curl component: " << curl_value_0 << "\n";
575  }
576  if (std::abs(vals(i,j,1) - curl_value_1) > INTREPID_TOL) {
577  errorFlag++;
578  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
579  // Output the multi-index of the value where the error is:
580  *outStream << " At multi-index { ";
581  *outStream << i << " ";*outStream << j << " ";*outStream << 1 << " ";
582  *outStream << "} computed curl component: " << vals(i,j,1)
583  << " but reference curl component: " << curl_value_1 << "\n";
584  }
585  }
586  }
587 
588  // Check D2 of basis function
589  int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
590  vals.resize(numFields, numPoints, D2cardinality);
591  quadBasis.getValues(vals, quadNodes, OPERATOR_D2);
592  for (int i = 0; i < numFields; i++) {
593  for (int j = 0; j < numPoints; j++) {
594  for (int k = 0; k < D2cardinality; k++) {
595 
596  // basisD2 is (F,P,Dk), compute offset:
597  int l = k + j * D2cardinality + i * D2cardinality * numPoints;
598  if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
599  errorFlag++;
600  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
601 
602  // Output the multi-index of the value where the error is:
603  *outStream << " At multi-index { ";
604  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
605  *outStream << "} computed D2 component: " << vals(i,j,k)
606  << " but reference D2 component: " << basisD2[l] << "\n";
607  }
608  }
609  }
610  }
611 
612 
613  // Check D3 of basis function
614  int D3cardinality = Intrepid::getDkCardinality(OPERATOR_D3, spaceDim);
615  vals.resize(numFields, numPoints, D3cardinality);
616  quadBasis.getValues(vals, quadNodes, OPERATOR_D3);
617  for (int i = 0; i < numFields; i++) {
618  for (int j = 0; j < numPoints; j++) {
619  for (int k = 0; k < D3cardinality; k++) {
620 
621  // basisD3 is (F,P,Dk), compute offset:
622  int l = k + j * D3cardinality + i * D3cardinality * numPoints;
623  if (std::abs(vals(i,j,k) - basisD3[l]) > INTREPID_TOL) {
624  errorFlag++;
625  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
626 
627  // Output the multi-index of the value where the error is:
628  *outStream << " At multi-index { ";
629  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
630  *outStream << "} computed D3 component: " << vals(i,j,k)
631  << " but reference D3 component: " << basisD2[l] << "\n";
632  }
633  }
634  }
635  }
636 
637  // Check D4 of basis function
638  int D4cardinality = Intrepid::getDkCardinality(OPERATOR_D4, spaceDim);
639  vals.resize(numFields, numPoints, D4cardinality);
640  quadBasis.getValues(vals, quadNodes, OPERATOR_D4);
641  for (int i = 0; i < numFields; i++) {
642  for (int j = 0; j < numPoints; j++) {
643  for (int k = 0; k < D4cardinality; k++) {
644 
645  // basisD4 is (F,P,Dk), compute offset:
646  int l = k + j * D4cardinality + i * D4cardinality * numPoints;
647  if (std::abs(vals(i,j,k) - basisD4[l]) > INTREPID_TOL) {
648  errorFlag++;
649  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
650 
651  // Output the multi-index of the value where the error is:
652  *outStream << " At multi-index { ";
653  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
654  *outStream << "} computed D4 component: " << vals(i,j,k)
655  << " but reference D4 component: " << basisD4[l] << "\n";
656  }
657  }
658  }
659  }
660 
661 
662  // Check all higher derivatives - must be zero.
663  for(EOperator op = OPERATOR_D5; op <= OPERATOR_D6; op++) {
664 
665  // The last dimension is the number of kth derivatives and needs to be resized for every Dk
666  int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
667  vals.resize(numFields, numPoints, DkCardin);
668 
669  quadBasis.getValues(vals, quadNodes, op);
670  for (int i = 0; i < vals.size(); i++) {
671  if (std::abs(vals[i]) > INTREPID_TOL) {
672  errorFlag++;
673  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
674 
675  // Get the multi-index of the value where the error is and the operator order
676  std::vector<int> myIndex;
677  vals.getMultiIndex(myIndex,i);
678  int ord = Intrepid::getOperatorOrder(op);
679  *outStream << " At multi-index { ";
680  for(int j = 0; j < vals.rank(); j++) {
681  *outStream << myIndex[j] << " ";
682  }
683  *outStream << "} computed D"<< ord <<" component: " << vals[i]
684  << " but reference D" << ord << " component: 0 \n";
685  }
686  }
687 }
688  }
689 
690  // Catch unexpected errors
691  catch (const std::logic_error & err) {
692  *outStream << err.what() << "\n\n";
693  errorFlag = -1000;
694  };
695 
696  if (errorFlag != 0)
697  std::cout << "End Result: TEST FAILED\n";
698  else
699  std::cout << "End Result: TEST PASSED\n";
700 
701  // reset format state of std::cout
702  std::cout.copyfmt(oldFormatState);
703 
704  return errorFlag;
705 }
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
Header file for utility class to provide point tools, such as barycentric coordinates, equispaced lattices, and warp-blend point distrubtions.
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 resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
Header file for the Intrepid::HGRAD_QUAD_Cn_FEM class.