src/MEDMEM/MEDMEM_Extractor.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   2 //
   3 //  This library is free software; you can redistribute it and/or
   4 //  modify it under the terms of the GNU Lesser General Public
   5 //  License as published by the Free Software Foundation; either
   6 //  version 2.1 of the License.
   7 //
   8 //  This library is distributed in the hope that it will be useful,
   9 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 //  Lesser General Public License for more details.
  12 //
  13 //  You should have received a copy of the GNU Lesser General Public
  14 //  License along with this library; if not, write to the Free Software
  15 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // File      : MEDMEM_Extractor.cxx
  20 // Created   : Thu Dec 18 11:10:11 2008
  21 // Author    : Edward AGAPOV (eap)
  22
  23 #include "MEDMEM_Extractor.hxx"
  24
  25 #include <MEDMEM_Field.hxx>
  26 #include <MEDMEM_Mesh.hxx>
  27 #include <MEDMEM_Meshing.hxx>
  28 #include <MEDMEM_Support.hxx>
  29
  30 #include <list>
  31
  32 #include <math.h>
  33
  34 using namespace MED_EN;
  35 using namespace std;
  36 using namespace MEDMEM;
  37
  38 namespace { // local tools
  39
  40   const int _POLYGON = -1; //!< map key to store connectivity of polygons
  41
  42   const double _TOLER = 1e-12;
  43
  44   //================================================================================
  45   /*!
  46    * \brief calculate cross product of two vectors
  47    */
  48   void crossProduct(const double* v1, const double* v2, double* res)
  49   {
  50     res[0] = v1[1] * v2[2] - v1[2] * v2[1];
  51     res[1] = v1[2] * v2[0] - v1[0] * v2[2];
  52     res[2] = v1[0] * v2[1] - v1[1] * v2[0];
  53   }
  54   //================================================================================
  55   /*!
  56    * \brief calculate dot product of two vectors
  57    */
  58   double dotProduct(const double* v1, const double* v2)
  59   {
  60     return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
  61   }
  62   //================================================================================
  63   /*!
  64    * \brief SUPPORT owning the mesh
  65    */
  66   struct TSupport : public SUPPORT
  67   {
  68     TSupport(MESH* mesh): SUPPORT(mesh) {}
  69     ~TSupport() { delete getMesh(); setMeshDirectly((MESH*)0); }
  70   };
  71   //================================================================================
  72   /*!
  73    * \brief Field owning the support
  74    */
  75   struct TField: public FIELD<double>
  76   {
  77     TField( const SUPPORT * Support, int nbComp ): FIELD<double>( Support, nbComp ) {}
  78     ~TField() { if ( _support ) delete _support; _support = (SUPPORT*)0; }
  79   };
  80   //================================================================================
  81   /*!
  82    * \brief Accessor to some ids. Provides operator[] and more-next access methods
  83    */
  84   class TIter
  85   {
  86     const int *_cur, *_end;
  87   public:
  88     TIter(const int* start, const int* end):_cur(start), _end(end) {}
  89     TIter(const int* start, int nb):_cur(start), _end(start+nb) {}
  90     int size() const { return _end - _cur; }
  91     int operator[] (int i) const { return _cur[i]; }
  92     bool more() const { return _cur < _end; }
  93     int next() { return *_cur++; }
  94     const int* begin() const { return _cur; }
  95     const int* end() const { return _end; }
  96   };
  97   //================================================================================
  98   /*!
  99    * \brief Edge linking two nodes, used to find equal edges and store edges in std containers
 100    */
 101   struct TEdge: public pair<int,int>
 102   {
 103     TEdge(const int n1=0, const int n2=0): pair<int,int>(n1,n2)
 104     { if ( n2 < n1 ) first = n2, second = n1; }
 105     TEdge(const TIter& nodes, int i )
 106     { *this = TEdge( nodes[i], nodes[ (i+1)%nodes.size() ]); }
 107     int node1() const { return first; }
 108     int node2() const { return second; }
 109   };
 110   //================================================================================
 111   /*!
 112    * \brief Tool providing iteration on edges of the cell of given classical type
 113    */
 114   struct TEdgeIterator
 115   {
 116     vector<TEdge>* _edges;
 117     TEdgeIterator(const medGeometryElement type);
 118     int getNbEdges() const { return _edges->size(); }
 119     TEdge getEdge(int i, const int* cellConn ) const
 120     { return TEdge( cellConn[(*_edges)[i].node1()], cellConn[(*_edges)[i].node2()]); }
 121     TEdge getEdge(int i, const TIter& cellNodes ) const
 122     { return TEdge( cellNodes[(*_edges)[i].node1()], cellNodes[(*_edges)[i].node2()]); }
 123   };
 124   //================================================================================
 125   /*!
 126    * \brief Comparator used to sort nodes of polygon
 127    */
 128   struct TNodeCompare
 129   {
 130     const double* _coords; //!< coordinates of mesh nodes in full interlace
 131     const double* _bc; //!< polygon barycentre
 132     int _i1, _i2; //!< indices of two ordinates used to compare points
 133
 134     TNodeCompare(const double* nodeCoords, int i1, int i2):
 135       _coords(nodeCoords),_i1(i1),_i2(i2) {}
 136
 137     void setBaryCenter(const double* bc) { _bc = bc; }
 138
 139     bool operator()(const int pt1, const int pt2) {
 140       const double* p1 = _coords + 3*(pt1-1);
 141       const double* p2 = _coords + 3*(pt2-1);
 142       // calculate angles of Op1 and Op2 with the i1 axis on plane (i1,i2)
 143       double ang1 = atan2(p1[_i1] - _bc[0], p1[_i2] - _bc[1]);
 144       double ang2 = atan2(p2[_i1] - _bc[0], p2[_i2] - _bc[1]);
 145       return ang1 > ang2;
 146     }
 147   };
 148   //================================================================================
 149   /*!
 150    * \brief Return tolerance to consider nodes coincident
 151    */
 152   double getTolerance(const MESH* mesh )
 153   {
 154     vector< vector<double> > bb = mesh->getBoundingBox();
 155     double diagonal = 0;
 156     for ( int j = 0; j < mesh->getSpaceDimension(); ++j ) {
 157       double dist = bb[0][j] - bb[1][j];
 158       diagonal += dist*dist;
 159     }
 160     return sqrt( diagonal ) * 1e-7;
 161   }
 162   //================================================================================
 163   /*!
 164    * \brief Line data and some methods
 165    */
 166   struct TLine
 167   {
 168     const double* _dir;
 169     const double* _coord;
 170     int           _maxDir; //!< index of maximal coordinate of _dir
 171     TLine( const double* dir, const double* coord): _dir(dir), _coord(coord) {
 172       _maxDir = 0;
 173       if ( fabs( _dir[_maxDir] ) < fabs( _dir[1] )) _maxDir = 1;
 174       if ( fabs( _dir[_maxDir] ) < fabs( _dir[2] )) _maxDir = 2;
 175     }
 176     //!< Check if intersection points coincide in _maxDir direction
 177     bool isSame( const double* p1, const double* p2 ) const {
 178       return fabs(p1[_maxDir] - p2[_maxDir]) < _TOLER;
 179     }
 180   };
 181   //================================================================================
 182   /*!
 183    * \brief Result of transfixing a face
 184    *
 185    * TIntersection's make up a chain via _prevSection field. Depending on whether
 186    * _prevSection is NULL or not, there are two types of TIntersection:
 187    * . _prevSection == NULL: TIntersection ending the chain. It is TIntersection from
 188    * which chain building starts.
 189    * . _prevSection != NULL: intermediate TIntersection, stores cells cut by a result segment.
 190    */
 191   struct TIntersection
 192   {
 193     double   _point[3];   //!< coordinates of itersection
 194     set<int> _cells;      //!< intersected cells
 195     int      _face;       //!< intersected face of a cell
 196     set<int> _startNodes; //!< nodes around which to look for next intersection
 197     TIntersection* _prevSection; //!< neighbor intersection
 198
 199     TIntersection(): _face(-1), _prevSection(NULL)
 200     {}
 201     ~TIntersection() {
 202       if ( _prevSection ) delete _prevSection; _prevSection=0;
 203     }
 204     void getRange( double& min, double& max, const int j ) const {
 205       if ( _point[j] < min ) min = _point[j];
 206       if ( _point[j] > max ) max = _point[j];
 207       if ( _prevSection ) _prevSection->getRange( min, max, j );
 208     }
 209     void reverse() {
 210       if ( _prevSection ) {
 211         _prevSection->reverse();
 212         _prevSection->_cells = _cells;
 213         _prevSection->_prevSection = this;
 214         _prevSection = NULL;
 215       }
 216     }
 217     int size() const { return 1 + ( _prevSection ? _prevSection->size() : 0 ); }
 218   };
 219   //================================================================================
 220   /*!
 221    * \brief Provider of comfortable access to connectivity data of MESH
 222    */
 223   struct TMeshData
 224   {
 225     int           _dim;
 226     double        _tolerance;
 227     const double* _coord;
 228     const int*    _cellConn;
 229     const int*    _cellConnIndex;
 230     const int*    _faceConn;
 231     const int*    _faceConnIndex;
 232     const int*    _face2Cell;
 233     const int*    _face2CellIndex;
 234     const int*    _cell2Face;
 235     const int*    _cell2FaceIndex;
 236     const int*    _node2Cell;
 237     const int*    _node2CellIndex;
 238     TMeshData(const MESH &mesh)
 239     {
 240       _tolerance      = getTolerance(&mesh);
 241       _dim            = mesh.getSpaceDimension();
 242       _coord          = mesh.getCoordinates(MED_FULL_INTERLACE);
 243       _cellConn       = mesh.getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
 244                                             MED_CELL, MED_ALL_ELEMENTS);
 245       _cellConnIndex  = mesh.getConnectivityIndex(MED_NODAL, MED_CELL);
 246       _cell2Face      = mesh.getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING,
 247                                             MED_CELL, MED_ALL_ELEMENTS);
 248       _cell2FaceIndex = mesh.getConnectivityIndex( MED_DESCENDING, MED_CELL );
 249       _face2Cell      = mesh.getReverseConnectivity( MED_DESCENDING );
 250       _face2CellIndex = mesh.getReverseConnectivityIndex( MED_DESCENDING );
 251       _faceConn       = mesh.getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
 252                                             MED_FACE, MED_ALL_ELEMENTS);
 253       _faceConnIndex  = mesh.getConnectivityIndex(MED_NODAL, MED_FACE);
 254       _node2Cell      = mesh.getReverseConnectivity( MED_NODAL );
 255       _node2CellIndex = mesh.getReverseConnectivityIndex( MED_NODAL );
 256     }
 257     double tolerance() const {
 258       return _tolerance; }
 259     const double* getNodeCoord( int node ) const {
 260       return _coord + _dim*(node-1); }
 261     TIter getFaces( int cell ) const {
 262       return TIter( _cell2Face+_cell2FaceIndex[cell-1]-1, _cell2Face+_cell2FaceIndex[cell]-1 ); }
 263     TIter getCellNodes( int cell ) const {
 264       return TIter( _cellConn+_cellConnIndex[cell-1]-1, _cellConn+_cellConnIndex[cell]-1 ); }
 265     TIter getFaceNodes( int face ) const {
 266       face = abs(face);
 267       return TIter( _faceConn+_faceConnIndex[face-1]-1, _faceConn+_faceConnIndex[face]-1 ); }
 268     TIter getCellsByNode( int node ) const {
 269       return TIter( _node2Cell+_node2CellIndex[node-1]-1, _node2Cell+_node2CellIndex[node]-1 ); }
 270     TIter getCellsByFace( int face ) const {
 271       face = abs(face);
 272       return TIter( _face2Cell+_face2CellIndex[face-1]-1, _face2Cell+_face2CellIndex[face]-1 ); }
 273     int isFreeFace( int face ) const {
 274       TIter cells = getCellsByFace( face );
 275       return ( cells[1] == 0 ) ? cells[0] : 0; }
 276   };
 277   //================================================================================
 278   /*!
 279    * \brief Calculates face normal
 280    *  \retval bool - false if face has zero area
 281    */
 282   bool calcFaceNormal( const int        face,
 283                        const TMeshData& meshData,
 284                        double*          norm)
 285   {
 286     TIter nodes = meshData.getFaceNodes( face );
 287     bool zeroArea = false;
 288     int i = 0;
 289     do {
 290       const double* p1 = meshData.getNodeCoord( nodes[i] );
 291       const double* p2 = meshData.getNodeCoord( nodes[i+1] );
 292       const double* p3 = meshData.getNodeCoord( nodes[i+2] );
 293       double p2p1[3] = { p1[0]-p2[0], p1[1]-p2[1], p1[2]-p2[2] };
 294       double p2p3[3] = { p3[0]-p2[0], p3[1]-p2[1], p3[2]-p2[2] };
 295       crossProduct( p2p3, p2p1, norm );
 296       double normSize2 = dotProduct( norm, norm );
 297       zeroArea = (normSize2 < DBL_MIN);
 298     }
 299     while ( zeroArea && i < 2 );
 300     return zeroArea;
 301   }
 302   //================================================================================
 303   /*!
 304    * \brief Fill set of cells sharing edge
 305    */
 306   void getCellsSharingEdge( const TEdge& edge, const TMeshData& meshData, set<int> & theCells )
 307   {
 308     TIter cells = meshData.getCellsByNode( edge.node1() );
 309     while ( cells.more() ) {
 310       int cell = cells.next();
 311       TIter nodes = meshData.getCellNodes( cell );
 312       TEdgeIterator edgeIter( medGeometryElement( 300 + nodes.size() ));
 313       for ( int i = 0, nb = edgeIter.getNbEdges(); i < nb; ++i ) {
 314         TEdge e = edgeIter.getEdge( i , nodes );
 315         if ( edge == e ) {
 316           theCells.insert( cell );
 317           break;
 318         }
 319       }
 320     }
 321   }
 322   bool canIntersect(const int        cell,
 323                     const TMeshData& meshData,
 324                     const TLine&     line);
 325
 326   TIntersection* intersect(const int        cell,
 327                            const TMeshData& meshData,
 328                            const TLine&     line,
 329                            set<int>&        checkedCells,
 330                            TIntersection*   prevInter=0);
 331
 332 } // noname namespace
 333
 334 namespace MEDMEM
 335 {
 336
 337   //================================================================================
 338   /*!
 339    * \brief Creates a tool
 340    *  \param inputField - input field
 341    *
 342    * The input field is supposed to comply with following conditions <ul>
 343    *<li>  it is constant by element (i.e. has 1 gauss point),</li>
 344    *<li>  it's support mesh does not contain poly elements,</li>
 345    *<li>  volumic elements have planar faces,</li>
 346    *<li>  surfasic elements have linear edges.</li></ul>
 347    */
 348   //================================================================================
 349
 350 Extractor::Extractor(const FIELD<double>& inputField) throw (MEDEXCEPTION)
 351 : _myInputField( & inputField )
 352 {
 353   const char* LOC = "Extractor::Extractor(inputField) :";
 354
 355   // Check if the input field complies with the conditions
 356
 357   if ( !inputField.getSupport() )
 358     throw MEDEXCEPTION(STRING(LOC) << "InputField has NULL support");
 359
 360   medEntityMesh entity = inputField.getSupport()->getEntity();
 361   if ( entity == MED_NODE || entity == MED_EDGE )
 362     throw MEDEXCEPTION(STRING(LOC) << "InputField has invalid supporting entity");
 363
 364   if ( inputField.getSupport()->getNumberOfElements(MED_ALL_ELEMENTS) == 0 )
 365     throw MEDEXCEPTION(STRING(LOC) << "InputField has support of zero size");
 366
 367   if ( inputField.getGaussPresence() && inputField.getNumberOfGaussPoints()[0] > 1 )
 368     throw MEDEXCEPTION(STRING(LOC) << "InputField is not constant be element");
 369
 370   MESH* mesh = inputField.getSupport()->getMesh();
 371   if ( !mesh )
 372     throw MEDEXCEPTION(STRING(LOC) << "InputField has support with NULL mesh");
 373
 374   if ( mesh->getSpaceDimension() < 2 )
 375     throw MEDEXCEPTION(STRING(LOC) << "InputField with 1D support not acceptable");
 376
 377   if ( mesh->getNumberOfPolygons() > 0 ||
 378        mesh->getNumberOfPolyhedron() > 0 )
 379     throw MEDEXCEPTION(STRING(LOC) << "InputField has supporting mesh with poly elements");
 380
 381   if ( mesh->getConnectivityptr()->getEntityDimension() < 2 )
 382     throw MEDEXCEPTION(STRING(LOC) << "Invalid entity dimension of connectivity");
 383 }
 384
 385 //================================================================================
 386   /*!
 387    * \brief Creates a field by cutting inputField by a plane
 388    *  \param coords - give a point to pass through by the plane
 389    *  \param normal - gives the plane normal
 390    *  \retval FIELD<double>* - resulting field holding ownership of its support,
 391    *                           which in its turn holds ownership of its mesh
 392    *
 393    * If the plane does not intersect the field, NULL is returned.
 394    */
 395 //================================================================================
 396
 397 FIELD<double>* Extractor::extractPlane(const double* coords, const double* normal)
 398   throw (MEDEXCEPTION)
 399 {
 400   const char* LOC = "Extractor::extractPlane(const double* coords, const double* normal) :";
 401
 402   // check agrument validity
 403   if ( !coords || !normal )
 404     throw MEDEXCEPTION(STRING(LOC) << "NULL argument");
 405
 406   double normalSize = sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
 407   if ( normalSize <= DBL_MIN )
 408     throw MEDEXCEPTION(STRING(LOC) << "normal has zero size");
 409
 410   if ( _myInputField->getSupport()->getMesh()->getSpaceDimension() < 3 )
 411     throw MEDEXCEPTION(STRING(LOC) << "Extraction by plane is possible in 3D space only ");
 412
 413   double norm[3] = { normal[0] / normalSize, normal[1] / normalSize, normal[2] / normalSize };
 414
 415   // cut mesh
 416   map<int,set<int> > new2oldCells;
 417   MESH* mesh = divideEdges( coords, norm, new2oldCells );
 418   if ( !mesh ) return 0;
 419
 420   return makeField( new2oldCells, mesh );
 421 }
 422
 423 //================================================================================
 424   /*!
 425    * \brief Creates a field by cutting inputField by a line
 426    *  \param coords - give a point to pass through by the line
 427    *  \param direction - gives a vector collinear to the line
 428    *  \retval FIELD<double>* - resulting field holding ownership of its support,
 429    *                           which in its turn holds ownership of its mesh
 430    *
 431    * If the line does not intersect the field, NULL is returned.
 432    */
 433 //================================================================================
 434
 435 FIELD<double>* Extractor::extractLine(const double* coords, const double* direction)
 436   throw (MEDEXCEPTION)
 437 {
 438   const char* LOC = "Extractor::extractLine(const double* coords, const double* direction) :";
 439
 440   // check agrument validity
 441   if ( !coords || !direction )
 442     throw MEDEXCEPTION(STRING(LOC) << "NULL argument");
 443
 444   double directionSize =
 445     sqrt(direction[0]*direction[0] + direction[1]*direction[1] + direction[2]*direction[2]);
 446   if ( directionSize <= DBL_MIN )
 447     throw MEDEXCEPTION(STRING(LOC) << "direction has zero size");
 448
 449   const SUPPORT* support = _myInputField->getSupport();
 450   const medGeometryElement* inTypes = support->getTypes();
 451   const int meshDim = inTypes[ support->getNumberOfTypes()-1 ] / 100;
 452
 453   if ( meshDim == 2 && support->getMesh()->getSpaceDimension() == 3 )
 454     throw MEDEXCEPTION(STRING(LOC) << "Extraction from 2D mesh not supported");
 455
 456   map<int,set<int> > new2oldCells;
 457   MESH* mesh = 0;
 458   if ( meshDim == 2 )
 459   {
 460     double norm[2] = { direction[1] / directionSize,
 461                        direction[0] / directionSize,  };
 462     // cut mesh
 463     mesh = divideEdges( coords, norm, new2oldCells );
 464   }
 465   else
 466   {
 467     double dir[3] = { direction[0] / directionSize,
 468                       direction[1] / directionSize,
 469                       direction[2] / directionSize };
 470     mesh = transfixFaces( coords, dir, new2oldCells );
 471   }
 472
 473   if ( !mesh ) return 0;
 474
 475   return makeField( new2oldCells, mesh );
 476 }
 477
 478 //================================================================================
 479 /*!
 480  * \brief Creates a new field and fill it from the input one
 481  *  \param new2oldCells - mapping of new to old cells
 482  *  \retval FIELD<double>* - resulting field
 483  */
 484 //================================================================================
 485
 486 FIELD<double>* Extractor::makeField( const map<int,set<int> >& new2oldCells,
 487                                      MESH*                     mesh ) const
 488 {
 489   // make new field
 490   int nbComp               = _myInputField->getNumberOfComponents();
 491   FIELD<double> * outField = new TField( new TSupport( mesh ), nbComp );
 492   double* outValues        = const_cast<double*>( outField->getValue() );
 493
 494   outField->setComponentsNames       ( _myInputField->getComponentsNames() );
 495   outField->setName                  ( STRING("Extracted from ")<< _myInputField->getName() );
 496   outField->setDescription           ( STRING("Created by MEDMEM::Extractor"));
 497   outField->setComponentsDescriptions( _myInputField->getComponentsDescriptions() );
 498   outField->setMEDComponentsUnits    ( _myInputField->getMEDComponentsUnits() );
 499
 500
 501   // put values to the new field
 502   map<int,set<int> >::const_iterator new_olds, noEnd = new2oldCells.end();
 503   for ( new_olds = new2oldCells.begin(); new_olds != noEnd; ++new_olds )
 504   {
 505     for ( int j = 0; j < nbComp; ++j )
 506     {
 507       int ind = ( new_olds->first - 1 ) * nbComp + j;
 508       outValues[ ind ] = 0.0;
 509       set<int>::const_iterator inCells = new_olds->second.begin();
 510       set<int>::const_iterator inEnd   = new_olds->second.end();
 511       for ( ; inCells != inEnd; ++inCells )
 512         outValues[ ind ] += _myInputField->getValueIJ( *inCells, j+1 );
 513
 514       outValues[ ind ] /= new_olds->second.size();
 515     }
 516   }
 517   return outField;
 518 }
 519
 520 //================================================================================
 521 /*!
 522  * \brief Makes a mesh by dividing edges of cells of the input mesh by plane
 523  *        in 3D or by line in 2D.
 524  *  \param coords - give a point to pass through by the plane or the line
 525  *  \param normal - gives the normal to plane or line
 526  *  \param new2oldCells - output map of new cells to cut input cell
 527  */
 528 //================================================================================
 529
 530 MESH* Extractor::divideEdges(const double*       coords,
 531                              const double*       normal,
 532                              map<int,set<int> >& new2oldCells)
 533 {
 534   const char* LOC = "Extractor::divideEdges()";
 535
 536   const SUPPORT* support            = _myInputField->getSupport();
 537   medEntityMesh entity              = support->getEntity();
 538   MESH* inMesh                      = support->getMesh();
 539   const medGeometryElement* inTypes = support->getTypes();
 540
 541   const int* inConn      = inMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
 542                                                    entity, MED_ALL_ELEMENTS);
 543   const int* inConnIndex = inMesh->getConnectivityIndex(MED_NODAL, entity);
 544   const int spaceDim     = inMesh->getSpaceDimension();
 545   const int meshDim      = inTypes[ support->getNumberOfTypes()-1 ] / 100;
 546
 547
 548   // connectivity of new cells by nb of nodes per cell
 549   map< int, vector< int > > newConnByNbNodes;
 550   int nbInputCells = support->getNumberOfElements( MED_ALL_ELEMENTS );
 551   int minNbNodesPerCell = 2, maxNbNodesPerCell = 2;
 552   if ( meshDim == 3 ) {
 553     newConnByNbNodes[ 3 ].reserve( nbInputCells/2 );
 554     newConnByNbNodes[ 4 ].reserve( nbInputCells/2 );
 555     minNbNodesPerCell = 3;
 556     maxNbNodesPerCell = int (MED_POLYGON); // polygones allowed
 557   }
 558   else {
 559     newConnByNbNodes[ 2 ].reserve( nbInputCells/2 );
 560   }
 561   list<int> nbNodesPerPolygon; // to make connectivity index of polygons
 562
 563   // new cells
 564   map< set<int>, int> oldNodes2newCell; //!< map connectivity of all old nodes to new cell id
 565
 566   // new nodes
 567   map< TEdge, int > cutEdge2newNodeId; // ids of nodes located between extremities of old edge
 568   map< int, int >   oldNode2newNodeId; // ids of nodes coincident with old node
 569
 570   // tolerance
 571   double tol = getTolerance( inMesh );
 572
 573
 574   // ----------------------------------------------
 575   // compute distances of nodes from plane or line
 576   // ----------------------------------------------
 577
 578   computeDistanceOfNodes(coords, normal);
 579
 580
 581   // ----------
 582   // Cut edges
 583   // ----------
 584
 585   int inCell = 0;
 586   for ( int iType = 0; iType < support->getNumberOfTypes(); ++iType) // loop on geom types
 587   {
 588     medGeometryElement type = inTypes[ iType ];
 589     TEdgeIterator edges( type );
 590
 591     const int* inCells = 0;
 592     if ( !support->isOnAllElements() )
 593       inCells = support->getNumber( type );
 594
 595     int nbInputCells = support->getNumberOfElements( type ); // loop on cells
 596     for ( int i = 0; i < nbInputCells; ++i )
 597     {
 598       int oldCell = inCells ? inCells[i] : ++inCell;
 599       const int* cellConn = inConn + (inConnIndex[ oldCell-1 ] - 1);
 600
 601       // Nodes of new mesh are either coincide with input nodes or lay on
 602       // edges cut by plane. If at least one edge of a cell is cut in the middle
 603       // then there will be the new cell, if the plane pass only through nodes of
 604       // input cell then the new cell is not necessarily created.
 605
 606       set<int> oldNodes, newNodes; //!< cut old nodes, new nodes on edges
 607       int nbEdgesOnPlane = 0;
 608       for ( int iEdge = 0; iEdge < edges.getNbEdges(); ++iEdge ) // loop on cell edges
 609       {
 610         // Analyse edge position in relation to the cutting plane or line
 611         const TEdge& edge = edges.getEdge( iEdge, cellConn );
 612         double dist1 = _myNodeDistance[ edge.node1()-1 ];
 613         double dist2 = _myNodeDistance[ edge.node2()-1 ];
 614         bool n1OnPlane = fabs( dist1 ) < tol;
 615         bool n2OnPlane = fabs( dist2 ) < tol;
 616         if ( n1OnPlane )
 617           oldNodes.insert( edge.node1() );
 618         if ( n2OnPlane )
 619           oldNodes.insert( edge.node2() );
 620         else if ( !n1OnPlane && dist1 * dist2 < 0 ) {
 621           // edge intersected
 622           int newNode = cutEdge2newNodeId.size() + oldNode2newNodeId.size() + 1;
 623           int node = cutEdge2newNodeId.insert( make_pair( edge, newNode )).first->second;
 624           newNodes.insert( node );
 625         }
 626         nbEdgesOnPlane += int( n1OnPlane && n2OnPlane );
 627       }
 628       int nbNodesInNewCell = oldNodes.size() + newNodes.size();
 629       if ( nbNodesInNewCell > maxNbNodesPerCell )
 630         throw MEDEXCEPTION(STRING(LOC) << "invalid input mesh");
 631
 632       if ( nbNodesInNewCell >= minNbNodesPerCell )
 633       {
 634         // Associate new and old cells
 635         int newCell = new2oldCells.size() + 1;
 636         // detect equal new cells on boundaries of old cells
 637         if ( newNodes.empty() && oldNodes.size() == nbEdgesOnPlane + int(meshDim==2)) {
 638           pair < map< set<int>, int>::iterator, bool > it_unique =
 639             oldNodes2newCell.insert( make_pair( oldNodes, newCell ));
 640           if ( !it_unique.second ) { // equal new faces
 641             int equalNewCell = it_unique.first->second;
 642             new2oldCells[ equalNewCell ].insert( oldCell );
 643             continue;
 644           }
 645         }
 646         set<int>& oldCells = // add a set of old cells to the end of new2oldCells
 647           new2oldCells.insert( new2oldCells.end(), make_pair(newCell, set<int>()))->second;
 648         oldCells.insert( oldCell );
 649
 650         // Store nodes
 651         vector< int >& connectivity =
 652           nbNodesInNewCell>4 ? newConnByNbNodes[_POLYGON] : newConnByNbNodes[nbNodesInNewCell];
 653         // nodes at edge intersection
 654         set<int>::iterator n, nEnd;
 655         for ( n = newNodes.begin(), nEnd = newNodes.end(); n != nEnd; ++n )
 656         {
 657           connectivity.push_back( *n );
 658         }
 659         // nodes coincident with input nodes
 660         for ( n = oldNodes.begin(), nEnd = oldNodes.end(); n != nEnd; ++n )
 661         {
 662           int newNode = cutEdge2newNodeId.size() + oldNode2newNodeId.size() + 1;
 663           int node = oldNode2newNodeId.insert( make_pair( *n, newNode )).first->second;
 664           connectivity.push_back( node );
 665         }
 666         if ( nbNodesInNewCell>4 )
 667           nbNodesPerPolygon.push_back( nbNodesInNewCell );
 668       }
 669     } // loop on cells, cutting thier edges
 670
 671   } // loop on geom types of input mesh
 672
 673
 674   // -----------------
 675   // Make coordinates
 676   // -----------------
 677
 678   int nbNodes = cutEdge2newNodeId.size() + oldNode2newNodeId.size();
 679   vector< double > resCoords( nbNodes * spaceDim );
 680
 681   const double* inCoords = inMesh->getCoordinates(MED_FULL_INTERLACE);
 682
 683   // nodes in the middle of edges
 684   map< TEdge, int >::iterator edge_node, enEnd = cutEdge2newNodeId.end();
 685   for ( edge_node = cutEdge2newNodeId.begin(); edge_node != enEnd; ++edge_node )
 686   {
 687     int newNode = edge_node->second;
 688     const TEdge& edge = edge_node->first;
 689
 690     double*     nodeCoords = & resCoords[ spaceDim * ( newNode-1 )];
 691     const double* n1Coords = inCoords +   spaceDim * ( edge.node1()-1 );
 692     const double* n2Coords = inCoords +   spaceDim * ( edge.node2()-1 );
 693
 694     double dist1 = _myNodeDistance[ edge.node1()-1 ];
 695     double dist2 = _myNodeDistance[ edge.node2()-1 ];
 696     double r1 = dist1 / ( dist1 - dist2 );
 697
 698     for ( int j = 0 ; j < spaceDim; ++j )
 699       nodeCoords[ j ] = ( 1.-r1 ) * n1Coords[ j ] + r1 * n2Coords[ j ];
 700   }
 701
 702   // nodes coincident with input nodes
 703   map< int, int >::iterator old_newNode, onEnd = oldNode2newNodeId.end();
 704   const size_t size = size_t( sizeof(double)*spaceDim );
 705   for ( old_newNode = oldNode2newNodeId.begin(); old_newNode != onEnd; ++old_newNode )
 706   {
 707     double*       newCoords = & resCoords[ spaceDim * ( old_newNode->second - 1 )];
 708     const double* oldCoords = inCoords + spaceDim * ( old_newNode->first - 1 );
 709     memcpy( newCoords, oldCoords, size );
 710   }
 711
 712   // --------------------
 713   // Sort nodes of cells
 714   // --------------------
 715
 716   sortNodes( newConnByNbNodes, &resCoords[0], coords, normal, nbNodesPerPolygon );
 717
 718   // ----------
 719   // Make mesh
 720   // ----------
 721
 722   // count classical types
 723   vector< medGeometryElement > types;
 724   vector< int > nbCellByType;
 725   map< int, vector< int > >::iterator nbNoConn, ncEnd =newConnByNbNodes.end();
 726   for ( nbNoConn = newConnByNbNodes.begin(); nbNoConn != ncEnd; ++nbNoConn )
 727   {
 728     int nbNodesPerCell = nbNoConn->first;
 729     int connSize = nbNoConn->second.size();
 730     if ( nbNodesPerCell >= 2 && nbNodesPerCell <= 4 && connSize > 0 ) {
 731       nbCellByType.push_back( connSize / nbNodesPerCell );
 732       types.push_back( medGeometryElement( (meshDim-1)*100 + nbNodesPerCell ));
 733     }
 734   }
 735   if ( types.empty() && newConnByNbNodes[_POLYGON].empty() )
 736     return 0;
 737
 738   MESHING* meshing = new MESHING();
 739
 740   meshing->setName(STRING("Cut of ") << inMesh->getName());
 741   meshing->setNumberOfTypes( types.size(), MED_CELL );
 742   meshing->setCoordinates( spaceDim, nbNodes, & resCoords[0],
 743                            inMesh->getCoordinatesSystem(), MED_FULL_INTERLACE );
 744   meshing->setTypes( &types[0], MED_CELL );
 745   meshing->setNumberOfElements( &nbCellByType[0], MED_CELL);
 746   for ( int i = 0; i < types.size(); ++i )
 747     meshing->setConnectivity( & newConnByNbNodes[ types[i]%100 ].front(), MED_CELL, types[i]);
 748   meshing->setMeshDimension( spaceDim /*meshDim-1*/ );
 749
 750   // polygons
 751   if ( !newConnByNbNodes[_POLYGON].empty() )
 752   {
 753     // make index
 754     vector<int> index;
 755     index.reserve( newConnByNbNodes[_POLYGON].size() / 5 );
 756     index.push_back( 1 );
 757     list<int>::iterator nbNodes = nbNodesPerPolygon.begin(), nnEnd = nbNodesPerPolygon.end();
 758     for ( ; nbNodes != nnEnd; ++nbNodes )
 759       index.push_back( index.back() + *nbNodes );
 760
 761     meshing->setPolygonsConnectivity( & index[0],
 762                                       & newConnByNbNodes[_POLYGON][0],
 763                                       index.size()-1,
 764                                       MED_CELL );
 765   }
 766   return meshing;
 767 }
 768
 769 //================================================================================
 770 /*!
 771  * \brief computes distance of each node from the plane or the line given by point and normal
 772  */
 773 //================================================================================
 774
 775 void Extractor::computeDistanceOfNodes(const double* point,
 776                                        const double* normal)
 777 {
 778   const MESH* mesh     = _myInputField->getSupport()->getMesh();
 779   const double * coord = mesh->getCoordinates(MED_FULL_INTERLACE);
 780   const int spaceDim   = mesh->getSpaceDimension();
 781
 782   _myNodeDistance.resize(mesh->getNumberOfNodes());
 783
 784   // compute dot product: normal * Vec(point,node)
 785   for ( int i = 0; i < mesh->getNumberOfNodes(); ++i )
 786   {
 787     _myNodeDistance[i] = 0.0;
 788     for ( int j = 0; j < spaceDim; ++j, ++coord )
 789     {
 790       _myNodeDistance[i] += normal[j] * (*coord - point[j]);
 791     }
 792   }
 793 }
 794
 795 //================================================================================
 796 /*!
 797  * \brief Orient elements correctly and sort nodes of polygons
 798  *  \param connByNbNodes - map of nb of nodes per cell to cell connectivity to rearrange
 799  *  \param nodeCoords - coordinates of nodes of a new mesh in full interlace
 800  *  \param point - point on plane or line
 801  *  \param normal - normal to plane or line
 802  *  \param nbNodesPerPolygon - index of connByNbNodes[_POLYGON] connectivity
 803  */
 804 //================================================================================
 805
 806 void Extractor::sortNodes( map< int, vector< int > >& connByNbNodes,
 807                            const double*              nodeCoords,
 808                            const double*              point,
 809                            const double*              normal,
 810                            const list<int> &          nbNodesPerPolygon)
 811 {
 812   const int spaceDim = _myInputField->getSupport()->getMesh()->getSpaceDimension();
 813
 814   if ( !connByNbNodes[2].empty() ) // 1D mesh
 815   {
 816     vector< int > & conn = connByNbNodes[2];
 817     if ( spaceDim == 2 )
 818     {
 819       // Orient edges along an coordinate axis,
 820       // select ordinate to check
 821       int ind = (fabs(normal[0]) < fabs(normal[1])) ? 1 : 0;
 822       // sorting
 823       for ( int i = 0; i < conn.size(); i += 2) {
 824         const double* p1 = nodeCoords + spaceDim*(conn[i]-1);
 825         const double* p2 = nodeCoords + spaceDim*(conn[i+1]-1);
 826         if ( p1[ind] > p2[ind] )
 827           std::swap( conn[i], conn[i+1] );
 828       }
 829     }
 830     else
 831     {
 832       // Reverse if necessary adjacent edges if they have equal nodes
 833       if ( conn.size() > 2 ) {
 834         if ( conn[0] == conn[2] || conn[0] == conn[3] )
 835           std::swap( conn[0], conn[1] );
 836         int i;
 837         for ( i = 2; i < conn.size()-2; i += 2) {
 838           if ( conn[i-1] == conn[i+1] )
 839             std::swap( conn[i], conn[i+1] );
 840           else if ( conn[i] == conn[i+2] || conn[i] == conn[i+3] )
 841             std::swap( conn[i], conn[i+1] );
 842         }
 843         if ( conn[i+1] == conn[i-1] )
 844           std::swap( conn[i], conn[i+1] );
 845       }
 846     }
 847     return;
 848   }
 849
 850   // Fix order of nodes
 851
 852   // select two ordinates for sorting
 853   int i1 = 0, i2 = 1, i3 = 2;
 854   if ( fabs( normal[i1] ) > fabs( normal[i3] )) swap( i1, i3 );
 855   if ( fabs( normal[i2] ) > fabs( normal[i3] )) swap( i2, i3 );
 856
 857   // comparator of nodes
 858   TNodeCompare nodeCompare( nodeCoords, i1, i2 );
 859
 860   map< int, vector< int > >::iterator nbN_conn = connByNbNodes.begin();
 861   for ( ; nbN_conn != connByNbNodes.end(); ++nbN_conn )
 862   {
 863     if ( nbN_conn->second.empty() ) continue;
 864
 865     int * conn    = & nbN_conn->second[0];
 866     int * connEnd = conn + nbN_conn->second.size();
 867
 868     int nbNodesPerFace = nbN_conn->first;
 869     list<int>::const_iterator nbPolyNodes, npEnd = nbNodesPerPolygon.end();
 870
 871     for ( nbPolyNodes = nbNodesPerPolygon.begin(); conn != connEnd; conn += nbNodesPerFace )
 872     {
 873       if ( nbPolyNodes != npEnd )
 874         nbNodesPerFace = *nbPolyNodes++;
 875
 876       // Sort nodes of polygons and quadrangles
 877
 878       if ( nbNodesPerFace > 3 )
 879       {
 880         // get barycenter
 881         double bary[2] = { 0., 0. };
 882         for ( int i = 0; i < nbNodesPerFace; ++i ) {
 883           const double* coord = nodeCoords + spaceDim*(conn[i]-1);
 884           bary[0] += coord[i1]; bary[1] += coord[i2];
 885         }
 886         bary[0] /= nbNodesPerFace; bary[1] /= nbNodesPerFace;
 887         nodeCompare.setBaryCenter( bary );
 888
 889         // sorting
 890         std::sort( conn, conn + nbNodesPerFace, nodeCompare);
 891       }
 892
 893       // Fix orientation of faces, orient them to have thier normal collinear with plane normal
 894
 895       // calculate cross product of two adjacent segments
 896       double dot = 0.;
 897       int i = 0;
 898       do {
 899         const double* p1 = nodeCoords + spaceDim*(conn[i+0]-1);
 900         const double* p2 = nodeCoords + spaceDim*(conn[i+1]-1);
 901         const double* p3 = nodeCoords + spaceDim*(conn[i+2]-1);
 902         double p2p1[2] = { p1[i1]-p2[i1] , p1[i2]-p2[i2] };
 903         double p2p3[2] = { p3[i1]-p2[i1] , p3[i2]-p2[i2] };
 904         dot = p2p3[1] * p2p1[0] - p2p3[0] * p2p1[1];
 905         ++i;
 906       }
 907       while ( dot == 0. && i+2 < nbNodesPerFace );
 908
 909       if ( dot * normal[i3] < 0. )
 910         std::reverse( conn, conn + nbNodesPerFace );
 911     }
 912   }
 913 }
 914 //================================================================================
 915 /*!
 916  * \brief Makes a 1D mesh by transfixing faces of 3D cells of input mesh by the line
 917  *  \param coords - give a point to pass through by the line
 918  *  \param direction - direction of the line
 919  *  \param new2oldCells - output map of new cells to cut input cell
 920  */
 921 //================================================================================
 922
 923 MESH* Extractor::transfixFaces( const double*       coords,
 924                                 const double*       direction,
 925                                 map<int,set<int> >& new2oldCells)
 926 {
 927   MESH* inMesh = _myInputField->getSupport()->getMesh();
 928   TMeshData inMeshData( *inMesh );
 929   TLine line( direction, coords );
 930
 931   const int nbFaces = inMesh->getNumberOfElements( MED_FACE, MED_ALL_ELEMENTS );
 932   const int dim = 3;
 933
 934   // Intersect 1st domain
 935
 936   vector< TIntersection*> chains;
 937   vector< pair< double, double > > ranges;
 938   int nbSegments = 0; // in the future mesh
 939
 940   set<int> checkedCells;
 941   checkedCells.insert(0); // 0 is returned by getCellsByFace() for boundary faces
 942
 943   int face = 1;
 944   for ( ; face <= nbFaces; ++face ) {
 945     if ( int soleCell = inMeshData.isFreeFace( face ))
 946       if ( checkedCells.insert(soleCell).second  &&
 947            canIntersect( soleCell, inMeshData, line ))
 948         if ( TIntersection* chain = intersect( soleCell, inMeshData, line, checkedCells )) {
 949           chains.push_back( chain );
 950           double min = DBL_MAX, max = -DBL_MAX;
 951           chain->getRange( min, max, line._maxDir );
 952           ranges.push_back( make_pair( min, max ));
 953           nbSegments += chain->size() - 1;
 954           break;
 955         }
 956   }
 957   if ( chains.empty() )
 958     return 0;
 959
 960   // Intersect the rest domains
 961
 962   for ( ; face <= nbFaces; ++face ) {
 963     if ( int soleCell = inMeshData.isFreeFace( face ))
 964       if ( checkedCells.insert(soleCell).second)
 965       {
 966         // check if at least one node of face is out of ranges of found chains
 967         TIter nodes = inMeshData.getFaceNodes( face );
 968         bool isOut = false;
 969         while ( nodes.more() && !isOut ) {
 970           double coord = inMeshData.getNodeCoord( nodes.next() )[ line._maxDir ];
 971           bool isIn = false;
 972           for ( int i = 0; i < ranges.size() && !isIn; ++i ) {
 973             const pair< double, double > & minMax = ranges[i];
 974             isIn = ( minMax.first < coord && coord < minMax.second );
 975           }
 976           isOut = !isIn;
 977         }
 978         // try to intersect
 979         if ( isOut && canIntersect( soleCell, inMeshData, line ))
 980           if ( TIntersection* chain = intersect( soleCell, inMeshData, line, checkedCells )) {
 981             chains.push_back( chain );
 982             double min = DBL_MAX, max = -DBL_MAX;
 983             chain->getRange( min, max, line._maxDir );
 984             ranges.push_back( make_pair( min, max ));
 985             nbSegments += chain->size() - 1;
 986           }
 987       }
 988   }
 989
 990   // Fill mesh data
 991
 992   int nbNodes = nbSegments + chains.size();
 993   vector< double > resCoords( nbNodes * dim );
 994   vector< int > resConn; resConn.reserve( nbSegments * 2 );
 995
 996   int iNode = 1, iSeg = 1;
 997   double* coord = & resCoords[0];
 998   const size_t cooSize = size_t( sizeof(double)*dim );
 999
1000   for ( int i = 0; i < chains.size(); ++i ) {
1001     TIntersection* section = chains[i];
1002     while ( section ) {
1003       memcpy( coord, section->_point, cooSize );
1004       coord += dim;
1005       if ( section->_prevSection ) {
1006         resConn.push_back( iNode++ );
1007         resConn.push_back( iNode );
1008         new2oldCells[ iSeg++ ] = section->_cells;
1009       }
1010       section = section->_prevSection;
1011     }
1012     iNode++;
1013     delete chains[i];
1014   }
1015
1016   // Create mesh
1017
1018   MESHING* meshing = new MESHING();
1019
1020   meshing->setName(STRING("Cut of ") << inMesh->getName());
1021   meshing->setNumberOfTypes( 1, MED_CELL );
1022   //meshing->setMeshDimension( dim );
1023   meshing->setCoordinates( dim, nbNodes, &resCoords[0],
1024                            inMesh->getCoordinatesSystem(), MED_FULL_INTERLACE );
1025   meshing->setTypes( &MED_SEG2, MED_CELL );
1026   meshing->setNumberOfElements( &nbSegments, MED_CELL);
1027   meshing->setConnectivity( & resConn[0], MED_CELL, MED_SEG2);
1028   meshing->setMeshDimension( dim );
1029
1030   return meshing;
1031 }
1032
1033 } // namespace MEDMEM
1034
1035
1036 //================================================================================
1037 /*!
1038  * \brief Constructs TEdgeIterator on given classical cell type
1039  */
1040 //================================================================================
1041
1042 TEdgeIterator::TEdgeIterator(const medGeometryElement type)
1043 {
1044   static map< medGeometryElement, vector<TEdge>* > _edgesByType;
1045   if ( _edgesByType.empty() ) {
1046     _edges = _edgesByType[ MED_TRIA3 ] = _edgesByType[ MED_TRIA6 ] = new vector<TEdge>();
1047     _edges->reserve( 3 );
1048     _edges->push_back( TEdge( 0, 1 ));
1049     _edges->push_back( TEdge( 1, 2 ));
1050     _edges->push_back( TEdge( 2, 0 ));
1051     _edges = _edgesByType[ MED_QUAD4 ] = _edgesByType[ MED_QUAD8 ] = new vector<TEdge>();
1052     _edges->reserve( 4 );
1053     _edges->push_back( TEdge( 0, 1 ));
1054     _edges->push_back( TEdge( 1, 2 ));
1055     _edges->push_back( TEdge( 2, 3 ));
1056     _edges->push_back( TEdge( 3, 0 ));
1057     _edges = _edgesByType[ MED_TETRA4 ] = _edgesByType[ MED_TETRA10 ] = new vector<TEdge>();
1058     _edges->reserve( 6 );
1059     _edges->push_back( TEdge( 0, 1 ));
1060     _edges->push_back( TEdge( 1, 2 ));
1061     _edges->push_back( TEdge( 2, 0 ));
1062     _edges->push_back( TEdge( 0, 3 ));
1063     _edges->push_back( TEdge( 1, 3 ));
1064     _edges->push_back( TEdge( 2, 3 ));
1065     _edges = _edgesByType[ MED_HEXA8 ] = _edgesByType[ MED_HEXA20 ] = new vector<TEdge>();
1066     _edges->reserve( 12 );
1067     _edges->push_back( TEdge( 0, 1 ));
1068     _edges->push_back( TEdge( 1, 2 ));
1069     _edges->push_back( TEdge( 2, 3 ));
1070     _edges->push_back( TEdge( 3, 0 ));
1071     _edges->push_back( TEdge( 4, 5 ));
1072     _edges->push_back( TEdge( 5, 6 ));
1073     _edges->push_back( TEdge( 6, 7 ));
1074     _edges->push_back( TEdge( 7, 4 ));
1075     _edges->push_back( TEdge( 0, 4 ));
1076     _edges->push_back( TEdge( 1, 5 ));
1077     _edges->push_back( TEdge( 2, 6 ));
1078     _edges->push_back( TEdge( 3, 7 ));
1079     _edges = _edgesByType[ MED_PYRA5 ] = _edgesByType[ MED_PYRA13 ] = new vector<TEdge>();
1080     _edges->reserve( 8 );
1081     _edges->push_back( TEdge( 0, 1 ));
1082     _edges->push_back( TEdge( 1, 2 ));
1083     _edges->push_back( TEdge( 2, 3 ));
1084     _edges->push_back( TEdge( 3, 0 ));
1085     _edges->push_back( TEdge( 0, 4 ));
1086     _edges->push_back( TEdge( 1, 4 ));
1087     _edges->push_back( TEdge( 2, 4 ));
1088     _edges->push_back( TEdge( 3, 4 ));
1089     _edges = _edgesByType[ MED_PENTA6 ] = _edgesByType[ MED_PENTA15 ] = new vector<TEdge>();
1090     _edges->reserve( 9 );
1091     _edges->push_back( TEdge( 0, 1 ));
1092     _edges->push_back( TEdge( 1, 2 ));
1093     _edges->push_back( TEdge( 2, 0 ));
1094     _edges->push_back( TEdge( 3, 4 ));
1095     _edges->push_back( TEdge( 4, 5 ));
1096     _edges->push_back( TEdge( 5, 3 ));
1097     _edges->push_back( TEdge( 0, 4 ));
1098     _edges->push_back( TEdge( 1, 5 ));
1099     _edges->push_back( TEdge( 2, 3 ));
1100     _edgesByType[ MED_NONE ]         = 0;
1101     _edgesByType[ MED_POINT1 ]       = 0;
1102     _edgesByType[ MED_SEG2 ]         = 0;
1103     _edgesByType[ MED_SEG3 ]         = 0;
1104     _edgesByType[ MED_POLYGON ]      = 0;
1105     _edgesByType[ MED_POLYHEDRA ]    = 0;
1106     _edgesByType[ MED_ALL_ELEMENTS ] = 0;
1107   }
1108   _edges = _edgesByType[type];
1109 }
1110
1111 namespace {
1112
1113 //================================================================================
1114 /*!
1115  * \brief Transfixes a cell with a line. Returns a head of chain of intersections
1116  *  \param cell - cell id to transfixe
1117  *  \param meshData - input 3D mesh
1118  *  \param line - cutting line
1119  *  \param checkedCells - set of cell already tried
1120  *  \param prevInter - previosly found intersection. Used to build a chain of
1121  *                     intersection via recursive call
1122  */
1123 //================================================================================
1124
1125 TIntersection* intersect(const int        cell,
1126                          const TMeshData& meshData,
1127                          const TLine&     line,
1128                          set<int>&        checkedCells,
1129                          TIntersection*   prevInter)
1130 {
1131   TIntersection* newSection = 0; // section to find
1132   TIntersection* auxSection = 0; // second section used when !prevInter
1133   list< TIntersection > bndSections;// list of intersection on edges
1134
1135   int avoidFace = prevInter ? prevInter->_face : -1;
1136
1137   TIter faces = meshData.getFaces( cell );
1138   while ( faces.more() ) /////////  loop on faces of the cell
1139   {
1140     int face = abs( faces.next() );
1141     if ( face == avoidFace )
1142       continue;
1143     TIter nodes = meshData.getFaceNodes( face );
1144
1145     // Get face normal
1146     // ----------------
1147
1148     double faceNormal[3];
1149     bool zeroArea = calcFaceNormal( face, meshData, faceNormal );
1150     if ( zeroArea )
1151       continue; // next face
1152
1153     // Is face parallel to line
1154     // -------------------------
1155
1156     double dirDotNorm = dotProduct( line._dir, faceNormal );
1157     const double* pFace0 = meshData.getNodeCoord( nodes[0] );
1158     if ( fabs( dirDotNorm ) < _TOLER )
1159     {
1160       // line || face, check if the line lays on the face
1161
1162       double lf[3] = { line._coord[0] - pFace0[0],
1163                        line._coord[1] - pFace0[1],
1164                        line._coord[2] - pFace0[2] };
1165       double lfDotNorm = dotProduct( lf, faceNormal );
1166       if ( fabs( lfDotNorm ) < _TOLER )
1167       {
1168         // =========================================================
1169         // Line lays on face. Intersect line with edges of the face
1170         // =========================================================
1171
1172         // Calculate distance of nodes from the line
1173         // ------------------------------------------
1174         double lineNormal[3];
1175         crossProduct( faceNormal, line._dir, lineNormal );
1176         vector<double> dist( nodes.size(), 0.0 );
1177         for ( int n = 0; n < nodes.size(); ++n )
1178         {
1179           const double* p = meshData.getNodeCoord( nodes[n] );
1180           for ( int j = 0; j < 3; ++j )
1181             dist[n] += lineNormal[j] * ( p[j] - line._coord[j] );
1182         }
1183         // Find intersections
1184         // -------------------
1185         vector<double> pCoords;  // intersection coordinates
1186         set<int>       pAtNodes; // intersected nodes of the face
1187         list<TEdge>    cutEdges; // intersected edges
1188         int nbPoints = 0;        // nb intersections
1189         pCoords.reserve(6);
1190         for ( int n = 0; n < nodes.size() && nbPoints < 2; ++n )
1191         {
1192           int n2 = (n+1) % nodes.size();
1193           double dist1 = dist[ n ];
1194           double dist2 = dist[ n2 ];
1195           bool n1OnLine = fabs( dist1 ) < meshData.tolerance();
1196           bool n2OnLine = fabs( dist2 ) < meshData.tolerance();
1197           if ( n1OnLine )
1198             pAtNodes.insert( nodes[n] );
1199           if ( n2OnLine )
1200             pAtNodes.insert( nodes[n2] );
1201           else if ( !n1OnLine && dist1 * dist2 < 0 ) {
1202             const double* p1 = meshData.getNodeCoord( nodes[n] );
1203             const double* p2 = meshData.getNodeCoord( nodes[n2] );
1204             double r1 = dist1 / ( dist1 - dist2 );
1205             for ( int j = 0 ; j < 3; ++j )
1206               pCoords.push_back(( 1.-r1 ) * p1[ j ] + r1 * p2[ j ]);
1207             cutEdges.push_back( TEdge( nodes[n], nodes[n2] ));
1208           }
1209           nbPoints = cutEdges.size() + pAtNodes.size();
1210         }
1211         // coords of intersection are stored in pCoords in order:
1212         // first go points on edges, then, points on nodes
1213         if ( nbPoints == 2 && !pAtNodes.empty() ) {
1214           set<int>::iterator n = pAtNodes.begin();
1215           while ( pCoords.size() != 6 ) { // there must be coords of two points
1216             const double* p = meshData.getNodeCoord( *n++ );
1217             for ( int j = 0 ; j < 3; ++j )
1218               pCoords.push_back( p[j] );
1219           }
1220         }
1221         // Store intersections
1222         // --------------------
1223         if ( nbPoints == 2 )
1224         {
1225           vector< TIntersection* > sections(nbPoints);
1226           const double* intCoord = & pCoords[0];
1227
1228           for ( int i = 0; i < nbPoints; ++i, intCoord += 3 )
1229           {
1230             // Set coords of intersection point
1231             sections[i] = new TIntersection;
1232             sections[i]->_point[0] = intCoord[0];
1233             sections[i]->_point[1] = intCoord[1];
1234             sections[i]->_point[2] = intCoord[2];
1235
1236             // Set intersected cells
1237             if ( cutEdges.empty() ) {
1238               // line can pass by edge shared by several cells
1239               TEdge cutEdge( *pAtNodes.begin(), *pAtNodes.rbegin() );
1240               getCellsSharingEdge( cutEdge, meshData, sections[i]->_cells );
1241             }
1242             if ( !cutEdges.empty() || sections[i]->_cells.empty() ) {
1243               // line pass by face between two cells
1244               TIter cells = meshData.getCellsByFace( face );
1245               while ( cells.more() )
1246                 if ( int elem = cells.next() )
1247                   sections[i]->_cells.insert( elem );
1248             }
1249             // Not to check the face at next search
1250             sections[i]->_face = face;
1251
1252             // Set nodes to start search of next intersection from
1253             if ( !cutEdges.empty() ) {
1254               sections[i]->_startNodes.insert( cutEdges.front().node1() );
1255               sections[i]->_startNodes.insert( cutEdges.front().node2() );
1256               cutEdges.pop_front();
1257             }
1258             else if ( pAtNodes.size() > 1 ) {
1259               set<int>::iterator p = pAtNodes.begin();
1260               if ( !line.isSame( intCoord, meshData.getNodeCoord( *p )))
1261                 ++p;
1262               sections[i]->_startNodes.insert( *p );
1263               pAtNodes.erase( p );
1264             }
1265             else {
1266               sections[i]->_startNodes.insert( *pAtNodes.begin() );
1267             }
1268           }
1269           if ( prevInter ) {
1270             // Only one point is needed, exclude already found intersection
1271             if ( line.isSame( prevInter->_point, sections[0]->_point ))
1272               std::swap( sections[0], sections[1] );
1273             delete sections[1];
1274             sections[1] = 0;
1275           }
1276           newSection = sections[0];
1277           auxSection = sections[1];
1278           if ( auxSection )
1279             auxSection->_cells = newSection->_cells;
1280
1281           bndSections.clear(); // remove already found intersections
1282
1283         } // if ( nbPoints == 2 )
1284
1285         break; // from loop on faces of cell
1286
1287       } // line lays on face
1288     }
1289     else
1290     {
1291       // ==================================
1292       // Line intersects the plane of face
1293       // ==================================
1294
1295       // Find distance of intersection point from line origin
1296       // t = faceNormal * ( pFace0 - line._coord ) / ( faceNormal * line._dir )
1297
1298       double pf0Coord[] = { pFace0[0] - line._coord[0],
1299                             pFace0[1] - line._coord[1],
1300                             pFace0[2] - line._coord[2] };
1301       double t = dotProduct( faceNormal, pf0Coord ) / dotProduct( faceNormal, line._dir );
1302
1303       // facePlane-line intersection point
1304       double ip[] = { line._coord[0] + line._dir[0] * t,
1305                       line._coord[1] + line._dir[1] * t,
1306                       line._coord[2] + line._dir[2] * t};
1307
1308       if ( prevInter && line.isSame( ip, prevInter->_point ))
1309         continue;
1310       if ( !bndSections.empty() && line.isSame( ip, bndSections.back()._point ))
1311         continue;
1312
1313       // Check if intersection point (ip) is inside the face.
1314       // ----------------------------------------------------
1315
1316       // do it in 2d, on the cartesian plane most normal to the face;
1317       // find indices on that plane: i1, i2
1318       int i1 = 0, i2 = 1, i3 = 2;
1319       if ( fabs( faceNormal[i1] ) > fabs( faceNormal[i3] )) swap( i1, i3 );
1320       if ( fabs( faceNormal[i2] ) > fabs( faceNormal[i3] )) swap( i2, i3 );
1321       if ( i2-i1 != 1 && i2 != 0 ) swap ( i1, i2 );
1322
1323       int inside = true, nbOnBoundary = 0;
1324       TEdge cutEdge;
1325       for ( int n = 0; n < nodes.size() && inside; ++n )
1326       {
1327         const double* p0 = meshData.getNodeCoord( nodes[n] );
1328         const double* p1 = meshData.getNodeCoord( nodes[ (n+1) % nodes.size() ] );
1329         double sign =
1330           faceNormal[i3]*((ip[i2] - p0[i2])*(p1[i1] - p0[i1]) - (ip[i1] - p0[i1])*(p1[i2] - p0[i2]));
1331         if ( sign < -DBL_MIN )
1332           inside = false;
1333         else if ( sign < DBL_MIN ) {
1334           nbOnBoundary++;
1335           cutEdge = TEdge( nodes, n );
1336         }
1337       }
1338
1339       // Store intersection point
1340       // -------------------------
1341       if ( inside )
1342       {
1343         TIntersection* section;
1344         if ( !nbOnBoundary )
1345           section = new TIntersection;
1346         else {
1347           if ( bndSections.size() >= 2 )
1348             continue; // go to next face
1349           bndSections.push_back( TIntersection() );
1350           section = & bndSections.back();
1351           // set nodes to start next searching from
1352           if ( nbOnBoundary == 1 ) {
1353             // edge is cut
1354             section->_startNodes.insert( cutEdge.node1() );
1355             section->_startNodes.insert( cutEdge.node2() );
1356           }
1357           else { // node is cut
1358             const double* p1 = meshData.getNodeCoord( cutEdge.node1() );
1359             if ( fabs( ip[i1]-p1[i1] ) < _TOLER && fabs( ip[i2]-p1[i2] ) < _TOLER  )
1360               section->_startNodes.insert( cutEdge.node1() );
1361             else
1362               section->_startNodes.insert( cutEdge.node2() );
1363           }
1364         }
1365         section->_point[0] = ip[0];
1366         section->_point[1] = ip[1];
1367         section->_point[2] = ip[2];
1368         section->_face     = face;
1369         section->_cells.insert( cell );
1370
1371         if ( !nbOnBoundary )
1372         {
1373           if ( !newSection )
1374             newSection = section;
1375           else
1376             auxSection = section;
1377           if ( prevInter || auxSection ) {
1378             bndSections.clear();
1379             break; // from loop on faces
1380           }
1381         }
1382       }
1383     }
1384   } // loop on faces of cell
1385
1386   // Care of intersections passing through edges
1387   // --------------------------------------------
1388
1389   if ( !bndSections.empty() )
1390   {
1391     if ( prevInter ) { // boundary section not equal to previous section
1392       if ( !newSection )
1393         newSection = new TIntersection( bndSections.front() );
1394     }
1395     else {
1396       if ( !newSection ) {
1397         newSection = new TIntersection( bndSections.front() );
1398         bndSections.pop_front();
1399       }
1400       if ( !auxSection && !bndSections.empty() ) {
1401         auxSection = new TIntersection( bndSections.front() );
1402       }
1403     }
1404   }
1405
1406   // Find the rest of chain starting from the found sections
1407   // --------------------------------------------------------
1408
1409   if ( newSection && ( prevInter || auxSection ))
1410   {
1411     TIntersection* chain[] = { newSection, auxSection };
1412     int chainLength[] = {0,0};
1413     for ( int i = 0; i < 2; ++i )
1414     {
1415       TIntersection* section = chain[i];
1416       if ( !section ) continue;
1417       // get cells to try to intersect next
1418       set<int> cellsToCheck;
1419       if ( !section->_startNodes.empty() ) {
1420         if ( section->_startNodes.size() == 1 ) {
1421           TIter cells = meshData.getCellsByNode( *section->_startNodes.begin() );
1422           cellsToCheck.insert( cells.begin(), cells.end() );
1423         }
1424         else {
1425           TEdge cutEdge( *section->_startNodes.begin(), *section->_startNodes.rbegin() );
1426           getCellsSharingEdge( cutEdge, meshData, cellsToCheck );
1427         }
1428       }
1429       else {
1430         TIter cells = meshData.getCellsByFace( section->_face );
1431         cellsToCheck.insert( cells.begin(), cells.end() );
1432       }
1433       // find the rest intersections
1434       chain[i] = 0;
1435       set<int>::iterator elem = cellsToCheck.begin(), elemEnd = cellsToCheck.end();
1436       for ( ; elem != elemEnd && !chain[i]; ++elem ) {
1437         if ( checkedCells.insert( *elem ).second &&
1438              section->_cells.find( *elem ) == section->_cells.end() )
1439         {
1440           chain[i] = intersect( *elem, meshData, line, checkedCells, section );
1441         }
1442       }
1443       if ( chain[i] )
1444         chainLength[i] = chain[i]->size();
1445     }
1446
1447     // Connect just found sections into a chain
1448     if ( prevInter ) {
1449       newSection->_prevSection = prevInter;
1450     }
1451     else {
1452       if ( chainLength[0] < chainLength[1] ) {
1453         std::swap( chain[0], chain[1] );
1454         std::swap( newSection, auxSection );
1455       }
1456       if ( chain[1] )
1457         chain[1]->reverse();
1458       newSection->_prevSection = auxSection;
1459     }
1460
1461     if ( chain[0] )
1462       return chain[0];
1463     return newSection;
1464   }
1465   else {
1466     delete newSection;
1467   }
1468
1469   return 0;
1470 }
1471
1472 //================================================================================
1473 /*!
1474  * \brief Evaluate if the line can intersect a cell
1475  */
1476 //================================================================================
1477
1478 bool canIntersect(const int        cell,
1479                   const TMeshData& meshData,
1480                   const TLine&     line)
1481 {
1482   // calculate bnd box of the cell
1483   double min[] = { DBL_MAX,DBL_MAX,DBL_MAX }, max[] = { -DBL_MAX,-DBL_MAX,-DBL_MAX };
1484
1485   TIter cellNodes = meshData.getCellNodes( cell );
1486   for ( int n = 0; n < cellNodes.size(); ++n ) {
1487     const double* node = meshData.getNodeCoord( cellNodes[n] );
1488     for ( int j = 0; j < 3; ++j ) {
1489       if ( node[j] < min[j] ) min[j] = node[j];
1490       if ( node[j] > max[j] ) max[j] = node[j];
1491     }
1492   }
1493   double xmin = 0, xmax = 0, ymin = 0, ymax = 0, zmin, zmax;
1494   double parmin, parmax, par1, par2;
1495   bool xToSet, yToSet;
1496   const double infinite = 1e100;
1497
1498   if (fabs(line._dir[0])>0.) {
1499     par1=(min[0]-line._coord[0])/line._dir[0];
1500     par2=(max[0]-line._coord[0])/line._dir[0];
1501     parmin=std::min(par1, par2);
1502     parmax=std::max(par1, par2);
1503     xToSet=true;
1504   }
1505   else {
1506     if (line._coord[0]<min[0] || max[0]<line._coord[0]) {
1507       return false;
1508     }
1509     xmin=line._coord[0];
1510     xmax=line._coord[0];
1511     parmin=-infinite;
1512     parmax=infinite;
1513     xToSet=false;
1514   }
1515
1516   if (fabs(line._dir[1])>0.) {
1517     par1=(min[1]-line._coord[1])/line._dir[1];
1518     par2=(max[1]-line._coord[1])/line._dir[1];
1519     if(parmax < std::min(par1,par2) || parmin > std::max(par1,par2))
1520       return false;
1521     parmin=std::max(parmin, std::min(par1,par2));
1522     parmax=std::min(parmax, std::max(par1,par2));
1523     yToSet=true;
1524   }
1525   else {
1526     if (line._coord[1]<min[1] || max[1]<line._coord[1]) {
1527       return false;
1528     }
1529     ymin=line._coord[1];
1530     ymax=line._coord[1];
1531     yToSet=false;
1532   }
1533
1534   if (fabs(line._dir[2])>0.) {
1535     par1=(min[2]-line._coord[2])/line._dir[2];
1536     par2=(max[2]-line._coord[2])/line._dir[2];
1537     if(parmax < std::min(par1,par2) || parmin > std::max(par1,par2))
1538       return false;
1539     parmin=std::max(parmin, std::min(par1,par2));
1540     parmax=std::min(parmax, std::max(par1,par2));
1541     par1=line._coord[2]+parmin*line._dir[2];
1542     par2=line._coord[2]+parmax*line._dir[2];
1543     zmin=std::min(par1, par2);
1544     zmax=std::max(par1, par2);
1545   }
1546   else {
1547     if (line._coord[2]<min[2] || max[2]<line._coord[2])
1548       return false;
1549     zmin=line._coord[2];
1550     zmax=line._coord[2];
1551   }
1552   if (zmax<min[2] || max[2]<zmin) return false;
1553
1554   if (xToSet) {
1555     par1=line._coord[0]+parmin*line._dir[0];
1556     par2=line._coord[0]+parmax*line._dir[0];
1557     xmin=std::min(par1, par2);
1558     xmax=std::max(par1, par2);
1559   }
1560   if (xmax<min[0] || max[0]<xmin) return false;
1561
1562   if (yToSet) {
1563     par1=line._coord[1]+parmin*line._dir[1];
1564     par2=line._coord[1]+parmax*line._dir[1];
1565     ymin=std::min(par1, par2);
1566     ymax=std::max(par1, par2);
1567   }
1568   if (ymax<min[1] || max[1]<ymin) return false;
1569
1570   return true;
1571 }
1572 } // unnamed namespace
1573