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"
57#include "Shards_CellTopology.hpp"
58
59#include <iostream>
60using namespace Intrepid;
61
66int main(int argc, char *argv[]) {
67
68 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
69
70 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
71 int iprint = argc - 1;
72
73 Teuchos::RCP<std::ostream> outStream;
74 Teuchos::oblackholestream bhs; // outputs nothing
75
76 if (iprint > 0)
77 outStream = Teuchos::rcp(&std::cout, false);
78 else
79 outStream = Teuchos::rcp(&bhs, false);
80
81 // Save the format state of the original std::cout.
82 Teuchos::oblackholestream oldFormatState;
83 oldFormatState.copyfmt(std::cout);
84
85 *outStream \
86 << "===============================================================================\n" \
87 << "| |\n" \
88 << "| Unit Test HCURL_TRI_In_FEM |\n" \
89 << "| |\n" \
90 << "| 1) Tests triangular Nedelec-Thomas basis |\n" \
91 << "| |\n" \
92 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
93 << "| Denis Ridzal (dridzal@sandia.gov) or |\n" \
94 << "| Robert Kirby (robert.c.kirby@ttu.edu) |\n" \
95 << "| |\n" \
96 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
97 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
98 << "| |\n" \
99 << "===============================================================================\n";
100
101 int errorFlag = 0;
102
103 // test for basis values, compare against H(div) basis, for they should be related by rotation
104 // in the lowest order case
105 try {
106 const int deg = 1;
107 Basis_HDIV_TRI_In_FEM<double,FieldContainer<double> > myBasisDIV( deg , POINTTYPE_EQUISPACED );
108 Basis_HCURL_TRI_In_FEM<double,FieldContainer<double> > myBasisCURL( deg , POINTTYPE_EQUISPACED );
109
110 // Get a lattice
111 const int np_lattice = PointTools::getLatticeSize( myBasisDIV.getBaseCellTopology() , deg , 0 );
112 FieldContainer<double> lattice( np_lattice , 2 );
113
114 FieldContainer<double> myBasisValuesDIV( myBasisDIV.getCardinality() , np_lattice , 2 );
115 FieldContainer<double> myBasisValuesCURL( myBasisDIV.getCardinality() , np_lattice , 2 );
116 PointTools::getLattice<double,FieldContainer<double> >( lattice ,
117 myBasisDIV.getBaseCellTopology() ,
118 deg ,
119 0 ,
120 POINTTYPE_EQUISPACED );
121
122 myBasisDIV.getValues( myBasisValuesDIV , lattice , OPERATOR_VALUE );
123 myBasisCURL.getValues( myBasisValuesCURL, lattice , OPERATOR_VALUE );
124
125 for (int i=0;i<myBasisValuesDIV.dimension(0);i++) {
126 for (int j=0;j<myBasisValuesDIV.dimension(1);j++) {
127 if (std::abs( myBasisValuesDIV(i,j,1) + myBasisValuesCURL(i,j,0) ) > INTREPID_TOL ) {
128 errorFlag++;
129 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
130 // Output the multi-index of the value where the error is:
131 *outStream << " At multi-index { ";
132 *outStream << i << " " << j << " and component 0";
133 *outStream << "} computed value: " << myBasisValuesCURL(i,j,0)
134 << " but correct value: " << -myBasisValuesDIV(i,j,1) << "\n";
135 *outStream << "Difference: " << std::abs( myBasisValuesCURL(i,j,0) + myBasisValuesDIV(i,j,1) ) << "\n";
136 }
137 if (std::abs( myBasisValuesDIV(i,j,0) - myBasisValuesCURL(i,j,1) ) > INTREPID_TOL ) {
138 errorFlag++;
139 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
140 // Output the multi-index of the value where the error is:
141 *outStream << " At multi-index { ";
142 *outStream << i << " " << j << " and component 1";
143 *outStream << "} computed value: " << myBasisValuesCURL(i,j,1)
144 << " but correct value: " << myBasisValuesDIV(i,j,0) << "\n";
145 *outStream << "Difference: " << std::abs( myBasisValuesCURL(i,j,1) - myBasisValuesDIV(i,j,1) ) << "\n";
146 }
147 }
148 }
149 }
150 catch (const std::exception & err) {
151 *outStream << err.what() << "\n\n";
152 errorFlag = -1000;
153 }
154
155 // next test: code-code comparison with FIAT
156 try {
157 const int deg = 2;
158 Basis_HCURL_TRI_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED );
159 const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
160 FieldContainer<double> lattice( np_lattice , 2 );
161 PointTools::getLattice<double,FieldContainer<double> >( lattice ,
162 myBasis.getBaseCellTopology() ,
163 deg ,
164 0 ,
165 POINTTYPE_EQUISPACED );
166 FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice , 2 );
167
168
169 myBasis.getValues( myBasisValues, lattice , OPERATOR_VALUE );
170
171 const double fiat_vals[] = {
1722.000000000000000e+00, 0.000000000000000e+00,
1735.000000000000000e-01, 2.500000000000000e-01,
174-1.000000000000000e+00, -1.000000000000000e+00,
1752.500000000000000e-01, 0.000000000000000e+00,
176-5.000000000000000e-01, -5.000000000000000e-01,
1770.000000000000000e+00, 0.000000000000000e+00,
178-1.000000000000000e+00, 0.000000000000000e+00,
1795.000000000000000e-01, 2.500000000000000e-01,
1802.000000000000000e+00, 2.000000000000000e+00,
181-5.000000000000000e-01, 0.000000000000000e+00,
1822.500000000000000e-01, 2.500000000000000e-01,
1830.000000000000000e+00, 0.000000000000000e+00,
1840.000000000000000e+00, 0.000000000000000e+00,
1850.000000000000000e+00, 2.500000000000000e-01,
1860.000000000000000e+00, 2.000000000000000e+00,
1875.000000000000000e-01, 0.000000000000000e+00,
188-2.500000000000000e-01, 2.500000000000000e-01,
1891.000000000000000e+00, 0.000000000000000e+00,
1900.000000000000000e+00, 0.000000000000000e+00,
1910.000000000000000e+00, -5.000000000000000e-01,
1920.000000000000000e+00, -1.000000000000000e+00,
193-2.500000000000000e-01, 0.000000000000000e+00,
194-2.500000000000000e-01, 2.500000000000000e-01,
195-2.000000000000000e+00, 0.000000000000000e+00,
1960.000000000000000e+00, 1.000000000000000e+00,
1970.000000000000000e+00, 5.000000000000000e-01,
1980.000000000000000e+00, 0.000000000000000e+00,
199-2.500000000000000e-01, -5.000000000000000e-01,
200-2.500000000000000e-01, -2.500000000000000e-01,
201-2.000000000000000e+00, -2.000000000000000e+00,
2020.000000000000000e+00, -2.000000000000000e+00,
2030.000000000000000e+00, -2.500000000000000e-01,
2040.000000000000000e+00, 0.000000000000000e+00,
205-2.500000000000000e-01, -5.000000000000000e-01,
2065.000000000000000e-01, 5.000000000000000e-01,
2071.000000000000000e+00, 1.000000000000000e+00,
2080.000000000000000e+00, 0.000000000000000e+00,
2090.000000000000000e+00, -7.500000000000000e-01,
2100.000000000000000e+00, 0.000000000000000e+00,
2111.500000000000000e+00, 0.000000000000000e+00,
2127.500000000000000e-01, 7.500000000000000e-01,
2130.000000000000000e+00, 0.000000000000000e+00,
2140.000000000000000e+00, 0.000000000000000e+00,
2150.000000000000000e+00, 1.500000000000000e+00,
2160.000000000000000e+00, 0.000000000000000e+00,
217-7.500000000000000e-01, 0.000000000000000e+00,
2187.500000000000000e-01, 7.500000000000000e-01,
2190.000000000000000e+00, 0.000000000000000e+00
220 };
221
222 int cur=0;
223 for (int i=0;i<myBasisValues.dimension(0);i++) {
224 for (int j=0;j<myBasisValues.dimension(1);j++) {
225 for (int k=0;k<myBasisValues.dimension(2);k++) {
226 if (std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) > INTREPID_TOL ) {
227 errorFlag++;
228 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
229
230 // Output the multi-index of the value where the error is:
231 *outStream << " At multi-index { ";
232 *outStream << i << " " << j << " " << k;
233 *outStream << "} computed value: " << myBasisValues(i,j,k)
234 << " but correct value: " << fiat_vals[cur] << "\n";
235 *outStream << "Difference: " << std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) << "\n";
236 }
237 cur++;
238 }
239 }
240 }
241 }
242 catch (const std::exception & err) {
243 *outStream << err.what() << "\n\n";
244 errorFlag = -1000;
245 }
246
247 try {
248 const int deg = 2;
249 Basis_HCURL_TRI_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED );
250 const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
251 FieldContainer<double> lattice( np_lattice , 2 );
252 PointTools::getLattice<double,FieldContainer<double> >( lattice ,
253 myBasis.getBaseCellTopology() ,
254 deg ,
255 0 ,
256 POINTTYPE_EQUISPACED );
257 FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice );
258
259
260 myBasis.getValues( myBasisValues, lattice , OPERATOR_CURL );
261
262 const double fiat_curls[] = {
2637.000000000000000e+00,
2642.500000000000000e+00,
265-2.000000000000000e+00,
2662.500000000000000e+00,
267-2.000000000000000e+00,
268-2.000000000000000e+00,
269-2.000000000000000e+00,
2702.500000000000000e+00,
2717.000000000000000e+00,
272-2.000000000000000e+00,
2732.500000000000000e+00,
274-2.000000000000000e+00,
275-2.000000000000000e+00,
2762.500000000000000e+00,
2777.000000000000000e+00,
278-2.000000000000000e+00,
2792.500000000000000e+00,
280-2.000000000000000e+00,
281-2.000000000000000e+00,
282-2.000000000000000e+00,
283-2.000000000000000e+00,
2842.500000000000000e+00,
2852.500000000000000e+00,
2867.000000000000000e+00,
287-2.000000000000000e+00,
288-2.000000000000000e+00,
289-2.000000000000000e+00,
2902.500000000000000e+00,
2912.500000000000000e+00,
2927.000000000000000e+00,
2937.000000000000000e+00,
2942.500000000000000e+00,
295-2.000000000000000e+00,
2962.500000000000000e+00,
297-2.000000000000000e+00,
298-2.000000000000000e+00,
299-9.000000000000000e+00,
300-4.500000000000000e+00,
3010.000000000000000e+00,
3020.000000000000000e+00,
3034.500000000000000e+00,
3049.000000000000000e+00,
3059.000000000000000e+00,
3060.000000000000000e+00,
307-9.000000000000000e+00,
3084.500000000000000e+00,
309-4.500000000000000e+00,
3100.000000000000000e+00
311 };
312
313 int cur=0;
314 for (int i=0;i<myBasisValues.dimension(0);i++) {
315 for (int j=0;j<myBasisValues.dimension(1);j++) {
316 if (std::abs( myBasisValues(i,j) - fiat_curls[cur] ) > INTREPID_TOL ) {
317 errorFlag++;
318 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
319
320 // Output the multi-index of the value where the error is:
321 *outStream << " At multi-index { ";
322 *outStream << i << " " << j;
323 *outStream << "} computed value: " << myBasisValues(i,j)
324 << " but correct value: " << fiat_curls[cur] << "\n";
325 *outStream << "Difference: " << std::abs( myBasisValues(i,j) - fiat_curls[cur] ) << "\n";
326 }
327 cur++;
328 }
329 }
330 }
331 catch (const std::exception & err) {
332 *outStream << err.what() << "\n\n";
333 errorFlag = -1000;
334 }
335
336 if (errorFlag != 0)
337 std::cout << "End Result: TEST FAILED\n";
338 else
339 std::cout << "End Result: TEST PASSED\n";
340
341 // reset format state of std::cout
342 std::cout.copyfmt(oldFormatState);
343
344 return errorFlag;
345}
int main(int argc, char *argv[])
Performs a code-code comparison to FIAT for Nedelec bases on triangles (values and curls)
Definition test_01.cpp:66
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HCURL_TRI_In_FEM class.
Header file for the Intrepid::HDIV_TRI_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...