Salome HOME
52499: Prismatic mesh is not computed on a prismatic shape
authoreap <eap@opencascade.com>
Wed, 24 Sep 2014 13:38:27 +0000 (17:38 +0400)
committereap <eap@opencascade.com>
Wed, 24 Sep 2014 13:38:27 +0000 (17:38 +0400)
22701: EDF SMESH: Crash when creating mesh

 1) Search of transformation using a least-square approximation
 2) Compute internal prism nodes using the found transformation

src/SMESHUtils/SMESH_TypeDefs.hxx
src/StdMeshers/StdMeshers_Prism_3D.cxx
src/StdMeshers/StdMeshers_Prism_3D.hxx
src/StdMeshers/StdMeshers_ProjectionUtils.cxx
src/StdMeshers/StdMeshers_ProjectionUtils.hxx
src/StdMeshers/StdMeshers_Projection_2D.cxx
src/StdMeshers/StdMeshers_Projection_2D.hxx
src/StdMeshers/StdMeshers_Projection_3D.cxx
src/StdMeshers/StdMeshers_RadialPrism_3D.cxx
src/StdMeshers/StdMeshers_ViscousLayers.cxx

index 9fec2e5..f5dedae 100644 (file)
@@ -43,7 +43,8 @@ typedef std::map<const SMDS_MeshElement*,
                  std::list<const SMDS_MeshElement*>, TIDCompare > TElemOfElemListMap;
 typedef std::map<const SMDS_MeshElement*,
                  std::list<const SMDS_MeshNode*>,    TIDCompare > TElemOfNodeListMap;
-typedef std::map<const SMDS_MeshNode*, const SMDS_MeshNode*>      TNodeNodeMap;
+typedef std::map<const SMDS_MeshNode*,
+                 const SMDS_MeshNode*,               TIDCompare>  TNodeNodeMap;
 
 //!< Set of elements sorted by ID, to be used to assure predictability of edition
 typedef std::set< const SMDS_MeshElement*, TIDCompare >      TIDSortedElemSet;
index afac9be..eb43e88 100644 (file)
@@ -77,7 +77,7 @@ using namespace std;
 #define SHOWYXZ(msg, xyz)
 #endif
 
-namespace TAssocTool = StdMeshers_ProjectionUtils;
+namespace NSProjUtils = StdMeshers_ProjectionUtils;
 
 typedef SMESH_Comment TCom;
 
@@ -157,6 +157,10 @@ namespace {
                                                             fatherAlgo->GetGen() );
       return algo;
     }
+    const NSProjUtils::TNodeNodeMap& GetNodesMap()
+    {
+      return _src2tgtNodes;
+    }
   };
   //=======================================================================
   /*!
@@ -539,7 +543,7 @@ StdMeshers_Prism_3D::StdMeshers_Prism_3D(int hypId, int studyId, SMESH_Gen* gen)
 {
   _name                    = "Prism_3D";
   _shapeType               = (1 << TopAbs_SOLID); // 1 bit per shape type
-  _onlyUnaryInput          = false; // accept all SOLIDs at once
+  _onlyUnaryInput          = false; // mesh all SOLIDs at once
   _requireDiscreteBoundary = false; // mesh FACEs and EDGEs by myself
   _supportSubmeshes        = true;  // "source" FACE must be meshed by other algo
   _neededLowerHyps[ 1 ]    = true;  // suppress warning on hiding a global 1D algo
@@ -578,12 +582,12 @@ bool StdMeshers_Prism_3D::CheckHypothesis(SMESH_Mesh&                          a
   for ( ; exp.More(); exp.Next() ) {
     ++nbFace;
     const TopoDS_Shape& face = exp.Current();
-    nbEdge = TAssocTool::Count( face, TopAbs_EDGE, 0 );
-    nbWire = TAssocTool::Count( face, TopAbs_WIRE, 0 );
+    nbEdge = NSProjUtils::Count( face, TopAbs_EDGE, 0 );
+    nbWire = NSProjUtils::Count( face, TopAbs_WIRE, 0 );
     if (  nbEdge!= 4 || nbWire!= 1 ) {
       if ( !notQuadFaces.empty() ) {
-        if ( TAssocTool::Count( notQuadFaces.back(), TopAbs_EDGE, 0 ) != nbEdge ||
-             TAssocTool::Count( notQuadFaces.back(), TopAbs_WIRE, 0 ) != nbWire )
+        if ( NSProjUtils::Count( notQuadFaces.back(), TopAbs_EDGE, 0 ) != nbEdge ||
+             NSProjUtils::Count( notQuadFaces.back(), TopAbs_WIRE, 0 ) != nbWire )
           RETURN_BAD_RESULT("Different not quad faces");
       }
       notQuadFaces.push_back( face );
@@ -595,7 +599,7 @@ bool StdMeshers_Prism_3D::CheckHypothesis(SMESH_Mesh&                          a
       RETURN_BAD_RESULT("Bad nb not quad faces: " << notQuadFaces.size());
 
     // check total nb faces
-    nbEdge = TAssocTool::Count( notQuadFaces.back(), TopAbs_EDGE, 0 );
+    nbEdge = NSProjUtils::Count( notQuadFaces.back(), TopAbs_EDGE, 0 );
     if ( nbFace != nbEdge + 2 )
       RETURN_BAD_RESULT("Bad nb of faces: " << nbFace << " but must be " << nbEdge + 2);
   }
@@ -1057,6 +1061,8 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
     return false;
 
   // Analyse mesh and geometry to find all block sub-shapes and submeshes
+  // (after fixing IPAL52499 myBlock is used only as a holder of boundary nodes
+  // and location of internal nodes is computed by StdMeshers_Sweeper)
   if ( !myBlock.Init( myHelper, thePrism ))
     return toSM( error( myBlock.GetError()));
 
@@ -1067,10 +1073,10 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
   // Try to get gp_Trsf to get all nodes from bottom ones
   vector<gp_Trsf> trsf;
   gp_Trsf bottomToTopTrsf;
-  if ( !myBlock.GetLayersTransformation( trsf, thePrism ))
-    trsf.clear();
-  else if ( !trsf.empty() )
-    bottomToTopTrsf = trsf.back();
+  // if ( !myBlock.GetLayersTransformation( trsf, thePrism ))
+  //   trsf.clear();
+  // else if ( !trsf.empty() )
+  //   bottomToTopTrsf = trsf.back();
 
   // To compute coordinates of a node inside a block, it is necessary to know
   // 1. normalized parameters of the node by which
@@ -1092,31 +1098,30 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
 
   // Create nodes inside the block
 
-  // try to use transformation (issue 0020680)
-  if ( !trsf.empty() )
+  // use transformation (issue 0020680, IPAL0052499)
+  StdMeshers_Sweeper sweeper;
+
+  // load boundary nodes
+  bool dummy;
+  list< TopoDS_Edge >::const_iterator edge = thePrism.myBottomEdges.begin();
+  for ( ; edge != thePrism.myBottomEdges.end(); ++edge )
   {
-    // loop on nodes inside the bottom face
-    TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin();
-    for ( ; bot_column != myBotToColumnMap.end(); ++bot_column )
-    {
-      const Prism_3D::TNode& tBotNode = bot_column->first; // bottom TNode
-      if ( tBotNode.GetPositionType() != SMDS_TOP_FACE )
-        continue; // node is not inside face 
+    int edgeID = meshDS->ShapeToIndex( *edge );
+    TParam2ColumnMap* u2col = const_cast<TParam2ColumnMap*>
+      ( myBlock.GetParam2ColumnMap( edgeID, dummy ));
+    TParam2ColumnMap::iterator u2colIt = u2col->begin();
+    for ( ; u2colIt != u2col->end(); ++u2colIt )
+      sweeper.myBndColumns.push_back( & u2colIt->second );
+  }
+  // load node columns inside the bottom face
+  TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin();
+  for ( ; bot_column != myBotToColumnMap.end(); ++bot_column )
+    sweeper.myIntColumns.push_back( & bot_column->second );
 
-      // column nodes; middle part of the column are zero pointers
-      TNodeColumn& column = bot_column->second;
-      TNodeColumn::iterator columnNodes = column.begin();
-      for ( int z = 0; columnNodes != column.end(); ++columnNodes, ++z)
-      {
-        const SMDS_MeshNode* & node = *columnNodes;
-        if ( node ) continue; // skip bottom or top node
+  const double tol = getSweepTolerance( thePrism );
 
-        gp_XYZ coords = tBotNode.GetCoords();
-        trsf[z-1].Transforms( coords );
-        node = meshDS->AddNode( coords.X(), coords.Y(), coords.Z() );
-        meshDS->SetNodeInVolume( node, volumeID );
-      }
-    } // loop on bottom nodes
+  if ( sweeper.ComputeNodes( *myHelper, tol ))
+  {
   }
   else // use block approach
   {
@@ -1127,7 +1132,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
     {
       const Prism_3D::TNode& tBotNode = bot_column->first; // bottom TNode
       if ( tBotNode.GetPositionType() != SMDS_TOP_FACE )
-        continue; // node is not inside the FACE 
+        continue; // node is not inside the FACE
 
       // column nodes; middle part of the column are zero pointers
       TNodeColumn& column = bot_column->second;
@@ -1853,8 +1858,8 @@ void StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
 bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf,
                                                  const Prism_3D::TPrismTopo& thePrism)
 {
-  SMESH_subMesh * botSM = myBlock.SubMesh( ID_BOT_FACE );
-  SMESH_subMesh * topSM = myBlock.SubMesh( ID_TOP_FACE );
+  SMESH_subMesh * botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom );
+  SMESH_subMesh * topSM = myHelper->GetMesh()->GetSubMesh( thePrism.myTop    );
 
   SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS();
   SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS();
@@ -1885,86 +1890,90 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf
                        <<" and #"<< topSM->GetId() << " seems different" ));
   ///RETURN_BAD_RESULT("Need to project but not allowed");
 
+  NSProjUtils::TNodeNodeMap n2nMap;
+  const NSProjUtils::TNodeNodeMap* n2nMapPtr = & n2nMap;
   if ( needProject )
   {
-    return projectBottomToTop( bottomToTopTrsf );
+    if ( !projectBottomToTop( bottomToTopTrsf, thePrism ))
+      return false;
+    n2nMapPtr = & TProjction2dAlgo::instance( this )->GetNodesMap();
   }
 
-  TopoDS_Face botFace = TopoDS::Face( myBlock.Shape( ID_BOT_FACE ));
-  TopoDS_Face topFace = TopoDS::Face( myBlock.Shape( ID_TOP_FACE ));
-  // associate top and bottom faces
-  TAssocTool::TShapeShapeMap shape2ShapeMap;
-  const bool sameTopo = 
-    TAssocTool::FindSubShapeAssociation( botFace, myBlock.Mesh(),
-                                         topFace, myBlock.Mesh(),
-                                         shape2ShapeMap);
-  if ( !sameTopo )
-    for ( size_t iQ = 0; iQ < thePrism.myWallQuads.size(); ++iQ )
-    {
-      const Prism_3D::TQuadList& quadList = thePrism.myWallQuads[iQ];
-      StdMeshers_FaceSidePtr      botSide = quadList.front()->side[ QUAD_BOTTOM_SIDE ];
-      StdMeshers_FaceSidePtr      topSide = quadList.back ()->side[ QUAD_TOP_SIDE ];
-      if ( botSide->NbEdges() == topSide->NbEdges() )
-      {
-        for ( int iE = 0; iE < botSide->NbEdges(); ++iE )
-        {
-          TAssocTool::InsertAssociation( botSide->Edge( iE ),
-                                         topSide->Edge( iE ), shape2ShapeMap );
-          TAssocTool::InsertAssociation( myHelper->IthVertex( 0, botSide->Edge( iE )),
-                                         myHelper->IthVertex( 0, topSide->Edge( iE )),
-                                         shape2ShapeMap );
-        }
-      }
-      else
+  if ( !n2nMapPtr || n2nMapPtr->size() < botSMDS->NbNodes() )
+  {
+    // associate top and bottom faces
+    NSProjUtils::TShapeShapeMap shape2ShapeMap;
+    const bool sameTopo =
+      NSProjUtils::FindSubShapeAssociation( thePrism.myBottom, myHelper->GetMesh(),
+                                           thePrism.myTop,    myHelper->GetMesh(),
+                                           shape2ShapeMap);
+    if ( !sameTopo )
+      for ( size_t iQ = 0; iQ < thePrism.myWallQuads.size(); ++iQ )
       {
-        TopoDS_Vertex vb, vt;
-        StdMeshers_FaceSidePtr sideB, sideT;
-        vb = myHelper->IthVertex( 0, botSide->Edge( 0 ));
-        vt = myHelper->IthVertex( 0, topSide->Edge( 0 ));
-        sideB = quadList.front()->side[ QUAD_LEFT_SIDE ];
-        sideT = quadList.back ()->side[ QUAD_LEFT_SIDE ];
-        if ( vb.IsSame( sideB->FirstVertex() ) &&
-             vt.IsSame( sideT->LastVertex() ))
+        const Prism_3D::TQuadList& quadList = thePrism.myWallQuads[iQ];
+        StdMeshers_FaceSidePtr      botSide = quadList.front()->side[ QUAD_BOTTOM_SIDE ];
+        StdMeshers_FaceSidePtr      topSide = quadList.back ()->side[ QUAD_TOP_SIDE ];
+        if ( botSide->NbEdges() == topSide->NbEdges() )
         {
-          TAssocTool::InsertAssociation( botSide->Edge( 0 ),
-                                         topSide->Edge( 0 ), shape2ShapeMap );
-          TAssocTool::InsertAssociation( vb, vt, shape2ShapeMap );
+          for ( int iE = 0; iE < botSide->NbEdges(); ++iE )
+          {
+            NSProjUtils::InsertAssociation( botSide->Edge( iE ),
+                                           topSide->Edge( iE ), shape2ShapeMap );
+            NSProjUtils::InsertAssociation( myHelper->IthVertex( 0, botSide->Edge( iE )),
+                                           myHelper->IthVertex( 0, topSide->Edge( iE )),
+                                           shape2ShapeMap );
+          }
         }
-        vb = myHelper->IthVertex( 1, botSide->Edge( botSide->NbEdges()-1 ));
-        vt = myHelper->IthVertex( 1, topSide->Edge( topSide->NbEdges()-1 ));
-        sideB = quadList.front()->side[ QUAD_RIGHT_SIDE ];
-        sideT = quadList.back ()->side[ QUAD_RIGHT_SIDE ];
-        if ( vb.IsSame( sideB->FirstVertex() ) &&
-             vt.IsSame( sideT->LastVertex() ))
+        else
         {
-          TAssocTool::InsertAssociation( botSide->Edge( botSide->NbEdges()-1 ),
-                                         topSide->Edge( topSide->NbEdges()-1 ),
-                                         shape2ShapeMap );
-          TAssocTool::InsertAssociation( vb, vt, shape2ShapeMap );
+          TopoDS_Vertex vb, vt;
+          StdMeshers_FaceSidePtr sideB, sideT;
+          vb = myHelper->IthVertex( 0, botSide->Edge( 0 ));
+          vt = myHelper->IthVertex( 0, topSide->Edge( 0 ));
+          sideB = quadList.front()->side[ QUAD_LEFT_SIDE ];
+          sideT = quadList.back ()->side[ QUAD_LEFT_SIDE ];
+          if ( vb.IsSame( sideB->FirstVertex() ) &&
+               vt.IsSame( sideT->LastVertex() ))
+          {
+            NSProjUtils::InsertAssociation( botSide->Edge( 0 ),
+                                           topSide->Edge( 0 ), shape2ShapeMap );
+            NSProjUtils::InsertAssociation( vb, vt, shape2ShapeMap );
+          }
+          vb = myHelper->IthVertex( 1, botSide->Edge( botSide->NbEdges()-1 ));
+          vt = myHelper->IthVertex( 1, topSide->Edge( topSide->NbEdges()-1 ));
+          sideB = quadList.front()->side[ QUAD_RIGHT_SIDE ];
+          sideT = quadList.back ()->side[ QUAD_RIGHT_SIDE ];
+          if ( vb.IsSame( sideB->FirstVertex() ) &&
+               vt.IsSame( sideT->LastVertex() ))
+          {
+            NSProjUtils::InsertAssociation( botSide->Edge( botSide->NbEdges()-1 ),
+                                           topSide->Edge( topSide->NbEdges()-1 ),
+                                           shape2ShapeMap );
+            NSProjUtils::InsertAssociation( vb, vt, shape2ShapeMap );
+          }
         }
       }
-    }
 
-  // Find matching nodes of top and bottom faces
-  TNodeNodeMap n2nMap;
-  if ( ! TAssocTool::FindMatchingNodesOnFaces( botFace, myBlock.Mesh(),
-                                               topFace, myBlock.Mesh(),
-                                               shape2ShapeMap, n2nMap ))
-  {
-    if ( sameTopo )
-      return toSM( error(TCom("Mesh on faces #") << botSM->GetId()
-                         <<" and #"<< topSM->GetId() << " seems different" ));
-    else
-      return toSM( error(TCom("Topology of faces #") << botSM->GetId()
-                         <<" and #"<< topSM->GetId() << " seems different" ));
+    // Find matching nodes of top and bottom faces
+    n2nMapPtr = & n2nMap;
+    if ( ! NSProjUtils::FindMatchingNodesOnFaces( thePrism.myBottom, myHelper->GetMesh(),
+                                                 thePrism.myTop,    myHelper->GetMesh(),
+                                                 shape2ShapeMap, n2nMap ))
+    {
+      if ( sameTopo )
+        return toSM( error(TCom("Mesh on faces #") << botSM->GetId()
+                           <<" and #"<< topSM->GetId() << " seems different" ));
+      else
+        return toSM( error(TCom("Topology of faces #") << botSM->GetId()
+                           <<" and #"<< topSM->GetId() << " seems different" ));
+    }
   }
 
   // Fill myBotToColumnMap
 
   int zSize = myBlock.VerticalSize();
-  //TNode prevTNode;
-  TNodeNodeMap::iterator bN_tN = n2nMap.begin();
-  for ( ; bN_tN != n2nMap.end(); ++bN_tN )
+  TNodeNodeMap::const_iterator bN_tN = n2nMapPtr->begin();
+  for ( ; bN_tN != n2nMapPtr->end(); ++bN_tN )
   {
     const SMDS_MeshNode* botNode = bN_tN->first;
     const SMDS_MeshNode* topNode = bN_tN->second;
@@ -1984,17 +1993,22 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf
 
 //================================================================================
 /*!
- * \brief Remove quadrangles from the top face and
- * create triangles there by projection from the bottom
+ * \brief Remove faces from the top face and re-create them by projection from the bottom
  * \retval bool - a success or not
  */
 //================================================================================
 
-bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf )
+bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf &             bottomToTopTrsf,
+                                              const Prism_3D::TPrismTopo& thePrism )
 {
-  SMESHDS_Mesh*  meshDS = myBlock.MeshDS();
-  SMESH_subMesh * botSM = myBlock.SubMesh( ID_BOT_FACE );
-  SMESH_subMesh * topSM = myBlock.SubMesh( ID_TOP_FACE );
+  if ( project2dMesh( thePrism.myBottom, thePrism.myTop ))
+  {
+    return true;
+  }
+
+  SMESHDS_Mesh*  meshDS = myHelper->GetMeshDS();
+  SMESH_subMesh * botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom );
+  SMESH_subMesh * topSM = myHelper->GetMesh()->GetSubMesh( thePrism.myTop );
 
   SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS();
   SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS();
@@ -2002,9 +2016,9 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf )
   if ( topSMDS && topSMDS->NbElements() > 0 )
     topSM->ComputeStateEngine( SMESH_subMesh::CLEAN );
 
-  const TopoDS_Face& botFace = TopoDS::Face( myBlock.Shape( ID_BOT_FACE )); // oriented within
-  const TopoDS_Face& topFace = TopoDS::Face( myBlock.Shape( ID_TOP_FACE )); //    the 3D SHAPE
-  int topFaceID = meshDS->ShapeToIndex( topFace );
+  const TopoDS_Face& botFace = thePrism.myBottom; // oriented within
+  const TopoDS_Face& topFace = thePrism.myTop;    //    the 3D SHAPE
+  int topFaceID = meshDS->ShapeToIndex( thePrism.myTop );
 
   SMESH_MesherHelper botHelper( *myHelper->GetMesh() );
   botHelper.SetSubShape( botFace );
@@ -2071,6 +2085,9 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf )
     column.resize( zSize );
     column.front() = botNode;
     column.back()  = topNode;
+
+    if ( _computeCanceled )
+      return toSM( error( SMESH_ComputeError::New(COMPERR_CANCELED)));
   }
 
   // Create top faces
@@ -2136,6 +2153,77 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf )
 }
 
 //=======================================================================
+//function : getSweepTolerance
+//purpose  : Compute tolerance to pass to StdMeshers_Sweeper
+//=======================================================================
+
+double StdMeshers_Prism_3D::getSweepTolerance( const Prism_3D::TPrismTopo& thePrism )
+{
+  SMESHDS_Mesh*    meshDS = myHelper->GetMeshDS();
+  SMESHDS_SubMesh * sm[2] = { meshDS->MeshElements( thePrism.myBottom ),
+                              meshDS->MeshElements( thePrism.myTop )    };
+  double minDist = 1e100;
+
+  vector< SMESH_TNodeXYZ > nodes;
+  for ( int iSM = 0; iSM < 2; ++iSM )
+  {
+    if ( !sm[ iSM ]) continue;
+
+    SMDS_ElemIteratorPtr fIt = sm[ iSM ]->GetElements();
+    while ( fIt->more() )
+    {
+      const SMDS_MeshElement* face = fIt->next();
+      const int            nbNodes = face->NbCornerNodes();
+      SMDS_ElemIteratorPtr     nIt = face->nodesIterator();
+
+      nodes.resize( nbNodes + 1 );
+      for ( int iN = 0; iN < nbNodes; ++iN )
+        nodes[ iN ] = nIt->next();
+      nodes.back() = nodes[0];
+      
+      // loop on links
+      double dist2;
+      for ( int iN = 0; iN < nbNodes; ++iN )
+      {
+        if ( nodes[ iN   ]._node->GetPosition()->GetDim() < 2 &&
+             nodes[ iN+1 ]._node->GetPosition()->GetDim() < 2 )
+        {
+          // it's a boundary link; measure distance of other
+          // nodes to this link
+          gp_XYZ linkDir = nodes[ iN ] - nodes[ iN+1 ];
+          double linkLen = linkDir.Modulus();
+          bool   isDegen = ( linkLen < numeric_limits<double>::min() );
+          if ( !isDegen ) linkDir /= linkLen;
+          for ( int iN2 = 0; iN2 < nbNodes; ++iN2 ) // loop on other nodes
+          {
+            if ( nodes[ iN2 ] == nodes[ iN ] ||
+                 nodes[ iN2 ] == nodes[ iN+1 ]) continue;
+            if ( isDegen )
+            {
+              dist2 = ( nodes[ iN ] - nodes[ iN2 ]).SquareModulus();
+            }
+            else
+            {
+              dist2 = linkDir.CrossSquareMagnitude( nodes[ iN ] - nodes[ iN2 ]);
+            }
+            if ( dist2 > numeric_limits<double>::min() )
+              minDist = Min ( minDist, dist2 );
+          }
+        }
+        // measure length link
+        else if ( nodes[ iN ]._node < nodes[ iN+1 ]._node ) // not to measure same link twice
+        {
+          dist2 = ( nodes[ iN ] - nodes[ iN+1 ]).SquareModulus();
+          if ( dist2 > numeric_limits<double>::min() )
+            minDist = Min ( minDist, dist2 );
+        }
+      }
+    }
+  }
+  return 0.1 * Sqrt ( minDist );
+}
+
+//=======================================================================
 //function : project2dMesh
 //purpose  : Project mesh faces from a source FACE of one prism (theSrcFace)
 //           to a source FACE of another prism (theTgtFace)
@@ -2261,6 +2349,12 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
         if ( E.IsSame( Edge( i ))) return i;
       return -1;
     }
+    bool IsSideFace( const TopoDS_Shape& face ) const
+    {
+      if ( _faces->Contains( face )) // avoid returning true for a prism top FACE
+        return ( !_face.IsNull() || !( face.IsSame( _faces->FindKey( _faces->Extent() ))));
+      return false;
+    }
   };
   //--------------------------------------------------------------------------------
   /*!
@@ -2392,12 +2486,11 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
 
       TEdgeWithNeighborsVec& botEdges = faceEdgesVec[ iF ];
       if ( botEdges.empty() )
-      {
         if ( !getEdges( botF, botEdges, /*noHoles=*/false ))
           break;
-        if ( allFaces.Extent()-1 <= (int) botEdges.size() )
-          continue; // all faces are adjacent to botF - no top FACE
-      }
+      if ( allFaces.Extent()-1 <= (int) botEdges.size() )
+        continue; // all faces are adjacent to botF - no top FACE
+
       // init data of side FACEs
       vector< PrismSide > sides( botEdges.size() );
       for ( int iS = 0; iS < botEdges.size(); ++iS )
@@ -2411,8 +2504,8 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
       }
 
       bool isOK = true; // ok for a current botF
-      bool isAdvanced = true;
-      int nbFoundSideFaces = 0;
+      bool isAdvanced = true; // is new data found in a current loop
+      int  nbFoundSideFaces = 0;
       for ( int iLoop = 0; isOK && isAdvanced; ++iLoop )
       {
         isAdvanced = false;
@@ -2420,7 +2513,8 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
         {
           PrismSide& side = sides[ iS ];
           if ( side._face.IsNull() )
-            continue;
+            continue; // probably the prism top face is the last of side._faces
+
           if ( side._topEdge.IsNull() )
           {
             // find vertical EDGEs --- EGDEs shared with neighbor side FACEs
@@ -2434,7 +2528,7 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
                 if ( side._isCheckedEdge[ iE ] ) continue;
                 const TopoDS_Edge&      vertE = side.Edge( iE );
                 const TopoDS_Shape& neighborF = getAnotherFace( side._face, vertE, facesOfEdge );
-                bool isEdgeShared = adjSide->_faces->Contains( neighborF );
+                bool isEdgeShared = adjSide->IsSideFace( neighborF );
                 if ( isEdgeShared )
                 {
                   isAdvanced = true;
@@ -2480,13 +2574,13 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
             {
               if ( side._leftSide->_faces->Contains( f ))
               {
-                stop = true;
+                stop = true; // probably f is the prism top face
                 side._leftSide->_face.Nullify();
                 side._leftSide->_topEdge.Nullify();
               }
               if ( side._rightSide->_faces->Contains( f ))
               {
-                stop = true;
+                stop = true; // probably f is the prism top face
                 side._rightSide->_face.Nullify();
                 side._rightSide->_topEdge.Nullify();
               }
@@ -2497,8 +2591,8 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
               side._topEdge.Nullify();
               continue;
             }
-            side._face = TopoDS::Face( f );
-            int faceID = allFaces.FindIndex( side._face );
+            side._face  = TopoDS::Face( f );
+            int faceID  = allFaces.FindIndex( side._face );
             side._edges = & faceEdgesVec[ faceID ];
             if ( side._edges->empty() )
               if ( !getEdges( side._face, * side._edges, /*noHoles=*/true ))
@@ -2805,8 +2899,9 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism,
   thePrism.myShape3D = shape3D;
   if ( thePrism.myBottom.IsNull() )
     thePrism.myBottom  = TopoDS::Face( botSM->GetSubShape() );
-  thePrism.myBottom.Orientation( myHelper->GetSubShapeOri( shape3D,
-                                                           thePrism.myBottom ));
+  thePrism.myBottom.Orientation( myHelper->GetSubShapeOri( shape3D, thePrism.myBottom ));
+  thePrism.myTop.   Orientation( myHelper->GetSubShapeOri( shape3D, thePrism.myTop ));
+
   // Get ordered bottom edges
   TopoDS_Face reverseBottom = // to have order of top EDGEs as in the top FACE
     TopoDS::Face( thePrism.myBottom.Reversed() );
@@ -4221,3 +4316,355 @@ gp_Pnt2d StdMeshers_PrismAsBlock::TPCurveOnHorFaceAdaptor::Value(const Standard_
   double r = ( U - i1->first ) / ( i2->first - i1->first );
   return i1->second * ( 1 - r ) + i2->second * r;
 }
+
+//================================================================================
+/*!
+ * \brief Projects internal nodes using transformation found by boundary nodes
+ */
+//================================================================================
+
+bool StdMeshers_Sweeper::projectIntPoints(const vector< gp_XYZ >&    fromBndPoints,
+                                          const vector< gp_XYZ >&    toBndPoints,
+                                          const vector< gp_XYZ >&    fromIntPoints,
+                                          vector< gp_XYZ >&          toIntPoints,
+                                          NSProjUtils::TrsfFinder3D& trsf,
+                                          vector< gp_XYZ > *         bndError)
+{
+  // find transformation
+  if ( trsf.IsIdentity() && !trsf.Solve( fromBndPoints, toBndPoints ))
+    return false;
+
+  // compute internal points using the found trsf
+  for ( size_t iP = 0; iP < fromIntPoints.size(); ++iP )
+  {
+    toIntPoints[ iP ] = trsf.Transform( fromIntPoints[ iP ]);
+  }
+
+  // compute boundary error
+  if ( bndError )
+  {
+    bndError->resize( fromBndPoints.size() );
+    gp_XYZ fromTrsf;
+    for ( size_t iP = 0; iP < fromBndPoints.size(); ++iP )
+    {
+      fromTrsf = trsf.Transform( fromBndPoints[ iP ] );
+      (*bndError)[ iP ]  = toBndPoints[ iP ] - fromTrsf;
+    }
+  }
+  return true;
+}
+
+//================================================================================
+/*!
+ * \brief Add boundary error to ineternal points
+ */
+//================================================================================
+
+void StdMeshers_Sweeper::applyBoundaryError(const vector< gp_XYZ >& bndPoints,
+                                            const vector< gp_XYZ >& bndError1,
+                                            const vector< gp_XYZ >& bndError2,
+                                            const double            r,
+                                            vector< gp_XYZ >&       intPoints,
+                                            vector< double >&       int2BndDist)
+{
+  // fix each internal point
+  const double eps = 1e-100;
+  for ( size_t iP = 0; iP < intPoints.size(); ++iP )
+  {
+    gp_XYZ & intPnt = intPoints[ iP ];
+
+    // compute distance from intPnt to each boundary node
+    double int2BndDistSum = 0;
+    for ( size_t iBnd = 0; iBnd < bndPoints.size(); ++iBnd )
+    {
+      int2BndDist[ iBnd ] = 1 / (( intPnt - bndPoints[ iBnd ]).SquareModulus() + eps );
+      int2BndDistSum += int2BndDist[ iBnd ];
+    }
+
+    // apply bndError
+    for ( size_t iBnd = 0; iBnd < bndPoints.size(); ++iBnd )
+    {
+      intPnt += bndError1[ iBnd ] * ( 1 - r ) * int2BndDist[ iBnd ] / int2BndDistSum;
+      intPnt += bndError2[ iBnd ] * r         * int2BndDist[ iBnd ] / int2BndDistSum;
+    }
+  }
+}
+
+//================================================================================
+/*!
+ * \brief Creates internal nodes of the prism
+ */
+//================================================================================
+
+bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper,
+                                       const double        tol)
+{
+  const size_t zSize = myBndColumns[0]->size();
+  const size_t zSrc = 0, zTgt = zSize-1;
+  if ( zSize < 3 ) return true;
+
+  vector< vector< gp_XYZ > > intPntsOfLayer( zSize ); // node coodinates to compute
+  // set coordinates of src and tgt nodes
+  for ( size_t z = 0; z < intPntsOfLayer.size(); ++z )
+    intPntsOfLayer[ z ].resize( myIntColumns.size() );
+  for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
+  {
+    intPntsOfLayer[ zSrc ][ iP ] = intPoint( iP, zSrc );
+    intPntsOfLayer[ zTgt ][ iP ] = intPoint( iP, zTgt );
+  }
+
+  // compute coordinates of internal nodes by projecting (transfroming) src and tgt
+  // nodes towards the central layer
+
+  vector< NSProjUtils::TrsfFinder3D > trsfOfLayer( zSize );
+  vector< vector< gp_XYZ > >          bndError( zSize );
+
+  // boundary points used to compute an affine transformation from a layer to a next one
+  vector< gp_XYZ > fromSrcBndPnts( myBndColumns.size() ), fromTgtBndPnts( myBndColumns.size() );
+  vector< gp_XYZ > toSrcBndPnts  ( myBndColumns.size() ), toTgtBndPnts  ( myBndColumns.size() );
+  for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
+  {
+    fromSrcBndPnts[ iP ] = bndPoint( iP, zSrc );
+    fromTgtBndPnts[ iP ] = bndPoint( iP, zTgt );
+  }
+
+  size_t zS = zSrc + 1;
+  size_t zT = zTgt - 1;
+  for ( ; zS < zT; ++zS, --zT ) // vertical loop on layers
+  {
+    for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
+    {
+      toSrcBndPnts[ iP ] = bndPoint( iP, zS );
+      toTgtBndPnts[ iP ] = bndPoint( iP, zT );
+    }
+    if (! projectIntPoints( fromSrcBndPnts, toSrcBndPnts,
+                            intPntsOfLayer[ zS-1 ], intPntsOfLayer[ zS ],
+                            trsfOfLayer   [ zS-1 ], & bndError[ zS-1 ]))
+      return false;
+    if (! projectIntPoints( fromTgtBndPnts, toTgtBndPnts,
+                            intPntsOfLayer[ zT+1 ], intPntsOfLayer[ zT ],
+                            trsfOfLayer   [ zT+1 ], & bndError[ zT+1 ]))
+      return false;
+
+    // if ( zT == zTgt - 1 )
+    // {
+    //   for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
+    //   {
+    //     gp_XYZ fromTrsf = trsfOfLayer   [ zT+1].Transform( fromTgtBndPnts[ iP ] );
+    //     cout << "mesh.AddNode( "
+    //          << fromTrsf.X() << ", "
+    //          << fromTrsf.Y() << ", "
+    //          << fromTrsf.Z() << ") " << endl;
+    //   }
+    //   for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
+    //     cout << "mesh.AddNode( "
+    //          << intPntsOfLayer[ zT ][ iP ].X() << ", "
+    //          << intPntsOfLayer[ zT ][ iP ].Y() << ", "
+    //          << intPntsOfLayer[ zT ][ iP ].Z() << ") " << endl;
+    // }
+
+    fromTgtBndPnts.swap( toTgtBndPnts );
+    fromSrcBndPnts.swap( toSrcBndPnts );
+  }
+
+  // Compute two projections of internal points to the central layer
+  // in order to evaluate an error of internal points
+
+  bool centerIntErrorIsSmall;
+  vector< gp_XYZ > centerSrcIntPnts( myIntColumns.size() );
+  vector< gp_XYZ > centerTgtIntPnts( myIntColumns.size() );
+
+  for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
+  {
+    toSrcBndPnts[ iP ] = bndPoint( iP, zS );
+    toTgtBndPnts[ iP ] = bndPoint( iP, zT );
+  }
+  if (! projectIntPoints( fromSrcBndPnts, toSrcBndPnts,
+                          intPntsOfLayer[ zS-1 ], centerSrcIntPnts,
+                          trsfOfLayer   [ zS-1 ], & bndError[ zS-1 ]))
+    return false;
+  if (! projectIntPoints( fromTgtBndPnts, toTgtBndPnts,
+                          intPntsOfLayer[ zT+1 ], centerTgtIntPnts,
+                          trsfOfLayer   [ zT+1 ], & bndError[ zT+1 ]))
+    return false;
+
+  // evaluate an error of internal points on the central layer
+  centerIntErrorIsSmall = true;
+  if ( zS == zT ) // odd zSize
+  {
+    for ( size_t iP = 0; ( iP < myIntColumns.size() && centerIntErrorIsSmall ); ++iP )
+      centerIntErrorIsSmall =
+        (centerSrcIntPnts[ iP ] - centerTgtIntPnts[ iP ]).SquareModulus() < tol*tol;
+  }
+  else // even zSize
+  {
+    for ( size_t iP = 0; ( iP < myIntColumns.size() && centerIntErrorIsSmall ); ++iP )
+      centerIntErrorIsSmall =
+        (intPntsOfLayer[ zS-1 ][ iP ] - centerTgtIntPnts[ iP ]).SquareModulus() < tol*tol;
+  }
+
+  // Evaluate an error of boundary points
+
+  bool bndErrorIsSmall = true;
+  for ( size_t iP = 0; ( iP < myBndColumns.size() && bndErrorIsSmall ); ++iP )
+  {
+    double sumError = 0;
+    for ( size_t z = 1; z < zS; ++z ) // loop on layers
+      sumError += ( bndError[ z-1     ][ iP ].Modulus() +
+                    bndError[ zSize-z ][ iP ].Modulus() );
+
+    bndErrorIsSmall = ( sumError < tol );
+  }
+
+  // compute final points on the central layer
+  std::vector< double > int2BndDist( myBndColumns.size() ); // work array of applyBoundaryError()
+  double r = zS / ( zSize - 1.);
+  if ( zS == zT )
+  {
+    for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
+    {
+      intPntsOfLayer[ zS ][ iP ] =
+        ( 1 - r ) * centerSrcIntPnts[ iP ] + r * centerTgtIntPnts[ iP ];
+    }
+    if ( !bndErrorIsSmall )
+    {
+      applyBoundaryError( toSrcBndPnts, bndError[ zS-1 ], bndError[ zS+1 ], r,
+                          intPntsOfLayer[ zS ], int2BndDist );
+    }
+  }
+  else
+  {
+    for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
+    {
+      intPntsOfLayer[ zS ][ iP ] =
+        r * intPntsOfLayer[ zS ][ iP ] + ( 1 - r ) * centerSrcIntPnts[ iP ];
+      intPntsOfLayer[ zT ][ iP ] =
+        r * intPntsOfLayer[ zT ][ iP ] + ( 1 - r ) * centerTgtIntPnts[ iP ];
+    }
+    if ( !bndErrorIsSmall )
+    {
+      applyBoundaryError( toSrcBndPnts, bndError[ zS-1 ], bndError[ zS+1 ], r,
+                          intPntsOfLayer[ zS ], int2BndDist );
+      applyBoundaryError( toTgtBndPnts, bndError[ zT+1 ], bndError[ zT-1 ], r,
+                          intPntsOfLayer[ zT ], int2BndDist );
+    }
+  }
+
+  //centerIntErrorIsSmall = true;
+  //bndErrorIsSmall = true;
+  if ( !centerIntErrorIsSmall )
+  {
+    // Compensate the central error; continue adding projection
+    // by going from central layer to the source and target ones
+
+    vector< gp_XYZ >& fromSrcIntPnts = centerSrcIntPnts;
+    vector< gp_XYZ >& fromTgtIntPnts = centerTgtIntPnts;
+    vector< gp_XYZ >  toSrcIntPnts( myIntColumns.size() );
+    vector< gp_XYZ >  toTgtIntPnts( myIntColumns.size() );
+    vector< gp_XYZ >  srcBndError( myBndColumns.size() );
+    vector< gp_XYZ >  tgtBndError( myBndColumns.size() );
+
+    fromTgtBndPnts.swap( toTgtBndPnts );
+    fromSrcBndPnts.swap( toSrcBndPnts );
+
+    for ( ++zS, --zT; zS < zTgt; ++zS, --zT ) // vertical loop on layers
+    {
+      // invert transformation
+      if ( !trsfOfLayer[ zS+1 ].Invert() )
+        trsfOfLayer[ zS+1 ] = NSProjUtils::TrsfFinder3D(); // to recompute
+      if ( !trsfOfLayer[ zT-1 ].Invert() )
+        trsfOfLayer[ zT-1 ] = NSProjUtils::TrsfFinder3D();
+
+      // project internal nodes and compute bnd error
+      for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
+      {
+        toSrcBndPnts[ iP ] = bndPoint( iP, zS );
+        toTgtBndPnts[ iP ] = bndPoint( iP, zT );
+      }
+      projectIntPoints( fromSrcBndPnts, toSrcBndPnts,
+                        fromSrcIntPnts, toSrcIntPnts,
+                        trsfOfLayer[ zS+1 ], & srcBndError );
+      projectIntPoints( fromTgtBndPnts, toTgtBndPnts,
+                        fromTgtIntPnts, toTgtIntPnts,
+                        trsfOfLayer[ zT-1 ], & tgtBndError );
+
+      // if ( zS == zTgt - 1 )
+      // {
+      //   cout << "mesh2 = smesh.Mesh()" << endl;
+      //   for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
+      //   {
+      //     gp_XYZ fromTrsf = trsfOfLayer   [ zS+1].Transform( fromSrcBndPnts[ iP ] );
+      //     cout << "mesh2.AddNode( "
+      //          << fromTrsf.X() << ", "
+      //          << fromTrsf.Y() << ", "
+      //          << fromTrsf.Z() << ") " << endl;
+      //   }
+      //   for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
+      //     cout << "mesh2.AddNode( "
+      //          << toSrcIntPnts[ iP ].X() << ", "
+      //          << toSrcIntPnts[ iP ].Y() << ", "
+      //          << toSrcIntPnts[ iP ].Z() << ") " << endl;
+      // }
+
+      // sum up 2 projections
+      r = zS / ( zSize - 1.);
+      vector< gp_XYZ >& zSIntPnts = intPntsOfLayer[ zS ];
+      vector< gp_XYZ >& zTIntPnts = intPntsOfLayer[ zT ];
+      for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
+      {
+        zSIntPnts[ iP ] = r * zSIntPnts[ iP ]  +  ( 1 - r ) * toSrcIntPnts[ iP ];
+        zTIntPnts[ iP ] = r * zTIntPnts[ iP ]  +  ( 1 - r ) * toTgtIntPnts[ iP ];
+      }
+
+      // compensate bnd error
+      if ( !bndErrorIsSmall )
+      {
+        applyBoundaryError( toSrcBndPnts, srcBndError, bndError[ zS+1 ], r,
+                            intPntsOfLayer[ zS ], int2BndDist );
+        applyBoundaryError( toTgtBndPnts, tgtBndError, bndError[ zT-1 ], r,
+                            intPntsOfLayer[ zT ], int2BndDist );
+      }
+
+      fromSrcBndPnts.swap( toSrcBndPnts );
+      fromSrcIntPnts.swap( toSrcIntPnts );
+      fromTgtBndPnts.swap( toTgtBndPnts );
+      fromTgtIntPnts.swap( toTgtIntPnts );
+    }
+  }  // if ( !centerIntErrorIsSmall )
+
+  else if ( !bndErrorIsSmall )
+  {
+    zS = zSrc + 1;
+    zT = zTgt - 1;
+    for ( ; zS < zT; ++zS, --zT ) // vertical loop on layers
+    {
+      for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
+      {
+        toSrcBndPnts[ iP ] = bndPoint( iP, zS );
+        toTgtBndPnts[ iP ] = bndPoint( iP, zT );
+      }
+      // compensate bnd error
+      applyBoundaryError( toSrcBndPnts, bndError[ zS-1 ], bndError[ zS-1 ], 0.5,
+                          intPntsOfLayer[ zS ], int2BndDist );
+      applyBoundaryError( toTgtBndPnts, bndError[ zT+1 ], bndError[ zT+1 ], 0.5,
+                          intPntsOfLayer[ zT ], int2BndDist );
+    }
+  }
+
+  // cout << "centerIntErrorIsSmall = " << centerIntErrorIsSmall<< endl;
+  // cout << "bndErrorIsSmall = " << bndErrorIsSmall<< endl;
+
+  // Create nodes
+  for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
+  {
+    vector< const SMDS_MeshNode* > & nodeCol = *myIntColumns[ iP ];
+    for ( size_t z = zSrc + 1; z < zTgt; ++z ) // vertical loop on layers
+    {
+      const gp_XYZ & xyz = intPntsOfLayer[ z ][ iP ];
+      if ( !( nodeCol[ z ] = helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() )))
+        return false;
+    }
+  }
+
+  return true;
+}
index f2929c3..6c7c52f 100644 (file)
@@ -56,6 +56,10 @@ namespace Prism_3D
   struct TNode;
   struct TPrismTopo;
 }
+namespace StdMeshers_ProjectionUtils
+{
+  class TrsfFinder3D;
+}
 class SMESHDS_SubMesh;
 class TopoDS_Edge;
 
@@ -398,7 +402,42 @@ private:
 
 }; // class StdMeshers_PrismAsBlock
 
-// =============================================
+// ===============================================
+/*!
+ * \brief Tool building internal nodes in a prism
+ */
+struct StdMeshers_Sweeper
+{
+  std::vector< TNodeColumn* > myBndColumns; // boundary nodes
+  std::vector< TNodeColumn* > myIntColumns; // internal nodes
+
+  bool ComputeNodes( SMESH_MesherHelper& helper,
+                     const double        tol );
+
+private:
+
+  gp_XYZ bndPoint( int iP, int z ) const
+  { return SMESH_TNodeXYZ( (*myBndColumns[ iP ])[ z ]); }
+
+  gp_XYZ intPoint( int iP, int z ) const
+  { return SMESH_TNodeXYZ( (*myIntColumns[ iP ])[ z ]); }
+
+  static bool projectIntPoints(const std::vector< gp_XYZ >& fromBndPoints,
+                               const std::vector< gp_XYZ >& toBndPoints,
+                               const std::vector< gp_XYZ >& fromIntPoints,
+                               std::vector< gp_XYZ >&       toIntPoints,
+                               StdMeshers_ProjectionUtils::TrsfFinder3D& trsf,
+                               std::vector< gp_XYZ > *      bndError);
+
+  static void applyBoundaryError(const std::vector< gp_XYZ >& bndPoints,
+                                 const std::vector< gp_XYZ >& bndError1,
+                                 const std::vector< gp_XYZ >& bndError2,
+                                 const double                 r,
+                                 std::vector< gp_XYZ >&       toIntPoints,
+                                 std::vector< double >&       int2BndDist);
+};
+
+// ===============================================
 /*!
  * \brief Algo building prisms on a prism shape
  */
@@ -480,7 +519,13 @@ public:
    *        create triangles there by projection from the bottom
     * \retval bool - a success or not
    */
-  bool projectBottomToTop( const gp_Trsf & bottomToTopTrsf );
+  bool projectBottomToTop( const gp_Trsf & bottomToTopTrsf,
+                           const Prism_3D::TPrismTopo& thePrism );
+
+  /*!
+   * \brief Compute tolerance to pass to StdMeshers_Sweeper
+   */
+  double getSweepTolerance( const Prism_3D::TPrismTopo& thePrism );
 
   /*!
    * \brief Project mesh faces from a source FACE of one prism to
index 92474f2..df96a9c 100644 (file)
@@ -66,6 +66,7 @@
 #include <TopoDS_Shape.hxx>
 #include <gp_Pnt.hxx>
 #include <gp_Vec.hxx>
+#include <math_Gauss.hxx>
 
 #include <numeric>
 #include <limits>
@@ -98,7 +99,7 @@ namespace HERE = StdMeshers_ProjectionUtils;
 
 namespace {
 
-  static SMESHDS_Mesh* theMeshDS[2] = { 0, 0 }; // used to debug only
+  static SMESHDS_Mesh* theMeshDS[2] = { 0, 0 }; // used for debug only
   long shapeIndex(const TopoDS_Shape& S)
   {
     if ( theMeshDS[0] && theMeshDS[1] )
@@ -2402,3 +2403,229 @@ void StdMeshers_ProjectionUtils::SetEventListener(SMESH_subMesh* subMesh,
     }
   }
 }
+
+namespace StdMeshers_ProjectionUtils
+{
+
+  //================================================================================
+  /*!
+   * \brief Computes transformation beween two sets of 2D points using
+   *        a least square approximation
+   *
+   * See "Surface Mesh Projection For Hexahedral Mesh Generation By Sweeping"
+   * by X.Roca, J.Sarrate, A.Huerta. (2.2)
+   */
+  //================================================================================
+
+  bool TrsfFinder2D::Solve( const vector< gp_XY >& srcPnts,
+                            const vector< gp_XY >& tgtPnts )
+  {
+    // find gravity centers
+    gp_XY srcGC( 0,0 ), tgtGC( 0,0 );
+    for ( size_t i = 0; i < srcPnts.size(); ++i )
+    {
+      srcGC += srcPnts[i];
+      tgtGC += tgtPnts[i];
+    }
+    srcGC /= srcPnts.size();
+    tgtGC /= tgtPnts.size();
+
+    // find trsf
+
+    math_Matrix mat (1,4,1,4, 0.);
+    math_Vector vec (1,4, 0.);
+
+    // cout << "m1 = smesh.Mesh('src')" << endl
+    //      << "m2 = smesh.Mesh('tgt')" << endl;
+    double xx = 0, xy = 0, yy = 0;
+    for ( size_t i = 0; i < srcPnts.size(); ++i )
+    {
+      gp_XY srcUV = srcPnts[i] - srcGC;
+      gp_XY tgtUV = tgtPnts[i] - tgtGC;
+      xx += srcUV.X() * srcUV.X();
+      yy += srcUV.Y() * srcUV.Y();
+      xy += srcUV.X() * srcUV.Y();
+      vec( 1 ) += srcUV.X() * tgtUV.X();
+      vec( 2 ) += srcUV.Y() * tgtUV.X();
+      vec( 3 ) += srcUV.X() * tgtUV.Y();
+      vec( 4 ) += srcUV.Y() * tgtUV.Y();
+      // cout << "m1.AddNode( " << srcUV.X() << ", " << srcUV.Y() << ", 0 )" << endl
+      //      << "m2.AddNode( " << tgtUV.X() << ", " << tgtUV.Y() << ", 0 )" << endl;
+    }
+    mat( 1,1 ) = mat( 3,3 ) = xx;
+    mat( 2,2 ) = mat( 4,4 ) = yy;
+    mat( 1,2 ) = mat( 2,1 ) = mat( 3,4 ) = mat( 4,3 ) = xy;
+
+    math_Gauss solver( mat );
+    if ( !solver.IsDone() )
+      return false;
+    solver.Solve( vec );
+    if ( vec.Norm2() < gp::Resolution() )
+      return false;
+    // cout << vec( 1 ) << "\t " << vec( 2 ) << endl
+    //      << vec( 3 ) << "\t " << vec( 4 ) << endl;
+
+    _trsf.SetTranslation( tgtGC );
+    _srcOrig = srcGC;
+
+    gp_Mat2d& M = const_cast< gp_Mat2d& >( _trsf.HVectorialPart());
+    M( 1,1 ) = vec( 1 );
+    M( 2,1 ) = vec( 2 );
+    M( 1,2 ) = vec( 3 );
+    M( 2,2 ) = vec( 4 );
+
+    return true;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Transforms a 2D points using a found transformation
+   */
+  //================================================================================
+
+  gp_XY TrsfFinder2D::Transform( const gp_Pnt2d& srcUV ) const
+  {
+    gp_XY uv = srcUV.XY() - _srcOrig ;
+    _trsf.Transforms( uv );
+    return uv;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Computes transformation beween two sets of 3D points using
+   *        a least square approximation
+   *
+   * See "Surface Mesh Projection For Hexahedral Mesh Generation By Sweeping"
+   * by X.Roca, J.Sarrate, A.Huerta. (2.4)
+   */
+  //================================================================================
+
+  bool TrsfFinder3D::Solve( const vector< gp_XYZ > & srcPnts,
+                            const vector< gp_XYZ > & tgtPnts )
+  {
+    // find gravity center
+    gp_XYZ srcGC( 0,0,0 ), tgtGC( 0,0,0 );
+    for ( size_t i = 0; i < srcPnts.size(); ++i )
+    {
+      srcGC += srcPnts[i];
+      tgtGC += tgtPnts[i];
+    }
+    srcGC /= srcPnts.size();
+    tgtGC /= tgtPnts.size();
+
+    gp_XYZ srcOrig = 2 * srcGC - tgtGC;
+    gp_XYZ tgtOrig = srcGC;
+
+    // find trsf
+
+    math_Matrix mat (1,9,1,9, 0.);
+    math_Vector vec (1,9, 0.);
+
+    double xx = 0, yy = 0, zz = 0;
+    double xy = 0, xz = 0, yz = 0;
+    for ( size_t i = 0; i < srcPnts.size(); ++i )
+    {
+      gp_XYZ src = srcPnts[i] - srcOrig;
+      gp_XYZ tgt = tgtPnts[i] - tgtOrig;
+      xx += src.X() * src.X();
+      yy += src.Y() * src.Y();
+      zz += src.Z() * src.Z();
+      xy += src.X() * src.Y();
+      xz += src.X() * src.Z();
+      yz += src.Y() * src.Z();
+      vec( 1 ) += src.X() * tgt.X();
+      vec( 2 ) += src.Y() * tgt.X();
+      vec( 3 ) += src.Z() * tgt.X();
+      vec( 4 ) += src.X() * tgt.Y();
+      vec( 5 ) += src.Y() * tgt.Y();
+      vec( 6 ) += src.Z() * tgt.Y();
+      vec( 7 ) += src.X() * tgt.Z();
+      vec( 8 ) += src.Y() * tgt.Z();
+      vec( 9 ) += src.Z() * tgt.Z();
+    }
+    mat( 1,1 ) = mat( 4,4 ) = mat( 7,7 ) = xx;
+    mat( 2,2 ) = mat( 5,5 ) = mat( 8,8 ) = yy;
+    mat( 3,3 ) = mat( 6,6 ) = mat( 9,9 ) = zz;
+    mat( 1,2 ) = mat( 2,1 ) = mat( 4,5 ) = mat( 5,4 ) = mat( 7,8 ) = mat( 8,7 ) = xy;
+    mat( 1,3 ) = mat( 3,1 ) = mat( 4,6 ) = mat( 6,4 ) = mat( 7,9 ) = mat( 9,7 ) = xz;
+    mat( 2,3 ) = mat( 3,2 ) = mat( 5,6 ) = mat( 6,5 ) = mat( 8,9 ) = mat( 9,8 ) = yz;
+
+    math_Gauss solver( mat );
+    if ( !solver.IsDone() )
+      return false;
+    solver.Solve( vec );
+    if ( vec.Norm2() < gp::Resolution() )
+      return false;
+    // cout << endl
+    //      << vec( 1 ) << "\t " << vec( 2 ) << "\t " << vec( 3 ) << endl
+    //      << vec( 4 ) << "\t " << vec( 5 ) << "\t " << vec( 6 ) << endl
+    //      << vec( 7 ) << "\t " << vec( 8 ) << "\t " << vec( 9 ) << endl;
+
+    _srcOrig = srcOrig;
+    _trsf.SetTranslation( tgtOrig );
+
+    gp_Mat& M = const_cast< gp_Mat& >( _trsf.HVectorialPart() );
+    M.SetRows( gp_XYZ( vec( 1 ), vec( 2 ), vec( 3 )),
+               gp_XYZ( vec( 4 ), vec( 5 ), vec( 6 )),
+               gp_XYZ( vec( 7 ), vec( 8 ), vec( 9 )));
+    return true;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Transforms a 3D point using a found transformation
+   */
+  //================================================================================
+
+  gp_XYZ TrsfFinder3D::Transform( const gp_Pnt& srcP ) const
+  {
+    gp_XYZ p = srcP.XYZ() - _srcOrig;
+    _trsf.Transforms( p );
+    return p;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Transforms a 3D vector using a found transformation
+   */
+  //================================================================================
+
+  gp_XYZ TrsfFinder3D::TransformVec( const gp_Vec& v ) const
+  {
+    return v.XYZ().Multiplied( _trsf.HVectorialPart() );
+  }
+  //================================================================================
+  /*!
+   * \brief Inversion
+   */
+  //================================================================================
+
+  bool TrsfFinder3D::Invert()
+  {
+    if (( _trsf.Form() == gp_Translation ) &&
+        ( _srcOrig.X() != 0 || _srcOrig.Y() != 0 || _srcOrig.Z() != 0 ))
+    {
+      // seems to be defined via Solve()
+      gp_XYZ newSrcOrig = _trsf.TranslationPart();
+      gp_Mat& M = const_cast< gp_Mat& >( _trsf.HVectorialPart() );
+      const double D = M.Determinant();
+      if ( D < 1e-3 * ( newSrcOrig - _srcOrig ).Modulus() )
+      {
+#ifdef _DEBUG_
+        cerr << "TrsfFinder3D::Invert()"
+             << "D " << M.Determinant() << " IsSingular " << M.IsSingular() << endl;
+#endif
+        return false;
+      }
+      gp_Mat Minv = M.Inverted();
+      _trsf.SetTranslation( _srcOrig );
+      _srcOrig = newSrcOrig;
+      M = Minv;
+    }
+    else
+    {
+      _trsf.Invert();
+    }
+    return true;
+  }
+}
index 31afb07..e7ea442 100644 (file)
 
 #include "SMESH_StdMeshers.hxx"
 
+#include "SMDS_MeshElement.hxx"
+
 #include <TopTools_DataMapOfShapeShape.hxx>
 #include <TopoDS_Edge.hxx>
-#include <TopoDS_Vertex.hxx>
 #include <TopoDS_Face.hxx>
+#include <TopoDS_Vertex.hxx>
+#include <gp_Trsf.hxx>
+#include <gp_Trsf2d.hxx>
 
 #include <list>
 #include <map>
@@ -77,7 +81,54 @@ namespace StdMeshers_ProjectionUtils
 {
   typedef StdMeshers_ShapeShapeBiDirectionMap                  TShapeShapeMap;
   typedef TopTools_IndexedDataMapOfShapeListOfShape            TAncestorMap;
-  typedef std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> TNodeNodeMap;
+  typedef std::map<const SMDS_MeshNode*, const SMDS_MeshNode*,
+                   TIDCompare>                                 TNodeNodeMap;
+
+
+  /*!
+   * \brief Finds transformation beween two sets of 2D points using
+   *        a least square approximation
+   */
+  class TrsfFinder2D
+  {
+    gp_Trsf2d _trsf;
+    gp_XY     _srcOrig;
+  public:
+    TrsfFinder2D(): _srcOrig(0,0) {}
+
+    void Set( const gp_Trsf2d& t ) { _trsf = t; } // it's an alternative to Solve()
+
+    bool Solve( const std::vector< gp_XY >& srcPnts,
+                const std::vector< gp_XY >& tgtPnts );
+
+    gp_XY Transform( const gp_Pnt2d& srcUV ) const;
+
+    bool IsIdentity() const { return ( _trsf.Form() == gp_Identity ); }
+  };
+  /*!
+   * \brief Finds transformation beween two sets of 3D points using
+   *        a least square approximation
+   */
+  class TrsfFinder3D
+  {
+    gp_Trsf _trsf;
+    gp_XYZ  _srcOrig;
+  public:
+    TrsfFinder3D(): _srcOrig(0,0,0) {}
+
+    void Set( const gp_Trsf& t ) { _trsf = t; } // it's an alternative to Solve()
+
+    bool Solve( const std::vector< gp_XYZ > & srcPnts,
+                const std::vector< gp_XYZ > & tgtPnts );
+
+    gp_XYZ Transform( const gp_Pnt& srcP ) const;
+
+    gp_XYZ TransformVec( const gp_Vec& v ) const;
+
+    bool IsIdentity() const { return ( _trsf.Form() == gp_Identity ); }
+
+    bool Invert();
+  };
 
   /*!
    * \brief Looks for association of all sub-shapes of two shapes
index b0118c7..3883290 100644 (file)
 
 #include "utilities.h"
 
+#include <BRepAdaptor_Surface.hxx>
 #include <BRep_Tool.hxx>
 #include <Bnd_B2d.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 #include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx>
@@ -374,123 +376,195 @@ namespace {
 
   //================================================================================
   /*!
+   * \brief Check if two consecutive EDGEs are connected in 2D
+   *  \param [in] E1 - a well oriented non-seam EDGE
+   *  \param [in] E2 - a possibly well oriented seam EDGE
+   *  \param [in] F - a FACE
+   *  \return bool - result
+   */
+  //================================================================================
+
+  bool are2dConnected( const TopoDS_Edge & E1,
+                       const TopoDS_Edge & E2,
+                       const TopoDS_Face & F )
+  {
+    double f,l;
+    Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( E1, F, f, l );
+    gp_Pnt2d uvLast1 = c1->Value( E1.Orientation() == TopAbs_REVERSED ? f : l );
+
+    Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( E2, F, f, l );
+    gp_Pnt2d uvFirst2 = c2->Value( f );
+    gp_Pnt2d uvLast2  = c2->Value( l );
+    double tol2 = 1e-5 * uvLast2.SquareDistance( uvFirst2 );
+
+    return (( uvLast1.SquareDistance( uvFirst2 ) < tol2 ) ||
+            ( uvLast1.SquareDistance( uvLast2 ) < tol2 ));
+  }
+
+  //================================================================================
+  /*!
+   * \brief Compose TSideVector for both FACEs keeping matching order of EDGEs
+   */
+  //================================================================================
+
+  bool getWires(const TopoDS_Face&                 tgtFace,
+                const TopoDS_Face&                 srcFace,
+                SMESH_Mesh *                       tgtMesh,
+                SMESH_Mesh *                       srcMesh,
+                const TAssocTool::TShapeShapeMap&  shape2ShapeMap,
+                TSideVector&                       srcWires,
+                TSideVector&                       tgtWires,
+                const bool                         is1DComputed)
+  {
+    // get ordered src EDGEs
+    TError err;
+    srcWires = StdMeshers_FaceSide::GetFaceWires( srcFace, *srcMesh,/*skipMediumNodes=*/0, err);
+    if ( err && !err->IsOK() )
+      return false;
+
+    // make corresponding sequence of tgt EDGEs
+    tgtWires.resize( srcWires.size() );
+    for ( size_t iW = 0; iW < srcWires.size(); ++iW )
+    {
+      list< TopoDS_Edge > tgtEdges;
+      StdMeshers_FaceSidePtr srcWire = srcWires[iW];
+      TopTools_IndexedMapOfShape edgeMap; // to detect seam edges
+      for ( int iE = 0; iE < srcWire->NbEdges(); ++iE )
+      {
+        TopoDS_Edge E = TopoDS::Edge( shape2ShapeMap( srcWire->Edge( iE ), /*isSrc=*/true));
+        // reverse a seam edge encountered for the second time
+        const int index = edgeMap.Add( E );
+        if ( index < edgeMap.Extent() ) // E is a seam
+        {
+          // check which of edges to reverse, E or one already being in tgtEdges
+          if ( are2dConnected( tgtEdges.back(), E, tgtFace ))
+          {
+            list< TopoDS_Edge >::iterator eIt = tgtEdges.begin();
+            std::advance( eIt, index-1 );
+            eIt->Reverse();
+          }
+          else
+          {
+            E.Reverse();
+          }
+        }
+        if ( srcWire->NbEdges() == 1 && tgtMesh == srcMesh ) // circle
+        {
+          // try to verify ori by propagation
+          pair<int,TopoDS_Edge> nE = StdMeshers_ProjectionUtils::GetPropagationEdge
+            ( srcMesh, E, srcWire->Edge( iE ));
+          if ( !nE.second.IsNull() )
+            E = nE.second;
+        }
+        tgtEdges.push_back( E );
+      }
+      tgtWires[ iW ].reset( new StdMeshers_FaceSide( tgtFace, tgtEdges, tgtMesh,
+                                                     /*theIsForward = */ true,
+                                                     /*theIgnoreMediumNodes = */false));
+      if ( is1DComputed &&
+           srcWires[iW]->GetUVPtStruct().size() !=
+           tgtWires[iW]->GetUVPtStruct().size())
+        return false;
+    }
+    return true;
+  }
+                             
+  //================================================================================
+  /*!
    * \brief Preform projection in case if tgtFace.IsPartner( srcFace ) and in case
-   * if projection by transformation is possible
+   * if projection by 3D transformation is possible
    */
   //================================================================================
 
-  bool projectPartner(const TopoDS_Face&                tgtFace,
-                      const TopoDS_Face&                srcFace,
-                      SMESH_Mesh *                      tgtMesh,
-                      SMESH_Mesh *                      srcMesh,
-                      const TAssocTool::TShapeShapeMap& shape2ShapeMap)
+  bool projectPartner(const TopoDS_Face&                 tgtFace,
+                      const TopoDS_Face&                 srcFace,
+                      const TSideVector&                 tgtWires,
+                      const TSideVector&                 srcWires,
+                      const TAssocTool::TShapeShapeMap&  shape2ShapeMap,
+                      TAssocTool::TNodeNodeMap&          src2tgtNodes)
   {
+    SMESH_Mesh *    tgtMesh = tgtWires[0]->GetMesh();
+    SMESH_Mesh *    srcMesh = srcWires[0]->GetMesh();
     SMESHDS_Mesh* tgtMeshDS = tgtMesh->GetMeshDS();
     SMESHDS_Mesh* srcMeshDS = srcMesh->GetMeshDS();
+    SMESH_MesherHelper helper( *tgtMesh );
 
     const double tol = 1.e-7 * srcMeshDS->getMaxDim();
 
-    gp_Trsf trsf; // transformation to get location of target nodes from source ones
+    // transformation to get location of target nodes from source ones
+    StdMeshers_ProjectionUtils::TrsfFinder3D trsf;
     if ( tgtFace.IsPartner( srcFace ))
     {
       gp_Trsf srcTrsf = srcFace.Location();
       gp_Trsf tgtTrsf = tgtFace.Location();
-      trsf = srcTrsf.Inverted() * tgtTrsf;
+      trsf.Set( srcTrsf.Inverted() * tgtTrsf );
     }
     else
     {
-      // Try to find the transformation
-
-      // make any local coord systems of src and tgt faces
-      vector<gp_Pnt> srcPP, tgtPP; // 3 points on face boundaries to make axes of CS
-      int tgtNbVert = SMESH_MesherHelper::Count( tgtFace, TopAbs_VERTEX, /*ignoreSame=*/true );
-      int srcNbVert = SMESH_MesherHelper::Count( srcFace, TopAbs_VERTEX, /*ignoreSame=*/true );
-      SMESH_subMesh *         srcSM = srcMesh->GetSubMesh( srcFace );
-      SMESH_subMeshIteratorPtr smIt = srcSM->getDependsOnIterator(/*includeSelf=*/false,false);
-      srcSM = smIt->next(); // sm of a vertex
-      while ( smIt->more() && srcPP.size() < 3 )
+      // Try to find the 3D transformation
+
+      const int totNbSeg = 50;
+      vector< gp_XYZ > srcPnts, tgtPnts;
+      srcPnts.reserve( totNbSeg );
+      tgtPnts.reserve( totNbSeg );
+      for ( size_t iW = 0; iW < srcWires.size(); ++iW )
       {
-        srcSM = smIt->next();
-        SMESHDS_SubMesh* srcSmds = srcSM->GetSubMeshDS();
-        if ( !srcSmds ) continue;
-        SMDS_NodeIteratorPtr nIt = srcSmds->GetNodes();
-        while ( nIt->more() )
+        const double minSegLen = srcWires[iW]->Length() / totNbSeg;
+        for ( int iE = 0; iE < srcWires[iW]->NbEdges(); ++iE )
         {
-          SMESH_TNodeXYZ p ( nIt->next());
-          bool pOK = false;
-          switch ( srcPP.size() )
-          {
-          case 0: pOK = true; break;
-
-          case 1: pOK = ( srcPP[0].SquareDistance( p ) > 10*tol ); break;
-            
-          case 2:
-            {
-              gp_Vec p0p1( srcPP[0], srcPP[1] ), p0p( srcPP[0], p );
-              // pOK = !p0p1.IsParallel( p0p, tol );
-              pOK = !p0p1.IsParallel( p0p, 3.14/20 ); // angle min 18 degrees
-              break;
-            }
-          }
-          if ( !pOK )
-            continue;
-
-          // find corresponding point on target shape
-          pOK = false;
-          gp_Pnt tgtP;
-          const TopoDS_Shape& tgtShape = shape2ShapeMap( srcSM->GetSubShape(), /*isSrc=*/true );
-          if ( tgtShape.ShapeType() == TopAbs_VERTEX )
+          int nbSeg    = Max( 1, int( srcWires[iW]->EdgeLength( iE ) / minSegLen ));
+          double srcU  = srcWires[iW]->FirstParameter( iE );
+          double tgtU  = tgtWires[iW]->FirstParameter( iE );
+          double srcDu = ( srcWires[iW]->LastParameter( iE )- srcU ) / nbSeg;
+          double tgtDu = ( tgtWires[iW]->LastParameter( iE )- tgtU ) / nbSeg;
+          for ( size_t i = 0; i < nbSeg; ++i  )
           {
-            tgtP = BRep_Tool::Pnt( TopoDS::Vertex( tgtShape ));
-            if ( srcNbVert == tgtNbVert || tgtPP.empty() )
-              pOK = true;
-            else
-              pOK = (( tgtP.Distance( tgtPP[0] ) > tol*tol ) &&
-                     ( tgtPP.size() == 1 || tgtP.Distance( tgtPP[1] ) > tol*tol ));
-            //cout << "V - nS " << p._node->GetID() << " - nT " << SMESH_Algo::VertexNode(TopoDS::Vertex( tgtShape),tgtMeshDS)->GetID() << endl;
+            srcPnts.push_back( srcWires[iW]->Value3d( srcU ).XYZ() );
+            tgtPnts.push_back( tgtWires[iW]->Value3d( tgtU ).XYZ() );
+            srcU += srcDu;
+            tgtU += tgtDu;
           }
-          else if ( tgtPP.size() > 0 )
-          {
-            if ( SMESHDS_SubMesh* tgtSmds = tgtMeshDS->MeshElements( tgtShape ))
-            {
-              double srcDist = srcPP[0].Distance( p );
-              double eTol = BRep_Tool::Tolerance( TopoDS::Edge( tgtShape ));
-              if (eTol < tol) eTol = tol;
-              SMDS_NodeIteratorPtr nItT = tgtSmds->GetNodes();
-              while ( nItT->more() && !pOK )
-              {
-                const SMDS_MeshNode* n = nItT->next();
-                tgtP = SMESH_TNodeXYZ( n );
-                pOK = ( fabs( srcDist - tgtPP[0].Distance( tgtP )) < 2*eTol );
-                //cout << "E - nS " << p._node->GetID() << " - nT " << n->GetID()<< " OK - " << pOK<< " " << fabs( srcDist - tgtPP[0].Distance( tgtP ))<< " tol " << eTol<< endl;
-              }
-            }
-          }
-          if ( !pOK )
-            continue;
-
-          srcPP.push_back( p );
-          tgtPP.push_back( tgtP );
         }
       }
-      if ( srcPP.size() != 3 )
+      if ( !trsf.Solve( srcPnts, tgtPnts ))
         return false;
 
-      // make transformation
-      gp_Trsf fromTgtCS, toSrcCS; // from/to global CS
-      gp_Ax2 srcCS( srcPP[0], gp_Vec( srcPP[0], srcPP[1] ), gp_Vec( srcPP[0], srcPP[2]));
-      gp_Ax2 tgtCS( tgtPP[0], gp_Vec( tgtPP[0], tgtPP[1] ), gp_Vec( tgtPP[0], tgtPP[2]));
-      toSrcCS  .SetTransformation( gp_Ax3( srcCS ));
-      fromTgtCS.SetTransformation( gp_Ax3( tgtCS ));
-      fromTgtCS.Invert();
+      // check trsf
 
-      trsf = fromTgtCS * toSrcCS;
+      bool trsfIsOK = true;
+      const int nbTestPnt = 20;
+      const size_t  iStep = Max( 1, int( srcPnts.size() / nbTestPnt ));
+      // check boundary
+      for ( size_t i = 0; ( i < srcPnts.size() && trsfIsOK ); i += iStep )
+      {
+        gp_Pnt trsfTgt = trsf.Transform( srcPnts[i] );
+        trsfIsOK = ( trsfTgt.SquareDistance( tgtPnts[i] ) < tol*tol );
+      }
+      // check an in-FACE point
+      if ( trsfIsOK )
+      {
+        BRepAdaptor_Surface srcSurf( srcFace );
+        gp_Pnt srcP =
+          srcSurf.Value( 0.5 * ( srcSurf.FirstUParameter() + srcSurf.LastUParameter() ),
+                         0.5 * ( srcSurf.FirstVParameter() + srcSurf.LastVParameter() ));
+        gp_Pnt tgtTrsfP = trsf.Transform( srcP );
+        TopLoc_Location loc;
+        GeomAPI_ProjectPointOnSurf& proj = helper.GetProjector( tgtFace, loc, 0.1*tol );
+        if ( !loc.IsIdentity() )
+          tgtTrsfP.Transform( loc.Transformation().Inverted() );
+        proj.Perform( tgtTrsfP );
+        trsfIsOK = ( proj.IsDone() &&
+                     proj.NbPoints() > 0 &&
+                     proj.LowerDistance() < tol );
+      }
+      if ( !trsfIsOK )
+        return false;
     }
 
     // Fill map of src to tgt nodes with nodes on edges
 
-    map<const SMDS_MeshNode* , const SMDS_MeshNode*> src2tgtNodes;
-    map<const SMDS_MeshNode* , const SMDS_MeshNode*>::iterator srcN_tgtN;
+    src2tgtNodes.clear();
+    TAssocTool::TNodeNodeMap::iterator srcN_tgtN;
 
     bool tgtEdgesMeshed = false;
     for ( TopExp_Explorer srcExp( srcFace, TopAbs_EDGE); srcExp.More(); srcExp.Next() )
@@ -518,31 +592,6 @@ namespace {
           )
         return false;
 
-      if ( !tgtEdge.IsPartner( srcEdge ))
-      {
-        if ( tgtNodes.empty() )
-          return false;
-        // check that transformation is OK by three nodes
-        gp_Pnt p0S = SMESH_TNodeXYZ( (srcNodes.begin())  ->second);
-        gp_Pnt p1S = SMESH_TNodeXYZ( (srcNodes.rbegin()) ->second);
-        gp_Pnt p2S = SMESH_TNodeXYZ( (++srcNodes.begin())->second);
-
-        gp_Pnt p0T = SMESH_TNodeXYZ( (tgtNodes.begin())  ->second);
-        gp_Pnt p1T = SMESH_TNodeXYZ( (tgtNodes.rbegin()) ->second);
-        gp_Pnt p2T = SMESH_TNodeXYZ( (++tgtNodes.begin())->second);
-
-        // transform source points, they must coincide with target ones
-        if ( p0T.SquareDistance( p0S.Transformed( trsf )) > tol ||
-             p1T.SquareDistance( p1S.Transformed( trsf )) > tol ||
-             p2T.SquareDistance( p2S.Transformed( trsf )) > tol )
-        {
-          //cout << "KO trsf, 3 dist: "
-          //<< p0T.SquareDistance( p0S.Transformed( trsf ))<< ", "
-          //<< p1T.SquareDistance( p1S.Transformed( trsf ))<< ", "
-          //<< p2T.SquareDistance( p2S.Transformed( trsf ))<< ", "<<endl;
-          return false;
-        }
-      }
       if ( !tgtNodes.empty() )
       {
         map< double, const SMDS_MeshNode* >::iterator u_tn = tgtNodes.begin();
@@ -563,16 +612,6 @@ namespace {
       if ( !tgtN || tgtV.ShapeType() != TopAbs_VERTEX )
         return false;
 
-      if ( !tgtV.IsPartner( srcV ))
-      {
-        // check that transformation is OK by three nodes
-        gp_Pnt p0S = SMESH_TNodeXYZ( srcN );
-        gp_Pnt p0T = SMESH_TNodeXYZ( tgtN );
-        if ( p0T.SquareDistance( p0S.Transformed( trsf )) > tol )
-        {
-          return false;
-        }
-      }
       src2tgtNodes.insert( make_pair( srcN, tgtN ));
     }
 
@@ -580,7 +619,6 @@ namespace {
     // Make new faces
 
     // prepare the helper to adding quadratic elements if necessary
-    SMESH_MesherHelper helper( *tgtMesh );
     helper.SetSubShape( tgtFace );
     helper.IsQuadraticSubMesh( tgtFace );
 
@@ -594,7 +632,7 @@ namespace {
     const SMDS_MeshNode* nullNode = 0;
 
     // indices of nodes to create properly oriented faces
-    bool isReverse = ( trsf.Form() != gp_Identity );
+    bool isReverse = ( !trsf.IsIdentity() );
     int tri1 = 1, tri2 = 2, quad1 = 1, quad3 = 3;
     if ( isReverse )
       std::swap( tri1, tri2 ), std::swap( quad1, quad3 );
@@ -614,7 +652,7 @@ namespace {
         if ( srcN_tgtN->second == nullNode )
         {
           // create a new node
-          gp_Pnt tgtP = gp_Pnt( SMESH_TNodeXYZ( srcNode )).Transformed( trsf );
+          gp_Pnt tgtP = trsf.Transform( SMESH_TNodeXYZ( srcNode ));
           SMDS_MeshNode* n = helper.AddNode( tgtP.X(), tgtP.Y(), tgtP.Z() );
           srcN_tgtN->second = n;
           switch ( srcNode->GetPosition()->GetTypeOfPosition() )
@@ -700,93 +738,27 @@ namespace {
 
   //================================================================================
   /*!
-   * \brief Check if two consecutive EDGEs are connected in 2D
-   *  \param [in] E1 - a well oriented non-seam EDGE
-   *  \param [in] E2 - a possibly well oriented seam EDGE
-   *  \param [in] F - a FACE
-   *  \return bool - result
-   */
-  //================================================================================
-
-  bool are2dConnected( const TopoDS_Edge & E1,
-                       const TopoDS_Edge & E2,
-                       const TopoDS_Face & F )
-  {
-    double f,l;
-    Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( E1, F, f, l );
-    gp_Pnt2d uvLast1 = c1->Value( E1.Orientation() == TopAbs_REVERSED ? f : l );
-
-    Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( E2, F, f, l );
-    gp_Pnt2d uvFirst2 = c2->Value( f );
-    gp_Pnt2d uvLast2  = c2->Value( l );
-    double tol2 = 1e-5 * uvLast2.SquareDistance( uvFirst2 );
-
-    return (( uvLast1.SquareDistance( uvFirst2 ) < tol2 ) ||
-            ( uvLast1.SquareDistance( uvLast2 ) < tol2 ));
-  }
-
-  //================================================================================
-  /*!
    * \brief Preform projection in case if the faces are similar in 2D space
    */
   //================================================================================
 
-  bool projectBy2DSimilarity(const TopoDS_Face&                tgtFace,
-                             const TopoDS_Face&                srcFace,
-                             SMESH_Mesh *                      tgtMesh,
-                             SMESH_Mesh *                      srcMesh,
-                             const TAssocTool::TShapeShapeMap& shape2ShapeMap,
-                             const bool                        is1DComputed)
+  bool projectBy2DSimilarity(const TopoDS_Face&                 tgtFace,
+                             const TopoDS_Face&                 srcFace,
+                             const TSideVector&                 tgtWires,
+                             const TSideVector&                 srcWires,
+                             const TAssocTool::TShapeShapeMap&  shape2ShapeMap,
+                             const bool                         is1DComputed,
+                             TAssocTool::TNodeNodeMap&          src2tgtNodes)
   {
-    // 1) Preparation
+    SMESH_Mesh * tgtMesh = tgtWires[0]->GetMesh();
+    SMESH_Mesh * srcMesh = srcWires[0]->GetMesh();
 
-    // get ordered src EDGEs
-    TError err;
-    TSideVector srcWires =
-      StdMeshers_FaceSide::GetFaceWires( srcFace, *srcMesh,/*ignoreMediumNodes = */false, err);
-    if ( err && !err->IsOK() )
-      return false;
+    // WARNING: we can have problems if the FACE is symmetrical in 2D,
+    // then the projection can be mirrored relating to what is expected
 
-    // make corresponding sequence of tgt EDGEs
-    TSideVector tgtWires( srcWires.size() );
-    for ( size_t iW = 0; iW < srcWires.size(); ++iW )
-    {
-      list< TopoDS_Edge > tgtEdges;
-      StdMeshers_FaceSidePtr srcWire = srcWires[iW];
-      TopTools_IndexedMapOfShape edgeMap; // to detect seam edges
-      for ( int iE = 0; iE < srcWire->NbEdges(); ++iE )
-      {
-        TopoDS_Edge E = TopoDS::Edge( shape2ShapeMap( srcWire->Edge( iE ), /*isSrc=*/true));
-        // reverse a seam edge encountered for the second time
-        const int index = edgeMap.Add( E );
-        if ( index < edgeMap.Extent() ) // E is a seam
-        {
-          // check which of edges to reverse, E or one already being in tgtEdges
-          if ( are2dConnected( tgtEdges.back(), E, tgtFace ))
-          {
-            list< TopoDS_Edge >::iterator eIt = tgtEdges.begin();
-            std::advance( eIt, index-1 );
-            eIt->Reverse();
-          }
-          else
-          {
-            E.Reverse();
-          }
-        }
-        tgtEdges.push_back( E );
-      }
-      tgtWires[ iW ].reset( new StdMeshers_FaceSide( tgtFace, tgtEdges, tgtMesh,
-                                                     /*theIsForward = */ true,
-                                                     /*theIgnoreMediumNodes = */false));
-      if ( is1DComputed &&
-           srcWires[iW]->GetUVPtStruct().size() !=
-           tgtWires[iW]->GetUVPtStruct().size())
-        return false;
-    }
-
-    // 2) Find transformation
+    // 1) Find 2D transformation
 
-    gp_Trsf2d trsf;
+    StdMeshers_ProjectionUtils::TrsfFinder2D trsf;
     {
       // get 2 pairs of corresponding UVs
       gp_Pnt2d srcP0 = srcWires[0]->Value2d(0.0);
@@ -801,34 +773,74 @@ namespace {
       toSrcCS  .SetTransformation( srcCS );
       fromTgtCS.SetTransformation( tgtCS );
       fromTgtCS.Invert();
-
-      trsf = fromTgtCS * toSrcCS;
+      trsf.Set( fromTgtCS * toSrcCS );
 
       // check transformation
+      bool trsfIsOK = true;
       const double tol = 1e-5 * gp_Vec2d( srcP0, srcP1 ).Magnitude();
-      for ( double u = 0.12; u < 1.; u += 0.1 )
+      for ( double u = 0.12; ( u < 1. && trsfIsOK ); u += 0.1 )
+      {
+        gp_Pnt2d srcUV  = srcWires[0]->Value2d( u );
+        gp_Pnt2d tgtUV  = tgtWires[0]->Value2d( u );
+        gp_Pnt2d tgtUV2 = trsf.Transform( srcUV );
+        trsfIsOK = ( tgtUV.Distance( tgtUV2 ) < tol );
+      }
+
+      // Find trsf using a least-square approximation
+      if ( !trsfIsOK )
       {
-        gp_Pnt2d srcUV = srcWires[0]->Value2d( u );
-        gp_Pnt2d tgtUV = tgtWires[0]->Value2d( u );
-        gp_Pnt2d tgtUV2 = srcUV.Transformed( trsf );
-        if ( tgtUV.Distance( tgtUV2 ) > tol )
+        // find trsf
+        const int totNbSeg = 50;
+        vector< gp_XY > srcPnts, tgtPnts;
+        srcPnts.resize( totNbSeg );
+        tgtPnts.resize( totNbSeg );
+        for ( size_t iW = 0; iW < srcWires.size(); ++iW )
+        {
+          const double minSegLen = srcWires[iW]->Length() / totNbSeg;
+          for ( int iE = 0; iE < srcWires[iW]->NbEdges(); ++iE )
+          {
+            int nbSeg    = Max( 1, int( srcWires[iW]->EdgeLength( iE ) / minSegLen ));
+            double srcU  = srcWires[iW]->FirstParameter( iE );
+            double tgtU  = tgtWires[iW]->FirstParameter( iE );
+            double srcDu = ( srcWires[iW]->LastParameter( iE )- srcU ) / nbSeg;
+            double tgtDu = ( tgtWires[iW]->LastParameter( iE )- tgtU ) / nbSeg;
+            for ( size_t i = 0; i < nbSeg; ++i, srcU += srcDu, tgtU += tgtDu  )
+            {
+              srcPnts.push_back( srcWires[iW]->Value2d( srcU ).XY() );
+              tgtPnts.push_back( tgtWires[iW]->Value2d( tgtU ).XY() );
+            }
+          }
+        }
+        if ( !trsf.Solve( srcPnts, tgtPnts ))
+          return false;
+
+        // check trsf
+
+        trsfIsOK = true;
+        const int nbTestPnt = 10;
+        const size_t  iStep = Max( 1, int( srcPnts.size() / nbTestPnt ));
+        for ( size_t i = 0; ( i < srcPnts.size() && trsfIsOK ); i += iStep )
+        {
+          gp_Pnt2d trsfTgt = trsf.Transform( srcPnts[i] );
+          trsfIsOK = ( trsfTgt.Distance( tgtPnts[i] ) < tol );
+        }
+        if ( !trsfIsOK )
           return false;
       }
-    }
+    } // "Find transformation" block
 
-    // 3) Projection
+    // 2) Projection
 
-    typedef map<const SMDS_MeshNode* , const SMDS_MeshNode*, TIDCompare> TN2NMap;
-    TN2NMap src2tgtNodes;
-    TN2NMap::iterator srcN_tgtN;
+    src2tgtNodes.clear();
+    TAssocTool::TNodeNodeMap::iterator srcN_tgtN;
 
     // fill src2tgtNodes in with nodes on EDGEs
-    for ( unsigned iW = 0; iW < srcWires.size(); ++iW )
+    for ( size_t iW = 0; iW < srcWires.size(); ++iW )
       if ( is1DComputed )
       {
         const vector<UVPtStruct>& srcUVs = srcWires[iW]->GetUVPtStruct();
         const vector<UVPtStruct>& tgtUVs = tgtWires[iW]->GetUVPtStruct();
-        for ( unsigned i = 0; i < srcUVs.size(); ++i )
+        for ( size_t i = 0; i < srcUVs.size(); ++i )
           src2tgtNodes.insert( make_pair( srcUVs[i].node, tgtUVs[i].node ));
       }
       else
@@ -880,8 +892,8 @@ namespace {
           // create a new node
           gp_Pnt2d srcUV = srcHelper.GetNodeUV( srcFace, srcNode,
                                                 elem->GetNode( helper.WrapIndex(i+1,nbN)), &uvOK);
-          gp_Pnt2d tgtUV = srcUV.Transformed( trsf );
-          gp_Pnt   tgtP  = tgtSurface->Value( tgtUV.X(), tgtUV.Y() );
+          gp_Pnt2d   tgtUV = trsf.Transform( srcUV );
+          gp_Pnt      tgtP = tgtSurface->Value( tgtUV.X(), tgtUV.Y() );
           SMDS_MeshNode* n = tgtMeshDS->AddNode( tgtP.X(), tgtP.Y(), tgtP.Z() );
           switch ( srcNode->GetPosition()->GetTypeOfPosition() )
           {
@@ -920,6 +932,74 @@ namespace {
 
   } // bool projectBy2DSimilarity(...)
 
+  //================================================================================
+  /*!
+   * \brief Fix bad faces by smoothing
+   */
+  //================================================================================
+
+  void fixDistortedFaces( SMESH_MesherHelper& helper )
+  {
+    // Detect bad faces
+
+    bool haveBadFaces = false;
+
+    const TopoDS_Face&  F = TopoDS::Face( helper.GetSubShape() );
+    SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( F );
+    if ( !smDS || smDS->NbElements() == 0 ) return;
+
+    SMDS_ElemIteratorPtr faceIt = smDS->GetElements();
+    double prevArea2D = 0;
+    vector< const SMDS_MeshNode* > nodes;
+    vector< gp_XY >                uv;
+    while ( faceIt->more() && !haveBadFaces )
+    {
+      const SMDS_MeshElement* face = faceIt->next();
+
+      // get nodes
+      nodes.resize( face->NbCornerNodes() );
+      SMDS_MeshElement::iterator n = face->begin_nodes();
+      for ( size_t i = 0; i < nodes.size(); ++n, ++i )
+        nodes[ i ] = *n;
+
+      // get UVs
+      const SMDS_MeshNode* inFaceNode = 0;
+      if ( helper.HasSeam() )
+        for ( size_t i = 0; ( i < nodes.size() && !inFaceNode ); ++i )
+          if ( !helper.IsSeamShape( nodes[ i ]->getshapeId() ))
+            inFaceNode = nodes[ i ];
+
+      uv.resize( nodes.size() );
+      for ( size_t i = 0; i < nodes.size(); ++i )
+        uv[ i ] = helper.GetNodeUV( F, nodes[ i ], inFaceNode );
+
+      // compare orientation of triangles
+      for ( int iT = 0, nbT = nodes.size()-2; iT < nbT; ++iT )
+      {
+        gp_XY v1 = uv[ iT+1 ] - uv[ 0 ];
+        gp_XY v2 = uv[ iT+2 ] - uv[ 0 ];
+        double area2D = v2 ^ v1;
+        if (( haveBadFaces = ( area2D * prevArea2D < 0 )))
+          break;
+        prevArea2D = area2D;
+      }
+    }
+
+    // Fix faces
+
+    if ( haveBadFaces )
+    {
+      SMESH_MeshEditor editor( helper.GetMesh() );
+
+      TIDSortedElemSet faces;
+      for ( faceIt = smDS->GetElements(); faceIt->more(); )
+        faces.insert( faces.end(), faceIt->next() );
+
+      set<const SMDS_MeshNode*> fixedNodes;
+      editor.Smooth( faces, fixedNodes, SMESH_MeshEditor::CENTROIDAL, 5 );
+    }
+  }
+  
 } // namespace
 
 
@@ -930,6 +1010,8 @@ namespace {
 
 bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theShape)
 {
+  _src2tgtNodes.clear();
+
   MESSAGE("Projection_2D Compute");
   if ( !_sourceHypo )
     return false;
@@ -1001,17 +1083,25 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape&
       is1DComputed = sm->IsMeshComputed();
   }
 
+  // get ordered src and tgt EDGEs
+  TSideVector srcWires, tgtWires;
+  if ( !getWires( tgtFace, srcFace, tgtMesh, srcMesh,
+                  shape2ShapeMap, srcWires, tgtWires, is1DComputed ))
+    return false;
+
   bool done = false;
 
 if ( !done )
+ if ( !done )
   {
     // try to project from the same face with different location
-    done = projectPartner( tgtFace, srcFace, tgtMesh, srcMesh, shape2ShapeMap );
+    done = projectPartner( tgtFace, srcFace, tgtWires, srcWires,
+                           shape2ShapeMap, _src2tgtNodes );
   }
   if ( !done )
   {
     // projection in case if the faces are similar in 2D space
-    done = projectBy2DSimilarity( tgtFace, srcFace, tgtMesh, srcMesh, shape2ShapeMap, is1DComputed);
+    done = projectBy2DSimilarity( tgtFace, srcFace, tgtWires, srcWires,
+                                  shape2ShapeMap, is1DComputed, _src2tgtNodes);
   }
 
   SMESH_MesherHelper helper( theMesh );
@@ -1019,6 +1109,7 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape&
 
   if ( !done )
   {
+    _src2tgtNodes.clear();
     // --------------------
     // Prepare to mapping 
     // --------------------
@@ -1100,7 +1191,7 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape&
       return error(COMPERR_BAD_INPUT_MESH,"Can't load mesh pattern from the source face");
 
     // --------------------
-    // Perform 2D mapping 
+    // Perform 2D mapping
     // --------------------
 
     // Compute mesh on a target face
@@ -1278,6 +1369,14 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape&
     if ( nbFaceBeforeMerge != nbFaceAtferMerge && !helper.HasDegeneratedEdges() )
       return error(COMPERR_BAD_INPUT_MESH, "Probably invalid node parameters on geom faces");
 
+
+    // ----------------------------------------------------------------
+    // The mapper can create distorted faces by placing nodes out of the FACE
+    // boundary -- fix bad faces by smoothing
+    // ----------------------------------------------------------------
+
+    fixDistortedFaces( helper );
+
     // ----------------------------------------------------------------
     // The mapper can't create quadratic elements, so convert if needed
     // ----------------------------------------------------------------
index 57e34db..878d218 100644 (file)
@@ -30,6 +30,7 @@
 #include "SMESH_StdMeshers.hxx"
 
 #include "SMESH_Algo.hxx"
+#include "StdMeshers_ProjectionUtils.hxx"
 
 class StdMeshers_ProjectionSource2D;
 
@@ -59,9 +60,12 @@ public:
    */
   virtual void SetEventListener(SMESH_subMesh* whenSetToSubMesh);
   
-protected:
 
-  const StdMeshers_ProjectionSource2D* _sourceHypo;
+ protected:
+
+  const StdMeshers_ProjectionSource2D*     _sourceHypo;
+
+  StdMeshers_ProjectionUtils::TNodeNodeMap _src2tgtNodes;
 
 };
 
index 3079c7b..93969df 100644 (file)
@@ -307,7 +307,7 @@ bool StdMeshers_Projection_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aS
       }
     }
     // Find matching nodes of tgt and src faces
-    TNodeNodeMap faceMatchingNodes;
+    TAssocTool::TNodeNodeMap faceMatchingNodes;
     if ( ! TAssocTool::FindMatchingNodesOnFaces( srcFace, srcMesh, tgtFace, tgtMesh, 
                                                  shape2ShapeMap, faceMatchingNodes ))
       return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #")
index 1741141..a06cf86 100644 (file)
@@ -239,7 +239,7 @@ bool StdMeshers_RadialPrism_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& a
     }
 
     // Find matching nodes of in and out faces
-    TNodeNodeMap nodeIn2OutMap;
+    ProjectionUtils::TNodeNodeMap nodeIn2OutMap;
     if ( ! ProjectionUtils::FindMatchingNodesOnFaces( inFace, &aMesh, outFace, &aMesh,
                                                       shape2ShapeMap, nodeIn2OutMap ))
       return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #")
index c04a632..eeec08f 100644 (file)
@@ -101,7 +101,7 @@ namespace VISCOUS_3D
   const double theSmoothThickToElemSizeRatio = 0.3;
 
   // what part of thickness is allowed till intersection
-  // defined by SALOME_TESTS/Grids/smesh/viscous_layers_00/A5
+  // (defined by SALOME_TESTS/Grids/smesh/viscous_layers_00/A5)
   const double theThickToIntersection = 1.5;
 
   bool needSmoothing( double cosin, double tgtThick, double elemSize )