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