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