Intrepid
test_01.cpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
49#include "Intrepid_HGRAD_HEX_C1_FEM.hpp"
50#include "Teuchos_oblackholestream.hpp"
51#include "Teuchos_RCP.hpp"
52#include "Teuchos_GlobalMPISession.hpp"
53
54using namespace std;
55using 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
71int 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
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
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
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
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}
Header file for utility class to provide multidimensional containers.
int getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
virtual int getDofOrdinal(const int subcDim, const int subcOrd, const int subcDofOrd)
DoF tag to ordinal lookup.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
virtual int getCardinality() const
Returns cardinality of the basis.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
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.