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_HGRAD_HEX_C2_FEM) |\n" \
94 << "| |\n" \
95 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
96 << "| 2) Basis values for VALUE, GRAD, and Dk operators |\n" \
97 << "| |\n" \
98 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
99 << "| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
100 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
101 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
102 << "| |\n" \
103 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
104 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
105 << "| |\n" \
106 << "===============================================================================\n"\
107 << "| TEST 1: Basis creation, exception testing |\n"\
108 << "===============================================================================\n";
109
110 // Define basis and error flag
111 // get points
112 const int deg=2;
113 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
114 FieldContainer<double> pts(PointTools::getLatticeSize(line,deg),1);
115 PointTools::getLattice<double,FieldContainer<double> >(pts,line,deg);
116
117 Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> > hexBasis(deg,deg,deg,pts,pts,pts);
118
119 int errorFlag = 0;
120
121 // Initialize throw counter for exception testing
122 int nException = 0;
123 int throwCounter = 0;
124
125 // Define arrayS containing the 27 nodes of hexahedron<27> topology
126 FieldContainer<double> hexNodes(27, 3);
127 // do it lexicographically as a lattice
128 hexNodes(0, 0) = -1.0; hexNodes(0, 1) = -1.0; hexNodes(0, 2) = -1.0;
129 hexNodes(1, 0) = 0.0; hexNodes(1, 1) = -1.0; hexNodes(1, 2) = -1.0;
130 hexNodes(2, 0) = 1.0; hexNodes(2, 1) = -1.0; hexNodes(2, 2) = -1.0;
131 hexNodes(3, 0) = -1.0; hexNodes(3, 1) = 0.0; hexNodes(3, 2) = -1.0;
132 hexNodes(4, 0) = 0.0; hexNodes(4, 1) = 0.0; hexNodes(4, 2) = -1.0;
133 hexNodes(5, 0) = 1.0; hexNodes(5, 1) = 0.0; hexNodes(5, 2) = -1.0;
134 hexNodes(6, 0) = -1.0; hexNodes(6, 1) = 1.0; hexNodes(6, 2) = -1.0;
135 hexNodes(7, 0) = 0.0; hexNodes(7, 1) = 1.0; hexNodes(7, 2) = -1.0;
136 hexNodes(8, 0) = 1.0; hexNodes(8, 1) = 1.0; hexNodes(8, 2) = -1.0;
137 hexNodes(9, 0) = -1.0; hexNodes(9, 1) = -1.0; hexNodes(9, 2) = 0.0;
138 hexNodes(10, 0) = 0.0; hexNodes(10, 1) = -1.0; hexNodes(10, 2) = 0.0;
139 hexNodes(11, 0) = 1.0; hexNodes(11, 1) = -1.0; hexNodes(11, 2) = 0.0;
140 hexNodes(12, 0) = -1.0; hexNodes(12, 1) = 0.0; hexNodes(12, 2) = 0.0;
141 hexNodes(13, 0) = 0.0; hexNodes(13, 1) = 0.0; hexNodes(13, 2) = 0.0;
142 hexNodes(14, 0) = 1.0; hexNodes(14, 1) = 0.0; hexNodes(14, 2) = 0.0;
143 hexNodes(15, 0) = -1.0; hexNodes(15, 1) = 1.0; hexNodes(15, 2) = 0.0;
144 hexNodes(16, 0) = 0.0; hexNodes(16, 1) = 1.0; hexNodes(16, 2) = 0.0;
145 hexNodes(17, 0) = 1.0; hexNodes(17, 1) = 1.0; hexNodes(17, 2) = 0.0;
146 hexNodes(18, 0) = -1.0; hexNodes(18, 1) = -1.0; hexNodes(18, 2) = 1.0;
147 hexNodes(19, 0) = 0.0; hexNodes(19, 1) = -1.0; hexNodes(19, 2) = 1.0;
148 hexNodes(20, 0) = 1.0; hexNodes(20, 1) = -1.0; hexNodes(20, 2) = 1.0;
149 hexNodes(21, 0) = -1.0; hexNodes(21, 1) = 0.0; hexNodes(21, 2) = 1.0;
150 hexNodes(22, 0) = 0.0; hexNodes(22, 1) = 0.0; hexNodes(22, 2) = 1.0;
151 hexNodes(23, 0) = 1.0; hexNodes(23, 1) = 0.0; hexNodes(23, 2) = 1.0;
152 hexNodes(24, 0) = -1.0; hexNodes(24, 1) = 1.0; hexNodes(24, 2) = 1.0;
153 hexNodes(25, 0) = 0.0; hexNodes(25, 1) = 1.0; hexNodes(25, 2) = 1.0;
154 hexNodes(26, 0) = 1.0; hexNodes(26, 1) = 1.0; hexNodes(26, 2) = 1.0;
155
156
157 // Generic array for the output values; needs to be properly resized depending on the operator type
158 FieldContainer<double> vals;
159
160 try{
161 // exception #1: CURL cannot be applied to scalar functions in 3D
162 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
163 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0), 4 );
164 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_CURL), throwCounter, nException );
165
166 // exception #2: DIV cannot be applied to scalar functions in 3D
167 // resize vals to rank-2 container with dimensions (num. basis functions, num. points)
168 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0) );
169 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_DIV), throwCounter, nException );
170
171 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
172 // getDofTag() to access invalid array elements thereby causing bounds check exception
173 // exception #3
174 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(3,10,0), throwCounter, nException );
175 // exception #4
176 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(1,2,1), throwCounter, nException );
177 // exception #5
178 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(0,4,1), throwCounter, nException );
179 // exception #6
180 INTREPID_TEST_COMMAND( hexBasis.getDofTag(28), throwCounter, nException );
181 // exception #7
182 INTREPID_TEST_COMMAND( hexBasis.getDofTag(-1), throwCounter, nException );
183
184#ifdef HAVE_INTREPID_DEBUG
185 // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays
186 // exception #8: input points array must be of rank-2
187 FieldContainer<double> badPoints1(4, 5, 3);
188 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
189
190 // exception #9 dimension 1 in the input point array must equal space dimension of the cell
191 FieldContainer<double> badPoints2(4, hexBasis.getBaseCellTopology().getDimension() - 1);
192 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
193
194 // exception #10 output values must be of rank-2 for OPERATOR_VALUE
195 FieldContainer<double> badVals1(4, 3, 1);
196 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
197
198 // exception #11 output values must be of rank-3 for OPERATOR_GRAD
199 FieldContainer<double> badVals2(4, 3);
200 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_GRAD), throwCounter, nException );
201
202 // exception #12 output values must be of rank-3 for OPERATOR_D1
203 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_D1), throwCounter, nException );
204
205 // exception #13 output values must be of rank-3 for OPERATOR_D2
206 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_D2), throwCounter, nException );
207
208 // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
209 FieldContainer<double> badVals3(hexBasis.getCardinality() + 1, hexNodes.dimension(0));
210 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException );
211
212 // exception #15 incorrect 1st dimension of output array (must equal number of points)
213 FieldContainer<double> badVals4(hexBasis.getCardinality(), hexNodes.dimension(0) + 1);
214 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException );
215
216 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
217 FieldContainer<double> badVals5(hexBasis.getCardinality(), hexNodes.dimension(0), hexBasis.getBaseCellTopology().getDimension() - 1);
218 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals5, hexNodes, OPERATOR_GRAD), throwCounter, nException );
219
220 // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D)
221 FieldContainer<double> badVals6(hexBasis.getCardinality(), hexNodes.dimension(0), 40);
222 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals6, hexNodes, OPERATOR_D2), throwCounter, nException );
223
224 // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D)
225 FieldContainer<double> badVals7(hexBasis.getCardinality(), hexNodes.dimension(0), 50);
226 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals7, hexNodes, OPERATOR_D3), throwCounter, nException );
227#endif
228
229 }
230 catch (const std::logic_error & err) {
231 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
232 *outStream << err.what() << '\n';
233 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
234 errorFlag = -1000;
235 };
236
237 // Check if number of thrown exceptions matches the one we expect
238 // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match.
239 if (throwCounter != nException) {
240 errorFlag++;
241 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
242 }
243
244 *outStream \
245 << "\n"
246 << "===============================================================================\n"\
247 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
248 << "===============================================================================\n";
249
250 try{
251 std::vector<std::vector<int> > allTags = hexBasis.getAllDofTags();
252
253 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
254 for (unsigned i = 0; i < allTags.size(); i++) {
255 int bfOrd = hexBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
256
257
258
259 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
260 if( !( (myTag[0] == allTags[i][0]) &&
261 (myTag[1] == allTags[i][1]) &&
262 (myTag[2] == allTags[i][2]) &&
263 (myTag[3] == allTags[i][3]) ) ) {
264 errorFlag++;
265 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
266 *outStream << " getDofOrdinal( {"
267 << allTags[i][0] << ", "
268 << allTags[i][1] << ", "
269 << allTags[i][2] << ", "
270 << allTags[i][3] << "}) = " << bfOrd <<" but \n";
271 *outStream << " getDofTag(" << bfOrd << ") = { "
272 << myTag[0] << ", "
273 << myTag[1] << ", "
274 << myTag[2] << ", "
275 << myTag[3] << "}\n";
276 }
277 }
278
279 // Now do the same but loop over basis functions
280 for( int bfOrd = 0; bfOrd < hexBasis.getCardinality(); bfOrd++) {
281 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
282 int myBfOrd = hexBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
283 if( bfOrd != myBfOrd) {
284 errorFlag++;
285 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
286 *outStream << " getDofTag(" << bfOrd << ") = { "
287 << myTag[0] << ", "
288 << myTag[1] << ", "
289 << myTag[2] << ", "
290 << myTag[3] << "} but getDofOrdinal({"
291 << myTag[0] << ", "
292 << myTag[1] << ", "
293 << myTag[2] << ", "
294 << myTag[3] << "} ) = " << myBfOrd << "\n";
295 }
296 }
297 }
298 catch (const std::logic_error & err){
299 *outStream << err.what() << "\n\n";
300 errorFlag = -1000;
301 };
302
303
304 *outStream \
305 << "\n"
306 << "===============================================================================\n"\
307 << "| TEST 3: correctness of basis function values |\n"\
308 << "===============================================================================\n";
309
310 outStream -> precision(20);
311
312 // VALUE: Each row gives the 27 correct basis set values at an evaluation point
313 double basisValues[] = {
314 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
315 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
316 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
317 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
318 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
319 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
320 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
321 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
322 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
323 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
324 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
325 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
326 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
327 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
328 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
329 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
330 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
331 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
332 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, \
333 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, \
334 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, \
335 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, \
336 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, \
337 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, \
338 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, 0, \
339 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000, 0, \
340 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.000 };
341
342
343 // GRAD, D1, D2, D3 and D4 test values are stored in files due to their large size
344 std::string fileName;
345 std::ifstream dataFile;
346
347 // GRAD and D1 values are stored in (F,P,D) format in a data file. Read file and do the test
348 std::vector<double> basisGrads; // Flat array for the gradient values.
349
350 fileName = "./testdata/HEX_C2_GradVals.dat";
351 dataFile.open(fileName.c_str());
352 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
353 ">>> ERROR (HGRAD_HEX_C2/test01): could not open GRAD values data file, test aborted.");
354 while (!dataFile.eof() ){
355 double temp;
356 string line; // string for one line of input file
357 std::getline(dataFile, line); // get next line from file
358 stringstream data_line(line); // convert to stringstream
359 while(data_line >> temp){ // extract value from line
360 basisGrads.push_back(temp); // push into vector
361 }
362 }
363 // It turns out that just closing and then opening the ifstream variable does not reset it
364 // and subsequent open() command fails. One fix is to explicitely clear the ifstream, or
365 // scope the variables.
366 dataFile.close();
367 dataFile.clear();
368
369
370 //D2: flat array with the values of D2 applied to basis functions. Multi-index is (F,P,D2cardinality)
371 std::vector<double> basisD2;
372 fileName = "./testdata/HEX_C2_D2Vals.dat";
373 dataFile.open(fileName.c_str());
374 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
375 ">>> ERROR (HGRAD_HEX_C2/test01): could not open D2 values data file, test aborted.");
376 while (!dataFile.eof() ){
377 double temp;
378 string line; // string for one line of input file
379 std::getline(dataFile, line); // get next line from file
380 stringstream data_line(line); // convert to stringstream
381 while(data_line >> temp){ // extract value from line
382 basisD2.push_back(temp); // push into vector
383 }
384 }
385 dataFile.close();
386 dataFile.clear();
387
388
389 //D3: flat array with the values of D3 applied to basis functions. Multi-index is (F,P,D3cardinality)
390 std::vector<double> basisD3;
391
392 fileName = "./testdata/HEX_C2_D3Vals.dat";
393 dataFile.open(fileName.c_str());
394 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
395 ">>> ERROR (HGRAD_HEX_C2/test01): could not open D3 values data file, test aborted.");
396
397 while (!dataFile.eof() ){
398 double temp;
399 string line; // string for one line of input file
400 std::getline(dataFile, line); // get next line from file
401 stringstream data_line(line); // convert to stringstream
402 while(data_line >> temp){ // extract value from line
403 basisD3.push_back(temp); // push into vector
404 }
405 }
406 dataFile.close();
407 dataFile.clear();
408
409
410 //D4: flat array with the values of D applied to basis functions. Multi-index is (F,P,D4cardinality)
411 std::vector<double> basisD4;
412
413 fileName = "./testdata/HEX_C2_D4Vals.dat";
414 dataFile.open(fileName.c_str());
415 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
416 ">>> ERROR (HGRAD_HEX_C2/test01): could not open D4 values data file, test aborted.");
417
418 while (!dataFile.eof() ){
419 double temp;
420 string line; // string for one line of input file
421 std::getline(dataFile, line); // get next line from file
422 stringstream data_line(line); // convert to stringstream
423 while(data_line >> temp){ // extract value from line
424 basisD4.push_back(temp); // push into vector
425 }
426 }
427 dataFile.close();
428 dataFile.clear();
429
430
431 try{
432
433 // Dimensions for the output arrays:
434 int numFields = hexBasis.getCardinality();
435 int numPoints = hexNodes.dimension(0);
436 int spaceDim = hexBasis.getBaseCellTopology().getDimension();
437
438 // Generic array for values, grads, curls, etc. that will be properly sized before each call
439 FieldContainer<double> vals;
440
441 // Check VALUE of basis functions: resize vals to rank-2 container:
442 vals.resize(numFields, numPoints);
443 hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE);
444 for (int i = 0; i < numFields; i++) {
445 for (int j = 0; j < numPoints; j++) {
446 int l = i + j * numFields;
447 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
448 errorFlag++;
449 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
450
451 // Output the multi-index of the value where the error is:
452 *outStream << " At multi-index { ";
453 *outStream << i << " ";*outStream << j << " ";
454 *outStream << "} computed value: " << vals(i,j)
455 << " but reference value: " << basisValues[l] << "\n";
456 }
457 }
458 }
459
460 // Check GRAD of basis function: resize vals to rank-3 container
461 vals.resize(numFields, numPoints, spaceDim);
462 hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD);
463 for (int i = 0; i < numFields; i++) {
464 for (int j = 0; j < numPoints; j++) {
465 for (int k = 0; k < spaceDim; k++) {
466
467 // basisGrads is (F,P,D), compute offset:
468 int l = k + j * spaceDim + i * spaceDim * numPoints;
469 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
470 errorFlag++;
471 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
472
473 // Output the multi-index of the value where the error is:
474 *outStream << " At multi-index { ";
475 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
476 *outStream << "} computed grad component: " << vals(i,j,k)
477 << " but reference grad component: " << basisGrads[l] << "\n";
478 }
479 }
480 }
481 }
482
483 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
484 hexBasis.getValues(vals, hexNodes, OPERATOR_D1);
485 for (int i = 0; i < numFields; i++) {
486 for (int j = 0; j < numPoints; j++) {
487 for (int k = 0; k < spaceDim; k++) {
488
489 // basisGrads is (F,P,D), compute offset:
490 int l = k + j * spaceDim + i * spaceDim * numPoints;
491 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
492 errorFlag++;
493 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
494
495 // Output the multi-index of the value where the error is:
496 *outStream << " At multi-index { ";
497 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
498 *outStream << "} computed D1 component: " << vals(i,j,k)
499 << " but reference D1 component: " << basisGrads[l] << "\n";
500 }
501 }
502 }
503 }
504
505
506 // Check D2 of basis function
507 int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
508 vals.resize(numFields, numPoints, D2cardinality);
509 hexBasis.getValues(vals, hexNodes, OPERATOR_D2);
510 for (int i = 0; i < numFields; i++) {
511 for (int j = 0; j < numPoints; j++) {
512 for (int k = 0; k < D2cardinality; k++) {
513
514 // basisD2 is (F,P,Dk), compute offset:
515 int l = k + j * D2cardinality + i * D2cardinality * numPoints;
516 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
517 errorFlag++;
518 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
519
520 // Output the multi-index of the value where the error is:
521 *outStream << " At multi-index { ";
522 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
523 *outStream << "} computed D2 component: " << vals(i,j,k)
524 << " but reference D2 component: " << basisD2[l] << "\n";
525 }
526 }
527 }
528 }
529
530
531 // Check D3 of basis function
532 int D3cardinality = Intrepid::getDkCardinality(OPERATOR_D3, spaceDim);
533 vals.resize(numFields, numPoints, D3cardinality);
534 hexBasis.getValues(vals, hexNodes, OPERATOR_D3);
535
536 for (int i = 0; i < numFields; i++) {
537 for (int j = 0; j < numPoints; j++) {
538 for (int k = 0; k < D3cardinality; k++) {
539
540 // basisD3 is (F,P,Dk), compute offset:
541 int l = k + j * D3cardinality + i * D3cardinality * numPoints;
542 if (std::abs(vals(i,j,k) - basisD3[l]) > INTREPID_TOL) {
543 errorFlag++;
544 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
545
546 // Output the multi-index of the value where the error is:
547 *outStream << " At multi-index { ";
548 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
549 *outStream << "} computed D3 component: " << vals(i,j,k)
550 << " but reference D3 component: " << basisD3[l] << "\n";
551 }
552 }
553 }
554 }
555
556
557 // Check D4 of basis function
558 int D4cardinality = Intrepid::getDkCardinality(OPERATOR_D4, spaceDim);
559 vals.resize(numFields, numPoints, D4cardinality);
560 hexBasis.getValues(vals, hexNodes, OPERATOR_D4);
561 for (int i = 0; i < numFields; i++) {
562 for (int j = 0; j < numPoints; j++) {
563 for (int k = 0; k < D4cardinality; k++) {
564
565 // basisD4 is (F,P,Dk), compute offset:
566 int l = k + j * D4cardinality + i * D4cardinality * numPoints;
567 if (std::abs(vals(i,j,k) - basisD4[l]) > INTREPID_TOL) {
568 errorFlag++;
569 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
570
571 // Output the multi-index of the value where the error is:
572 *outStream << " At multi-index { ";
573 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
574 *outStream << "} computed D4 component: " << vals(i,j,k)
575 << " but reference D4 component: " << basisD2[l] << "\n";
576 }
577 }
578 }
579 }
580
581
582
583 // Check D7 to D10 - must be zero. This basis does not cover D5 and D6
584 for(EOperator op = OPERATOR_D7; op < OPERATOR_MAX; op++) {
585
586 // The last dimension is the number of kth derivatives and needs to be resized for every Dk
587 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
588 vals.resize(numFields, numPoints, DkCardin);
589
590 hexBasis.getValues(vals, hexNodes, op);
591 for (int i = 0; i < vals.size(); i++) {
592 if (std::abs(vals[i]) > INTREPID_TOL) {
593 errorFlag++;
594 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
595
596 // Get the multi-index of the value where the error is and the operator order
597 std::vector<int> myIndex;
598 vals.getMultiIndex(myIndex,i);
599 int ord = Intrepid::getOperatorOrder(op);
600 *outStream << " At multi-index { ";
601 for(int j = 0; j < vals.rank(); j++) {
602 *outStream << myIndex[j] << " ";
603 }
604 *outStream << "} computed D"<< ord <<" component: " << vals[i]
605 << " but reference D" << ord << " component: 0 \n";
606 }
607 }
608 }
609 }
610 // Catch unexpected errors
611 catch (const std::logic_error & err) {
612 *outStream << err.what() << "\n\n";
613 errorFlag = -1000;
614 };
615
616 if (errorFlag != 0)
617 std::cout << "End Result: TEST FAILED\n";
618 else
619 std::cout << "End Result: TEST PASSED\n";
620
621 // reset format state of std::cout
622 std::cout.copyfmt(oldFormatState);
623
624 return errorFlag;
625}
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Header file for utility class to provide point tools, such as barycentric coordinates,...
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.
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...