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
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_TET_C2_FEM) |\n" \
93 << "| |\n" \
94 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 << "| 2) Basis values for VALUE, GRAD, and Dk operators |\n" \
96 << "| |\n" \
97 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
98 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
99 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
100 << "| |\n" \
101 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
102 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
103 << "| |\n" \
104 << "===============================================================================\n"\
105 << "| TEST 1: Basis creation, exception testing |\n"\
106 << "===============================================================================\n";
107
108 // Define basis and error flag
109 Basis_HGRAD_TET_C2_FEM<double, FieldContainer<double> > tetBasis;
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 10 nodes of Tetrahedron<10>: 4 vertices + 6 edge midpoints
117 FieldContainer<double> tetNodes(10, 3);
118 tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0;
119 tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0;
120 tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0;
121 tetNodes(3,0) = 0.0; tetNodes(3,1) = 0.0; tetNodes(3,2) = 1.0;
122
123 tetNodes(4,0) = 0.5; tetNodes(4,1) = 0.0; tetNodes(4,2) = 0.0;
124 tetNodes(5,0) = 0.5; tetNodes(5,1) = 0.5; tetNodes(5,2) = 0.0;
125 tetNodes(6,0) = 0.0; tetNodes(6,1) = 0.5; tetNodes(6,2) = 0.0;
126 tetNodes(7,0) = 0.0; tetNodes(7,1) = 0.0; tetNodes(7,2) = 0.5;
127 tetNodes(8,0) = 0.5; tetNodes(8,1) = 0.0; tetNodes(8,2) = 0.5;
128 tetNodes(9,0) = 0.0; tetNodes(9,1) = 0.5; tetNodes(9,2) = 0.5;
129
130
131 // Generic array for the output values; needs to be properly resized depending on the operator type
132 FieldContainer<double> vals;
133
134 try{
135 // exception #1: CURL cannot be applied to scalar functions
136 // resize vals to rank-3 container with dimensions (num. points, num. basis functions, arbitrary)
137 vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0), 4 );
138 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException );
139
140 // exception #2: DIV cannot be applied to scalar functions
141 // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
142 vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0) );
143 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_DIV), throwCounter, nException );
144
145 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
146 // getDofTag() to access invalid array elements thereby causing bounds check exception
147 // exception #3
148 INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(3,0,0), throwCounter, nException );
149 // exception #4
150 INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(1,1,1), throwCounter, nException );
151 // exception #5
152 INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(0,4,0), throwCounter, nException );
153 // exception #6
154 INTREPID_TEST_COMMAND( tetBasis.getDofTag(10), throwCounter, nException );
155 // exception #7
156 INTREPID_TEST_COMMAND( tetBasis.getDofTag(-1), throwCounter, nException );
157
158#ifdef HAVE_INTREPID_DEBUG
159 // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays
160 // exception #8: input points array must be of rank-2
161 FieldContainer<double> badPoints1(4, 5, 3);
162 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
163
164 // exception #9 dimension 1 in the input point array must equal space dimension of the cell
165 FieldContainer<double> badPoints2(4, tetBasis.getBaseCellTopology().getDimension() - 1);
166 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
167
168 // exception #10 output values must be of rank-2 for OPERATOR_VALUE
169 FieldContainer<double> badVals1(4, 3, 1);
170 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
171
172 // exception #11 output values must be of rank-3 for OPERATOR_GRAD
173 FieldContainer<double> badVals2(4, 3);
174 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_GRAD), throwCounter, nException );
175
176 // exception #12 output values must be of rank-3 for OPERATOR_D1
177 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_D1), throwCounter, nException );
178
179 // exception #13 output values must be of rank-3 for OPERATOR_D2
180 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_D2), throwCounter, nException );
181
182 // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
183 FieldContainer<double> badVals3(tetBasis.getCardinality() + 1, tetNodes.dimension(0));
184 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
185
186 // exception #15 incorrect 1st dimension of output array (must equal number of points)
187 FieldContainer<double> badVals4(tetBasis.getCardinality(), tetNodes.dimension(0) + 1);
188 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals4, tetNodes, OPERATOR_VALUE), throwCounter, nException );
189
190 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
191 FieldContainer<double> badVals5(tetBasis.getCardinality(), tetNodes.dimension(0), tetBasis.getBaseCellTopology().getDimension() + 1);
192 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals5, tetNodes, OPERATOR_GRAD), throwCounter, nException );
193
194 // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D)
195 FieldContainer<double> badVals6(tetBasis.getCardinality(), tetNodes.dimension(0), 40);
196 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_D1), throwCounter, nException );
197
198 // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D)
199 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_D2), throwCounter, nException );
200#endif
201
202 }
203 catch (const std::logic_error & err) {
204 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
205 *outStream << err.what() << '\n';
206 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
207 errorFlag = -1000;
208 };
209
210 // Check if number of thrown exceptions matches the one we expect
211 if (throwCounter != nException) {
212 errorFlag++;
213 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
214 }
215
216 *outStream \
217 << "\n"
218 << "===============================================================================\n"\
219 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
220 << "===============================================================================\n";
221
222 try{
223 std::vector<std::vector<int> > allTags = tetBasis.getAllDofTags();
224
225 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
226 for (unsigned i = 0; i < allTags.size(); i++) {
227 int bfOrd = tetBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
228
229 std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
230 if( !( (myTag[0] == allTags[i][0]) &&
231 (myTag[1] == allTags[i][1]) &&
232 (myTag[2] == allTags[i][2]) &&
233 (myTag[3] == allTags[i][3]) ) ) {
234 errorFlag++;
235 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
236 *outStream << " getDofOrdinal( {"
237 << allTags[i][0] << ", "
238 << allTags[i][1] << ", "
239 << allTags[i][2] << ", "
240 << allTags[i][3] << "}) = " << bfOrd <<" but \n";
241 *outStream << " getDofTag(" << bfOrd << ") = { "
242 << myTag[0] << ", "
243 << myTag[1] << ", "
244 << myTag[2] << ", "
245 << myTag[3] << "}\n";
246 }
247 }
248
249 // Now do the same but loop over basis functions
250 for( int bfOrd = 0; bfOrd < tetBasis.getCardinality(); bfOrd++) {
251 std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
252 int myBfOrd = tetBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
253 if( bfOrd != myBfOrd) {
254 errorFlag++;
255 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
256 *outStream << " getDofTag(" << bfOrd << ") = { "
257 << myTag[0] << ", "
258 << myTag[1] << ", "
259 << myTag[2] << ", "
260 << myTag[3] << "} but getDofOrdinal({"
261 << myTag[0] << ", "
262 << myTag[1] << ", "
263 << myTag[2] << ", "
264 << myTag[3] << "} ) = " << myBfOrd << "\n";
265 }
266 }
267 }
268 catch (const std::logic_error & err){
269 *outStream << err.what() << "\n\n";
270 errorFlag = -1000;
271 };
272
273 *outStream \
274 << "\n"
275 << "===============================================================================\n"\
276 << "| TEST 3: correctness of basis function values |\n"\
277 << "===============================================================================\n";
278
279 outStream -> precision(20);
280
281 // VALUE: in (F,P) format
282 double basisValues[] = {
283 1.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00000, 0, 0, 0, 0, 0, 0, 0, \
284 0, 0, 0, 1.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00000, 0, 0, 0, 0, \
285 0, 0, 0, 0, 0, 0, 1.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00000, 0, \
286 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
287 1.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00000, 0, 0, 0, 0, 0, 0, 0, \
288 0, 0, 0, 1.00000 };
289
290 // GRAD and D1: in (F,P,D) format
291 double basisGrads[] = {
292 -3.00000, -3.00000, -3.00000, 1.00000, 1.00000, 1.00000, 1.00000, \
293 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, -1.00000, -1.00000, \
294 -1.00000, 1.00000, 1.00000, 1.00000, -1.00000, -1.00000, -1.00000, \
295 -1.00000, -1.00000, -1.00000, 1.00000, 1.00000, 1.00000, 1.00000, \
296 1.00000, 1.00000, -1.00000, 0, 0, 3.00000, 0, 0, -1.00000, 0, 0, \
297 -1.00000, 0, 0, 1.00000, 0, 0, 1.00000, 0, 0, -1.00000, 0, 0, \
298 -1.00000, 0, 0, 1.00000, 0, 0, -1.00000, 0, 0, 0, -1.00000, 0, 0, \
299 -1.00000, 0, 0, 3.00000, 0, 0, -1.00000, 0, 0, -1.00000, 0, 0, \
300 1.00000, 0, 0, 1.00000, 0, 0, -1.00000, 0, 0, -1.00000, 0, 0, \
301 1.00000, 0, 0, 0, -1.00000, 0, 0, -1.00000, 0, 0, -1.00000, 0, 0, \
302 3.00000, 0, 0, -1.00000, 0, 0, -1.00000, 0, 0, -1.00000, 0, 0, \
303 1.00000, 0, 0, 1.00000, 0, 0, 1.00000, 4.00000, 0, 0, -4.00000, \
304 -4.00000, -4.00000, 0, 0, 0, 0, 0, 0, 0, -2.00000, -2.00000, \
305 -2.00000, -2.00000, -2.00000, 2.00000, 0, 0, 2.00000, 0, 0, -2.00000, \
306 -2.00000, -2.00000, 0, 0, 0, 0, 0, 0, 0, 4.00000, 0, 4.00000, 0, 0, \
307 0, 0, 0, 0, 2.00000, 0, 2.00000, 2.00000, 0, 2.00000, 0, 0, 0, 0, 0, \
308 0, 2.00000, 0, 2.00000, 0, 0, 0, 4.00000, 0, 0, 0, 0, -4.00000, \
309 -4.00000, -4.00000, 0, 0, 0, 0, 2.00000, 0, -2.00000, -2.00000, \
310 -2.00000, -2.00000, 0, -2.00000, 0, 2.00000, 0, 0, 0, 0, -2.00000, \
311 -2.00000, -2.00000, 0, 0, 4.00000, 0, 0, 0, 0, 0, 0, -4.00000, \
312 -4.00000, -4.00000, 0, 0, 2.00000, 0, 0, 0, 0, 0, 2.00000, -2.00000, \
313 -2.00000, 0, -2.00000, -2.00000, -2.00000, -2.00000, -2.00000, \
314 -2.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 4.00000, 0, 0, 0, 0, \
315 2.00000, 0, 0, 2.00000, 0, 0, 0, 2.00000, 0, 0, 2.00000, 0, 2.00000, \
316 2.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.00000, 0, 4.00000, 0, 0, 0, \
317 0, 0, 0, 2.00000, 0, 0, 2.00000, 0, 2.00000, 0, 0, 2.00000, 0, 0, \
318 2.00000, 2.00000};
319
320 // D2 values in (F,P, Dk) format
321 double basisD2[]={
322 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
323 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
324 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
325 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
326 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
327 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
328 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
329 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, \
330 4.00000, 4.00000, 4.00000, 4.00000, 4.00000, 0, 0, 0, 0, 0, 4.00000, \
331 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, \
332 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, \
333 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, \
334 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, \
335 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, \
336 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, \
337 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 0, 0, 4.00000, \
338 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, \
339 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, \
340 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, \
341 0, 0, 4.00000, -8.00000, -4.00000, -4.00000, 0, 0, 0, -8.00000, \
342 -4.00000, -4.00000, 0, 0, 0, -8.00000, -4.00000, -4.00000, 0, 0, 0, \
343 -8.00000, -4.00000, -4.00000, 0, 0, 0, -8.00000, -4.00000, -4.00000, \
344 0, 0, 0, -8.00000, -4.00000, -4.00000, 0, 0, 0, -8.00000, -4.00000, \
345 -4.00000, 0, 0, 0, -8.00000, -4.00000, -4.00000, 0, 0, 0, -8.00000, \
346 -4.00000, -4.00000, 0, 0, 0, -8.00000, -4.00000, -4.00000, 0, 0, 0, \
347 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, \
348 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, \
349 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, \
350 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, -4.00000, 0, -8.00000, -4.00000, \
351 0, 0, -4.00000, 0, -8.00000, -4.00000, 0, 0, -4.00000, 0, -8.00000, \
352 -4.00000, 0, 0, -4.00000, 0, -8.00000, -4.00000, 0, 0, -4.00000, 0, \
353 -8.00000, -4.00000, 0, 0, -4.00000, 0, -8.00000, -4.00000, 0, 0, \
354 -4.00000, 0, -8.00000, -4.00000, 0, 0, -4.00000, 0, -8.00000, \
355 -4.00000, 0, 0, -4.00000, 0, -8.00000, -4.00000, 0, 0, -4.00000, 0, \
356 -8.00000, -4.00000, 0, 0, 0, -4.00000, 0, -4.00000, -8.00000, 0, 0, \
357 -4.00000, 0, -4.00000, -8.00000, 0, 0, -4.00000, 0, -4.00000, \
358 -8.00000, 0, 0, -4.00000, 0, -4.00000, -8.00000, 0, 0, -4.00000, 0, \
359 -4.00000, -8.00000, 0, 0, -4.00000, 0, -4.00000, -8.00000, 0, 0, \
360 -4.00000, 0, -4.00000, -8.00000, 0, 0, -4.00000, 0, -4.00000, \
361 -8.00000, 0, 0, -4.00000, 0, -4.00000, -8.00000, 0, 0, -4.00000, 0, \
362 -4.00000, -8.00000, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, \
363 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, \
364 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, \
365 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 0, 0, \
366 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, \
367 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, \
368 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, 0, 0, 0, 4.00000, 0, 0, \
369 0, 0, 0, 4.00000, 0
370 };
371
372
373 try{
374
375 // Dimensions for the output arrays:
376 int numFields = tetBasis.getCardinality();
377 int numPoints = tetNodes.dimension(0);
378 int spaceDim = tetBasis.getBaseCellTopology().getDimension();
379
380 // Generic array for values, grads, curls, etc. that will be properly sized before each call
381 FieldContainer<double> vals;
382
383 // Check VALUE of basis functions: resize vals to rank-2 container:
384 vals.resize(numFields, numPoints);
385 tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE);
386 for (int i = 0; i < numFields; i++) {
387 for (int j = 0; j < numPoints; j++) {
388 int l = i + j * numFields;
389 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
390 errorFlag++;
391 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
392
393 // Output the multi-index of the value where the error is:
394 *outStream << " At multi-index { ";
395 *outStream << i << " ";*outStream << j << " ";
396 *outStream << "} computed value: " << vals(i,j)
397 << " but reference value: " << basisValues[l] << "\n";
398 }
399 }
400 }
401
402 // Check GRAD of basis function: resize vals to rank-3 container
403 vals.resize(numFields, numPoints, spaceDim);
404 tetBasis.getValues(vals, tetNodes, OPERATOR_GRAD);
405 for (int i = 0; i < numFields; i++) {
406 for (int j = 0; j < numPoints; j++) {
407 for (int k = 0; k < spaceDim; k++) {
408
409 // basisGrads is (F,P,D), compute offset:
410 int l = k + j * spaceDim + i * spaceDim * numPoints;
411 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
412 errorFlag++;
413 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
414
415 // Output the multi-index of the value where the error is:
416 *outStream << " At multi-index { ";
417 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
418 *outStream << "} computed grad component: " << vals(i,j,k)
419 << " but reference grad component: " << basisGrads[l] << "\n";
420 }
421 }
422 }
423 }
424
425 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
426 tetBasis.getValues(vals, tetNodes, OPERATOR_D1);
427 for (int i = 0; i < numFields; i++) {
428 for (int j = 0; j < numPoints; j++) {
429 for (int k = 0; k < spaceDim; k++) {
430
431 // basisGrads is (F,P,D), compute offset:
432 int l = k + j * spaceDim + i * spaceDim * numPoints;
433 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
434 errorFlag++;
435 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
436
437 // Output the multi-index of the value where the error is:
438 *outStream << " At multi-index { ";
439 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
440 *outStream << "} computed D1 component: " << vals(i,j,k)
441 << " but reference D1 component: " << basisGrads[l] << "\n";
442 }
443 }
444 }
445 }
446
447 // Check D2 of basis function
448 int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
449 vals.resize(numFields, numPoints, D2cardinality);
450 tetBasis.getValues(vals, tetNodes, OPERATOR_D2);
451 for (int i = 0; i < numFields; i++) {
452 for (int j = 0; j < numPoints; j++) {
453 for (int k = 0; k < D2cardinality; k++) {
454
455 // basisD2 is (F,P,Dk), compute offset:
456 int l = k + j * D2cardinality + i * D2cardinality * numPoints;
457 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
458 errorFlag++;
459 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
460
461 // Output the multi-index of the value where the error is:
462 *outStream << " At multi-index { ";
463 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
464 *outStream << "} computed D2 component: " << vals(i,j,k)
465 << " but reference D2 component: " << basisD2[l] << "\n";
466 }
467 }
468 }
469 }
470
471
472 // Check all higher derivatives - must be zero.
473 for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
474
475 // The last dimension is the number of kth derivatives and needs to be resized for every Dk
476 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
477 vals.resize(numFields, numPoints, DkCardin);
478
479 tetBasis.getValues(vals, tetNodes, op);
480 for (int i = 0; i < vals.size(); i++) {
481 if (std::abs(vals[i]) > INTREPID_TOL) {
482 errorFlag++;
483 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
484
485 // Get the multi-index of the value where the error is and the operator order
486 std::vector<int> myIndex;
487 vals.getMultiIndex(myIndex,i);
488 int ord = Intrepid::getOperatorOrder(op);
489 *outStream << " At multi-index { ";
490 for(int j = 0; j < vals.rank(); j++) {
491 *outStream << myIndex[j] << " ";
492 }
493 *outStream << "} computed D"<< ord <<" component: " << vals[i]
494 << " but reference D" << ord << " component: 0 \n";
495 }
496 }
497 }
498 }
499
500 // Catch unexpected errors
501 catch (const std::logic_error & err) {
502 *outStream << err.what() << "\n\n";
503 errorFlag = -1000;
504 };
505
506 if (errorFlag != 0)
507 std::cout << "End Result: TEST FAILED\n";
508 else
509 std::cout << "End Result: TEST PASSED\n";
510
511 // reset format state of std::cout
512 std::cout.copyfmt(oldFormatState);
513
514 return errorFlag;
515}
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_TET_C2_FEM class.
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.