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
54#include "Teuchos_oblackholestream.hpp"
55#include "Teuchos_RCP.hpp"
56#include "Teuchos_GlobalMPISession.hpp"
57
58using namespace std;
59using namespace Intrepid;
60
61#define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
62{ \
63 ++nException; \
64 try { \
65 S ; \
66 } \
67 catch (const std::logic_error & err) { \
68 ++throwCounter; \
69 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
70 *outStream << err.what() << '\n'; \
71 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
72 }; \
73}
74
75
76int main(int argc, char *argv[]) {
77
78 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
79
80 // This little trick lets us print to std::cout only if
81 // a (dummy) command-line argument is provided.
82 int iprint = argc - 1;
83 Teuchos::RCP<std::ostream> outStream;
84 Teuchos::oblackholestream bhs; // outputs nothing
85 if (iprint > 0)
86 outStream = Teuchos::rcp(&std::cout, false);
87 else
88 outStream = Teuchos::rcp(&bhs, false);
89
90 // Save the format state of the original std::cout.
91 Teuchos::oblackholestream oldFormatState;
92 oldFormatState.copyfmt(std::cout);
93
94 *outStream \
95 << "===============================================================================\n" \
96 << "| |\n" \
97 << "| Unit Test (Basis_HGRAD_LINE_Cn_FEM_JACOBI) |\n" \
98 << "| |\n" \
99 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
100 << "| 2) Basis values for VALUE, GRAD, CURL, and Dk operators |\n" \
101 << "| |\n" \
102 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
103 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
104 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
105 << "| |\n" \
106 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
107 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
108 << "| |\n" \
109 << "===============================================================================\n"\
110 << "| TEST 1: Basis creation, exception testing |\n"\
111 << "===============================================================================\n";
112
113
114 // Define basis and error flag
115 double alpha = 0.0, beta = 0.0;
116 Basis_HGRAD_LINE_Cn_FEM_JACOBI<double, FieldContainer<double> > lineBasis(5, alpha, beta);
117 int errorFlag = 0;
118
119 // Initialize throw counter for exception testing
120 int nException = 0;
121 int throwCounter = 0;
122
123 // Define array containing vertices of the reference Line and a few other points
124 int numIntervals = 100;
125 FieldContainer<double> lineNodes(numIntervals+1, 1);
126 for (int i=0; i<numIntervals+1; i++) {
127 lineNodes(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
128 }
129
130 // Generic array for the output values; needs to be properly resized depending on the operator type
131 FieldContainer<double> vals;
132
133 try{
134 // Exceptions 1-5: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
135 // getDofTag() to access invalid array elements thereby causing bounds check exception
136 // exception #1
137 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(2,0,0), throwCounter, nException );
138 // exception #2
139 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,1,1), throwCounter, nException );
140 // exception #3
141 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,0,7), throwCounter, nException );
142 // not an exception
143 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,0,5), throwCounter, nException ); --nException;
144 // exception #4
145 INTREPID_TEST_COMMAND( lineBasis.getDofTag(6), throwCounter, nException );
146 // exception #5
147 INTREPID_TEST_COMMAND( lineBasis.getDofTag(-1), throwCounter, nException );
148 // not an exception
149 INTREPID_TEST_COMMAND( lineBasis.getDofTag(5), throwCounter, nException ); --nException;
150#ifdef HAVE_INTREPID_DEBUG
151 // Exceptions 6-16 test exception handling with incorrectly dimensioned input/output arrays
152 // exception #6: input points array must be of rank-2
153 FieldContainer<double> badPoints1(4, 5, 3);
154 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
155
156 // exception #7: dimension 1 in the input point array must equal space dimension of the cell
157 FieldContainer<double> badPoints2(4, 3);
158 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
159
160 // exception #8: output values must be of rank-2 for OPERATOR_VALUE
161 FieldContainer<double> badVals1(4, 3, 1);
162 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals1, lineNodes, OPERATOR_VALUE), throwCounter, nException );
163
164 // exception #9: output values must be of rank-3 for OPERATOR_GRAD
165 FieldContainer<double> badVals2(4, 3);
166 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_GRAD), throwCounter, nException );
167
168 // exception #10: output values must be of rank-3 for OPERATOR_CURL
169 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_CURL), throwCounter, nException );
170
171 // exception #11: output values must be of rank-2 for OPERATOR_DIV
172 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_DIV), throwCounter, nException );
173
174 // exception #12: output values must be of rank-2 for OPERATOR_D1
175 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_D1), throwCounter, nException );
176
177 // exception #13: incorrect 0th dimension of output array (must equal number of basis functions)
178 FieldContainer<double> badVals3(lineBasis.getCardinality() + 1, lineNodes.dimension(0));
179 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals3, lineNodes, OPERATOR_VALUE), throwCounter, nException );
180
181 // exception #14: incorrect 1st dimension of output array (must equal number of points)
182 FieldContainer<double> badVals4(lineBasis.getCardinality(), lineNodes.dimension(0) + 1);
183 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals4, lineNodes, OPERATOR_VALUE), throwCounter, nException );
184
185 // exception #15: incorrect 2nd dimension of output array (must equal spatial dimension)
186 FieldContainer<double> badVals5(lineBasis.getCardinality(), lineNodes.dimension(0), 2);
187 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals5, lineNodes, OPERATOR_GRAD), throwCounter, nException );
188
189 // not an exception
190 FieldContainer<double> goodVals2(lineBasis.getCardinality(), lineNodes.dimension(0));
191 INTREPID_TEST_COMMAND( lineBasis.getValues(goodVals2, lineNodes, OPERATOR_VALUE), throwCounter, nException ); --nException;
192#endif
193
194 }
195 catch (const std::logic_error & err) {
196 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
197 *outStream << err.what() << '\n';
198 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
199 errorFlag = -1000;
200 };
201
202 // Check if number of thrown exceptions matches the one we expect
203 if (throwCounter != nException) {
204 errorFlag++;
205 *outStream << std::setw(70) << "FAILURE! Incorrect number of exceptions." << "\n";
206 }
207
208
209 *outStream \
210 << "\n"
211 << "===============================================================================\n"\
212 << "| TEST 3: orthogonality of basis functions |\n"\
213 << "===============================================================================\n";
214
215 outStream -> precision(20);
216
217 try {
218
219 // Check orthogonality property for Legendre polynomials.
220 int maxorder = 10;
221
222 DefaultCubatureFactory<double> cubFactory; // create factory
223 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >()); // create cell topology
224 for (int ordi=0; ordi < maxorder; ordi++) {
225 //create left basis
226 Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasisLeft =
227 Teuchos::rcp(new Basis_HGRAD_LINE_Cn_FEM_JACOBI<double,FieldContainer<double> >(ordi) );
228
229 for (int ordj=0; ordj < maxorder; ordj++) {
230
231 //create right basis
232 Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasisRight =
233 Teuchos::rcp(new Basis_HGRAD_LINE_Cn_FEM_JACOBI<double,FieldContainer<double> >(ordj) );
234
235 // get cubature points and weights
236 Teuchos::RCP<Cubature<double> > lineCub = cubFactory.create(line, ordi+ordj);
237 int numPoints = lineCub->getNumPoints();
238 FieldContainer<double> cubPoints (numPoints, lineCub->getDimension());
239 FieldContainer<double> cubWeights(numPoints);
240 FieldContainer<double> cubWeightsC(1, numPoints);
241 lineCub->getCubature(cubPoints, cubWeights);
242 // "reshape" weights
243 for (int i=0; i<numPoints; i++) { cubWeightsC(0,i) = cubWeights(i); }
244
245
246 // get basis values
247 int numFieldsLeft = lineBasisLeft ->getCardinality();
248 int numFieldsRight = lineBasisRight->getCardinality();
249 FieldContainer<double> valsLeft(numFieldsLeft,numPoints),
250 valsRight(numFieldsRight,numPoints);
251 lineBasisLeft ->getValues(valsLeft, cubPoints, OPERATOR_VALUE);
252 lineBasisRight->getValues(valsRight, cubPoints, OPERATOR_VALUE);
253
254 // reshape by cloning and integrate
255 FieldContainer<double> valsLeftC(1, numFieldsLeft,numPoints),
256 valsRightC(1, numFieldsRight,numPoints),
257 massMatrix(1, numFieldsLeft, numFieldsRight);
258 ArrayTools::cloneFields<double>(valsLeftC, valsLeft);
259 ArrayTools::cloneFields<double>(valsRightC, valsRight);
260 ArrayTools::scalarMultiplyDataField<double>(valsRightC, cubWeightsC, valsRightC);
261 FunctionSpaceTools::integrate<double>(massMatrix, valsLeftC, valsRightC, COMP_CPP);
262
263 // check orthogonality property
264 for (int i=0; i<numFieldsLeft; i++) {
265 for (int j=0; j<numFieldsRight; j++) {
266
267 if (i==j) {
268 if ( std::abs(massMatrix(0,i,j)-(double)(2.0/(2.0*j+1.0))) > INTREPID_TOL ) {
269 *outStream << "Incorrect ii (\"diagonal\") value for i=" << i << ", j=" << j << ": "
270 << massMatrix(0,i,j) << " != " << "2/(2*" << j << "+1)\n\n";
271 errorFlag++;
272 }
273 }
274 else {
275 if ( std::abs(massMatrix(0,i,j)) > INTREPID_TOL ) {
276 *outStream << "Incorrect ij (\"off-diagonal\") value for i=" << i << ", j=" << j << ": "
277 << massMatrix(0,i,j) << " != " << "0\n\n";
278 errorFlag++;
279 }
280 }
281 }
282 }
283
284 }
285 }
286
287 }
288 // Catch unexpected errors
289 catch (const std::logic_error & err) {
290 *outStream << err.what() << "\n\n";
291 errorFlag = -1000;
292 };
293
294 *outStream \
295 << "\n"
296 << "===============================================================================\n"\
297 << "| TEST 4: correctness of basis function derivatives |\n"\
298 << "===============================================================================\n";
299
300 outStream -> precision(20);
301
302 // function values stored by bf, then pt
303 double basisValues[] = {
304 1.000000000000000, 1.000000000000000, 1.000000000000000, \
305 1.000000000000000, -1.000000000000000, -0.3333333333333333, \
306 0.3333333333333333, 1.000000000000000, 1.000000000000000, \
307 -0.3333333333333333, -0.3333333333333333, 1.000000000000000, \
308 -1.000000000000000, 0.4074074074074074, -0.4074074074074074, \
309 1.000000000000000};
310
311 double basisD1Values[] =
312 {0, 0, 0, 0, 1.000000000000000, 1.000000000000000, 1.000000000000000, \
313 1.000000000000000, -3.000000000000000, -1.000000000000000, \
314 1.000000000000000, 3.000000000000000, 6.000000000000000, \
315 -0.6666666666666667, -0.6666666666666667, 6.000000000000000};
316
317 double basisD2Values[] =
318 {0, 0, 0, 0, 0, 0, 0, 0, 3.000000000000000, 3.000000000000000, \
319 3.000000000000000, 3.000000000000000, -15.00000000000000, \
320 -5.000000000000000, 5.000000000000000, 15.00000000000000};
321
322 double basisD3Values[] =
323 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.00000000000000, \
324 15.00000000000000, 15.00000000000000, 15.00000000000000};
325
326
327
328 try {
329 Basis_HGRAD_LINE_Cn_FEM_JACOBI<double, FieldContainer<double> > lineBasis3(3, alpha, beta);
330 int numIntervals = 3;
331 FieldContainer<double> lineNodes3(numIntervals+1, 1);
332 FieldContainer<double> vals;
333 for (int i=0; i<numIntervals+1; i++) {
334 lineNodes3(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
335 }
336 int numFields = lineBasis3.getCardinality();
337 int numPoints = lineNodes3.dimension(0);
338
339 // test basis values
340 vals.resize(numFields, numPoints);
341 lineBasis3.getValues(vals,lineNodes3,OPERATOR_VALUE);
342 for (int i = 0; i < numFields; i++) {
343 for (int j = 0; j < numPoints; j++) {
344
345 // Compute offset for (F,P) container
346 int l = j + i * numPoints;
347 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
348 errorFlag++;
349 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
350
351 // Output the multi-index of the value where the error is:
352 *outStream << " At multi-index { ";
353 *outStream << i << " ";*outStream << j << " ";
354 *outStream << "} computed value: " << vals(i,j)
355 << " but reference value: " << basisValues[l] << "\n";
356 }
357 }
358 }
359
360 // test basis derivatives
361 vals.resize(numFields, numPoints,1);
362 lineBasis3.getValues(vals,lineNodes3,OPERATOR_D1);
363 for (int i = 0; i < numFields; i++) {
364 for (int j = 0; j < numPoints; j++) {
365
366 // Compute offset for (F,P) container
367 int l = j + i * numPoints;
368 if (std::abs(vals(i,j,0) - basisD1Values[l]) > INTREPID_TOL) {
369 errorFlag++;
370 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
371
372 // Output the multi-index of the value where the error is:
373 *outStream << " At multi-index { ";
374 *outStream << i << " ";*outStream << j << " ";
375 *outStream << "} computed value: " << vals(i,j,0)
376 << " but reference value: " << basisD1Values[l] << "\n";
377 }
378 }
379 }
380
381 vals.resize(numFields, numPoints,1);
382 lineBasis3.getValues(vals,lineNodes3,OPERATOR_D2);
383 for (int i = 0; i < numFields; i++) {
384 for (int j = 0; j < numPoints; j++) {
385
386 // Compute offset for (F,P) container
387 int l = j + i * numPoints;
388 if (std::abs(vals(i,j,0) - basisD2Values[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,0)
396 << " but reference value: " << basisD2Values[l] << "\n";
397 }
398 }
399 }
400
401 vals.resize(numFields, numPoints,1);
402 lineBasis3.getValues(vals,lineNodes3,OPERATOR_D3);
403 for (int i = 0; i < numFields; i++) {
404 for (int j = 0; j < numPoints; j++) {
405
406 // Compute offset for (F,P) container
407 int l = j + i * numPoints;
408 if (std::abs(vals(i,j,0) - basisD3Values[l]) > INTREPID_TOL) {
409 errorFlag++;
410 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
411
412 // Output the multi-index of the value where the error is:
413 *outStream << " At multi-index { ";
414 *outStream << i << " ";*outStream << j << " ";
415 *outStream << "} computed value: " << vals(i,j,0)
416 << " but reference value: " << basisD3Values[l] << "\n";
417 }
418 }
419 }
420 }
421 // Catch unexpected errors
422 catch (const std::logic_error & err) {
423 *outStream << err.what() << "\n\n";
424 errorFlag = -1000;
425 };
426
427
428 if (errorFlag != 0)
429 std::cout << "End Result: TEST FAILED\n";
430 else
431 std::cout << "End Result: TEST PASSED\n";
432
433 // reset format state of std::cout
434 std::cout.copyfmt(oldFormatState);
435
436 return errorFlag;
437}
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for the Intrepid::HGRAD_LINE_Cn_FEM class.
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1,...