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_HDIV_TET_I1_FEM) |\n" \
93 << "| |\n" \
94 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 << "| 3) Exception tests on input arguments and invalid operators |\n" \
96 << "| 2) Basis values for VALUE and DIV operators |\n" \
97 << "| |\n" \
98 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
99 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
100 << "| Kara Peterson (kjpeter@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
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 4 vertices of the reference TET and its 6 edge midpoints.
118 FieldContainer<double> tetNodes(10, 3);
119 tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0;
120 tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0;
121 tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0;
122 tetNodes(3,0) = 0.0; tetNodes(3,1) = 0.0; tetNodes(3,2) = 1.0;
123 tetNodes(4,0) = 0.5; tetNodes(4,1) = 0.0; tetNodes(4,2) = 0.0;
124 tetNodes(5,0) = 0.5; tetNodes(5,1) = 0.5; tetNodes(5,2) = 0.0;
125 tetNodes(6,0) = 0.0; tetNodes(6,1) = 0.5; tetNodes(6,2) = 0.0;
126 tetNodes(7,0) = 0.0; tetNodes(7,1) = 0.0; tetNodes(7,2) = 0.5;
127 tetNodes(8,0) = 0.5; tetNodes(8,1) = 0.0; tetNodes(8,2) = 0.5;
128 tetNodes(9,0) = 0.0; tetNodes(9,1) = 0.5; tetNodes(9,2) = 0.5;
129
130 // Generic array for the output values; needs to be properly resized depending on the operator type
132
133 try{
134 // exception #1: GRAD cannot be applied to HDIV functions
135 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
136 vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0), 3 );
137 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_GRAD), throwCounter, nException );
138
139 // exception #2: CURL cannot be applied to HDIV functions
140 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException );
141
142 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
143 // getDofTag() to access invalid array elements thereby causing bounds check exception
144 // exception #3
145 INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(3,0,0), throwCounter, nException );
146 // exception #4
147 INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(1,1,1), throwCounter, nException );
148 // exception #5
149 INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(0,4,1), throwCounter, nException );
150 // exception #6
151 INTREPID_TEST_COMMAND( tetBasis.getDofTag(7), throwCounter, nException );
152 // exception #7
153 INTREPID_TEST_COMMAND( tetBasis.getDofTag(-1), throwCounter, nException );
154
155#ifdef HAVE_INTREPID_DEBUG
156 // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays
157 // exception #8: input points array must be of rank-2
158 FieldContainer<double> badPoints1(4, 5, 3);
159 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
160
161 // exception #9 dimension 1 in the input point array must equal space dimension of the cell
162 FieldContainer<double> badPoints2(4, 2);
163 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
164
165 // exception #10 output values must be of rank-3 for OPERATOR_VALUE
166 FieldContainer<double> badVals1(4, 3);
167 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
168
169 // exception #11 output values must be of rank-2 for OPERATOR_DIV
170 FieldContainer<double> badVals2(4, 3, 1);
171 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_VALUE), throwCounter, nException );
172
173 // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
174 FieldContainer<double> badVals3(tetBasis.getCardinality() + 1, tetNodes.dimension(0), 3);
175 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
176
177 // exception #13 incorrect 0th dimension of output array (must equal number of basis functions)
178 FieldContainer<double> badVals4(tetBasis.getCardinality() + 1, tetNodes.dimension(0));
179 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals4, tetNodes, OPERATOR_DIV), throwCounter, nException );
180
181 // exception #14 incorrect 1st dimension of output array (must equal number of points)
182 FieldContainer<double> badVals5(tetBasis.getCardinality(), tetNodes.dimension(0) + 1, 3);
183 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals5, tetNodes, OPERATOR_VALUE), throwCounter, nException );
184
185 // exception #15 incorrect 1st dimension of output array (must equal number of points)
186 FieldContainer<double> badVals6(tetBasis.getCardinality(), tetNodes.dimension(0) + 1);
187 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_DIV), throwCounter, nException );
188
189 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
190 FieldContainer<double> badVals7(tetBasis.getCardinality(), tetNodes.dimension(0), 4);
191 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals7, tetNodes, OPERATOR_VALUE), throwCounter, nException );
192#endif
193
194 }
195 catch (const std::logic_error & err) {
196 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
197 *outStream << err.what() << '\n';
198 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
199 errorFlag = -1000;
200 };
201
202 // Check if number of thrown exceptions matches the one we expect
203 if (throwCounter != nException) {
204 errorFlag++;
205 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
206 }
207
208 *outStream \
209 << "\n"
210 << "===============================================================================\n"\
211 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
212 << "===============================================================================\n";
213
214 try{
215 std::vector<std::vector<int> > allTags = tetBasis.getAllDofTags();
216
217 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
218 for (unsigned i = 0; i < allTags.size(); i++) {
219 int bfOrd = tetBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
220
221 std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
222 if( !( (myTag[0] == allTags[i][0]) &&
223 (myTag[1] == allTags[i][1]) &&
224 (myTag[2] == allTags[i][2]) &&
225 (myTag[3] == allTags[i][3]) ) ) {
226 errorFlag++;
227 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
228 *outStream << " getDofOrdinal( {"
229 << allTags[i][0] << ", "
230 << allTags[i][1] << ", "
231 << allTags[i][2] << ", "
232 << allTags[i][3] << "}) = " << bfOrd <<" but \n";
233 *outStream << " getDofTag(" << bfOrd << ") = { "
234 << myTag[0] << ", "
235 << myTag[1] << ", "
236 << myTag[2] << ", "
237 << myTag[3] << "}\n";
238 }
239 }
240
241 // Now do the same but loop over basis functions
242 for( int bfOrd = 0; bfOrd < tetBasis.getCardinality(); bfOrd++) {
243 std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
244 int myBfOrd = tetBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
245 if( bfOrd != myBfOrd) {
246 errorFlag++;
247 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
248 *outStream << " getDofTag(" << bfOrd << ") = { "
249 << myTag[0] << ", "
250 << myTag[1] << ", "
251 << myTag[2] << ", "
252 << myTag[3] << "} but getDofOrdinal({"
253 << myTag[0] << ", "
254 << myTag[1] << ", "
255 << myTag[2] << ", "
256 << myTag[3] << "} ) = " << myBfOrd << "\n";
257 }
258 }
259 }
260 catch (const std::logic_error & err){
261 *outStream << err.what() << "\n\n";
262 errorFlag = -1000;
263 };
264
265 *outStream \
266 << "\n"
267 << "===============================================================================\n"\
268 << "| TEST 3: correctness of basis function values |\n"\
269 << "===============================================================================\n";
270
271 outStream -> precision(20);
272
273 // VALUE: Correct basis values in (P,F,D) format: each row gives the 4x3 correct basis set values
274 // at an evaluation point. Note that getValues returns results as an (F,P,D) array.
275 double basisValues[] = {
276 // 4 vertices
277 0.,-2.0,0., 0.,0.,0., -2.0,0.,0., 0.,0.,-2.0,
278 2.0,-2.0,0., 2.0,0.,0., 0.,0.,0., 2.0,0.,-2.0,
279 0.,0.,0., 0.,2.0,0., -2.0,2.0,0., 0,2.0,-2.0,
280 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, 0.,0.,0.,
281 // 6 edge midpoints
282 1.0,-2.0,0., 1.0,0.,0., -1.0,0.,0., 1.0,0.,-2.0,
283 1.0,-1.0,0., 1.0,1.0,0., -1.0,1.0,0., 1.0,1.0,-2.0,
284 0.,-1.0,0., 0.,1.0,0., -2.0,1.0,0., 0.,1.0,-2.0,
285 0.,-2.0,1.0, 0.,0.,1.0, -2.0,0.,1.0, 0.,0.,-1.0,
286 1.0,-2.0,1.0, 1.0,0.,1.0, -1.0,0.,1.0, 1.0,0.,-1.0,
287 0.,-1.0,1.0, 0.,1.0,1.0, -2.0,1.0,1.0, 0.,1.0,-1.0
288 // bf0 bf1 bf2 bf3
289 };
290
291 // DIV: each row gives the 4 correct values of the divergence of the 4 basis functions
292 double basisDivs[] = {
293 // 4 vertices
294 6.0, 6.0, 6.0, 6.0,
295 6.0, 6.0, 6.0, 6.0,
296 6.0, 6.0, 6.0, 6.0,
297 6.0, 6.0, 6.0, 6.0,
298 // 6 edge midpoints
299 6.0, 6.0, 6.0, 6.0,
300 6.0, 6.0, 6.0, 6.0,
301 6.0, 6.0, 6.0, 6.0,
302 6.0, 6.0, 6.0, 6.0,
303 6.0, 6.0, 6.0, 6.0,
304 6.0, 6.0, 6.0, 6.0
305 };
306
307 try{
308
309 // Dimensions for the output arrays:
310 int numFields = tetBasis.getCardinality();
311 int numPoints = tetNodes.dimension(0);
312 int spaceDim = tetBasis.getBaseCellTopology().getDimension();
313
314 // Generic array for values and curls that will be properly sized before each call
316
317 // Check VALUE of basis functions: resize vals to rank-3 container:
318 vals.resize(numFields, numPoints, spaceDim);
319 tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE);
320 for (int i = 0; i < numFields; i++) {
321 for (int j = 0; j < numPoints; j++) {
322 for (int k = 0; k < spaceDim; k++) {
323 // basisValues is (P,F,D) array so its multiindex is (j,i,k) and not (i,j,k)!
324 int l = k + i * spaceDim + j * spaceDim * numFields;
325 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
326 errorFlag++;
327 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
328
329 // Output the multi-index of the value where the error is:
330 *outStream << " At (Field,Point,Dim) multi-index { ";
331 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
332 *outStream << "} computed value: " << vals(i,j,k)
333 << " but reference value: " << basisValues[l] << "\n";
334 }
335 }
336 }
337 }
338
339
340 // Check DIV of basis function: resize vals to rank-2 container
341 vals.resize(numFields, numPoints);
342 tetBasis.getValues(vals, tetNodes, OPERATOR_DIV);
343 for (int i = 0; i < numFields; i++) {
344 for (int j = 0; j < numPoints; j++) {
345 int l = i + j * numFields;
346 if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) {
347 errorFlag++;
348 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
349
350 // Output the multi-index of the value where the error is:
351 *outStream << " At multi-index { ";
352 *outStream << i << " ";*outStream << j << " ";
353 *outStream << "} computed divergence component: " << vals(i,j)
354 << " but reference divergence component: " << basisDivs[l] << "\n";
355 }
356 }
357 }
358
359 }
360
361 // Catch unexpected errors
362 catch (const std::logic_error & err) {
363 *outStream << err.what() << "\n\n";
364 errorFlag = -1000;
365 };
366
367 *outStream \
368 << "\n"
369 << "===============================================================================\n"\
370 << "| TEST 4: correctness of DoF locations |\n"\
371 << "===============================================================================\n";
372
373 try{
374 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
375 Teuchos::rcp(new Basis_HDIV_TET_I1_FEM<double, FieldContainer<double> >);
376 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
377 Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
378
379 int spaceDim = 3;
381 FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality(),spaceDim); // last dimension is spatial dim
382
383 // Check exceptions.
384#ifdef HAVE_INTREPID_DEBUG
385 cvals.resize(1,2,3);
386 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
387 cvals.resize(3,2);
388 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
389 cvals.resize(4,2);
390 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
391#endif
392 cvals.resize(4,spaceDim);
393 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
394 // Check if number of thrown exceptions matches the one we expect
395 if (throwCounter != nException) {
396 errorFlag++;
397 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
398 }
399
400 // Check mathematical correctness
401 FieldContainer<double> normals(basis->getCardinality(),spaceDim); // normals at each point basis point
402 normals(0,0) = 0.0; normals(0,1) = -0.5; normals(0,2) = 0.0;
403 normals(1,0) = 0.5; normals(1,1) = 0.5; normals(1,2) = 0.5;
404 normals(2,0) = -0.5; normals(2,1) = 0.0; normals(2,2) = 0.0;
405 normals(3,0) = 0.0; normals(3,1) = 0.0; normals(3,2) = -0.5;
406
407 basis->getValues(bvals, cvals, OPERATOR_VALUE);
408 char buffer[120];
409 for (int i=0; i<bvals.dimension(0); i++) {
410 for (int j=0; j<bvals.dimension(1); j++) {
411
412 double normal = 0.0;
413 for(int d=0;d<spaceDim;d++)
414 normal += bvals(i,j,d)*normals(j,d);
415
416 if ((i != j) && (std::abs(normal - 0.0) > INTREPID_TOL)) {
417 errorFlag++;
418 sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), normal, 0.0);
419 *outStream << buffer;
420 }
421 else if ((i == j) && (std::abs(normal - 1.0) > INTREPID_TOL)) {
422 errorFlag++;
423 sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), normal, 1.0);
424 *outStream << buffer;
425 }
426 }
427 }
428
429 }
430 catch (const std::logic_error & err){
431 *outStream << err.what() << "\n\n";
432 errorFlag = -1000;
433 };
434
435 if (errorFlag != 0)
436 std::cout << "End Result: TEST FAILED\n";
437 else
438 std::cout << "End Result: TEST PASSED\n";
439
440 // reset format state of std::cout
441 std::cout.copyfmt(oldFormatState);
442
443 return errorFlag;
444}
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.
Header file for the Intrepid::HDIV_TET_I1_FEM class.
Implementation of the default H(div)-compatible FEM basis of degree 1 on Tetrahedron cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Tetrahedron 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.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.