Salome HOME
PR: relax constraints on node distances on StdMeshers_import_1D
[modules/smesh.git] / src / StdMeshers / StdMeshers_Prism_3D.cxx
index 00a53610546f39a9b6cf983e50dba7ed1accbc78..cac08865fea808fdd0de25153d6d43e93fd9fe6c 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -68,9 +68,15 @@ using namespace std;
 // gp_Pnt p (xyz); \
 // cout << msg << " ("<< p.X() << "; " <<p.Y() << "; " <<p.Z() << ") " <<endl;\
 // }
+#ifdef _DEBUG_
+#define DBGOUT(msg) //cout << msg << endl;
+#else
+#define DBGOUT(msg)
+#endif
+
+namespace TAssocTool = StdMeshers_ProjectionUtils;
 
-typedef StdMeshers_ProjectionUtils TAssocTool;
-typedef SMESH_Comment              TCom;
+typedef SMESH_Comment TCom;
 
 enum { ID_BOT_FACE = SMESH_Block::ID_Fxy0,
        ID_TOP_FACE = SMESH_Block::ID_Fxy1,
@@ -457,20 +463,13 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh
   if ( nbSolids < 1 )
     return true;
 
-  Prism_3D::TPrismTopo prism;
-
-  if ( nbSolids == 1 )
-  {
-    return ( initPrism( prism, TopExp_Explorer( theShape, TopAbs_SOLID ).Current() ) &&
-             compute( prism ));
-  }
-
   TopTools_IndexedDataMapOfShapeListOfShape faceToSolids;
   TopExp::MapShapesAndAncestors( theShape, TopAbs_FACE, TopAbs_SOLID, faceToSolids );
 
   // look for meshed FACEs ("source" FACEs) that must be prism bottoms
-  list< TopoDS_Face > meshedFaces;//, notQuadMeshedFaces, notQuadFaces;
+  list< TopoDS_Face > meshedFaces, notQuadMeshedFaces, notQuadFaces;
   const bool meshHasQuads = ( theMesh.NbQuadrangles() > 0 );
+  //StdMeshers_Quadrangle_2D* quadAlgo = TQuadrangleAlgo::instance( this );
   for ( int iF = 1; iF < faceToSolids.Extent(); ++iF )
   {
     const TopoDS_Face& face = TopoDS::Face( faceToSolids.FindKey( iF ));
@@ -479,17 +478,36 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh
     {
       if ( !meshHasQuads ||
            !helper.IsSameElemGeometry( faceSM->GetSubMeshDS(), SMDSGeom_QUADRANGLE ) ||
-           !helper.IsStructured( faceSM ))
-        // notQuadMeshedFaces are of higher priority
+           !helper.IsStructured( faceSM )
+           )
+        notQuadMeshedFaces.push_front( face );
+      else if ( myHelper->Count( face, TopAbs_EDGE, /*ignoreSame=*/false ) != 4 )
         meshedFaces.push_front( face );
       else
         meshedFaces.push_back( face );
     }
+    // not add not quadrilateral FACE as we can't compute it
+    // else if ( !quadAlgo->CheckNbEdges( theMesh, face ))
+    // // not add not quadrilateral FACE as it can be a prism side
+    // // else if ( myHelper->Count( face, TopAbs_EDGE, /*ignoreSame=*/false ) != 4 )
+    // {
+    //   notQuadFaces.push_back( face );
+    // }
   }
-  //meshedFaces.splice( meshedFaces.begin(), notQuadMeshedFaces );
+  // notQuadFaces are of medium priority, put them before ordinary meshed faces
+  meshedFaces.splice( meshedFaces.begin(), notQuadFaces );
+  // notQuadMeshedFaces are of highest priority, put them before notQuadFaces
+  meshedFaces.splice( meshedFaces.begin(), notQuadMeshedFaces );
 
-  // if ( meshedFaces.empty() )
-  //   return error( COMPERR_BAD_INPUT_MESH, "No meshed source faces found" );
+  Prism_3D::TPrismTopo prism;
+
+  if ( nbSolids == 1 )
+  {
+    if ( !meshedFaces.empty() )
+      prism.myBottom = meshedFaces.front();
+    return ( initPrism( prism, TopExp_Explorer( theShape, TopAbs_SOLID ).Current() ) &&
+             compute( prism ));
+  }
 
   TopTools_MapOfShape meshedSolids;
   list< Prism_3D::TPrismTopo > meshedPrism;
@@ -619,6 +637,27 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh
     // TODO. there are other ways to find out the source FACE:
     // propagation, topological similarity, ect.
 
+    // simply try to mesh all not meshed SOLIDs
+    if ( meshedFaces.empty() )
+    {
+      for ( TopExp_Explorer solid( theShape, TopAbs_SOLID ); solid.More(); solid.Next() )
+      {
+        mySetErrorToSM = false;
+        prism.Clear();
+        if ( !meshedSolids.Contains( solid.Current() ) &&
+             initPrism( prism, solid.Current() ))
+        {
+          mySetErrorToSM = true;
+          if ( !compute( prism ))
+            return false;
+          meshedFaces.push_front( prism.myTop );
+          meshedFaces.push_front( prism.myBottom );
+          meshedPrism.push_back( prism );
+          meshedSolids.Add( solid.Current() );
+        }
+        mySetErrorToSM = true;
+      }
+    }
 
     if ( meshedFaces.empty() ) // set same error to 10 not-computed solids
     {
@@ -633,7 +672,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh
           SMESH_subMesh* sm = theMesh.GetSubMesh( solid.Current() );
           sm->GetComputeError() = err;
         }
-      return false;
+      return error( err );
     }
   }
   return true;
@@ -708,7 +747,7 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
   // -------------------------
 
   // Compose a vector of indixes of right neighbour FACE for each wall FACE
-  // that is not so evident in case of several WIREs
+  // that is not so evident in case of several WIREs in the bottom FACE
   thePrism.myRightQuadIndex.clear();
   for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
     thePrism.myRightQuadIndex.push_back( i+1 );
@@ -1013,6 +1052,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
 {
   SMESH_Mesh*     mesh = myHelper->GetMesh();
   SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
+  DBGOUT( endl << "COMPUTE Prism " << meshDS->ShapeToIndex( thePrism.myShape3D ));
 
   TProjction1dAlgo* projector1D = TProjction1dAlgo::instance( this );
   StdMeshers_Quadrangle_2D* quadAlgo = TQuadrangleAlgo::instance( this, myHelper );
@@ -1078,6 +1118,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
         const TopoDS_Edge& srcE = lftSide->Edge(i);
         SMESH_subMesh*    srcSM = mesh->GetSubMesh( srcE );
         if ( !srcSM->IsMeshComputed() ) {
+          DBGOUT( "COMPUTE V edge " << srcSM->GetId() );
           srcSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
           srcSM->ComputeStateEngine       ( SMESH_subMesh::COMPUTE );
           if ( !srcSM->IsMeshComputed() )
@@ -1112,7 +1153,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
         }
         continue;
       }
-      // Compute
+      // Compute 'vertical projection'
       if ( nbTgtMeshed == 0 )
       {
         // compute nodes on target VERTEXes
@@ -1131,8 +1172,9 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
         }
 
         // compute nodes on target EDGEs
+        DBGOUT( "COMPUTE V edge (proj) " << shapeID( lftSide->Edge(0)));
         rgtSide->Reverse(); // direct it same as the lftSide
-        myHelper->SetElementsOnShape( false );
+        myHelper->SetElementsOnShape( false ); // myHelper holds the prism shape
         TopoDS_Edge tgtEdge;
         for ( size_t iN = 1; iN < srcNodeStr.size()-1; ++iN ) // add nodes
         {
@@ -1143,25 +1185,24 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
         }
         for ( size_t iN = 1; iN < srcNodeStr.size(); ++iN ) // add segments
         {
-          SMDS_MeshElement* newEdge = myHelper->AddEdge( newNodes[ iN-1 ], newNodes[ iN ] );
+          // find an EDGE to set a new segment
           std::pair<int, TopAbs_ShapeEnum> id2type = 
             myHelper->GetMediumPos( newNodes[ iN-1 ], newNodes[ iN ] );
-          if ( id2type.second == TopAbs_EDGE )
-          {
-            meshDS->SetMeshElementOnShape( newEdge, id2type.first );
-          }
-          else // new nodes are on different EDGEs; put one of them on VERTEX
+          if ( id2type.second != TopAbs_EDGE )
           {
+            // new nodes are on different EDGEs; put one of them on VERTEX
             const int      edgeIndex = rgtSide->EdgeIndex( srcNodeStr[ iN-1 ].normParam );
             const double vertexParam = rgtSide->LastParameter( edgeIndex );
             const gp_Pnt           p = BRep_Tool::Pnt( rgtSide->LastVertex( edgeIndex ));
             const int         isPrev = ( Abs( srcNodeStr[ iN-1 ].normParam - vertexParam ) <
                                          Abs( srcNodeStr[ iN   ].normParam - vertexParam ));
-            meshDS->SetMeshElementOnShape( newEdge, newNodes[ iN-(1-isPrev) ]->getshapeId() );
             meshDS->UnSetNodeOnShape( newNodes[ iN-isPrev ] );
             meshDS->SetNodeOnVertex ( newNodes[ iN-isPrev ], rgtSide->LastVertex( edgeIndex ));
-            meshDS->MoveNode( newNodes[ iN-isPrev ], p.X(), p.Y(), p.Z() );
+            meshDS->MoveNode        ( newNodes[ iN-isPrev ], p.X(), p.Y(), p.Z() );
+            id2type.first = newNodes[ iN-(1-isPrev) ]->getshapeId();
           }
+          SMDS_MeshElement* newEdge = myHelper->AddEdge( newNodes[ iN-1 ], newNodes[ iN ] );
+          meshDS->SetMeshElementOnShape( newEdge, id2type.first );
         }
         myHelper->SetElementsOnShape( true );
         for ( int i = 0; i < rgtSide->NbEdges(); ++i ) // update state of sub-meshes
@@ -1198,28 +1239,43 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
       // Top EDGEs must be projections from the bottom ones
       // to compute stuctured quad mesh on wall FACEs
       // ---------------------------------------------------
-      const TopoDS_Edge& botE = (*quad)->side[ QUAD_BOTTOM_SIDE ]->Edge(0);
-      const TopoDS_Edge& topE = (*quad)->side[ QUAD_TOP_SIDE    ]->Edge(0);
-
-      projector1D->myHyp.SetSourceEdge( botE );
-
-      SMESH_subMesh* tgtEdgeSm = mesh->GetSubMesh( topE );
-      if ( !tgtEdgeSm->IsMeshComputed() )
       {
-        // compute nodes on VERTEXes
-        tgtEdgeSm->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
-        // project segments
-        projector1D->InitComputeError();
-        bool ok = projector1D->Compute( *mesh, topE );
-        if ( !ok )
+        const TopoDS_Edge& botE = (*quad)->side[ QUAD_BOTTOM_SIDE ]->Edge(0);
+        const TopoDS_Edge& topE = (*quad)->side[ QUAD_TOP_SIDE    ]->Edge(0);
+        SMESH_subMesh*    botSM = mesh->GetSubMesh( botE );
+        SMESH_subMesh*    topSM = mesh->GetSubMesh( topE );
+        SMESH_subMesh*    srcSM = botSM;
+        SMESH_subMesh*    tgtSM = topSM;
+        if ( !srcSM->IsMeshComputed() && topSM->IsMeshComputed() )
+          std::swap( srcSM, tgtSM );
+
+        if ( !srcSM->IsMeshComputed() )
         {
-          SMESH_ComputeErrorPtr err = projector1D->GetComputeError();
-          if ( err->IsOK() ) err->myName = COMPERR_ALGO_FAILED;
-          tgtEdgeSm->GetComputeError() = err;
-          return false;
+          DBGOUT( "COMPUTE H edge " << srcSM->GetId());
+          srcSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE ); // nodes on VERTEXes
+          srcSM->ComputeStateEngine( SMESH_subMesh::COMPUTE );        // segments on the EDGE
         }
+        srcSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+
+        if ( !tgtSM->IsMeshComputed() )
+        {
+          // compute nodes on VERTEXes
+          tgtSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
+          // project segments
+          DBGOUT( "COMPUTE H edge (proj) " << tgtSM->GetId());
+          projector1D->myHyp.SetSourceEdge( TopoDS::Edge( srcSM->GetSubShape() ));
+          projector1D->InitComputeError();
+          bool ok = projector1D->Compute( *mesh, tgtSM->GetSubShape() );
+          if ( !ok )
+          {
+            SMESH_ComputeErrorPtr err = projector1D->GetComputeError();
+            if ( err->IsOK() ) err->myName = COMPERR_ALGO_FAILED;
+            tgtSM->GetComputeError() = err;
+            return false;
+          }
+        }
+        tgtSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
       }
-      tgtEdgeSm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
 
       // Compute quad mesh on wall FACEs
       // -------------------------------
@@ -1234,6 +1290,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
                               "Not all edges have valid algorithm and hypothesis"));
         // mesh the <face>
         quadAlgo->InitComputeError();
+        DBGOUT( "COMPUTE Quad face " << fSM->GetId());
         bool ok = quadAlgo->Compute( *mesh, face );
         fSM->GetComputeError() = quadAlgo->GetComputeError();
         if ( !ok )
@@ -1511,7 +1568,12 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top()
   SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS();
 
   if ( !botSMDS || botSMDS->NbElements() == 0 )
-    return toSM( error(TCom("No elememts on face #") << botSM->GetId() ));
+  {
+    _gen->Compute( *myHelper->GetMesh(), botSM->GetSubShape() );
+    botSMDS = botSM->GetSubMeshDS();
+    if ( !botSMDS || botSMDS->NbElements() == 0 )
+      return toSM( error(TCom("No elements on face #") << botSM->GetId() ));
+  }
 
   bool needProject = !topSM->IsMeshComputed();
   if ( !needProject && 
@@ -1646,7 +1708,7 @@ bool StdMeshers_Prism_3D::projectBottomToTop()
   // if the bottom faces is orienetd OK then top faces must be reversed
   bool reverseTop = true;
   if ( myHelper->NbAncestors( botFace, *myBlock.Mesh(), TopAbs_SOLID ) > 1 )
-    reverseTop = ! SMESH_Algo::IsReversedSubMesh( TopoDS::Face( botFace ), meshDS );
+    reverseTop = ! myHelper->IsReversedSubMesh( TopoDS::Face( botFace ));
   int iFrw, iRev, *iPtr = &( reverseTop ? iRev : iFrw );
 
   // loop on bottom mesh faces
@@ -2160,8 +2222,10 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
     {
       if ( len2edgeMap.size() != nbEdges )
         RETURN_BAD_RESULT("Uniqueness of edge lengths not assured");
-      map< double, int >::reverse_iterator maxLen_i = len2edgeMap.rbegin();
-      map< double, int >::reverse_iterator midLen_i = ++len2edgeMap.rbegin();
+
+      multimap< double, int >::reverse_iterator maxLen_i = len2edgeMap.rbegin();
+      multimap< double, int >::reverse_iterator midLen_i = ++len2edgeMap.rbegin();
+
       double maxLen = maxLen_i->first;
       double midLen = ( len2edgeMap.size() == 1 ) ? 0 : midLen_i->first;
       switch ( nbEdges ) {