Intrepid
test_01.cpp
1 #include "Teuchos_Array.hpp"
2 #include "Teuchos_RCP.hpp"
3 #include "Teuchos_oblackholestream.hpp"
4 #include "Teuchos_GlobalMPISession.hpp"
9 #include "Intrepid_Utils.hpp"
10 #include "Intrepid_Types.hpp"
11 
12 
13 using Teuchos::Array;
15 using Intrepid::Basis;
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 }
Contains definitions of custom data types in Intrepid.
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Intrepid utilities.
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin, Aeronautics, Imperial College London) within Intrepid.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...
Header file for the Intrepid::HGRAD_QUAD_Cn_FEM class.
An abstract base class that defines interface for bases that are tensor products of simpler bases...
Header file for the Intrepid::CubaturePolylib class.
Header file for the Intrepid::TensorProductSpaceTools class.