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
51#include "Teuchos_oblackholestream.hpp"
52#include "Teuchos_GlobalMPISession.hpp"
53#include "Shards_CellTopology.hpp"
54
55using namespace std;
56using namespace Intrepid;
57
58
59#define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
60{ \
61 ++nException; \
62 try { \
63 S ; \
64 } \
65 catch (const std::logic_error & err) { \
66 ++throwCounter; \
67 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
68 *outStream << err.what() << '\n'; \
69 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
70 }; \
71}
72
73
74int main(int argc, char *argv[]) {
75
76 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
77
78 // This little trick lets us print to std::cout only if
79 // a (dummy) command-line argument is provided.
80 int iprint = argc - 1;
81 Teuchos::RCP<std::ostream> outStream;
82 Teuchos::oblackholestream bhs; // outputs nothing
83 if (iprint > 0)
84 outStream = Teuchos::rcp(&std::cout, false);
85 else
86 outStream = Teuchos::rcp(&bhs, false);
87
88 // Save the format state of the original std::cout.
89 Teuchos::oblackholestream oldFormatState;
90 oldFormatState.copyfmt(std::cout);
91
92 *outStream \
93 << "===============================================================================\n" \
94 << "| |\n" \
95 << "| Unit Test (PointTools) |\n" \
96 << "| |\n" \
97 << "| 1) Construction of equispaced and warped lattices on simplices |\n" \
98 << "| |\n" \
99 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
100 << "| Denis Ridzal (dridzal@sandia.gov) or |\n" \
101 << "| Robert Kirby (robert.c.kirby@ttu.edu) |\n" \
102 << "| |\n" \
103 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
104 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
105 << "| |\n" \
106 << "===============================================================================\n";
107
108
109
110 int errorFlag = 0;
111
112 shards::CellTopology myLine_2( shards::getCellTopologyData< shards::Line<2> >() );
113 shards::CellTopology myTri_3( shards::getCellTopologyData< shards::Triangle<3> >() );
114 shards::CellTopology myTet_4( shards::getCellTopologyData< shards::Tetrahedron<4> >() );
115
116
117 *outStream \
118 << "\n"
119 << "===============================================================================\n"\
120 << "| TEST 1: size of lattices |\n"\
121 << "===============================================================================\n";
122
123 try{
124 // first try the lattices with offset = 0. This is a spot-check
125
126 if (PointTools::getLatticeSize( myLine_2 , 4 , 0 ) != 5) {
127 errorFlag++;
128 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
129 *outStream << " size of 4th order lattice on a line with no offset: " << PointTools::getLatticeSize( myLine_2 , 4 , 0 ) << "\n";
130 *outStream << " should be 5\n";
131 }
132
133
134 if (PointTools::getLatticeSize( myTri_3 , 3 , 0 ) != 10) {
135 errorFlag++;
136 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
137 *outStream << " size of 3rd order lattice on a line with no offset: " << PointTools::getLatticeSize( myTri_3 , 3 , 0 ) << "\n";
138 *outStream << " should be 10\n";
139 }
140
141
142 if (PointTools::getLatticeSize( myTet_4 , 3 , 0 ) != 20) {
143 errorFlag++;
144 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
145 *outStream << " size of 3rd order lattice on a tet with no offset: " << PointTools::getLatticeSize( myTet_4 , 3 , 0 ) << "\n";
146 *outStream << " should be 20\n";
147 }
148
149
150 // check with the offset = 1
151 if (PointTools::getLatticeSize( myLine_2 , 3 , 1 ) != 2) {
152 errorFlag++;
153 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
154 *outStream << " size of 3rd order lattice on a line with 1 offset: " << PointTools::getLatticeSize( myLine_2 , 3 , 1 ) << "\n";
155 *outStream << " should be 2\n";
156 }
157
158 if (PointTools::getLatticeSize( myTri_3 , 4 , 1 ) != 3) {
159 errorFlag++;
160 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
161 *outStream << " size of 4th order lattice on a triangle with 1 offset: " << PointTools::getLatticeSize( myTri_3 , 4 , 1 ) << "\n";
162 *outStream << " should be 3\n";
163 }
164
165 if (PointTools::getLatticeSize( myTet_4 , 5 , 1 ) != 4) {
166 errorFlag++;
167 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
168 *outStream << " size of 5th order lattice on a tetrahedron with 1 offset: " << PointTools::getLatticeSize( myTet_4 , 5 , 1 ) << "\n";
169 *outStream << " should be 4\n";
170 }
171
172 }
173 catch (std::exception &err) {
174 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
175 *outStream << err.what() << '\n';
176 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
177 errorFlag = -1000;
178 };
179
180 // Now verify that we throw an exception on some of the non-supported cell types.
181
182 *outStream \
183 << "\n"
184 << "===============================================================================\n"\
185 << "| TEST 2: check for unsupported cell types \n"\
186 << "===============================================================================\n";
187 try{
188 try {
189 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Quadrilateral<4> >() , 3 , 0 );
190 }
191 catch (const std::invalid_argument & err) {
192 *outStream << err.what() << "\n";
193 }
194
195 try {
196 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Hexahedron<8> >() , 3 , 0 );
197 }
198 catch (const std::invalid_argument & err) {
199 *outStream << err.what() << "\n";
200 }
201 }
202 catch (std::exception &err) {
203 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
204 *outStream << err.what() << '\n';
205 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
206 errorFlag = -1000;
207 };
208
209 // Check for malformed point arrays
210
211
212 *outStream \
213 << "\n"
214 << "===============================================================================\n"\
215 << "| TEST 2: malformed point arrays \n"\
216 << "===============================================================================\n";
217 try{
218 // line: not enough points allocated
219 try {
220 FieldContainer<double> pts(3,1);
221 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Line<2> >() , 5 , 0 );
222 }
223 catch (const std::invalid_argument & err) {
224 *outStream << err.what() << "\n";
225 }
226 // line: wrong dimension for points
227 try {
228 FieldContainer<double> pts(6,2);
229 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Line<2> >() , 5 , 0 );
230 }
231 catch (const std::invalid_argument & err) {
232 *outStream << err.what() << "\n";
233 }
234 // triangle: too many points allocated
235 try {
236 FieldContainer<double> pts(4,2);
237 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Triangle<3> >() , 3 , 1 );
238 }
239 catch (const std::invalid_argument & err) {
240 *outStream << err.what() << "\n";
241 }
242 // triangle: wrong dimension for points
243 try {
244 FieldContainer<double> pts(6,1);
245 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Triangle<3> >() , 3 , 0 );
246 }
247 catch (const std::invalid_argument & err) {
248 *outStream << err.what() << "\n";
249 }
250 // tetrahedron: not enough points allocated
251 try {
252 FieldContainer<double> pts(4,2);
253 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Tetrahedron<4> >() , 2 , 0 );
254 }
255 catch (const std::invalid_argument & err) {
256 *outStream << err.what() << "\n";
257 }
258 // tetrahedron: wrong dimension for points
259 try {
260 FieldContainer<double> pts(4,2);
261 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Tetrahedron<4> >() , 1 , 0 );
262 }
263 catch (const std::invalid_argument & err) {
264 *outStream << err.what() << "\n";
265 }
266
267 }
268 catch (std::exception &err) {
269 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
270 *outStream << err.what() << '\n';
271 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
272 errorFlag = -1000;
273 };
274
275
276 *outStream \
277 << "\n"
278 << "===============================================================================\n"\
279 << "| TEST 3: check values of triangular lattice compared to Warburton's code \n"\
280 << "===============================================================================\n";
281 try {
282 // triangle case
283 const int order = 4;
284 const int offset = 0;
285 int numPts = PointTools::getLatticeSize( myTri_3 , order , offset );
286 int ptDim = 2;
287 FieldContainer<double> warpBlendPts( numPts , ptDim );
288 PointTools::getLattice<double,FieldContainer<double> >( warpBlendPts ,
289 myTri_3 ,
290 order ,
291 offset ,
292 POINTTYPE_WARPBLEND );
293 FieldContainer<double> verts( 1, 3 , 2 );
294 verts(0,0,0) = -1.0;
295 verts(0,0,1) = -1.0/sqrt(3.0);
296 verts(0,1,0) = 1.0;
297 verts(0,1,1) = -1.0/sqrt(3.0);
298 verts(0,2,0) = 0.0;
299 verts(0,2,1) = 2.0/sqrt(3.0);
300
301 // holds points on the equilateral triangle
302 FieldContainer<double> warpBlendMappedPts( numPts , ptDim );
303
304 CellTools<double>::mapToPhysicalFrame(warpBlendMappedPts ,
305 warpBlendPts ,
306 verts ,
307 myTri_3 ,
308 0 );
309
310 // Values from MATLAB code
311 double points[] = { -1.000000000000000 , -0.577350269189626 ,
312 -0.654653670707977 , -0.577350269189626 ,
313 -0.000000000000000 , -0.577350269189626 ,
314 0.654653670707977 , -0.577350269189626 ,
315 1.000000000000000 , -0.577350269189626 ,
316 -0.827326835353989 , -0.278271574919028 ,
317 -0.327375261332958 , -0.189010195256608 ,
318 0.327375261332958 , -0.189010195256608 ,
319 0.827326835353989, -0.278271574919028,
320 -0.500000000000000, 0.288675134594813,
321 0.000000000000000, 0.378020390513215,
322 0.500000000000000, 0.288675134594813,
323 -0.172673164646011, 0.855621844108654,
324 0.172673164646011, 0.855621844108654,
325 0, 1.154700538379252 };
326
327 // compare
328 for (int i=0;i<numPts;i++) {
329 for (int j=0;j<2;j++) {
330 int l = 2*i+j;
331 if (std::abs(warpBlendMappedPts(i,j) - points[l]) > INTREPID_TOL) {
332 errorFlag++;
333 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
334
335 // Output the multi-index of the value where the error is:
336 *outStream << " At multi-index { ";
337 *outStream << i << " ";*outStream << j << " ";
338 *outStream << "} computed value: " << warpBlendMappedPts(i,j)
339 << " but correct value: " << points[l] << "\n";
340 }
341 }
342 }
343
344
345 }
346 catch (std::exception &err) {
347 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
348 *outStream << err.what() << '\n';
349 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
350 errorFlag = -1000;
351 };
352
353
354 *outStream \
355 << "\n"
356 << "===============================================================================\n"\
357 << "| TEST 4: check values of tetrahedral lattice compared to Warburton's code \n"\
358 << "===============================================================================\n";
359 try {
360 // triangle case
361 const int order = 6;
362 const int offset = 0;
363 int numPts = PointTools::getLatticeSize( myTet_4 , order , offset );
364 int ptDim = 3;
365 FieldContainer<double> warpBlendPts( numPts , ptDim );
366 PointTools::getLattice<double,FieldContainer<double> >( warpBlendPts ,
367 myTet_4 ,
368 order ,
369 offset ,
370 POINTTYPE_WARPBLEND );
371
372 FieldContainer<double> verts(1,4,3);
373 verts(0,0,0) = -1.0;
374 verts(0,0,1) = -1.0/sqrt(3.0);
375 verts(0,0,2) = -1.0/sqrt(6.0);
376 verts(0,1,0) = 1.0;
377 verts(0,1,1) = -1.0/sqrt(3.0);
378 verts(0,1,2) = -1.0/sqrt(6.0);
379 verts(0,2,0) = 0.0;
380 verts(0,2,1) = 2.0 / sqrt(3.0);
381 verts(0,2,2) = -1.0/sqrt(6.0);
382 verts(0,3,0) = 0.0;
383 verts(0,3,1) = 0.0;
384 verts(0,3,2) = 3.0 / sqrt(6.0);
385
386
387 // points on the equilateral tet
388 FieldContainer<double> warpBlendMappedPts( numPts , ptDim );
389
390 CellTools<double>::mapToPhysicalFrame(warpBlendMappedPts ,
391 warpBlendPts ,
392 verts ,
393 myTet_4 ,
394 0 );
395
396 // Values from MATLAB code
397 double points[] = { -1.000000000000000, -0.577350269189626, -0.408248290463863,
398 -0.830223896278567, -0.577350269189626, -0.408248290463863,
399 -0.468848793470714, -0.577350269189626, -0.408248290463863,
400 -0.000000000000000, -0.577350269189626, -0.408248290463863,
401 0.468848793470714, -0.577350269189626, -0.408248290463863,
402 0.830223896278567, -0.577350269189626, -0.408248290463863,
403 1.000000000000000, -0.577350269189626, -0.408248290463863,
404 -0.915111948139283, -0.430319850411323, -0.408248290463863,
405 -0.660434383303427, -0.381301968982318, -0.408248290463863,
406 -0.239932664820086, -0.368405260495326, -0.408248290463863,
407 0.239932664820086, -0.368405260495326, -0.408248290463863,
408 0.660434383303426, -0.381301968982318, -0.408248290463863,
409 0.915111948139283, -0.430319850411323, -0.408248290463863,
410 -0.734424396735357, -0.117359831084509, -0.408248290463863,
411 -0.439014646886819, -0.023585152684228, -0.408248290463863,
412 -0.000000000000000, -0.000000000000000, -0.408248290463863,
413 0.439014646886819, -0.023585152684228, -0.408248290463863,
414 0.734424396735357, -0.117359831084509, -0.408248290463863,
415 -0.500000000000000, 0.288675134594813, -0.408248290463863,
416 -0.199081982066733, 0.391990413179555, -0.408248290463863,
417 0.199081982066733, 0.391990413179555, -0.408248290463863,
418 0.500000000000000, 0.288675134594813, -0.408248290463863,
419 -0.265575603264643, 0.694710100274135, -0.408248290463863,
420 -0.000000000000000, 0.762603937964635, -0.408248290463863,
421 0.265575603264643, 0.694710100274135, -0.408248290463863,
422 -0.084888051860716, 1.007670119600949, -0.408248290463863,
423 0.084888051860716, 1.007670119600949, -0.408248290463863,
424 0, 1.154700538379252, -0.408248290463863,
425 -0.915111948139284, -0.528340129596858, -0.269626682252082,
426 -0.660434383303427, -0.512000835787190, -0.223412180441618,
427 -0.239932664820086, -0.507701932958193, -0.211253047073435,
428 0.239932664820086, -0.507701932958193, -0.211253047073435,
429 0.660434383303426, -0.512000835787190, -0.223412180441618,
430 0.915111948139284, -0.528340129596858, -0.269626682252082,
431 -0.773622922202284, -0.315952535579882, -0.223412180441618,
432 -0.421605613935553, -0.243414114697549, -0.172119771139157,
433 -0.000000000000000, -0.224211101329670, -0.158541190167514,
434 0.421605613935553, -0.243414114697549, -0.172119771139157,
435 0.773622922202284, -0.315952535579882, -0.223412180441618,
436 -0.559649103902302, 0.046063183547205, -0.211253047073435,
437 -0.194172509561981, 0.112105550664835, -0.158541190167514,
438 0.194172509561981, 0.112105550664835, -0.158541190167514,
439 0.559649103902302, 0.046063183547205, -0.211253047073435,
440 -0.319716439082216, 0.461638749410988, -0.211253047073435,
441 -0.000000000000000, 0.486828229395098, -0.172119771139157,
442 0.319716439082216, 0.461638749410988, -0.211253047073435,
443 -0.113188538898858, 0.827953371367071, -0.223412180441618,
444 0.113188538898858, 0.827953371367071, -0.223412180441618,
445 -0.000000000000000, 1.056680259193716, -0.269626682252082,
446 -0.734424396735357, -0.424020123154587, 0.025434853622935,
447 -0.439014646886819, -0.392761897021160, 0.113846468290170,
448 -0.000000000000000, -0.384900179459751, 0.136082763487954,
449 0.439014646886819, -0.392761897021160, 0.113846468290170,
450 0.734424396735357, -0.424020123154587, 0.025434853622935,
451 -0.559649103902302, -0.183816888326860, 0.113846468290170,
452 -0.194172509561981, -0.112105550664835, 0.158541190167514,
453 0.194172509561981, -0.112105550664835, 0.158541190167514,
454 0.559649103902302, -0.183816888326860, 0.113846468290170,
455 -0.333333333333333, 0.192450089729875, 0.136082763487954,
456 -0.000000000000000, 0.224211101329670, 0.158541190167514,
457 0.333333333333333, 0.192450089729875, 0.136082763487954,
458 -0.120634457015483, 0.576578785348020, 0.113846468290170,
459 0.120634457015482, 0.576578785348020, 0.113846468290170,
460 -0.000000000000000, 0.848040246309174, 0.025434853622935,
461 -0.500000000000000, -0.288675134594813, 0.408248290463863,
462 -0.199081982066733, -0.254236708399899, 0.505654869247127,
463 0.199081982066733, -0.254236708399899, 0.505654869247127,
464 0.500000000000000, -0.288675134594813, 0.408248290463863,
465 -0.319716439082216, -0.045291699705599, 0.505654869247127,
466 -0.000000000000000, -0.000000000000000, 0.516359313417471,
467 0.319716439082216, -0.045291699705599, 0.505654869247127,
468 -0.120634457015483, 0.299528408105498, 0.505654869247127,
469 0.120634457015483, 0.299528408105498, 0.505654869247127,
470 -0.000000000000000, 0.577350269189626, 0.408248290463863,
471 -0.265575603264643, -0.153330146035039, 0.791061727304791,
472 -0.000000000000000, -0.130698866804872, 0.855072651347100,
473 0.265575603264643, -0.153330146035039, 0.791061727304791,
474 -0.113188538898858, 0.065349433402436, 0.855072651347100,
475 0.113188538898858, 0.065349433402436, 0.855072651347099,
476 0, 0.306660292070078, 0.791061727304791,
477 -0.084888051860717, -0.049010139592768, 1.086123263179808,
478 0.084888051860717, -0.049010139592768, 1.086123263179808,
479 0.000000000000000, 0.098020279185535, 1.086123263179808,
480 0, 0, 1.224744871391589
481 };
482
483 // compare
484 for (int i=0;i<numPts;i++) {
485 for (int j=0;j<ptDim;j++) {
486 int l = ptDim*i+j;
487 if (std::abs(warpBlendMappedPts(i,j) - points[l]) > INTREPID_TOL) {
488 errorFlag++;
489 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
490
491 // Output the multi-index of the value where the error is:
492 *outStream << " At multi-index { ";
493 *outStream << i << " ";*outStream << j << " ";
494 *outStream << "} computed value: " << warpBlendMappedPts(i,j)
495 << " but correct value: " << points[l] << "\n";
496 }
497 }
498 }
499
500
501 }
502 catch (std::exception &err) {
503 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
504 *outStream << err.what() << '\n';
505 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
506 errorFlag = -1000;
507 };
508
509
510
511 if (errorFlag != 0)
512 std::cout << "End Result: TEST FAILED\n";
513 else
514 std::cout << "End Result: TEST PASSED\n";
515
516 // reset format state of std::cout
517 std::cout.copyfmt(oldFormatState);
518
519 return errorFlag;
520}
Header file for utility class to provide multidimensional containers.
Header file for utility class to provide point tools, such as barycentric coordinates,...
static void mapToPhysicalFrame(ArrayPhysPoint &physPoints, const ArrayRefPoint &refPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computes F, the reference-to-physical frame map.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
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...