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"
9#include "Intrepid_Utils.hpp"
10#include "Intrepid_Types.hpp"
11
12
13using Teuchos::Array;
15using 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
31int 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}
Header file for the Intrepid::CubaturePolylib class.
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Header file for the Intrepid::HGRAD_QUAD_Cn_FEM class.
Header file for the Intrepid::TensorProductSpaceTools class.
Contains definitions of custom data types in Intrepid.
static const double INTREPID_TOL
General purpose tolerance in, e.g., internal Newton's method to invert ref to phys maps.
Intrepid utilities.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin,...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
An abstract base class that defines interface for bases that are tensor products of simpler bases.