Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/example/CellTools/example_01.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00044 
00049 #include "Intrepid_FieldContainer.hpp"
00050 #include "Intrepid_CellTools.hpp"
00051 #include "Intrepid_RealSpaceTools.hpp"
00052 #include "Shards_CellTopology.hpp"
00053 #include "Teuchos_GlobalMPISession.hpp"
00054 
00055 using namespace std;
00056 using namespace Intrepid;
00057 using namespace shards;
00058 
00059 int main(int argc, char *argv[]) {
00060 
00061   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00062 
00063   std::cout \
00064   << "===============================================================================\n" \
00065   << "|                                                                             |\n" \
00066   << "|                   Example use of the CellTools class                        |\n" \
00067   << "|                                                                             |\n" \
00068   << "|     1) Using shards::CellTopology to get cell types and topology            |\n" \
00069   << "|     2) Using CellTools to get cell Jacobians and their inverses and dets    |\n" \
00070   << "|     3) Testing points for inclusion in reference and physical cells         |\n" \
00071   << "|     4) Mapping points to and from reference cells with base and extended    |\n" \
00072   << "|        topologies.                                                          |\n" \
00073   << "|                                                                             |\n" \
00074   << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00075   << "|                      Denis Ridzal (dridzal@sandia.gov)                      |\n" \
00076   << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00077   << "|                                                                             |\n" \
00078   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00079   << "|  Shards's website:   http://trilinos.sandia.gov/packages/shards             |\n" \
00080   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00081   << "|                                                                             |\n" \
00082   << "===============================================================================\n"\
00083   << "| EXAMPLE 1: Query of cell types and topology                                 |\n"\
00084   << "===============================================================================\n";
00085   
00086   
00087   typedef CellTools<double>       CellTools;
00088   typedef RealSpaceTools<double>  RealSpaceTools;
00089   typedef shards::CellTopology    CellTopology;
00090   
00091   
00092   // Vector to hold cell topologies
00093   std::vector<CellTopology> shardsTopologies;
00094 
00095   
00096   // All 2D cell topologies in Shards:
00097   int cellDim = 2;
00098   shards::getTopologies(shardsTopologies, cellDim);
00099   std::cout << "Number of all " << cellDim 
00100     << "D Shards cell topologies = " << shardsTopologies.size() << "\n\n";
00101   
00102   for(unsigned i = 0; i < shardsTopologies.size(); i++){
00103     std::cout << shardsTopologies[i].getName() << "\n"; 
00104   }
00105   std::cout <<"\n";
00106 
00107   
00108   // All standard 2D cell topologies in Shards:
00109   shards::getTopologies(shardsTopologies, cellDim, 
00110                         shards::STANDARD_CELL);
00111   std::cout << "Number of all " 
00112     << cellDim << "D standard Shards cell topologies = " << shardsTopologies.size() << "\n\n";
00113   
00114   for(unsigned i = 0; i < shardsTopologies.size(); i++){
00115     std::cout << shardsTopologies[i].getName() << "\n"; 
00116   }
00117   std::cout <<"\n";
00118   
00119   
00120   // All standard 2D cells with base topologies in Shards:
00121   shards::getTopologies(shardsTopologies, cellDim, 
00122                         shards::STANDARD_CELL, 
00123                         shards::BASE_TOPOLOGY);
00124   std::cout << "Number of all " << cellDim 
00125     << "D standard Shards cells with base topologies = " << shardsTopologies.size() << "\n\n";
00126   
00127   for(unsigned i = 0; i < shardsTopologies.size(); i++){
00128     std::cout << shardsTopologies[i].getName() << "\n"; 
00129   }
00130   std::cout <<"\n";
00131   
00132   
00133   // All standard 2D cells with extended topologies in Shards:
00134   shards::getTopologies(shardsTopologies, cellDim, 
00135                         shards::STANDARD_CELL, 
00136                         shards::EXTENDED_TOPOLOGY);
00137   std::cout << "Number of all " << cellDim 
00138     << "D standard Shards cells with extended topologies = " << shardsTopologies.size() << "\n\n";
00139   
00140   for(unsigned i = 0; i < shardsTopologies.size(); i++){
00141     std::cout << shardsTopologies[i].getName() << "\n"; 
00142   }
00143   std::cout <<"\n";
00144 
00145   
00146   // All non-standard 2D cells with base topologies in Shards:
00147   shards::getTopologies(shardsTopologies, cellDim, 
00148                         shards::NONSTANDARD_CELL, 
00149                         shards::BASE_TOPOLOGY);
00150   std::cout << "Number of all " << cellDim 
00151     << "D non-standard Shards cells with base topologies = " << shardsTopologies.size() << "\n\n";
00152   
00153   for(unsigned i = 0; i < shardsTopologies.size(); i++){
00154     std::cout << shardsTopologies[i].getName() << "\n"; 
00155   }
00156   std::cout <<"\n";
00157   
00158  
00159   // All non-standard 2D cells with extended topologies in Shards:
00160   shards::getTopologies(shardsTopologies, cellDim, 
00161                         shards::NONSTANDARD_CELL, 
00162                         shards::EXTENDED_TOPOLOGY);
00163   std::cout << "Number of all " << cellDim 
00164     << "D non-standard Shards cells with extended topologies = " << shardsTopologies.size() << "\n\n";
00165   
00166   for(unsigned i = 0; i < shardsTopologies.size(); i++){
00167     std::cout << shardsTopologies[i].getName() << "\n"; 
00168   }
00169   std::cout <<"\n";
00170   
00171   // Finally, get all shards cell topologies and print them!
00172   shards::getTopologies(shardsTopologies);
00173   std::cout << "Number of all Shards cell topologies = " << shardsTopologies.size() << "\n\n";
00174   for(unsigned i = 0; i < shardsTopologies.size(); i++){
00175     std::cout << shardsTopologies[i] << "\n"; 
00176   }
00177   std::cout <<"\n";
00178   
00179   
00180   
00181 cout \
00182 << "===============================================================================\n"\
00183 << "| EXAMPLE 2: Using CellTools to get Jacobian, Jacobian inverse & Jacobian det |\n"\
00184 << "===============================================================================\n";
00185 
00186 // 4 triangles with basic triangle topology: number of nodes = number of vertices
00187 CellTopology triangle_3(shards::getCellTopologyData<Triangle<3> >() );
00188 int numCells = 4;
00189 int numNodes = triangle_3.getNodeCount();
00190 int spaceDim = triangle_3.getDimension();
00191 
00192 // Rank-3 array with dimensions (C,N,D) for the node coordinates of 3 traingle cells
00193 FieldContainer<double> triNodes(numCells, numNodes, spaceDim);
00194 
00195 // Initialize node data: accessor is (cellOrd, nodeOrd, coordinateOrd)
00196 triNodes(0, 0, 0) = 0.0;  triNodes(0, 0, 1) = 0.0;    // 1st triangle =  the reference tri
00197 triNodes(0, 1, 0) = 1.0;  triNodes(0, 1, 1) = 0.0;
00198 triNodes(0, 2, 0) = 0.0;  triNodes(0, 2, 1) = 1.0;
00199 
00200 triNodes(1, 0, 0) = 1.0;  triNodes(1, 0, 1) = 1.0;    // 2nd triangle = 1st shifted by (1,1)
00201 triNodes(1, 1, 0) = 2.0;  triNodes(1, 1, 1) = 1.0;
00202 triNodes(1, 2, 0) = 1.0;  triNodes(1, 2, 1) = 2.0;
00203 
00204 triNodes(2, 0, 0) = 0.0;  triNodes(2, 0, 1) = 0.0;    // 3rd triangle = flip 1st vertical
00205 triNodes(2, 1, 0) =-1.0;  triNodes(2, 1, 1) = 0.0;
00206 triNodes(2, 2, 0) = 0.0;  triNodes(2, 2, 1) = 1.0;
00207 
00208 triNodes(3, 0, 0) = 2.0;  triNodes(3, 0, 1) = 1.0;    // 4th triangle = just a triangle
00209 triNodes(3, 1, 0) = 3.0;  triNodes(3, 1, 1) = 0.5;
00210 triNodes(3, 2, 0) = 3.5;  triNodes(3, 2, 1) = 2.0;
00211 
00212 
00213 // Rank-2 array with dimensions (P,D) for some points on the reference triangle
00214 int numRefPoints = 2;
00215 FieldContainer<double> refPoints(numRefPoints, spaceDim);
00216 refPoints(0,0) = 0.0;   refPoints(0,1) = 0.0;
00217 refPoints(1,0) = 0.5;   refPoints(1,1) = 0.5;
00218 
00219 
00220 // Rank-4 array (C,P,D,D) for the Jacobian and its inverse and Rank-2 array (C,P) for its determinant
00221 FieldContainer<double> triJacobian(numCells, numRefPoints, spaceDim, spaceDim);
00222 FieldContainer<double> triJacobInv(numCells, numRefPoints, spaceDim, spaceDim);
00223 FieldContainer<double> triJacobDet(numCells, numRefPoints);
00224 
00225 // Rank-4 and Rank-2 auxiliary arrays
00226 FieldContainer<double> rank4Aux (numCells, numRefPoints, spaceDim, spaceDim);
00227 FieldContainer<double> rank2Aux (numCells, numRefPoints);
00228 
00229 
00230 // Methods to compute cell Jacobians, their inverses and their determinants
00231 CellTools::setJacobian(triJacobian, refPoints, triNodes, triangle_3); 
00232 CellTools::setJacobianInv(triJacobInv, triJacobian );
00233 CellTools::setJacobianDet(triJacobDet, triJacobian );
00234 
00235 // Checks: compute det(Inv(DF)) and Inv(Inv(DF))
00236 RealSpaceTools::det(rank2Aux, triJacobInv);
00237 RealSpaceTools::inverse(rank4Aux, triJacobInv);
00238 
00239 // Print data
00240 std::cout 
00241 << std::scientific<< std::setprecision(4)
00242 << std::right << std::setw(16) << "DF(P)" 
00243 << std::right << std::setw(30) << "Inv(DF(P))"
00244 << std::right << std::setw(30) << "Inv(Inv(DF(P)))\n";
00245 
00246 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){
00247   std::cout 
00248   << "===============================================================================\n"
00249   << "Cell " << cellOrd << "\n";
00250   for(int pointOrd = 0; pointOrd < numRefPoints; pointOrd++){
00251     std::cout << "Point = ("
00252     << std::setw(4) << std::right << refPoints(pointOrd, 0) << ","<< std::setw(4) << std::right<< refPoints(pointOrd,1) << ")\n";
00253     for(int row = 0; row < spaceDim; row++){
00254       std::cout 
00255       << std::setw(11) << std::right << triJacobian(cellOrd, pointOrd, row, 0) << " "
00256       << std::setw(11) << std::right << triJacobian(cellOrd, pointOrd, row, 1) 
00257       //
00258       << std::setw(16) << std::right << triJacobInv(cellOrd, pointOrd, row, 0) << " "
00259       << std::setw(11) << std::right << triJacobInv(cellOrd, pointOrd, row, 1)
00260       //
00261       << std::setw(16) << std::right << rank4Aux(cellOrd, pointOrd, row, 0) << " "
00262       << std::setw(11) << std::right << rank4Aux(cellOrd, pointOrd, row, 1) << "\n";
00263     }
00264     std::cout 
00265       << setw(5)<<std::left<< "Determinant:\n"
00266       << std::setw(11) << std::right << triJacobDet(cellOrd, pointOrd)
00267       << std::setw(28) << std::right << rank2Aux(cellOrd, pointOrd)
00268       << std::setw(28) << std::right << " product = " << triJacobDet(cellOrd, pointOrd)*rank2Aux(cellOrd, pointOrd);
00269     std::cout<< "\n\n";
00270   }
00271 }
00272 
00273 
00274 // 2 Quadrilateral cells with base topology 
00275 CellTopology quad_4(shards::getCellTopologyData<Quadrilateral<4> >() );
00276 numCells = 2;
00277 numNodes = quad_4.getNodeCount();
00278 spaceDim = quad_4.getDimension();
00279 
00280 FieldContainer<double> quadNodes(numCells, numNodes, spaceDim);
00281 // 1st QUAD
00282 quadNodes(0,0,0) = 1.00;  quadNodes(0,0,1) = 1.00;               
00283 quadNodes(0,1,0) = 2.00;  quadNodes(0,1,1) = 0.75;
00284 quadNodes(0,2,0) = 1.75;  quadNodes(0,2,1) = 2.00;  
00285 quadNodes(0,3,0) = 1.25;  quadNodes(0,3,1) = 2.00, 
00286 // 2ND QUAD
00287 quadNodes(1,0,0) = 2.00;  quadNodes(1,0,1) = 0.75;               
00288 quadNodes(1,1,0) = 3.00;  quadNodes(1,1,1) = 1.25;
00289 quadNodes(1,2,0) = 2.75;  quadNodes(1,2,1) = 2.25;
00290 quadNodes(1,3,0) = 1.75;  quadNodes(1,3,1) = 2.00;
00291 
00292 
00293 
00294 // 1 Hexahedron cell with base topology: number of nodes = number of vertices
00295 CellTopology hex_8(shards::getCellTopologyData<Hexahedron<8> >() );
00296 numCells = 1;
00297 numNodes = hex_8.getNodeCount();
00298 spaceDim = hex_8.getDimension();
00299 
00300 FieldContainer<double> hexNodes(numCells, numNodes, spaceDim);
00301 hexNodes(0,0,0) =
00302 // bottom face vertices
00303 hexNodes(0,0,0) = 1.00;   hexNodes(0,0,1) = 1.00;   hexNodes(0,0,2) = 0.00;          
00304 hexNodes(0,1,0) = 2.00;   hexNodes(0,1,1) = 0.75;   hexNodes(0,1,2) =-0.25;
00305 hexNodes(0,2,0) = 1.75;   hexNodes(0,2,1) = 2.00;   hexNodes(0,2,2) = 0.00;
00306 hexNodes(0,3,0) = 1.25;   hexNodes(0,3,1) = 2.00;   hexNodes(0,3,1) = 0.25;
00307 // top face vertices
00308 hexNodes(0,4,0) = 1.25;   hexNodes(0,4,1) = 0.75;   hexNodes(0,4,2) = 0.75;          
00309 hexNodes(0,5,0) = 1.75;   hexNodes(0,5,1) = 1.00;   hexNodes(0,5,2) = 1.00;
00310 hexNodes(0,6,0) = 2.00;   hexNodes(0,6,1) = 2.00;   hexNodes(0,6,2) = 1.25;
00311 hexNodes(0,7,0) = 1.00;   hexNodes(0,7,1) = 2.00;   hexNodes(0,7,2) = 1.00;
00312 
00313 
00314 
00315 std::cout << std::setprecision(16) << "\n" \
00316 << "===============================================================================\n"\
00317 << "| EXAMPLE 3: Using single point inclusion test method                         |\n"\
00318 << "===============================================================================\n";
00319 
00320 // Define cell topologies
00321 CellTopology edge3(shards::getCellTopologyData<Line<3> >() );
00322 CellTopology tri6 (shards::getCellTopologyData<Triangle<6> >() );
00323 CellTopology quad9(shards::getCellTopologyData<Quadrilateral<9> >() );
00324 CellTopology tet4 (shards::getCellTopologyData<Tetrahedron<> >() );
00325 CellTopology hex27(shards::getCellTopologyData<Hexahedron<27> >() );
00326 CellTopology wedge(shards::getCellTopologyData<Wedge<> >() );
00327 CellTopology pyr  (shards::getCellTopologyData<Pyramid<> >() );
00328 
00329 // Points that are close to the boundaries of their reference cells
00330 double point_in_edge[1]    = {1.0-INTREPID_EPSILON};
00331 double point_in_quad[2]    = {1.0,                  1.0-INTREPID_EPSILON};
00332 double point_in_tri[2]     = {0.5-INTREPID_EPSILON, 0.5-INTREPID_EPSILON};
00333 double point_in_tet[3]     = {0.5-INTREPID_EPSILON, 0.5-INTREPID_EPSILON, 2.0*INTREPID_EPSILON};
00334 double point_in_hex[3]     = {1.0-INTREPID_EPSILON, 1.0-INTREPID_EPSILON, 1.0-INTREPID_EPSILON};
00335 double point_in_wedge[3]   = {0.5,                  0.25,                 1.0-INTREPID_EPSILON};
00336 double point_in_pyramid[3] = {-INTREPID_EPSILON,    INTREPID_EPSILON,     1.0-INTREPID_EPSILON};
00337 
00338 // Run the inclusion test for each point and print results
00339 int in_edge     = CellTools::checkPointInclusion( point_in_edge,    1, edge3, INTREPID_THRESHOLD );
00340 int in_quad     = CellTools::checkPointInclusion( point_in_quad,    2, quad9 );
00341 int in_tri      = CellTools::checkPointInclusion( point_in_tri,     2, tri6 );
00342 int in_tet      = CellTools::checkPointInclusion( point_in_tet,     3, tet4 );
00343 int in_hex      = CellTools::checkPointInclusion( point_in_hex,     3, hex27 );
00344 int in_prism    = CellTools::checkPointInclusion( point_in_wedge,   3, wedge );
00345 int in_pyramid  = CellTools::checkPointInclusion( point_in_pyramid, 3, pyr );
00346 
00347 if(in_edge) {
00348   std::cout << "(" << point_in_edge[0] << ")" 
00349   << " is inside reference Line " << endl;
00350 }
00351 if(in_quad) {
00352   std::cout << "(" << point_in_quad[0] << "," << point_in_quad[1]  << ")" 
00353   << " is inside reference Quadrilateral " << endl;
00354 }
00355 if(in_tri) {
00356   std::cout << "(" << point_in_tri[0] << "," << point_in_tri[1] << ")"  
00357   << " is inside reference Triangle " << endl;
00358 }
00359 if(in_tet) {
00360   std::cout << "(" << point_in_tet[0] << "," << point_in_tet[1] << "," << point_in_tet[2]<<")" 
00361   << " is inside reference Tetrahedron " << endl;
00362 }
00363 if(in_hex) {
00364   std::cout << "(" << point_in_hex[0] << "," << point_in_hex[1] << "," << point_in_hex[2]<<")" 
00365   << " is inside reference Hexahedron " << endl;
00366 }
00367 if(in_prism) {
00368   std::cout << "(" << point_in_wedge[0] << "," << point_in_wedge[1] << "," << point_in_wedge[2]<<")" 
00369   << " is inside reference Wedge " << endl;
00370 }
00371 if(in_pyramid) {
00372   std::cout << "(" << point_in_pyramid[0] << "," << point_in_pyramid[1] << "," << point_in_pyramid[2]<<")" 
00373   << " is inside reference Pyramid " << endl;
00374 }
00375 
00376 // Change the points to be outside their reference cells.
00377 double small = 2.0*INTREPID_THRESHOLD;
00378 point_in_edge[0] += small;
00379 
00380 point_in_tri[0] += small;
00381 point_in_tri[1] += small;
00382 
00383 point_in_pyramid[0] += small;
00384 point_in_pyramid[1] += small;
00385 point_in_pyramid[2] += small;
00386 
00387 in_edge     = CellTools::checkPointInclusion(point_in_edge,    1,   edge3);
00388 in_tri      = CellTools::checkPointInclusion(point_in_tri,     2,   tri6);
00389 in_pyramid  = CellTools::checkPointInclusion(point_in_pyramid, 3,   pyr);
00390 
00391 std::cout << "\nChecking if perturbed Points belong to reference cell: " << endl;
00392 if(!in_edge) {
00393   std::cout << "(" << point_in_edge[0] << ")" << " is NOT inside reference Line " << endl;
00394 }
00395 if(!in_tri) {
00396   std::cout << "(" << point_in_tri[0] << "," << point_in_tri[1] << ")"  << " is NOT inside reference Triangle " << endl;
00397 }
00398 if(!in_pyramid) {
00399   std::cout << "(" << point_in_pyramid[0] << "," << point_in_pyramid[1] << "," << point_in_pyramid[2]<<")" 
00400   << " is NOT inside reference Pyramid " << endl;
00401 }
00402 
00403 
00404 
00405 std::cout << std::setprecision(6) << "\n" \
00406 << "===============================================================================\n"\
00407 << "| EXAMPLE 4-A: Using pointwise inclusion test method for reference cells      |\n"\
00408 << "===============================================================================\n";
00409 
00410 
00411 // Rank-1 array for one 2D reference point and rank-1 array for the test result
00412 FieldContainer<double> onePoint(2);
00413 FieldContainer<int> testOnePoint(1);
00414 
00415 onePoint(0) = 0.2;   onePoint(1) = 0.3;
00416 
00417 std::cout <<"\t Pointwise inclusion test for Triangle<6>: rank-1 array with a single 2D point: \n";
00418 
00419 CellTools::checkPointwiseInclusion(testOnePoint, onePoint, tri6);  
00420 
00421 std::cout << "point(" 
00422 << std::setw(13) << std::right << onePoint(0) << "," 
00423 << std::setw(13) << std::right << onePoint(1) << ") ";
00424 if( testOnePoint(0) ) {
00425   std::cout << " is inside. \n";
00426 }
00427 else{
00428   std::cout << " is not inside. \n";
00429 }
00430 std::cout << "\n";
00431 
00432 
00433 // Rank-2 array for 4 2D reference points (vector of points) and rank-1 array for the test result
00434 FieldContainer<double>  fourPoints(4, 2);
00435 FieldContainer<int> testFourPoints(4);
00436 
00437 fourPoints(0,0) = 0.5;   fourPoints(0,1) = 0.5;
00438 fourPoints(1,0) = 1.0;   fourPoints(1,1) = 1.1;
00439 fourPoints(2,0) =-1.0;   fourPoints(2,1) =-1.1;
00440 fourPoints(3,0) =-1.0;   fourPoints(3,1) = 0.5;
00441 
00442 std::cout <<"\t  Pointwise inclusion test for Quadrilateral<9>: rank-2 array with 4 2D points: \n";
00443 
00444 CellTools::checkPointwiseInclusion(testFourPoints, fourPoints, quad9);
00445 
00446 for(int i1 = 0; i1 < fourPoints.dimension(0); i1++) {
00447   std::cout << " point(" << i1 << ") = (" 
00448   << std::setw(13) << std::right << fourPoints(i1, 0) << "," 
00449   << std::setw(13) << std::right << fourPoints(i1, 1) << ") ";
00450   if( testFourPoints(i1) ) {
00451     std::cout << " is inside. \n";
00452   }
00453   else{
00454     std::cout << " is not inside. \n";
00455   }
00456 }
00457 std::cout << "\n";
00458 
00459 
00460 // Rank-3 array for 6 2D points and rank-2 array for the test result
00461 FieldContainer<double>  sixPoints(2, 3, 2);
00462 FieldContainer<int> testSixPoints(2, 3);
00463 
00464 sixPoints(0,0,0) = -1.0;   sixPoints(0,0,1) =  1.0;
00465 sixPoints(0,1,0) =  1.0;   sixPoints(0,1,1) =  0.0;
00466 sixPoints(0,2,0) =  0.0;   sixPoints(0,2,1) =  1.0;
00467 sixPoints(1,0,0) = -1.0;   sixPoints(1,0,1) = -1.0;
00468 sixPoints(1,1,0) =  0.1;   sixPoints(1,1,1) =  0.2;
00469 sixPoints(1,2,0) =  0.2;   sixPoints(1,2,1) =  0.3;
00470 
00471 std::cout <<"\t  Pointwise inclusion test for Triangle<6>: rank-3 array with six 2D points: \n";
00472 
00473 CellTools::checkPointwiseInclusion(testSixPoints, sixPoints, tri6);
00474 
00475 for(int i0 = 0; i0 < sixPoints.dimension(0); i0++){
00476   for(int i1 = 0; i1 < sixPoints.dimension(1); i1++) {
00477     std::cout << " point(" << i0 << "," << i1 << ") = (" 
00478     << std::setw(13) << std::right << sixPoints(i0, i1, 0) << "," 
00479     << std::setw(13) << std::right << sixPoints(i0, i1, 1) << ") ";
00480     if( testSixPoints(i0, i1) ) {
00481       std::cout << " is inside. \n";
00482     }
00483     else{
00484       std::cout << " is not inside. \n";
00485     }
00486     
00487   }
00488 }
00489 std::cout << "\n";
00490 
00491 
00492 // Rank-3 array for 6 3D reference points and rank-2 array for the test results
00493 FieldContainer<double> six3DPoints(2, 3, 3);
00494 FieldContainer<int> testSix3DPoints(2, 3);
00495 
00496 six3DPoints(0,0,0) = -1.0;   six3DPoints(0,0,1) =  1.0;   six3DPoints(0,0,2) =  1.0;
00497 six3DPoints(0,1,0) =  1.0;   six3DPoints(0,1,1) =  1.0;   six3DPoints(0,1,2) = -1.0;
00498 six3DPoints(0,2,0) =  0.0;   six3DPoints(0,2,1) =  1.1;   six3DPoints(0,2,2) =  1.0;
00499 six3DPoints(1,0,0) = -1.1;   six3DPoints(1,0,1) = -1.0;   six3DPoints(1,0,2) = -1.0;
00500 six3DPoints(1,1,0) =  0.1;   six3DPoints(1,1,1) =  0.2;   six3DPoints(1,1,2) =  0.2;
00501 six3DPoints(1,2,0) =  1.1;   six3DPoints(1,2,1) =  0.3;   six3DPoints(1,2,2) =  0.3;
00502 
00503 std::cout <<"\t  Pointwise inclusion test for Hexahedron<27>: rank-3 array with six 3D points: \n";
00504 
00505 CellTools::checkPointwiseInclusion(testSix3DPoints, six3DPoints, hex27);
00506 
00507 
00508 
00509 for(int i0 = 0; i0 < six3DPoints.dimension(0); i0++){
00510   for(int i1 = 0; i1 < six3DPoints.dimension(1); i1++) {
00511     std::cout << " point(" << i0 << "," << i1 << ") = (" 
00512     << std::setw(13) << std::right << six3DPoints(i0, i1, 0) << "," 
00513     << std::setw(13) << std::right << six3DPoints(i0, i1, 1) << "," 
00514     << std::setw(13) << std::right << six3DPoints(i0, i1, 2) << ") ";
00515     if( testSix3DPoints(i0, i1) ) {
00516       std::cout << " is inside. \n";
00517     }
00518     else{
00519       std::cout << " is not inside. \n";
00520     }
00521   }
00522 }
00523 
00524 
00525 std::cout << std::setprecision(6) << "\n" \
00526 << "===============================================================================\n"\
00527 << "| EXAMPLE 4-B: Pointwise inclusion test for a physical cell in cell workset   |\n"\
00528 << "===============================================================================\n";
00529 
00530 // Rank-2 array for a single set of 5 2D physical points and rank-1 array for the test results
00531 FieldContainer<double>  fivePoints(5,2);
00532 FieldContainer<int> testFivePoints(5);
00533 
00534 // These points will be tested for inclusion in the last Triangle cell specified by triNodes
00535 fivePoints(0, 0) = 2.1 ;   fivePoints(0, 1) = 1.0 ;       // in
00536 fivePoints(1, 0) = 3.0 ;   fivePoints(1, 1) = 0.75;       // in
00537 fivePoints(2, 0) = 3.5 ;   fivePoints(2, 1) = 1.9 ;       // out
00538 fivePoints(3, 0) = 2.5 ;   fivePoints(3, 1) = 1.0 ;       // in
00539 fivePoints(4, 0) = 2.75;   fivePoints(4, 1) = 2.0 ;       // out
00540 
00541 CellTools::checkPointwiseInclusion(testFivePoints, fivePoints, triNodes, triangle_3,  3);
00542 
00543 std::cout << " Vertices of Triangle #3: \n" 
00544 << "\t(" << triNodes(3, 0, 0) << ", " << triNodes(3, 0, 1) << ")\n"
00545 << "\t(" << triNodes(3, 1, 0) << ", " << triNodes(3, 1, 1) << ")\n"
00546 << "\t(" << triNodes(3, 1, 0) << ", " << triNodes(3, 1, 1) << ")\n"
00547 << " Inclusion test results for the physical points: \n\n";
00548 
00549 for(int i1 = 0; i1 < fivePoints.dimension(0); i1++) {
00550   std::cout << " point(" << i1 << ") = (" 
00551   << std::setw(13) << std::right << fivePoints(i1, 0) << "," 
00552   << std::setw(13) << std::right << fivePoints(i1, 1) << ") ";
00553   if( testFivePoints(i1) ) {
00554     std::cout << " is inside. \n";
00555   }
00556   else{
00557     std::cout << " is not inside. \n";
00558   }
00559 }
00560 std::cout << "\n";
00561 
00562 
00563 
00564 std::cout << std::setprecision(6) << "\n" \
00565 << "===============================================================================\n"\
00566 << "| EXAMPLE 4-C: Pointwise inclusion test for all physical cells in cell workset|\n"\
00567 << "===============================================================================\n";
00568 
00569 // Rank-3 array for 4 sets of 2 2D physical points and rank-2 array for the test results
00570 FieldContainer<double>  fourPointSets(4, 2, 2);
00571 FieldContainer<int> testFourSets(4, 2);
00572 
00573 // 1st point set - will be tested for inclusion in Triangle #0
00574 fourPointSets(0, 0, 0) = 0.25;     fourPointSets(0, 0, 1) = 0.75;       // in
00575 fourPointSets(0, 1, 0) = 0.00;     fourPointSets(0, 1, 1) =-0.01;       // out
00576 
00577 // 2nd point set - will be tested for inclusion in Triangle #1
00578 fourPointSets(1, 0, 0) = 1.50;     fourPointSets(1, 0, 1) = 1.50;       // in
00579 fourPointSets(1, 1, 0) = 0.99;     fourPointSets(1, 1, 1) = 1.50;       // out
00580 
00581 // 3rd point set - will be tested for inclusion in Triangle #2
00582 fourPointSets(2, 0, 0) =-0.25;     fourPointSets(2, 0, 1) = 0.70;       // in
00583 fourPointSets(2, 1, 0) = 0.0001;   fourPointSets(2, 1, 1) = 0.50;       // out
00584 
00585 // 4th point set - will be tested for inclusion in Triangle #3
00586 fourPointSets(3, 0, 0) = 3.00;     fourPointSets(3, 0, 1) = 1.00;       // in
00587 fourPointSets(3, 1, 0) = 3.50;     fourPointSets(3, 1, 1) = 0.50;       // out
00588 
00589 CellTools::checkPointwiseInclusion(testFourSets, fourPointSets, triNodes, triangle_3);
00590 
00591 for(int cell = 0; cell < triNodes.dimension(0); cell++){
00592 
00593   std::cout << " Testing point set inclusion for Triangle #" << cell << " with vertices \n" 
00594   << "\t(" << triNodes(cell, 0, 0) << ", " << triNodes(cell, 0, 1) << ")\n"
00595   << "\t(" << triNodes(cell, 1, 0) << ", " << triNodes(cell, 1, 1) << ")\n"
00596   << "\t(" << triNodes(cell, 1, 0) << ", " << triNodes(cell, 1, 1) << ")\n"
00597   << " Results for physical point set indexed by the cell ordinal  " << cell << " \n\n";
00598   
00599   for(int i1 = 0; i1 < fourPointSets.dimension(1); i1++) {
00600     std::cout << " point(" << i1 << ") = (" 
00601     << std::setw(13) << std::right << fourPointSets(cell, i1, 0) << "," 
00602     << std::setw(13) << std::right << fourPointSets(cell, i1, 1) << ") ";
00603     if( testFourSets(cell, i1) ) {
00604       std::cout << " is inside  Triangle #" << cell << " \n";
00605     }
00606     else{
00607       std::cout << " is outside Triangle #" << cell << " \n";
00608     }
00609   }
00610   if(cell < triNodes.dimension(0) - 1) {
00611     std::cout << "-------------------------------------------------------------------------------\n";
00612   }
00613   else{
00614     std::cout <<" \n"; 
00615   }
00616 }// cell
00617 
00618 
00619 
00620 std::cout << "\n" \
00621 << "===============================================================================\n"\
00622 << "| EXAMPLE 5: Using point set inclusion test method                            |\n"\
00623 << "===============================================================================\n";
00624 
00625 std::cout <<"\t  Point set inclusion test for Triangle<6>: rank-2 array with four 2D point: \n";
00626 
00627 if( CellTools::checkPointsetInclusion(fourPoints, tri6) ) {
00628   std::cout << "\t - All points are inside the reference Triangle<6> cell. \n\n ";
00629 }
00630 else{
00631   std::cout << "\t - At least one point is not inside the reference Triangle<6> cell. \n\n";
00632 }
00633 
00634 std::cout <<"\t  Point set inclusion test for Hexahedron<27>: rank-3 array with six 3D point: \n";
00635 
00636 if( CellTools::checkPointsetInclusion(six3DPoints, hex27) ) {
00637   std::cout << "\t - All points are inside the reference Hexahedron<27> cell. \n\n ";
00638 }
00639 else{
00640   std::cout << "\t - At least one point is not inside the reference Hexahedron<27> cell. \n\n";
00641 }
00642 
00643 
00644 
00645 std::cout << std::setprecision(4) << "\n" \
00646 << "===============================================================================\n"\
00647 << "| EXAMPLE 6: mapping a single point set to physical cells with base topology  |\n"\
00648 << "===============================================================================\n";
00649 
00650 // Rank-3 array with dimensions (P, D) for points on the reference triangle
00651 FieldContainer<double> refTriPoints(3,triangle_3.getDimension() );
00652 refTriPoints(0,0) = 0.2;  refTriPoints(0,1) = 0.0;      // on edge 0
00653 refTriPoints(1,0) = 0.4;  refTriPoints(1,1) = 0.6;      // on edge 1
00654 refTriPoints(2,0) = 0.0;  refTriPoints(2,1) = 0.8;      // on edge 2
00655 
00656 // Reference points will be mapped to physical cells with vertices in triNodes: define the appropriate array
00657 FieldContainer<double> physTriPoints(triNodes.dimension(0),      // cell count
00658                                      refTriPoints.dimension(0),     // point count
00659                                      refTriPoints.dimension(1));    // point dimension (=2)
00660 
00661 CellTools::mapToPhysicalFrame(physTriPoints, refTriPoints, triNodes, triangle_3);
00662 
00663 for(int cell = 0; cell < triNodes.dimension(0); cell++){
00664   std::cout << "====== Triangle " << cell << " ====== \n";
00665   for(int pt = 0; pt < refTriPoints.dimension(0); pt++){
00666     std::cout 
00667     <<  "(" 
00668     << std::setw(13) << std::right << refTriPoints(pt,0) << "," 
00669     << std::setw(13) << std::right << refTriPoints(pt,1) << ") -> "
00670     <<  "(" 
00671     << std::setw(13) << std::right << physTriPoints(cell, pt, 0) << "," 
00672     << std::setw(13) << std::right << physTriPoints(cell, pt, 1) << ") \n"; 
00673   }
00674   std::cout << "\n";
00675 }
00676 
00677 
00678 std::cout << std::setprecision(4) << "\n" \
00679 << "===============================================================================\n"\
00680 << "| EXAMPLE 7: mapping a single point set to physical cells with ext. topology  |\n"\
00681 << "===============================================================================\n";
00682 /*
00683  *  This example illustrates reference-to-physical mapping for a triangle with curved sides. The
00684  *  physical curved triangle is defined as follows:
00685  *
00686  *    - edge 0: starts at (1, -1/2) ends at (1,1)     and is a vertical line segment
00687  *    - edge 1: starts at (1,1)     ends at (0,0)     and is parabolic segment lying on y=x^2
00688  *    - edge 2: starts at (0,0)     ends at (1,-1/2)  and is parabolic segment lying on y = -1/2 x^2
00689  * 
00690  *  The triangle is uniquely specified by its 3 vertices and 3 edge points. Therefore, to compute the
00691  *  map from reference to physical coordinates we specify Triangle<6> topology.
00692  */
00693 
00694 // A (C,V,D) array with the 6 nodes of the physical Triangle<6>
00695 FieldContainer<double> tri6Nodes(1,6,2);
00696 tri6Nodes(0,0,0) = 1.0;    tri6Nodes(0,0,1) = -0.5;
00697 tri6Nodes(0,1,0) = 1.0;    tri6Nodes(0,1,1) =  1.0;
00698 tri6Nodes(0,2,0) = 0.0;    tri6Nodes(0,2,1) =  0.0;
00699 
00700 tri6Nodes(0,3,0) = 1.0;    tri6Nodes(0,3,1) =  0.0;
00701 tri6Nodes(0,4,0) = 0.5;    tri6Nodes(0,4,1) =  0.25;
00702 tri6Nodes(0,5,0) = 0.5;    tri6Nodes(0,5,1) = -0.125;
00703 
00704 // A (P, D) array with the 6 nodes of the reference Triangle<6> plus some extra points
00705 FieldContainer<double> refTri6Points(9,2);
00706 refTri6Points(0,0) = 0.0;    refTri6Points(0,1) =  0.0;
00707 refTri6Points(1,0) = 1.0;    refTri6Points(1,1) =  0.0;
00708 refTri6Points(2,0) = 0.0;    refTri6Points(2,1) =  1.0;
00709 
00710 refTri6Points(3,0) = 0.5;    refTri6Points(3,1) =  0.0;
00711 refTri6Points(4,0) = 0.5;    refTri6Points(4,1) =  0.5;
00712 refTri6Points(5,0) = 0.0;    refTri6Points(5,1) =  0.5;
00713 
00714 refTri6Points(6,0) = 0.75;   refTri6Points(6,1) =  0.0;               // on ref edge 0
00715 refTri6Points(7,0) = 0.25;   refTri6Points(7,1) =  0.75;              // on ref edge 1
00716 refTri6Points(8,0) = 0.00;   refTri6Points(8,1) =  0.25;              // on ref edge 2
00717 
00718 // A (C,P,D) array for the images of the reference points in physical frame
00719 FieldContainer<double> physTri6Points(tri6Nodes.dimension(0),             // cell count
00720                                       refTri6Points.dimension(0),         // point count
00721                                       refTri6Points.dimension(1));        // point dimension (=2)
00722 
00723 // Define the cell topology: use the extended Triangle<6> topology to fit the curved edges
00724 CellTopology triangle_6(shards::getCellTopologyData<Triangle<6> >() );
00725 
00726 // 
00727 CellTools::mapToPhysicalFrame(physTri6Points, refTri6Points, tri6Nodes, triangle_6);
00728 
00729 for(int cell = 0; cell < tri6Nodes.dimension(0); cell++){
00730   std::cout << "====== Triangle " << cell << " ====== \n";
00731   for(int pt = 0; pt < refTri6Points.dimension(0); pt++){
00732     std::cout 
00733     <<  "(" 
00734     << std::setw(13) << std::right << refTri6Points(pt,0) << "," 
00735     << std::setw(13) << std::right << refTri6Points(pt,1) << ") -> "
00736     <<  "(" 
00737     << std::setw(13) << std::right << physTri6Points(cell, pt, 0) << "," 
00738     << std::setw(13) << std::right << physTri6Points(cell, pt, 1) << ") \n"; 
00739   }
00740   std::cout << "\n";
00741 }
00742 
00743 
00744 
00745 std::cout << "\n" \
00746 << "===============================================================================\n"\
00747 << "| EXAMPLE 8: mapping a single physical point set to reference frame           |\n"\
00748 << "===============================================================================\n";
00749 /*
00750  * This example shows use of mapToReferenceFrame with rank-2 array (P,D) of points in physical frame.
00751  * Points are mapped back to reference frame using the mapping for one of the physical cells whose
00752  * nodes are passed as an argument. Therefore, this use case requires a valid cell ordinal (relative
00753  * to the nodes array) to be specified. The output array for the images of the points is also rank-2 (P,D).
00754  *
00755  */
00756 // Rank-2 arrays with dimensions (P, D) for physical points and their preimages
00757 FieldContainer<double> physPoints(5, triangle_3.getDimension() ); 
00758 FieldContainer<double> preImages (5, triangle_3.getDimension() ); 
00759 
00760 // First 3 points are the vertices of the last triangle (ordinal = 3) stored in triNodes
00761 physPoints(0,0) = triNodes(3, 0, 0);   physPoints(0,1) = triNodes(3, 0, 1) ;
00762 physPoints(1,0) = triNodes(3, 1, 0);   physPoints(1,1) = triNodes(3, 1, 1);
00763 physPoints(2,0) = triNodes(3, 2, 0);   physPoints(2,1) = triNodes(3, 2, 1);
00764 
00765 // last 2 points are just some arbitrary points contained in the last triangle
00766 physPoints(3,0) = 3.0;                    physPoints(3,1) = 1.0;
00767 physPoints(4,0) = 2.5;                    physPoints(4,1) = 1.1;
00768 
00769 
00770 // Map physical points from triangle 3 to the reference frame
00771 int whichCell = 3;
00772 CellTools::mapToReferenceFrame(preImages, physPoints, triNodes, triangle_3, whichCell  );
00773 
00774 std::cout << " Mapping from Triangle #"<< whichCell << " with vertices: \n" 
00775 << "\t(" << triNodes(whichCell, 0, 0) << ", " << triNodes(whichCell, 0, 1) << ")\n"
00776 << "\t(" << triNodes(whichCell, 1, 0) << ", " << triNodes(whichCell, 1, 1) << ")\n"
00777 << "\t(" << triNodes(whichCell, 1, 0) << ", " << triNodes(whichCell, 1, 1) << ")\n\n"
00778 << " Physical points and their reference cell preimages: \n";
00779 
00780 for(int pt = 0; pt < physPoints.dimension(0); pt++){
00781   std::cout 
00782   <<  "(" << std::setw(13) << std::right << physPoints(pt,0) << "," 
00783   << std::setw(13) << std::right << physPoints(pt,1) << ") -> "
00784   <<  "(" 
00785   << std::setw(13) << std::right << preImages(pt, 0) << "," 
00786   << std::setw(13) << std::right << preImages(pt, 1) << ") \n"; 
00787 }
00788 std::cout << "\n";
00789 
00790 // As a check, map pre-images back to Triangle #3
00791 FieldContainer<double> images(5, triangle_3.getDimension() );
00792 CellTools::mapToPhysicalFrame(images, preImages, triNodes, triangle_3, whichCell);
00793 
00794 std::cout << " Check: map preimages back to Triangle #3: \n"; 
00795 for(int pt = 0; pt < images.dimension(0); pt++){
00796   std::cout 
00797   <<  "(" << std::setw(13) << std::right << preImages(pt,0) << "," 
00798   << std::setw(13) << std::right << preImages(pt,1) << ") -> "
00799   <<  "(" 
00800   << std::setw(13) << std::right << images(pt, 0) << "," 
00801   << std::setw(13) << std::right << images(pt, 1) << ") \n"; 
00802 }
00803 std::cout << "\n";
00804 
00805 
00806 std::cout << "\n" \
00807 << "===============================================================================\n"\
00808 << "| EXAMPLE 9: mapping physical point sets on a cell workset to reference frame |\n"\
00809 << "===============================================================================\n";
00810 /*
00811  * This example shows use of mapToReferenceFrame with rank-3 array (C,P,D) of points in physical frame.
00812  * For each cell ordinal (relative to the nodes array), the associated point set is mapped back to 
00813  * reference frame using the mapping corresponding to the physical cell with that ordinal. This use case
00814  * requires the default value -1 for the cell ordinal. The output array for the images of the points 
00815  * is also rank-3 (C,P,D).
00816  *
00817  */
00818 // Rank-3 arrays with dimensions (C, P, D) for physical points and their preimages
00819 FieldContainer<double> physPointSets(triNodes.dimension(0), 2, triangle_3.getDimension() ); 
00820 FieldContainer<double> preImageSets (triNodes.dimension(0), 2, triangle_3.getDimension() ); 
00821 
00822 // Point set on Triangle #0
00823 physPointSets(0,0,0) = 0.25; physPointSets(0,0,1) = 0.25; 
00824 physPointSets(0,1,0) = 0.50; physPointSets(0,1,1) = 0.50; 
00825 
00826 // Point set on Triangle #1
00827 physPointSets(1,0,0) = 1.00; physPointSets(1,0,1) = 1.00; 
00828 physPointSets(1,1,0) = 1.25; physPointSets(1,1,1) = 1.75; 
00829 
00830 // Point set on Triangle #2
00831 physPointSets(2,0,0) =-0.25; physPointSets(2,0,1) = 0.25; 
00832 physPointSets(2,1,0) =-0.50; physPointSets(2,1,1) = 0.50; 
00833 
00834 // Point set on Triangle #3
00835 physPointSets(3,0,0) = 3.00; physPointSets(3,0,1) = 1.00; 
00836 physPointSets(3,1,0) = 2.00; physPointSets(3,1,1) = 1.00; 
00837 
00838 
00839 // Map physical point sets to reference frame: requires default value of cell ordinal (skip last arg)
00840 CellTools::mapToReferenceFrame(preImageSets, physPointSets, triNodes, triangle_3);
00841 
00842 for(int cell = 0; cell < triNodes.dimension(0); cell++){
00843   std::cout << "====== Triangle " << cell << " ====== \n";
00844   for(int pt = 0; pt < physPointSets.dimension(1); pt++){
00845     std::cout 
00846     <<  "(" 
00847     << std::setw(13) << std::right << physPointSets(cell, pt, 0) << "," 
00848     << std::setw(13) << std::right << physPointSets(cell, pt, 1) << ") -> "
00849     <<  "(" 
00850     << std::setw(13) << std::right << preImageSets(cell, pt, 0) << "," 
00851     << std::setw(13) << std::right << preImageSets(cell, pt, 1) << ") \n"; 
00852   }
00853   std::cout << "\n";
00854 }
00855   
00856 // As a check, map preImageSets back to physical frame
00857 FieldContainer<double> postImageSets(triNodes.dimension(0), 2, triangle_3.getDimension() );
00858 
00859 CellTools::mapToPhysicalFrame(postImageSets, preImageSets, triNodes, triangle_3);
00860 
00861 std::cout << " Check: map preimages back to Triangles: \n"; 
00862 for(int cell = 0; cell < triNodes.dimension(0); cell++){
00863   std::cout << "====== Triangle " << cell << " ====== \n";
00864   for(int pt = 0; pt < preImageSets.dimension(1); pt++){
00865     std::cout 
00866     <<  "(" 
00867     << std::setw(13) << std::right << preImageSets(cell, pt, 0) << "," 
00868     << std::setw(13) << std::right << preImageSets(cell, pt, 1) << ") -> "
00869     <<  "(" 
00870     << std::setw(13) << std::right << postImageSets(cell, pt, 0) << "," 
00871     << std::setw(13) << std::right << postImageSets(cell, pt, 1) << ") \n"; 
00872   }
00873   std::cout << "\n";
00874 }
00875 
00876 return 0;
00877 }
00878 
00879 
00880 
00881 
00882 
00883 
00884 
00885 
00886 
00887 
00888 
00889 
00890 
00891 
00892 
00893 
00894 
00895 
00896 
00897 
00898 
00899 
00900 
00901 
00902 
00903 
00904 
00905