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 
50 #include "Intrepid_FieldContainer.hpp"
51 #include "Teuchos_oblackholestream.hpp"
52 #include "Teuchos_RCP.hpp"
53 #include "Teuchos_GlobalMPISession.hpp"
54 #include "Intrepid_PointTools.hpp"
55 #include "Intrepid_HCURL_TET_In_FEM.hpp"
56 #include "Shards_CellTopology.hpp"
57 
58 #include <iostream>
59 using namespace Intrepid;
60 
65 int 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[] = {
119 1.000000000000001e+00, -2.498001805406602e-16, -1.665334536937735e-16,
120 9.999999999999998e-01, 1.000000000000000e+00, 1.000000000000000e+00,
121 5.828670879282072e-16, 1.110223024625157e-16, 2.498001805406602e-16,
122 7.771561172376096e-16, 8.326672684688674e-17, 1.110223024625157e-16,
123 2.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,
126 2.220446049250313e-16, -5.551115123125783e-17, 1.084013877651281e-16,
127 3.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,
130 5.551115123125783e-17, -2.220446049250313e-16, -8.326672684688674e-17,
131 -2.220446049250313e-16, -5.551115123125783e-17, 9.999999999999999e-01,
132 1.665334536937735e-16, 1.110223024625157e-16, -6.383782391594650e-16,
133 1.110223024625157e-16, 1.110223024625157e-16, -1.110223024625157e-16,
134 9.999999999999990e-01, 9.999999999999994e-01, 9.999999999999996e-01,
135 1.387778780781446e-16, -2.496931404305374e-17, -1.665334536937735e-16,
136 -2.498001805406602e-16, -2.149987498083074e-16, 1.000000000000000e+00,
137 8.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,
142 2.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 }
Intrepid::Basis::getCardinality
virtual int getCardinality() const
Returns cardinality of the basis.
Definition: Intrepid_BasisDef.hpp:100
Intrepid_HCURL_TET_In_FEM.hpp
Header file for the Intrepid::HCURL_TET_In_FEM class.
Intrepid_PointTools.hpp
Header file for utility class to provide point tools, such as barycentric coordinates, equispaced lattices, and warp-blend point distrubtions.
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid::Basis_HCURL_TET_In_FEM
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Te...
Definition: Intrepid_HCURL_TET_In_FEM.hpp:95
Intrepid::Basis::getBaseCellTopology
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
Definition: Intrepid_BasisDef.hpp:112
Intrepid::Basis_HCURL_TET_In_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Tetrahedron cell.
Definition: Intrepid_HCURL_TET_In_FEMDef.hpp:433
Intrepid::PointTools::getLatticeSize
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...
Definition: Intrepid_PointTools.hpp:214
Generated by   doxygen 1.8.5