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
51#include "Teuchos_oblackholestream.hpp"
52#include "Teuchos_RCP.hpp"
53#include "Teuchos_GlobalMPISession.hpp"
54
55using namespace std;
56using namespace Intrepid;
57
58#define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
59{ \
60 ++nException; \
61 try { \
62 S ; \
63 } \
64 catch (const std::logic_error & err) { \
65 ++throwCounter; \
66 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
67 *outStream << err.what() << '\n'; \
68 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
69 }; \
70}
71
72int main(int argc, char *argv[]) {
73
74 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
75
76 // This little trick lets us print to std::cout only if
77 // a (dummy) command-line argument is provided.
78 int iprint = argc - 1;
79 Teuchos::RCP<std::ostream> outStream;
80 Teuchos::oblackholestream bhs; // outputs nothing
81 if (iprint > 0)
82 outStream = Teuchos::rcp(&std::cout, false);
83 else
84 outStream = Teuchos::rcp(&bhs, false);
85
86 // Save the format state of the original std::cout.
87 Teuchos::oblackholestream oldFormatState;
88 oldFormatState.copyfmt(std::cout);
89
90 *outStream \
91 << "===============================================================================\n" \
92 << "| |\n" \
93 << "| Unit Test (Basis_HCURL_HEX_In_FEM) |\n" \
94 << "| |\n" \
95 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
96 << "| 2) Basis values for VALUE and CURL 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
110 const int deg = 1;
111 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
112 FieldContainer<double> closedPts(PointTools::getLatticeSize(line,deg),1);
113 FieldContainer<double> openPts(PointTools::getLatticeSize(line,deg+1,1),1);
114 PointTools::getLattice<double,FieldContainer<double> >(closedPts,line,deg);
115 PointTools::getLattice<double,FieldContainer<double> >(openPts,line,deg+1,1);
116
117 Basis_HCURL_HEX_In_FEM<double, FieldContainer<double> > hexBasis(deg,closedPts,openPts);
118 int errorFlag = 0;
119
120 // Initialize throw counter for exception testing
121 int nException = 0;
122 int throwCounter = 0;
123
124 // Define array containing the 8 vertices of the reference HEX, its center and 6 face centers
125 FieldContainer<double> hexNodes(8, 3);
126 hexNodes(0,0) = -1.0; hexNodes(0,1) = -1.0; hexNodes(0,2) = -1.0;
127 hexNodes(1,0) = 1.0; hexNodes(1,1) = -1.0; hexNodes(1,2) = -1.0;
128 hexNodes(2,0) = -1.0; hexNodes(2,1) = 1.0; hexNodes(2,2) = -1.0;
129 hexNodes(3,0) = 1.0; hexNodes(3,1) = 1.0; hexNodes(3,2) = -1.0;
130 hexNodes(4,0) = -1.0; hexNodes(4,1) = -1.0; hexNodes(4,2) = 1.0;
131 hexNodes(5,0) = 1.0; hexNodes(5,1) = -1.0; hexNodes(5,2) = 1.0;
132 hexNodes(6,0) = -1.0; hexNodes(6,1) = 1.0; hexNodes(6,2) = 1.0;
133 hexNodes(7,0) = 1.0; hexNodes(7,1) = 1.0; hexNodes(7,2) = 1.0;
134
135
136 try{
137 // Generic array for the output values; needs to be properly resized depending on the operator type
138 FieldContainer<double> vals;
139
140 // exception #1: GRAD cannot be applied to HCURL functions
141 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
142 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0), 4 );
143 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD), throwCounter, nException );
144
145 // exception #2: DIV cannot be applied to HCURL functions
146 // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
147 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0) );
148 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_DIV), throwCounter, nException );
149
150 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
151 // getDofTag() to access invalid array elements thereby causing bounds check exception
152 // exception #3
153 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(3,0,0), throwCounter, nException );
154 // exception #4
155 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(1,1,1), throwCounter, nException );
156 // exception #5
157 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(0,4,1), throwCounter, nException );
158 // exception #6
159 INTREPID_TEST_COMMAND( hexBasis.getDofTag(12), throwCounter, nException );
160 // exception #7
161 INTREPID_TEST_COMMAND( hexBasis.getDofTag(-1), throwCounter, nException );
162
163#ifdef HAVE_INTREPID_DEBUG
164 // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays
165 // exception #8: input points array must be of rank-2
166 FieldContainer<double> badPoints1(4, 5, 3);
167 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
168
169 // exception #9 dimension 1 in the input point array must equal space dimension of the cell
170 FieldContainer<double> badPoints2(4, 2);
171 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
172
173 // exception #10 output values must be of rank-3 for OPERATOR_VALUE
174 FieldContainer<double> badVals1(4, 3);
175 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
176
177 // exception #11 output values must be of rank-3 for OPERATOR_CURL
178 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_CURL), throwCounter, nException );
179
180 // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
181 FieldContainer<double> badVals2(hexBasis.getCardinality() + 1, hexNodes.dimension(0), 3);
182 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
183
184 // exception #13 incorrect 1st dimension of output array (must equal number of points)
185 FieldContainer<double> badVals3(hexBasis.getCardinality(), hexNodes.dimension(0) + 1, 3);
186 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
187
188 // exception #14: incorrect 2nd dimension of output array (must equal the space dimension)
189 FieldContainer<double> badVals4(hexBasis.getCardinality(), hexNodes.dimension(0), 4);
190 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
191
192 // exception #15: incorrect 2nd dimension of output array (must equal the space dimension)
193 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_CURL), throwCounter, nException ) ;
194
195 // exception #16: D2 cannot be applied to HCURL functions
196 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
197 vals.resize(hexBasis.getCardinality(),
198 hexNodes.dimension(0),
199 Intrepid::getDkCardinality(OPERATOR_D2, hexBasis.getBaseCellTopology().getDimension()));
200 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_D2), throwCounter, nException );
201
202#endif
203
204 }
205 catch (const std::logic_error & err) {
206 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
207 *outStream << err.what() << '\n';
208 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
209 errorFlag = -1000;
210 };
211
212 // Check if number of thrown exceptions matches the one we expect
213 // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match.
214 if (throwCounter != nException) {
215 errorFlag++;
216 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
217 }
218//#endif
219
220 *outStream \
221 << "\n"
222 << "===============================================================================\n"\
223 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
224 << "===============================================================================\n";
225
226 try{
227 std::vector<std::vector<int> > allTags = hexBasis.getAllDofTags();
228
229 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
230 for (unsigned i = 0; i < allTags.size(); i++) {
231 int bfOrd = hexBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
232
233// for (unsigned j=0;j<4;j++) std::cout << allTags[i][j] << " "; std::cout << std::endl;
234
235 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
236 if( !( (myTag[0] == allTags[i][0]) &&
237 (myTag[1] == allTags[i][1]) &&
238 (myTag[2] == allTags[i][2]) &&
239 (myTag[3] == allTags[i][3]) ) ) {
240 errorFlag++;
241 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
242 *outStream << " getDofOrdinal( {"
243 << allTags[i][0] << ", "
244 << allTags[i][1] << ", "
245 << allTags[i][2] << ", "
246 << allTags[i][3] << "}) = " << bfOrd <<" but \n";
247 *outStream << " getDofTag(" << bfOrd << ") = { "
248 << myTag[0] << ", "
249 << myTag[1] << ", "
250 << myTag[2] << ", "
251 << myTag[3] << "}\n";
252 }
253 }
254
255 // Now do the same but loop over basis functions
256 for( int bfOrd = 0; bfOrd < hexBasis.getCardinality(); bfOrd++) {
257 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
258 int myBfOrd = hexBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
259 if( bfOrd != myBfOrd) {
260 errorFlag++;
261 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
262 *outStream << " getDofTag(" << bfOrd << ") = { "
263 << myTag[0] << ", "
264 << myTag[1] << ", "
265 << myTag[2] << ", "
266 << myTag[3] << "} but getDofOrdinal({"
267 << myTag[0] << ", "
268 << myTag[1] << ", "
269 << myTag[2] << ", "
270 << myTag[3] << "} ) = " << myBfOrd << "\n";
271 }
272 }
273 }
274 catch (const std::logic_error & err){
275 *outStream << err.what() << "\n\n";
276 errorFlag = -1000;
277 };
278
279 *outStream \
280 << "\n"
281 << "===============================================================================\n"\
282 << "| TEST 3: correctness of basis function values |\n"\
283 << "===============================================================================\n";
284
285 outStream -> precision(20);
286
287 // VALUE: Each row pair gives the 12x3 correct basis set values at an evaluation point: (P,F,D) layout
288 double basisValues[] = {
289 1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
290 0,0,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
291 0,0,0, 0,0,0, 0,0,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0,
292 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 1,0,0, 1,0,0,
293 0,1,0, 0,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
294 0,0,0, 0,1,0, 0,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
295 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,1,0, 0,0,0, 0,1,0, 0,0,0,
296 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,1,0, 0,0,0, 0,1,0,
297 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0,
298 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0,
299 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0,
300 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1
301 };
302
303 // CURL
304
305 double basisCurls[] = {
306 0,-0.5,0.5, 0,-0.5,0.5, 0,0,0.5, 0,0,0.5, 0,-0.5,0, 0,-0.5,0, 0,0,0, 0,0,0,
307 0,0,-0.5, 0,0,-0.5, 0,-0.5,-0.5, 0,-0.5,-0.5, 0,0,0, 0,0,0, 0,-0.5,0, 0,-0.5,0,
308 0,0.5,0, 0,0.5,0, 0,0,0, 0,0,0, 0,0.5,0.5, 0,0.5,0.5, 0,0,0.5, 0,0,0.5,
309 0,0,0, 0,0,0, 0,0.5,0, 0,0.5,0, 0,0,-0.5, 0,0,-0.5, 0,0.5,-0.5, 0,0.5,-0.5,
310 // y-component basis functions
311 // first y-component bf is (0,1/4(1-x)(1-z),0)
312 // curl is (1/4(1-x),0,-1/4(1-z))
313 0.5,0,-0.5, 0,0,-0.5, 0.5,0,-0.5, 0,0,-0.5, 0.5,0,0, 0,0,0, 0.5,0,0, 0,0,0,
314 // second y-component bf is (0,1/4(1+x)(1-z),0)
315 // curl is (1/4(1+x),0,1/4(1-z))
316 0,0,0.5, 0.5,0,0.5, 0,0,0.5, 0.5,0,0.5, 0,0,0, 0.5,0,0, 0,0,0, 0.5,0,0,
317 // third y-component bf is (0,1/4(1-x)(1+z),0)
318 // curl is (-1/4(1-x),0,-1/4(1+z))
319 -0.5,0,0, 0,0,0, -0.5,0,0, 0,0,0, -0.5,0,-0.5, 0,0,-0.5, -0.5,0,-0.5, 0,0,-0.5,
320 // fourth y-component bf is (0,1/4(1+x)(1+z),0)
321 // curl is (-1/4(1+x),0,1/4(1+z))
322 0.0,0,0, -0.5,0,0, 0.0,0,0, -0.5,0,0, 0.0,0,0.5, -0.5,0,0.5, 0.0,0,0.5, -0.5,0,0.5,
323 // first z-component bf is (0,0,1/4(1-x)(1-y))
324 // curl is (-1/4(1-x),1/4(1-y),0)
325 -0.5,0.5,0, 0,0.5,0, -0.5,0,0, 0,0,0, -0.5,0.5,0, 0,0.5,0, -0.5,0,0, 0,0,0,
326 // second z-component bf is (0,0,1/4(1+x)(1-y))
327 // curl is (-1/4(1+x),1/4(1-y),0)
328 0.0,-0.5,0, -0.5,-0.5,0, 0,0,0, -0.5,0,0, 0,-0.5,0, -0.5,-0.5,0, 0,0,0, -0.5,0,0,
329 // third z-component bf is (0,0,1/4(1-x)(1+y))
330 // curl is (1/4(1-x),1/4(1+y),0)
331 0.5,0,0, 0,0,0, 0.5,0.5,0, 0,0.5,0, 0.5,0,0, 0,0,0, 0.5,0.5,0, 0,0.5,0,
332 // fourth z-component bf is (0,0,1/4(1+x)(1+y))
333 // curl is (1/4(1+x),-1/4(1+y),0)
334 0,0,0, 0.5,0,0, 0,-0.5,0, 0.5,-0.5,0, 0,0,0, 0.5,0,0, 0,-0.5,0, 0.5,-0.5,0
335 };
336
337
338 try{
339
340 // Dimensions for the output arrays:
341 int numFields = hexBasis.getCardinality();
342 int numPoints = hexNodes.dimension(0);
343 int spaceDim = hexBasis.getBaseCellTopology().getDimension();
344
345 // Generic array for values and curls that will be properly sized before each call
346 FieldContainer<double> vals;
347
348 // Check VALUE of basis functions: resize vals to rank-3 container:
349 vals.resize(numFields, numPoints, spaceDim);
350 hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE);
351 for (int i = 0; i < numFields; i++) {
352 for (int j = 0; j < numPoints; j++) {
353 for (int k = 0; k < spaceDim; k++) {
354
355 // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k
356 int l = k + i * spaceDim * numPoints + j * spaceDim;
357 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
358 errorFlag++;
359 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
360
361 // Output the multi-index of the value where the error is:
362 *outStream << " At multi-index { ";
363 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
364 *outStream << "} computed value: " << vals(i,j,k)
365 << " but reference value: " << basisValues[l] << "\n";
366 }
367 }
368 }
369 }
370
371 // Check CURL of basis function: resize vals to rank-3 container
372 vals.resize(numFields, numPoints, spaceDim);
373 hexBasis.getValues(vals, hexNodes, OPERATOR_CURL);
374 for (int i = 0; i < numFields; i++) {
375 for (int j = 0; j < numPoints; j++) {
376 for (int k = 0; k < spaceDim; k++) {
377 // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k
378 int l = k + i * spaceDim * numPoints + j * spaceDim;
379 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
380 errorFlag++;
381 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
382
383 // Output the multi-index of the value where the error is:
384 *outStream << " At multi-index { ";
385 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
386 *outStream << "} computed curl component: " << vals(i,j,k)
387 << " but reference curl component: " << basisCurls[l] << "\n";
388 }
389 }
390 }
391 }
392
393 }
394
395 // Catch unexpected errors
396 catch (const std::logic_error & err) {
397 *outStream << err.what() << "\n\n";
398 errorFlag = -1000;
399 };
400
401 if (errorFlag != 0)
402 std::cout << "End Result: TEST FAILED\n";
403 else
404 std::cout << "End Result: TEST PASSED\n";
405
406 // reset format state of std::cout
407 std::cout.copyfmt(oldFormatState);
408
409 return errorFlag;
410}
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HCURL_HEX_In_FEM class.
Header file for utility class to provide point tools, such as barycentric coordinates,...
static int getLatticeSize(const shards::CellTopology &cellType, const int order, const int offset=0)
Computes the number of points in a lattice of a given order on a simplex (currently disabled for othe...