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}
71
72int main(int argc, char *argv[]) {
73 Teuchos::oblackholestream oldFormatState;
74 oldFormatState.copyfmt(std::cout);
75
76 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
77
78 // This little trick lets us print to std::cout only if
79 // a (dummy) command-line argument is provided.
80 int iprint = argc - 1;
81 Teuchos::RCP<std::ostream> outStream;
82 Teuchos::oblackholestream bhs; // outputs nothing
83 if (iprint > 0)
84 outStream = Teuchos::rcp(&std::cout, false);
85 else
86 outStream = Teuchos::rcp(&bhs, false);
87
88 // Save the format state of the original std::cout.
89
90 *outStream \
91 << "===============================================================================\n" \
92 << "| |\n" \
93 << "| Unit Test (Basis_HGRAD_LINE_C1_FEM) |\n" \
94 << "| |\n" \
95 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
96 << "| 2) Basis values for VALUE, GRAD, CURL, DIV, and Dk operators |\n" \
97 << "| |\n" \
98 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
99 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
100 << "| Kara Peterson (dridzal@sandia.gov). |\n" \
101 << "| |\n" \
102 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
103 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
104 << "| |\n" \
105 << "===============================================================================\n"\
106 << "| TEST 1: Basis creation, exception testing |\n"\
107 << "===============================================================================\n";
108
109 // Define basis and error flag
110 Basis_HGRAD_LINE_C1_FEM<double, FieldContainer<double> > lineBasis;
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 2 vertices of the reference Line, its center and another point
118 FieldContainer<double> lineNodes(4, 1);
119 lineNodes(0,0) = -1.0;
120 lineNodes(1,0) = 1.0;
121 lineNodes(2,0) = 0.0;
122 lineNodes(3,0) = 0.5;
123
124 // Generic array for the output values; needs to be properly resized depending on the operator type
125 FieldContainer<double> vals;
126
127 try{
128 // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
129 vals.resize(lineBasis.getCardinality(), lineNodes.dimension(0) );
130
131 // Exceptions 1-5: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
132 // getDofTag() to access invalid array elements thereby causing bounds check exception
133 // exception #1
134 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(2,0,0), throwCounter, nException );
135 // exception #2
136 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,1,1), throwCounter, nException );
137 // exception #3
138 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(0,4,0), throwCounter, nException );
139 // exception #4
140 INTREPID_TEST_COMMAND( lineBasis.getDofTag(5), throwCounter, nException );
141 // exception #5
142 INTREPID_TEST_COMMAND( lineBasis.getDofTag(-1), throwCounter, nException );
143
144#ifdef HAVE_INTREPID_DEBUG
145 // Exceptions 6-17 test exception handling with incorrectly dimensioned input/output arrays
146 // exception #6: input points array must be of rank-2
147 FieldContainer<double> badPoints1(4, 5, 3);
148 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
149
150 // exception #7 dimension 1 in the input point array must equal space dimension of the cell
151 FieldContainer<double> badPoints2(4, 2);
152 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
153
154 // exception #8 output values must be of rank-2 for OPERATOR_VALUE
155 FieldContainer<double> badVals1(4, 3, 1);
156 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals1, lineNodes, OPERATOR_VALUE), throwCounter, nException );
157
158 // exception #9 output values must be of rank-3 for OPERATOR_GRAD
159 FieldContainer<double> badVals2(4, 3);
160 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_GRAD), throwCounter, nException );
161
162 // exception #10 output values must be of rank-3 for OPERATOR_DIV
163 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_DIV), throwCounter, nException );
164
165 // exception #11 output values must be of rank-3 for OPERATOR_CURL
166 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_CURL), throwCounter, nException );
167
168 // exception #12 output values must be of rank-3 for OPERATOR_D2
169 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_D2), throwCounter, nException );
170
171 // exception #13 incorrect 1st dimension of output array (must equal number of basis functions)
172 FieldContainer<double> badVals3(lineBasis.getCardinality() + 1, lineNodes.dimension(0));
173 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals3, lineNodes, OPERATOR_VALUE), throwCounter, nException );
174
175 // exception #14 incorrect 0th dimension of output array (must equal number of points)
176 FieldContainer<double> badVals4(lineBasis.getCardinality(), lineNodes.dimension(0) + 1);
177 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals4, lineNodes, OPERATOR_VALUE), throwCounter, nException );
178
179 // exception #15: incorrect 2nd dimension of output array (must equal the space dimension)
180 FieldContainer<double> badVals5(lineBasis.getCardinality(), lineNodes.dimension(0), 2);
181 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals5, lineNodes, OPERATOR_GRAD), throwCounter, nException );
182
183 // exception #16: incorrect 2nd dimension of output array (must equal D2 cardinality in 1D)
184 FieldContainer<double> badVals6(lineBasis.getCardinality(), lineNodes.dimension(0), 40);
185 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals6, lineNodes, OPERATOR_D2), throwCounter, nException );
186
187 // exception #17: incorrect 2nd dimension of output array (must equal D3 cardinality in 1D)
188 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals6, lineNodes, OPERATOR_D3), throwCounter, nException );
189#endif
190
191 }
192 catch (const std::logic_error & err) {
193 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
194 *outStream << err.what() << '\n';
195 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
196 errorFlag = -1000;
197 };
198
199 // Check if number of thrown exceptions matches the one we expect
200 if (throwCounter != nException) {
201 errorFlag++;
202 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
203 }
204
205 *outStream \
206 << "\n"
207 << "===============================================================================\n"\
208 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
209 << "===============================================================================\n";
210
211 try{
212 std::vector<std::vector<int> > allTags = lineBasis.getAllDofTags();
213
214 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
215 for (unsigned i = 0; i < allTags.size(); i++) {
216 int bfOrd = lineBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
217
218 std::vector<int> myTag = lineBasis.getDofTag(bfOrd);
219 if( !( (myTag[0] == allTags[i][0]) &&
220 (myTag[1] == allTags[i][1]) &&
221 (myTag[2] == allTags[i][2]) &&
222 (myTag[3] == allTags[i][3]) ) ) {
223 errorFlag++;
224 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
225 *outStream << " getDofOrdinal( {"
226 << allTags[i][0] << ", "
227 << allTags[i][1] << ", "
228 << allTags[i][2] << ", "
229 << allTags[i][3] << "}) = " << bfOrd <<" but \n";
230 *outStream << " getDofTag(" << bfOrd << ") = { "
231 << myTag[0] << ", "
232 << myTag[1] << ", "
233 << myTag[2] << ", "
234 << myTag[3] << "}\n";
235 }
236 }
237
238 // Now do the same but loop over basis functions
239 for( int bfOrd = 0; bfOrd < lineBasis.getCardinality(); bfOrd++) {
240 std::vector<int> myTag = lineBasis.getDofTag(bfOrd);
241 int myBfOrd = lineBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
242 if( bfOrd != myBfOrd) {
243 errorFlag++;
244 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
245 *outStream << " getDofTag(" << bfOrd << ") = { "
246 << myTag[0] << ", "
247 << myTag[1] << ", "
248 << myTag[2] << ", "
249 << myTag[3] << "} but getDofOrdinal({"
250 << myTag[0] << ", "
251 << myTag[1] << ", "
252 << myTag[2] << ", "
253 << myTag[3] << "} ) = " << myBfOrd << "\n";
254 }
255 }
256 }
257 catch (const std::logic_error & err){
258 *outStream << err.what() << "\n\n";
259 errorFlag = -1000;
260 };
261
262 *outStream \
263 << "\n"
264 << "===============================================================================\n"\
265 << "| TEST 3: correctness of basis function values |\n"\
266 << "===============================================================================\n";
267
268 outStream -> precision(20);
269
270 // VALUE: Each row gives the 2 correct basis set values at an evaluation point
271 double basisValues[] = {
272 1.0, 0.0,
273 0.0, 1.0,
274 0.5, 0.5,
275 0.25, 0.75
276 };
277
278 // GRAD, DIV, CURL and D1: each row gives the 2 correct values of the gradients of the 2 basis functions
279 double basisDerivs[] = {
280 -0.5, 0.5,
281 -0.5, 0.5,
282 -0.5, 0.5,
283 -0.5, 0.5
284 };
285
286 try{
287
288 // Dimensions for the output arrays:
289 int numFields = lineBasis.getCardinality();
290 int numPoints = lineNodes.dimension(0);
291 int spaceDim = lineBasis.getBaseCellTopology().getDimension();
292
293 // Generic array for values, grads, curls, etc. that will be properly sized before each call
294 FieldContainer<double> vals;
295
296 // Check VALUE of basis functions: resize vals to rank-2 container:
297 vals.resize(numFields, numPoints);
298 lineBasis.getValues(vals, lineNodes, OPERATOR_VALUE);
299 for (int i = 0; i < numFields; i++) {
300 for (int j = 0; j < numPoints; j++) {
301 int l = i + j * numFields;
302 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
303 errorFlag++;
304 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
305
306 // Output the multi-index of the value where the error is:
307 *outStream << " At multi-index { ";
308 *outStream << i << " ";*outStream << j << " ";
309 *outStream << "} computed value: " << vals(i,j)
310 << " but reference value: " << basisValues[l] << "\n";
311 }
312 }
313 }
314
315 // Check derivatives of basis function: resize vals to rank-3 container
316 vals.resize(numFields, numPoints, spaceDim);
317 lineBasis.getValues(vals, lineNodes, OPERATOR_GRAD);
318 for (int i = 0; i < numFields; i++) {
319 for (int j = 0; j < numPoints; j++) {
320 for (int k = 0; k < spaceDim; k++) {
321 int l = k + i * spaceDim + j * spaceDim * numFields;
322 if (std::abs(vals(i,j,k) - basisDerivs[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 << " ";*outStream << k << " ";
329 *outStream << "} computed grad component: " << vals(i,j,k)
330 << " but reference grad component: " << basisDerivs[l] << "\n";
331 }
332 }
333 }
334 }
335
336 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
337 lineBasis.getValues(vals, lineNodes, OPERATOR_D1);
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) - basisDerivs[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 D1 component: " << vals(i,j,k)
350 << " but reference D1 component: " << basisDerivs[l] << "\n";
351 }
352 }
353 }
354 }
355
356
357 // Check CURL of basis function: resize vals just for illustration!
358 vals.resize(numFields, numPoints, spaceDim);
359 lineBasis.getValues(vals, lineNodes, OPERATOR_CURL);
360 for (int i = 0; i < numFields; i++) {
361 for (int j = 0; j < numPoints; j++) {
362 for (int k = 0; k < spaceDim; k++) {
363 int l = k + i * spaceDim + j * spaceDim * numFields;
364 if (std::abs(vals(i,j,k) - basisDerivs[l]) > INTREPID_TOL) {
365 errorFlag++;
366 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
367
368 // Output the multi-index of the value where the error is:
369 *outStream << " At multi-index { ";
370 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
371 *outStream << "} computed curl component: " << vals(i,j,k)
372 << " but reference curl component: " << basisDerivs[l] << "\n";
373 }
374 }
375 }
376 }
377
378
379 // Check DIV of basis function (do not resize vals because it has the correct size: DIV = CURL)
380 lineBasis.getValues(vals, lineNodes, OPERATOR_DIV);
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) - basisDerivs[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 D1 component: " << vals(i,j,k)
393 << " but reference D1 component: " << basisDerivs[l] << "\n";
394 }
395 }
396 }
397 }
398
399
400 // Check all higher derivatives - must be zero.
401 for(EOperator op = OPERATOR_D2; op < OPERATOR_MAX; op++) {
402
403 // The last dimension is the number of kth derivatives and needs to be resized for every Dk
404 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
405 vals.resize(numFields, numPoints, DkCardin);
406
407 lineBasis.getValues(vals, lineNodes, op);
408 for (int i = 0; i < vals.size(); i++) {
409 if (std::abs(vals[i]) > INTREPID_TOL) {
410 errorFlag++;
411 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
412
413 // Get the multi-index of the value where the error is and the operator order
414 std::vector<int> myIndex;
415 vals.getMultiIndex(myIndex,i);
416 int ord = Intrepid::getOperatorOrder(op);
417 *outStream << " At multi-index { ";
418 for(int j = 0; j < vals.rank(); j++) {
419 *outStream << myIndex[j] << " ";
420 }
421 *outStream << "} computed D"<< ord <<" component: " << vals[i]
422 << " but reference D" << ord << " component: 0 \n";
423 }
424 }
425 }
426 }
427
428 // Catch unexpected errors
429 catch (const std::logic_error & err) {
430 *outStream << err.what() << "\n\n";
431 errorFlag = -1000;
432 };
433
434 if (errorFlag != 0)
435 std::cout << "End Result: TEST FAILED\n";
436 else
437 std::cout << "End Result: TEST PASSED\n";
438
439 // reset format state of std::cout
440 std::cout.copyfmt(oldFormatState);
441 return errorFlag;
442}
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_LINE_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.