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
51#include "Teuchos_oblackholestream.hpp"
52#include "Teuchos_RCP.hpp"
53#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}
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_TRI_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 (dridzal@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
109 // Define basis and error flag
111 int errorFlag = 0;
112
113 // Initialize throw counter for exception testing
114 int nException = 0;
115 int throwCounter = 0;
116
117 // Define array containing the 3 vertices of the reference Triangle, its center and another point
118 FieldContainer<double> triNodes(5, 2);
119 triNodes(0,0) = 0.0; triNodes(0,1) = 0.0;
120 triNodes(1,0) = 1.0; triNodes(1,1) = 0.0;
121 triNodes(2,0) = 0.0; triNodes(2,1) = 1.0;
122 triNodes(3,0) = 0.5; triNodes(3,1) = 0.5;
123 triNodes(4,0) = 0.0; triNodes(4,1) = 0.75;
124
125
126 try{
127 // Generic array for the output values; needs to be properly resized depending on the operator type
129
130 // exception #1: DIV cannot be applied to scalar functions
131 // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
132 vals.resize(triBasis.getCardinality(), triNodes.dimension(0) );
133 INTREPID_TEST_COMMAND( triBasis.getValues(vals, triNodes, OPERATOR_DIV), throwCounter, nException );
134
135 // Exceptions 2-6: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
136 // getDofTag() to access invalid array elements thereby causing bounds check exception
137 // exception #2
138 INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(2,0,0), throwCounter, nException );
139 // exception #3
140 INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(1,1,1), throwCounter, nException );
141 // exception #4
142 INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(0,4,0), throwCounter, nException );
143 // exception #5
144 INTREPID_TEST_COMMAND( triBasis.getDofTag(5), throwCounter, nException );
145 // exception #6
146 INTREPID_TEST_COMMAND( triBasis.getDofTag(-1), throwCounter, nException );
147
148#ifdef HAVE_INTREPID_DEBUG
149 // Exceptions 7-17 test exception handling with incorrectly dimensioned input/output arrays
150 // exception #7: input points array must be of rank-2
151 FieldContainer<double> badPoints1(4, 5, 3);
152 INTREPID_TEST_COMMAND( triBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
153
154 // exception #8 dimension 1 in the input point array must equal space dimension of the cell
155 FieldContainer<double> badPoints2(4, 3);
156 INTREPID_TEST_COMMAND( triBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
157
158 // exception #9 output values must be of rank-2 for OPERATOR_VALUE
159 FieldContainer<double> badVals1(4, 3, 1);
160 INTREPID_TEST_COMMAND( triBasis.getValues(badVals1, triNodes, OPERATOR_VALUE), throwCounter, nException );
161
162 // exception #10 output values must be of rank-3 for OPERATOR_GRAD
163 FieldContainer<double> badVals2(4, 3);
164 INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_GRAD), throwCounter, nException );
165
166 // exception #11 output values must be of rank-3 for OPERATOR_CURL
167 INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_CURL), throwCounter, nException );
168
169 // exception #12 output values must be of rank-3 for OPERATOR_D2
170 INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_D2), throwCounter, nException );
171
172 // exception #13 incorrect 1st dimension of output array (must equal number of basis functions)
173 FieldContainer<double> badVals3(triBasis.getCardinality() + 1, triNodes.dimension(0));
174 INTREPID_TEST_COMMAND( triBasis.getValues(badVals3, triNodes, OPERATOR_VALUE), throwCounter, nException );
175
176 // exception #14 incorrect 0th dimension of output array (must equal number of points)
177 FieldContainer<double> badVals4(triBasis.getCardinality(), triNodes.dimension(0) + 1);
178 INTREPID_TEST_COMMAND( triBasis.getValues(badVals4, triNodes, OPERATOR_VALUE), throwCounter, nException );
179
180 // exception #15: incorrect 2nd dimension of output array (must equal the space dimension)
181 FieldContainer<double> badVals5(triBasis.getCardinality(), triNodes.dimension(0), 4);
182 INTREPID_TEST_COMMAND( triBasis.getValues(badVals5, triNodes, OPERATOR_GRAD), throwCounter, nException );
183
184 // exception #16: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D)
185 FieldContainer<double> badVals6(triBasis.getCardinality(), triNodes.dimension(0), 40);
186 INTREPID_TEST_COMMAND( triBasis.getValues(badVals6, triNodes, OPERATOR_D2), throwCounter, nException );
187
188 // exception #17: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D)
189 INTREPID_TEST_COMMAND( triBasis.getValues(badVals6, triNodes, OPERATOR_D3), throwCounter, nException );
190#endif
191
192 }
193 catch (const std::logic_error & err) {
194 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
195 *outStream << err.what() << '\n';
196 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
197 errorFlag = -1000;
198 };
199
200 // Check if number of thrown exceptions matches the one we expect
201 if (throwCounter != nException) {
202 errorFlag++;
203 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
204 }
205
206 *outStream \
207 << "\n"
208 << "===============================================================================\n"\
209 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
210 << "===============================================================================\n";
211
212 try{
213 std::vector<std::vector<int> > allTags = triBasis.getAllDofTags();
214
215 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
216 for (unsigned i = 0; i < allTags.size(); i++) {
217 int bfOrd = triBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
218
219 std::vector<int> myTag = triBasis.getDofTag(bfOrd);
220 if( !( (myTag[0] == allTags[i][0]) &&
221 (myTag[1] == allTags[i][1]) &&
222 (myTag[2] == allTags[i][2]) &&
223 (myTag[3] == allTags[i][3]) ) ) {
224 errorFlag++;
225 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
226 *outStream << " getDofOrdinal( {"
227 << allTags[i][0] << ", "
228 << allTags[i][1] << ", "
229 << allTags[i][2] << ", "
230 << allTags[i][3] << "}) = " << bfOrd <<" but \n";
231 *outStream << " getDofTag(" << bfOrd << ") = { "
232 << myTag[0] << ", "
233 << myTag[1] << ", "
234 << myTag[2] << ", "
235 << myTag[3] << "}\n";
236 }
237 }
238
239 // Now do the same but loop over basis functions
240 for( int bfOrd = 0; bfOrd < triBasis.getCardinality(); bfOrd++) {
241 std::vector<int> myTag = triBasis.getDofTag(bfOrd);
242 int myBfOrd = triBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
243 if( bfOrd != myBfOrd) {
244 errorFlag++;
245 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
246 *outStream << " getDofTag(" << bfOrd << ") = { "
247 << myTag[0] << ", "
248 << myTag[1] << ", "
249 << myTag[2] << ", "
250 << myTag[3] << "} but getDofOrdinal({"
251 << myTag[0] << ", "
252 << myTag[1] << ", "
253 << myTag[2] << ", "
254 << myTag[3] << "} ) = " << myBfOrd << "\n";
255 }
256 }
257 }
258 catch (const std::logic_error & err){
259 *outStream << err.what() << "\n\n";
260 errorFlag = -1000;
261 };
262
263 *outStream \
264 << "\n"
265 << "===============================================================================\n"\
266 << "| TEST 3: correctness of basis function values |\n"\
267 << "===============================================================================\n";
268
269 outStream -> precision(20);
270
271 // VALUE: Each row gives the 3 correct basis set values at an evaluation point
272 double basisValues[] = {
273 1.0, 0.0, 0.0,
274 0.0, 1.0, 0.0,
275 0.0, 0.0, 1.0,
276 0.0, 0.5, 0.5,
277 0.25,0.0, 0.75
278 };
279
280 // GRAD and D1: each row gives the 6 correct values of the gradients of the 3 basis functions
281 double basisGrads[] = {
282 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
283 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
284 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
285 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
286 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
287 };
288
289 // CURL: each row gives the 6 correct values of the curls of the 3 basis functions
290 double basisCurls[] = {
291 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0,
292 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0,
293 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0,
294 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0,
295 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0
296 };
297
298 try{
299
300 // Dimensions for the output arrays:
301 int numFields = triBasis.getCardinality();
302 int numPoints = triNodes.dimension(0);
303 int spaceDim = triBasis.getBaseCellTopology().getDimension();
304
305 // Generic array for values, grads, curls, etc. that will be properly sized before each call
307
308 // Check VALUE of basis functions: resize vals to rank-2 container:
309 vals.resize(numFields, numPoints);
310 triBasis.getValues(vals, triNodes, OPERATOR_VALUE);
311 for (int i = 0; i < numFields; i++) {
312 for (int j = 0; j < numPoints; j++) {
313 int l = i + j * numFields;
314 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
315 errorFlag++;
316 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
317
318 // Output the multi-index of the value where the error is:
319 *outStream << " At multi-index { ";
320 *outStream << i << " ";*outStream << j << " ";
321 *outStream << "} computed value: " << vals(i,j)
322 << " but reference value: " << basisValues[l] << "\n";
323 }
324 }
325 }
326
327 // Check GRAD of basis function: resize vals to rank-3 container
328 vals.resize(numFields, numPoints, spaceDim);
329 triBasis.getValues(vals, triNodes, OPERATOR_GRAD);
330 for (int i = 0; i < numFields; i++) {
331 for (int j = 0; j < numPoints; j++) {
332 for (int k = 0; k < spaceDim; k++) {
333 int l = k + i * spaceDim + j * spaceDim * numFields;
334 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
335 errorFlag++;
336 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
337
338 // Output the multi-index of the value where the error is:
339 *outStream << " At multi-index { ";
340 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
341 *outStream << "} computed grad component: " << vals(i,j,k)
342 << " but reference grad component: " << basisGrads[l] << "\n";
343 }
344 }
345 }
346 }
347
348 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
349 triBasis.getValues(vals, triNodes, OPERATOR_D1);
350 for (int i = 0; i < numFields; i++) {
351 for (int j = 0; j < numPoints; j++) {
352 for (int k = 0; k < spaceDim; k++) {
353 int l = k + i * spaceDim + j * spaceDim * numFields;
354 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
355 errorFlag++;
356 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
357
358 // Output the multi-index of the value where the error is:
359 *outStream << " At multi-index { ";
360 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
361 *outStream << "} computed D1 component: " << vals(i,j,k)
362 << " but reference D1 component: " << basisGrads[l] << "\n";
363 }
364 }
365 }
366 }
367
368
369 // Check CURL of basis function: resize vals just for illustration!
370 vals.resize(numFields, numPoints, spaceDim);
371 triBasis.getValues(vals, triNodes, OPERATOR_CURL);
372 for (int i = 0; i < numFields; i++) {
373 for (int j = 0; j < numPoints; j++) {
374 for (int k = 0; k < spaceDim; k++) {
375 int l = k + i * spaceDim + j * spaceDim * numFields;
376 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
377 errorFlag++;
378 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
379
380 // Output the multi-index of the value where the error is:
381 *outStream << " At multi-index { ";
382 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
383 *outStream << "} computed curl component: " << vals(i,j,k)
384 << " but reference curl component: " << basisCurls[l] << "\n";
385 }
386 }
387 }
388 }
389
390 // Check all higher derivatives - must be zero.
391 for(EOperator op = OPERATOR_D2; op < OPERATOR_MAX; op++) {
392
393 // The last dimension is the number of kth derivatives and needs to be resized for every Dk
394 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
395 vals.resize(numFields, numPoints, DkCardin);
396
397 triBasis.getValues(vals, triNodes, op);
398 for (int i = 0; i < vals.size(); i++) {
399 if (std::abs(vals[i]) > INTREPID_TOL) {
400 errorFlag++;
401 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
402
403 // Get the multi-index of the value where the error is and the operator order
404 std::vector<int> myIndex;
405 vals.getMultiIndex(myIndex,i);
406 int ord = Intrepid::getOperatorOrder(op);
407 *outStream << " At multi-index { ";
408 for(int j = 0; j < vals.rank(); j++) {
409 *outStream << myIndex[j] << " ";
410 }
411 *outStream << "} computed D"<< ord <<" component: " << vals[i]
412 << " but reference D" << ord << " component: 0 \n";
413 }
414 }
415 }
416 }
417
418 // Catch unexpected errors
419 catch (const std::logic_error & err) {
420 *outStream << err.what() << "\n\n";
421 errorFlag = -1000;
422 };
423
424 *outStream \
425 << "\n"
426 << "===============================================================================\n"\
427 << "| TEST 4: correctness of DoF locations |\n"\
428 << "===============================================================================\n";
429
430 try{
431 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
432 Teuchos::rcp(new Basis_HGRAD_TRI_C1_FEM<double, FieldContainer<double> >);
433 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
434 Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
435
437 FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality());
438
439 // Check exceptions.
440#ifdef HAVE_INTREPID_DEBUG
441 cvals.resize(1,2,3);
442 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
443 cvals.resize(4,2);
444 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
445 cvals.resize(4,3);
446 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
447#endif
448 cvals.resize(3,2);
449 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
450 // Check if number of thrown exceptions matches the one we expect
451 if (throwCounter != nException) {
452 errorFlag++;
453 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
454 }
455
456 // Check mathematical correctness.
457 basis->getValues(bvals, cvals, OPERATOR_VALUE);
458 char buffer[120];
459 for (int i=0; i<bvals.dimension(0); i++) {
460 for (int j=0; j<bvals.dimension(1); j++) {
461 if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) {
462 errorFlag++;
463 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);
464 *outStream << buffer;
465 }
466 else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) {
467 errorFlag++;
468 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);
469 *outStream << buffer;
470 }
471 }
472 }
473
474 }
475 catch (const std::logic_error & err){
476 *outStream << err.what() << "\n\n";
477 errorFlag = -1000;
478 };
479
480 if (errorFlag != 0)
481 std::cout << "End Result: TEST FAILED\n";
482 else
483 std::cout << "End Result: TEST PASSED\n";
484
485 // reset format state of std::cout
486 std::cout.copyfmt(oldFormatState);
487
488 return errorFlag;
489}
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_TRI_C1_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.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Triangle cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
virtual int getDofOrdinal(const int subcDim, const int subcOrd, const int subcDofOrd)
DoF tag to ordinal lookup.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
virtual int getCardinality() const
Returns cardinality of the basis.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
void getMultiIndex(int &i0, const int valueEnum) const
Returns the multi-index of a value, based on its enumeration, as a list, for rank-1 containers.
int rank() const
Return rank of the FieldContainer = number of indices used to tag the multi-indexed value.