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_HDIV_TRI_In_FEM.hpp"
56 #include "Intrepid_HCURL_TRI_In_FEM.hpp"
57 #include "Shards_CellTopology.hpp"
58 
59 #include <iostream>
60 using namespace Intrepid;
61 
66 int 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[] = {
172 2.000000000000000e+00, 0.000000000000000e+00,
173 5.000000000000000e-01, 2.500000000000000e-01,
174 -1.000000000000000e+00, -1.000000000000000e+00,
175 2.500000000000000e-01, 0.000000000000000e+00,
176 -5.000000000000000e-01, -5.000000000000000e-01,
177 0.000000000000000e+00, 0.000000000000000e+00,
178 -1.000000000000000e+00, 0.000000000000000e+00,
179 5.000000000000000e-01, 2.500000000000000e-01,
180 2.000000000000000e+00, 2.000000000000000e+00,
181 -5.000000000000000e-01, 0.000000000000000e+00,
182 2.500000000000000e-01, 2.500000000000000e-01,
183 0.000000000000000e+00, 0.000000000000000e+00,
184 0.000000000000000e+00, 0.000000000000000e+00,
185 0.000000000000000e+00, 2.500000000000000e-01,
186 0.000000000000000e+00, 2.000000000000000e+00,
187 5.000000000000000e-01, 0.000000000000000e+00,
188 -2.500000000000000e-01, 2.500000000000000e-01,
189 1.000000000000000e+00, 0.000000000000000e+00,
190 0.000000000000000e+00, 0.000000000000000e+00,
191 0.000000000000000e+00, -5.000000000000000e-01,
192 0.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,
196 0.000000000000000e+00, 1.000000000000000e+00,
197 0.000000000000000e+00, 5.000000000000000e-01,
198 0.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,
202 0.000000000000000e+00, -2.000000000000000e+00,
203 0.000000000000000e+00, -2.500000000000000e-01,
204 0.000000000000000e+00, 0.000000000000000e+00,
205 -2.500000000000000e-01, -5.000000000000000e-01,
206 5.000000000000000e-01, 5.000000000000000e-01,
207 1.000000000000000e+00, 1.000000000000000e+00,
208 0.000000000000000e+00, 0.000000000000000e+00,
209 0.000000000000000e+00, -7.500000000000000e-01,
210 0.000000000000000e+00, 0.000000000000000e+00,
211 1.500000000000000e+00, 0.000000000000000e+00,
212 7.500000000000000e-01, 7.500000000000000e-01,
213 0.000000000000000e+00, 0.000000000000000e+00,
214 0.000000000000000e+00, 0.000000000000000e+00,
215 0.000000000000000e+00, 1.500000000000000e+00,
216 0.000000000000000e+00, 0.000000000000000e+00,
217 -7.500000000000000e-01, 0.000000000000000e+00,
218 7.500000000000000e-01, 7.500000000000000e-01,
219 0.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[] = {
263 7.000000000000000e+00,
264 2.500000000000000e+00,
265 -2.000000000000000e+00,
266 2.500000000000000e+00,
267 -2.000000000000000e+00,
268 -2.000000000000000e+00,
269 -2.000000000000000e+00,
270 2.500000000000000e+00,
271 7.000000000000000e+00,
272 -2.000000000000000e+00,
273 2.500000000000000e+00,
274 -2.000000000000000e+00,
275 -2.000000000000000e+00,
276 2.500000000000000e+00,
277 7.000000000000000e+00,
278 -2.000000000000000e+00,
279 2.500000000000000e+00,
280 -2.000000000000000e+00,
281 -2.000000000000000e+00,
282 -2.000000000000000e+00,
283 -2.000000000000000e+00,
284 2.500000000000000e+00,
285 2.500000000000000e+00,
286 7.000000000000000e+00,
287 -2.000000000000000e+00,
288 -2.000000000000000e+00,
289 -2.000000000000000e+00,
290 2.500000000000000e+00,
291 2.500000000000000e+00,
292 7.000000000000000e+00,
293 7.000000000000000e+00,
294 2.500000000000000e+00,
295 -2.000000000000000e+00,
296 2.500000000000000e+00,
297 -2.000000000000000e+00,
298 -2.000000000000000e+00,
299 -9.000000000000000e+00,
300 -4.500000000000000e+00,
301 0.000000000000000e+00,
302 0.000000000000000e+00,
303 4.500000000000000e+00,
304 9.000000000000000e+00,
305 9.000000000000000e+00,
306 0.000000000000000e+00,
307 -9.000000000000000e+00,
308 4.500000000000000e+00,
309 -4.500000000000000e+00,
310 0.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 }
Intrepid::Basis::getCardinality
virtual int getCardinality() const
Returns cardinality of the basis.
Definition: Intrepid_BasisDef.hpp:100
Intrepid_PointTools.hpp
Header file for utility class to provide point tools, such as barycentric coordinates, equispaced lattices, and warp-blend point distrubtions.
Intrepid::Basis_HDIV_TRI_In_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
Definition: Intrepid_HDIV_TRI_In_FEMDef.hpp:285
Intrepid::Basis_HCURL_TRI_In_FEM
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Tr...
Definition: Intrepid_HCURL_TRI_In_FEM.hpp:83
Intrepid_HDIV_TRI_In_FEM.hpp
Header file for the Intrepid::HDIV_TRI_In_FEM class.
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid::Basis_HDIV_TRI_In_FEM
Implementation of the default H(div)-compatible Raviart-Thomas basis of arbitrary degree on Triangle ...
Definition: Intrepid_HDIV_TRI_In_FEM.hpp:93
Intrepid_HCURL_TRI_In_FEM.hpp
Header file for the Intrepid::HCURL_TRI_In_FEM class.
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_TRI_In_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
Definition: Intrepid_HCURL_TRI_In_FEMDef.hpp:288
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