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
48#include "Intrepid_FieldContainer.hpp"
49#include "Intrepid_HCURL_HEX_In_FEM.hpp"
50#include "Intrepid_PointTools.hpp"
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
137 // Generic array for the output values; needs to be properly resized depending on the operator type
138 FieldContainer<double> vals;
139
140 try{
141 // exception #1: GRAD cannot be applied to HCURL functions
142 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
143 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0), 4 );
144 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD), throwCounter, nException );
145
146 // exception #2: DIV cannot be applied to HCURL functions
147 // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
148 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0) );
149 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_DIV), throwCounter, nException );
150
151 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
152 // getDofTag() to access invalid array elements thereby causing bounds check exception
153 // exception #3
154 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(3,0,0), throwCounter, nException );
155 // exception #4
156 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(1,1,1), throwCounter, nException );
157 // exception #5
158 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(0,4,1), throwCounter, nException );
159 // exception #6
160 INTREPID_TEST_COMMAND( hexBasis.getDofTag(12), throwCounter, nException );
161 // exception #7
162 INTREPID_TEST_COMMAND( hexBasis.getDofTag(-1), throwCounter, nException );
163
164#ifdef HAVE_INTREPID_DEBUG
165 // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays
166 // exception #8: input points array must be of rank-2
167 FieldContainer<double> badPoints1(4, 5, 3);
168 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
169
170 // exception #9 dimension 1 in the input point array must equal space dimension of the cell
171 FieldContainer<double> badPoints2(4, 2);
172 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
173
174 // exception #10 output values must be of rank-3 for OPERATOR_VALUE
175 FieldContainer<double> badVals1(4, 3);
176 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
177
178 // exception #11 output values must be of rank-3 for OPERATOR_CURL
179 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_CURL), throwCounter, nException );
180
181 // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
182 FieldContainer<double> badVals2(hexBasis.getCardinality() + 1, hexNodes.dimension(0), 3);
183 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
184
185 // exception #13 incorrect 1st dimension of output array (must equal number of points)
186 FieldContainer<double> badVals3(hexBasis.getCardinality(), hexNodes.dimension(0) + 1, 3);
187 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
188
189 // exception #14: incorrect 2nd dimension of output array (must equal the space dimension)
190 FieldContainer<double> badVals4(hexBasis.getCardinality(), hexNodes.dimension(0), 4);
191 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
192
193 // exception #15: incorrect 2nd dimension of output array (must equal the space dimension)
194 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_CURL), throwCounter, nException ) ;
195
196 // exception #16: D2 cannot be applied to HCURL functions
197 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
198 vals.resize(hexBasis.getCardinality(),
199 hexNodes.dimension(0),
200 Intrepid::getDkCardinality(OPERATOR_D2, hexBasis.getBaseCellTopology().getDimension()));
201 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_D2), throwCounter, nException );
202
203#endif
204
205 }
206 catch (const std::logic_error & err) {
207 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
208 *outStream << err.what() << '\n';
209 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
210 errorFlag = -1000;
211 };
212
213 // Check if number of thrown exceptions matches the one we expect
214 // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match.
215 if (throwCounter != nException) {
216 errorFlag++;
217 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
218 }
219//#endif
220
221 *outStream \
222 << "\n"
223 << "===============================================================================\n"\
224 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
225 << "===============================================================================\n";
226
227 try{
228 std::vector<std::vector<int> > allTags = hexBasis.getAllDofTags();
229
230 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
231 for (unsigned i = 0; i < allTags.size(); i++) {
232 int bfOrd = hexBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
233
234// for (unsigned j=0;j<4;j++) std::cout << allTags[i][j] << " "; std::cout << std::endl;
235
236 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
237 if( !( (myTag[0] == allTags[i][0]) &&
238 (myTag[1] == allTags[i][1]) &&
239 (myTag[2] == allTags[i][2]) &&
240 (myTag[3] == allTags[i][3]) ) ) {
241 errorFlag++;
242 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
243 *outStream << " getDofOrdinal( {"
244 << allTags[i][0] << ", "
245 << allTags[i][1] << ", "
246 << allTags[i][2] << ", "
247 << allTags[i][3] << "}) = " << bfOrd <<" but \n";
248 *outStream << " getDofTag(" << bfOrd << ") = { "
249 << myTag[0] << ", "
250 << myTag[1] << ", "
251 << myTag[2] << ", "
252 << myTag[3] << "}\n";
253 }
254 }
255
256 // Now do the same but loop over basis functions
257 for( int bfOrd = 0; bfOrd < hexBasis.getCardinality(); bfOrd++) {
258 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
259 int myBfOrd = hexBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
260 if( bfOrd != myBfOrd) {
261 errorFlag++;
262 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
263 *outStream << " getDofTag(" << bfOrd << ") = { "
264 << myTag[0] << ", "
265 << myTag[1] << ", "
266 << myTag[2] << ", "
267 << myTag[3] << "} but getDofOrdinal({"
268 << myTag[0] << ", "
269 << myTag[1] << ", "
270 << myTag[2] << ", "
271 << myTag[3] << "} ) = " << myBfOrd << "\n";
272 }
273 }
274 }
275 catch (const std::logic_error & err){
276 *outStream << err.what() << "\n\n";
277 errorFlag = -1000;
278 };
279
280 *outStream \
281 << "\n"
282 << "===============================================================================\n"\
283 << "| TEST 3: correctness of basis function values |\n"\
284 << "===============================================================================\n";
285
286 outStream -> precision(20);
287
288 // VALUE: Each row pair gives the 12x3 correct basis set values at an evaluation point: (P,F,D) layout
289 double basisValues[] = {
290 1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
291 0,0,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
292 0,0,0, 0,0,0, 0,0,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0,
293 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 1,0,0, 1,0,0,
294 0,1,0, 0,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
295 0,0,0, 0,1,0, 0,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
296 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,1,0, 0,0,0, 0,1,0, 0,0,0,
297 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,1,0, 0,0,0, 0,1,0,
298 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0,
299 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0,
300 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0,
301 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1
302 };
303
304 // CURL
305
306 double basisCurls[] = {
307 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,
308 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,
309 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,
310 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,
311 // y-component basis functions
312 // first y-component bf is (0,1/4(1-x)(1-z),0)
313 // curl is (1/4(1-x),0,-1/4(1-z))
314 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,
315 // second y-component bf is (0,1/4(1+x)(1-z),0)
316 // curl is (1/4(1+x),0,1/4(1-z))
317 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,
318 // third y-component bf is (0,1/4(1-x)(1+z),0)
319 // curl is (-1/4(1-x),0,-1/4(1+z))
320 -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,
321 // fourth y-component bf is (0,1/4(1+x)(1+z),0)
322 // curl is (-1/4(1+x),0,1/4(1+z))
323 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,
324 // first z-component bf is (0,0,1/4(1-x)(1-y))
325 // curl is (-1/4(1-x),1/4(1-y),0)
326 -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,
327 // second z-component bf is (0,0,1/4(1+x)(1-y))
328 // curl is (-1/4(1+x),1/4(1-y),0)
329 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,
330 // third z-component bf is (0,0,1/4(1-x)(1+y))
331 // curl is (1/4(1-x),1/4(1+y),0)
332 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,
333 // fourth z-component bf is (0,0,1/4(1+x)(1+y))
334 // curl is (1/4(1+x),-1/4(1+y),0)
335 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
336 };
337
338
339 try{
340
341 // Dimensions for the output arrays:
342 int numFields = hexBasis.getCardinality();
343 int numPoints = hexNodes.dimension(0);
344 int spaceDim = hexBasis.getBaseCellTopology().getDimension();
345
346 // Generic array for values and curls that will be properly sized before each call
347 FieldContainer<double> vals;
348
349 // Check VALUE of basis functions: resize vals to rank-3 container:
350 vals.resize(numFields, numPoints, spaceDim);
351 hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE);
352 for (int i = 0; i < numFields; i++) {
353 for (int j = 0; j < numPoints; j++) {
354 for (int k = 0; k < spaceDim; k++) {
355
356 // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k
357 int l = k + i * spaceDim * numPoints + j * spaceDim;
358 if (std::abs(vals(i,j,k) - basisValues[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 value: " << vals(i,j,k)
366 << " but reference value: " << basisValues[l] << "\n";
367 }
368 }
369 }
370 }
371
372 // Check CURL of basis function: resize vals to rank-3 container
373 vals.resize(numFields, numPoints, spaceDim);
374 hexBasis.getValues(vals, hexNodes, OPERATOR_CURL);
375 for (int i = 0; i < numFields; i++) {
376 for (int j = 0; j < numPoints; j++) {
377 for (int k = 0; k < spaceDim; k++) {
378 // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k
379 int l = k + i * spaceDim * numPoints + j * spaceDim;
380 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
381 errorFlag++;
382 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
383
384 // Output the multi-index of the value where the error is:
385 *outStream << " At multi-index { ";
386 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
387 *outStream << "} computed curl component: " << vals(i,j,k)
388 << " but reference curl component: " << basisCurls[l] << "\n";
389 }
390 }
391 }
392 }
393
394 }
395
396 // Catch unexpected errors
397 catch (const std::logic_error & err) {
398 *outStream << err.what() << "\n\n";
399 errorFlag = -1000;
400 };
401
402 if (errorFlag != 0)
403 std::cout << "End Result: TEST FAILED\n";
404 else
405 std::cout << "End Result: TEST PASSED\n";
406
407 // reset format state of std::cout
408 std::cout.copyfmt(oldFormatState);
409
410 return errorFlag;
411}
main
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
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid_HCURL_HEX_In_FEM.hpp
Header file for the Intrepid::HCURL_HEX_In_FEM class.
Intrepid_PointTools.hpp
Header file for utility class to provide point tools, such as barycentric coordinates,...
Intrepid::Basis_HCURL_HEX_In_FEM
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedral cell.
Definition: Intrepid_HCURL_HEX_In_FEM.hpp:71
Intrepid::FieldContainer::resize
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
Definition: Intrepid_FieldContainerDef.hpp:1137
Generated by doxygen 1.9.4