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
49#include "Intrepid_HGRAD_WEDGE_C1_FEM.hpp"
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_C1_FEM) |\n" \
93 << "| |\n" \
94 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 << "| 2) Basis values for VALUE, GRAD, CURL, 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_C1_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 // Define array containing the 6 vertices of the reference TRIPRISM and some other points.
117 FieldContainer<double> wedgeNodes(12, 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.25; wedgeNodes(6,1) = 0.5; wedgeNodes(6,2) = -1.0;
126 wedgeNodes(7,0) = 0.5; wedgeNodes(7,1) = 0.25; wedgeNodes(7,2) = 0.0;
127 wedgeNodes(8,0) = 0.25; wedgeNodes(8,1) = 0.30; wedgeNodes(8,2) = 1.0;
128 wedgeNodes(9,0) = 0.3; wedgeNodes(9,1) = 0.0; wedgeNodes(9,2) = 0.75;
129 wedgeNodes(10,0)= 0.0; wedgeNodes(10,1)= 0.44; wedgeNodes(10,2)= -0.23;
130 wedgeNodes(11,0)= 0.4; wedgeNodes(11,1)= 0.6; wedgeNodes(11,2)= 0.0;
131
132
133 // Generic array for the output values; needs to be properly resized depending on the operator type
134 FieldContainer<double> vals;
135
136 try{
137 // exception #1: CURL cannot be applied to scalar functions
138 // resize vals to rank-3 container with dimensions (num. points, num. basis functions)
139 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 );
140 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
141
142 // exception #2: DIV cannot be applied to scalar functions
143 // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
144 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) );
145 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
146
147 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
148 // getDofTag() to access invalid array elements thereby causing bounds check exception
149 // exception #3
150 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException );
151 // exception #4
152 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException );
153 // exception #5
154 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,6,0), throwCounter, nException );
155 // exception #6
156 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(7), throwCounter, nException );
157 // exception #7
158 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException );
159
160#ifdef HAVE_INTREPID_DEBUG
161 // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays
162 // exception #8: input points array must be of rank-2
163 FieldContainer<double> badPoints1(4, 5, 3);
164 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
165
166 // exception #9 dimension 1 in the input point array must equal space dimension of the cell
167 FieldContainer<double> badPoints2(4, 2);
168 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
169
170 // exception #10 output values must be of rank-2 for OPERATOR_VALUE
171 FieldContainer<double> badVals1(4, 3, 1);
172 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
173
174 // exception #11 output values must be of rank-3 for OPERATOR_GRAD
175 FieldContainer<double> badVals2(4, 3);
176 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
177
178 // exception #12 output values must be of rank-3 for OPERATOR_D1
179 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D1), throwCounter, nException );
180
181 // exception #13 output values must be of rank-3 for OPERATOR_D2
182 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D2), throwCounter, nException );
183
184 // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
185 FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0));
186 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
187
188 // exception #15 incorrect 1st dimension of output array (must equal number of points)
189 FieldContainer<double> badVals4(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1);
190 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
191
192 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
193 FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 4);
194 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
195
196 // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D)
197 FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 40);
198 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D2), throwCounter, nException );
199
200 // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D)
201 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D3), throwCounter, nException );
202#endif
203
204 }
205 catch (const std::logic_error & err) {
206 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
207 *outStream << err.what() << '\n';
208 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
209 errorFlag = -1000;
210 };
211
212 // Check if number of thrown exceptions matches the one we expect - 18
213 if (throwCounter != nException) {
214 errorFlag++;
215 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
216 }
217
218 *outStream \
219 << "\n"
220 << "===============================================================================\n"\
221 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
222 << "===============================================================================\n";
223
224 try{
225 std::vector<std::vector<int> > allTags = wedgeBasis.getAllDofTags();
226
227 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
228 for (unsigned i = 0; i < allTags.size(); i++) {
229 int bfOrd = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
230
231 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
232 if( !( (myTag[0] == allTags[i][0]) &&
233 (myTag[1] == allTags[i][1]) &&
234 (myTag[2] == allTags[i][2]) &&
235 (myTag[3] == allTags[i][3]) ) ) {
236 errorFlag++;
237 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
238 *outStream << " getDofOrdinal( {"
239 << allTags[i][0] << ", "
240 << allTags[i][1] << ", "
241 << allTags[i][2] << ", "
242 << allTags[i][3] << "}) = " << bfOrd <<" but \n";
243 *outStream << " getDofTag(" << bfOrd << ") = { "
244 << myTag[0] << ", "
245 << myTag[1] << ", "
246 << myTag[2] << ", "
247 << myTag[3] << "}\n";
248 }
249 }
250
251 // Now do the same but loop over basis functions
252 for( int bfOrd = 0; bfOrd < wedgeBasis.getCardinality(); bfOrd++) {
253 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
254 int myBfOrd = wedgeBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
255 if( bfOrd != myBfOrd) {
256 errorFlag++;
257 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
258 *outStream << " getDofTag(" << bfOrd << ") = { "
259 << myTag[0] << ", "
260 << myTag[1] << ", "
261 << myTag[2] << ", "
262 << myTag[3] << "} but getDofOrdinal({"
263 << myTag[0] << ", "
264 << myTag[1] << ", "
265 << myTag[2] << ", "
266 << myTag[3] << "} ) = " << myBfOrd << "\n";
267 }
268 }
269 }
270 catch (const std::logic_error & err){
271 *outStream << err.what() << "\n\n";
272 errorFlag = -1000;
273 };
274
275 *outStream \
276 << "\n"
277 << "===============================================================================\n"\
278 << "| TEST 3: correctness of basis function values |\n"\
279 << "===============================================================================\n";
280
281 outStream -> precision(20);
282
283 // VALUE: Each row gives the 4 correct basis set values at an evaluation point
284 double basisValues[] = {
285 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
286 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
287 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
288 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
289 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
290 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
291 //
292 0.25, 0.25, 0.5, 0., 0., 0.,
293 0.125, 0.25, 0.125, 0.125, 0.25, 0.125,
294 0., 0., 0., 0.45, 0.25, 0.3,
295 0.0875, 0.0375, 0, 0.6125, 0.2625, 0,
296 0.3444, 0, 0.2706, 0.2156, 0, 0.1694,
297 0., 0.2, 0.3, 0., 0.2, 0.3
298 };
299
300 // GRAD and D1: each row gives the 3 x 4 correct values of the gradients of the 4 basis functions
301 double basisGrads[] = {
302 -1., -1., -0.5, 1., 0, 0, 0, 1., 0, 0., 0., 0.5, 0., 0, 0, 0, 0., 0, \
303 -1., -1., 0., 1., 0, -0.5, 0, 1., 0, 0., 0., 0., 0., 0, 0.5, 0, 0., \
304 0, -1., -1., 0., 1., 0, 0, 0, 1., -0.5, 0., 0., 0., 0., 0, 0, 0, 0., \
305 0.5, 0., 0., -0.5, 0., 0, 0, 0, 0., 0, -1., -1., 0.5, 1., 0, 0, 0, \
306 1., 0, 0., 0., 0., 0., 0, -0.5, 0, 0., 0, -1., -1., 0., 1., 0, 0.5, \
307 0, 1., 0, 0., 0., 0., 0., 0, 0, 0, 0., -0.5, -1., -1., 0., 1., 0, 0, \
308 0, 1., 0.5, -1., -1., -0.125, 1., 0, -0.125, 0, 1., -0.25, 0., 0., \
309 0.125, 0., 0, 0.125, 0, 0., 0.25, -0.5, -0.5, -0.125, 0.5, 0, -0.25, \
310 0, 0.5, -0.125, -0.5, -0.5, 0.125, 0.5, 0, 0.25, 0, 0.5, 0.125, 0., \
311 0., -0.225, 0., 0, -0.125, 0, 0., -0.15, -1., -1., 0.225, 1., 0, \
312 0.125, 0, 1., 0.15, -0.125, -0.125, -0.35, 0.125, 0, -0.15, 0, 0.125, \
313 0, -0.875, -0.875, 0.35, 0.875, 0, 0.15, 0, 0.875, 0, -0.615, -0.615, \
314 -0.28, 0.615, 0, 0, 0, 0.615, -0.22, -0.385, -0.385, 0.28, 0.385, 0, \
315 0, 0, 0.385, 0.22, -0.5, -0.5, 0., 0.5, 0, -0.2, 0, 0.5, -0.3, -0.5, \
316 -0.5, 0., 0.5, 0, 0.2, 0, 0.5, 0.3
317 };
318 //D2: flat array with the values of D2 applied to basis functions. Multi-index is (P,F,K)
319 double basisD2[] = {
320 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, \
321 -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, \
322 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, \
323 -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, \
324 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, \
325 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, \
326 -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, \
327 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, \
328 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, \
329 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, \
330 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, \
331 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, \
332 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, \
333 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, \
334 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, \
335 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, \
336 -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, \
337 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, \
338 -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, \
339 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, \
340 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, \
341 -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, \
342 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0
343 };
344 try{
345
346 // Dimensions for the output arrays:
347 int numFields = wedgeBasis.getCardinality();
348 int numPoints = wedgeNodes.dimension(0);
349 int spaceDim = wedgeBasis.getBaseCellTopology().getDimension();
350 int D2Cardin = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
351
352 // Generic array for values, grads, curls, etc. that will be properly sized before each call
353 FieldContainer<double> vals;
354
355 // Check VALUE of basis functions: resize vals to rank-2 container:
356 vals.resize(numFields, numPoints);
357 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE);
358 for (int i = 0; i < numFields; i++) {
359 for (int j = 0; j < numPoints; j++) {
360 int l = i + j * numFields;
361 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
362 errorFlag++;
363 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
364
365 // Output the multi-index of the value where the error is:
366 *outStream << " At multi-index { ";
367 *outStream << i << " ";*outStream << j << " ";
368 *outStream << "} computed value: " << vals(i,j)
369 << " but reference value: " << basisValues[l] << "\n";
370 }
371 }
372 }
373
374 // Check GRAD of basis function: resize vals to rank-3 container
375 vals.resize(numFields, numPoints, spaceDim);
376 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD);
377 for (int i = 0; i < numFields; i++) {
378 for (int j = 0; j < numPoints; j++) {
379 for (int k = 0; k < spaceDim; k++) {
380 int l = k + i * spaceDim + j * spaceDim * numFields;
381 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
382 errorFlag++;
383 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
384
385 // Output the multi-index of the value where the error is:
386 *outStream << " At multi-index { ";
387 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
388 *outStream << "} computed grad component: " << vals(i,j,k)
389 << " but reference grad component: " << basisGrads[l] << "\n";
390 }
391 }
392 }
393 }
394
395 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
396 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D1);
397 for (int i = 0; i < numFields; i++) {
398 for (int j = 0; j < numPoints; j++) {
399 for (int k = 0; k < spaceDim; k++) {
400 int l = k + i * spaceDim + j * spaceDim * numFields;
401 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
402 errorFlag++;
403 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
404
405 // Output the multi-index of the value where the error is:
406 *outStream << " At multi-index { ";
407 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
408 *outStream << "} computed D1 component: " << vals(i,j,k)
409 << " but reference D1 component: " << basisGrads[l] << "\n";
410 }
411 }
412 }
413 }
414
415 // Check D2 of basis function
416 vals.resize(numFields, numPoints, D2Cardin);
417 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D2);
418 for (int i = 0; i < numFields; i++) {
419 for (int j = 0; j < numPoints; j++) {
420 for (int k = 0; k < D2Cardin; k++) {
421 int l = k + i * D2Cardin + j * D2Cardin * numFields;
422 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
423 errorFlag++;
424 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
425
426 // Output the multi-index of the value where the error is:
427 *outStream << " At multi-index { ";
428 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
429 *outStream << "} computed D2 component: " << vals(i,j,k)
430 << " but reference D2 component: " << basisD2[l] << "\n";
431 }
432 }
433 }
434 }
435
436 // Check all higher derivatives - must be zero.
437 for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
438
439 // The last dimension is the number of kth derivatives and needs to be resized for every Dk
440 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
441 vals.resize(numFields, numPoints, DkCardin);
442
443 wedgeBasis.getValues(vals, wedgeNodes, op);
444 for (int i = 0; i < vals.size(); i++) {
445 if (std::abs(vals[i]) > INTREPID_TOL) {
446 errorFlag++;
447 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
448
449 // Get the multi-index of the value where the error is and the operator order
450 std::vector<int> myIndex;
451 vals.getMultiIndex(myIndex,i);
452 int ord = Intrepid::getOperatorOrder(op);
453 *outStream << " At multi-index { ";
454 for(int j = 0; j < vals.rank(); j++) {
455 *outStream << myIndex[j] << " ";
456 }
457 *outStream << "} computed D"<< ord <<" component: " << vals[i]
458 << " but reference D" << ord << " component: 0 \n";
459 }
460 }
461 }
462 }
463
464 // Catch unexpected errors
465 catch (const std::logic_error & err) {
466 *outStream << err.what() << "\n\n";
467 errorFlag = -1000;
468 };
469
470 if (errorFlag != 0)
471 std::cout << "End Result: TEST FAILED\n";
472 else
473 std::cout << "End Result: TEST PASSED\n";
474
475 // reset format state of std::cout
476 std::cout.copyfmt(oldFormatState);
477
478 return errorFlag;
479}
Header file for utility class to provide multidimensional containers.
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.