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