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         newNodes [ i ] = 0;
 311
 312         it_isnew = n2n->insert( make_pair( *node, (SMDS_MeshNode*)0 ));
 313         n2nIt    = it_isnew.first;
 314
 315         const SMDS_MeshNode* & newNode = n2nIt->second;
 316         if ( !it_isnew.second && !newNode )
 317           break; // a node is mapped to NULL - it is OUT of the FACE
 318
 319         if ( newNode )
 320         {
 321           if ( !subShapeIDs.count( newNode->getshapeId() ))
 322             break; // node is Imported onto other FACE
 323           if ( !isIn && bndNodes.count( *node ))
 324             nodeState[ i ] = TopAbs_ON;
 325         }
 326         else
 327         {
 328           // find a node pre-existing on EDGE or VERTEX
 329           dist2foundNodes.clear();
 330           existingNodeOcTr.NodesAround( nXYZ, dist2foundNodes, groupTol );
 331           if ( !dist2foundNodes.empty() )
 332           {
 333             newNode = dist2foundNodes.begin()->second;
 334             nodeState[ i ] = TopAbs_ON;
 335           }
 336         }
 337
 338         if ( !newNode )
 339         {
 340           // find out if node lies on the surface of theShape
 341           gp_XY uv( Precision::Infinite(), 0 );
 342           isOut = ( !helper.CheckNodeUV( geomFace, *node, uv, groupTol, /*force=*/true ) ||
 343                     bndBox2d.IsOut( uv ));
 344           if ( !isOut && !isIn ) // classify
 345           {
 346             classifier.Perform( geomFace, uv, clsfTol );
 347             nodeState[i] = classifier.State();
 348             isOut = ( nodeState[i] == TopAbs_OUT );
 349           }
 350           if ( !isOut ) // create a new node
 351           {
 352             newNode = tgtMesh->AddNode( nXYZ.X(), nXYZ.Y(), nXYZ.Z());
 353             tgtMesh->SetNodeOnFace( newNode, shapeID, uv.X(), uv.Y() );
 354             nbCreatedNodes++;
 355             if ( nodeState[i] == TopAbs_ON )
 356               bndNodes.insert( *node );
 357             else
 358               isIn = true;
 359           }
 360         }
 361         if ( !(newNodes[i] = newNode ) || isOut )
 362           break;
 363       }
 364
 365       if ( !newNodes.back() )
 366         continue; // not all nodes of the face lie on theShape
 367
 368       if ( !isIn ) // if all nodes are on FACE boundary, a mesh face can be OUT
 369       {
 370         // check state of nodes created for other faces
 371         for ( size_t i = 0; i < nodeState.size() && !isIn; ++i )
 372         {
 373           if ( nodeState[i] != TopAbs_UNKNOWN ) continue;
 374           gp_XY uv = helper.GetNodeUV( geomFace, newNodes[i] );
 375           classifier.Perform( geomFace, uv, clsfTol );
 376           nodeState[i] = classifier.State();
 377           isIn = ( nodeState[i] == TopAbs_IN );
 378         }
 379         if ( !isIn ) // classify face center
 380         {
 381           gp_XYZ gc( 0., 0., 0 );
 382           for ( size_t i = 0; i < newNodes.size(); ++i )
 383             gc += SMESH_TNodeXYZ( newNodes[i] );
 384           gc /= newNodes.size();
 385
 386           TopLoc_Location loc;
 387           GeomAPI_ProjectPointOnSurf& proj = helper.GetProjector( geomFace,
 388                                                                   loc,
 389                                                                   helper.MaxTolerance( geomFace ));
 390           if ( !loc.IsIdentity() ) loc.Transformation().Inverted().Transforms( gc );
 391           proj.Perform( gc );
 392           if ( !proj.IsDone() || proj.NbPoints() < 1 )
 393             continue;
 394           Quantity_Parameter U,V;
 395           proj.LowerDistanceParameters(U,V);
 396           gp_XY uv( U,V );
 397           classifier.Perform( geomFace, uv, clsfTol );
 398           if ( classifier.State() != TopAbs_IN )
 399             continue;
 400         }
 401       }
 402
 403       // try to find already created face
 404       SMDS_MeshElement * newFace = 0;
 405       if ( nbCreatedNodes == 0 &&
 406            tgtMesh->FindElement(newNodes, SMDSAbs_Face, /*noMedium=*/false))
 407         continue; // repeated face in source groups already created
 408
 409       // check future face orientation
 410       const int nbCorners = face->NbCornerNodes();
 411       const bool isQuad   = ( nbCorners != (int) newNodes.size() );
 412       if ( toCheckOri )
 413       {
 414         int iNode = -1;
 415         gp_Vec geomNorm;
 416         do
 417         {
 418           gp_XY uv = helper.GetNodeUV( geomFace, newNodes[++iNode] );
 419           surface->D1( uv.X(),uv.Y(), p, du,dv );
 420           geomNorm = reverse ? dv^du : du^dv;
 421         }
 422         while ( geomNorm.SquareMagnitude() < 1e-6 && iNode+1 < nbCorners );
 423
 424         int iNext = helper.WrapIndex( iNode+1, nbCorners );
 425         int iPrev = helper.WrapIndex( iNode-1, nbCorners );
 426
 427         SMESH_TNodeXYZ prevNode( newNodes[iPrev] );
 428         SMESH_TNodeXYZ curNode ( newNodes[iNode] );
 429         SMESH_TNodeXYZ nextNode( newNodes[iNext] );
 430         gp_Vec n1n0( prevNode - curNode);
 431         gp_Vec n1n2( nextNode - curNode );
 432         gp_Vec meshNorm = n1n2 ^ n1n0;
 433
 434         if ( geomNorm * meshNorm < 0 )
 435           SMDS_MeshCell::applyInterlace
 436             ( SMDS_MeshCell::reverseSmdsOrder( face->GetEntityType() ), newNodes );
 437       }
 438
 439       // make a new face
 440       if ( face->IsPoly() )
 441         newFace = tgtMesh->AddPolygonalFace( newNodes );
 442       else
 443         switch ( newNodes.size() )
 444         {
 445         case 3:
 446           newFace = tgtMesh->AddFace( newNodes[0], newNodes[1], newNodes[2] );
 447           break;
 448         case 4:
 449           newFace = tgtMesh->AddFace( newNodes[0], newNodes[1], newNodes[2], newNodes[3] );
 450           break;
 451         case 6:
 452           newFace = tgtMesh->AddFace( newNodes[0], newNodes[1], newNodes[2],
 453                                       newNodes[3], newNodes[4], newNodes[5]);
 454           break;
 455         case 8:
 456           newFace = tgtMesh->AddFace( newNodes[0], newNodes[1], newNodes[2], newNodes[3],
 457                                       newNodes[4], newNodes[5], newNodes[6], newNodes[7]);
 458           break;
 459         default: continue;
 460         }
 461       tgtMesh->SetMeshElementOnShape( newFace, shapeID );
 462       e2e->insert( make_pair( face, newFace ));
 463
 464       // collect links
 465       const SMDS_MeshNode* medium = 0;
 466       for ( int i = 0; i < nbCorners; ++i )
 467       {
 468         const SMDS_MeshNode* n1 = newNodes[i];
 469         const SMDS_MeshNode* n2 = newNodes[ (i+1)%nbCorners ];
 470         if ( isQuad ) // quadratic face
 471           medium = newNodes[i+nbCorners];
 472         link2Nb = linkCount.insert( make_pair( TLink( n1, n2, medium ), 0)).first;
 473         ++link2Nb->second;
 474         // if ( link2Nb->second == 1 )
 475         // {
 476         //   // measure link length
 477         //   double len2 = SMESH_TNodeXYZ( n1 ).SquareDistance( n2 );
 478         //   if ( len2 < minGroupTol )
 479         //     minGroupTol = len2;
 480         // }
 481       }
 482     }
 483     // Remove OUT nodes from n2n map
 484     for ( n2nIt = n2n->begin(); n2nIt != n2n->end(); )
 485       if ( !n2nIt->second )
 486         n2n->erase( n2nIt++ );
 487       else
 488         ++n2nIt;
 489   }
 490
 491
 492   // ==========================================================
 493   // Put nodes on geom edges and create edges on them;
 494   // check if the whole geom face is covered by imported faces
 495   // ==========================================================
 496
 497   // use large tolerance for projection of nodes to edges because of
 498   // BLSURF mesher specifics (issue 0020918, Study2.hdf)
 499   const double projTol = minGroupTol;
 500
 501   bool isFaceMeshed = false;
 502   SMESHDS_SubMesh* tgtFaceSM = tgtMesh->MeshElements( theShape );
 503   if ( tgtFaceSM )
 504   {
 505     // the imported mesh is valid if all external links (encountered once)
 506     // lie on geom edges
 507     subShapeIDs.erase( shapeID ); // to contain edges and vertices only
 508     double u, f, l;
 509     for ( link2Nb = linkCount.begin(); link2Nb != linkCount.end(); ++link2Nb)
 510     {
 511       const TLink& link = (*link2Nb).first;
 512       int nbFaces = link2Nb->second;
 513       if ( nbFaces == 1 )
 514       {
 515         // check if a not shared link lies on face boundary
 516         bool nodesOnBoundary = true;
 517         list< TopoDS_Shape > bndShapes;
 518         for ( int is1stN = 0; is1stN < 2 && nodesOnBoundary; ++is1stN )
 519         {
 520           const SMDS_MeshNode* n = is1stN ? link.node1() : link.node2();
 521           if ( !subShapeIDs.count( n->getshapeId() )) // n is assigned to FACE
 522           {
 523             for ( size_t iE = 0; iE < edges.size(); ++iE )
 524               if ( helper.CheckNodeU( edges[iE], n, u=0, projTol, /*force=*/true ))
 525               {
 526                 BRep_Tool::Range(edges[iE],f,l);
 527                 if ( Abs(u-f) < 2 * faceTol || Abs(u-l) < 2 * faceTol )
 528                   // duplicated node on vertex
 529                   return error("Source elements overlap one another");
 530                 tgtFaceSM->RemoveNode( n, /*isNodeDeleted=*/false );
 531                 tgtMesh->SetNodeOnEdge( n, edges[iE], u );
 532                 break;
 533               }
 534             nodesOnBoundary = subShapeIDs.count( n->getshapeId());
 535           }
 536           if ( nodesOnBoundary )
 537           {
 538             TopoDS_Shape s = helper.GetSubShapeByNode( n, tgtMesh );
 539             if ( s.ShapeType() == TopAbs_VERTEX )
 540               bndShapes.push_front( s ); // vertex first
 541             else
 542               bndShapes.push_back( s ); // edges last
 543           }
 544         }
 545         if ( !nodesOnBoundary )
 546         {
 547           error("free internal link"); // just for an easier debug
 548           break;
 549         }
 550         if ( bndShapes.front().ShapeType() == TopAbs_EDGE && // all link nodes are on EDGEs
 551              bndShapes.front() != bndShapes.back() )
 552           // link nodes on different geom edges
 553           return error(COMPERR_BAD_INPUT_MESH, "Source nodes mismatch target vertices");
 554
 555         // find geom edge the link is on
 556         if ( bndShapes.back().ShapeType() != TopAbs_EDGE ) // all link nodes are on VERTEXes
 557         {
 558           // find geom edge by two vertices
 559           TopoDS_Shape geomEdge = helper.GetCommonAncestor( bndShapes.back(),
 560                                                             bndShapes.front(),
 561                                                             theMesh, TopAbs_EDGE );
 562           if ( geomEdge.IsNull() )
 563           {
 564             error("free internal link");
 565             break; // vertices belong to different edges
 566           }
 567           bndShapes.push_back( geomEdge );
 568         }
 569
 570         // create an edge if not yet exists
 571         newNodes.resize(2);
 572         newNodes[0] = link.node1(), newNodes[1] = link.node2();
 573         const SMDS_MeshElement* edge = tgtMesh->FindElement( newNodes, SMDSAbs_Edge );
 574         if ( edge ) continue;
 575
 576         if ( link._reversed ) std::swap( newNodes[0], newNodes[1] );
 577         if ( link._medium )
 578         {
 579           edge = tgtMesh->AddEdge( newNodes[0], newNodes[1], link._medium );
 580
 581           TopoDS_Edge geomEdge = TopoDS::Edge(bndShapes.back());
 582           helper.CheckNodeU( geomEdge, link._medium, u, projTol, /*force=*/true );
 583           tgtFaceSM->RemoveNode( link._medium, /*isNodeDeleted=*/false );
 584           tgtMesh->SetNodeOnEdge( (SMDS_MeshNode*)link._medium, geomEdge, u );
 585         }
 586         else
 587         {
 588           edge = tgtMesh->AddEdge( newNodes[0], newNodes[1]);
 589         }
 590         if ( !edge )
 591           return false;
 592
 593         tgtMesh->SetMeshElementOnShape( edge, bndShapes.back() );
 594       }
 595       else if ( nbFaces > 2 )
 596       {
 597         return error( COMPERR_BAD_INPUT_MESH, "Non-manifold source mesh");
 598       }
 599     }
 600     isFaceMeshed = ( link2Nb == linkCount.end() && !linkCount.empty());
 601     if ( isFaceMeshed )
 602     {
 603       // check that source faces do not overlap:
 604       // there must be only two edges sharing each vertex and bound to sub-edges of theShape
 605       SMESH_MeshEditor editor( &theMesh );
 606       set<int>::iterator subID = subShapeIDs.begin();
 607       for ( ; subID != subShapeIDs.end(); ++subID )
 608       {
 609         const TopoDS_Shape& s = tgtMesh->IndexToShape( *subID );
 610         if ( s.ShapeType() != TopAbs_VERTEX ) continue;
 611         const SMDS_MeshNode* n = SMESH_Algo::VertexNode( TopoDS::Vertex(s), tgtMesh );
 612         SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator(SMDSAbs_Edge);
 613         int nbEdges = 0;
 614         while ( eIt->more() )
 615         {
 616           const SMDS_MeshElement* edge = eIt->next();
 617           int sId = editor.FindShape( edge );
 618           nbEdges += subShapeIDs.count( sId );
 619         }
 620         if ( nbEdges < 2 && !helper.IsRealSeam( s ))
 621           return false; // weird
 622         if ( nbEdges > 2 )
 623           return error( COMPERR_BAD_INPUT_MESH, "Source elements overlap one another");
 624       }
 625     }
 626   }
 627   if ( !isFaceMeshed )
 628     return error( COMPERR_BAD_INPUT_MESH,
 629                   "Source elements don't cover totally the geometrical face" );
 630
 631   if ( helper.HasSeam() )
 632   {
 633     // links on seam edges are shared by two faces, so no edges were created on them
 634     // by the previous detection of 2D mesh boundary
 635     for ( size_t iE = 0; iE < edges.size(); ++iE )
 636     {
 637       if ( !helper.IsRealSeam( edges[iE] )) continue;
 638       const TopoDS_Edge& seamEdge = edges[iE];
 639       // to find nodes lying on the seamEdge we check nodes of mesh faces sharing a node on one
 640       // of its vertices; after finding another node on seamEdge we continue the same way
 641       // until finding all nodes.
 642       TopoDS_Vertex      seamVertex = helper.IthVertex( 0, seamEdge );
 643       const SMDS_MeshNode* vertNode = SMESH_Algo::VertexNode( seamVertex, tgtMesh );
 644       set< const SMDS_MeshNode* > checkedNodes; checkedNodes.insert( vertNode );
 645       set< const SMDS_MeshElement* > checkedFaces;
 646       // as a face can have more than one node on the seamEdge, there is a difficulty in selecting
 647       // one of those nodes to treat next; so we simply find all nodes on the seamEdge and
 648       // then sort them by U on edge
 649       typedef list< pair< double, const SMDS_MeshNode* > > TUNodeList;
 650       TUNodeList nodesOnSeam;
 651       double u = helper.GetNodeU( seamEdge, vertNode );
 652       nodesOnSeam.push_back( make_pair( u, vertNode ));
 653       TUNodeList::iterator u2nIt = nodesOnSeam.begin();
 654       for ( ; u2nIt != nodesOnSeam.end(); ++u2nIt )
 655       {
 656         const SMDS_MeshNode* startNode = (*u2nIt).second;
 657         SMDS_ElemIteratorPtr faceIt = startNode->GetInverseElementIterator( SMDSAbs_Face );
 658         while ( faceIt->more() )
 659         {
 660           const SMDS_MeshElement* face = faceIt->next();
 661           if ( !checkedFaces.insert( face ).second ) continue;
 662           for ( int i = 0, nbNodes = face->NbCornerNodes(); i < nbNodes; ++i )
 663           {
 664             const SMDS_MeshNode* n = face->GetNode( i );
 665             if ( n == startNode || !checkedNodes.insert( n ).second ) continue;
 666             if ( helper.CheckNodeU( seamEdge, n, u=0, projTol, /*force=*/true ))
 667               nodesOnSeam.push_back( make_pair( u, n ));
 668           }
 669         }
 670       }
 671       // sort the found nodes by U on the seamEdge; most probably they are in a good order,
 672       // so we can use the hint to spead-up map filling
 673       map< double, const SMDS_MeshNode* > u2nodeMap;
 674       for ( u2nIt = nodesOnSeam.begin(); u2nIt != nodesOnSeam.end(); ++u2nIt )
 675         u2nodeMap.insert( u2nodeMap.end(), *u2nIt );
 676
 677       // create edges
 678       {
 679         SMESH_MesherHelper seamHelper( theMesh );
 680         seamHelper.SetSubShape( edges[ iE ]);
 681         seamHelper.SetElementsOnShape( true );
 682
 683         if ( !checkedFaces.empty() && (*checkedFaces.begin())->IsQuadratic() )
 684           for ( set< const SMDS_MeshElement* >::iterator fIt = checkedFaces.begin();
 685                 fIt != checkedFaces.end(); ++fIt )
 686             seamHelper.AddTLinks( static_cast<const SMDS_MeshFace*>( *fIt ));
 687
 688         map< double, const SMDS_MeshNode* >::iterator n1, n2, u2nEnd = u2nodeMap.end();
 689         for ( n2 = u2nodeMap.begin(), n1 = n2++; n2 != u2nEnd; ++n1, ++n2 )
 690         {
 691           const SMDS_MeshNode* node1 = n1->second;
 692           const SMDS_MeshNode* node2 = n2->second;
 693           seamHelper.AddEdge( node1, node2 );
 694           if ( node2->getshapeId() == helper.GetSubShapeID() )
 695           {
 696             tgtFaceSM->RemoveNode( node2, /*isNodeDeleted=*/false );
 697             tgtMesh->SetNodeOnEdge( const_cast<SMDS_MeshNode*>( node2 ), seamEdge, n2->first );
 698           }
 699         }
 700       }
 701     } // loop on edges to find seam ones
 702   } // if ( helper.HasSeam() )
 703
 704   // notify sub-meshes of edges on computation
 705   for ( size_t iE = 0; iE < edges.size(); ++iE )
 706   {
 707     SMESH_subMesh * sm = theMesh.GetSubMesh( edges[iE] );
 708     // if ( SMESH_Algo::isDegenerated( edges[iE] ))
 709     //   sm->SetIsAlwaysComputed( true );
 710     sm->ComputeStateEngine(SMESH_subMesh::CHECK_COMPUTE_STATE);
 711     if ( sm->GetComputeState() != SMESH_subMesh::COMPUTE_OK )
 712       return error(SMESH_Comment("Failed to create segments on the edge #") << sm->GetId());
 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 }