Back to index

salome-med  6.5.0
MEDMEM_Extractor.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 // File      : MEDMEM_Extractor.cxx
00021 // Created   : Thu Dec 18 11:10:11 2008
00022 // Author    : Edward AGAPOV (eap)
00023 //
00024 #include "MEDMEM_Extractor.hxx"
00025 
00026 #include <MEDMEM_Field.hxx>
00027 #include <MEDMEM_Mesh.hxx>
00028 #include <MEDMEM_Meshing.hxx>
00029 #include <MEDMEM_Support.hxx>
00030 
00031 #include <list>
00032 
00033 #include <math.h>
00034 
00035 using namespace MED_EN;
00036 using namespace std;
00037 using namespace MEDMEM;
00038 
00039 namespace { // local tools
00040 
00041   const int _POLYGON = -1; 
00042 
00043   const double _TOLER = 1e-12;
00044 
00045   //================================================================================
00049   void crossProduct(const double* v1, const double* v2, double* res)
00050   {
00051     res[0] = v1[1] * v2[2] - v1[2] * v2[1];
00052     res[1] = v1[2] * v2[0] - v1[0] * v2[2];
00053     res[2] = v1[0] * v2[1] - v1[1] * v2[0];
00054   }
00055   //================================================================================
00059   double dotProduct(const double* v1, const double* v2)
00060   {
00061     return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
00062   }
00063   //================================================================================
00067   class TIter
00068   {
00069     const int *_cur, *_end;
00070   public:
00071     TIter(const int* start, const int* end):_cur(start), _end(end) {}
00072     TIter(const int* start, int nb):_cur(start), _end(start+nb) {}
00073     int size() const { return _end - _cur; }
00074     int operator[] (int i) const { return _cur[i]; }
00075     bool more() const { return _cur < _end; }
00076     int next() { return *_cur++; }
00077     const int* begin() const { return _cur; }
00078     const int* end() const { return _end; }
00079   };
00080   //================================================================================
00084   struct TEdge: public pair<int,int>
00085   {
00086     TEdge(const int n1=0, const int n2=0): pair<int,int>(n1,n2)
00087     { if ( n2 < n1 ) first = n2, second = n1; }
00088     TEdge(const TIter& nodes, int i )
00089     { *this = TEdge( nodes[i], nodes[ (i+1)%nodes.size() ]); }
00090     int node1() const { return first; }
00091     int node2() const { return second; }
00092   };
00093   //================================================================================
00097   struct TEdgeIterator
00098   {
00099     vector<TEdge>* _edges;
00100     TEdgeIterator(const medGeometryElement type);
00101     int getNbEdges() const { return _edges->size(); }
00102     TEdge getEdge(int i, const int* cellConn ) const
00103     { return TEdge( cellConn[(*_edges)[i].node1()], cellConn[(*_edges)[i].node2()]); }
00104     TEdge getEdge(int i, const TIter& cellNodes ) const
00105     { return TEdge( cellNodes[(*_edges)[i].node1()], cellNodes[(*_edges)[i].node2()]); }
00106   };
00107   //================================================================================
00111   struct TNodeCompare
00112   {
00113     const double* _coords; 
00114     const double* _bc; 
00115     int _i1, _i2; 
00116 
00117     TNodeCompare(const double* nodeCoords, int i1, int i2):
00118       _coords(nodeCoords),_i1(i1),_i2(i2) {}
00119 
00120     void setBaryCenter(const double* bc) { _bc = bc; }
00121 
00122     bool operator()(const int pt1, const int pt2) {
00123       const double* p1 = _coords + 3*(pt1-1);
00124       const double* p2 = _coords + 3*(pt2-1);
00125       // calculate angles of Op1 and Op2 with the i1 axis on plane (i1,i2)
00126       double ang1 = atan2(p1[_i1] - _bc[0], p1[_i2] - _bc[1]);
00127       double ang2 = atan2(p2[_i1] - _bc[0], p2[_i2] - _bc[1]);
00128       return ang1 > ang2;
00129     }
00130   };
00131   //================================================================================
00135   double getTolerance(const MESH* mesh )
00136   {
00137     vector< vector<double> > bb = mesh->getBoundingBox();
00138     double diagonal = 0;
00139     for ( int j = 0; j < mesh->getSpaceDimension(); ++j ) {
00140       double dist = bb[0][j] - bb[1][j];
00141       diagonal += dist*dist;
00142     }
00143     return sqrt( diagonal ) * 1e-7;
00144   }
00145   //================================================================================
00149   struct TLine
00150   {
00151     const double* _dir;
00152     const double* _coord;
00153     int           _maxDir; 
00154     TLine( const double* dir, const double* coord): _dir(dir), _coord(coord) {
00155       _maxDir = 0;
00156       if ( fabs( _dir[_maxDir] ) < fabs( _dir[1] )) _maxDir = 1;
00157       if ( fabs( _dir[_maxDir] ) < fabs( _dir[2] )) _maxDir = 2;
00158     }
00160     bool isSame( const double* p1, const double* p2 ) const {
00161       return fabs(p1[_maxDir] - p2[_maxDir]) < _TOLER;
00162     }
00163   };
00164   //================================================================================
00174   struct TIntersection
00175   {
00176     double   _point[3];   
00177     set<int> _cells;      
00178     int      _face;       
00179     set<int> _startNodes; 
00180     TIntersection* _prevSection; 
00181 
00182     TIntersection(): _face(-1), _prevSection(NULL)
00183     {}
00184     ~TIntersection() {
00185       if ( _prevSection ) delete _prevSection; _prevSection=0;
00186     }
00187     void getRange( double& min, double& max, const int j ) const {
00188       if ( _point[j] < min ) min = _point[j];
00189       if ( _point[j] > max ) max = _point[j];
00190       if ( _prevSection ) _prevSection->getRange( min, max, j );
00191     }
00192     void reverse() {
00193       if ( _prevSection ) {
00194         _prevSection->reverse();
00195         _prevSection->_cells = _cells;
00196         _prevSection->_prevSection = this;
00197         _prevSection = NULL;
00198       }
00199     }
00200     int size() const { return 1 + ( _prevSection ? _prevSection->size() : 0 ); }
00201   };
00202   //================================================================================
00206   struct TMeshData
00207   {
00208     int           _dim;
00209     double        _tolerance;
00210     const double* _coord;
00211     const int*    _cellConn;
00212     const int*    _cellConnIndex;
00213     const int*    _faceConn;
00214     const int*    _faceConnIndex;
00215     const int*    _face2Cell;
00216     const int*    _face2CellIndex;
00217     const int*    _cell2Face;
00218     const int*    _cell2FaceIndex;
00219     const int*    _node2Cell;
00220     const int*    _node2CellIndex;
00221     TMeshData(const MESH &mesh)
00222     {
00223       _tolerance      = getTolerance(&mesh);
00224       _dim            = mesh.getSpaceDimension();
00225       _coord          = mesh.getCoordinates(MED_FULL_INTERLACE);
00226       _cellConn       = mesh.getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS);
00227       _cellConnIndex  = mesh.getConnectivityIndex(MED_NODAL, MED_CELL);
00228       _cell2Face      = mesh.getConnectivity( MED_DESCENDING, MED_CELL, MED_ALL_ELEMENTS);
00229       _cell2FaceIndex = mesh.getConnectivityIndex( MED_DESCENDING, MED_CELL );
00230       _face2Cell      = mesh.getReverseConnectivity( MED_DESCENDING );
00231       _face2CellIndex = mesh.getReverseConnectivityIndex( MED_DESCENDING );
00232       _faceConn       = mesh.getConnectivity( MED_NODAL, MED_FACE, MED_ALL_ELEMENTS);
00233       _faceConnIndex  = mesh.getConnectivityIndex(MED_NODAL, MED_FACE);
00234       _node2Cell      = mesh.getReverseConnectivity( MED_NODAL );
00235       _node2CellIndex = mesh.getReverseConnectivityIndex( MED_NODAL );
00236     }
00237     double tolerance() const {
00238       return _tolerance; }
00239     const double* getNodeCoord( int node ) const {
00240       return _coord + _dim*(node-1); }
00241     TIter getFaces( int cell ) const {
00242       return TIter( _cell2Face+_cell2FaceIndex[cell-1]-1, _cell2Face+_cell2FaceIndex[cell]-1 ); }
00243     TIter getCellNodes( int cell ) const {
00244       return TIter( _cellConn+_cellConnIndex[cell-1]-1, _cellConn+_cellConnIndex[cell]-1 ); }
00245     TIter getFaceNodes( int face ) const {
00246       face = abs(face);
00247       return TIter( _faceConn+_faceConnIndex[face-1]-1, _faceConn+_faceConnIndex[face]-1 ); }
00248     TIter getCellsByNode( int node ) const {
00249       return TIter( _node2Cell+_node2CellIndex[node-1]-1, _node2Cell+_node2CellIndex[node]-1 ); }
00250     TIter getCellsByFace( int face ) const {
00251       face = abs(face);
00252       return TIter( _face2Cell+_face2CellIndex[face-1]-1, _face2Cell+_face2CellIndex[face]-1 ); }
00253     int isFreeFace( int face ) const {
00254       TIter cells = getCellsByFace( face );
00255       return ( cells[1] == 0 ) ? cells[0] : 0; }
00256   };
00257   //================================================================================
00262   bool calcFaceNormal( const int        face,
00263                        const TMeshData& meshData,
00264                        double*          norm)
00265   {
00266     TIter nodes = meshData.getFaceNodes( face );
00267     bool zeroArea = false;
00268     int i = 0;
00269     do {
00270       const double* p1 = meshData.getNodeCoord( nodes[i] );
00271       const double* p2 = meshData.getNodeCoord( nodes[i+1] );
00272       const double* p3 = meshData.getNodeCoord( nodes[i+2] );
00273       double p2p1[3] = { p1[0]-p2[0], p1[1]-p2[1], p1[2]-p2[2] };
00274       double p2p3[3] = { p3[0]-p2[0], p3[1]-p2[1], p3[2]-p2[2] };
00275       crossProduct( p2p3, p2p1, norm );
00276       double normSize2 = dotProduct( norm, norm );
00277       zeroArea = (normSize2 < DBL_MIN);
00278     }
00279     while ( zeroArea && i < 2 );
00280     return zeroArea;
00281   }
00282   //================================================================================
00286   void getCellsSharingEdge( const TEdge& edge, const TMeshData& meshData, set<int> & theCells )
00287   {
00288     TIter cells = meshData.getCellsByNode( edge.node1() );
00289     while ( cells.more() ) {
00290       int cell = cells.next();
00291       TIter nodes = meshData.getCellNodes( cell );
00292       TEdgeIterator edgeIter( medGeometryElement( 300 + nodes.size() ));
00293       for ( int i = 0, nb = edgeIter.getNbEdges(); i < nb; ++i ) {
00294         TEdge e = edgeIter.getEdge( i , nodes );
00295         if ( edge == e ) {
00296           theCells.insert( cell );
00297           break;
00298         }
00299       }
00300     }
00301   }
00302   bool canIntersect(const int        cell,
00303                     const TMeshData& meshData,
00304                     const TLine&     line);
00305 
00306   TIntersection* intersect(const int        cell,
00307                            const TMeshData& meshData,
00308                            const TLine&     line,
00309                            set<int>&        checkedCells,
00310                            TIntersection*   prevInter=0);
00311 
00312 } // noname namespace
00313 
00314 namespace MEDMEM
00315 {
00316 
00317   //================================================================================
00328   //================================================================================
00329 
00330 Extractor::Extractor(const FIELD<double>& inputField) throw (MEDEXCEPTION)
00331 : _myInputField( & inputField )    
00332 {
00333   const char* LOC = "Extractor::Extractor(inputField) :";
00334 
00335   // Check if the input field complies with the conditions
00336 
00337   if ( !inputField.getSupport() )
00338     throw MEDEXCEPTION(STRING(LOC) << "InputField has NULL support");
00339 
00340   medEntityMesh entity = inputField.getSupport()->getEntity();
00341   if ( entity == MED_NODE || entity == MED_EDGE )
00342     throw MEDEXCEPTION(STRING(LOC) << "InputField has invalid supporting entity");
00343 
00344   if ( inputField.getSupport()->getNumberOfElements(MED_ALL_ELEMENTS) == 0 )
00345     throw MEDEXCEPTION(STRING(LOC) << "InputField has support of zero size");
00346 
00347   if ( inputField.getGaussPresence() && inputField.getNumberOfGaussPoints()[0] > 1 )
00348     throw MEDEXCEPTION(STRING(LOC) << "InputField is not constant be element");
00349 
00350   const GMESH* mesh = inputField.getSupport()->getMesh();
00351   if ( !mesh )
00352       throw MEDEXCEPTION(STRING(LOC) << "InputField has support with NULL mesh");
00353 
00354   if ( mesh->getSpaceDimension() < 2 )
00355       throw MEDEXCEPTION(STRING(LOC) << "InputField with 1D support not acceptable");
00356 
00357   if ( mesh->getNumberOfElements(MED_CELL, MED_POLYGON) > 0 ||
00358        mesh->getNumberOfElements(MED_CELL, MED_POLYHEDRA) > 0 )
00359       throw MEDEXCEPTION(STRING(LOC) << "InputField has supporting mesh with poly elements");
00360 
00361   if ( mesh->getMeshDimension() < 2 )
00362       throw MEDEXCEPTION(STRING(LOC) << "Invalid entity dimension of connectivity");
00363 
00364   _myInputField->addReference();
00365   _myInputMesh = mesh->convertInMESH();
00366 }
00367 
00368 Extractor::~Extractor()
00369 {
00370   _myInputField->removeReference();
00371   _myInputMesh->removeReference();
00372 }
00373 
00374 //================================================================================
00384 //================================================================================
00385 
00386 FIELD<double>* Extractor::extractPlane(const double* coords, const double* normal)
00387   throw (MEDEXCEPTION)
00388 {
00389   const char* LOC = "Extractor::extractPlane(const double* coords, const double* normal) :";
00390 
00391   // check agrument validity
00392   if ( !coords || !normal )
00393     throw MEDEXCEPTION(STRING(LOC) << "NULL argument");
00394 
00395   double normalSize = sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
00396   if ( normalSize <= DBL_MIN )
00397     throw MEDEXCEPTION(STRING(LOC) << "normal has zero size");
00398 
00399   if ( _myInputField->getSupport()->getMesh()->getSpaceDimension() < 3 )
00400     throw MEDEXCEPTION(STRING(LOC) << "Extraction by plane is possible in 3D space only ");
00401 
00402   double norm[3] = { normal[0] / normalSize, normal[1] / normalSize, normal[2] / normalSize };
00403 
00404   // cut mesh
00405   map<int,set<int> > new2oldCells;
00406   MESH* mesh = divideEdges( coords, norm, new2oldCells );
00407   if ( !mesh ) return 0;
00408 
00409   FIELD<double>* ret=makeField( new2oldCells, mesh );
00410   return ret;
00411 }
00412 
00413 //================================================================================
00423 //================================================================================
00424 
00425 FIELD<double>* Extractor::extractLine(const double* coords, const double* direction)
00426   throw (MEDEXCEPTION)
00427 {
00428   const char* LOC = "Extractor::extractLine(const double* coords, const double* direction) :";
00429 
00430   // check agrument validity
00431   if ( !coords || !direction )
00432     throw MEDEXCEPTION(STRING(LOC) << "NULL argument");
00433 
00434   double directionSize =
00435     sqrt(direction[0]*direction[0] + direction[1]*direction[1] + direction[2]*direction[2]);
00436   if ( directionSize <= DBL_MIN )
00437     throw MEDEXCEPTION(STRING(LOC) << "direction has zero size");
00438 
00439   const SUPPORT* support = _myInputField->getSupport();
00440   const medGeometryElement* inTypes = support->getTypes();
00441   const int meshDim = inTypes[ support->getNumberOfTypes()-1 ] / 100;
00442 
00443   if ( meshDim == 2 && support->getMesh()->getSpaceDimension() == 3 )
00444     throw MEDEXCEPTION(STRING(LOC) << "Extraction from 2D mesh not supported");
00445 
00446   map<int,set<int> > new2oldCells;
00447   MESH* mesh = 0;
00448   if ( meshDim == 2 )
00449   {
00450     double norm[2] = { direction[1] / directionSize,
00451                        direction[0] / directionSize,  };
00452     // cut mesh
00453     mesh = divideEdges( coords, norm, new2oldCells );
00454   }
00455   else
00456   {
00457     double dir[3] = { direction[0] / directionSize,
00458                       direction[1] / directionSize,
00459                       direction[2] / directionSize };
00460     mesh = transfixFaces( coords, dir, new2oldCells );
00461   }
00462 
00463   if ( !mesh ) return 0;
00464     
00465   FIELD<double>*ret=makeField( new2oldCells, mesh );
00466   return ret;
00467 }
00468 
00469 //================================================================================
00475 //================================================================================
00476 
00477 FIELD<double>* Extractor::makeField( const map<int,set<int> >& new2oldCells,
00478                                      MESH*                     mesh ) const
00479 {
00480   // make new field
00481   int nbComp = _myInputField->getNumberOfComponents();
00482   const SUPPORT *sup = mesh->getSupportOnAll( MED_CELL );
00483   FIELD<double> * outField = new FIELD<double>( sup, nbComp );
00484 
00485   outField->setComponentsNames       ( _myInputField->getComponentsNames() );
00486   outField->setName                  ( STRING("Extracted from ")<< _myInputField->getName() );
00487   outField->setDescription           ( STRING("Created by MEDMEM::Extractor"));
00488   outField->setComponentsDescriptions( _myInputField->getComponentsDescriptions() );
00489   outField->setMEDComponentsUnits    ( _myInputField->getMEDComponentsUnits() );
00490 
00491   sup->removeReference(); // to delete mesh as soon as outField dies
00492 
00493   // put values to the new field
00494 
00495   double* outValues = const_cast<double*>( outField->getValue() );
00496 
00497   map<int,set<int> >::const_iterator new_olds, noEnd = new2oldCells.end();
00498   for ( new_olds = new2oldCells.begin(); new_olds != noEnd; ++new_olds )
00499   {
00500     for ( int j = 0; j < nbComp; ++j )
00501     {
00502       int ind = ( new_olds->first - 1 ) * nbComp + j;
00503       outValues[ ind ] = 0.0;
00504       set<int>::const_iterator inCells = new_olds->second.begin();
00505       set<int>::const_iterator inEnd   = new_olds->second.end();    
00506       for ( ; inCells != inEnd; ++inCells )
00507         outValues[ ind ] += _myInputField->getValueIJ( *inCells, j+1 );
00508 
00509       outValues[ ind ] /= new_olds->second.size();
00510     }
00511   }
00512   return outField;
00513 }
00514 
00515 //================================================================================
00523 //================================================================================
00524 
00525 MESH* Extractor::divideEdges(const double*       coords,
00526                              const double*       normal,
00527                              map<int,set<int> >& new2oldCells)
00528 {
00529   const char* LOC = "Extractor::divideEdges()";
00530 
00531   const SUPPORT* support            = _myInputField->getSupport();
00532   medEntityMesh entity              = support->getEntity();
00533   const MESH* inMesh                = _myInputMesh;//support->getMesh();
00534   const medGeometryElement* inTypes = support->getTypes();
00535 
00536   const int* inConn      = inMesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
00537   const int* inConnIndex = inMesh->getConnectivityIndex(MED_NODAL, entity);
00538   const int spaceDim     = inMesh->getSpaceDimension();
00539   const int meshDim      = inTypes[ support->getNumberOfTypes()-1 ] / 100;
00540 
00541 
00542   // connectivity of new cells by nb of nodes per cell
00543   map< int, vector< int > > newConnByNbNodes;
00544   int nbInputCells = support->getNumberOfElements( MED_ALL_ELEMENTS );
00545   int minNbNodesPerCell = 2, maxNbNodesPerCell = 2;
00546   if ( meshDim == 3 ) {
00547     newConnByNbNodes[ 3 ].reserve( nbInputCells/2 );
00548     newConnByNbNodes[ 4 ].reserve( nbInputCells/2 );
00549     minNbNodesPerCell = 3;
00550     maxNbNodesPerCell = int (MED_POLYGON); // polygones allowed
00551   }
00552   else {
00553     newConnByNbNodes[ 2 ].reserve( nbInputCells/2 );
00554   }
00555   list<int> nbNodesPerPolygon; // to make connectivity index of polygons
00556 
00557   // new cells
00558   map< set<int>, int> oldNodes2newCell; 
00559 
00560   // new nodes
00561   map< TEdge, int > cutEdge2newNodeId; // ids of nodes located between extremities of old edge
00562   map< int, int >   oldNode2newNodeId; // ids of nodes coincident with old node
00563 
00564   // tolerance
00565   double tol = getTolerance( inMesh );
00566 
00567 
00568   // ----------------------------------------------
00569   // compute distances of nodes from plane or line
00570   // ----------------------------------------------
00571 
00572   computeDistanceOfNodes(coords, normal);
00573 
00574 
00575   // ----------
00576   // Cut edges
00577   // ----------
00578 
00579   int inCell = 0;
00580   for ( int iType = 0; iType < support->getNumberOfTypes(); ++iType) // loop on geom types
00581   {
00582     medGeometryElement type = inTypes[ iType ];
00583     TEdgeIterator edges( type );
00584 
00585     const int* inCells = 0;
00586     if ( !support->isOnAllElements() )
00587       inCells = support->getNumber( type );
00588 
00589     int nbInputCells = support->getNumberOfElements( type ); // loop on cells
00590     for ( int i = 0; i < nbInputCells; ++i )
00591     {
00592       int oldCell = inCells ? inCells[i] : ++inCell;
00593       const int* cellConn = inConn + (inConnIndex[ oldCell-1 ] - 1);
00594 
00595       // Nodes of new mesh are either coincide with input nodes or lay on
00596       // edges cut by plane. If at least one edge of a cell is cut in the middle
00597       // then there will be the new cell, if the plane pass only through nodes of
00598       // input cell then the new cell is not necessarily created.
00599 
00600       set<int> oldNodes, newNodes; 
00601       int nbEdgesOnPlane = 0;
00602       for ( int iEdge = 0; iEdge < edges.getNbEdges(); ++iEdge ) // loop on cell edges
00603       {
00604         // Analyse edge position in relation to the cutting plane or line
00605         const TEdge& edge = edges.getEdge( iEdge, cellConn );
00606         double dist1 = _myNodeDistance[ edge.node1()-1 ];
00607         double dist2 = _myNodeDistance[ edge.node2()-1 ];
00608         bool n1OnPlane = fabs( dist1 ) < tol;
00609         bool n2OnPlane = fabs( dist2 ) < tol;
00610         if ( n1OnPlane )
00611           oldNodes.insert( edge.node1() );
00612         if ( n2OnPlane )
00613           oldNodes.insert( edge.node2() );
00614         else if ( !n1OnPlane && dist1 * dist2 < 0 ) {
00615           // edge intersected
00616           int newNode = cutEdge2newNodeId.size() + oldNode2newNodeId.size() + 1;
00617           int node = cutEdge2newNodeId.insert( make_pair( edge, newNode )).first->second;
00618           newNodes.insert( node );
00619         }
00620         nbEdgesOnPlane += int( n1OnPlane && n2OnPlane );
00621       }
00622       int nbNodesInNewCell = oldNodes.size() + newNodes.size();
00623       if ( nbNodesInNewCell > maxNbNodesPerCell )
00624         throw MEDEXCEPTION(STRING(LOC) << "invalid input mesh");
00625 
00626       if ( nbNodesInNewCell >= minNbNodesPerCell )
00627       {
00628         // Associate new and old cells
00629         int newCell = new2oldCells.size() + 1;
00630         // detect equal new cells on boundaries of old cells
00631         if ( newNodes.empty() && (int)oldNodes.size() == nbEdgesOnPlane + int(meshDim==2)) {
00632           pair < map< set<int>, int>::iterator, bool > it_unique =
00633             oldNodes2newCell.insert( make_pair( oldNodes, newCell ));
00634           if ( !it_unique.second ) { // equal new faces
00635             int equalNewCell = it_unique.first->second;
00636             new2oldCells[ equalNewCell ].insert( oldCell );
00637             continue;
00638           }
00639         }
00640         set<int>& oldCells = // add a set of old cells to the end of new2oldCells
00641           new2oldCells.insert( new2oldCells.end(), make_pair(newCell, set<int>()))->second;
00642         oldCells.insert( oldCell );
00643 
00644         // Store nodes
00645         vector< int >& connectivity =
00646           nbNodesInNewCell>4 ? newConnByNbNodes[_POLYGON] : newConnByNbNodes[nbNodesInNewCell];
00647         // nodes at edge intersection
00648         set<int>::iterator n, nEnd;
00649         for ( n = newNodes.begin(), nEnd = newNodes.end(); n != nEnd; ++n )
00650         {
00651           connectivity.push_back( *n );
00652         }
00653         // nodes coincident with input nodes
00654         for ( n = oldNodes.begin(), nEnd = oldNodes.end(); n != nEnd; ++n )
00655         {
00656           int newNode = cutEdge2newNodeId.size() + oldNode2newNodeId.size() + 1;
00657           int node = oldNode2newNodeId.insert( make_pair( *n, newNode )).first->second;
00658           connectivity.push_back( node );
00659         }
00660         if ( nbNodesInNewCell>4 )
00661           nbNodesPerPolygon.push_back( nbNodesInNewCell );
00662       }
00663     } // loop on cells, cutting thier edges
00664 
00665   } // loop on geom types of input mesh
00666 
00667 
00668   // -----------------
00669   // Make coordinates
00670   // -----------------
00671 
00672   int nbNodes = cutEdge2newNodeId.size() + oldNode2newNodeId.size();
00673   vector< double > resCoords( nbNodes * spaceDim );
00674 
00675   const double* inCoords = inMesh->getCoordinates(MED_FULL_INTERLACE);
00676 
00677   // nodes in the middle of edges
00678   map< TEdge, int >::iterator edge_node, enEnd = cutEdge2newNodeId.end();
00679   for ( edge_node = cutEdge2newNodeId.begin(); edge_node != enEnd; ++edge_node )
00680   {
00681     int newNode = edge_node->second;
00682     const TEdge& edge = edge_node->first;
00683 
00684     double*     nodeCoords = & resCoords[ spaceDim * ( newNode-1 )];
00685     const double* n1Coords = inCoords +   spaceDim * ( edge.node1()-1 );
00686     const double* n2Coords = inCoords +   spaceDim * ( edge.node2()-1 );
00687 
00688     double dist1 = _myNodeDistance[ edge.node1()-1 ];
00689     double dist2 = _myNodeDistance[ edge.node2()-1 ];
00690     double r1 = dist1 / ( dist1 - dist2 );
00691 
00692     for ( int j = 0 ; j < spaceDim; ++j )
00693       nodeCoords[ j ] = ( 1.-r1 ) * n1Coords[ j ] + r1 * n2Coords[ j ];
00694   }
00695 
00696   // nodes coincident with input nodes
00697   map< int, int >::iterator old_newNode, onEnd = oldNode2newNodeId.end();
00698   const size_t size = size_t( sizeof(double)*spaceDim );
00699   for ( old_newNode = oldNode2newNodeId.begin(); old_newNode != onEnd; ++old_newNode )
00700   {
00701     double*       newCoords = & resCoords[ spaceDim * ( old_newNode->second - 1 )];
00702     const double* oldCoords = inCoords + spaceDim * ( old_newNode->first - 1 );
00703     memcpy( newCoords, oldCoords, size );
00704   }
00705 
00706   // --------------------
00707   // Sort nodes of cells
00708   // --------------------
00709   if ( nbNodes > 0 )
00710     sortNodes( newConnByNbNodes, &resCoords[0], coords, normal, nbNodesPerPolygon );
00711 
00712   // ----------
00713   // Make mesh
00714   // ----------
00715 
00716   // count types
00717   vector< medGeometryElement > types;
00718   vector< int > nbCellByType;
00719   map< int, vector< int > >::iterator nbNoConn, ncEnd =newConnByNbNodes.end();
00720   for ( nbNoConn = newConnByNbNodes.begin(); nbNoConn != ncEnd; ++nbNoConn )
00721   {
00722     int nbNodesPerCell = nbNoConn->first;
00723     int connSize = nbNoConn->second.size();
00724     if ( connSize == 0 ) continue;
00725     if ( nbNodesPerCell >= 2 && nbNodesPerCell <= 4 )
00726     {
00727       nbCellByType.push_back( connSize / nbNodesPerCell );
00728       types.push_back( medGeometryElement( (meshDim-1)*100 + nbNodesPerCell ));
00729     }
00730     else
00731     {
00732       nbCellByType.push_back( nbNodesPerPolygon.size() );
00733       types.push_back( MED_POLYGON );
00734     }
00735   }
00736   if ( types.empty() )
00737     return 0;
00738 
00739   MESHING* meshing = new MESHING();
00740 
00741   meshing->setName(STRING("Cut of ") << inMesh->getName());
00742   meshing->setNumberOfTypes( types.size(), MED_CELL );
00743   meshing->setCoordinates( spaceDim, nbNodes, & resCoords[0],
00744                            inMesh->getCoordinatesSystem(), MED_FULL_INTERLACE );
00745   meshing->setTypes( &types[0], MED_CELL );
00746   meshing->setNumberOfElements( &nbCellByType[0], MED_CELL);
00747   for ( unsigned i = 0; i < types.size(); ++i )
00748     if ( types[i] != MED_POLYGON )
00749     {
00750       meshing->setConnectivity( MED_CELL, types[i], & newConnByNbNodes[ types[i]%100 ].front());
00751     }
00752     else
00753     {
00754       // make index
00755       vector<int> index;
00756       index.reserve( nbNodesPerPolygon.size()+1 );
00757       index.push_back( 1 );
00758       list<int>::iterator nbNodes = nbNodesPerPolygon.begin(), nnEnd = nbNodesPerPolygon.end();
00759       for ( ; nbNodes != nnEnd; ++nbNodes )
00760         index.push_back( index.back() + *nbNodes );
00761     
00762       meshing->setConnectivity( MED_CELL, types[i], & newConnByNbNodes[ _POLYGON ].front(),
00763                                 & index[0]);
00764     }
00765 
00766   return meshing;
00767 }
00768 
00769 //================================================================================
00773 //================================================================================
00774 
00775 void Extractor::computeDistanceOfNodes(const double* point,
00776                                        const double* normal)
00777 {
00778   const MESH* mesh     = _myInputMesh;
00779   const double * coord = mesh->getCoordinates(MED_FULL_INTERLACE);
00780   const int spaceDim   = mesh->getSpaceDimension();
00781 
00782   _myNodeDistance.resize(mesh->getNumberOfNodes());
00783 
00784   // compute dot product: normal * Vec(point,node)
00785   for ( int i = 0; i < mesh->getNumberOfNodes(); ++i )
00786   {
00787     _myNodeDistance[i] = 0.0;
00788     for ( int j = 0; j < spaceDim; ++j, ++coord )
00789     {
00790       _myNodeDistance[i] += normal[j] * (*coord - point[j]);
00791     }
00792   }
00793 }
00794 
00795 //================================================================================
00804 //================================================================================
00805 
00806 void Extractor::sortNodes( map< int, vector< int > >& connByNbNodes,
00807                            const double*              nodeCoords,
00808                            const double*              point,
00809                            const double*              normal,
00810                            const list<int> &          nbNodesPerPolygon)
00811 {
00812   const int spaceDim = _myInputField->getSupport()->getMesh()->getSpaceDimension();
00813 
00814   if ( !connByNbNodes[2].empty() ) // 1D mesh
00815   {
00816     vector< int > & conn = connByNbNodes[2];
00817     if ( spaceDim == 2 )
00818     {
00819       // Orient edges along an coordinate axis,
00820       // select ordinate to check
00821       int ind = (fabs(normal[0]) < fabs(normal[1])) ? 1 : 0;
00822       // sorting
00823       for ( unsigned i = 0; i < conn.size(); i += 2) {
00824         const double* p1 = nodeCoords + spaceDim*(conn[i]-1);
00825         const double* p2 = nodeCoords + spaceDim*(conn[i+1]-1);
00826         if ( p1[ind] > p2[ind] )
00827           std::swap( conn[i], conn[i+1] );
00828       }
00829     }
00830     else
00831     {
00832       // Reverse if necessary adjacent edges if they have equal nodes
00833       if ( conn.size() > 2 ) {
00834         if ( conn[0] == conn[2] || conn[0] == conn[3] )
00835           std::swap( conn[0], conn[1] );
00836         int i;
00837         for ( i = 2; i < (int)conn.size()-2; i += 2) {
00838           if ( conn[i-1] == conn[i+1] )
00839             std::swap( conn[i], conn[i+1] );
00840           else if ( conn[i] == conn[i+2] || conn[i] == conn[i+3] )
00841             std::swap( conn[i], conn[i+1] );
00842         }
00843         if ( conn[i+1] == conn[i-1] )
00844           std::swap( conn[i], conn[i+1] );
00845       }
00846     }
00847     return;
00848   }
00849 
00850   // Fix order of nodes
00851 
00852   // select two ordinates for sorting
00853   int i1 = 0, i2 = 1, i3 = 2;
00854   if ( fabs( normal[i1] ) > fabs( normal[i3] )) swap( i1, i3 );
00855   if ( fabs( normal[i2] ) > fabs( normal[i3] )) swap( i2, i3 );
00856 
00857   // comparator of nodes
00858   TNodeCompare nodeCompare( nodeCoords, i1, i2 );
00859 
00860   map< int, vector< int > >::iterator nbN_conn = connByNbNodes.begin();
00861   for ( ; nbN_conn != connByNbNodes.end(); ++nbN_conn )
00862   {
00863     if ( nbN_conn->second.empty() ) continue;
00864 
00865     int * conn    = & nbN_conn->second[0];
00866     int * connEnd = conn + nbN_conn->second.size();
00867 
00868     int nbNodesPerFace = nbN_conn->first;
00869     list<int>::const_iterator nbPolyNodes, npEnd = nbNodesPerPolygon.end();
00870 
00871     for ( nbPolyNodes = nbNodesPerPolygon.begin(); conn != connEnd; conn += nbNodesPerFace )
00872     {
00873       if ( nbPolyNodes != npEnd )
00874         nbNodesPerFace = *nbPolyNodes++;
00875 
00876       // Sort nodes of polygons and quadrangles
00877       
00878       if ( nbNodesPerFace > 3 )
00879       {
00880         // get barycenter
00881         double bary[2] = { 0., 0. };
00882         for ( int i = 0; i < nbNodesPerFace; ++i ) {
00883           const double* coord = nodeCoords + spaceDim*(conn[i]-1);
00884           bary[0] += coord[i1]; bary[1] += coord[i2];
00885         }
00886         bary[0] /= nbNodesPerFace; bary[1] /= nbNodesPerFace;
00887         nodeCompare.setBaryCenter( bary );
00888 
00889         // sorting
00890         std::sort( conn, conn + nbNodesPerFace, nodeCompare);
00891       }
00892 
00893       // Fix orientation of faces, orient them to have thier normal collinear with plane normal
00894 
00895       // calculate cross product of two adjacent segments
00896       double dot = 0.;
00897       int i = 0;
00898       do {
00899         const double* p1 = nodeCoords + spaceDim*(conn[i+0]-1);
00900         const double* p2 = nodeCoords + spaceDim*(conn[i+1]-1);
00901         const double* p3 = nodeCoords + spaceDim*(conn[i+2]-1);
00902         double p2p1[2] = { p1[i1]-p2[i1] , p1[i2]-p2[i2] };
00903         double p2p3[2] = { p3[i1]-p2[i1] , p3[i2]-p2[i2] };
00904         dot = p2p3[1] * p2p1[0] - p2p3[0] * p2p1[1];
00905         ++i;
00906       }
00907       while ( dot == 0. && i+2 < nbNodesPerFace );
00908 
00909       if ( dot * normal[i3] < 0. )
00910         std::reverse( conn, conn + nbNodesPerFace );
00911     }
00912   }
00913 }
00914 //================================================================================
00921 //================================================================================
00922 
00923 MESH* Extractor::transfixFaces( const double*       coords,
00924                                 const double*       direction,
00925                                 map<int,set<int> >& new2oldCells)
00926 {
00927   const MESH* inMesh = _myInputMesh;
00928   TMeshData inMeshData( *inMesh );
00929   TLine line( direction, coords );
00930 
00931   const int nbFaces = inMesh->getNumberOfElements( MED_FACE, MED_ALL_ELEMENTS );
00932   const int dim = 3;
00933 
00934   // Intersect 1st domain
00935 
00936   vector< TIntersection*> chains;
00937   vector< pair< double, double > > ranges;
00938   int nbSegments = 0; // in the future mesh
00939 
00940   set<int> checkedCells;
00941   checkedCells.insert(0); // 0 is returned by getCellsByFace() for boundary faces
00942 
00943   int face = 1;
00944   for ( ; face <= nbFaces; ++face ) {
00945     if ( int soleCell = inMeshData.isFreeFace( face ))
00946       if ( checkedCells.insert(soleCell).second  &&
00947            canIntersect( soleCell, inMeshData, line ))
00948         if ( TIntersection* chain = intersect( soleCell, inMeshData, line, checkedCells )) {
00949           chains.push_back( chain );
00950           double min = DBL_MAX, max = -DBL_MAX;
00951           chain->getRange( min, max, line._maxDir );
00952           ranges.push_back( make_pair( min, max ));
00953           nbSegments += chain->size() - 1;
00954           break;
00955         }
00956   }
00957   if ( chains.empty() )
00958     return 0;
00959 
00960   // Intersect the rest domains
00961 
00962   for ( ; face <= nbFaces; ++face ) {
00963     if ( int soleCell = inMeshData.isFreeFace( face ))
00964       if ( checkedCells.insert(soleCell).second)
00965       {
00966         // check if at least one node of face is out of ranges of found chains
00967         TIter nodes = inMeshData.getFaceNodes( face );
00968         bool isOut = false;
00969         while ( nodes.more() && !isOut ) {
00970           double coord = inMeshData.getNodeCoord( nodes.next() )[ line._maxDir ];
00971           bool isIn = false;
00972           for ( unsigned i = 0; i < ranges.size() && !isIn; ++i ) {
00973             const pair< double, double > & minMax = ranges[i];
00974             isIn = ( minMax.first < coord && coord < minMax.second );
00975           }
00976           isOut = !isIn;
00977         }
00978         // try to intersect
00979         if ( isOut && canIntersect( soleCell, inMeshData, line ))
00980           if ( TIntersection* chain = intersect( soleCell, inMeshData, line, checkedCells )) {
00981             chains.push_back( chain );
00982             double min = DBL_MAX, max = -DBL_MAX;
00983             chain->getRange( min, max, line._maxDir );
00984             ranges.push_back( make_pair( min, max ));
00985             nbSegments += chain->size() - 1;
00986           }
00987       }
00988   }
00989 
00990   // Fill mesh data
00991 
00992   int nbNodes = nbSegments + chains.size();
00993   vector< double > resCoords( nbNodes * dim );
00994   vector< int > resConn; resConn.reserve( nbSegments * 2 );
00995 
00996   int iNode = 1, iSeg = 1;
00997   double* coord = & resCoords[0];
00998   const size_t cooSize = size_t( sizeof(double)*dim );
00999 
01000   for ( unsigned i = 0; i < chains.size(); ++i ) {
01001     TIntersection* section = chains[i];
01002     while ( section ) {
01003       memcpy( coord, section->_point, cooSize );
01004       coord += dim;
01005       if ( section->_prevSection ) {
01006         resConn.push_back( iNode++ );
01007         resConn.push_back( iNode );
01008         new2oldCells[ iSeg++ ] = section->_cells;
01009       }
01010       section = section->_prevSection;
01011     }
01012     iNode++;
01013     delete chains[i];
01014   }
01015 
01016   // Create mesh
01017 
01018   MESHING* meshing = new MESHING();
01019 
01020   meshing->setName(STRING("Cut of ") << inMesh->getName());
01021   meshing->setNumberOfTypes( 1, MED_CELL );
01022   meshing->setCoordinates( dim, nbNodes, &resCoords[0],
01023                            inMesh->getCoordinatesSystem(), MED_FULL_INTERLACE );
01024   meshing->setTypes( &MED_SEG2, MED_CELL );
01025   meshing->setNumberOfElements( &nbSegments, MED_CELL);
01026   meshing->setConnectivity( MED_CELL, MED_SEG2, & resConn[0]);
01027 
01028   return meshing;
01029 }
01030 
01031 } // namespace MEDMEM
01032 
01033 
01034 class MapGeoEdge : public map< medGeometryElement, vector<TEdge>* >
01035 {
01036 public:
01037   MapGeoEdge();
01038   ~MapGeoEdge();
01039 };
01040 
01041 MapGeoEdge::MapGeoEdge()
01042 {
01043   std::vector<TEdge> *edges=(*this)[MED_TRIA3]=(*this)[MED_TRIA6]=new vector<TEdge>();
01044   edges->reserve( 3 );
01045   edges->push_back( TEdge( 0, 1 ));
01046   edges->push_back( TEdge( 1, 2 ));
01047   edges->push_back( TEdge( 2, 0 ));
01048   edges=(*this)[MED_QUAD4]=(*this)[MED_QUAD8]=new vector<TEdge>();
01049   edges->reserve( 4 );
01050   edges->push_back( TEdge( 0, 1 ));
01051   edges->push_back( TEdge( 1, 2 ));
01052   edges->push_back( TEdge( 2, 3 ));
01053   edges->push_back( TEdge( 3, 0 ));
01054   edges=(*this)[MED_TETRA4]=(*this)[MED_TETRA10]=new vector<TEdge>();
01055   edges->reserve( 6 );
01056   edges->push_back( TEdge( 0, 1 ));
01057   edges->push_back( TEdge( 1, 2 ));
01058   edges->push_back( TEdge( 2, 0 ));
01059   edges->push_back( TEdge( 0, 3 ));
01060   edges->push_back( TEdge( 1, 3 ));
01061   edges->push_back( TEdge( 2, 3 ));
01062   edges=(*this)[MED_HEXA8]=(*this)[MED_HEXA20]=new vector<TEdge>();
01063   edges->reserve( 12 );
01064   edges->push_back( TEdge( 0, 1 ));
01065   edges->push_back( TEdge( 1, 2 ));
01066   edges->push_back( TEdge( 2, 3 ));
01067   edges->push_back( TEdge( 3, 0 ));
01068   edges->push_back( TEdge( 4, 5 ));
01069   edges->push_back( TEdge( 5, 6 ));
01070   edges->push_back( TEdge( 6, 7 ));
01071   edges->push_back( TEdge( 7, 4 ));
01072   edges->push_back( TEdge( 0, 4 ));
01073   edges->push_back( TEdge( 1, 5 ));
01074   edges->push_back( TEdge( 2, 6 ));
01075   edges->push_back( TEdge( 3, 7 ));
01076   edges=(*this)[MED_PYRA5]=(*this)[MED_PYRA13]=new vector<TEdge>();
01077   edges->reserve( 8 );
01078   edges->push_back( TEdge( 0, 1 ));
01079   edges->push_back( TEdge( 1, 2 ));
01080   edges->push_back( TEdge( 2, 3 ));
01081   edges->push_back( TEdge( 3, 0 ));
01082   edges->push_back( TEdge( 0, 4 ));
01083   edges->push_back( TEdge( 1, 4 ));
01084   edges->push_back( TEdge( 2, 4 ));
01085   edges->push_back( TEdge( 3, 4 ));
01086   edges=(*this)[MED_PENTA6]=(*this)[MED_PENTA15]=new vector<TEdge>();
01087   edges->reserve( 9 );
01088   edges->push_back( TEdge( 0, 1 ));
01089   edges->push_back( TEdge( 1, 2 ));
01090   edges->push_back( TEdge( 2, 0 ));
01091   edges->push_back( TEdge( 3, 4 ));
01092   edges->push_back( TEdge( 4, 5 ));
01093   edges->push_back( TEdge( 5, 3 ));
01094   edges->push_back( TEdge( 0, 4 ));
01095   edges->push_back( TEdge( 1, 5 ));
01096   edges->push_back( TEdge( 2, 3 ));
01097   (*this)[MED_NONE]         = 0;
01098   (*this)[MED_POINT1]       = 0;
01099   (*this)[MED_SEG2]         = 0;
01100   (*this)[MED_SEG3]         = 0;
01101   (*this)[MED_POLYGON]      = 0;
01102   (*this)[MED_POLYHEDRA]    = 0;
01103   (*this)[MED_ALL_ELEMENTS] = 0;
01104 }
01105 
01106 MapGeoEdge::~MapGeoEdge()
01107 {
01108   delete (*this)[MED_TRIA6];
01109   delete (*this)[MED_QUAD8];
01110   delete (*this)[MED_TETRA10];
01111   delete (*this)[MED_HEXA20];
01112   delete (*this)[MED_PYRA13];
01113   delete (*this)[MED_PENTA15];
01114 }
01115 
01116 //================================================================================
01120 //================================================================================
01121 
01122 TEdgeIterator::TEdgeIterator(const medGeometryElement type)
01123 {
01124   static MapGeoEdge _edgesByType;
01125   _edges = _edgesByType[type];
01126 }
01127 
01128 namespace {
01129 
01130 //================================================================================
01140 //================================================================================
01141 
01142 TIntersection* intersect(const int        cell,
01143                          const TMeshData& meshData,
01144                          const TLine&     line,
01145                          set<int>&        checkedCells,
01146                          TIntersection*   prevInter)
01147 {
01148   TIntersection* newSection = 0; // section to find
01149   TIntersection* auxSection = 0; // second section used when !prevInter
01150   list< TIntersection > bndSections;// list of intersection on edges
01151 
01152   int avoidFace = prevInter ? prevInter->_face : -1;
01153 
01154   TIter faces = meshData.getFaces( cell );
01155   while ( faces.more() ) 
01156   {
01157     int face = abs( faces.next() );
01158     if ( face == avoidFace )
01159       continue;
01160     TIter nodes = meshData.getFaceNodes( face );
01161 
01162     // Get face normal
01163     // ----------------
01164 
01165     double faceNormal[3];
01166     bool zeroArea = calcFaceNormal( face, meshData, faceNormal );
01167     if ( zeroArea )
01168       continue; // next face
01169 
01170     // Is face parallel to line
01171     // -------------------------
01172 
01173     double dirDotNorm = dotProduct( line._dir, faceNormal );
01174     const double* pFace0 = meshData.getNodeCoord( nodes[0] );
01175     if ( fabs( dirDotNorm ) < _TOLER )
01176     {
01177       // line || face, check if the line lays on the face
01178 
01179       double lf[3] = { line._coord[0] - pFace0[0],
01180                        line._coord[1] - pFace0[1],
01181                        line._coord[2] - pFace0[2] };
01182       double lfDotNorm = dotProduct( lf, faceNormal );
01183       if ( fabs( lfDotNorm ) < _TOLER )
01184       {
01185         // =========================================================
01186         // Line lays on face. Intersect line with edges of the face
01187         // =========================================================
01188 
01189         // Calculate distance of nodes from the line
01190         // ------------------------------------------
01191         double lineNormal[3];
01192         crossProduct( faceNormal, line._dir, lineNormal );
01193         vector<double> dist( nodes.size(), 0.0 );
01194         for ( int n = 0; n < nodes.size(); ++n )
01195         {
01196           const double* p = meshData.getNodeCoord( nodes[n] );
01197           for ( int j = 0; j < 3; ++j )
01198             dist[n] += lineNormal[j] * ( p[j] - line._coord[j] );
01199         }
01200         // Find intersections
01201         // -------------------
01202         vector<double> pCoords;  // intersection coordinates
01203         set<int>       pAtNodes; // intersected nodes of the face
01204         list<TEdge>    cutEdges; // intersected edges
01205         int nbPoints = 0;        // nb intersections
01206         pCoords.reserve(6);
01207         for ( int n = 0; n < nodes.size() && nbPoints < 2; ++n )
01208         {
01209           int n2 = (n+1) % nodes.size();
01210           double dist1 = dist[ n ];
01211           double dist2 = dist[ n2 ];
01212           bool n1OnLine = fabs( dist1 ) < meshData.tolerance();
01213           bool n2OnLine = fabs( dist2 ) < meshData.tolerance();
01214           if ( n1OnLine )
01215             pAtNodes.insert( nodes[n] );
01216           if ( n2OnLine )
01217             pAtNodes.insert( nodes[n2] );
01218           else if ( !n1OnLine && dist1 * dist2 < 0 ) {
01219             const double* p1 = meshData.getNodeCoord( nodes[n] );
01220             const double* p2 = meshData.getNodeCoord( nodes[n2] );
01221             double r1 = dist1 / ( dist1 - dist2 );
01222             for ( int j = 0 ; j < 3; ++j )
01223               pCoords.push_back(( 1.-r1 ) * p1[ j ] + r1 * p2[ j ]);
01224             cutEdges.push_back( TEdge( nodes[n], nodes[n2] ));
01225           }
01226           nbPoints = cutEdges.size() + pAtNodes.size();
01227         }
01228         // coords of intersection are stored in pCoords in order:
01229         // first go points on edges, then, points on nodes
01230         if ( nbPoints == 2 && !pAtNodes.empty() ) {
01231           set<int>::iterator n = pAtNodes.begin();
01232           while ( pCoords.size() != 6 ) { // there must be coords of two points
01233             const double* p = meshData.getNodeCoord( *n++ );
01234             for ( int j = 0 ; j < 3; ++j )
01235               pCoords.push_back( p[j] );
01236           }
01237         }
01238         // Store intersections
01239         // --------------------
01240         if ( nbPoints == 2 )
01241         {
01242           vector< TIntersection* > sections(nbPoints);
01243           const double* intCoord = & pCoords[0];
01244 
01245           for ( int i = 0; i < nbPoints; ++i, intCoord += 3 )
01246           {
01247             // Set coords of intersection point
01248             sections[i] = new TIntersection;
01249             sections[i]->_point[0] = intCoord[0];
01250             sections[i]->_point[1] = intCoord[1];
01251             sections[i]->_point[2] = intCoord[2];
01252 
01253             // Set intersected cells
01254             if ( cutEdges.empty() ) {
01255               // line can pass by edge shared by several cells
01256               TEdge cutEdge( *pAtNodes.begin(), *pAtNodes.rbegin() );
01257               getCellsSharingEdge( cutEdge, meshData, sections[i]->_cells );
01258             }
01259             if ( !cutEdges.empty() || sections[i]->_cells.empty() ) {
01260               // line pass by face between two cells
01261               TIter cells = meshData.getCellsByFace( face );
01262               while ( cells.more() )
01263                 if ( int elem = cells.next() )
01264                   sections[i]->_cells.insert( elem );
01265             }
01266             // Not to check the face at next search
01267             sections[i]->_face = face;
01268 
01269             // Set nodes to start search of next intersection from
01270             if ( !cutEdges.empty() ) {
01271               sections[i]->_startNodes.insert( cutEdges.front().node1() );
01272               sections[i]->_startNodes.insert( cutEdges.front().node2() );
01273               cutEdges.pop_front();
01274             }
01275             else if ( pAtNodes.size() > 1 ) {
01276               set<int>::iterator p = pAtNodes.begin();
01277               if ( !line.isSame( intCoord, meshData.getNodeCoord( *p )))
01278                 ++p;
01279               sections[i]->_startNodes.insert( *p );
01280               pAtNodes.erase( p );
01281             }
01282             else {
01283               sections[i]->_startNodes.insert( *pAtNodes.begin() );
01284             }
01285           }
01286           if ( prevInter ) {
01287             // Only one point is needed, exclude already found intersection
01288             if ( line.isSame( prevInter->_point, sections[0]->_point ))
01289               std::swap( sections[0], sections[1] );
01290             delete sections[1];
01291             sections[1] = 0;
01292           }
01293           newSection = sections[0];
01294           auxSection = sections[1];
01295           if ( auxSection )
01296             auxSection->_cells = newSection->_cells;
01297 
01298           bndSections.clear(); // remove already found intersections
01299 
01300         } // if ( nbPoints == 2 )
01301 
01302         break; // from loop on faces of cell
01303 
01304       } // line lays on face
01305     }
01306     else
01307     {
01308       // ==================================
01309       // Line intersects the plane of face
01310       // ==================================
01311 
01312       // Find distance of intersection point from line origin
01313       // t = faceNormal * ( pFace0 - line._coord ) / ( faceNormal * line._dir )
01314 
01315       double pf0Coord[] = { pFace0[0] - line._coord[0],
01316                             pFace0[1] - line._coord[1],
01317                             pFace0[2] - line._coord[2] };
01318       double t = dotProduct( faceNormal, pf0Coord ) / dotProduct( faceNormal, line._dir );
01319 
01320       // facePlane-line intersection point
01321       double ip[] = { line._coord[0] + line._dir[0] * t,
01322                       line._coord[1] + line._dir[1] * t,
01323                       line._coord[2] + line._dir[2] * t};
01324 
01325       if ( prevInter && line.isSame( ip, prevInter->_point ))
01326         continue;
01327       if ( !bndSections.empty() && line.isSame( ip, bndSections.back()._point ))
01328         continue;
01329 
01330       // Check if intersection point (ip) is inside the face.
01331       // ----------------------------------------------------
01332 
01333       // do it in 2d, on the cartesian plane most normal to the face;
01334       // find indices on that plane: i1, i2
01335       int i1 = 0, i2 = 1, i3 = 2;
01336       if ( fabs( faceNormal[i1] ) > fabs( faceNormal[i3] )) swap( i1, i3 );
01337       if ( fabs( faceNormal[i2] ) > fabs( faceNormal[i3] )) swap( i2, i3 );
01338       if ( i2-i1 != 1 && i2 != 0 ) swap ( i1, i2 );
01339 
01340       int inside = true, nbOnBoundary = 0;
01341       TEdge cutEdge;
01342       for ( int n = 0; n < nodes.size() && inside; ++n )
01343       {
01344         const double* p0 = meshData.getNodeCoord( nodes[n] );
01345         const double* p1 = meshData.getNodeCoord( nodes[ (n+1) % nodes.size() ] );
01346         double sign =
01347           faceNormal[i3]*((ip[i2] - p0[i2])*(p1[i1] - p0[i1]) - (ip[i1] - p0[i1])*(p1[i2] - p0[i2]));
01348         if ( sign < -DBL_MIN )
01349           inside = false;
01350         else if ( sign < DBL_MIN ) {
01351           nbOnBoundary++;
01352           cutEdge = TEdge( nodes, n );
01353         }
01354       }
01355 
01356       // Store intersection point
01357       // -------------------------
01358       if ( inside )
01359       {
01360         TIntersection* section;
01361         if ( !nbOnBoundary )
01362           section = new TIntersection;
01363         else {
01364           if ( bndSections.size() >= 2 )
01365             continue; // go to next face
01366           bndSections.push_back( TIntersection() );
01367           section = & bndSections.back();
01368           // set nodes to start next searching from
01369           if ( nbOnBoundary == 1 ) {
01370             // edge is cut
01371             section->_startNodes.insert( cutEdge.node1() );
01372             section->_startNodes.insert( cutEdge.node2() );
01373           }
01374           else { // node is cut
01375             const double* p1 = meshData.getNodeCoord( cutEdge.node1() );
01376             if ( fabs( ip[i1]-p1[i1] ) < _TOLER && fabs( ip[i2]-p1[i2] ) < _TOLER  )
01377               section->_startNodes.insert( cutEdge.node1() );
01378             else
01379               section->_startNodes.insert( cutEdge.node2() );
01380           }
01381         }
01382         section->_point[0] = ip[0];
01383         section->_point[1] = ip[1];
01384         section->_point[2] = ip[2];
01385         section->_face     = face;
01386         section->_cells.insert( cell );
01387 
01388         if ( !nbOnBoundary )
01389         {
01390           if ( !newSection )
01391             newSection = section;
01392           else
01393             auxSection = section;
01394           if ( prevInter || auxSection ) {
01395             bndSections.clear();
01396             break; // from loop on faces
01397           }
01398         }
01399       }
01400     }
01401   } // loop on faces of cell
01402 
01403   // Care of intersections passing through edges
01404   // --------------------------------------------
01405 
01406   if ( !bndSections.empty() )
01407   {
01408     if ( prevInter ) { // boundary section not equal to previous section
01409       if ( !newSection )
01410         newSection = new TIntersection( bndSections.front() );
01411     }
01412     else {
01413       if ( !newSection ) {
01414         newSection = new TIntersection( bndSections.front() );
01415         bndSections.pop_front();
01416       }
01417       if ( !auxSection && !bndSections.empty() ) {
01418         auxSection = new TIntersection( bndSections.front() );
01419       }
01420     }
01421   }
01422 
01423   // Find the rest of chain starting from the found sections
01424   // --------------------------------------------------------
01425 
01426   if ( newSection && ( prevInter || auxSection ))
01427   {
01428     TIntersection* chain[] = { newSection, auxSection };
01429     int chainLength[] = {0,0};
01430     for ( int i = 0; i < 2; ++i )
01431     {
01432       TIntersection* section = chain[i];
01433       if ( !section ) continue;
01434       // get cells to try to intersect next
01435       set<int> cellsToCheck;
01436       if ( !section->_startNodes.empty() ) {
01437         if ( section->_startNodes.size() == 1 ) {
01438           TIter cells = meshData.getCellsByNode( *section->_startNodes.begin() );
01439           cellsToCheck.insert( cells.begin(), cells.end() );
01440         }
01441         else {
01442           TEdge cutEdge( *section->_startNodes.begin(), *section->_startNodes.rbegin() );
01443           getCellsSharingEdge( cutEdge, meshData, cellsToCheck );
01444         }
01445       }
01446       else {
01447         TIter cells = meshData.getCellsByFace( section->_face );
01448         cellsToCheck.insert( cells.begin(), cells.end() );
01449       }
01450       // find the rest intersections
01451       chain[i] = 0;
01452       set<int>::iterator elem = cellsToCheck.begin(), elemEnd = cellsToCheck.end();
01453       for ( ; elem != elemEnd && !chain[i]; ++elem ) {
01454         if ( checkedCells.insert( *elem ).second &&
01455              section->_cells.find( *elem ) == section->_cells.end() )
01456         {
01457           chain[i] = intersect( *elem, meshData, line, checkedCells, section );
01458         }
01459       }
01460       if ( chain[i] )
01461         chainLength[i] = chain[i]->size();
01462     }
01463 
01464     // Connect just found sections into a chain
01465     if ( prevInter ) {
01466       newSection->_prevSection = prevInter;
01467     }
01468     else {
01469       if ( chainLength[0] < chainLength[1] ) {
01470         std::swap( chain[0], chain[1] );
01471         std::swap( newSection, auxSection );
01472       }
01473       if ( chain[1] )
01474         chain[1]->reverse();
01475       newSection->_prevSection = auxSection;
01476     }
01477 
01478     if ( chain[0] )
01479       return chain[0];
01480     return newSection;
01481   }
01482   else {
01483     delete newSection;
01484   }
01485 
01486   return 0;
01487 }
01488 
01489 //================================================================================
01493 //================================================================================
01494 
01495 bool canIntersect(const int        cell,
01496                   const TMeshData& meshData,
01497                   const TLine&     line)
01498 {
01499   // calculate bnd box of the cell
01500   double min[] = { DBL_MAX,DBL_MAX,DBL_MAX }, max[] = { -DBL_MAX,-DBL_MAX,-DBL_MAX };
01501 
01502   TIter cellNodes = meshData.getCellNodes( cell );
01503   for ( int n = 0; n < cellNodes.size(); ++n ) {
01504     const double* node = meshData.getNodeCoord( cellNodes[n] );
01505     for ( int j = 0; j < 3; ++j ) {
01506       if ( node[j] < min[j] ) min[j] = node[j];
01507       if ( node[j] > max[j] ) max[j] = node[j];
01508     }
01509   }
01510   double xmin = 0, xmax = 0, ymin = 0, ymax = 0, zmin, zmax;
01511   double parmin, parmax, par1, par2;
01512   bool xToSet, yToSet;
01513   const double infinite = 1e100;
01514 
01515   if (fabs(line._dir[0])>0.) {
01516     par1=(min[0]-line._coord[0])/line._dir[0];
01517     par2=(max[0]-line._coord[0])/line._dir[0];
01518     parmin=std::min(par1, par2);
01519     parmax=std::max(par1, par2);
01520     xToSet=true;
01521   }
01522   else {
01523     if (line._coord[0]<min[0] || max[0]<line._coord[0]) {
01524       return false;
01525     }
01526     xmin=line._coord[0];
01527     xmax=line._coord[0];
01528     parmin=-infinite;
01529     parmax=infinite;
01530     xToSet=false;
01531   }
01532 
01533   if (fabs(line._dir[1])>0.) {
01534     par1=(min[1]-line._coord[1])/line._dir[1];
01535     par2=(max[1]-line._coord[1])/line._dir[1];
01536     if(parmax < std::min(par1,par2) || parmin > std::max(par1,par2))
01537       return false;
01538     parmin=std::max(parmin, std::min(par1,par2));
01539     parmax=std::min(parmax, std::max(par1,par2));
01540     yToSet=true;
01541   }
01542   else {
01543     if (line._coord[1]<min[1] || max[1]<line._coord[1]) {
01544       return false;
01545     }
01546     ymin=line._coord[1];
01547     ymax=line._coord[1];
01548     yToSet=false;
01549   }
01550 
01551   if (fabs(line._dir[2])>0.) {
01552     par1=(min[2]-line._coord[2])/line._dir[2];
01553     par2=(max[2]-line._coord[2])/line._dir[2];
01554     if(parmax < std::min(par1,par2) || parmin > std::max(par1,par2))
01555       return false;
01556     parmin=std::max(parmin, std::min(par1,par2));
01557     parmax=std::min(parmax, std::max(par1,par2));
01558     par1=line._coord[2]+parmin*line._dir[2];
01559     par2=line._coord[2]+parmax*line._dir[2];
01560     zmin=std::min(par1, par2);
01561     zmax=std::max(par1, par2);
01562   }
01563   else {
01564     if (line._coord[2]<min[2] || max[2]<line._coord[2])
01565       return false;
01566     zmin=line._coord[2];
01567     zmax=line._coord[2];
01568   }
01569   if (zmax<min[2] || max[2]<zmin) return false;
01570 
01571   if (xToSet) {
01572     par1=line._coord[0]+parmin*line._dir[0];
01573     par2=line._coord[0]+parmax*line._dir[0];
01574     xmin=std::min(par1, par2);
01575     xmax=std::max(par1, par2);
01576   }
01577   if (xmax<min[0] || max[0]<xmin) return false;
01578 
01579   if (yToSet) {
01580     par1=line._coord[1]+parmin*line._dir[1];
01581     par2=line._coord[1]+parmax*line._dir[1];
01582     ymin=std::min(par1, par2);
01583     ymax=std::max(par1, par2);
01584   }
01585   if (ymax<min[1] || max[1]<ymin) return false;
01586 
01587   return true;
01588 }
01589 } // unnamed namespace
01590