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
44
51#include "Teuchos_oblackholestream.hpp"
52#include "Teuchos_RCP.hpp"
53#include "Teuchos_GlobalMPISession.hpp"
56#include "Shards_CellTopology.hpp"
57
58#include <iostream>
59using namespace Intrepid;
60
65int main(int argc, char *argv[]) {
66
67 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
68
69 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
70 int iprint = argc - 1;
71
72 Teuchos::RCP<std::ostream> outStream;
73 Teuchos::oblackholestream bhs; // outputs nothing
74
75 if (iprint > 0)
76 outStream = Teuchos::rcp(&std::cout, false);
77 else
78 outStream = Teuchos::rcp(&bhs, false);
79
80 // Save the format state of the original std::cout.
81 Teuchos::oblackholestream oldFormatState;
82 oldFormatState.copyfmt(std::cout);
83
84 *outStream \
85 << "===============================================================================\n" \
86 << "| |\n" \
87 << "| Unit Test HCURL_TET_In_FEM |\n" \
88 << "| |\n" \
89 << "| 1) Tests tetrahedral Nedelec basis |\n" \
90 << "| |\n" \
91 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
92 << "| Denis Ridzal (dridzal@sandia.gov) or |\n" \
93 << "| Robert Kirby (robert.c.kirby@ttu.edu) |\n" \
94 << "| |\n" \
95 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
96 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
97 << "| |\n" \
98 << "===============================================================================\n";
99
100 int errorFlag = 0;
101
102 // code-code comparison with FIAT
103 try {
104 const int deg = 1;
105 Basis_HCURL_TET_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED );
106
107 const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
108 FieldContainer<double> lattice( np_lattice , 3 );
109 FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice , 3 );
110 PointTools::getLattice<double,FieldContainer<double> >( lattice ,
111 myBasis.getBaseCellTopology() ,
112 deg ,
113 0 ,
114 POINTTYPE_EQUISPACED );
115
116 myBasis.getValues( myBasisValues , lattice , OPERATOR_VALUE );
117
118 const double fiat_vals[] = {
1191.000000000000001e+00, -2.498001805406602e-16, -1.665334536937735e-16,
1209.999999999999998e-01, 1.000000000000000e+00, 1.000000000000000e+00,
1215.828670879282072e-16, 1.110223024625157e-16, 2.498001805406602e-16,
1227.771561172376096e-16, 8.326672684688674e-17, 1.110223024625157e-16,
1232.081668171172169e-16, -2.914335439641036e-16, 1.280865063236792e-16,
124-3.191891195797325e-16, 1.000000000000000e+00, -4.293998586504916e-17,
125-9.999999999999994e-01, 2.081668171172169e-16, 2.400576428367544e-16,
1262.220446049250313e-16, -5.551115123125783e-17, 1.084013877651281e-16,
1273.469446951953614e-16, -1.000000000000000e+00, 1.387778780781446e-16,
128-1.804112415015879e-16, 1.942890293094024e-16, -1.387778780781446e-16,
129-9.999999999999993e-01, -9.999999999999996e-01, -9.999999999999998e-01,
1305.551115123125783e-17, -2.220446049250313e-16, -8.326672684688674e-17,
131-2.220446049250313e-16, -5.551115123125783e-17, 9.999999999999999e-01,
1321.665334536937735e-16, 1.110223024625157e-16, -6.383782391594650e-16,
1331.110223024625157e-16, 1.110223024625157e-16, -1.110223024625157e-16,
1349.999999999999990e-01, 9.999999999999994e-01, 9.999999999999996e-01,
1351.387778780781446e-16, -2.496931404305374e-17, -1.665334536937735e-16,
136-2.498001805406602e-16, -2.149987498083074e-16, 1.000000000000000e+00,
1378.326672684688674e-17, -3.769887250591415e-17, 8.326672684688674e-17,
138-9.999999999999994e-01, 1.556977698723022e-16, 2.220446049250313e-16,
139-9.422703950001342e-18, 1.665334536937735e-16, -2.359223927328458e-16,
140-9.422703950001268e-18, -8.326672684688674e-17, 1.387778780781446e-17,
141-7.525083148581445e-17, 2.775557561562891e-17, 1.000000000000000e+00,
1422.789513560273035e-16, -9.999999999999998e-01, -5.551115123125783e-17
143 };
144
145 int cur=0;
146 for (int i=0;i<myBasisValues.dimension(0);i++) {
147 for (int j=0;j<myBasisValues.dimension(1);j++) {
148 for (int k=0;k<myBasisValues.dimension(2);k++) {
149 if (std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) > 10.0*INTREPID_TOL ) {
150 errorFlag++;
151 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
152
153 // Output the multi-index of the value where the error is:
154 *outStream << " At multi-index { ";
155 *outStream << i << " " << j << " " << k;
156 *outStream << "} computed value: " << myBasisValues(i,j,k)
157 << " but correct value: " << fiat_vals[cur] << "\n";
158 *outStream << "Difference: " << std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) << "\n";
159 }
160 cur++;
161 }
162 }
163 }
164 }
165 catch (const std::exception & err) {
166 *outStream << err.what() << "\n\n";
167 errorFlag = -1000;
168 }
169
170 try {
171 const int deg = 1;
172 Basis_HCURL_TET_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED );
173
174 const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
175 FieldContainer<double> lattice( np_lattice , 3 );
176 FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice , 3 );
177 PointTools::getLattice<double,FieldContainer<double> >( lattice ,
178 myBasis.getBaseCellTopology() ,
179 deg ,
180 0 ,
181 POINTTYPE_EQUISPACED );
182
183 myBasis.getValues( myBasisValues , lattice , OPERATOR_CURL );
184
185 const double fiat_curls[] = {
186 -5.551115123125783e-16, -2.000000000000000e+00, 2.000000000000000e+00,
187 -5.551115123125783e-16, -2.000000000000000e+00, 2.000000000000000e+00,
188 -5.551115123125783e-16, -2.000000000000000e+00, 2.000000000000000e+00,
189 -5.551115123125783e-16, -2.000000000000000e+00, 2.000000000000000e+00,
190 -4.440892098500626e-16, -2.775557561562891e-16, 2.000000000000000e+00,
191 -4.440892098500626e-16, -2.775557561562891e-16, 2.000000000000000e+00,
192 -4.440892098500626e-16, -2.775557561562891e-16, 2.000000000000000e+00,
193 -4.440892098500626e-16, -2.775557561562891e-16, 2.000000000000000e+00,
194 -2.000000000000000e+00, 5.551115123125783e-17, 2.000000000000000e+00,
195 -2.000000000000000e+00, 5.551115123125783e-17, 2.000000000000000e+00,
196 -2.000000000000000e+00, 5.551115123125783e-17, 2.000000000000000e+00,
197 -2.000000000000000e+00, 5.551115123125783e-17, 2.000000000000000e+00,
198 -2.000000000000000e+00, 2.000000000000000e+00, 9.861075762086680e-17,
199 -2.000000000000000e+00, 2.000000000000000e+00, 9.861075762086680e-17,
200 -2.000000000000000e+00, 2.000000000000000e+00, 9.861075762086680e-17,
201 -2.000000000000000e+00, 2.000000000000000e+00, 9.861075762086680e-17,
202 -2.775557561562891e-17, -2.000000000000000e+00, 4.287451790760826e-16,
203 -2.775557561562891e-17, -2.000000000000000e+00, 4.287451790760826e-16,
204 -2.775557561562891e-17, -2.000000000000000e+00, 4.287451790760826e-16,
205 -2.775557561562891e-17, -2.000000000000000e+00, 4.287451790760826e-16,
206 2.000000000000000e+00, -2.185751579730777e-16, 1.526556658859590e-16,
207 2.000000000000000e+00, -2.185751579730777e-16, 1.526556658859590e-16,
208 2.000000000000000e+00, -2.185751579730777e-16, 1.526556658859590e-16,
209 2.000000000000000e+00, -2.185751579730777e-16, 1.526556658859590e-16
210 };
211
212 int cur=0;
213 for (int i=0;i<myBasisValues.dimension(0);i++) {
214 for (int j=0;j<myBasisValues.dimension(1);j++) {
215 for (int k=0;k<myBasisValues.dimension(2);k++) {
216 if (std::abs( myBasisValues(i,j,k) - fiat_curls[cur] ) > 10.0*INTREPID_TOL ) {
217 errorFlag++;
218 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
219
220 // Output the multi-index of the value where the error is:
221 *outStream << " At multi-index { ";
222 *outStream << i << " " << j << " " << k;
223 *outStream << "} computed value: " << myBasisValues(i,j,k)
224 << " but correct value: " << fiat_curls[cur] << "\n";
225 *outStream << "Difference: " << std::abs( myBasisValues(i,j,k) - fiat_curls[cur] ) << "\n";
226 }
227 cur++;
228 }
229 }
230 }
231 }
232 catch (const std::exception & err) {
233 *outStream << err.what() << "\n\n";
234 errorFlag = -1000;
235 }
236
237
238 if (errorFlag != 0)
239 std::cout << "End Result: TEST FAILED\n";
240 else
241 std::cout << "End Result: TEST PASSED\n";
242
243 // reset format state of std::cout
244 std::cout.copyfmt(oldFormatState);
245
246 return errorFlag;
247}
int main(int argc, char *argv[])
Performs a code-code comparison to FIAT for Nedelec bases on tets (values and curls)
Definition test_01.cpp:65
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HCURL_TET_In_FEM class.
Header file for utility class to provide point tools, such as barycentric coordinates,...
static int getLatticeSize(const shards::CellTopology &cellType, const int order, const int offset=0)
Computes the number of points in a lattice of a given order on a simplex (currently disabled for othe...