Intrepid
test_01.cpp
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38 // Denis Ridzal (dridzal@sandia.gov), or
39 // Kara Peterson (kjpeter@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
49 #include "Intrepid_HGRAD_QUAD_C2_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_QUAD_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 (kjpeter@sandia.gov). |\n" \
100  << "| |\n" \
101  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
102  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
103  << "| |\n" \
104  << "===============================================================================\n"\
105  << "| TEST 1: Basis creation, exception testing |\n"\
106  << "===============================================================================\n";
107 
108  // Define basis and error flag
110  int errorFlag = 0;
111 
112  // Initialize throw counter for exception testing
113  int nException = 0;
114  int throwCounter = 0;
115 
116  // Array with the 4 vertices, 4 edge midpoints, center of the reference QUAD and a random point.
117  FieldContainer<double> quadNodes(10, 2);
118  quadNodes(0,0) = -1.0; quadNodes(0,1) = -1.0;
119  quadNodes(1,0) = 1.0; quadNodes(1,1) = -1.0;
120  quadNodes(2,0) = 1.0; quadNodes(2,1) = 1.0;
121  quadNodes(3,0) = -1.0; quadNodes(3,1) = 1.0;
122  // edge midpoints
123  quadNodes(4,0) = 0.0; quadNodes(4,1) = -1.0;
124  quadNodes(5,0) = 1.0; quadNodes(5,1) = 0.0;
125  quadNodes(6,0) = 0.0; quadNodes(6,1) = 1.0;
126  quadNodes(7,0) = -1.0; quadNodes(7,1) = 0.0;
127  // center & random point
128  quadNodes(8,0) = 0.0; quadNodes(8,1) = 0.0;
129  quadNodes(9,0) =1./3.; quadNodes(9,1) =-3./5.;
130 
131 
132  // Generic array for the output values; needs to be properly resized depending on the operator type
134 
135  try{
136  // exception #1: DIV cannot be applied to scalar functions
137  // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
138  vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0));
139  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException );
140 
141  // Exceptions 2-6: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
142  // getDofTag() to access invalid array elements thereby causing bounds check exception
143  // exception #2
144  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException );
145  // exception #3
146  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException );
147  // exception #4
148  INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,0), throwCounter, nException );
149  // exception #5
150  INTREPID_TEST_COMMAND( quadBasis.getDofTag(10), throwCounter, nException );
151  // exception #6
152  INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException );
153 
154 #ifdef HAVE_INTREPID_DEBUG
155  // Exceptions 7- test exception handling with incorrectly dimensioned input/output arrays
156  // exception #7: input points array must be of rank-2
157  FieldContainer<double> badPoints1(4, 5, 3);
158  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
159 
160  // exception #8 dimension 1 in the input point array must equal space dimension of the cell
161  FieldContainer<double> badPoints2(4, quadBasis.getBaseCellTopology().getDimension() + 1);
162  INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
163 
164  // exception #9 output values must be of rank-2 for OPERATOR_VALUE
165  FieldContainer<double> badVals1(4, 3, 1);
166  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
167 
168  // exception #10 output values must be of rank-3 for OPERATOR_GRAD
169  FieldContainer<double> badVals2(4, 3);
170  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_GRAD), throwCounter, nException );
171 
172  // exception #11 output values must be of rank-3 for OPERATOR_CURL
173  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_CURL), throwCounter, nException );
174 
175  // exception #12 output values must be of rank-3 for OPERATOR_D2
176  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_D2), throwCounter, nException );
177 
178  // exception #13 incorrect 0th dimension of output array (must equal number of basis functions)
179  FieldContainer<double> badVals3(quadBasis.getCardinality() + 1, quadNodes.dimension(0));
180  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException );
181 
182  // exception #14 incorrect 1st dimension of output array (must equal number of points in quadNodes)
183  FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0) + 1);
184  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException );
185 
186  // exception #15: incorrect 2nd dimension of output array (must equal the space dimension)
187  FieldContainer<double> badVals5(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() + 1);
188  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals5, quadNodes, OPERATOR_GRAD), throwCounter, nException );
189 
190  // exception #16: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D)
191  FieldContainer<double> badVals6(quadBasis.getCardinality(), quadNodes.dimension(0), 40);
192  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D2), throwCounter, nException );
193 
194  // exception #17: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D)
195  INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D3), throwCounter, nException );
196 #endif
197 
198  }
199  catch (const std::logic_error & err) {
200  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
201  *outStream << err.what() << '\n';
202  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
203  errorFlag = -1000;
204  };
205 
206  // Check if number of thrown exceptions matches the one we expect
207  if (throwCounter != nException) {
208  errorFlag++;
209  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
210  }
211 
212  *outStream \
213  << "\n"
214  << "===============================================================================\n"\
215  << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
216  << "===============================================================================\n";
217 
218  try{
219  std::vector<std::vector<int> > allTags = quadBasis.getAllDofTags();
220 
221  // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
222  for (unsigned i = 0; i < allTags.size(); i++) {
223  int bfOrd = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
224 
225  std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
226  if( !( (myTag[0] == allTags[i][0]) &&
227  (myTag[1] == allTags[i][1]) &&
228  (myTag[2] == allTags[i][2]) &&
229  (myTag[3] == allTags[i][3]) ) ) {
230  errorFlag++;
231  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
232  *outStream << " getDofOrdinal( {"
233  << allTags[i][0] << ", "
234  << allTags[i][1] << ", "
235  << allTags[i][2] << ", "
236  << allTags[i][3] << "}) = " << bfOrd <<" but \n";
237  *outStream << " getDofTag(" << bfOrd << ") = { "
238  << myTag[0] << ", "
239  << myTag[1] << ", "
240  << myTag[2] << ", "
241  << myTag[3] << "}\n";
242  }
243  }
244 
245  // Now do the same but loop over basis functions
246  for( int bfOrd = 0; bfOrd < quadBasis.getCardinality(); bfOrd++) {
247  std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
248  int myBfOrd = quadBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
249  if( bfOrd != myBfOrd) {
250  errorFlag++;
251  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
252  *outStream << " getDofTag(" << bfOrd << ") = { "
253  << myTag[0] << ", "
254  << myTag[1] << ", "
255  << myTag[2] << ", "
256  << myTag[3] << "} but getDofOrdinal({"
257  << myTag[0] << ", "
258  << myTag[1] << ", "
259  << myTag[2] << ", "
260  << myTag[3] << "} ) = " << myBfOrd << "\n";
261  }
262  }
263  }
264  catch (const std::logic_error & err){
265  *outStream << err.what() << "\n\n";
266  errorFlag = -1000;
267  };
268 
269  *outStream \
270  << "\n"
271  << "===============================================================================\n"\
272  << "| TEST 3: correctness of basis function values |\n"\
273  << "===============================================================================\n";
274 
275  outStream -> precision(20);
276 
277  // VALUE: Correct basis values in (F,P) format:
278  double basisValues[] = {
279  1.000000000000000, 0, 0, 0, 0, 0, 0, 0, 0, -0.05333333333333334, 0, \
280  1.000000000000000, 0, 0, 0, 0, 0, 0, 0, 0.1066666666666667, 0, 0, \
281  1.000000000000000, 0, 0, 0, 0, 0, 0, -0.02666666666666666, 0, 0, 0, \
282  1.000000000000000, 0, 0, 0, 0, 0, 0.01333333333333333, 0, 0, 0, 0, \
283  1.000000000000000, 0, 0, 0, 0, 0.4266666666666667, 0, 0, 0, 0, 0, \
284  1.000000000000000, 0, 0, 0, 0.1422222222222222, 0, 0, 0, 0, 0, 0, \
285  1.000000000000000, 0, 0, -0.1066666666666667, 0, 0, 0, 0, 0, 0, 0, \
286  1.000000000000000, 0, -0.07111111111111112, 0, 0, 0, 0, 0, 0, 0, 0, \
287  1.000000000000000, 0.5688888888888890};
288 
289  // GRAD and D1: Correct gradients and D1 in (F,P,D) format: each row contains 6x2 values of
290  double basisGrads[] = {
291  -1.500000000000000, -1.500000000000000, 0.5000000000000000, 0, 0, 0, \
292  0, 0.5000000000000000, -0.5000000000000000, 0, 0, 0, 0, 0, 0, \
293  -0.5000000000000000, 0, 0, -0.08000000000000002, 0.1222222222222222, \
294  -0.5000000000000000, 0, 1.500000000000000, -1.500000000000000, 0, \
295  0.5000000000000000, 0, 0, 0.5000000000000000, 0, 0, \
296  -0.5000000000000000, 0, 0, 0, 0, 0, 0, 0.3999999999999999, \
297  -0.2444444444444444, 0, 0, 0, -0.5000000000000000, 1.500000000000000, \
298  1.500000000000000, -0.5000000000000000, 0, 0, 0, 0, \
299  0.5000000000000000, 0.5000000000000000, 0, 0, 0, 0, 0, \
300  -0.09999999999999998, -0.02222222222222221, 0, -0.5000000000000000, \
301  0, 0, 0.5000000000000000, 0, -1.500000000000000, 1.500000000000000, \
302  0, 0, 0, 0, -0.5000000000000000, 0, 0, 0.5000000000000000, 0, 0, \
303  0.02000000000000000, 0.01111111111111112, 2.000000000000000, 0, \
304  -2.000000000000000, 0, 0, 0, 0, 0, 0, -1.500000000000000, 0, 0, 0, \
305  0.5000000000000000, 0, 0, 0, -0.5000000000000000, \
306  -0.3199999999999999, -0.9777777777777779, 0, 0, 0, 2.000000000000000, \
307  0, -2.000000000000000, 0, 0, 0, 0, 1.500000000000000, 0, 0, 0, \
308  -0.5000000000000000, 0, 0.5000000000000000, 0, 0.5333333333333334, \
309  0.2666666666666666, 0, 0, 0, 0, -2.000000000000000, 0, \
310  2.000000000000000, 0, 0, -0.5000000000000000, 0, 0, 0, \
311  1.500000000000000, 0, 0, 0, 0.5000000000000000, 0.07999999999999997, \
312  -0.08888888888888888, 0, 2.000000000000000, 0, 0, 0, 0, 0, \
313  -2.000000000000000, 0, 0, 0.5000000000000000, 0, 0, 0, \
314  -1.500000000000000, 0, -0.5000000000000000, 0, -0.1066666666666667, \
315  -0.1333333333333333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.000000000000000, \
316  -2.000000000000000, 0, 0, -2.000000000000000, 2.000000000000000, 0, \
317  0, 0, -0.4266666666666667, 1.066666666666667};
318 
319  // D2: Correct multiset of second order partials in (F,P,Dk) format. D2 cardinality = 3 for 2D
320  double basisD2[] = {
321  1.000000000000000, 2.250000000000000, 1.000000000000000, \
322  1.000000000000000, -0.7500000000000000, 0, 0, 0.2500000000000000, 0, \
323  0, -0.7500000000000000, 1.000000000000000, 1.000000000000000, \
324  0.7500000000000000, 0, 0, -0.2500000000000000, 0, 0, \
325  -0.2500000000000000, 0, 0, 0.7500000000000000, 1.000000000000000, 0, \
326  0.2500000000000000, 0, 0.4800000000000000, 0.1833333333333334, \
327  -0.1111111111111111, 1.000000000000000, 0.7500000000000000, 0, \
328  1.000000000000000, -2.250000000000000, 1.000000000000000, 0, \
329  0.7500000000000000, 1.000000000000000, 0, -0.2500000000000000, 0, \
330  1.000000000000000, -0.7500000000000000, 0, 0, -0.7500000000000000, \
331  1.000000000000000, 0, 0.2500000000000000, 0, 0, 0.2500000000000000, \
332  0, 0, -0.2500000000000000, 0, 0.4800000000000000, \
333  -0.9166666666666666, 0.2222222222222222, 0, 0.2500000000000000, 0, 0, \
334  -0.7500000000000000, 1.000000000000000, 1.000000000000000, \
335  2.250000000000000, 1.000000000000000, 1.000000000000000, \
336  -0.7500000000000000, 0, 0, -0.2500000000000000, 0, 0, \
337  0.7500000000000000, 1.000000000000000, 1.000000000000000, \
338  0.7500000000000000, 0, 0, -0.2500000000000000, 0, 0, \
339  0.2500000000000000, 0, -0.1200000000000000, -0.08333333333333331, \
340  0.2222222222222222, 0, 0.7500000000000000, 1.000000000000000, 0, \
341  -0.2500000000000000, 0, 1.000000000000000, 0.7500000000000000, 0, \
342  1.000000000000000, -2.250000000000000, 1.000000000000000, 0, \
343  0.2500000000000000, 0, 0, 0.2500000000000000, 0, 1.000000000000000, \
344  -0.7500000000000000, 0, 0, -0.7500000000000000, 1.000000000000000, 0, \
345  -0.2500000000000000, 0, -0.1200000000000000, 0.01666666666666666, \
346  -0.1111111111111111, -2.000000000000000, -3.000000000000000, 0, \
347  -2.000000000000000, 3.000000000000000, 0, 0, -1.000000000000000, 0, \
348  0, 1.000000000000000, 0, -2.000000000000000, 0, 1.000000000000000, 0, \
349  1.000000000000000, 0, 0, 0, 1.000000000000000, 0, -1.000000000000000, \
350  0, 0, 0, 1.000000000000000, -0.9600000000000000, 0.7333333333333332, \
351  0.8888888888888890, 0, -1.000000000000000, 0, 0, 3.000000000000000, \
352  -2.000000000000000, 0, -3.000000000000000, -2.000000000000000, 0, \
353  1.000000000000000, 0, 0, 1.000000000000000, 0, 1.000000000000000, 0, \
354  -2.000000000000000, 0, -1.000000000000000, 0, 1.000000000000000, 0, \
355  0, 1.000000000000000, 0, 0, 0.6400000000000001, 1.000000000000000, \
356  -0.4444444444444444, 0, -1.000000000000000, 0, 0, 1.000000000000000, \
357  0, -2.000000000000000, -3.000000000000000, 0, -2.000000000000000, \
358  3.000000000000000, 0, 0, 0, 1.000000000000000, 0, -1.000000000000000, \
359  0, -2.000000000000000, 0, 1.000000000000000, 0, 1.000000000000000, 0, \
360  0, 0, 1.000000000000000, 0.2400000000000000, 0.06666666666666665, \
361  0.8888888888888890, 0, -3.000000000000000, -2.000000000000000, 0, \
362  1.000000000000000, 0, 0, -1.000000000000000, 0, 0, 3.000000000000000, \
363  -2.000000000000000, 0, -1.000000000000000, 0, 1.000000000000000, 0, \
364  0, 0, 1.000000000000000, 0, 1.000000000000000, 0, -2.000000000000000, \
365  1.000000000000000, 0, 0, 0.6400000000000001, -0.2000000000000001, \
366  0.2222222222222222, 0, 4.000000000000000, 0, 0, -4.000000000000000, \
367  0, 0, 4.000000000000000, 0, 0, -4.000000000000000, 0, 0, 0, \
368  -2.000000000000000, -2.000000000000000, 0, 0, 0, 0, \
369  -2.000000000000000, -2.000000000000000, 0, 0, -2.000000000000000, 0, \
370  -2.000000000000000, -1.280000000000000, -0.7999999999999998, \
371  -1.777777777777778};
372 
373  //D3: Correct multiset of second order partials in (F,P,Dk) format. D3 cardinality = 4 for 2D
374  double basisD3[] = {
375  0, -1.500000000000000, -1.500000000000000, 0, 0, -1.500000000000000, \
376  0.5000000000000000, 0, 0, 0.5000000000000000, 0.5000000000000000, 0, \
377  0, 0.5000000000000000, -1.500000000000000, 0, 0, -1.500000000000000, \
378  -0.5000000000000000, 0, 0, -0.5000000000000000, 0.5000000000000000, \
379  0, 0, 0.5000000000000000, -0.5000000000000000, 0, 0, \
380  -0.5000000000000000, -1.500000000000000, 0, 0, -0.5000000000000000, \
381  -0.5000000000000000, 0, 0, -1.100000000000000, -0.1666666666666667, \
382  0, 0, -1.500000000000000, -0.5000000000000000, 0, 0, \
383  -1.500000000000000, 1.500000000000000, 0, 0, 0.5000000000000000, \
384  1.500000000000000, 0, 0, 0.5000000000000000, -0.5000000000000000, 0, \
385  0, -1.500000000000000, 0.5000000000000000, 0, 0, -0.5000000000000000, \
386  1.500000000000000, 0, 0, 0.5000000000000000, 0.5000000000000000, 0, \
387  0, -0.5000000000000000, -0.5000000000000000, 0, 0, \
388  -0.5000000000000000, 0.5000000000000000, 0, 0, -1.100000000000000, \
389  0.8333333333333333, 0, 0, -0.5000000000000000, -0.5000000000000000, \
390  0, 0, -0.5000000000000000, 1.500000000000000, 0, 0, \
391  1.500000000000000, 1.500000000000000, 0, 0, 1.500000000000000, \
392  -0.5000000000000000, 0, 0, -0.5000000000000000, 0.5000000000000000, \
393  0, 0, 0.5000000000000000, 1.500000000000000, 0, 0, 1.500000000000000, \
394  0.5000000000000000, 0, 0, 0.5000000000000000, -0.5000000000000000, 0, \
395  0, 0.5000000000000000, 0.5000000000000000, 0, 0, \
396  -0.09999999999999998, 0.8333333333333333, 0, 0, -0.5000000000000000, \
397  -1.500000000000000, 0, 0, -0.5000000000000000, 0.5000000000000000, 0, \
398  0, 1.500000000000000, 0.5000000000000000, 0, 0, 1.500000000000000, \
399  -1.500000000000000, 0, 0, -0.5000000000000000, -0.5000000000000000, \
400  0, 0, 0.5000000000000000, 0.5000000000000000, 0, 0, \
401  1.500000000000000, -0.5000000000000000, 0, 0, 0.5000000000000000, \
402  -1.500000000000000, 0, 0, 0.5000000000000000, -0.5000000000000000, 0, \
403  0, -0.09999999999999998, -0.1666666666666667, 0, 0, \
404  3.000000000000000, 2.000000000000000, 0, 0, 3.000000000000000, \
405  -2.000000000000000, 0, 0, -1.000000000000000, -2.000000000000000, 0, \
406  0, -1.000000000000000, 2.000000000000000, 0, 0, 3.000000000000000, 0, \
407  0, 0, 1.000000000000000, -2.000000000000000, 0, 0, \
408  -1.000000000000000, 0, 0, 0, 1.000000000000000, 2.000000000000000, 0, \
409  0, 1.000000000000000, 0, 0, 0, 2.200000000000000, \
410  -0.6666666666666665, 0, 0, 2.000000000000000, 1.000000000000000, 0, \
411  0, 2.000000000000000, -3.000000000000000, 0, 0, -2.000000000000000, \
412  -3.000000000000000, 0, 0, -2.000000000000000, 1.000000000000000, 0, \
413  0, 2.000000000000000, -1.000000000000000, 0, 0, 0, \
414  -3.000000000000000, 0, 0, -2.000000000000000, -1.000000000000000, 0, \
415  0, 0, 1.000000000000000, 0, 0, 0, -1.000000000000000, 0, 0, \
416  1.200000000000000, -1.666666666666667, 0, 0, 1.000000000000000, \
417  2.000000000000000, 0, 0, 1.000000000000000, -2.000000000000000, 0, 0, \
418  -3.000000000000000, -2.000000000000000, 0, 0, -3.000000000000000, \
419  2.000000000000000, 0, 0, 1.000000000000000, 0, 0, 0, \
420  -1.000000000000000, -2.000000000000000, 0, 0, -3.000000000000000, 0, \
421  0, 0, -1.000000000000000, 2.000000000000000, 0, 0, \
422  -1.000000000000000, 0, 0, 0, 0.2000000000000000, -0.6666666666666665, \
423  0, 0, 2.000000000000000, 3.000000000000000, 0, 0, 2.000000000000000, \
424  -1.000000000000000, 0, 0, -2.000000000000000, -1.000000000000000, 0, \
425  0, -2.000000000000000, 3.000000000000000, 0, 0, 2.000000000000000, \
426  1.000000000000000, 0, 0, 0, -1.000000000000000, 0, 0, \
427  -2.000000000000000, 1.000000000000000, 0, 0, 0, 3.000000000000000, 0, \
428  0, 0, 1.000000000000000, 0, 0, 1.200000000000000, 0.3333333333333334, \
429  0, 0, -4.000000000000000, -4.000000000000000, 0, 0, \
430  -4.000000000000000, 4.000000000000000, 0, 0, 4.000000000000000, \
431  4.000000000000000, 0, 0, 4.000000000000000, -4.000000000000000, 0, 0, \
432  -4.000000000000000, 0, 0, 0, 0, 4.000000000000000, 0, 0, \
433  4.000000000000000, 0, 0, 0, 0, -4.000000000000000, 0, 0, 0, 0, 0, 0, \
434  -2.400000000000000, 1.333333333333333, 0};
435  //D4: Correct multiset of second order partials in (F,P,Dk) format. D4 cardinality = 5 for 2D
436  double basisD4[] = {
437  0, 0, 1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
438  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
439  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
440  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
441  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
442  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
443  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
444  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
445  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
446  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
447  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
448  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
449  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
450  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
451  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
452  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
453  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
454  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
455  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
456  1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
457  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
458  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
459  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
460  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
461  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
462  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
463  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
464  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
465  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
466  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
467  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
468  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
469  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
470  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
471  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
472  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
473  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
474  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
475  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
476  -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
477  4.000000000000000, 0, 0, 0, 0, 4.000000000000000, 0, 0, 0, 0, \
478  4.000000000000000, 0, 0, 0, 0, 4.000000000000000, 0, 0, 0, 0, \
479  4.000000000000000, 0, 0, 0, 0, 4.000000000000000, 0, 0, 0, 0, \
480  4.000000000000000, 0, 0, 0, 0, 4.000000000000000, 0, 0, 0, 0, \
481  4.000000000000000, 0, 0, 0, 0, 4.000000000000000, 0, 0};
482 
483  try{
484 
485  // Dimensions for the output arrays:
486  int numFields = quadBasis.getCardinality();
487  int numPoints = quadNodes.dimension(0);
488  int spaceDim = quadBasis.getBaseCellTopology().getDimension();
489 
490  // Generic array for values, grads, curls, etc. that will be properly sized before each call
492 
493  // Check VALUE of basis functions: resize vals to rank-2 container:
494  vals.resize(numFields, numPoints);
495  quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE);
496  for (int i = 0; i < numFields; i++) {
497  for (int j = 0; j < numPoints; j++) {
498 
499  // Compute offset for (F,P) container
500  int l = j + i * numPoints;
501  if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
502  errorFlag++;
503  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
504 
505  // Output the multi-index of the value where the error is:
506  *outStream << " At multi-index { ";
507  *outStream << i << " ";*outStream << j << " ";
508  *outStream << "} computed value: " << vals(i,j)
509  << " but reference value: " << basisValues[l] << "\n";
510  }
511  }
512  }
513 
514  // Check GRAD of basis function: resize vals to rank-3 container
515  vals.resize(numFields, numPoints, spaceDim);
516  quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD);
517  for (int i = 0; i < numFields; i++) {
518  for (int j = 0; j < numPoints; j++) {
519  for (int k = 0; k < spaceDim; k++) {
520 
521  // basisGrads is (F,P,D), compute offset:
522  int l = k + j * spaceDim + i * spaceDim * numPoints;
523  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
524  errorFlag++;
525  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
526 
527  // Output the multi-index of the value where the error is:
528  *outStream << " At multi-index { ";
529  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
530  *outStream << "} computed grad component: " << vals(i,j,k)
531  << " but reference grad component: " << basisGrads[l] << "\n";
532  }
533  }
534  }
535  }
536 
537 
538  // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
539  quadBasis.getValues(vals, quadNodes, OPERATOR_D1);
540  for (int i = 0; i < numFields; i++) {
541  for (int j = 0; j < numPoints; j++) {
542  for (int k = 0; k < spaceDim; k++) {
543 
544  // basisGrads is (F,P,D), compute offset:
545  int l = k + j * spaceDim + i * spaceDim * numPoints;
546  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
547  errorFlag++;
548  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
549 
550  // Output the multi-index of the value where the error is:
551  *outStream << " At multi-index { ";
552  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
553  *outStream << "} computed D1 component: " << vals(i,j,k)
554  << " but reference D1 component: " << basisGrads[l] << "\n";
555  }
556  }
557  }
558  }
559 
560 
561  // Check CURL of basis function: resize vals just for illustration!
562  vals.resize(numFields, numPoints, spaceDim);
563  quadBasis.getValues(vals, quadNodes, OPERATOR_CURL);
564  for (int i = 0; i < numFields; i++) {
565  for (int j = 0; j < numPoints; j++) {
566  // We will use "rotated" basisGrads to check CURL: get offsets to extract (u_y, -u_x)
567  int curl_0 = 1 + j * spaceDim + i * spaceDim * numPoints; // position of y-derivative
568  int curl_1 = 0 + j * spaceDim + i * spaceDim * numPoints; // position of x-derivative
569 
570  double curl_value_0 = basisGrads[curl_0];
571  double curl_value_1 =-basisGrads[curl_1];
572  if (std::abs(vals(i,j,0) - curl_value_0) > INTREPID_TOL) {
573  errorFlag++;
574  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
575  // Output the multi-index of the value where the error is:
576  *outStream << " At multi-index { ";
577  *outStream << i << " ";*outStream << j << " ";*outStream << 0 << " ";
578  *outStream << "} computed curl component: " << vals(i,j,0)
579  << " but reference curl component: " << curl_value_0 << "\n";
580  }
581  if (std::abs(vals(i,j,1) - curl_value_1) > INTREPID_TOL) {
582  errorFlag++;
583  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
584  // Output the multi-index of the value where the error is:
585  *outStream << " At multi-index { ";
586  *outStream << i << " ";*outStream << j << " ";*outStream << 1 << " ";
587  *outStream << "} computed curl component: " << vals(i,j,1)
588  << " but reference curl component: " << curl_value_1 << "\n";
589  }
590  }
591  }
592 
593  // Check D2 of basis function
594  int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
595  vals.resize(numFields, numPoints, D2cardinality);
596  quadBasis.getValues(vals, quadNodes, OPERATOR_D2);
597  for (int i = 0; i < numFields; i++) {
598  for (int j = 0; j < numPoints; j++) {
599  for (int k = 0; k < D2cardinality; k++) {
600 
601  // basisD2 is (F,P,Dk), compute offset:
602  int l = k + j * D2cardinality + i * D2cardinality * numPoints;
603  if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
604  errorFlag++;
605  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
606 
607  // Output the multi-index of the value where the error is:
608  *outStream << " At multi-index { ";
609  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
610  *outStream << "} computed D2 component: " << vals(i,j,k)
611  << " but reference D2 component: " << basisD2[l] << "\n";
612  }
613  }
614  }
615  }
616 
617 
618  // Check D3 of basis function
619  int D3cardinality = Intrepid::getDkCardinality(OPERATOR_D3, spaceDim);
620  vals.resize(numFields, numPoints, D3cardinality);
621  quadBasis.getValues(vals, quadNodes, OPERATOR_D3);
622  for (int i = 0; i < numFields; i++) {
623  for (int j = 0; j < numPoints; j++) {
624  for (int k = 0; k < D3cardinality; k++) {
625 
626  // basisD3 is (F,P,Dk), compute offset:
627  int l = k + j * D3cardinality + i * D3cardinality * numPoints;
628  if (std::abs(vals(i,j,k) - basisD3[l]) > INTREPID_TOL) {
629  errorFlag++;
630  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
631 
632  // Output the multi-index of the value where the error is:
633  *outStream << " At multi-index { ";
634  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
635  *outStream << "} computed D3 component: " << vals(i,j,k)
636  << " but reference D3 component: " << basisD2[l] << "\n";
637  }
638  }
639  }
640  }
641 
642 
643  // Check D4 of basis function
644  int D4cardinality = Intrepid::getDkCardinality(OPERATOR_D4, spaceDim);
645  vals.resize(numFields, numPoints, D4cardinality);
646  quadBasis.getValues(vals, quadNodes, OPERATOR_D4);
647  for (int i = 0; i < numFields; i++) {
648  for (int j = 0; j < numPoints; j++) {
649  for (int k = 0; k < D4cardinality; k++) {
650 
651  // basisD4 is (F,P,Dk), compute offset:
652  int l = k + j * D4cardinality + i * D4cardinality * numPoints;
653  if (std::abs(vals(i,j,k) - basisD4[l]) > INTREPID_TOL) {
654  errorFlag++;
655  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
656 
657  // Output the multi-index of the value where the error is:
658  *outStream << " At multi-index { ";
659  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
660  *outStream << "} computed D4 component: " << vals(i,j,k)
661  << " but reference D4 component: " << basisD2[l] << "\n";
662  }
663  }
664  }
665  }
666 
667 
668  // Check all higher derivatives - must be zero.
669  for(EOperator op = OPERATOR_D5; op < OPERATOR_MAX; op++) {
670 
671  // The last dimension is the number of kth derivatives and needs to be resized for every Dk
672  int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
673  vals.resize(numFields, numPoints, DkCardin);
674 
675  quadBasis.getValues(vals, quadNodes, op);
676  for (int i = 0; i < vals.size(); i++) {
677  if (std::abs(vals[i]) > INTREPID_TOL) {
678  errorFlag++;
679  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
680 
681  // Get the multi-index of the value where the error is and the operator order
682  std::vector<int> myIndex;
683  vals.getMultiIndex(myIndex,i);
684  int ord = Intrepid::getOperatorOrder(op);
685  *outStream << " At multi-index { ";
686  for(int j = 0; j < vals.rank(); j++) {
687  *outStream << myIndex[j] << " ";
688  }
689  *outStream << "} computed D"<< ord <<" component: " << vals[i]
690  << " but reference D" << ord << " component: 0 \n";
691  }
692  }
693  }
694  }
695  // Catch unexpected errors
696  catch (const std::logic_error & err) {
697  *outStream << err.what() << "\n\n";
698  errorFlag = -1000;
699  };
700 
701 
702  *outStream \
703  << "\n"
704  << "===============================================================================\n"\
705  << "| TEST 4: correctness of DoF locations |\n"\
706  << "===============================================================================\n";
707 
708  try{
709  Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
710  Teuchos::rcp(new Basis_HGRAD_QUAD_C2_FEM<double, FieldContainer<double> >);
711  Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
712  Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
713 
715  FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality());
716 
717  // Check exceptions.
718 #ifdef HAVE_INTREPID_DEBUG
719  cvals.resize(1,2,3);
720  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
721  cvals.resize(8,2);
722  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
723  cvals.resize(9,1);
724  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
725 #endif
726  cvals.resize(9,2);
727  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
728  // Check if number of thrown exceptions matches the one we expect
729  if (throwCounter != nException) {
730  errorFlag++;
731  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
732  }
733 
734  // Check mathematical correctness.
735  basis->getValues(bvals, cvals, OPERATOR_VALUE);
736  char buffer[120];
737  for (int i=0; i<bvals.dimension(0); i++) {
738  for (int j=0; j<bvals.dimension(1); j++) {
739  if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) {
740  errorFlag++;
741  sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), bvals(i,j), 0.0);
742  *outStream << buffer;
743  }
744  else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) {
745  errorFlag++;
746  sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), bvals(i,j), 1.0);
747  *outStream << buffer;
748  }
749  }
750  }
751 
752  }
753  catch (const std::logic_error & err){
754  *outStream << err.what() << "\n\n";
755  errorFlag = -1000;
756  };
757 
758  if (errorFlag != 0)
759  std::cout << "End Result: TEST FAILED\n";
760  else
761  std::cout << "End Result: TEST PASSED\n";
762 
763  // reset format state of std::cout
764  std::cout.copyfmt(oldFormatState);
765 
766  return errorFlag;
767 }
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.
This is an interface class for bases whose degrees of freedom can be associated with spatial location...
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...
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Quadrilateral cell...
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...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Quadrilateral cell.