51 #include "Teuchos_oblackholestream.hpp"
52 #include "Teuchos_RCP.hpp"
53 #include "Teuchos_GlobalMPISession.hpp"
54 #include "Intrepid_HGRAD_TRI_Cn_FEM_ORTH.hpp"
57 #include "Shards_CellTopology.hpp"
59 using namespace Intrepid;
66 int main(
int argc,
char *argv[]) {
68 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
71 int iprint = argc - 1;
73 Teuchos::RCP<std::ostream> outStream;
74 Teuchos::oblackholestream bhs;
77 outStream = Teuchos::rcp(&std::cout,
false);
79 outStream = Teuchos::rcp(&bhs,
false);
82 Teuchos::oblackholestream oldFormatState;
83 oldFormatState.copyfmt(std::cout);
86 <<
"===============================================================================\n" \
88 <<
"| Unit Test OrthogonalBases |\n" \
90 <<
"| 1) Tests orthogonality of triangular orthogonal basis (Dubiner) |\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" \
96 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
97 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
99 <<
"===============================================================================\n";
118 myBasis.
getValues( basisAtCubPts , cubPts , OPERATOR_VALUE );
121 for (
int i=0;i<polydim;i++) {
122 for (
int j=i;j<polydim;j++) {
125 cur += cubWts(k) * basisAtCubPts( i , k ) * basisAtCubPts( j , k );
127 if (i != j && fabs( cur ) > INTREPID_TOL) {
130 else if (i == j && fabs( cur ) < INTREPID_TOL ) {
137 catch (
const std::exception & err) {
138 *outStream << err.what() <<
"\n\n";
146 shards::CellTopology myTri_3( shards::getCellTopologyData< shards::Triangle<3> >() );
149 PointTools::getLattice<double,FieldContainer<double> >( lattice ,
153 POINTTYPE_EQUISPACED );
157 myBasis.
getValues( dBasisAtLattice , lattice , OPERATOR_D1 );
159 double fiat_vals[] = {
160 0.000000000000000e+00, 0.000000000000000e+00,
161 0.000000000000000e+00, 0.000000000000000e+00,
162 0.000000000000000e+00, 0.000000000000000e+00,
163 0.000000000000000e+00, 0.000000000000000e+00,
164 0.000000000000000e+00, 0.000000000000000e+00,
165 0.000000000000000e+00, 0.000000000000000e+00,
166 0.000000000000000e+00, 0.000000000000000e+00,
167 0.000000000000000e+00, 0.000000000000000e+00,
168 0.000000000000000e+00, 0.000000000000000e+00,
169 0.000000000000000e+00, 0.000000000000000e+00,
170 3.464101615137754e+00, 1.732050807568877e+00,
171 3.464101615137754e+00, 1.732050807568877e+00,
172 3.464101615137754e+00, 1.732050807568877e+00,
173 3.464101615137754e+00, 1.732050807568877e+00,
174 3.464101615137754e+00, 1.732050807568877e+00,
175 3.464101615137754e+00, 1.732050807568877e+00,
176 3.464101615137754e+00, 1.732050807568877e+00,
177 3.464101615137754e+00, 1.732050807568877e+00,
178 3.464101615137754e+00, 1.732050807568877e+00,
179 3.464101615137754e+00, 1.732050807568877e+00,
180 0.000000000000000e+00, 3.000000000000000e+00,
181 0.000000000000000e+00, 3.000000000000000e+00,
182 0.000000000000000e+00, 3.000000000000000e+00,
183 0.000000000000000e+00, 3.000000000000000e+00,
184 0.000000000000000e+00, 3.000000000000000e+00,
185 0.000000000000000e+00, 3.000000000000000e+00,
186 0.000000000000000e+00, 3.000000000000000e+00,
187 0.000000000000000e+00, 3.000000000000000e+00,
188 0.000000000000000e+00, 3.000000000000000e+00,
189 0.000000000000000e+00, 3.000000000000000e+00,
190 -1.643167672515498e+01, -5.477225575051661e+00,
191 -5.477225575051661e+00, 0.000000000000000e+00,
192 5.477225575051660e+00, 5.477225575051660e+00,
193 1.643167672515498e+01, 1.095445115010332e+01,
194 -1.095445115010332e+01, -3.651483716701107e+00,
195 -9.121412916732176e-16, 1.825741858350553e+00,
196 1.095445115010332e+01, 7.302967433402213e+00,
197 -5.477225575051661e+00, -1.825741858350554e+00,
198 5.477225575051660e+00, 3.651483716701107e+00,
199 0.000000000000000e+00, 0.000000000000000e+00,
200 -4.242640687119285e+00, -1.272792206135786e+01,
201 -4.242640687119285e+00, -5.656854249492381e+00,
202 -4.242640687119285e+00, 1.414213562373094e+00,
203 -4.242640687119285e+00, 8.485281374238570e+00,
204 2.828427124746189e+00, -5.656854249492381e+00,
205 2.828427124746189e+00, 1.414213562373094e+00,
206 2.828427124746189e+00, 8.485281374238568e+00,
207 9.899494936611664e+00, 1.414213562373094e+00,
208 9.899494936611664e+00, 8.485281374238568e+00,
209 1.697056274847714e+01, 8.485281374238570e+00,
210 0.000000000000000e+00, -9.797958971132712e+00,
211 0.000000000000000e+00, -9.797958971132712e+00,
212 0.000000000000000e+00, -9.797958971132710e+00,
213 0.000000000000000e+00, -9.797958971132712e+00,
214 0.000000000000000e+00, -1.632993161855452e+00,
215 0.000000000000000e+00, -1.632993161855452e+00,
216 0.000000000000000e+00, -1.632993161855452e+00,
217 0.000000000000000e+00, 6.531972647421806e+00,
218 0.000000000000000e+00, 6.531972647421806e+00,
219 0.000000000000000e+00, 1.469693845669907e+01,
220 4.489988864128730e+01, 1.122497216032182e+01,
221 -4.988876515698587e+00, -6.236095644623235e+00,
222 -4.988876515698591e+00, 1.247219128924645e+00,
223 4.489988864128730e+01, 3.367491648096547e+01,
224 1.995550606279436e+01, 4.988876515698590e+00,
225 -4.988876515698589e+00, -2.494438257849295e+00,
226 1.995550606279435e+01, 1.496662954709576e+01,
227 4.988876515698590e+00, 1.247219128924648e+00,
228 4.988876515698586e+00, 3.741657386773940e+00,
229 0.000000000000000e+00, 0.000000000000000e+00,
230 1.897366596101028e+01, 2.846049894151541e+01,
231 6.324555320336759e+00, -7.378647873726218e+00,
232 -6.324555320336757e+00, -1.370320319406298e+01,
233 -1.897366596101028e+01, 9.486832980505138e+00,
234 -1.686548085423136e+01, 4.216370213557840e+00,
235 -1.404333387430680e-15, -2.108185106778921e+00,
236 1.686548085423135e+01, 2.108185106778919e+01,
237 -2.319003617456811e+01, -5.270462766947300e+00,
238 2.319003617456811e+01, 1.791957340762081e+01,
239 0.000000000000000e+00, 0.000000000000000e+00,
240 4.898979485566356e+00, 3.184336665618131e+01,
241 4.898979485566356e+00, 1.224744871391589e+01,
242 4.898979485566356e+00, -7.348469228349532e+00,
243 4.898979485566356e+00, -2.694438717061496e+01,
244 -3.265986323710904e+00, -4.898979485566356e+00,
245 -3.265986323710904e+00, -1.632993161855452e+00,
246 -3.265986323710904e+00, 1.632993161855451e+00,
247 1.143095213298816e+01, -7.348469228349537e+00,
248 1.143095213298816e+01, 1.877942136133769e+01,
249 4.898979485566356e+01, 2.449489742783178e+01,
250 0.000000000000000e+00, 2.121320343559643e+01,
251 0.000000000000000e+00, 2.121320343559643e+01,
252 0.000000000000000e+00, 2.121320343559642e+01,
253 0.000000000000000e+00, 2.121320343559643e+01,
254 0.000000000000000e+00, -4.714045207910317e+00,
255 0.000000000000000e+00, -4.714045207910317e+00,
256 0.000000000000000e+00, -4.714045207910317e+00,
257 0.000000000000000e+00, 2.357022603955157e+00,
258 0.000000000000000e+00, 2.357022603955157e+00,
259 0.000000000000000e+00, 4.242640687119285e+01
262 int fiat_index_cur = 0;
263 for (
int i=0;i<polydim;i++) {
264 for (
int j=0;j<np_lattice;j++) {
265 for (
int k=0;k<2;k++) {
266 if (std::abs( dBasisAtLattice(i,j,k) - fiat_vals[fiat_index_cur] ) > INTREPID_TOL ) {
268 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
271 *outStream <<
" At multi-index { ";
272 *outStream << i <<
" " << j <<
" " << k;
273 *outStream <<
"} computed value: " << dBasisAtLattice(i,j,k)
274 <<
" but correct value: " << fiat_vals[fiat_index_cur] <<
"\n";
275 *outStream <<
"Difference: " << std::abs( dBasisAtLattice(i,j,k) - fiat_vals[fiat_index_cur] ) <<
"\n";
282 catch (
const std::exception & err) {
283 *outStream << err.what() <<
"\n\n";
290 shards::CellTopology myTri_3( shards::getCellTopologyData< shards::Triangle<3> >() );
295 PointTools::getLattice<double,FieldContainer<double> >( lattice ,
299 POINTTYPE_EQUISPACED );
303 myBasis.
getValues( dBasisAtLattice , lattice , OPERATOR_D2 );
305 const double fiat_vals[] = {
306 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
307 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
308 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
309 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
310 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
311 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
312 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
313 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
314 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
315 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
316 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
317 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
318 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
319 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
320 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
321 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
322 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
323 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
324 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
325 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
326 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
327 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
328 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
329 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
330 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
331 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
332 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
333 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
334 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
335 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
336 3.286335345030997e+01, 1.643167672515498e+01, 5.477225575051661e+00,
337 3.286335345030997e+01, 1.643167672515498e+01, 5.477225575051661e+00,
338 3.286335345030997e+01, 1.643167672515498e+01, 5.477225575051661e+00,
339 3.286335345030997e+01, 1.643167672515498e+01, 5.477225575051661e+00,
340 3.286335345030997e+01, 1.643167672515498e+01, 5.477225575051661e+00,
341 3.286335345030997e+01, 1.643167672515498e+01, 5.477225575051661e+00,
342 3.286335345030997e+01, 1.643167672515498e+01, 5.477225575051661e+00,
343 3.286335345030997e+01, 1.643167672515498e+01, 5.477225575051661e+00,
344 3.286335345030997e+01, 1.643167672515498e+01, 5.477225575051661e+00,
345 3.286335345030997e+01, 1.643167672515498e+01, 5.477225575051661e+00,
346 0.000000000000000e+00, 2.121320343559642e+01, 2.121320343559642e+01,
347 0.000000000000000e+00, 2.121320343559642e+01, 2.121320343559642e+01,
348 0.000000000000000e+00, 2.121320343559642e+01, 2.121320343559642e+01,
349 0.000000000000000e+00, 2.121320343559642e+01, 2.121320343559642e+01,
350 0.000000000000000e+00, 2.121320343559642e+01, 2.121320343559642e+01,
351 0.000000000000000e+00, 2.121320343559642e+01, 2.121320343559642e+01,
352 0.000000000000000e+00, 2.121320343559642e+01, 2.121320343559642e+01,
353 0.000000000000000e+00, 2.121320343559642e+01, 2.121320343559642e+01,
354 0.000000000000000e+00, 2.121320343559642e+01, 2.121320343559642e+01,
355 0.000000000000000e+00, 2.121320343559642e+01, 2.121320343559642e+01,
356 0.000000000000000e+00, 0.000000000000000e+00, 2.449489742783178e+01,
357 0.000000000000000e+00, 0.000000000000000e+00, 2.449489742783178e+01,
358 0.000000000000000e+00, 0.000000000000000e+00, 2.449489742783177e+01,
359 0.000000000000000e+00, 0.000000000000000e+00, 2.449489742783178e+01,
360 0.000000000000000e+00, 0.000000000000000e+00, 2.449489742783178e+01,
361 0.000000000000000e+00, 0.000000000000000e+00, 2.449489742783178e+01,
362 0.000000000000000e+00, 0.000000000000000e+00, 2.449489742783178e+01,
363 0.000000000000000e+00, 0.000000000000000e+00, 2.449489742783177e+01,
364 0.000000000000000e+00, 0.000000000000000e+00, 2.449489742783177e+01,
365 0.000000000000000e+00, 0.000000000000000e+00, 2.449489742783178e+01,
366 -2.244994432064365e+02, -8.979977728257460e+01, -2.244994432064365e+01,
367 -7.483314773547885e+01, -1.496662954709577e+01, 7.483314773547880e+00,
368 7.483314773547882e+01, 5.986651818838305e+01, 3.741657386773942e+01,
369 2.244994432064365e+02, 1.346996659238619e+02, 6.734983296193094e+01,
370 -1.496662954709577e+02, -5.986651818838306e+01, -1.496662954709577e+01,
371 -1.246222254316567e-14, 1.496662954709576e+01, 1.496662954709576e+01,
372 1.496662954709576e+02, 8.979977728257458e+01, 4.489988864128730e+01,
373 -7.483314773547885e+01, -2.993325909419154e+01, -7.483314773547884e+00,
374 7.483314773547882e+01, 4.489988864128729e+01, 2.244994432064365e+01,
375 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00,
376 -3.794733192202055e+01, -1.517893276880822e+02, -9.486832980505139e+01,
377 -3.794733192202055e+01, -6.324555320336759e+01, -6.324555320336759e+00,
378 -3.794733192202055e+01, 2.529822128134703e+01, 8.221921916437785e+01,
379 -3.794733192202055e+01, 1.138419957660617e+02, 1.707629936490925e+02,
380 5.059644256269407e+01, -6.324555320336759e+01, -5.059644256269407e+01,
381 5.059644256269407e+01, 2.529822128134703e+01, 3.794733192202055e+01,
382 5.059644256269407e+01, 1.138419957660616e+02, 1.264911064067352e+02,
383 1.391402170474087e+02, 2.529822128134704e+01, -6.324555320336762e+00,
384 1.391402170474087e+02, 1.138419957660617e+02, 8.221921916437786e+01,
385 2.276839915321233e+02, 1.138419957660617e+02, 3.794733192202055e+01,
386 0.000000000000000e+00, -5.878775382679627e+01, -1.616663230236897e+02,
387 0.000000000000000e+00, -5.878775382679627e+01, -9.308061022576078e+01,
388 0.000000000000000e+00, -5.878775382679627e+01, -2.449489742783179e+01,
389 0.000000000000000e+00, -5.878775382679627e+01, 4.409081537009720e+01,
390 0.000000000000000e+00, 9.797958971132708e+00, -5.878775382679628e+01,
391 0.000000000000000e+00, 9.797958971132708e+00, 9.797958971132703e+00,
392 0.000000000000000e+00, 9.797958971132708e+00, 7.838367176906168e+01,
393 0.000000000000000e+00, 7.838367176906168e+01, 4.409081537009718e+01,
394 0.000000000000000e+00, 7.838367176906168e+01, 1.126765281680262e+02,
395 0.000000000000000e+00, 1.469693845669907e+02, 1.469693845669907e+02,
396 0.000000000000000e+00, 0.000000000000000e+00, -1.272792206135786e+02,
397 0.000000000000000e+00, 0.000000000000000e+00, -1.272792206135786e+02,
398 0.000000000000000e+00, 0.000000000000000e+00, -1.272792206135785e+02,
399 0.000000000000000e+00, 0.000000000000000e+00, -1.272792206135786e+02,
400 0.000000000000000e+00, 0.000000000000000e+00, -2.828427124746190e+01,
401 0.000000000000000e+00, 0.000000000000000e+00, -2.828427124746190e+01,
402 0.000000000000000e+00, 0.000000000000000e+00, -2.828427124746190e+01,
403 0.000000000000000e+00, 0.000000000000000e+00, 7.071067811865474e+01,
404 0.000000000000000e+00, 0.000000000000000e+00, 7.071067811865474e+01,
405 0.000000000000000e+00, 0.000000000000000e+00, 1.697056274847714e+02
408 int fiat_index_cur = 0;
409 for (
int i=0;i<polydim;i++) {
410 for (
int j=0;j<np_lattice;j++) {
411 for (
int k=0;k<3;k++) {
412 if (std::abs( dBasisAtLattice(i,j,k) - fiat_vals[fiat_index_cur] ) > 10.0*INTREPID_TOL ) {
414 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
417 *outStream <<
" At multi-index { ";
418 *outStream << i <<
" " << j <<
" " << k;
419 *outStream <<
"} computed value: " << dBasisAtLattice(i,j,k)
420 <<
" but correct value: " << fiat_vals[fiat_index_cur] <<
"\n";
421 *outStream <<
"Difference: " << std::abs( dBasisAtLattice(i,j,k) - fiat_vals[fiat_index_cur] ) <<
"\n";
428 catch (
const std::exception & err) {
429 *outStream << err.what() <<
"\n\n";
435 std::cout <<
"End Result: TEST FAILED\n";
437 std::cout <<
"End Result: TEST PASSED\n";
440 std::cout.copyfmt(oldFormatState);
Implementation of the default H(grad)-compatible orthogonal basis (Dubiner) of arbitrary degree on tr...
virtual int getCardinality() const
Returns cardinality of the basis.
virtual void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
Header file for utility class to provide multidimensional containers.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
virtual int getNumPoints() const
Returns the number of cubature points.
Header file for the Intrepid::CubatureDirectTriDefault class.
Defines direct integration rules on a triangle.