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_HEX_C1_FEM.hpp"
50 #include "Teuchos_oblackholestream.hpp"
51 #include "Teuchos_RCP.hpp"
52 #include "Teuchos_GlobalMPISession.hpp"
53 
54 using namespace std;
55 using namespace Intrepid;
56 
57 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
58 { \
59  ++nException; \
60  try { \
61  S ; \
62  } \
63  catch (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_HEX_C1_FEM) |\n" \
93  << "| |\n" \
94  << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95  << "| 2) Basis values for VALUE, GRAD, CURL, and Dk operators |\n" \
96  << "| |\n" \
97  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
98  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
99  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
100  << "| |\n" \
101  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
102  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
103  << "| |\n" \
104  << "===============================================================================\n"\
105  << "| TEST 1: Basis creation, exception testing |\n"\
106  << "===============================================================================\n";
107 
108  // Define basis and error flag
109  Basis_HGRAD_HEX_C1_FEM<double, FieldContainer<double> > hexBasis;
110  int errorFlag = 0;
111 
112  // Initialize throw counter for exception testing
113  int nException = 0;
114  int throwCounter = 0;
115 
116  // Define array containing the 8 vertices of the reference HEX, its center and 6 face centers
117  FieldContainer<double> hexNodes(15, 3);
118  hexNodes(0,0) = -1.0; hexNodes(0,1) = -1.0; hexNodes(0,2) = -1.0;
119  hexNodes(1,0) = 1.0; hexNodes(1,1) = -1.0; hexNodes(1,2) = -1.0;
120  hexNodes(2,0) = 1.0; hexNodes(2,1) = 1.0; hexNodes(2,2) = -1.0;
121  hexNodes(3,0) = -1.0; hexNodes(3,1) = 1.0; hexNodes(3,2) = -1.0;
122 
123  hexNodes(4,0) = -1.0; hexNodes(4,1) = -1.0; hexNodes(4,2) = 1.0;
124  hexNodes(5,0) = 1.0; hexNodes(5,1) = -1.0; hexNodes(5,2) = 1.0;
125  hexNodes(6,0) = 1.0; hexNodes(6,1) = 1.0; hexNodes(6,2) = 1.0;
126  hexNodes(7,0) = -1.0; hexNodes(7,1) = 1.0; hexNodes(7,2) = 1.0;
127 
128  hexNodes(8,0) = 0.0; hexNodes(8,1) = 0.0; hexNodes(8,2) = 0.0;
129 
130  hexNodes(9,0) = 1.0; hexNodes(9,1) = 0.0; hexNodes(9,2) = 0.0;
131  hexNodes(10,0)= -1.0; hexNodes(10,1)= 0.0; hexNodes(10,2)= 0.0;
132 
133  hexNodes(11,0)= 0.0; hexNodes(11,1)= 1.0; hexNodes(11,2)= 0.0;
134  hexNodes(12,0)= 0.0; hexNodes(12,1)= -1.0; hexNodes(12,2)= 0.0;
135 
136  hexNodes(13,0)= 0.0; hexNodes(13,1)= 0.0; hexNodes(13,2)= 1.0;
137  hexNodes(14,0)= 0.0; hexNodes(14,1)= 0.0; hexNodes(14,2)= -1.0;
138 
139 
140  // Generic array for the output values; needs to be properly resized depending on the operator type
141  FieldContainer<double> vals;
142 
143  try{
144  // exception #1: CURL cannot be applied to scalar functions in 3D
145  // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
146  vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0), 4 );
147  INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_CURL), throwCounter, nException );
148 
149  // exception #2: DIV cannot be applied to scalar functions in 3D
150  // resize vals to rank-2 container with dimensions (num. basis functions, num. points)
151  vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0) );
152  INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_DIV), throwCounter, nException );
153 
154  // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
155  // getDofTag() to access invalid array elements thereby causing bounds check exception
156  // exception #3
157  INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(3,0,0), throwCounter, nException );
158  // exception #4
159  INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(1,1,1), throwCounter, nException );
160  // exception #5
161  INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(0,4,1), throwCounter, nException );
162  // exception #6
163  INTREPID_TEST_COMMAND( hexBasis.getDofTag(8), throwCounter, nException );
164  // exception #7
165  INTREPID_TEST_COMMAND( hexBasis.getDofTag(-1), throwCounter, nException );
166 
167 #ifdef HAVE_INTREPID_DEBUG
168  // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays
169  // exception #8: input points array must be of rank-2
170  FieldContainer<double> badPoints1(4, 5, 3);
171  INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
172 
173  // exception #9 dimension 1 in the input point array must equal space dimension of the cell
174  FieldContainer<double> badPoints2(4, 2);
175  INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
176 
177  // exception #10 output values must be of rank-2 for OPERATOR_VALUE
178  FieldContainer<double> badVals1(4, 3, 1);
179  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
180 
181  // exception #11 output values must be of rank-3 for OPERATOR_GRAD
182  FieldContainer<double> badVals2(4, 3);
183  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_GRAD), throwCounter, nException );
184 
185  // exception #12 output values must be of rank-3 for OPERATOR_D1
186  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_D1), throwCounter, nException );
187 
188  // exception #13 output values must be of rank-3 for OPERATOR_D2
189  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_D2), throwCounter, nException );
190 
191  // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
192  FieldContainer<double> badVals3(hexBasis.getCardinality() + 1, hexNodes.dimension(0));
193  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException );
194 
195  // exception #15 incorrect 1st dimension of output array (must equal number of points)
196  FieldContainer<double> badVals4(hexBasis.getCardinality(), hexNodes.dimension(0) + 1);
197  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException );
198 
199  // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
200  FieldContainer<double> badVals5(hexBasis.getCardinality(), hexNodes.dimension(0), 4);
201  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals5, hexNodes, OPERATOR_GRAD), throwCounter, nException );
202 
203  // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D)
204  FieldContainer<double> badVals6(hexBasis.getCardinality(), hexNodes.dimension(0), 40);
205  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals6, hexNodes, OPERATOR_D2), throwCounter, nException );
206 
207  // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D)
208  FieldContainer<double> badVals7(hexBasis.getCardinality(), hexNodes.dimension(0), 50);
209  INTREPID_TEST_COMMAND( hexBasis.getValues(badVals7, hexNodes, OPERATOR_D3), throwCounter, nException );
210 #endif
211 
212  }
213  catch (const std::logic_error & err) {
214  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
215  *outStream << err.what() << '\n';
216  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
217  errorFlag = -1000;
218  };
219 
220  // Check if number of thrown exceptions matches the one we expect
221  // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match.
222  if (throwCounter != nException) {
223  errorFlag++;
224  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
225  }
226 
227  *outStream \
228  << "\n"
229  << "===============================================================================\n"\
230  << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
231  << "===============================================================================\n";
232 
233  try{
234  std::vector<std::vector<int> > allTags = hexBasis.getAllDofTags();
235 
236  // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
237  for (unsigned i = 0; i < allTags.size(); i++) {
238  int bfOrd = hexBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
239 
240  std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
241  if( !( (myTag[0] == allTags[i][0]) &&
242  (myTag[1] == allTags[i][1]) &&
243  (myTag[2] == allTags[i][2]) &&
244  (myTag[3] == allTags[i][3]) ) ) {
245  errorFlag++;
246  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
247  *outStream << " getDofOrdinal( {"
248  << allTags[i][0] << ", "
249  << allTags[i][1] << ", "
250  << allTags[i][2] << ", "
251  << allTags[i][3] << "}) = " << bfOrd <<" but \n";
252  *outStream << " getDofTag(" << bfOrd << ") = { "
253  << myTag[0] << ", "
254  << myTag[1] << ", "
255  << myTag[2] << ", "
256  << myTag[3] << "}\n";
257  }
258  }
259 
260  // Now do the same but loop over basis functions
261  for( int bfOrd = 0; bfOrd < hexBasis.getCardinality(); bfOrd++) {
262  std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
263  int myBfOrd = hexBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
264  if( bfOrd != myBfOrd) {
265  errorFlag++;
266  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
267  *outStream << " getDofTag(" << bfOrd << ") = { "
268  << myTag[0] << ", "
269  << myTag[1] << ", "
270  << myTag[2] << ", "
271  << myTag[3] << "} but getDofOrdinal({"
272  << myTag[0] << ", "
273  << myTag[1] << ", "
274  << myTag[2] << ", "
275  << myTag[3] << "} ) = " << myBfOrd << "\n";
276  }
277  }
278  }
279  catch (const std::logic_error & err){
280  *outStream << err.what() << "\n\n";
281  errorFlag = -1000;
282  };
283 
284  *outStream \
285  << "\n"
286  << "===============================================================================\n"\
287  << "| TEST 3: correctness of basis function values |\n"\
288  << "===============================================================================\n";
289 
290  outStream -> precision(20);
291 
292  // VALUE: Each row gives the 8 correct basis set values at an evaluation point
293  double basisValues[] = {
294  // bottom 4 vertices
295  1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
296  0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
297  0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
298  0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
299  // top 4 vertices
300  0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
301  0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
302  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
303  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
304  // center {0, 0, 0}
305  0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
306  // faces { 1, 0, 0} and {-1, 0, 0}
307  0.0, 0.25, 0.25, 0.0, 0.0, 0.25, 0.25, 0.0,
308  0.25, 0.0, 0.0, 0.25, 0.25, 0.0, 0.0, 0.25,
309  // faces { 0, 1, 0} and { 0,-1, 0}
310  0.0, 0.0, 0.25, 0.25, 0.0, 0.0, 0.25, 0.25,
311  0.25, 0.25, 0.0, 0.0, 0.25, 0.25, 0.0, 0.0,
312  // faces {0, 0, 1} and {0, 0, -1}
313  0.0, 0.0, 0.0, 0.0, 0.25, 0.25, 0.25, 0.25,
314  0.25, 0.25, 0.25, 0.25, 0.0, 0.0, 0.0, 0.0,
315  };
316 
317  // GRAD and D1: each row gives the 3x8 correct values of the gradients of the 8 basis functions
318  double basisGrads[] = {
319  // points 0-3
320  -0.5,-0.5,-0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
321  -0.5, 0.0, 0.0, 0.5,-0.5,-0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
322  0.0, 0.0, 0.0, 0.0,-0.5, 0.0, 0.5, 0.5,-0.5, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0,
323  0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, -0.5, 0.5,-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5,
324  // points 4-7
325  0.0, 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.5,-0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0,
326  0.0, 0.0, 0.0, 0.0, 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.5, 0.0, 0.0, 0.5,-0.5, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0,
327  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.5, 0.0, 0.5, 0.5, 0.5, -0.5, 0.0, 0.0,
328  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.5, 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, -0.5, 0.5, 0.5,
329  // point 8
330  -0.125,-0.125,-0.125, 0.125,-0.125,-0.125, 0.125, 0.125,-0.125, \
331  -0.125, 0.125,-0.125, -0.125,-0.125, 0.125, 0.125,-0.125, 0.125, \
332  0.125, 0.125, 0.125, -0.125, 0.125, 0.125,
333  // point 9
334  -0.125, 0.0, 0.0, 0.125,-0.25, -0.25, 0.125, 0.25, -0.25, -0.125, 0.0, 0.0, \
335  -0.125, 0.0, 0.0, 0.125,-0.25, 0.25, 0.125, 0.25, 0.25, -0.125, 0.0, 0.0,
336  // point 10
337  -0.125,-0.25, -0.25, 0.125, 0.0, 0.0, 0.125, 0.0, 0.0, -0.125, 0.25, -0.25,\
338  -0.125,-0.25, 0.25, 0.125, 0.0, 0.0, 0.125, 0.0, 0.0, -0.125, 0.25, 0.25,
339  // point 11
340  0.0, -0.125, 0.0, 0.0, -0.125, 0.0, 0.25, 0.125,-0.25, -0.25, 0.125,-0.25,\
341  0.0, -0.125, 0.0, 0.0, -0.125, 0.0, 0.25, 0.125, 0.25, -0.25, 0.125, 0.25,
342  // point 12
343  -0.25, -0.125,-0.25, 0.25, -0.125,-0.25, 0.0, 0.125, 0.0, 0.0, 0.125, 0.0, \
344  -0.25, -0.125, 0.25, 0.25, -0.125, 0.25, 0.0, 0.125, 0.0, 0.0, 0.125, 0.0,
345  // point 13
346  0.0, 0.0, -0.125, 0.0, 0.0, -0.125, 0.0, 0.0, -0.125, 0.0, 0.0, -0.125, \
347  -0.25, -0.25, 0.125, 0.25, -0.25, 0.125, 0.25, 0.25, 0.125, -0.25, 0.25, 0.125,
348  // point 14
349  -0.25, -0.25, -0.125, 0.25, -0.25, -0.125, 0.25, 0.25, -0.125, -0.25, 0.25, -0.125, \
350  0.0, 0.0, 0.125, 0.0, 0.0, 0.125, 0.0, 0.0, 0.125, 0.0, 0.0, 0.125
351  };
352 
353  //D2: flat array with the values of D2 applied to basis functions. Multi-index is (P,F,K)
354  double basisD2[] = {
355  // point 0
356  0, 0.25, 0.25, 0, 0.25, 0, 0, -0.25, -0.25, 0, 0., 0, 0, 0.25, 0., 0, \
357  0., 0, 0, -0.25, 0., 0, -0.25, 0, 0, 0., -0.25, 0, -0.25, 0, 0, 0., \
358  0.25, 0, 0., 0, 0, 0., 0., 0, 0., 0, 0, 0., 0., 0, 0.25, 0., \
359  // point 1
360  0, 0.25, 0.25, 0, 0., 0, 0, -0.25, -0.25, 0, 0.25, 0, 0, 0.25, 0., 0, \
361  -0.25, 0, 0, -0.25, 0., 0, 0., 0, 0, 0., -0.25, 0, 0., 0, 0, 0., \
362  0.25, 0, -0.25, 0, 0, 0., 0., 0, 0.25, 0, 0, 0., 0., 0, 0., 0., \
363  // Point 2
364  0, 0.25, 0., 0, 0., 0, 0, -0.25, 0., 0, 0.25, 0, 0, 0.25, -0.25, 0, \
365  -0.25, 0, 0, -0.25, 0.25, 0, 0., 0, 0, 0., 0., 0, 0., 0, 0, 0., 0., \
366  0, -0.25, 0, 0, 0., 0.25, 0, 0.25, 0, 0, 0., -0.25, 0, 0., 0., \
367  // Point 3
368  0, 0.25, 0., 0, 0.25, 0, 0, -0.25, 0., 0, 0., 0, 0, 0.25, -0.25, 0, \
369  0., 0, 0, -0.25, 0.25, 0, -0.25, 0, 0, 0., 0., 0, -0.25, 0, 0, 0., \
370  0., 0, 0., 0, 0, 0., 0.25, 0, 0., 0, 0, 0., -0.25, 0, 0.25, 0.,\
371  // Point 4
372  0, 0., 0.25, 0, 0.25, 0, 0, 0., -0.25, 0, 0., 0, 0, 0., 0., 0, 0., 0, \
373  0, 0., 0., 0, -0.25, 0, 0, 0.25, -0.25, 0, -0.25, 0, 0, -0.25, 0.25, \
374  0, 0., 0, 0, 0.25, 0., 0, 0., 0, 0, -0.25, 0., 0, 0.25, 0., \
375  // Point 5
376  0, 0., 0.25, 0, 0., 0, 0, 0., -0.25, 0, 0.25, 0, 0, 0., 0., 0, -0.25, \
377  0, 0, 0., 0., 0, 0., 0, 0, 0.25, -0.25, 0, 0., 0, 0, -0.25, 0.25, 0, \
378  -0.25, 0, 0, 0.25, 0., 0, 0.25, 0, 0, -0.25, 0., 0, 0., 0., \
379  // Point 6
380  0, 0., 0., 0, 0., 0, 0, 0., 0., 0, 0.25, 0, 0, 0., -0.25, 0, -0.25, \
381  0, 0, 0., 0.25, 0, 0., 0, 0, 0.25, 0., 0, 0., 0, 0, -0.25, 0., 0, \
382  -0.25, 0, 0, 0.25, 0.25, 0, 0.25, 0, 0, -0.25, -0.25, 0, 0., 0., \
383  // Point 7
384  0, 0., 0., 0, 0.25, 0, 0, 0., 0., 0, 0., 0, 0, 0., -0.25, 0, 0., 0, \
385  0, 0., 0.25, 0, -0.25, 0, 0, 0.25, 0., 0, -0.25, 0, 0, -0.25, 0., 0, \
386  0., 0, 0, 0.25, 0.25, 0, 0., 0, 0, -0.25, -0.25, 0, 0.25, 0., \
387  // Point 8
388  0, 0.125, 0.125, 0, 0.125, 0, 0, -0.125, -0.125, 0, 0.125, 0, 0, \
389  0.125, -0.125, 0, -0.125, 0, 0, -0.125, 0.125, 0, -0.125, 0, 0, \
390  0.125, -0.125, 0, -0.125, 0, 0, -0.125, 0.125, 0, -0.125, 0, 0, \
391  0.125, 0.125, 0, 0.125, 0, 0, -0.125, -0.125, 0, 0.125, 0., \
392  // Point 9
393  0, 0.125, 0.125, 0, 0., 0, 0, -0.125, -0.125, 0, 0.25, 0, 0, 0.125, \
394  -0.125, 0, -0.25, 0, 0, -0.125, 0.125, 0, 0., 0, 0, 0.125, -0.125, 0, \
395  0., 0, 0, -0.125, 0.125, 0, -0.25, 0, 0, 0.125, 0.125, 0, 0.25, 0, 0, \
396  -0.125, -0.125, 0, 0., 0., \
397  // Point 10
398  0, 0.125, 0.125, 0, 0.25, 0, 0, -0.125, -0.125, 0, 0., 0, 0, 0.125, \
399  -0.125, 0, 0., 0, 0, -0.125, 0.125, 0, -0.25, 0, 0, 0.125, -0.125, 0, \
400  -0.25, 0, 0, -0.125, 0.125, 0, 0., 0, 0, 0.125, 0.125, 0, 0., 0, 0, \
401  -0.125, -0.125, 0, 0.25, 0., \
402  // Point 11
403  0, 0.125, 0., 0, 0.125, 0, 0, -0.125, 0., 0, 0.125, 0, 0, 0.125, \
404  -0.25, 0, -0.125, 0, 0, -0.125, 0.25, 0, -0.125, 0, 0, 0.125, 0., 0, \
405  -0.125, 0, 0, -0.125, 0., 0, -0.125, 0, 0, 0.125, 0.25, 0, 0.125, 0, \
406  0, -0.125, -0.25, 0, 0.125, 0., \
407  // Point 12
408  0, 0.125, 0.25, 0, 0.125, 0, 0, -0.125, -0.25, 0, 0.125, 0, 0, 0.125, \
409  0., 0, -0.125, 0, 0, -0.125, 0., 0, -0.125, 0, 0, 0.125, -0.25, 0, \
410  -0.125, 0, 0, -0.125, 0.25, 0, -0.125, 0, 0, 0.125, 0., 0, 0.125, 0, \
411  0, -0.125, 0., 0, 0.125, 0., \
412  // Point 13
413  0, 0., 0.125, 0, 0.125, 0, 0, 0., -0.125, 0, 0.125, 0, 0, 0., -0.125, \
414  0, -0.125, 0, 0, 0., 0.125, 0, -0.125, 0, 0, 0.25, -0.125, 0, -0.125, \
415  0, 0, -0.25, 0.125, 0, -0.125, 0, 0, 0.25, 0.125, 0, 0.125, 0, 0, \
416  -0.25, -0.125, 0, 0.125, 0., \
417  // Point 14
418  0, 0.25, 0.125, 0, 0.125, 0, 0, -0.25, -0.125, 0, 0.125, 0, 0, 0.25, \
419  -0.125, 0, -0.125, 0, 0, -0.25, 0.125, 0, -0.125, 0, 0, 0., -0.125, \
420  0, -0.125, 0, 0, 0., 0.125, 0, -0.125, 0, 0, 0., 0.125, 0, 0.125, 0, \
421  0, 0., -0.125, 0, 0.125, 0.
422  };
423 
424  try{
425 
426  // Dimensions for the output arrays:
427  int numFields = hexBasis.getCardinality();
428  int numPoints = hexNodes.dimension(0);
429  int spaceDim = hexBasis.getBaseCellTopology().getDimension();
430  int D2Cardin = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
431 
432  // Generic array for values, grads, curls, etc. that will be properly sized before each call
433  FieldContainer<double> vals;
434 
435  // Check VALUE of basis functions: resize vals to rank-2 container:
436  vals.resize(numFields, numPoints);
437  hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE);
438  for (int i = 0; i < numFields; i++) {
439  for (int j = 0; j < numPoints; j++) {
440  int l = i + j * numFields;
441  if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
442  errorFlag++;
443  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
444 
445  // Output the multi-index of the value where the error is:
446  *outStream << " At multi-index { ";
447  *outStream << i << " ";*outStream << j << " ";
448  *outStream << "} computed value: " << vals(i,j)
449  << " but reference value: " << basisValues[l] << "\n";
450  }
451  }
452  }
453 
454  // Check GRAD of basis function: resize vals to rank-3 container
455  vals.resize(numFields, numPoints, spaceDim);
456  hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD);
457  for (int i = 0; i < numFields; i++) {
458  for (int j = 0; j < numPoints; j++) {
459  for (int k = 0; k < spaceDim; k++) {
460  int l = k + i * spaceDim + j * spaceDim * numFields;
461  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
462  errorFlag++;
463  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
464 
465  // Output the multi-index of the value where the error is:
466  *outStream << " At multi-index { ";
467  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
468  *outStream << "} computed grad component: " << vals(i,j,k)
469  << " but reference grad component: " << basisGrads[l] << "\n";
470  }
471  }
472  }
473  }
474 
475  // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
476  hexBasis.getValues(vals, hexNodes, OPERATOR_D1);
477  for (int i = 0; i < numFields; i++) {
478  for (int j = 0; j < numPoints; j++) {
479  for (int k = 0; k < spaceDim; k++) {
480  int l = k + i * spaceDim + j * spaceDim * numFields;
481  if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
482  errorFlag++;
483  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
484 
485  // Output the multi-index of the value where the error is:
486  *outStream << " At multi-index { ";
487  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
488  *outStream << "} computed D1 component: " << vals(i,j,k)
489  << " but reference D1 component: " << basisGrads[l] << "\n";
490  }
491  }
492  }
493  }
494 
495 
496  // Check D2 of basis function
497  vals.resize(numFields, numPoints, D2Cardin);
498  hexBasis.getValues(vals, hexNodes, OPERATOR_D2);
499  for (int i = 0; i < numFields; i++) {
500  for (int j = 0; j < numPoints; j++) {
501  for (int k = 0; k < D2Cardin; k++) {
502  int l = k + i * D2Cardin + j * D2Cardin * numFields;
503  if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
504  errorFlag++;
505  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
506 
507  // Output the multi-index of the value where the error is:
508  *outStream << " At multi-index { ";
509  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
510  *outStream << "} computed D2 component: " << vals(i,j,k)
511  << " but reference D2 component: " << basisD2[l] << "\n";
512  }
513  }
514  }
515  }
516 
517  // Check all higher derivatives - must be zero.
518  for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
519 
520  // The last dimension is the number of kth derivatives and needs to be resized for every Dk
521  int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
522  vals.resize(numFields, numPoints, DkCardin);
523 
524  hexBasis.getValues(vals, hexNodes, op);
525  for (int i = 0; i < vals.size(); i++) {
526  if (std::abs(vals[i]) > INTREPID_TOL) {
527  errorFlag++;
528  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
529 
530  // Get the multi-index of the value where the error is and the operator order
531  std::vector<int> myIndex;
532  vals.getMultiIndex(myIndex,i);
533  int ord = Intrepid::getOperatorOrder(op);
534  *outStream << " At multi-index { ";
535  for(int j = 0; j < vals.rank(); j++) {
536  *outStream << myIndex[j] << " ";
537  }
538  *outStream << "} computed D"<< ord <<" component: " << vals[i]
539  << " but reference D" << ord << " component: 0 \n";
540  }
541  }
542  }
543  }
544 
545  // Catch unexpected errors
546  catch (const std::logic_error & err) {
547  *outStream << err.what() << "\n\n";
548  errorFlag = -1000;
549  };
550 
551  *outStream \
552  << "\n"
553  << "===============================================================================\n"\
554  << "| TEST 4: correctness of DoF locations |\n"\
555  << "===============================================================================\n";
556 
557  try{
558  Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
559  Teuchos::rcp(new Basis_HGRAD_HEX_C1_FEM<double, FieldContainer<double> >);
560  Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
561  Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
562 
563  FieldContainer<double> cvals;
564  FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality());
565 
566  // Check exceptions.
567 #ifdef HAVE_INTREPID_DEBUG
568  cvals.resize(1,2,3);
569  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
570  cvals.resize(3,2);
571  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
572  cvals.resize(8,2);
573  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
574 #endif
575  cvals.resize(8,3);
576  INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
577  // Check if number of thrown exceptions matches the one we expect
578  if (throwCounter != nException) {
579  errorFlag++;
580  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
581  }
582 
583  // Check mathematical correctness.
584  basis->getValues(bvals, cvals, OPERATOR_VALUE);
585  char buffer[120];
586  for (int i=0; i<bvals.dimension(0); i++) {
587  for (int j=0; j<bvals.dimension(1); j++) {
588  if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) {
589  errorFlag++;
590  sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), bvals(i,j), 0.0);
591  *outStream << buffer;
592  }
593  else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) {
594  errorFlag++;
595  sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), bvals(i,j), 1.0);
596  *outStream << buffer;
597  }
598  }
599  }
600 
601  }
602  catch (const std::logic_error & err){
603  *outStream << err.what() << "\n\n";
604  errorFlag = -1000;
605  };
606 
607  if (errorFlag != 0)
608  std::cout << "End Result: TEST FAILED\n";
609  else
610  std::cout << "End Result: TEST PASSED\n";
611 
612  // reset format state of std::cout
613  std::cout.copyfmt(oldFormatState);
614 
615  return errorFlag;
616 }
Intrepid::FieldContainer::size
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
Definition: Intrepid_FieldContainerDef.hpp:601
Intrepid::Basis::getDofTag
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
Definition: Intrepid_BasisDef.hpp:78
Intrepid::Basis::getCardinality
virtual int getCardinality() const
Returns cardinality of the basis.
Definition: Intrepid_BasisDef.hpp:100
Intrepid::DofCoordsInterface
This is an interface class for bases whose degrees of freedom can be associated with spatial location...
Definition: Intrepid_Basis.hpp:353
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::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_HGRAD_HEX_C1_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell...
Definition: Intrepid_HGRAD_HEX_C1_FEM.hpp:91
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
Intrepid::Basis_HGRAD_HEX_C1_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
Definition: Intrepid_HGRAD_HEX_C1_FEMDef.hpp:100
Generated by   doxygen 1.8.5