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