23562: EDF 17098 - problem with Extrusion 3D

[modules/smesh.git] / src / StdMeshers / StdMeshers_QuadToTriaAdaptor.cxx

diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx
index f133994fa1d1041d097d51ad43a7dfe833175ff5..fd3a5bb10f8376a3dbb7789039c964905b7facd3 100644 (file)
--- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx
+++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx
@@ -296,12 +296,13 @@ namespace
       SMESH_ComputeErrorPtr& err = sm->GetComputeError();
       if ( !err || err->IsOK() )
       {
-        err = SMESH_ComputeError::New( COMPERR_BAD_INPUT_MESH, msg, sm->GetAlgo() );
-        err->myBadElements.push_back( face1 );
-        err->myBadElements.push_back( face2 );
+        SMESH_BadInputElements* badElems =
+          new SMESH_BadInputElements( mesh.GetMeshDS(),COMPERR_BAD_INPUT_MESH, msg, sm->GetAlgo() );
+        badElems->add( face1 );
+        badElems->add( face2 );
+        err.reset( badElems );
       }
     }
-    //throw SALOME_Exception( msg.c_str() );
 
     return false; // == "algo fails"
   }
@@ -485,9 +486,6 @@ static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
   gp_XYZ vert1 = P2.XYZ();
   gp_XYZ vert2 = P3.XYZ();
 
-  /* calculate distance from vert0 to ray origin */
-  gp_XYZ  tvec = orig - vert0;
-
   gp_XYZ edge1 = vert1 - vert0;
   gp_XYZ edge2 = vert2 - vert0;
 
@@ -497,9 +495,13 @@ static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
   /* if determinant is near zero, ray lies in plane of triangle */
   double det = edge1 * pvec;
 
-  if (det > -EPSILON && det < EPSILON)
+  const double ANGL_EPSILON = 1e-12;
+  if ( det > -ANGL_EPSILON && det < ANGL_EPSILON )
     return false;
 
+  /* calculate distance from vert0 to ray origin */
+  gp_XYZ  tvec = orig - vert0;
+
   /* calculate U parameter and test bounds */
   double u = ( tvec * pvec ) / det;
   //if (u < 0.0 || u > 1.0)
@@ -615,7 +617,7 @@ bool StdMeshers_QuadToTriaAdaptor::LimitHeight (gp_Pnt&
         }
       }
     }
-    if ( height < 1e-5 * idealHeight && intFace )
+    if ( height < 1e-2 * idealHeight && intFace )
       return overlapError( aMesh, NotCheckedFace, intFace, Shape );
   }
 
@@ -1012,6 +1014,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
       groupDS = 0;
   }
 
+  const bool toFindVolumes = aMesh.NbVolumes() > 0;
+
   vector<const SMDS_MeshElement*> myPyramids;
   SMESH_MesherHelper helper(aMesh);
   helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh());
@@ -1023,13 +1027,16 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
   if ( !myElemSearcher )
     myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS );
   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
+  SMESHUtils::Deleter<SMESH_ElementSearcher>
+    volSearcher( SMESH_MeshAlgos::GetElementSearcher( *meshDS ));
+  vector< const SMDS_MeshElement* > suspectFaces, foundVolumes;
 
   TColgp_Array1OfPnt PN(1,5);
   TColgp_Array1OfVec VN(1,4);
   vector<const SMDS_MeshNode*> FNodes(5);
   TColgp_SequenceOfPnt aContour;
 
-  SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(/*idInceasingOrder=*/true);
+  SMDS_FaceIteratorPtr fIt = meshDS->facesIterator();
   while( fIt->more())
   {
     const SMDS_MeshElement* face = fIt->next();
@@ -1042,7 +1049,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
     if ( what == NOT_QUAD )
       continue;
     if ( volumes[0] && volumes[1] )
-      continue; // face is shared by two volumes - no space for a pyramid
+      continue; // face is shared by two volumes - no room for a pyramid
 
     if ( what == DEGEN_QUAD )
     {
@@ -1143,6 +1150,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
     }
 
     // Restrict pyramid height by intersection with other faces
+
     gp_Vec tmpDir(PC,PCbest); tmpDir.Normalize();
     double tmp = PN(1).Distance(PN(3)) + PN(2).Distance(PN(4));
     // far points: in (PC, PCbest) direction and vice-versa
@@ -1154,8 +1162,24 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
     gp_Pnt intPnt     [2];
     int    intFaceInd [2] = { 0, 0 };
 
+    if ( toFindVolumes && 0 ) // non-conformal mesh is not suitable for any mesher so far
+    {
+      // there are volumes in the mesh, in a non-conformal mesh an neighbor
+      // volume can be not found yet
+      for ( int isRev = 0; isRev < 2; ++isRev )
+      {
+        if ( volumes[isRev] ) continue;
+        gp_Pnt testPnt = PC.XYZ() + tmpDir.XYZ() * height * ( isRev ? -0.1: 0.1 );
+        foundVolumes.clear();
+        if ( volSearcher->FindElementsByPoint( testPnt, SMDSAbs_Volume, foundVolumes ))
+          volumes[isRev] = foundVolumes[0];
+      }
+      if ( volumes[0] && volumes[1] )
+        continue; // no room for a pyramid
+    }
+
     gp_Ax1 line( PC, tmpDir );
-    vector< const SMDS_MeshElement* > suspectFaces;
+    suspectFaces.clear();
     searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces);
 
     for ( size_t iF = 0; iF < suspectFaces.size(); ++iF )
@@ -1208,7 +1232,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
       if ( !intersected[isRev] ) continue;
       double pyramidH = Min( height, dist2int[isRev]/3. );
       gp_Pnt    Papex = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH);
-      if ( pyramidH < 1e-5 * height )
+      if ( pyramidH < 1e-2 * height )
         return overlapError( aMesh, face, suspectFaces[ intFaceInd[isRev] ] );
 
       if ( !LimitHeight( Papex, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/false ))