Intrepid
test_01.cpp
Go to the documentation of this file.
1 #include "Teuchos_Array.hpp"
2 #include "Teuchos_RCP.hpp"
3 #include "Teuchos_oblackholestream.hpp"
4 #include "Teuchos_GlobalMPISession.hpp"
5 #include "Intrepid_TensorProductSpaceTools.hpp"
6 #include "Intrepid_HGRAD_QUAD_Cn_FEM.hpp"
7 #include "Intrepid_HGRAD_HEX_Cn_FEM.hpp"
8 #include "Intrepid_CubaturePolylib.hpp"
9 #include "Intrepid_Utils.hpp"
10 #include "Intrepid_Types.hpp"
11 
12 
13 using Teuchos::Array;
14 using Intrepid::FieldContainer;
15 using Intrepid::Basis;
16 using Intrepid::TensorBasis;
17 
18 #define INTREPID_TEST_COMMAND( S ) \
19 { \
20  try { \
21  S ; \
22  } \
23  catch (const std::logic_error & err) { \
24  *outStream << "Expected Error ----------------------------------------------------------------\n"; \
25  *outStream << err.what() << '\n'; \
26  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
27  }; \
28 }
29 
30 
31 int main( int argc , char **argv )
32 {
33  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
34 
35  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
36  int iprint = argc - 1;
37 
38  Teuchos::RCP<std::ostream> outStream;
39  Teuchos::oblackholestream bhs; // outputs nothing
40 
41  if (iprint > 0)
42  outStream = Teuchos::rcp(&std::cout, false);
43  else
44  outStream = Teuchos::rcp(&bhs, false);
45 
46  // Save the format state of the original std::cout.
47  Teuchos::oblackholestream oldFormatState;
48  oldFormatState.copyfmt(std::cout);
49 
50 
51  *outStream \
52  << "===============================================================================\n" \
53  << "| |\n" \
54  << "| Unit Test TensorProductSpace Tools |\n" \
55  << "| |\n" \
56  << "| Tests sum-factored polynomial evaluation and integration |\n" \
57  << "| |\n" \
58  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
59  << "| Denis Ridzal (dridzal@sandia.gov) or |\n" \
60  << "| Robert Kirby (robert.c.kirby@ttu.edu) |\n" \
61  << "| |\n" \
62  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
63  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
64  << "| |\n" \
65  << "===============================================================================\n";
66 
67 
68 
69  int errorFlag = 0;
70 
71  Array<RCP<TensorBasis<double,FieldContainer<double> > > > basesByDim(4);
72  Array<RCP<FieldContainer<double> > > cubPtsByDim(4);
73 
74  Intrepid::CubaturePolylib<double> cpl(2,Intrepid::PL_GAUSS_LOBATTO);
75 
76  FieldContainer<double> cubPoints( cpl.getNumPoints() ,1 );
77  FieldContainer<double> cubWeights( cpl.getNumPoints() );
78 
79  cpl.getCubature( cubPoints, cubWeights );
80 
81  basesByDim[2] = Teuchos::rcp( new Intrepid::Basis_HGRAD_QUAD_Cn_FEM<double,FieldContainer<double> >( 2 , Intrepid::POINTTYPE_SPECTRAL ) );
82  basesByDim[3] = Teuchos::rcp( new Intrepid::Basis_HGRAD_HEX_Cn_FEM<double,FieldContainer<double> >( 2 , Intrepid::POINTTYPE_SPECTRAL ) );
83 
84 
85  // get points
86  FieldContainer<double> quadPts( cpl.getNumPoints() * cpl.getNumPoints() , 2 );
87  for (int j=0;j<cpl.getNumPoints();j++)
88  {
89  for (int i=0;i<cpl.getNumPoints();i++)
90  {
91  int index = j*cpl.getNumPoints() + i;
92  quadPts(index,0) = cubPoints(i,0);
93  quadPts(index,1) = cubPoints(j,0);
94  }
95  }
96 
97  FieldContainer<double> cubPts( cpl.getNumPoints() * cpl.getNumPoints() * cpl.getNumPoints() , 3 );
98  for (int k=0;k<cpl.getNumPoints();k++)
99  {
100  for (int j=0;j<cpl.getNumPoints();j++)
101  {
102  for (int i=0;i<cpl.getNumPoints();i++)
103  {
104  int index = k* cpl.getNumPoints() * cpl.getNumPoints() + j*cpl.getNumPoints() + i;
105  cubPts(index,0) = cubPoints(i,0);
106  cubPts(index,1) = cubPoints(j,0);
107  cubPts(index,2) = cubPoints(k,0);
108  }
109  }
110  }
111 
112  cubPtsByDim[2] = Teuchos::rcp( &quadPts , false );
113  cubPtsByDim[3] = Teuchos::rcp( &cubPts , false );
114 
115  int space_dim = 2;
116 
117  Array<Array<RCP<Basis<double,FieldContainer<double> > > > > &bases = basesByDim[space_dim]->getBases();
118 
119  FieldContainer<double> coeff(1,1,basesByDim[space_dim]->getCardinality());
120 
121 
122 
123  Array<RCP<FieldContainer<double> > > pts( space_dim );
124  pts[0] = Teuchos::rcp( &cubPoints, false );
125  for (int i=1;i<space_dim;i++)
126  {
127  pts[i] = pts[0];
128  }
129 
130  Array<RCP<FieldContainer<double> > > wts(space_dim);
131  wts[0] = Teuchos::rcp( &cubWeights , false );
132  for (int i=1;i<space_dim;i++)
133  {
134  wts[i] = wts[0];
135  }
136 
137  FieldContainer<double> Phix(bases[0][0]->getCardinality(),
138  cpl.getNumPoints() );
139  FieldContainer<double> Phiy(bases[0][1]->getCardinality(),
140  cpl.getNumPoints() );
141  FieldContainer<double> DPhix(bases[0][0]->getCardinality(),
142  cpl.getNumPoints(), 1 );
143  FieldContainer<double> DPhiy(bases[0][1]->getCardinality(),
144  cpl.getNumPoints(), 1 );
145 
146  bases[0][0]->getValues( Phix , cubPoints, Intrepid::OPERATOR_VALUE );
147  bases[0][1]->getValues( Phiy , cubPoints, Intrepid::OPERATOR_VALUE );
148  bases[0][0]->getValues( DPhix , cubPoints, Intrepid::OPERATOR_D1 );
149  bases[0][1]->getValues( DPhiy , cubPoints, Intrepid::OPERATOR_D1 );
150 
151  Array<RCP<FieldContainer<double> > > basisVals(2);
152  basisVals[0] = Teuchos::rcp( &Phix , false );
153  basisVals[1] = Teuchos::rcp( &Phiy , false );
154 
155  Array<RCP<FieldContainer<double> > > basisDVals(2);
156  basisDVals[0] = Teuchos::rcp( &DPhix , false );
157  basisDVals[1] = Teuchos::rcp( &DPhiy , false );
158 
159  FieldContainer<double> vals(1,1,pts[0]->size() * pts[1]->size() );
160 
161  // first basis function is the polynomial.
162  coeff(0,0,0) = 1.0;
163 
164  Intrepid::TensorProductSpaceTools::evaluate<double,FieldContainer<double>,FieldContainer<double>,FieldContainer<double> >( vals , coeff , basisVals );
165 
166  FieldContainer<double> grads( 1 , 1, pts[0]->size() * pts[1]->size() , 2 );
167 
168  Intrepid::TensorProductSpaceTools::evaluateGradient<double,FieldContainer<double>,FieldContainer<double>,FieldContainer<double> >( grads ,
169  coeff ,
170  basisVals ,
171  basisDVals );
172 
173  // confirm by comparing to actual gradients
174  FieldContainer<double> fullVals(basesByDim[space_dim]->getCardinality(),
175  basesByDim[space_dim]->getCardinality());
176  FieldContainer<double> fullGrads(basesByDim[space_dim]->getCardinality(),
177  basesByDim[space_dim]->getCardinality(),
178  space_dim );
179 
180 
181  basesByDim[space_dim]->getValues( fullVals ,
182  quadPts ,
183  Intrepid::OPERATOR_VALUE );
184  basesByDim[space_dim]->getValues( fullGrads ,
185  quadPts ,
186  Intrepid::OPERATOR_GRAD );
187 
188  for (int i=0;i<fullVals.dimension(1);i++)
189  {
190  if (std::abs( fullVals(0,i) - vals(0,0,i) ) > Intrepid::INTREPID_TOL )
191  {
192  errorFlag++;
193  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
194 
195  // Output the multi-index of the value where the error is:
196  *outStream << " Evaluating first bf At multi-index { ";
197  *outStream << i;
198  *outStream << "} brute force value: " << fullVals(0,i)
199  << " but tensor-product value: " << vals(0,i) << "\n";
200  *outStream << "Difference: " << std::abs( fullVals(0,i) - vals(0,i) ) << "\n";
201  }
202  }
203 
204  for (int i=0;i<fullGrads.dimension(1);i++)
205  {
206  for (int j=0;j<fullGrads.dimension(2);j++)
207  {
208  if (std::abs( fullGrads(0,i,j) - grads(0,0,i,j) ) > Intrepid::INTREPID_TOL )
209  {
210  errorFlag++;
211  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
212 
213  // Output the multi-index of the value where the error is:
214  *outStream << " Evaluating first bf At multi-index { ";
215  *outStream << i << " " << j;
216  *outStream << "} brute force value: " << fullGrads(0,i,j)
217  << " but tensor-product value: " << grads(0,0,i,j) << "\n";
218  *outStream << "Difference: " << std::abs( fullGrads(0,i,j) - grads(0,0,i,j) ) << "\n";
219  }
220  }
221  }
222 
223 
224  // now test moments.
225  // I've already evaluated the first basis function at the quadrature points.
226  // why not use it?
227 
228  FieldContainer<double> momentsNaive(1,basesByDim[2]->getCardinality());
229  for (int i=0;i<basesByDim[2]->getCardinality();i++)
230  {
231  momentsNaive(0,i) = 0.0;
232  for (int qpty=0;qpty<cubPoints.dimension(0);qpty++)
233  {
234  for (int qptx=0;qptx<cubPoints.dimension(0);qptx++)
235  {
236  momentsNaive(0,i) += cubWeights(qpty) * cubWeights(qptx) *
237  vals( 0, 0, qpty*cubPoints.dimension(0)+qptx )
238  * fullVals(i,qpty*cubPoints.dimension(0)+qptx);
239  }
240  }
241  }
242 
243  FieldContainer<double> momentsClever(1,1,basesByDim[space_dim]->getCardinality());
244  Intrepid::TensorProductSpaceTools::moments<double,FieldContainer<double>,FieldContainer<double>,FieldContainer<double>,FieldContainer<double> >( momentsClever ,
245  vals ,
246  basisVals ,
247  wts );
248  for (int j=0;j<momentsClever.dimension(0);j++)
249  {
250  if (std::abs( momentsClever(0,0,j) - momentsNaive(0,j) ) > Intrepid::INTREPID_TOL )
251  {
252  errorFlag++;
253  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
254 
255  // Output the multi-index of the value where the error is:
256  *outStream << " At multi-index { ";
257  *outStream << " " << j;
258  *outStream << "} brute force value: " << momentsNaive(0,j)
259  << " but sum-factored value: " << momentsClever(0,0,j) << "\n";
260  *outStream << "Difference: " << std::abs( momentsNaive(0,j) - momentsClever(0,0,j) ) << "\n";
261  }
262  }
263 
264  if (errorFlag != 0)
265  {
266  std::cout << "End Result: TEST FAILED\n";
267  }
268  else
269  {
270  std::cout << "End Result: TEST PASSED\n";
271  }
272 
273  // reset format state of std::cout
274  std::cout.copyfmt(oldFormatState);
275 
276  return errorFlag;
277 
278 }
main
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
Intrepid_CubaturePolylib.hpp
Header file for the Intrepid::CubaturePolylib class.
Intrepid_HGRAD_HEX_Cn_FEM.hpp
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Intrepid_HGRAD_QUAD_Cn_FEM.hpp
Header file for the Intrepid::HGRAD_QUAD_Cn_FEM class.
Intrepid_TensorProductSpaceTools.hpp
Header file for the Intrepid::TensorProductSpaceTools class.
Intrepid_Types.hpp
Contains definitions of custom data types in Intrepid.
Intrepid::INTREPID_TOL
static const double INTREPID_TOL
General purpose tolerance in, e.g., internal Newton's method to invert ref to phys maps.
Definition: Intrepid_Types.hpp:122
Intrepid_Utils.hpp
Intrepid utilities.
Intrepid::Basis_HGRAD_HEX_Cn_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell.
Definition: Intrepid_HGRAD_HEX_Cn_FEM.hpp:75
Intrepid::Basis_HGRAD_QUAD_Cn_FEM
Definition: Intrepid_HGRAD_QUAD_Cn_FEM.hpp:69
Intrepid::Basis
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Definition: Intrepid_Basis.hpp:89
Intrepid::CubaturePolylib
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin,...
Definition: Intrepid_CubaturePolylib.hpp:74
Intrepid::FieldContainer
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
Definition: Intrepid_FieldContainer.hpp:78
Intrepid::TensorBasis
An abstract base class that defines interface for bases that are tensor products of simpler bases.
Definition: Intrepid_TensorBasis.hpp:68
Generated by doxygen 1.9.1