src/MEDMEM/MEDMEM_Extractor.cxx

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