Intrepid
test_02.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
44
52#include "Shards_CellTopology.hpp"
53
54#include "Teuchos_oblackholestream.hpp"
55#include "Teuchos_RCP.hpp"
56#include "Teuchos_GlobalMPISession.hpp"
57#include "Teuchos_ScalarTraits.hpp"
58
59using namespace std;
60using namespace Intrepid;
61using namespace shards;
62
63#define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
64{ \
65 ++nException; \
66 try { \
67 S ; \
68 } \
69 catch (const std::logic_error & err) { \
70 ++throwCounter; \
71 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
72 *outStream << err.what() << '\n'; \
73 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
74 }; \
75}
76
77
78
79int main(int argc, char *argv[]) {
80
81 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
82
84 // typedef shards::CellTopology CellTopology;
85
86 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
87 int iprint = argc - 1;
88
89 Teuchos::RCP<std::ostream> outStream;
90 Teuchos::oblackholestream bhs; // outputs nothing
91
92 if (iprint > 0)
93 outStream = Teuchos::rcp(&std::cout, false);
94 else
95 outStream = Teuchos::rcp(&bhs, false);
96
97 // Save the format state of the original std::cout.
98 Teuchos::oblackholestream oldFormatState;
99 oldFormatState.copyfmt(std::cout);
100
101 *outStream \
102 << "===============================================================================\n" \
103 << "| |\n" \
104 << "| Unit Test CellTools |\n" \
105 << "| |\n" \
106 << "| 1) Mapping to and from reference cells with base and extended topologies|\n" \
107 << "| using default initial guesses when computing the inverse F^{-1} |\n" \
108 << "| 2) Repeat all tests from 1) using user-defined initial guess for F^{-1} |\n" \
109 << "| 3) Exception testing |\n" \
110 << "| |\n" \
111 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \
112 << "| Denis Ridzal (dridzal@sandia.gov), or |\n" \
113 << "| Kara Peterson (kjpeter@sandia.gov) |\n" \
114 << "| |\n" \
115 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
116 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
117 << "| |\n" \
118 << "===============================================================================\n";
119
120 int errorFlag = 0;
121
122 // Collect all supported cell topologies
123 std::vector<shards::CellTopology> supportedTopologies;
124 supportedTopologies.push_back(shards::getCellTopologyData<Triangle<3> >() );
125 supportedTopologies.push_back(shards::getCellTopologyData<Triangle<6> >() );
126 supportedTopologies.push_back(shards::getCellTopologyData<Quadrilateral<4> >() );
127 supportedTopologies.push_back(shards::getCellTopologyData<Quadrilateral<9> >() );
128 supportedTopologies.push_back(shards::getCellTopologyData<Tetrahedron<4> >() );
129 supportedTopologies.push_back(shards::getCellTopologyData<Tetrahedron<10> >() );
130 supportedTopologies.push_back(shards::getCellTopologyData<Hexahedron<8> >() );
131 supportedTopologies.push_back(shards::getCellTopologyData<Hexahedron<20> >() );
132 supportedTopologies.push_back(shards::getCellTopologyData<Hexahedron<27> >() );
133 supportedTopologies.push_back(shards::getCellTopologyData<Wedge<6> >() );
134 supportedTopologies.push_back(shards::getCellTopologyData<Wedge<15> >() );
135 supportedTopologies.push_back(shards::getCellTopologyData<Wedge<18> >() );
136 supportedTopologies.push_back(shards::getCellTopologyData<Pyramid<5> >() );
137 supportedTopologies.push_back(shards::getCellTopologyData<Pyramid<13> >() );
138
139 // Declare iterator to loop over the cell topologies
140 std::vector<shards::CellTopology>::iterator topo_iterator;
141
142 // Test 1 scope
143 try{
144
145 *outStream \
146 << "\n"
147 << "===============================================================================\n"\
148 << "| Test 1: computing F(x) and F^{-1}(x) using default initial guesses. |\n"\
149 << "===============================================================================\n\n";
150 /*
151 * Test summary:
152 *
153 * A reference point set is mapped to physical frame and then back to reference frame.
154 * Test passes if the final set of points matches the first set of points. The cell workset
155 * is generated by perturbing randomly the cellWorkset of a reference cell with the specified
156 * cell topology.
157 *
158 */
159 // Declare arrays for cell workset and point sets. We will have 10 cells in the wset. and 10 pts per pt. set
160 FieldContainer<double> cellWorkset; // physical cell workset
161 FieldContainer<double> refPoints; // reference point set(s)
162 FieldContainer<double> physPoints; // physical point set(s)
163 FieldContainer<double> controlPoints; // preimages: physical points mapped back to ref. frame
164
165 // We will use cubature factory to get some points on the reference cells. Declare necessary arrays
167 FieldContainer<double> cubPoints;
168 FieldContainer<double> cubWeights;
169
170 // Initialize number of cells in the cell workset
171 int numCells = 10;
172
173
174 // Loop over cell topologies, make cell workset for each one by perturbing the cellWorkset & test methods
175 for(topo_iterator = supportedTopologies.begin(); topo_iterator != supportedTopologies.end(); ++topo_iterator){
176
177 // 1. Define a single reference point set using cubature factory with order 4 cubature
178 Teuchos::RCP<Cubature<double> > cellCubature = cubFactory.create( (*topo_iterator), 4);
179 int cubDim = cellCubature -> getDimension();
180 int numPts = cellCubature -> getNumPoints();
181 cubPoints.resize(numPts, cubDim);
182 cubWeights.resize(numPts);
183 cellCubature -> getCubature(cubPoints, cubWeights);
184
185 // 2. Define a cell workset by perturbing the cellWorkset of the reference cell with the specified topology
186 // 2.1 Resize dimensions of the rank-3 (C,N,D) cell workset array for the current topology
187 int numNodes = (*topo_iterator).getNodeCount();
188 int cellDim = (*topo_iterator).getDimension();
189 cellWorkset.resize(numCells, numNodes, cellDim);
190
191 // 2.2 Copy cellWorkset of the reference cell with the same topology to temp rank-2 (N,D) array
192 FieldContainer<double> refCellNodes(numNodes, cellDim );
193 CellTools::getReferenceSubcellNodes(refCellNodes, cellDim, 0, (*topo_iterator) );
194
195 // 2.3 Create randomly perturbed version of the reference cell and save in the cell workset array
196 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){
197
198 // Move vertices +/-0.125 along their axes. Gives nondegenerate cells for base and extended topologies
199 for(int nodeOrd = 0; nodeOrd < numNodes; nodeOrd++){
200 for(int d = 0; d < cellDim; d++){
201 double delta = Teuchos::ScalarTraits<double>::random()/16.0;
202 cellWorkset(cellOrd, nodeOrd, d) = refCellNodes(nodeOrd, d) + delta;
203 } // d
204 }// nodeOrd
205 }// cellOrd
206 /*
207 * 3.1 Test 1: single point set to single physical cell: map ref. point set in rank-2 (P,D) array
208 * to a physical point set in rank-2 (P,D) array for a specified cell ordinal. Use the cub.
209 * points array for this test. Resize physPoints and controlPoints to rank-2 (P,D) arrays.
210 */
211 physPoints.resize(numPts, cubDim);
212 controlPoints.resize(numPts, cubDim);
213
214 *outStream
215 << " Mapping a set of " << numPts << " points to one cell in a workset of " << numCells << " "
216 << (*topo_iterator).getName() << " cells. \n";
217
218 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){
219
220 // Forward map:: requires cell ordinal
221 CellTools::mapToPhysicalFrame(physPoints, cubPoints, cellWorkset, (*topo_iterator), cellOrd);
222 // Inverse map: requires cell ordinal
223 CellTools::mapToReferenceFrame(controlPoints, physPoints, cellWorkset, (*topo_iterator), cellOrd);
224
225 // Points in controlPoints should match the originals in cubPoints up to a tolerance
226 for(int pt = 0; pt < numPts; pt++){
227 for(int d = 0; d < cellDim; d++){
228
229 if( abs( controlPoints(pt, d) - cubPoints(pt, d) ) > 100.0*INTREPID_TOL ){
230 errorFlag++;
231 *outStream
232 << std::setw(70) << "^^^^----FAILURE!" << "\n"
233 << " Mapping a single point set to a single physical cell in a workset failed for: \n"
234 << " Cell Topology = " << (*topo_iterator).getName() << "\n"
235 << " Physical cell ordinal in workset = " << cellOrd << "\n"
236 << " Reference point ordinal = " << setprecision(12) << pt << "\n"
237 << " At reference point coordinate = " << setprecision(12) << d << "\n"
238 << " Original value = " << cubPoints(pt, d) << "\n"
239 << " F^{-1}F(P_d) = " << controlPoints(pt, d) <<"\n";
240 }
241 }// d
242 }// pt
243 }// cellOrd
244 /*
245 * 3.2 Test 2: single point set to multiple physical cells: map ref. point set in rank-2 (P,D) array
246 * to a physical point set in rank-3 (C, P,D) array for all cell ordinals. Use the cub.
247 * points array for this test. Resize physPoints and controlPoints to rank-3 (C,P,D) arrays.
248 */
249 physPoints.clear();
250 controlPoints.clear();
251 physPoints.resize(numCells, numPts, cubDim);
252 controlPoints.resize(numCells, numPts, cubDim);
253
254 *outStream
255 << " Mapping a set of " << numPts << " points to all cells in workset of " << numCells << " "
256 << (*topo_iterator).getName() << " cells. \n";
257
258 // Forward map: do not specify cell ordinal
259 CellTools::mapToPhysicalFrame(physPoints, cubPoints, cellWorkset, (*topo_iterator));
260 // Inverse map: do not specify cell ordinal
261 CellTools::mapToReferenceFrame(controlPoints, physPoints, cellWorkset, (*topo_iterator));
262
263 // Check: points in controlPoints should match the originals in cubPoints up to a tolerance
264 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){
265 for(int pt = 0; pt < numPts; pt++){
266 for(int d = 0; d < cellDim; d++){
267
268 if( abs( controlPoints(cellOrd, pt, d) - cubPoints(pt, d) ) > 100.0*INTREPID_TOL ){
269 errorFlag++;
270 *outStream
271 << std::setw(70) << "^^^^----FAILURE!" << "\n"
272 << " Mapping a single point set to all physical cells in a workset failed for: \n"
273 << " Cell Topology = " << (*topo_iterator).getName() << "\n"
274 << " Physical cell ordinal in workset = " << cellOrd << "\n"
275 << " Reference point ordinal = " << setprecision(12) << pt << "\n"
276 << " At reference point coordinate = " << setprecision(12) << d << "\n"
277 << " Original value = " << cubPoints(pt, d) << "\n"
278 << " F^{-1}F(P_d) = " << controlPoints(cellOrd, pt, d) <<"\n";
279 }
280 }// d
281 }// pt
282 }// cellOrd
283 /*
284 * 3.3 Test 3: multiple point sets to multiple physical cells: map ref. point sets in rank-3 (C,P,D) array
285 * to physical point sets in rank-3 (C, P,D) array for all cell ordinals. The (C,P,D) array
286 * with reference point sets is obtained by cloning the cubature point array.
287 */
288 physPoints.clear();
289 controlPoints.clear();
290 refPoints.resize(numCells, numPts, cubDim);
291 physPoints.resize(numCells, numPts, cubDim);
292 controlPoints.resize(numCells, numPts, cubDim);
293
294 // Clone cubature points in refPoints:
295 for(int c = 0; c < numCells; c++){
296 for(int pt = 0; pt < numPts; pt++){
297 for(int d = 0; d < cellDim; d++){
298 refPoints(c, pt, d) = cubPoints(pt, d);
299 }// d
300 }// pt
301 }// c
302
303 *outStream
304 << " Mapping " << numCells << " sets of " << numPts << " points to corresponding cells in workset of " << numCells << " "
305 << (*topo_iterator).getName() << " cells. \n";
306
307 // Forward map: do not specify cell ordinal
308 CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator));
309 // Inverse map: do not specify cell ordinal
310 CellTools::mapToReferenceFrame(controlPoints, physPoints, cellWorkset, (*topo_iterator));
311
312 // Check: points in controlPoints should match the originals in cubPoints up to a tolerance
313 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){
314 for(int pt = 0; pt < numPts; pt++){
315 for(int d = 0; d < cellDim; d++){
316
317 if( abs( controlPoints(cellOrd, pt, d) - cubPoints(pt, d) ) > 100.0*INTREPID_TOL ){
318 errorFlag++;
319 *outStream
320 << std::setw(70) << "^^^^----FAILURE!" << "\n"
321 << " Mapping multiple point sets to corresponding physical cells in a workset failed for: \n"
322 << " Cell Topology = " << (*topo_iterator).getName() << "\n"
323 << " Physical cell ordinal in workset = " << cellOrd << "\n"
324 << " Reference point ordinal = " << setprecision(12) << pt << "\n"
325 << " At reference point coordinate = " << setprecision(12) << d << "\n"
326 << " Original value = " << refPoints(cellOrd, pt, d) << "\n"
327 << " F^{-1}F(P_d) = " << controlPoints(cellOrd, pt, d) <<"\n";
328 }
329 }// d
330 }// pt
331 }// cellOrd
332 }// topo_iterator
333 }// try using default initial guesses branch
334
335 /*************************************************************************************************
336 * Wrap up test: check if the test broke down unexpectedly due to an exception *
337 ************************************************************************************************/
338
339 catch (const std::logic_error & err) {
340 *outStream << err.what() << "\n";
341 errorFlag = -1000;
342 };
343
344
345 // Test 2: repeat all of the above using the original points as user-defined initial guess
346 try{
347
348 *outStream \
349 << "\n"
350 << "===============================================================================\n"\
351 << "| Test 2: computing F(x) and F^{-1}(x) using user-defined initial guess. |\n"\
352 << "===============================================================================\n\n";
353 /*
354 * Test summary:
355 *
356 * Repeats all parts of Test 1 using user-defined initial guesses in the computation of F^{-1}.
357 * The guesses are simply the exact solutions (the original set of points we started with)
358 * and so, this tests runs much faster because Newton converges in a single iteration.
359 *
360 */
361
362 // Declare arrays for cell workset and point sets. We will have 10 cells in the wset. and 10 pts per pt. set
363 FieldContainer<double> cellWorkset; // physical cell workset
364 FieldContainer<double> physPoints; // physical point set(s)
365 FieldContainer<double> controlPoints; // preimages: physical points mapped back to ref. frame
366 FieldContainer<double> initialGuess; // User-defined initial guesses for F^{-1}
367
368 // We will use cubature factory to get some points on the reference cells. Declare necessary arrays
370 FieldContainer<double> cubPoints;
371 FieldContainer<double> cubWeights;
372
373 // Initialize number of cells in the cell workset
374 int numCells = 10;
375
376
377 // Loop over cell topologies, make cell workset for each one by perturbing the cellWorkset & test methods
378 for(topo_iterator = supportedTopologies.begin(); topo_iterator != supportedTopologies.end(); ++topo_iterator){
379
380 // 1. Define a single reference point set using cubature factory with order 6 cubature
381 Teuchos::RCP<Cubature<double> > cellCubature = cubFactory.create( (*topo_iterator), 4);
382 int cubDim = cellCubature -> getDimension();
383 int numPts = cellCubature -> getNumPoints();
384 cubPoints.resize(numPts, cubDim);
385 cubWeights.resize(numPts);
386 cellCubature -> getCubature(cubPoints, cubWeights);
387
388 // 2. Define a cell workset by perturbing the cellWorkset of the reference cell with the specified topology
389 // 2.1 Resize dimensions of the rank-3 (C,N,D) cell workset array for the current topology
390 int numNodes = (*topo_iterator).getNodeCount();
391 int cellDim = (*topo_iterator).getDimension();
392 cellWorkset.resize(numCells, numNodes, cellDim);
393
394 // 2.2 Copy cellWorkset of the reference cell with the same topology to temp rank-2 (N,D) array
395 FieldContainer<double> refCellNodes(numNodes, cellDim );
396 CellTools::getReferenceSubcellNodes(refCellNodes, cellDim, 0, (*topo_iterator) );
397
398 // 2.3 Create randomly perturbed version of the reference cell and save in the cell workset array
399 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){
400
401 // Move vertices +/-0.125 along their axes. Gives nondegenerate cells for base and extended topologies
402 for(int nodeOrd = 0; nodeOrd < numNodes; nodeOrd++){
403 for(int d = 0; d < cellDim; d++){
404 double delta = Teuchos::ScalarTraits<double>::random()/16.0;
405 cellWorkset(cellOrd, nodeOrd, d) = refCellNodes(nodeOrd, d) + delta;
406 } // d
407 }// nodeOrd
408 }// cellOrd
409 /*
410 * 3.1 Test 1: single point set to single physical cell: map ref. point set in rank-2 (P,D) array
411 * to a physical point set in rank-2 (P,D) array for a specified cell ordinal. Use the cub.
412 * points array for this test. Resize physPoints and controlPoints to rank-2 (P,D) arrays.
413 */
414 physPoints.resize(numPts, cubDim);
415 controlPoints.resize(numPts, cubDim);
416
417 *outStream
418 << " Mapping a set of " << numPts << " points to one cell in a workset of " << numCells << " "
419 << (*topo_iterator).getName() << " cells. \n";
420
421 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){
422
423 // Forward map:: requires cell ordinal
424 CellTools::mapToPhysicalFrame(physPoints, cubPoints, cellWorkset, (*topo_iterator), cellOrd);
425 // Inverse map: requires cell ordinal. Use cubPoints as initial guess
426 CellTools::mapToReferenceFrameInitGuess(controlPoints, cubPoints, physPoints, cellWorkset, (*topo_iterator), cellOrd);
427
428 // Points in controlPoints should match the originals in cubPoints up to a tolerance
429 for(int pt = 0; pt < numPts; pt++){
430 for(int d = 0; d < cellDim; d++){
431
432 if( abs( controlPoints(pt, d) - cubPoints(pt, d) ) > 100.0*INTREPID_TOL ){
433 errorFlag++;
434 *outStream
435 << std::setw(70) << "^^^^----FAILURE!" << "\n"
436 << " Mapping a single point set to a single physical cell in a workset failed for: \n"
437 << " Cell Topology = " << (*topo_iterator).getName() << "\n"
438 << " Physical cell ordinal in workset = " << cellOrd << "\n"
439 << " Reference point ordinal = " << setprecision(12) << pt << "\n"
440 << " At reference point coordinate = " << setprecision(12) << d << "\n"
441 << " Original value = " << cubPoints(pt, d) << "\n"
442 << " F^{-1}F(P_d) = " << controlPoints(pt, d) <<"\n";
443 }
444 }// d
445 }// pt
446 }// cellOrd
447 /*
448 * 3.2 Test 2: single point set to multiple physical cells: map ref. point set in rank-2 (P,D) array
449 * to a physical point set in rank-3 (C, P,D) array for all cell ordinals. Use the cub.
450 * points array for this test. Resize physPoints and controlPoints to rank-3 (C,P,D) arrays.
451 */
452 physPoints.clear();
453 controlPoints.clear();
454 physPoints.resize(numCells, numPts, cubDim);
455 controlPoints.resize(numCells, numPts, cubDim);
456
457 // Clone cubature points in initialGuess:
458 initialGuess.resize(numCells, numPts, cubDim);
459 for(int c = 0; c < numCells; c++){
460 for(int pt = 0; pt < numPts; pt++){
461 for(int d = 0; d < cellDim; d++){
462 initialGuess(c, pt, d) = cubPoints(pt, d);
463 }// d
464 }// pt
465 }// c
466
467 *outStream
468 << " Mapping a set of " << numPts << " points to all cells in workset of " << numCells << " "
469 << (*topo_iterator).getName() << " cells. \n";
470
471 // Forward map: do not specify cell ordinal
472 CellTools::mapToPhysicalFrame(physPoints, cubPoints, cellWorkset, (*topo_iterator));
473 // Inverse map: do not specify cell ordinal
474 CellTools::mapToReferenceFrameInitGuess(controlPoints, initialGuess, physPoints, cellWorkset, (*topo_iterator));
475
476 // Check: points in controlPoints should match the originals in cubPoints up to a tolerance
477 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){
478 for(int pt = 0; pt < numPts; pt++){
479 for(int d = 0; d < cellDim; d++){
480
481 if( abs( controlPoints(cellOrd, pt, d) - cubPoints(pt, d) ) > 100.0*INTREPID_TOL ){
482 errorFlag++;
483 *outStream
484 << std::setw(70) << "^^^^----FAILURE!" << "\n"
485 << " Mapping a single point set to all physical cells in a workset failed for: \n"
486 << " Cell Topology = " << (*topo_iterator).getName() << "\n"
487 << " Physical cell ordinal in workset = " << cellOrd << "\n"
488 << " Reference point ordinal = " << setprecision(12) << pt << "\n"
489 << " At reference point coordinate = " << setprecision(12) << d << "\n"
490 << " Original value = " << cubPoints(pt, d) << "\n"
491 << " F^{-1}F(P_d) = " << controlPoints(cellOrd, pt, d) <<"\n";
492 }
493 }// d
494 }// pt
495 }// cellOrd
496 /*
497 * 3.3 Test 3: multiple point sets to multiple physical cells: map ref. point sets in rank-3 (C,P,D) array
498 * to physical point sets in rank-3 (C, P,D) array for all cell ordinals. The initialGuess
499 * array from last test is used as the required (C,P,D) array of reference points for the
500 * forward map and as the user-defined initial guess array for the inverse map
501 */
502 physPoints.clear();
503 controlPoints.clear();
504 physPoints.resize(numCells, numPts, cubDim);
505 controlPoints.resize(numCells, numPts, cubDim);
506
507 *outStream
508 << " Mapping " << numCells << " sets of " << numPts << " points to corresponding cells in workset of " << numCells << " "
509 << (*topo_iterator).getName() << " cells. \n";
510
511 // Forward map: do not specify cell ordinal
512 CellTools::mapToPhysicalFrame(physPoints, initialGuess, cellWorkset, (*topo_iterator));
513 // Inverse map: do not specify cell ordinal
514 CellTools::mapToReferenceFrameInitGuess(controlPoints, initialGuess, physPoints, cellWorkset, (*topo_iterator));
515
516 // Check: points in controlPoints should match the originals in cubPoints up to a tolerance
517 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){
518 for(int pt = 0; pt < numPts; pt++){
519 for(int d = 0; d < cellDim; d++){
520
521 if( abs( controlPoints(cellOrd, pt, d) - cubPoints(pt, d) ) > 100.0*INTREPID_TOL ){
522 errorFlag++;
523 *outStream
524 << std::setw(70) << "^^^^----FAILURE!" << "\n"
525 << " Mapping multiple point sets to corresponding physical cells in a workset failed for: \n"
526 << " Cell Topology = " << (*topo_iterator).getName() << "\n"
527 << " Physical cell ordinal in workset = " << cellOrd << "\n"
528 << " Reference point ordinal = " << setprecision(12) << pt << "\n"
529 << " At reference point coordinate = " << setprecision(12) << d << "\n"
530 << " Original value = " << initialGuess(cellOrd, pt, d) << "\n"
531 << " F^{-1}F(P_d) = " << controlPoints(cellOrd, pt, d) <<"\n";
532 }
533 }// d
534 }// pt
535 }// cellOrd
536 } //topo-iterator
537 }// try user-defined initial guess
538
539 /*************************************************************************************************
540 * Wrap up test: check if the test broke down unexpectedly due to an exception *
541 ************************************************************************************************/
542
543 catch (const std::logic_error & err) {
544 *outStream << err.what() << "\n";
545 errorFlag = -1000;
546 };
547
548 *outStream \
549 << "\n"
550 << "===============================================================================\n"\
551 << "| Test 3: Exception testing - only when HAVE_INTREPID_DEBUG is defined. |\n"\
552 << "===============================================================================\n\n";
553 /*
554 * Test summary:
555 * Calls methods of CellTools class with incorrectly configured arguments. This test is run only
556 * in debug mode because many of the exceptions are checked only in that mode.
557 *
558 */
559
560 // Initialize throw counter for exception testing
561 int nException = 0;
562 int throwCounter = 0;
563
564 try {
565
566#ifdef HAVE_INTREPID_DEBUG
567 // Some arbitrary dimensions
568 int C = 10;
569 int P = 21;
570 int N;
571 int D;
572 int V;
573
574 // Array arguments
575 FieldContainer<double> jacobian;
576 FieldContainer<double> jacobianInv;
577 FieldContainer<double> jacobianDet;
579 FieldContainer<double> cellWorkset;
580 FieldContainer<double> physPoints;
581 FieldContainer<double> refPoints;
582 FieldContainer<double> initGuess;
583
584 /***********************************************************************************************
585 * Exception tests for setJacobian method *
586 **********************************************************************************************/
587
588 // Use the second cell topology for these tests (Triangle<6>)
589 topo_iterator = supportedTopologies.begin() + 1;
590 D = (*topo_iterator).getDimension();
591 N = (*topo_iterator).getNodeCount();
592 V = (*topo_iterator).getVertexCount();
593
594 // 1. incorrect jacobian rank
595 jacobian.resize(C, P, D);
596 points.resize(P, D);
597 cellWorkset.resize(C, N, D);
598 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ),
599 throwCounter, nException );
600
601 // 2. Incorrect cellWorkset rank
602 cellWorkset.resize(C, D);
603 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ),
604 throwCounter, nException );
605
606 // 3. Incorrect points rank
607 cellWorkset.resize(C, N, D);
608 points.resize(D);
609 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ),
610 throwCounter, nException );
611
612 // 4. points rank incompatible with whichCell = valid cell ordinal
613 points.resize(C, P, D);
614 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator), 0 ),
615 throwCounter, nException );
616
617 // 5. Non-matching dim
618 jacobian.resize(C, P, D, D);
619 points.resize(C, P, D - 1);
620 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ),
621 throwCounter, nException );
622
623 // 6. Non-matching dim
624 jacobian.resize(C, P, D, D);
625 points.resize(C, P - 1, D);
626 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ),
627 throwCounter, nException );
628
629 // 7. Non-matching dim
630 jacobian.resize(C, P, D, D);
631 points.resize(C - 1, P, D);
632 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ),
633 throwCounter, nException );
634
635 // 8. Non-matching dim
636 jacobian.resize(C, P, D, D);
637 points.resize(C, P, D);
638 cellWorkset.resize(C, N, D - 1);
639 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ),
640 throwCounter, nException );
641
642 // 9. Non-matching dim
643 jacobian.resize(C, P, D, D);
644 points.resize(C, P, D);
645 cellWorkset.resize(C - 1, N, D);
646 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ),
647 throwCounter, nException );
648
649 // 10. Incompatible ranks
650 jacobian.resize(C, D, D);
651 points.resize(C, P, D);
652 cellWorkset.resize(C, N, D);
653 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ),
654 throwCounter, nException );
655
656 /***********************************************************************************************
657 * Exception tests for setJacobianInv method *
658 **********************************************************************************************/
659
660 // 11. incompatible ranks
661 jacobian.resize(C, P, D, D);
662 jacobianInv.resize(P, D, D);
663 INTREPID_TEST_COMMAND( CellTools::setJacobianInv(jacobianInv, jacobian),
664 throwCounter, nException );
665
666 // 12. incorrect ranks
667 jacobian.resize(D, D);
668 jacobianInv.resize(D, D);
669 INTREPID_TEST_COMMAND( CellTools::setJacobianInv(jacobianInv, jacobian),
670 throwCounter, nException );
671
672 // 13. nonmatching dimensions
673 jacobian.resize(C, P, D, D - 1);
674 jacobianInv.resize(C, P, D, D);
675 INTREPID_TEST_COMMAND( CellTools::setJacobianInv(jacobianInv, jacobian),
676 throwCounter, nException );
677
678 // 14. nonmatching dimensions
679 jacobian.resize(C, P, D - 1, D);
680 jacobianInv.resize(C, P, D, D);
681 INTREPID_TEST_COMMAND( CellTools::setJacobianInv(jacobianInv, jacobian),
682 throwCounter, nException );
683
684 // 15. nonmatching dimensions
685 jacobian.resize(C, P - 1, D, D);
686 jacobianInv.resize(C, P, D, D);
687 INTREPID_TEST_COMMAND( CellTools::setJacobianInv(jacobianInv, jacobian),
688 throwCounter, nException );
689
690 // 16. nonmatching dimensions
691 jacobian.resize(C - 1, P, D, D);
692 jacobianInv.resize(C, P, D, D);
693 INTREPID_TEST_COMMAND( CellTools::setJacobianInv(jacobianInv, jacobian),
694 throwCounter, nException );
695
696 /***********************************************************************************************
697 * Exception tests for setJacobianDet method *
698 **********************************************************************************************/
699
700 // 17. Incompatible ranks
701 jacobian.resize(C, P, D, D);
702 jacobianDet.resize(C, P, D);
703 INTREPID_TEST_COMMAND( CellTools::setJacobianDet(jacobianDet, jacobian),
704 throwCounter, nException );
705
706 // 18. Incompatible ranks
707 jacobian.resize(P, D, D);
708 jacobianDet.resize(C, P);
709 INTREPID_TEST_COMMAND( CellTools::setJacobianDet(jacobianDet, jacobian),
710 throwCounter, nException );
711
712 // 19. Incorrect rank
713 jacobian.resize(D, D);
714 jacobianDet.resize(C, P);
715 INTREPID_TEST_COMMAND( CellTools::setJacobianDet(jacobianDet, jacobian),
716 throwCounter, nException );
717
718 // 20. Incorrect rank
719 jacobian.resize(C, P, D, D);
720 jacobianDet.resize(C);
721 INTREPID_TEST_COMMAND( CellTools::setJacobianDet(jacobianDet, jacobian),
722 throwCounter, nException );
723
724 // 21. Incorrect dimension
725 jacobian.resize(C, P, D, D);
726 jacobianDet.resize(C, P-1);
727 INTREPID_TEST_COMMAND( CellTools::setJacobianDet(jacobianDet, jacobian),
728 throwCounter, nException );
729
730 // 22. Incorrect dimension
731 jacobian.resize(C - 1, P, D, D);
732 jacobianDet.resize(C, P);
733 INTREPID_TEST_COMMAND( CellTools::setJacobianDet(jacobianDet, jacobian),
734 throwCounter, nException );
735
736 /***********************************************************************************************
737 * Exception tests for mapToPhysicalFrame method *
738 **********************************************************************************************/
739
740 // 23. Incorrect refPoint rank
741 refPoints.resize(P);
742 physPoints.resize(P, D);
743 cellWorkset.resize(C, N, D);
744 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator) ),
745 throwCounter, nException );
746 // 24. Incorrect workset rank
747 cellWorkset.resize(P, D);
748 refPoints.resize(P, D);
749 physPoints.resize(C, P, D);
750 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator) ),
751 throwCounter, nException );
752
753 // 25. Incompatible ranks
754 refPoints.resize(C, P, D);
755 physPoints.resize(P, D);
756 cellWorkset.resize(C, N, D);
757 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator) ),
758 throwCounter, nException );
759
760 // 26. Incompatible dimensions
761 refPoints.resize(C, P, D);
762 physPoints.resize(C, P, D - 1);
763 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator) ),
764 throwCounter, nException );
765
766 // 27. Incompatible dimensions
767 refPoints.resize(C, P, D);
768 physPoints.resize(C, P - 1, D);
769 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator) ),
770 throwCounter, nException );
771
772 // 28. Incompatible dimensions
773 refPoints.resize(C, P, D);
774 physPoints.resize(C - 1, P, D);
775 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator) ),
776 throwCounter, nException );
777
778 // 29. Incorrect physPoints rank when whichCell is valid cell ordinal
779 refPoints.resize(P, D);
780 physPoints.resize(C, P, D);
781 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator), 0 ),
782 throwCounter, nException );
783
784 /***********************************************************************************************
785 * Exception tests for mapToReferenceFrame method (with default initial guesses) *
786 **********************************************************************************************/
787
788 // 30. incompatible ranks
789 refPoints.resize(C, P, D);
790 physPoints.resize(P, D);
791 cellWorkset.resize(C, N, D);
792 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator) ),
793 throwCounter, nException );
794
795 // 31. Incompatible ranks with whichCell = valid cell ordinal
796 refPoints.resize(C, P, D);
797 physPoints.resize(C, P, D);
798 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator), 0 ),
799 throwCounter, nException );
800
801 // 32. Incompatible ranks with whichCell = -1 (default)
802 refPoints.resize(P, D);
803 physPoints.resize(P, D);
804 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator)),
805 throwCounter, nException );
806
807 // 33. Nonmatching dimensions
808 refPoints.resize(C, P, D - 1);
809 physPoints.resize(C, P, D);
810 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator)),
811 throwCounter, nException );
812
813 // 34. Nonmatching dimensions
814 refPoints.resize(C, P - 1, D);
815 physPoints.resize(C, P, D);
816 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator)),
817 throwCounter, nException );
818
819 // 35. Nonmatching dimensions
820 refPoints.resize(C - 1, P, D);
821 physPoints.resize(C, P, D);
822 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator)),
823 throwCounter, nException );
824
825 // 36. Incorrect rank for cellWorkset
826 refPoints.resize(C, P, D);
827 physPoints.resize(C, P, D);
828 cellWorkset.resize(C, N);
829 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator)),
830 throwCounter, nException );
831
832 /***********************************************************************************************
833 * Exception tests for mapToReferenceFrameInitGuess method (initial guess is a parameter) *
834 **********************************************************************************************/
835
836 // 37. Incompatible ranks
837 refPoints.resize(C, P, D);
838 physPoints.resize(C, P, D);
839 initGuess.resize(P, D);
840 cellWorkset.resize(C, N, D);
841 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, (*topo_iterator) ),
842 throwCounter, nException );
843
844 // 38. Incompatible ranks when whichCell is valid ordinal
845 refPoints.resize(P, D);
846 physPoints.resize(P, D);
847 initGuess.resize(C, P, D);
848 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, (*topo_iterator), 0),
849 throwCounter, nException );
850
851 // 39. Nonmatching dimensions
852 refPoints.resize(C, P, D);
853 physPoints.resize(C, P, D);
854 initGuess.resize(C, P, D - 1);
855 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, (*topo_iterator)),
856 throwCounter, nException );
857
858 // 40. Nonmatching dimensions
859 initGuess.resize(C, P - 1, D);
860 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, (*topo_iterator)),
861 throwCounter, nException );
862
863 // 41. Nonmatching dimensions
864 initGuess.resize(C - 1, P, D);
865 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, (*topo_iterator)),
866 throwCounter, nException );
867
868 /***********************************************************************************************
869 * Exception tests for mapToReferenceSubcell method *
870 **********************************************************************************************/
871
872 FieldContainer<double> refSubcellPoints;
873 FieldContainer<double> paramPoints;
874 int subcellDim = 2;
875 int subcellOrd = 0;
876
877 // This should set cell topology to Tetrahedron<10> so that we have real edges and faces.
878 topo_iterator += 5;
879 D = (*topo_iterator).getDimension();
880
881 // 42. Incorrect array rank
882 refSubcellPoints.resize(P,3);
883 paramPoints.resize(P);
884 INTREPID_TEST_COMMAND( CellTools::mapToReferenceSubcell(refSubcellPoints, paramPoints, subcellDim, subcellOrd, (*topo_iterator)),
885 throwCounter, nException );
886
887 // 43. Incorrect array rank
888 refSubcellPoints.resize(P);
889 paramPoints.resize(P, 2);
890 INTREPID_TEST_COMMAND( CellTools::mapToReferenceSubcell(refSubcellPoints, paramPoints, subcellDim, subcellOrd, (*topo_iterator)),
891 throwCounter, nException );
892
893 // 44. Incorrect array dimension for face of 3D cell (should be 3)
894 refSubcellPoints.resize(P, 2);
895 paramPoints.resize(P, 2);
896 INTREPID_TEST_COMMAND( CellTools::mapToReferenceSubcell(refSubcellPoints, paramPoints, subcellDim, subcellOrd, (*topo_iterator)),
897 throwCounter, nException );
898
899 // 45. Incorrect array dimension for parametrization domain of a face of 3D cell (should be 2)
900 refSubcellPoints.resize(P, 3);
901 paramPoints.resize(P, 3);
902 INTREPID_TEST_COMMAND( CellTools::mapToReferenceSubcell(refSubcellPoints, paramPoints, subcellDim, subcellOrd, (*topo_iterator)),
903 throwCounter, nException );
904
905 /***********************************************************************************************
906 * Exception tests for getReferenceEdgeTangent method *
907 **********************************************************************************************/
908
909 FieldContainer<double> refEdgeTangent;
910
911 // 46. Incorrect rank
912 refEdgeTangent.resize(C,P,D);
913 INTREPID_TEST_COMMAND( CellTools::getReferenceEdgeTangent(refEdgeTangent, 0, (*topo_iterator)),
914 throwCounter, nException );
915
916 // 47. Incorrect dimension D for Tet<10> cell
917 refEdgeTangent.resize(2);
918 INTREPID_TEST_COMMAND( CellTools::getReferenceEdgeTangent(refEdgeTangent, 0, (*topo_iterator)),
919 throwCounter, nException );
920
921 // 48. Invalid edge ordinal for Tet<10>
922 refEdgeTangent.resize(C,P,D);
923 INTREPID_TEST_COMMAND( CellTools::getReferenceEdgeTangent(refEdgeTangent, 10, (*topo_iterator)),
924 throwCounter, nException );
925
926 /***********************************************************************************************
927 * Exception tests for getReferenceFaceTangents method *
928 **********************************************************************************************/
929
930 FieldContainer<double> refFaceTanU;
931 FieldContainer<double> refFaceTanV;
932
933 // 49. Incorrect rank
934 refFaceTanU.resize(P, D);
935 refFaceTanV.resize(D);
936 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceTangents(refFaceTanU, refFaceTanV, 0, (*topo_iterator)),
937 throwCounter, nException );
938
939 // 50. Incorrect rank
940 refFaceTanU.resize(D);
941 refFaceTanV.resize(P, D);
942 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceTangents(refFaceTanU, refFaceTanV, 0, (*topo_iterator)),
943 throwCounter, nException );
944
945 // 51. Incorrect dimension for 3D cell
946 refFaceTanU.resize(D - 1);
947 refFaceTanV.resize(D);
948 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceTangents(refFaceTanU, refFaceTanV, 0, (*topo_iterator)),
949 throwCounter, nException );
950
951 // 52. Incorrect dimension for 3D cell
952 refFaceTanU.resize(D);
953 refFaceTanV.resize(D - 1);
954 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceTangents(refFaceTanU, refFaceTanV, 0, (*topo_iterator)),
955 throwCounter, nException );
956
957 // 53. Invalid face ordinal
958 refFaceTanU.resize(D);
959 refFaceTanV.resize(D);
960 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceTangents(refFaceTanU, refFaceTanV, 10, (*topo_iterator)),
961 throwCounter, nException );
962
963 /***********************************************************************************************
964 * Exception tests for getReferenceSide/FaceNormal methods *
965 **********************************************************************************************/
966
967 FieldContainer<double> refSideNormal;
968
969 // 54-55. Incorrect rank
970 refSideNormal.resize(C,P);
971 INTREPID_TEST_COMMAND( CellTools::getReferenceSideNormal(refSideNormal, 0, (*topo_iterator)),
972 throwCounter, nException );
973 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceNormal(refSideNormal, 0, (*topo_iterator)),
974 throwCounter, nException );
975
976 // 56-57. Incorrect dimension for 3D cell
977 refSideNormal.resize(D - 1);
978 INTREPID_TEST_COMMAND( CellTools::getReferenceSideNormal(refSideNormal, 0, (*topo_iterator)),
979 throwCounter, nException );
980 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceNormal(refSideNormal, 0, (*topo_iterator)),
981 throwCounter, nException );
982
983 // 58-59. Invalid side ordinal for Tet<10>
984 refSideNormal.resize(D);
985 INTREPID_TEST_COMMAND( CellTools::getReferenceSideNormal(refSideNormal, 10, (*topo_iterator)),
986 throwCounter, nException );
987 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceNormal(refSideNormal, 10, (*topo_iterator)),
988 throwCounter, nException );
989
990 // 60. Incorrect dimension for 2D cell: reset topo_iterator to the first cell in supportedTopologies which is Tri<3>
991 topo_iterator = supportedTopologies.begin();
992 D = (*topo_iterator).getDimension();
993 refSideNormal.resize(D - 1);
994 INTREPID_TEST_COMMAND( CellTools::getReferenceSideNormal(refSideNormal, 0, (*topo_iterator)),
995 throwCounter, nException );
996
997 // 61. Invalid side ordinal for Tri<3>
998 refSideNormal.resize(D);
999 INTREPID_TEST_COMMAND( CellTools::getReferenceSideNormal(refSideNormal, 10, (*topo_iterator)),
1000 throwCounter, nException );
1001
1002 // 62. Cannot call the "face" method for 2D cells
1003 refSideNormal.resize(D);
1004 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceNormal(refSideNormal, 0, (*topo_iterator)),
1005 throwCounter, nException );
1006
1007 /***********************************************************************************************
1008 * Exception tests for checkPoint/Pointset/PointwiseInclusion methods *
1009 **********************************************************************************************/
1010 points.resize(2,3,3,4);
1011 FieldContainer<int> inCell;
1012
1013 // 63. Point dimension does not match cell topology
1014 double * point = 0;
1015 INTREPID_TEST_COMMAND(CellTools::checkPointInclusion(point, (*topo_iterator).getDimension() + 1, (*topo_iterator) ),
1016 throwCounter, nException );
1017
1018 // 64. Invalid cell topology
1019 CellTopology pentagon_5(shards::getCellTopologyData<shards::Pentagon<> >() );
1020 INTREPID_TEST_COMMAND(CellTools::checkPointInclusion(point, pentagon_5.getDimension(), pentagon_5 ),
1021 throwCounter, nException );
1022
1023 // 65. Incorrect spatial dimension of points
1024 points.resize(10, 10, (*topo_iterator).getDimension() + 1);
1025 INTREPID_TEST_COMMAND(CellTools::checkPointsetInclusion(points, (*topo_iterator) ),
1026 throwCounter, nException );
1027
1028 // 66. Incorrect rank of input array
1029 points.resize(10,10,10,3);
1030 INTREPID_TEST_COMMAND(CellTools::checkPointsetInclusion(points, (*topo_iterator) ),
1031 throwCounter, nException );
1032
1033 // 67. Incorrect rank of output array
1034 points.resize(10,10,(*topo_iterator).getDimension() );
1035 inCell.resize(10);
1036 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ),
1037 throwCounter, nException );
1038
1039 // 68. Incorrect rank of output array
1040 points.resize(10, (*topo_iterator).getDimension() );
1041 inCell.resize(10, 10);
1042 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ),
1043 throwCounter, nException );
1044
1045 // 69. Incorrect rank of output array
1046 points.resize((*topo_iterator).getDimension() );
1047 inCell.resize(10, 10);
1048 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ),
1049 throwCounter, nException );
1050
1051 // 70. Incorrect dimension of output array
1052 points.resize(10, 10, (*topo_iterator).getDimension() );
1053 inCell.resize(10, 9);
1054 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ),
1055 throwCounter, nException );
1056
1057 // 71. Incorrect dimension of output array
1058 points.resize(10, 10, (*topo_iterator).getDimension() );
1059 inCell.resize(9, 10);
1060 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ),
1061 throwCounter, nException );
1062
1063 // 72. Incorrect dimension of output array
1064 points.resize(10, (*topo_iterator).getDimension() );
1065 inCell.resize(9);
1066 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ),
1067 throwCounter, nException );
1068
1069 // 73. Incorrect spatial dimension of input array
1070 points.resize(10, 10, (*topo_iterator).getDimension() + 1);
1071 inCell.resize(10, 10);
1072 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ),
1073 throwCounter, nException );
1074
1075 // 74. Incorrect rank of input array.
1076 points.resize(10,10,10,3);
1077 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ),
1078 throwCounter, nException );
1079
1080
1081 physPoints.resize(C, P, D);
1082 inCell.resize(C, P);
1083 // 75. Invalid rank of cellWorkset
1084 cellWorkset.resize(C, N, D, D);
1085 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator) ),
1086 throwCounter, nException );
1087
1088 // 76. Invalid dimension 1 (node count) of cellWorkset
1089 cellWorkset.resize(C, N + 1, D);
1090 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator) ),
1091 throwCounter, nException );
1092
1093 // 77. Invalid dimension 2 (spatial dimension) of cellWorkset
1094 cellWorkset.resize(C, N, D + 1);
1095 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator) ),
1096 throwCounter, nException );
1097
1098 // 78. Invalid whichCell value (exceeds cell count in the workset)
1099 cellWorkset.resize(C, N, D);
1100 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator), C + 1 ),
1101 throwCounter, nException );
1102
1103 // 79. Invalid whichCell for rank-3 physPoints (must be -1, here it is valid cell ordinal)
1104 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator), 0 ),
1105 throwCounter, nException );
1106
1107 // 80. Invalid whichCell for rank-2 physPoints (must be a valid cell ordinal, here it is the default -1)
1108 physPoints.resize(P, D);
1109 inCell.resize(P);
1110 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator) ),
1111 throwCounter, nException );
1112
1113 // 81. Incompatible ranks of I/O arrays
1114 physPoints.resize(C, P, D);
1115 inCell.resize(P);
1116 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator)),
1117 throwCounter, nException );
1118
1119 // 82. Incompatible ranks of I/O arrays
1120 physPoints.resize(P, D);
1121 inCell.resize(C, P);
1122 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator), 0),
1123 throwCounter, nException );
1124
1125 // 83. Incompatible dimensions of I/O arrays
1126 physPoints.resize(C, P, D);
1127 inCell.resize(C, P + 1);
1128 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator)),
1129 throwCounter, nException );
1130
1131 // 84. Incompatible dimensions of I/O arrays: rank-3 Input
1132 physPoints.resize(C + 1, P, D);
1133 inCell.resize(C, P);
1134 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator)),
1135 throwCounter, nException );
1136
1137 // 85. Incompatible dimensions of I/O arrays: rank-2 Input
1138 physPoints.resize(P, D);
1139 inCell.resize(P + 1);
1140 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator), 0 ),
1141 throwCounter, nException );
1142
1143
1144 /***********************************************************************************************
1145 * Exception tests for getReferenceVertex/vertices/Node/Nodes methods *
1146 **********************************************************************************************/
1147
1148 FieldContainer<double> subcellNodes;
1149
1150 // 86-89. Cell does not have reference cell
1151 INTREPID_TEST_COMMAND(CellTools::getReferenceVertex(pentagon_5, 0), throwCounter, nException);
1152 INTREPID_TEST_COMMAND(CellTools::getReferenceNode(pentagon_5, 0), throwCounter, nException);
1153 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellVertices(subcellNodes, 0, 0, pentagon_5), throwCounter, nException);
1154 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellNodes(subcellNodes, 0, 0, pentagon_5), throwCounter, nException);
1155
1156 // Use last cell topology (Wedge<18>) for these tests
1157 topo_iterator = supportedTopologies.end() - 1;
1158 D = (*topo_iterator).getDimension();
1159 int subcDim = D - 1;
1160 int S = (*topo_iterator).getSubcellCount(subcDim);
1161 V = (*topo_iterator).getVertexCount(subcDim, S - 1);
1162 subcellNodes.resize(V, D);
1163 // 90. subcell ordinal out of range
1164 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellVertices(subcellNodes, subcDim, S + 1, (*topo_iterator)),
1165 throwCounter, nException);
1166
1167 // 91. subcell dim out of range
1168 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellVertices(subcellNodes, D + 1, S, (*topo_iterator)),
1169 throwCounter, nException);
1170
1171 // 92. Incorrect rank for subcellNodes
1172 subcellNodes.resize(V, D, D);
1173 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellVertices(subcellNodes, subcDim, S - 1, (*topo_iterator)),
1174 throwCounter, nException);
1175
1176 // 93. Incorrect dimension for subcellNodes
1177 subcellNodes.resize(V - 1, D);
1178 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellVertices(subcellNodes, subcDim, S - 1, (*topo_iterator)),
1179 throwCounter, nException);
1180
1181 // 94. Incorrect dimension for subcellNodes
1182 subcellNodes.resize(V, D - 1);
1183 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellVertices(subcellNodes, subcDim, S - 1, (*topo_iterator)),
1184 throwCounter, nException);
1185
1186
1187 N = (*topo_iterator).getNodeCount(subcDim, S - 1);
1188 subcellNodes.resize(N, D);
1189 // 95. subcell ordinal out of range
1190 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellNodes(subcellNodes, subcDim, S + 1, (*topo_iterator)),
1191 throwCounter, nException);
1192
1193 // 96. subcell dim out of range
1194 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellNodes(subcellNodes, D + 1, S, (*topo_iterator)),
1195 throwCounter, nException);
1196
1197 // 97. Incorrect rank for subcellNodes
1198 subcellNodes.resize(N, D, D);
1199 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellNodes(subcellNodes, subcDim, S - 1, (*topo_iterator)),
1200 throwCounter, nException);
1201
1202 // 98. Incorrect dimension for subcellNodes
1203 subcellNodes.resize(N - 1, D);
1204 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellNodes(subcellNodes, subcDim, S - 1, (*topo_iterator)),
1205 throwCounter, nException);
1206
1207 // 99. Incorrect dimension for subcellNodes
1208 subcellNodes.resize(N, D - 1);
1209 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellNodes(subcellNodes, subcDim, S - 1, (*topo_iterator)),
1210 throwCounter, nException);
1211
1212#endif
1213 } // try exception testing
1214
1215 /*************************************************************************************************
1216 * Wrap up test: check if the test broke down unexpectedly due to an exception *
1217 ************************************************************************************************/
1218
1219 catch (const std::logic_error & err) {
1220 *outStream << err.what() << "\n";
1221 errorFlag = -1000;
1222 }
1223
1224 // Check if number of thrown exceptions matches the one we expect
1225 if (throwCounter != nException) {
1226 errorFlag++;
1227 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
1228 }
1229
1230
1231 if (errorFlag != 0)
1232 std::cout << "End Result: TEST FAILED\n";
1233 else
1234 std::cout << "End Result: TEST PASSED\n";
1235
1236 // reset format state of std::cout
1237 std::cout.copyfmt(oldFormatState);
1238
1239 return errorFlag;
1240}
1241
1242
1243
1244
1245
1246
1247
Header file for the Intrepid::CellTools class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
A stateless class for operations on cell data. Provides methods for:
static const double * getReferenceVertex(const shards::CellTopology &cell, const int vertexOrd)
Retrieves the Cartesian coordinates of a reference cell vertex.
static int checkPointsetInclusion(const ArrayPoint &points, const shards::CellTopology &cellTopo, const double &threshold=INTREPID_THRESHOLD)
Checks if a set of points belongs to a reference cell.
static void getReferenceEdgeTangent(ArrayEdgeTangent &refEdgeTangent, const int &edgeOrd, const shards::CellTopology &parentCell)
Computes constant tangent vectors to edges of 2D or 3D reference cells.
static void mapToReferenceSubcell(ArraySubcellPoint &refSubcellPoints, const ArrayParamPoint &paramPoints, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
static void getReferenceSideNormal(ArraySideNormal &refSideNormal, const int &sideOrd, const shards::CellTopology &parentCell)
Computes constant normal vectors to sides of 2D or 3D reference cells.
static void getReferenceSubcellNodes(ArraySubcellNode &subcellNodes, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Retrieves the Cartesian coordinates of all nodes of a reference subcell.
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.
static void getReferenceFaceTangents(ArrayFaceTangentU &refFaceTanU, ArrayFaceTangentV &refFaceTanV, const int &faceOrd, const shards::CellTopology &parentCell)
Computes pairs of constant tangent vectors to faces of a 3D reference cells.
static void setJacobianDet(ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F.
static const double * getReferenceNode(const shards::CellTopology &cell, const int nodeOrd)
Retrieves the Cartesian coordinates of a reference cell node.
static void mapToReferenceFrame(ArrayRefPoint &refPoints, const ArrayPhysPoint &physPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computes , the inverse of the reference-to-physical frame map using a default initial guess.
static void checkPointwiseInclusion(ArrayIncl &inRefCell, const ArrayPoint &points, const shards::CellTopology &cellTopo, const double &threshold=INTREPID_THRESHOLD)
Checks every point in a set for inclusion in a reference cell.
static void mapToReferenceFrameInitGuess(ArrayRefPoint &refPoints, const ArrayInitGuess &initGuess, const ArrayPhysPoint &physPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computation of , the inverse of the reference-to-physical frame map using user-supplied initial guess...
static void getReferenceFaceNormal(ArrayFaceNormal &refFaceNormal, const int &faceOrd, const shards::CellTopology &parentCell)
Computes constant normal vectors to faces of 3D reference cell.
static void setJacobianInv(ArrayJacInv &jacobianInv, const ArrayJac &jacobian)
Computes the inverse of the Jacobian matrix DF of the reference-to-physical frame map F.
static int checkPointInclusion(const Scalar *point, const int pointDim, const shards::CellTopology &cellTopo, const double &threshold=INTREPID_THRESHOLD)
Checks if a point belongs to a reference cell.
static void getReferenceSubcellVertices(ArraySubcellVert &subcellVertices, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Retrieves the Cartesian coordinates of all vertices of a reference subcell.
A factory class that generates specific instances of cubatures.
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > &degree)
Factory method.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
void clear()
Clears FieldContainer to trivial container (one with rank = 0 and size = 0)