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