src/StdMeshers/StdMeshers_Import_1D2D.cxx

   1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 //  SMESH SMESH : implementaion of SMESH idl descriptions
  24 //  File   : StdMeshers_Import_1D2D.cxx
  25 //  Module : SMESH
  26 //
  27 #include "StdMeshers_Import_1D2D.hxx"
  28
  29 #include "StdMeshers_Import_1D.hxx"
  30 #include "StdMeshers_ImportSource.hxx"
  31
  32 #include "SMDS_MeshElement.hxx"
  33 #include "SMDS_MeshNode.hxx"
  34 #include "SMESHDS_Group.hxx"
  35 #include "SMESHDS_Mesh.hxx"
  36 #include "SMESH_Comment.hxx"
  37 #include "SMESH_Gen.hxx"
  38 #include "SMESH_Group.hxx"
  39 #include "SMESH_Mesh.hxx"
  40 #include "SMESH_MesherHelper.hxx"
  41 #include "SMESH_OctreeNode.hxx"
  42 #include "SMESH_subMesh.hxx"
  43
  44 #include "Utils_SALOME_Exception.hxx"
  45 #include "utilities.h"
  46
  47 #include <BRepClass_FaceClassifier.hxx>
  48 #include <BRepTools.hxx>
  49 #include <BRep_Builder.hxx>
  50 #include <BRep_Tool.hxx>
  51 #include <Bnd_B2d.hxx>
  52 #include <GeomAPI_ProjectPointOnSurf.hxx>
  53 #include <GeomAdaptor_Surface.hxx>
  54 #include <Precision.hxx>
  55 #include <TopExp.hxx>
  56 #include <TopExp_Explorer.hxx>
  57 #include <TopoDS.hxx>
  58 #include <TopoDS_Compound.hxx>
  59 #include <TopoDS_Edge.hxx>
  60 #include <TopoDS_Vertex.hxx>
  61
  62 #include <numeric>
  63
  64 using namespace std;
  65
  66 namespace
  67 {
  68   double getMinElemSize2( const SMESHDS_GroupBase* srcGroup )
  69   {
  70     double minSize2 = 1e100;
  71     SMDS_ElemIteratorPtr srcElems = srcGroup->GetElements();
  72     while ( srcElems->more() ) // loop on group contents
  73     {
  74       const SMDS_MeshElement* face = srcElems->next();
  75       int nbN = face->NbCornerNodes();
  76
  77       SMESH_TNodeXYZ prevN( face->GetNode( nbN-1 ));
  78       for ( int i = 0; i < nbN; ++i )
  79       {
  80         SMESH_TNodeXYZ n( face->GetNode( i ) );
  81         double size2 = ( n - prevN ).SquareModulus();
  82         minSize2 = std::min( minSize2, size2 );
  83         prevN = n;
  84       }
  85     }
  86     return minSize2;
  87   }
  88 }
  89
  90 //=============================================================================
  91 /*!
  92  * Creates StdMeshers_Import_1D2D
  93  */
  94 //=============================================================================
  95
  96 StdMeshers_Import_1D2D::StdMeshers_Import_1D2D(int hypId, int studyId, SMESH_Gen * gen)
  97   :SMESH_2D_Algo(hypId, studyId, gen), _sourceHyp(0)
  98 {
  99   MESSAGE("StdMeshers_Import_1D2D::StdMeshers_Import_1D2D");
 100   _name = "Import_1D2D";
 101   _shapeType = (1 << TopAbs_FACE);
 102
 103   _compatibleHypothesis.push_back("ImportSource2D");
 104   _requireDiscreteBoundary = false;
 105   _supportSubmeshes = true;
 106 }
 107
 108 //=============================================================================
 109 /*!
 110  * Check presence of a hypothesis
 111  */
 112 //=============================================================================
 113
 114 bool StdMeshers_Import_1D2D::CheckHypothesis
 115                          (SMESH_Mesh&                          aMesh,
 116                           const TopoDS_Shape&                  aShape,
 117                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 118 {
 119   _sourceHyp = 0;
 120
 121   const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
 122   if ( hyps.size() == 0 )
 123   {
 124     aStatus = SMESH_Hypothesis::HYP_MISSING;
 125     return false;  // can't work with no hypothesis
 126   }
 127
 128   if ( hyps.size() > 1 )
 129   {
 130     aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
 131     return false;
 132   }
 133
 134   const SMESHDS_Hypothesis *theHyp = hyps.front();
 135
 136   string hypName = theHyp->GetName();
 137
 138   if (hypName == _compatibleHypothesis.front())
 139   {
 140     _sourceHyp = (StdMeshers_ImportSource1D *)theHyp;
 141     aStatus = SMESH_Hypothesis::HYP_OK;
 142     return true;
 143   }
 144
 145   aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 146   return true;
 147 }
 148
 149 namespace
 150 {
 151   /*!
 152    * \brief OrientedLink additionally storing a medium node
 153    */
 154   struct TLink : public SMESH_OrientedLink
 155   {
 156     const SMDS_MeshNode* _medium;
 157     TLink( const SMDS_MeshNode* n1,
 158            const SMDS_MeshNode* n2,
 159            const SMDS_MeshNode* medium=0)
 160       : SMESH_OrientedLink( n1,n2 ), _medium( medium ) {}
 161   };
 162 }
 163
 164 //=============================================================================
 165 /*!
 166  * Import elements from the other mesh
 167  */
 168 //=============================================================================
 169
 170 bool StdMeshers_Import_1D2D::Compute(SMESH_Mesh & theMesh, const TopoDS_Shape & theShape)
 171 {
 172   if ( !_sourceHyp ) return false;
 173
 174   const vector<SMESH_Group*>& srcGroups = _sourceHyp->GetGroups(/*loaded=*/true);
 175   if ( srcGroups.empty() )
 176     return error("Invalid source groups");
 177
 178   bool allGroupsEmpty = true;
 179   for ( size_t iG = 0; iG < srcGroups.size() && allGroupsEmpty; ++iG )
 180     allGroupsEmpty = srcGroups[iG]->GetGroupDS()->IsEmpty();
 181   if ( allGroupsEmpty )
 182     return error("No faces in source groups");
 183
 184   SMESH_MesherHelper helper(theMesh);
 185   helper.SetSubShape(theShape);
 186   SMESHDS_Mesh* tgtMesh = theMesh.GetMeshDS();
 187
 188   const TopoDS_Face& geomFace = TopoDS::Face( theShape );
 189   const double faceTol = helper.MaxTolerance( geomFace );
 190   const int shapeID = tgtMesh->ShapeToIndex( geomFace );
 191   const bool toCheckOri = (helper.NbAncestors( geomFace, theMesh, TopAbs_SOLID ) == 1 );
 192
 193   Handle(Geom_Surface) surface = BRep_Tool::Surface( geomFace );
 194   const bool reverse =
 195     ( helper.GetSubShapeOri( tgtMesh->ShapeToMesh(), geomFace ) == TopAbs_REVERSED );
 196   gp_Pnt p; gp_Vec du, dv;
 197
 198   BRepClass_FaceClassifier classifier;
 199   Bnd_B2d bndBox2d;
 200   {
 201     Standard_Real umin,umax,vmin,vmax;
 202     BRepTools::UVBounds(geomFace,umin,umax,vmin,vmax);
 203     gp_XY pmin( umin,vmin ), pmax( umax,vmax );
 204     bndBox2d.Add( pmin );
 205     bndBox2d.Add( pmax );
 206     if ( helper.HasSeam() )
 207     {
 208       const int i = helper.GetPeriodicIndex();
 209       pmin.SetCoord( i, helper.GetOtherParam( pmin.Coord( i )));
 210       pmax.SetCoord( i, helper.GetOtherParam( pmax.Coord( i )));
 211       bndBox2d.Add( pmin );
 212       bndBox2d.Add( pmax );
 213     }
 214     bndBox2d.Enlarge( 1e-2 * Sqrt( bndBox2d.SquareExtent() ));
 215   }
 216
 217   set<int> subShapeIDs;
 218   subShapeIDs.insert( shapeID );
 219
 220   // nodes already existing on sub-shapes of the FACE
 221   TIDSortedNodeSet existingNodes;
 222
 223   // get/make nodes on vertices and add them to existingNodes
 224   TopExp_Explorer exp( theShape, TopAbs_VERTEX );
 225   for ( ; exp.More(); exp.Next() )
 226   {
 227     const TopoDS_Vertex& v = TopoDS::Vertex( exp.Current() );
 228     if ( !subShapeIDs.insert( tgtMesh->ShapeToIndex( v )).second )
 229       continue;
 230     const SMDS_MeshNode* n = SMESH_Algo::VertexNode( v, tgtMesh );
 231     if ( !n )
 232     {
 233       _gen->Compute(theMesh,v,/*anUpward=*/true);
 234       n = SMESH_Algo::VertexNode( v, tgtMesh );
 235       if ( !n ) return false; // very strange
 236     }
 237     existingNodes.insert( n );
 238   }
 239
 240   // get EDGEs and their ids and get existing nodes on EDGEs
 241   vector< TopoDS_Edge > edges;
 242   for ( exp.Init( theShape, TopAbs_EDGE ); exp.More(); exp.Next() )
 243   {
 244     const TopoDS_Edge & edge = TopoDS::Edge( exp.Current() );
 245     if ( !SMESH_Algo::isDegenerated( edge ))
 246       if ( subShapeIDs.insert( tgtMesh->ShapeToIndex( edge )).second )
 247       {
 248         edges.push_back( edge );
 249         if ( SMESHDS_SubMesh* eSM = tgtMesh->MeshElements( edge ))
 250         {
 251           typedef SMDS_StdIterator< const SMDS_MeshNode*, SMDS_NodeIteratorPtr > iterator;
 252           existingNodes.insert( iterator( eSM->GetNodes() ), iterator() );
 253         }
 254       }
 255   }
 256   // octree to find existing nodes
 257   SMESH_OctreeNode existingNodeOcTr( existingNodes );
 258   std::map<double, const SMDS_MeshNode*> dist2foundNodes;
 259
 260   // to count now many times a link between nodes encounters
 261   map<TLink, int> linkCount;
 262   map<TLink, int>::iterator link2Nb;
 263   double minGroupTol = Precision::Infinite();
 264
 265   // =========================
 266   // Import faces from groups
 267   // =========================
 268
 269   StdMeshers_Import_1D::TNodeNodeMap* n2n;
 270   StdMeshers_Import_1D::TElemElemMap* e2e;
 271   vector<TopAbs_State>         nodeState;
 272   vector<const SMDS_MeshNode*> newNodes; // of a face
 273   set   <const SMDS_MeshNode*> bndNodes; // nodes classified ON
 274
 275   for ( size_t iG = 0; iG < srcGroups.size(); ++iG )
 276   {
 277     const SMESHDS_GroupBase* srcGroup = srcGroups[iG]->GetGroupDS();
 278
 279     const int meshID = srcGroup->GetMesh()->GetPersistentId();
 280     const SMESH_Mesh* srcMesh = GetMeshByPersistentID( meshID );
 281     if ( !srcMesh ) continue;
 282     StdMeshers_Import_1D::getMaps( srcMesh, &theMesh, n2n, e2e );
 283
 284     const double groupTol = 0.5 * sqrt( getMinElemSize2( srcGroup ));
 285     minGroupTol = std::min( groupTol, minGroupTol );
 286
 287     GeomAdaptor_Surface S( surface );
 288     const double clsfTol = Min( S.UResolution( 0.1 * groupTol ),
 289                                 S.VResolution( 0.1 * groupTol ));
 290
 291     StdMeshers_Import_1D::TNodeNodeMap::iterator n2nIt;
 292     pair< StdMeshers_Import_1D::TNodeNodeMap::iterator, bool > it_isnew;
 293
 294     SMDS_ElemIteratorPtr srcElems = srcGroup->GetElements();
 295     while ( srcElems->more() ) // loop on group contents
 296     {
 297       const SMDS_MeshElement* face = srcElems->next();
 298
 299       // find or create nodes of a new face
 300       nodeState.resize( face->NbNodes() );
 301       newNodes.resize( nodeState.size() );
 302       newNodes.back() = 0;
 303       int nbCreatedNodes = 0;
 304       bool isOut = false, isIn = false; // if at least one node isIn - do not classify other nodes
 305       SMDS_MeshElement::iterator node = face->begin_nodes();
 306       for ( size_t i = 0; i < newNodes.size(); ++i, ++node )
 307       {
 308         SMESH_TNodeXYZ nXYZ = *node;
 309         nodeState[ i ] = TopAbs_UNKNOWN;
 310
 311         it_isnew = n2n->insert( make_pair( *node, (SMDS_MeshNode*)0 ));
 312         n2nIt    = it_isnew.first;
 313
 314         const SMDS_MeshNode* & newNode = n2nIt->second;
 315         if ( !it_isnew.second && !newNode )
 316           break; // a node is mapped to NULL - it is OUT of the FACE
 317
 318         if ( newNode )
 319         {
 320           if ( !subShapeIDs.count( newNode->getshapeId() ))
 321             break; // node is Imported onto other FACE
 322           if ( !isIn && bndNodes.count( *node ))
 323             nodeState[ i ] = TopAbs_ON;
 324         }
 325         else
 326         {
 327           // find a node pre-existing on EDGE or VERTEX
 328           dist2foundNodes.clear();
 329           existingNodeOcTr.NodesAround( nXYZ, dist2foundNodes, groupTol );
 330           if ( !dist2foundNodes.empty() )
 331           {
 332             newNode = dist2foundNodes.begin()->second;
 333             nodeState[ i ] = TopAbs_ON;
 334           }
 335         }
 336
 337         if ( !newNode )
 338         {
 339           // find out if node lies on the surface of theShape
 340           gp_XY uv( Precision::Infinite(), 0 );
 341           isOut = ( !helper.CheckNodeUV( geomFace, *node, uv, groupTol, /*force=*/true ) ||
 342                     bndBox2d.IsOut( uv ));
 343           if ( !isOut && !isIn ) // classify
 344           {
 345             classifier.Perform( geomFace, uv, clsfTol );
 346             nodeState[i] = classifier.State();
 347             isOut = ( nodeState[i] == TopAbs_OUT );
 348           }
 349           if ( !isOut ) // create a new node
 350           {
 351             newNode = tgtMesh->AddNode( nXYZ.X(), nXYZ.Y(), nXYZ.Z());
 352             tgtMesh->SetNodeOnFace( newNode, shapeID, uv.X(), uv.Y() );
 353             nbCreatedNodes++;
 354             if ( nodeState[i] == TopAbs_ON )
 355               bndNodes.insert( *node );
 356             else
 357               isIn = true;
 358           }
 359         }
 360         if ( !(newNodes[i] = newNode ) || isOut )
 361           break;
 362       }
 363
 364       if ( !newNodes.back() )
 365         continue; // not all nodes of the face lie on theShape
 366
 367       if ( !isIn ) // if all nodes are on FACE boundary, a mesh face can be OUT
 368       {
 369         // check state of nodes created for other faces
 370         for ( size_t i = 0; i < nodeState.size() && !isIn; ++i )
 371         {
 372           if ( nodeState[i] != TopAbs_UNKNOWN ) continue;
 373           gp_XY uv = helper.GetNodeUV( geomFace, newNodes[i] );
 374           classifier.Perform( geomFace, uv, clsfTol );
 375           nodeState[i] = classifier.State();
 376           isIn = ( nodeState[i] == TopAbs_IN );
 377         }
 378         if ( !isIn ) // classify face center
 379         {
 380           gp_XYZ gc( 0., 0., 0 );
 381           for ( size_t i = 0; i < newNodes.size(); ++i )
 382             gc += SMESH_TNodeXYZ( newNodes[i] );
 383           gc /= newNodes.size();
 384
 385           TopLoc_Location loc;
 386           GeomAPI_ProjectPointOnSurf& proj = helper.GetProjector( geomFace,
 387                                                                   loc,
 388                                                                   helper.MaxTolerance( geomFace ));
 389           if ( !loc.IsIdentity() ) loc.Transformation().Inverted().Transforms( gc );
 390           proj.Perform( gc );
 391           if ( !proj.IsDone() || proj.NbPoints() < 1 )
 392             continue;
 393           Quantity_Parameter U,V;
 394           proj.LowerDistanceParameters(U,V);
 395           gp_XY uv( U,V );
 396           classifier.Perform( geomFace, uv, clsfTol );
 397           if ( classifier.State() != TopAbs_IN )
 398             continue;
 399         }
 400       }
 401
 402       // try to find already created face
 403       SMDS_MeshElement * newFace = 0;
 404       if ( nbCreatedNodes == 0 &&
 405            tgtMesh->FindElement(newNodes, SMDSAbs_Face, /*noMedium=*/false))
 406         continue; // repeated face in source groups already created
 407
 408       // check future face orientation
 409       const int nbCorners = face->NbCornerNodes();
 410       const bool isQuad   = ( nbCorners != (int) newNodes.size() );
 411       if ( toCheckOri )
 412       {
 413         int iNode = -1;
 414         gp_Vec geomNorm;
 415         do
 416         {
 417           gp_XY uv = helper.GetNodeUV( geomFace, newNodes[++iNode] );
 418           surface->D1( uv.X(),uv.Y(), p, du,dv );
 419           geomNorm = reverse ? dv^du : du^dv;
 420         }
 421         while ( geomNorm.SquareMagnitude() < 1e-6 && iNode+1 < nbCorners );
 422
 423         int iNext = helper.WrapIndex( iNode+1, nbCorners );
 424         int iPrev = helper.WrapIndex( iNode-1, nbCorners );
 425
 426         SMESH_TNodeXYZ prevNode( newNodes[iPrev] );
 427         SMESH_TNodeXYZ curNode ( newNodes[iNode] );
 428         SMESH_TNodeXYZ nextNode( newNodes[iNext] );
 429         gp_Vec n1n0( prevNode - curNode);
 430         gp_Vec n1n2( nextNode - curNode );
 431         gp_Vec meshNorm = n1n2 ^ n1n0;
 432
 433         if ( geomNorm * meshNorm < 0 )
 434           SMDS_MeshCell::applyInterlace
 435             ( SMDS_MeshCell::reverseSmdsOrder( face->GetEntityType() ), newNodes );
 436       }
 437
 438       // make a new face
 439       if ( face->IsPoly() )
 440         newFace = tgtMesh->AddPolygonalFace( newNodes );
 441       else
 442         switch ( newNodes.size() )
 443         {
 444         case 3:
 445           newFace = tgtMesh->AddFace( newNodes[0], newNodes[1], newNodes[2] );
 446           break;
 447         case 4:
 448           newFace = tgtMesh->AddFace( newNodes[0], newNodes[1], newNodes[2], newNodes[3] );
 449           break;
 450         case 6:
 451           newFace = tgtMesh->AddFace( newNodes[0], newNodes[1], newNodes[2],
 452                                       newNodes[3], newNodes[4], newNodes[5]);
 453           break;
 454         case 8:
 455           newFace = tgtMesh->AddFace( newNodes[0], newNodes[1], newNodes[2], newNodes[3],
 456                                       newNodes[4], newNodes[5], newNodes[6], newNodes[7]);
 457           break;
 458         default: continue;
 459         }
 460       tgtMesh->SetMeshElementOnShape( newFace, shapeID );
 461       e2e->insert( make_pair( face, newFace ));
 462
 463       // collect links
 464       const SMDS_MeshNode* medium = 0;
 465       for ( int i = 0; i < nbCorners; ++i )
 466       {
 467         const SMDS_MeshNode* n1 = newNodes[i];
 468         const SMDS_MeshNode* n2 = newNodes[ (i+1)%nbCorners ];
 469         if ( isQuad ) // quadratic face
 470           medium = newNodes[i+nbCorners];
 471         link2Nb = linkCount.insert( make_pair( TLink( n1, n2, medium ), 0)).first;
 472         ++link2Nb->second;
 473         // if ( link2Nb->second == 1 )
 474         // {
 475         //   // measure link length
 476         //   double len2 = SMESH_TNodeXYZ( n1 ).SquareDistance( n2 );
 477         //   if ( len2 < minGroupTol )
 478         //     minGroupTol = len2;
 479         // }
 480       }
 481     }
 482     // Remove OUT nodes from n2n map
 483     for ( n2nIt = n2n->begin(); n2nIt != n2n->end(); )
 484       if ( !n2nIt->second )
 485         n2n->erase( n2nIt++ );
 486       else
 487         ++n2nIt;
 488   }
 489
 490
 491   // ==========================================================
 492   // Put nodes on geom edges and create edges on them;
 493   // check if the whole geom face is covered by imported faces
 494   // ==========================================================
 495
 496   // use large tolerance for projection of nodes to edges because of
 497   // BLSURF mesher specifics (issue 0020918, Study2.hdf)
 498   const double projTol = minGroupTol;
 499
 500   bool isFaceMeshed = false;
 501   SMESHDS_SubMesh* tgtFaceSM = tgtMesh->MeshElements( theShape );
 502   if ( tgtFaceSM )
 503   {
 504     // the imported mesh is valid if all external links (encountered once)
 505     // lie on geom edges
 506     subShapeIDs.erase( shapeID ); // to contain edges and vertices only
 507     double u, f, l;
 508     for ( link2Nb = linkCount.begin(); link2Nb != linkCount.end(); ++link2Nb)
 509     {
 510       const TLink& link = (*link2Nb).first;
 511       int nbFaces = link2Nb->second;
 512       if ( nbFaces == 1 )
 513       {
 514         // check if a not shared link lies on face boundary
 515         bool nodesOnBoundary = true;
 516         list< TopoDS_Shape > bndShapes;
 517         for ( int is1stN = 0; is1stN < 2 && nodesOnBoundary; ++is1stN )
 518         {
 519           const SMDS_MeshNode* n = is1stN ? link.node1() : link.node2();
 520           if ( !subShapeIDs.count( n->getshapeId() )) // n is assigned to FACE
 521           {
 522             for ( size_t iE = 0; iE < edges.size(); ++iE )
 523               if ( helper.CheckNodeU( edges[iE], n, u=0, projTol, /*force=*/true ))
 524               {
 525                 BRep_Tool::Range(edges[iE],f,l);
 526                 if ( Abs(u-f) < 2 * faceTol || Abs(u-l) < 2 * faceTol )
 527                   // duplicated node on vertex
 528                   return error("Source elements overlap one another");
 529                 tgtFaceSM->RemoveNode( n, /*isNodeDeleted=*/false );
 530                 tgtMesh->SetNodeOnEdge( n, edges[iE], u );
 531                 break;
 532               }
 533             nodesOnBoundary = subShapeIDs.count( n->getshapeId());
 534           }
 535           if ( nodesOnBoundary )
 536           {
 537             TopoDS_Shape s = helper.GetSubShapeByNode( n, tgtMesh );
 538             if ( s.ShapeType() == TopAbs_VERTEX )
 539               bndShapes.push_front( s ); // vertex first
 540             else
 541               bndShapes.push_back( s ); // edges last
 542           }
 543         }
 544         if ( !nodesOnBoundary )
 545         {
 546           error("free internal link"); // just for an easier debug
 547           break;
 548         }
 549         if ( bndShapes.front().ShapeType() == TopAbs_EDGE && // all link nodes are on EDGEs
 550              bndShapes.front() != bndShapes.back() )
 551           // link nodes on different geom edges
 552           return error(COMPERR_BAD_INPUT_MESH, "Source nodes mismatch target vertices");
 553
 554         // find geom edge the link is on
 555         if ( bndShapes.back().ShapeType() != TopAbs_EDGE ) // all link nodes are on VERTEXes
 556         {
 557           // find geom edge by two vertices
 558           TopoDS_Shape geomEdge = helper.GetCommonAncestor( bndShapes.back(),
 559                                                             bndShapes.front(),
 560                                                             theMesh, TopAbs_EDGE );
 561           if ( geomEdge.IsNull() )
 562           {
 563             error("free internal link");
 564             break; // vertices belong to different edges
 565           }
 566           bndShapes.push_back( geomEdge );
 567         }
 568
 569         // create an edge if not yet exists
 570         newNodes.resize(2);
 571         newNodes[0] = link.node1(), newNodes[1] = link.node2();
 572         const SMDS_MeshElement* edge = tgtMesh->FindElement( newNodes, SMDSAbs_Edge );
 573         if ( edge ) continue;
 574
 575         if ( link._reversed ) std::swap( newNodes[0], newNodes[1] );
 576         if ( link._medium )
 577         {
 578           edge = tgtMesh->AddEdge( newNodes[0], newNodes[1], link._medium );
 579
 580           TopoDS_Edge geomEdge = TopoDS::Edge(bndShapes.back());
 581           helper.CheckNodeU( geomEdge, link._medium, u, projTol, /*force=*/true );
 582           tgtFaceSM->RemoveNode( link._medium, /*isNodeDeleted=*/false );
 583           tgtMesh->SetNodeOnEdge( (SMDS_MeshNode*)link._medium, geomEdge, u );
 584         }
 585         else
 586         {
 587           edge = tgtMesh->AddEdge( newNodes[0], newNodes[1]);
 588         }
 589         if ( !edge )
 590           return false;
 591
 592         tgtMesh->SetMeshElementOnShape( edge, bndShapes.back() );
 593       }
 594       else if ( nbFaces > 2 )
 595       {
 596         return error( COMPERR_BAD_INPUT_MESH, "Non-manifold source mesh");
 597       }
 598     }
 599     isFaceMeshed = ( link2Nb == linkCount.end() && !linkCount.empty());
 600     if ( isFaceMeshed )
 601     {
 602       // check that source faces do not overlap:
 603       // there must be only two edges sharing each vertex and bound to sub-edges of theShape
 604       SMESH_MeshEditor editor( &theMesh );
 605       set<int>::iterator subID = subShapeIDs.begin();
 606       for ( ; subID != subShapeIDs.end(); ++subID )
 607       {
 608         const TopoDS_Shape& s = tgtMesh->IndexToShape( *subID );
 609         if ( s.ShapeType() != TopAbs_VERTEX ) continue;
 610         const SMDS_MeshNode* n = SMESH_Algo::VertexNode( TopoDS::Vertex(s), tgtMesh );
 611         SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator(SMDSAbs_Edge);
 612         int nbEdges = 0;
 613         while ( eIt->more() )
 614         {
 615           const SMDS_MeshElement* edge = eIt->next();
 616           int sId = editor.FindShape( edge );
 617           nbEdges += subShapeIDs.count( sId );
 618         }
 619         if ( nbEdges < 2 && !helper.IsRealSeam( s ))
 620           return false; // weird
 621         if ( nbEdges > 2 )
 622           return error( COMPERR_BAD_INPUT_MESH, "Source elements overlap one another");
 623       }
 624     }
 625   }
 626   if ( !isFaceMeshed )
 627     return error( COMPERR_BAD_INPUT_MESH,
 628                   "Source elements don't cover totally the geometrical face" );
 629
 630   if ( helper.HasSeam() )
 631   {
 632     // links on seam edges are shared by two faces, so no edges were created on them
 633     // by the previous detection of 2D mesh boundary
 634     for ( size_t iE = 0; iE < edges.size(); ++iE )
 635     {
 636       if ( !helper.IsRealSeam( edges[iE] )) continue;
 637       const TopoDS_Edge& seamEdge = edges[iE];
 638       // to find nodes lying on the seamEdge we check nodes of mesh faces sharing a node on one
 639       // of its vertices; after finding another node on seamEdge we continue the same way
 640       // until finding all nodes.
 641       TopoDS_Vertex      seamVertex = helper.IthVertex( 0, seamEdge );
 642       const SMDS_MeshNode* vertNode = SMESH_Algo::VertexNode( seamVertex, tgtMesh );
 643       set< const SMDS_MeshNode* > checkedNodes; checkedNodes.insert( vertNode );
 644       set< const SMDS_MeshElement* > checkedFaces;
 645       // as a face can have more than one node on the seamEdge, there is a difficulty in selecting
 646       // one of those nodes to treat next; so we simply find all nodes on the seamEdge and
 647       // then sort them by U on edge
 648       typedef list< pair< double, const SMDS_MeshNode* > > TUNodeList;
 649       TUNodeList nodesOnSeam;
 650       double u = helper.GetNodeU( seamEdge, vertNode );
 651       nodesOnSeam.push_back( make_pair( u, vertNode ));
 652       TUNodeList::iterator u2nIt = nodesOnSeam.begin();
 653       for ( ; u2nIt != nodesOnSeam.end(); ++u2nIt )
 654       {
 655         const SMDS_MeshNode* startNode = (*u2nIt).second;
 656         SMDS_ElemIteratorPtr faceIt = startNode->GetInverseElementIterator( SMDSAbs_Face );
 657         while ( faceIt->more() )
 658         {
 659           const SMDS_MeshElement* face = faceIt->next();
 660           if ( !checkedFaces.insert( face ).second ) continue;
 661           for ( int i = 0, nbNodes = face->NbCornerNodes(); i < nbNodes; ++i )
 662           {
 663             const SMDS_MeshNode* n = face->GetNode( i );
 664             if ( n == startNode || !checkedNodes.insert( n ).second ) continue;
 665             if ( helper.CheckNodeU( seamEdge, n, u=0, projTol, /*force=*/true ))
 666               nodesOnSeam.push_back( make_pair( u, n ));
 667           }
 668         }
 669       }
 670       // sort the found nodes by U on the seamEdge; most probably they are in a good order,
 671       // so we can use the hint to spead-up map filling
 672       map< double, const SMDS_MeshNode* > u2nodeMap;
 673       for ( u2nIt = nodesOnSeam.begin(); u2nIt != nodesOnSeam.end(); ++u2nIt )
 674         u2nodeMap.insert( u2nodeMap.end(), *u2nIt );
 675
 676       // create edges
 677       {
 678         SMESH_MesherHelper seamHelper( theMesh );
 679         seamHelper.SetSubShape( edges[ iE ]);
 680         seamHelper.SetElementsOnShape( true );
 681
 682         if ( (*checkedFaces.begin())->IsQuadratic() )
 683           for ( set< const SMDS_MeshElement* >::iterator fIt = checkedFaces.begin();
 684                 fIt != checkedFaces.end(); ++fIt )
 685             seamHelper.AddTLinks( static_cast<const SMDS_MeshFace*>( *fIt ));
 686
 687         map< double, const SMDS_MeshNode* >::iterator n1, n2, u2nEnd = u2nodeMap.end();
 688         for ( n2 = u2nodeMap.begin(), n1 = n2++; n2 != u2nEnd; ++n1, ++n2 )
 689         {
 690           const SMDS_MeshNode* node1 = n1->second;
 691           const SMDS_MeshNode* node2 = n2->second;
 692           seamHelper.AddEdge( node1, node2 );
 693           if ( node2->getshapeId() == helper.GetSubShapeID() )
 694           {
 695             tgtFaceSM->RemoveNode( node2, /*isNodeDeleted=*/false );
 696             tgtMesh->SetNodeOnEdge( const_cast<SMDS_MeshNode*>( node2 ), seamEdge, n2->first );
 697           }
 698         }
 699       }
 700     } // loop on edges to find seam ones
 701   } // if ( helper.HasSeam() )
 702
 703   // notify sub-meshes of edges on computation
 704   for ( size_t iE = 0; iE < edges.size(); ++iE )
 705   {
 706     SMESH_subMesh * sm = theMesh.GetSubMesh( edges[iE] );
 707     // if ( SMESH_Algo::isDegenerated( edges[iE] ))
 708     //   sm->SetIsAlwaysComputed( true );
 709     sm->ComputeStateEngine(SMESH_subMesh::CHECK_COMPUTE_STATE);
 710     if ( sm->GetComputeState() != SMESH_subMesh::COMPUTE_OK )
 711       return error(SMESH_Comment("Failed to create segments on the edge ")
 712                    << tgtMesh->ShapeToIndex( edges[iE ]));
 713   }
 714
 715   // ============
 716   // Copy meshes
 717   // ============
 718
 719   vector<SMESH_Mesh*> srcMeshes = _sourceHyp->GetSourceMeshes();
 720   for ( size_t i = 0; i < srcMeshes.size(); ++i )
 721     StdMeshers_Import_1D::importMesh( srcMeshes[i], theMesh, _sourceHyp, theShape );
 722
 723   return true;
 724 }
 725
 726 //=============================================================================
 727 /*!
 728  * \brief Set needed event listeners and create a submesh for a copied mesh
 729  *
 730  * This method is called only if a submesh has HYP_OK algo_state.
 731  */
 732 //=============================================================================
 733
 734 void StdMeshers_Import_1D2D::SetEventListener(SMESH_subMesh* subMesh)
 735 {
 736   if ( !_sourceHyp )
 737   {
 738     const TopoDS_Shape& tgtShape = subMesh->GetSubShape();
 739     SMESH_Mesh*         tgtMesh  = subMesh->GetFather();
 740     Hypothesis_Status aStatus;
 741     CheckHypothesis( *tgtMesh, tgtShape, aStatus );
 742   }
 743   StdMeshers_Import_1D::setEventListener( subMesh, _sourceHyp );
 744 }
 745 void StdMeshers_Import_1D2D::SubmeshRestored(SMESH_subMesh* subMesh)
 746 {
 747   SetEventListener(subMesh);
 748 }
 749
 750 //=============================================================================
 751 /*!
 752  * Predict nb of mesh entities created by Compute()
 753  */
 754 //=============================================================================
 755
 756 bool StdMeshers_Import_1D2D::Evaluate(SMESH_Mesh &         theMesh,
 757                                       const TopoDS_Shape & theShape,
 758                                       MapShapeNbElems&     aResMap)
 759 {
 760   if ( !_sourceHyp ) return false;
 761
 762   const vector<SMESH_Group*>& srcGroups = _sourceHyp->GetGroups();
 763   if ( srcGroups.empty() )
 764     return error("Invalid source groups");
 765
 766   vector<int> aVec(SMDSEntity_Last,0);
 767
 768   bool toCopyMesh, toCopyGroups;
 769   _sourceHyp->GetCopySourceMesh(toCopyMesh, toCopyGroups);
 770   if ( toCopyMesh ) // the whole mesh is copied
 771   {
 772     vector<SMESH_Mesh*> srcMeshes = _sourceHyp->GetSourceMeshes();
 773     for ( unsigned i = 0; i < srcMeshes.size(); ++i )
 774     {
 775       SMESH_subMesh* sm = StdMeshers_Import_1D::getSubMeshOfCopiedMesh( theMesh, *srcMeshes[i]);
 776       if ( !sm || aResMap.count( sm )) continue; // already counted
 777       const SMDS_MeshInfo& aMeshInfo = srcMeshes[i]->GetMeshDS()->GetMeshInfo();
 778       for (int i = 0; i < SMDSEntity_Last; i++)
 779         aVec[i] = aMeshInfo.NbEntities((SMDSAbs_EntityType)i);
 780     }
 781   }
 782   else
 783   {
 784     // std-like iterator used to get coordinates of nodes of mesh element
 785     typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
 786
 787     SMESH_MesherHelper helper(theMesh);
 788     helper.SetSubShape(theShape);
 789
 790     const TopoDS_Face& geomFace = TopoDS::Face( theShape );
 791
 792     // take into account nodes on vertices
 793     TopExp_Explorer exp( theShape, TopAbs_VERTEX );
 794     for ( ; exp.More(); exp.Next() )
 795       theMesh.GetSubMesh( exp.Current())->Evaluate( aResMap );
 796
 797     // to count now many times a link between nodes encounters,
 798     // negative nb additionally means that a link is quadratic
 799     map<SMESH_TLink, int> linkCount;
 800     map<SMESH_TLink, int>::iterator link2Nb;
 801
 802     // count faces and nodes imported from groups
 803     set<const SMDS_MeshNode* > allNodes;
 804     gp_XY uv;
 805     double minGroupTol = 1e100;
 806     for ( int iG = 0; iG < srcGroups.size(); ++iG )
 807     {
 808       const SMESHDS_GroupBase* srcGroup = srcGroups[iG]->GetGroupDS();
 809       const double groupTol = 0.5 * sqrt( getMinElemSize2( srcGroup ));
 810       minGroupTol = std::min( groupTol, minGroupTol );
 811       SMDS_ElemIteratorPtr srcElems = srcGroup->GetElements();
 812       SMDS_MeshNode *tmpNode =helper.AddNode(0,0,0);
 813       while ( srcElems->more() ) // loop on group contents
 814       {
 815         const SMDS_MeshElement* face = srcElems->next();
 816         // find out if face is located on geomEdge by projecting
 817         // a gravity center of face to geomFace
 818         gp_XYZ gc(0,0,0);
 819         gc = accumulate( TXyzIterator(face->nodesIterator()), TXyzIterator(), gc)/face->NbNodes();
 820         tmpNode->setXYZ( gc.X(), gc.Y(), gc.Z());
 821         if ( helper.CheckNodeUV( geomFace, tmpNode, uv, groupTol, /*force=*/true ))
 822         {
 823           ++aVec[ face->GetEntityType() ];
 824
 825           // collect links
 826           int nbConers = face->NbCornerNodes();
 827           for ( int i = 0; i < face->NbNodes(); ++i )
 828           {
 829             const SMDS_MeshNode* n1 = face->GetNode(i);
 830             allNodes.insert( n1 );
 831             if ( i < nbConers )
 832             {
 833               const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbConers );
 834               link2Nb = linkCount.insert( make_pair( SMESH_TLink( n1, n2 ), 0)).first;
 835               if ( (*link2Nb).second )
 836                 link2Nb->second += (link2Nb->second < 0 ) ? -1 : 1;
 837               else
 838                 link2Nb->second += ( face->IsQuadratic() ) ? -1 : 1;
 839             }
 840           }
 841         }
 842       }
 843       helper.GetMeshDS()->RemoveNode(tmpNode);
 844     }
 845
 846     int nbNodes = allNodes.size();
 847     allNodes.clear();
 848
 849     // count nodes and edges on geom edges
 850
 851     double u;
 852     for ( exp.Init(theShape, TopAbs_EDGE); exp.More(); exp.Next() )
 853     {
 854       TopoDS_Edge geomEdge = TopoDS::Edge( exp.Current() );
 855       SMESH_subMesh* sm = theMesh.GetSubMesh( geomEdge );
 856       vector<int>& edgeVec = aResMap[sm];
 857       if ( edgeVec.empty() )
 858       {
 859         edgeVec.resize(SMDSEntity_Last,0);
 860         for ( link2Nb = linkCount.begin(); link2Nb != linkCount.end(); )
 861         {
 862           const SMESH_TLink& link = (*link2Nb).first;
 863           int nbFacesOfLink = Abs( link2Nb->second );
 864           bool eraseLink = ( nbFacesOfLink != 1 );
 865           if ( nbFacesOfLink == 1 )
 866           {
 867             if ( helper.CheckNodeU( geomEdge, link.node1(), u, minGroupTol, /*force=*/true )&&
 868                  helper.CheckNodeU( geomEdge, link.node2(), u, minGroupTol, /*force=*/true ))
 869             {
 870               bool isQuadratic = ( link2Nb->second < 0 );
 871               ++edgeVec[ isQuadratic ? SMDSEntity_Quad_Edge : SMDSEntity_Edge ];
 872               ++edgeVec[ SMDSEntity_Node ];
 873               --nbNodes;
 874               eraseLink = true;
 875             }
 876           }
 877           if ( eraseLink )
 878             linkCount.erase(link2Nb++);
 879           else
 880             link2Nb++;
 881         }
 882         if ( edgeVec[ SMDSEntity_Node] > 0 )
 883           --edgeVec[ SMDSEntity_Node ]; // for one node on vertex
 884       }
 885       else if ( !helper.IsSeamShape( geomEdge ) ||
 886                 geomEdge.Orientation() == TopAbs_FORWARD )
 887       {
 888         nbNodes -= 1+edgeVec[ SMDSEntity_Node ];
 889       }
 890     }
 891
 892     aVec[SMDSEntity_Node] = nbNodes;
 893   }
 894
 895   SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
 896   aResMap.insert(make_pair(sm,aVec));
 897
 898   return true;
 899 }