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_C1_FEM.hpp"
50#include "Teuchos_oblackholestream.hpp"
51#include "Teuchos_RCP.hpp"
52#include "Teuchos_GlobalMPISession.hpp"
53
54using namespace std;
55using namespace Intrepid;
56
57#define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
58{ \
59 ++nException; \
60 try { \
61 S ; \
62 } \
63 catch (const std::logic_error & err) { \
64 ++throwCounter; \
65 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
66 *outStream << err.what() << '\n'; \
67 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
68 }; \
69}
70
71int main(int argc, char *argv[]) {
72
73 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
74
75 // This little trick lets us print to std::cout only if
76 // a (dummy) command-line argument is provided.
77 int iprint = argc - 1;
78 Teuchos::RCP<std::ostream> outStream;
79 Teuchos::oblackholestream bhs; // outputs nothing
80 if (iprint > 0)
81 outStream = Teuchos::rcp(&std::cout, false);
82 else
83 outStream = Teuchos::rcp(&bhs, false);
84
85 // Save the format state of the original std::cout.
86 Teuchos::oblackholestream oldFormatState;
87 oldFormatState.copyfmt(std::cout);
88
89 *outStream \
90 << "===============================================================================\n" \
91 << "| |\n" \
92 << "| Unit Test (Basis_HGRAD_QUAD_C1_FEM) |\n" \
93 << "| |\n" \
94 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 << "| 2) Basis values for VALUE, GRAD, CURL, and Dk operators |\n" \
96 << "| |\n" \
97 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
98 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
99 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
100 << "| |\n" \
101 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
102 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
103 << "| |\n" \
104 << "===============================================================================\n"\
105 << "| TEST 1: Basis creation, exception testing |\n"\
106 << "===============================================================================\n";
107
108 // Define basis and error flag
110 int errorFlag = 0;
111
112 // Initialize throw counter for exception testing
113 int nException = 0;
114 int throwCounter = 0;
115
116 // Define array containing the 4 vertices of the reference QUAD and its center.
117 FieldContainer<double> quadNodes(5, 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 quadNodes(4,0) = 0.0; quadNodes(4,1) = 0.0;
123
124 // Generic array for the output values; needs to be properly resized depending on the operator type
126
127 try{
128 // exception #1: DIV cannot be applied to scalar functions
129 // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
130 vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0));
131 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException );
132
133 // Exceptions 2-6: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
134 // getDofTag() to access invalid array elements thereby causing bounds check exception
135 // exception #2
136 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException );
137 // exception #3
138 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException );
139 // exception #4
140 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,0), throwCounter, nException );
141 // exception #5
142 INTREPID_TEST_COMMAND( quadBasis.getDofTag(5), throwCounter, nException );
143 // exception #6
144 INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException );
145
146#ifdef HAVE_INTREPID_DEBUG
147 // Exceptions 7- test exception handling with incorrectly dimensioned input/output arrays
148 // exception #7: input points array must be of rank-2
149 FieldContainer<double> badPoints1(4, 5, 3);
150 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
151
152 // exception #8 dimension 1 in the input point array must equal space dimension of the cell
153 FieldContainer<double> badPoints2(4, 3);
154 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
155
156 // exception #9 output values must be of rank-2 for OPERATOR_VALUE
157 FieldContainer<double> badVals1(4, 3, 1);
158 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
159
160 // exception #10 output values must be of rank-3 for OPERATOR_GRAD
161 FieldContainer<double> badVals2(4, 3);
162 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_GRAD), throwCounter, nException );
163
164 // exception #11 output values must be of rank-3 for OPERATOR_CURL
165 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_CURL), throwCounter, nException );
166
167 // exception #12 output values must be of rank-3 for OPERATOR_D2
168 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_D2), throwCounter, nException );
169
170 // exception #13 incorrect 0th dimension of output array (must equal number of basis functions)
171 FieldContainer<double> badVals3(quadBasis.getCardinality() + 1, quadNodes.dimension(0));
172 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException );
173
174 // exception #14 incorrect 1st dimension of output array (must equal number of points)
175 FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0) + 1);
176 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException );
177
178 // exception #15: incorrect 2nd dimension of output array (must equal the space dimension)
179 FieldContainer<double> badVals5(quadBasis.getCardinality(), quadNodes.dimension(0), 4);
180 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals5, quadNodes, OPERATOR_GRAD), throwCounter, nException );
181
182 // exception #16: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D)
183 FieldContainer<double> badVals6(quadBasis.getCardinality(), quadNodes.dimension(0), 40);
184 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D2), throwCounter, nException );
185
186 // exception #17: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D)
187 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D3), throwCounter, nException );
188#endif
189
190 }
191 catch (const std::logic_error & err) {
192 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
193 *outStream << err.what() << '\n';
194 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
195 errorFlag = -1000;
196 };
197
198 // Check if number of thrown exceptions matches the one we expect
199 if (throwCounter != nException) {
200 errorFlag++;
201 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
202 }
203
204 *outStream \
205 << "\n"
206 << "===============================================================================\n"\
207 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
208 << "===============================================================================\n";
209
210 try{
211 std::vector<std::vector<int> > allTags = quadBasis.getAllDofTags();
212
213 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
214 for (unsigned i = 0; i < allTags.size(); i++) {
215 int bfOrd = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
216
217 std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
218 if( !( (myTag[0] == allTags[i][0]) &&
219 (myTag[1] == allTags[i][1]) &&
220 (myTag[2] == allTags[i][2]) &&
221 (myTag[3] == allTags[i][3]) ) ) {
222 errorFlag++;
223 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
224 *outStream << " getDofOrdinal( {"
225 << allTags[i][0] << ", "
226 << allTags[i][1] << ", "
227 << allTags[i][2] << ", "
228 << allTags[i][3] << "}) = " << bfOrd <<" but \n";
229 *outStream << " getDofTag(" << bfOrd << ") = { "
230 << myTag[0] << ", "
231 << myTag[1] << ", "
232 << myTag[2] << ", "
233 << myTag[3] << "}\n";
234 }
235 }
236
237 // Now do the same but loop over basis functions
238 for( int bfOrd = 0; bfOrd < quadBasis.getCardinality(); bfOrd++) {
239 std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
240 int myBfOrd = quadBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
241 if( bfOrd != myBfOrd) {
242 errorFlag++;
243 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
244 *outStream << " getDofTag(" << bfOrd << ") = { "
245 << myTag[0] << ", "
246 << myTag[1] << ", "
247 << myTag[2] << ", "
248 << myTag[3] << "} but getDofOrdinal({"
249 << myTag[0] << ", "
250 << myTag[1] << ", "
251 << myTag[2] << ", "
252 << myTag[3] << "} ) = " << myBfOrd << "\n";
253 }
254 }
255 }
256 catch (const std::logic_error & err){
257 *outStream << err.what() << "\n\n";
258 errorFlag = -1000;
259 };
260
261 *outStream \
262 << "\n"
263 << "===============================================================================\n"\
264 << "| TEST 3: correctness of basis function values |\n"\
265 << "===============================================================================\n";
266
267 outStream -> precision(20);
268
269 // VALUE: Each row gives the 4 correct basis set values at an evaluation point
270 double basisValues[] = {
271 1.0, 0.0, 0.0, 0.0,
272 0.0, 1.0, 0.0, 0.0,
273 0.0, 0.0, 1.0, 0.0,
274 0.0, 0.0, 0.0, 1.0,
275 0.25,0.25,0.25,0.25
276 };
277
278 // GRAD and D1: each row gives the 8 correct values of the gradients of the 4 basis functions
279 double basisGrads[] = {
280 -0.5, -0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5,
281 -0.5, 0.0, 0.5, -0.5, 0.0, 0.5, 0.0, 0.0,
282 0.0, 0.0, 0.0, -0.5, 0.5, 0.5, -0.5, 0.0,
283 0.0, -0.5, 0.0, 0.0, 0.5, 0.0, -0.5, 0.5,
284 -0.25,-0.25, 0.25,-0.25, 0.25, 0.25, -0.25, 0.25
285 };
286
287 // CURL: each row gives the 8 correct values of the curls of the 4 basis functions
288 double basisCurls[] = {
289 -0.5, 0.5, 0.0, -0.5, 0.0, 0.0, 0.5, 0.0,
290 0.0, 0.5, -0.5, -0.5, 0.5, 0.0, 0.0, 0.0,
291 0.0, 0.0, -0.5, 0.0, 0.5, -0.5, 0.0, 0.5,
292 -0.5, 0.0, 0.0, 0.0, 0.0, -0.5, 0.5, 0.5,
293 -0.25, 0.25, -0.25,-0.25, 0.25,-0.25, 0.25, 0.25
294 };
295
296 //D2: each row gives the 12 correct values of all 2nd derivatives of the 4 basis functions
297 double basisD2[] = {
298 0.0, 0.25, 0.0, 0.0,-0.25, 0.0, 0.0, 0.25, 0.0, 0.0,-0.25, 0.0,
299 0.0, 0.25, 0.0, 0.0,-0.25, 0.0, 0.0, 0.25, 0.0, 0.0,-0.25, 0.0,
300 0.0, 0.25, 0.0, 0.0,-0.25, 0.0, 0.0, 0.25, 0.0, 0.0,-0.25, 0.0,
301 0.0, 0.25, 0.0, 0.0,-0.25, 0.0, 0.0, 0.25, 0.0, 0.0,-0.25, 0.0,
302 0.0, 0.25, 0.0, 0.0,-0.25, 0.0, 0.0, 0.25, 0.0, 0.0,-0.25, 0.0,
303 };
304
305 try{
306
307 // Dimensions for the output arrays:
308 int numFields = quadBasis.getCardinality();
309 int numPoints = quadNodes.dimension(0);
310 int spaceDim = quadBasis.getBaseCellTopology().getDimension();
311 int D2Cardin = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
312
313 // Generic array for values, grads, curls, etc. that will be properly sized before each call
315
316 // Check VALUE of basis functions: resize vals to rank-2 container:
317 vals.resize(numFields, numPoints);
318 quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE);
319 for (int i = 0; i < numFields; i++) {
320 for (int j = 0; j < numPoints; j++) {
321 int l = i + j * numFields;
322 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
323 errorFlag++;
324 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
325
326 // Output the multi-index of the value where the error is:
327 *outStream << " At multi-index { ";
328 *outStream << i << " ";*outStream << j << " ";
329 *outStream << "} computed value: " << vals(i,j)
330 << " but reference value: " << basisValues[l] << "\n";
331 }
332 }
333 }
334
335 // Check GRAD of basis function: resize vals to rank-3 container
336 vals.resize(numFields, numPoints, spaceDim);
337 quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD);
338 for (int i = 0; i < numFields; i++) {
339 for (int j = 0; j < numPoints; j++) {
340 for (int k = 0; k < spaceDim; k++) {
341 int l = k + i * spaceDim + j * spaceDim * numFields;
342 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
343 errorFlag++;
344 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
345
346 // Output the multi-index of the value where the error is:
347 *outStream << " At multi-index { ";
348 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
349 *outStream << "} computed grad component: " << vals(i,j,k)
350 << " but reference grad component: " << basisGrads[l] << "\n";
351 }
352 }
353 }
354 }
355
356
357 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
358 quadBasis.getValues(vals, quadNodes, OPERATOR_D1);
359 for (int i = 0; i < numFields; i++) {
360 for (int j = 0; j < numPoints; j++) {
361 for (int k = 0; k < spaceDim; k++) {
362 int l = k + i * spaceDim + j * spaceDim * numFields;
363 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
364 errorFlag++;
365 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
366
367 // Output the multi-index of the value where the error is:
368 *outStream << " At multi-index { ";
369 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
370 *outStream << "} computed D1 component: " << vals(i,j,k)
371 << " but reference D1 component: " << basisGrads[l] << "\n";
372 }
373 }
374 }
375 }
376
377
378 // Check CURL of basis function: resize vals just for illustration!
379 vals.resize(numFields, numPoints, spaceDim);
380 quadBasis.getValues(vals, quadNodes, OPERATOR_CURL);
381 for (int i = 0; i < numFields; i++) {
382 for (int j = 0; j < numPoints; j++) {
383 for (int k = 0; k < spaceDim; k++) {
384 int l = k + i * spaceDim + j * spaceDim * numFields;
385 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
386 errorFlag++;
387 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
388
389 // Output the multi-index of the value where the error is:
390 *outStream << " At multi-index { ";
391 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
392 *outStream << "} computed curl component: " << vals(i,j,k)
393 << " but reference curl component: " << basisCurls[l] << "\n";
394 }
395 }
396 }
397 }
398
399 // Check D2 of basis function
400 vals.resize(numFields, numPoints, D2Cardin);
401 quadBasis.getValues(vals, quadNodes, OPERATOR_D2);
402 for (int i = 0; i < numFields; i++) {
403 for (int j = 0; j < numPoints; j++) {
404 for (int k = 0; k < D2Cardin; k++) {
405 int l = k + i * D2Cardin + j * D2Cardin * numFields;
406 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
407 errorFlag++;
408 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
409
410 // Output the multi-index of the value where the error is:
411 *outStream << " At multi-index { ";
412 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
413 *outStream << "} computed D2 component: " << vals(i,j,k)
414 << " but reference D2 component: " << basisD2[l] << "\n";
415 }
416 }
417 }
418 }
419
420
421 // Check all higher derivatives - must be zero.
422 for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
423
424 // The last dimension is the number of kth derivatives and needs to be resized for every Dk
425 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
426 vals.resize(numFields, numPoints, DkCardin);
427
428 quadBasis.getValues(vals, quadNodes, op);
429 for (int i = 0; i < vals.size(); i++) {
430 if (std::abs(vals[i]) > INTREPID_TOL) {
431 errorFlag++;
432 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
433
434 // Get the multi-index of the value where the error is and the operator order
435 std::vector<int> myIndex;
436 vals.getMultiIndex(myIndex,i);
437 int ord = Intrepid::getOperatorOrder(op);
438 *outStream << " At multi-index { ";
439 for(int j = 0; j < vals.rank(); j++) {
440 *outStream << myIndex[j] << " ";
441 }
442 *outStream << "} computed D"<< ord <<" component: " << vals[i]
443 << " but reference D" << ord << " component: 0 \n";
444 }
445 }
446 }
447 }
448 // Catch unexpected errors
449 catch (const std::logic_error & err) {
450 *outStream << err.what() << "\n\n";
451 errorFlag = -1000;
452 };
453
454
455 *outStream \
456 << "\n"
457 << "===============================================================================\n"\
458 << "| TEST 4: correctness of DoF locations |\n"\
459 << "===============================================================================\n";
460
461 try{
462 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
463 Teuchos::rcp(new Basis_HGRAD_QUAD_C1_FEM<double, FieldContainer<double> >);
464 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
465 Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
466
468 FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality());
469
470 // Check exceptions.
471#ifdef HAVE_INTREPID_DEBUG
472 cvals.resize(1,2,3);
473 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
474 cvals.resize(3,2);
475 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
476 cvals.resize(4,3);
477 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
478#endif
479 cvals.resize(4,2);
480 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
481 // Check if number of thrown exceptions matches the one we expect
482 if (throwCounter != nException) {
483 errorFlag++;
484 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
485 }
486
487 // Check mathematical correctness.
488 basis->getValues(bvals, cvals, OPERATOR_VALUE);
489 char buffer[120];
490 for (int i=0; i<bvals.dimension(0); i++) {
491 for (int j=0; j<bvals.dimension(1); j++) {
492 if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) {
493 errorFlag++;
494 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);
495 *outStream << buffer;
496 }
497 else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) {
498 errorFlag++;
499 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);
500 *outStream << buffer;
501 }
502 }
503 }
504
505 }
506 catch (const std::logic_error & err){
507 *outStream << err.what() << "\n\n";
508 errorFlag = -1000;
509 };
510
511 if (errorFlag != 0)
512 std::cout << "End Result: TEST FAILED\n";
513 else
514 std::cout << "End Result: TEST PASSED\n";
515
516 // reset format state of std::cout
517 std::cout.copyfmt(oldFormatState);
518
519 return errorFlag;
520}
int main(int argc, char *argv[])
Performs a code-code comparison to FIAT for Nedelec bases on tets (values and curls)
Definition: test_01.cpp:65
Header file for utility class to provide multidimensional containers.
int getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Quadrilateral cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Quadrilateral cell.
virtual int getDofOrdinal(const int subcDim, const int subcOrd, const int subcDofOrd)
DoF tag to ordinal lookup.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
virtual int getCardinality() const
Returns cardinality of the basis.
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.