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
49#include "Intrepid_HGRAD_TET_C1_FEM.hpp"
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_HGRAD_TET_C1_FEM) |\n" \
93 << "| |\n" \
94 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 << "| 2) Basis values for VALUE, GRAD, CURL, and Dk 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 4 vertices of the reference TET and its center.
117 FieldContainer<double> tetNodes(8, 3);
118 tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0;
119 tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0;
120 tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0;
121 tetNodes(3,0) = 0.0; tetNodes(3,1) = 0.0; tetNodes(3,2) = 1.0;
122 tetNodes(4,0) = 0.25; tetNodes(4,1) = 0.5; tetNodes(4,2) = 0.1;
123 tetNodes(5,0) = 0.5; tetNodes(5,1) = 0.5; tetNodes(5,2) = 0.0;
124 tetNodes(6,0) = 0.5; tetNodes(6,1) = 0.0; tetNodes(6,2) = 0.5;
125 tetNodes(7,0) = 0.0; tetNodes(7,1) = 0.5; tetNodes(7,2) = 0.5;
126
127
128 // Generic array for the output values; needs to be properly resized depending on the operator type
130
131 try{
132 // exception #1: CURL cannot be applied to scalar functions
133 // resize vals to rank-3 container with dimensions (num. points, num. basis functions, arbitrary)
134 vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0), 4 );
135 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException );
136
137 // exception #2: DIV cannot be applied to scalar functions
138 // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
139 vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0) );
140 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_DIV), 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,0), throwCounter, nException );
150 // exception #6
151 INTREPID_TEST_COMMAND( tetBasis.getDofTag(5), throwCounter, nException );
152 // exception #7
153 INTREPID_TEST_COMMAND( tetBasis.getDofTag(-1), throwCounter, nException );
154
155#ifdef HAVE_INTREPID_DEBUG
156 // Exceptions 8-18 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, tetBasis.getBaseCellTopology().getDimension() - 1);
163 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
164
165 // exception #10 output values must be of rank-2 for OPERATOR_VALUE
166 FieldContainer<double> badVals1(4, 3, 1);
167 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
168
169 // exception #11 output values must be of rank-3 for OPERATOR_GRAD
170 FieldContainer<double> badVals2(4, 3);
171 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_GRAD), throwCounter, nException );
172
173 // exception #12 output values must be of rank-3 for OPERATOR_D1
174 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_D1), throwCounter, nException );
175
176 // exception #13 output values must be of rank-3 for OPERATOR_D2
177 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_D2), throwCounter, nException );
178
179 // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
180 FieldContainer<double> badVals3(tetBasis.getCardinality() + 1, tetNodes.dimension(0));
181 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
182
183 // exception #15 incorrect 1st dimension of output array (must equal number of points)
184 FieldContainer<double> badVals4(tetBasis.getCardinality(), tetNodes.dimension(0) + 1);
185 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals4, tetNodes, OPERATOR_VALUE), throwCounter, nException );
186
187 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
188 FieldContainer<double> badVals5(tetBasis.getCardinality(), tetNodes.dimension(0), tetBasis.getBaseCellTopology().getDimension() + 1);
189 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals5, tetNodes, OPERATOR_GRAD), throwCounter, nException );
190
191 // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D)
192 FieldContainer<double> badVals6(tetBasis.getCardinality(), tetNodes.dimension(0), 40);
193 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_D1), throwCounter, nException );
194
195 // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D)
196 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_D2), throwCounter, nException );
197#endif
198
199 }
200 catch (const std::logic_error & err) {
201 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
202 *outStream << err.what() << '\n';
203 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
204 errorFlag = -1000;
205 };
206
207 // Check if number of thrown exceptions matches the one we expect
208 if (throwCounter != nException) {
209 errorFlag++;
210 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
211 }
212
213 *outStream \
214 << "\n"
215 << "===============================================================================\n"\
216 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
217 << "===============================================================================\n";
218
219 try{
220 std::vector<std::vector<int> > allTags = tetBasis.getAllDofTags();
221
222 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
223 for (unsigned i = 0; i < allTags.size(); i++) {
224 int bfOrd = tetBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
225
226 std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
227 if( !( (myTag[0] == allTags[i][0]) &&
228 (myTag[1] == allTags[i][1]) &&
229 (myTag[2] == allTags[i][2]) &&
230 (myTag[3] == allTags[i][3]) ) ) {
231 errorFlag++;
232 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
233 *outStream << " getDofOrdinal( {"
234 << allTags[i][0] << ", "
235 << allTags[i][1] << ", "
236 << allTags[i][2] << ", "
237 << allTags[i][3] << "}) = " << bfOrd <<" but \n";
238 *outStream << " getDofTag(" << bfOrd << ") = { "
239 << myTag[0] << ", "
240 << myTag[1] << ", "
241 << myTag[2] << ", "
242 << myTag[3] << "}\n";
243 }
244 }
245
246 // Now do the same but loop over basis functions
247 for( int bfOrd = 0; bfOrd < tetBasis.getCardinality(); bfOrd++) {
248 std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
249 int myBfOrd = tetBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
250 if( bfOrd != myBfOrd) {
251 errorFlag++;
252 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
253 *outStream << " getDofTag(" << bfOrd << ") = { "
254 << myTag[0] << ", "
255 << myTag[1] << ", "
256 << myTag[2] << ", "
257 << myTag[3] << "} but getDofOrdinal({"
258 << myTag[0] << ", "
259 << myTag[1] << ", "
260 << myTag[2] << ", "
261 << myTag[3] << "} ) = " << myBfOrd << "\n";
262 }
263 }
264 }
265 catch (const std::logic_error & err){
266 *outStream << err.what() << "\n\n";
267 errorFlag = -1000;
268 };
269
270 *outStream \
271 << "\n"
272 << "===============================================================================\n"\
273 << "| TEST 3: correctness of basis function values |\n"\
274 << "===============================================================================\n";
275
276 outStream -> precision(20);
277
278 // VALUE: Each row gives the 4 correct basis set values at an evaluation point
279 double basisValues[] = {
280 1.0, 0.0, 0.0, 0.0,
281 0.0, 1.0, 0.0, 0.0,
282 0.0, 0.0, 1.0, 0.0,
283 0.0, 0.0, 0.0, 1.0,
284 0.15,0.25,0.5, 0.1,
285 0.0, 0.5, 0.5, 0.0,
286 0.0, 0.5, 0.0, 0.5,
287 0.0, 0.0, 0.5, 0.5
288 };
289
290 // GRAD and D1: each row gives the 3 x 4 correct values of the gradients of the 4 basis functions
291 double basisGrads[] = {
292 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
293 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
294 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
295 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
296 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
297 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
298 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
299 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
300 };
301
302 try{
303
304 // Dimensions for the output arrays:
305 int numFields = tetBasis.getCardinality();
306 int numPoints = tetNodes.dimension(0);
307 int spaceDim = tetBasis.getBaseCellTopology().getDimension();
308
309 // Generic array for values, grads, curls, etc. that will be properly sized before each call
311
312 // Check VALUE of basis functions: resize vals to rank-2 container:
313 vals.resize(numFields, numPoints);
314 tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE);
315 for (int i = 0; i < numFields; i++) {
316 for (int j = 0; j < numPoints; j++) {
317 int l = i + j * numFields;
318 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
319 errorFlag++;
320 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
321
322 // Output the multi-index of the value where the error is:
323 *outStream << " At multi-index { ";
324 *outStream << i << " ";*outStream << j << " ";
325 *outStream << "} computed value: " << vals(i,j)
326 << " but reference value: " << basisValues[l] << "\n";
327 }
328 }
329 }
330
331 // Check GRAD of basis function: resize vals to rank-3 container
332 vals.resize(numFields, numPoints, spaceDim);
333 tetBasis.getValues(vals, tetNodes, OPERATOR_GRAD);
334 for (int i = 0; i < numFields; i++) {
335 for (int j = 0; j < numPoints; j++) {
336 for (int k = 0; k < spaceDim; k++) {
337 int l = k + i * spaceDim + j * spaceDim * numFields;
338 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
339 errorFlag++;
340 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
341
342 // Output the multi-index of the value where the error is:
343 *outStream << " At multi-index { ";
344 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
345 *outStream << "} computed grad component: " << vals(i,j,k)
346 << " but reference grad component: " << basisGrads[l] << "\n";
347 }
348 }
349 }
350 }
351
352 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
353 tetBasis.getValues(vals, tetNodes, OPERATOR_D1);
354 for (int i = 0; i < numFields; i++) {
355 for (int j = 0; j < numPoints; j++) {
356 for (int k = 0; k < spaceDim; k++) {
357 int l = k + i * spaceDim + j * spaceDim * numFields;
358 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
359 errorFlag++;
360 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
361
362 // Output the multi-index of the value where the error is:
363 *outStream << " At multi-index { ";
364 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
365 *outStream << "} computed D1 component: " << vals(i,j,k)
366 << " but reference D1 component: " << basisGrads[l] << "\n";
367 }
368 }
369 }
370 }
371
372 // Check all higher derivatives - must be zero.
373 for(EOperator op = OPERATOR_D2; op < OPERATOR_MAX; op++) {
374
375 // The last dimension is the number of kth derivatives and needs to be resized for every Dk
376 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
377 vals.resize(numFields, numPoints, DkCardin);
378
379 tetBasis.getValues(vals, tetNodes, op);
380 for (int i = 0; i < vals.size(); i++) {
381 if (std::abs(vals[i]) > INTREPID_TOL) {
382 errorFlag++;
383 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
384
385 // Get the multi-index of the value where the error is and the operator order
386 std::vector<int> myIndex;
387 vals.getMultiIndex(myIndex,i);
388 int ord = Intrepid::getOperatorOrder(op);
389 *outStream << " At multi-index { ";
390 for(int j = 0; j < vals.rank(); j++) {
391 *outStream << myIndex[j] << " ";
392 }
393 *outStream << "} computed D"<< ord <<" component: " << vals[i]
394 << " but reference D" << ord << " component: 0 \n";
395 }
396 }
397 }
398 }
399
400 // Catch unexpected errors
401 catch (const std::logic_error & err) {
402 *outStream << err.what() << "\n\n";
403 errorFlag = -1000;
404 };
405
406 if (errorFlag != 0)
407 std::cout << "End Result: TEST FAILED\n";
408 else
409 std::cout << "End Result: TEST PASSED\n";
410
411 // reset format state of std::cout
412 std::cout.copyfmt(oldFormatState);
413
414 return errorFlag;
415}
Header file for utility class to provide multidimensional containers.
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.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Tetrahedron cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation 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....
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
void getMultiIndex(int &i0, const int valueEnum) const
Returns the multi-index of a value, based on its enumeration, as a list, for rank-1 containers.
int rank() const
Return rank of the FieldContainer = number of indices used to tag the multi-indexed value.