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 
49 #include "Intrepid_CellTools.hpp"
50 #include "Intrepid_FieldContainer.hpp"
51 #include "Intrepid_DefaultCubatureFactory.hpp"
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 
59 using namespace std;
60 using namespace Intrepid;
61 using namespace shards;
62 
63 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
64 { \
65  ++nException; \
66  try { \
67  S ; \
68  } \
69  catch (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 
79 int main(int argc, char *argv[]) {
80 
81  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
82 
83  typedef CellTools<double> CellTools;
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
166  DefaultCubatureFactory<double> cubFactory;
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 (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
369  DefaultCubatureFactory<double> cubFactory;
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 (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;
578  FieldContainer<double> points;
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(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 
Intrepid::FieldContainer::clear
void clear()
Clears FieldContainer to trivial container (one with rank = 0 and size = 0)
Definition: Intrepid_FieldContainerDef.hpp:1048
Intrepid_CellTools.hpp
Header file for the Intrepid::CellTools class.
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid_DefaultCubatureFactory.hpp
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Intrepid::FieldContainer::resize
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
Definition: Intrepid_FieldContainerDef.hpp:1137
Intrepid::DefaultCubatureFactory
A factory class that generates specific instances of cubatures.
Definition: Intrepid_DefaultCubatureFactory.hpp:77
Intrepid::DefaultCubatureFactory::create
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > &degree)
Factory method.
Definition: Intrepid_DefaultCubatureFactoryDef.hpp:53
Intrepid::CellTools
A stateless class for operations on cell data. Provides methods for:
Definition: Intrepid_CellTools.hpp:111
Generated by   doxygen 1.8.5