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