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_WEDGE_I1_FEM) |\n" \
93 << "| |\n" \
94 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 << "| 2) Basis values for VALUE and DIV 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
110 int errorFlag = 0;
111
112 // Initialize throw counter for exception testing
113 int nException = 0;
114 int throwCounter = 0;
115
116 // Define array containing the 6 vertices of the reference WEDGE and 6 other points.
117 FieldContainer<double> wedgeNodes(12, 3);
118 wedgeNodes(0,0) = 0.0; wedgeNodes(0,1) = 0.0; wedgeNodes(0,2) = -1.0;
119 wedgeNodes(1,0) = 1.0; wedgeNodes(1,1) = 0.0; wedgeNodes(1,2) = -1.0;
120 wedgeNodes(2,0) = 0.0; wedgeNodes(2,1) = 1.0; wedgeNodes(2,2) = -1.0;
121 wedgeNodes(3,0) = 0.0; wedgeNodes(3,1) = 0.0; wedgeNodes(3,2) = 1.0;
122 wedgeNodes(4,0) = 1.0; wedgeNodes(4,1) = 0.0; wedgeNodes(4,2) = 1.0;
123 wedgeNodes(5,0) = 0.0; wedgeNodes(5,1) = 1.0; wedgeNodes(5,2) = 1.0;
124
125 wedgeNodes(6,0) = 0.25; wedgeNodes(6,1) = 0.5; wedgeNodes(6,2) = -1.0;
126 wedgeNodes(7,0) = 0.5; wedgeNodes(7,1) = 0.25; wedgeNodes(7,2) = 0.0;
127 wedgeNodes(8,0) = 0.25; wedgeNodes(8,1) = 0.25; wedgeNodes(8,2) = 1.0;
128 wedgeNodes(9,0) = 0.25; wedgeNodes(9,1) = 0.0; wedgeNodes(9,2) = 0.75;
129 wedgeNodes(10,0)= 0.0; wedgeNodes(10,1)= 0.5; wedgeNodes(10,2)= -0.25;
130 wedgeNodes(11,0)= 0.5; wedgeNodes(11,1)= 0.5; wedgeNodes(11,2)= 0.0;
131
132
133 try{
134 // Generic array for the output values; needs to be properly resized depending on the operator type
136
137 // exception #1: GRAD cannot be applied to HDIV functions
138 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
139 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 );
140 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
141
142 // exception #2: CURL cannot be applied to HDIV functions
143 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_CURL), throwCounter, nException );
144
145 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
146 // getDofTag() to access invalid array elements thereby causing bounds check exception
147 // exception #3
148 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException );
149 // exception #4
150 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException );
151 // exception #5
152 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,4,1), throwCounter, nException );
153 // exception #6
154 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(11), throwCounter, nException );
155 // exception #7
156 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException );
157
158#ifdef HAVE_INTREPID_DEBUG
159 // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays
160 // exception #8: input points array must be of rank-2
161 FieldContainer<double> badPoints1(4, 5, 3);
162 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
163
164 // exception #9 dimension 1 in the input point array must equal space dimension of the cell
165 FieldContainer<double> badPoints2(4, 2);
166 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
167
168 // exception #10 output values must be of rank-3 for OPERATOR_VALUE
169 FieldContainer<double> badVals1(4, 3);
170 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
171
172 // exception #11 output values must be of rank-2 for OPERATOR_DIV
173 FieldContainer<double> badVals2(4, 3, 1);
174 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
175
176 // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
177 FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0), 3);
178 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
179
180 // exception #13 incorrect 0th dimension of output array (must equal number of basis functions)
181 FieldContainer<double> badVals4(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0));
182 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
183
184 // exception #14 incorrect 1st dimension of output array (must equal number of points)
185 FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1, 3);
186 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
187
188 // exception #15 incorrect 1st dimension of output array (must equal number of points)
189 FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1);
190 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
191
192 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
193 FieldContainer<double> badVals7(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 4);
194 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals7, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
195#endif
196
197 }
198 catch (const std::logic_error & err) {
199 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
200 *outStream << err.what() << '\n';
201 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
202 errorFlag = -1000;
203 };
204
205 // Check if number of thrown exceptions matches the one we expect
206 if (throwCounter != nException) {
207 errorFlag++;
208 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
209 }
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 = wedgeBasis.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 = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
223
224 std::vector<int> myTag = wedgeBasis.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 < wedgeBasis.getCardinality(); bfOrd++) {
246 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
247 int myBfOrd = wedgeBasis.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 5x3 correct basis set values at an evaluation point
277 double basisValues[] = {
278 0, -0.500000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, -2.00000, 0, 0, 0, \
279 0.500000, -0.500000, 0, 0.500000, 0, 0, 0, 0, 0, 0, 0, -2.00000, 0, \
280 0, 0, 0, 0, 0, 0, 0.500000, 0, -0.500000, 0.500000, 0, 0, 0, \
281 -2.00000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, \
282 0, 0, 0, 2.00000, 0.500000, -0.500000, 0, 0.500000, 0, 0, 0, 0, 0, 0, \
283 0, 0, 0, 0, 2.00000, 0, 0, 0, 0, 0.500000, 0, -0.500000, 0.500000, 0, \
284 0, 0, 0, 0, 0, 2.00000, 0.125000, -0.250000, 0, 0.125000, 0.250000, \
285 0, -0.375000, 0.250000, 0, 0, 0, -2.00000, 0, 0, 0, 0.250000, \
286 -0.375000, 0, 0.250000, 0.125000, 0, -0.250000, 0.125000, 0, 0, 0, \
287 -1.00000, 0, 0, 1.00000, 0.125000, -0.375000, 0, 0.125000, 0.125000, \
288 0, -0.375000, 0.125000, 0, 0, 0, 0, 0, 0, 2.00000, 0.125000, \
289 -0.500000, 0, 0.125000, 0, 0, -0.375000, 0, 0, 0, 0, -0.250000, 0, 0, \
290 1.75000, 0, -0.250000, 0, 0, 0.250000, 0, -0.500000, 0.250000, 0, 0, \
291 0, -1.25000, 0, 0, 0.750000, 0.250000, -0.250000, 0, 0.250000, \
292 0.250000, 0, -0.250000, 0.250000, 0, 0, 0, -1.00000, 0, 0, 1.00000};
293
294 // DIV: each row pair gives the 5 correct values of the divergence of the 5 basis functions
295 double basisDivs[] = {
296 // 6 vertices
297 1.0, 1.0, 1.0, 1.0, 1.0,
298 1.0, 1.0, 1.0, 1.0, 1.0,
299 1.0, 1.0, 1.0, 1.0, 1.0,
300 1.0, 1.0, 1.0, 1.0, 1.0,
301 1.0, 1.0, 1.0, 1.0, 1.0,
302 1.0, 1.0, 1.0, 1.0, 1.0,
303 // 6 other points
304 1.0, 1.0, 1.0, 1.0, 1.0,
305 1.0, 1.0, 1.0, 1.0, 1.0,
306 1.0, 1.0, 1.0, 1.0, 1.0,
307 1.0, 1.0, 1.0, 1.0, 1.0,
308 1.0, 1.0, 1.0, 1.0, 1.0,
309 1.0, 1.0, 1.0, 1.0, 1.0
310 };
311
312 try{
313
314 // Dimensions for the output arrays:
315 int numFields = wedgeBasis.getCardinality();
316 int numPoints = wedgeNodes.dimension(0);
317 int spaceDim = wedgeBasis.getBaseCellTopology().getDimension();
318
319 // Generic array for values and curls that will be properly sized before each call
321
322 // Check VALUE of basis functions: resize vals to rank-3 container:
323 vals.resize(numFields, numPoints, spaceDim);
324 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE);
325 for (int i = 0; i < numFields; i++) {
326 for (int j = 0; j < numPoints; j++) {
327 for (int k = 0; k < spaceDim; k++) {
328 int l = k + i * spaceDim + j * spaceDim * numFields;
329 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
330 errorFlag++;
331 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
332
333 // Output the multi-index of the value where the error is:
334 *outStream << " At multi-index { ";
335 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
336 *outStream << "} computed value: " << vals(i,j,k)
337 << " but reference value: " << basisValues[l] << "\n";
338 }
339 }
340 }
341 }
342
343
344 // Check DIV of basis function: resize vals to rank-2 container
345 vals.resize(numFields, numPoints);
346 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV);
347 for (int i = 0; i < numFields; i++) {
348 for (int j = 0; j < numPoints; j++) {
349 int l = i + j * numFields;
350 if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) {
351 errorFlag++;
352 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
353
354 // Output the multi-index of the value where the error is:
355 *outStream << " At multi-index { ";
356 *outStream << i << " ";*outStream << j << " ";
357 *outStream << "} computed divergence component: " << vals(i,j)
358 << " but reference divergence component: " << basisDivs[l] << "\n";
359 }
360 }
361 }
362
363 }
364
365 // Catch unexpected errors
366 catch (const std::logic_error & err) {
367 *outStream << err.what() << "\n\n";
368 errorFlag = -1000;
369 };
370
371 if (errorFlag != 0)
372 std::cout << "End Result: TEST FAILED\n";
373 else
374 std::cout << "End Result: TEST PASSED\n";
375
376 // reset format state of std::cout
377 std::cout.copyfmt(oldFormatState);
378
379 return errorFlag;
380}
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HDIV_WEDGE_I1_FEM class.
Implementation of the default H(div)-compatible FEM basis of degree 1 on Wedge cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Wedge 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.