Intrepid
test_01.cpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
50#include "Teuchos_oblackholestream.hpp"
51#include "Teuchos_RCP.hpp"
52#include "Teuchos_GlobalMPISession.hpp"
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_HCURL_QUAD_I1_FEM) |\n" \
93 << "| |\n" \
94 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 << "| 2) Basis values for VALUE and CURL 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_HCURL_QUAD_I1_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 of the reference Quadrilateral, its center and 4 more points
117 FieldContainer<double> quadNodes(9, 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
123 quadNodes(4,0) = 0.0; quadNodes(4,1) = 0.0;
124 quadNodes(5,0) = 0.0; quadNodes(5,1) = -0.5;
125 quadNodes(6,0) = 0.0; quadNodes(6,1) = 0.5;
126 quadNodes(7,0) = -0.5; quadNodes(7,1) = 0.0;
127 quadNodes(8,0) = 0.5; quadNodes(8,1) = 0.0;
128
129
130 try{
131 // Generic array for the output values; needs to be properly resized depending on the operator type
132 FieldContainer<double> vals;
133
134 // exception #1: GRAD cannot be applied to HCURL functions
135 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
136 vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0), 4 );
137 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD), throwCounter, nException );
138
139 // exception #2: DIV cannot be applied to HCURL functions
140 // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
141 vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0) );
142 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException );
143
144 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
145 // getDofTag() to access invalid array elements thereby causing bounds check exception
146 // exception #3
147 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException );
148 // exception #4
149 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException );
150 // exception #5
151 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,1), throwCounter, nException );
152 // exception #6
153 INTREPID_TEST_COMMAND( quadBasis.getDofTag(12), throwCounter, nException );
154 // exception #7
155 INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException );
156
157#ifdef HAVE_INTREPID_DEBUG
158 // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays
159 // exception #8: input points array must be of rank-2
160 FieldContainer<double> badPoints1(4, 5, 3);
161 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
162
163 // exception #9 dimension 1 in the input point array must equal space dimension of the cell
164 FieldContainer<double> badPoints2(4, quadBasis.getBaseCellTopology().getDimension() + 1);
165 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
166
167 // exception #10 output values must be of rank-3 for OPERATOR_VALUE in 2D
168 FieldContainer<double> badVals1(4, 3);
169 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
170
171 FieldContainer<double> badCurls1(4,3,2);
172 // exception #11 output values must be of rank-2 for OPERATOR_CURL
173 INTREPID_TEST_COMMAND( quadBasis.getValues(badCurls1, quadNodes, OPERATOR_CURL), throwCounter, nException );
174
175 // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
176 FieldContainer<double> badVals2(quadBasis.getCardinality() + 1, quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension());
177 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
178
179 // exception #13 incorrect 1st dimension of output array (must equal number of points)
180 FieldContainer<double> badVals3(quadBasis.getCardinality(), quadNodes.dimension(0) + 1, quadBasis.getBaseCellTopology().getDimension() );
181 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
182
183 // exception #14: incorrect 2nd dimension of output array for VALUE (must equal the space dimension)
184 FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() - 1);
185 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
186
187 // exception #15: D2 cannot be applied to HCURL functions
188 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
189 vals.resize(quadBasis.getCardinality(),
190 quadNodes.dimension(0),
191 Intrepid::getDkCardinality(OPERATOR_D2, quadBasis.getBaseCellTopology().getDimension()));
192 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_D2), throwCounter, nException );
193#endif
194
195 }
196 catch (const std::logic_error & err) {
197 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
198 *outStream << err.what() << '\n';
199 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
200 errorFlag = -1000;
201 };
202
203 // Check if number of thrown exceptions matches the one we expect
204 // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match.
205 if (throwCounter != nException) {
206 errorFlag++;
207 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
208 }
209//#endif
210
211 *outStream \
212 << "\n"
213 << "===============================================================================\n"\
214 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
215 << "===============================================================================\n";
216
217 try{
218 std::vector<std::vector<int> > allTags = quadBasis.getAllDofTags();
219
220 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
221 for (unsigned i = 0; i < allTags.size(); i++) {
222 int bfOrd = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
223
224 std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
225 if( !( (myTag[0] == allTags[i][0]) &&
226 (myTag[1] == allTags[i][1]) &&
227 (myTag[2] == allTags[i][2]) &&
228 (myTag[3] == allTags[i][3]) ) ) {
229 errorFlag++;
230 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
231 *outStream << " getDofOrdinal( {"
232 << allTags[i][0] << ", "
233 << allTags[i][1] << ", "
234 << allTags[i][2] << ", "
235 << allTags[i][3] << "}) = " << bfOrd <<" but \n";
236 *outStream << " getDofTag(" << bfOrd << ") = { "
237 << myTag[0] << ", "
238 << myTag[1] << ", "
239 << myTag[2] << ", "
240 << myTag[3] << "}\n";
241 }
242 }
243
244 // Now do the same but loop over basis functions
245 for( int bfOrd = 0; bfOrd < quadBasis.getCardinality(); bfOrd++) {
246 std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
247 int myBfOrd = quadBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
248 if( bfOrd != myBfOrd) {
249 errorFlag++;
250 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
251 *outStream << " getDofTag(" << bfOrd << ") = { "
252 << myTag[0] << ", "
253 << myTag[1] << ", "
254 << myTag[2] << ", "
255 << myTag[3] << "} but getDofOrdinal({"
256 << myTag[0] << ", "
257 << myTag[1] << ", "
258 << myTag[2] << ", "
259 << myTag[3] << "} ) = " << myBfOrd << "\n";
260 }
261 }
262 }
263 catch (const std::logic_error & err){
264 *outStream << err.what() << "\n\n";
265 errorFlag = -1000;
266 };
267
268 *outStream \
269 << "\n"
270 << "===============================================================================\n"\
271 << "| TEST 3: correctness of basis function values |\n"\
272 << "===============================================================================\n";
273
274 outStream -> precision(20);
275
276 // VALUE: Each row pair gives the 4x2 correct basis set values at an evaluation point: (P,F,D) layout
277 double basisValues[] = {
278 0.500000, 0, 0, 0, 0, 0, 0, -0.500000, 0.500000, 0, 0, 0.500000, 0, \
279 0, 0, 0, 0, 0, 0, 0.500000, -0.500000, 0, 0, 0, 0, 0, 0, 0, \
280 -0.500000, 0, 0, -0.500000, 0.250000, 0, 0, 0.250000, -0.250000, 0, \
281 0, -0.250000, 0.375000, 0, 0, 0.250000, -0.125000, 0, 0, -0.250000, \
282 0.125000, 0, 0, 0.250000, -0.375000, 0, 0, -0.250000, 0.250000, 0, 0, \
283 0.125000, -0.250000, 0, 0, -0.375000, 0.250000, 0, 0, 0.375000, \
284 -0.250000, 0, 0, -0.125000
285 };
286
287 // CURL: correct values in (F,P) format
288 double basisCurls[] = {
289 0.25, 0.25, 0.25, 0.25,
290 0.25, 0.25, 0.25, 0.25,
291 0.25, 0.25, 0.25, 0.25,
292 0.25, 0.25, 0.25, 0.25,
293 0.25, 0.25, 0.25, 0.25,
294 0.25, 0.25, 0.25, 0.25,
295 0.25, 0.25, 0.25, 0.25,
296 0.25, 0.25, 0.25, 0.25,
297 0.25, 0.25, 0.25, 0.25
298 };
299
300
301 try{
302
303 // Dimensions for the output arrays:
304 int numFields = quadBasis.getCardinality();
305 int numPoints = quadNodes.dimension(0);
306 int spaceDim = quadBasis.getBaseCellTopology().getDimension();
307
308 // Generic array for values and curls that will be properly sized before each call
309 FieldContainer<double> vals;
310
311 // Check VALUE of basis functions: resize vals to rank-3 container:
312 vals.resize(numFields, numPoints, spaceDim);
313 quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE);
314 for (int i = 0; i < numFields; i++) {
315 for (int j = 0; j < numPoints; j++) {
316 for (int k = 0; k < spaceDim; k++) {
317
318 // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k
319 int l = k + i * spaceDim + j * spaceDim * numFields;
320 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
321 errorFlag++;
322 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
323
324 // Output the multi-index of the value where the error is:
325 *outStream << " At multi-index { ";
326 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
327 *outStream << "} computed value: " << vals(i,j,k)
328 << " but reference value: " << basisValues[l] << "\n";
329 }
330 }
331 }
332 }
333
334 // Check CURL of basis function: resize vals to rank-2 container
335 vals.resize(numFields, numPoints);
336 quadBasis.getValues(vals, quadNodes, OPERATOR_CURL);
337 for (int i = 0; i < numFields; i++) {
338 for (int j = 0; j < numPoints; j++) {
339 int l = i + j * numFields;
340 if (std::abs(vals(i,j) - basisCurls[l]) > INTREPID_TOL) {
341 errorFlag++;
342 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
343
344 // Output the multi-index of the value where the error is:
345 *outStream << " At multi-index { ";
346 *outStream << i << " ";*outStream << j << " ";
347 *outStream << "} computed curl component: " << vals(i,j)
348 << " but reference curl component: " << basisCurls[l] << "\n";
349 }
350 }
351 }
352 }
353
354 // Catch unexpected errors
355 catch (const std::logic_error & err) {
356 *outStream << err.what() << "\n\n";
357 errorFlag = -1000;
358 };
359
360 *outStream \
361 << "\n"
362 << "===============================================================================\n"\
363 << "| TEST 4: correctness of DoF locations |\n"\
364 << "===============================================================================\n";
365
366 try{
367 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
368 Teuchos::rcp(new Basis_HCURL_QUAD_I1_FEM<double, FieldContainer<double> >);
369 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
370 Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
371
372 int spaceDim = 2;
373 FieldContainer<double> cvals;
374 FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality(),2); // last dimension is spatial dim
375
376 // Check exceptions.
377#ifdef HAVE_INTREPID_DEBUG
378 cvals.resize(1,2,3);
379 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
380 cvals.resize(3,2);
381 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
382 cvals.resize(4,3);
383 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
384#endif
385 cvals.resize(4,spaceDim);
386 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
387 // Check if number of thrown exceptions matches the one we expect
388 if (throwCounter != nException) {
389 errorFlag++;
390 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
391 }
392
393 // Check mathematical correctness
394 FieldContainer<double> tangents(basis->getCardinality(),spaceDim); // tangents at each point basis point
395 tangents(0,0) = 2.0; tangents(0,1) = 0.0;
396 tangents(1,0) = 0.0; tangents(1,1) = 2.0;
397 tangents(2,0) = -2.0; tangents(2,1) = 0.0;
398 tangents(3,0) = 0.0; tangents(3,1) = -2.0;
399
400 basis->getValues(bvals, cvals, OPERATOR_VALUE);
401 char buffer[120];
402 for (int i=0; i<bvals.dimension(0); i++) {
403 for (int j=0; j<bvals.dimension(1); j++) {
404
405 double tangent = 0.0;
406 for(int d=0;d<spaceDim;d++)
407 tangent += bvals(i,j,d)*tangents(j,d);
408
409 if ((i != j) && (std::abs(tangent - 0.0) > INTREPID_TOL)) {
410 errorFlag++;
411 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), tangent, 0.0);
412 *outStream << buffer;
413 }
414 else if ((i == j) && (std::abs(tangent - 1.0) > INTREPID_TOL)) {
415 errorFlag++;
416 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), tangent, 1.0);
417 *outStream << buffer;
418 }
419 }
420 }
421
422 }
423 catch (const std::logic_error & err){
424 *outStream << err.what() << "\n\n";
425 errorFlag = -1000;
426 };
427
428
429 if (errorFlag != 0)
430 std::cout << "End Result: TEST FAILED\n";
431 else
432 std::cout << "End Result: TEST PASSED\n";
433
434 // reset format state of std::cout
435 std::cout.copyfmt(oldFormatState);
436
437 return errorFlag;
438}
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HCURL_QUAD_I1_FEM class.
Implementation of the default H(curl)-compatible FEM basis of degree 1 on Quadrilateral cell.