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