Salome HOME
0021543: EDF 1978 SMESH: Viscous layer for 2D meshes
[modules/smesh.git] / src / StdMeshers / StdMeshers_Regular_1D.cxx
index f4fa4a18e96020bc26912fcb8ae98d0190737bc6..6ea520001c195e66985a5b9e329d9c9ff0e0c099 100644 (file)
@@ -1,30 +1,30 @@
-//  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2012  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
+// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 //
-//  This library is free software; you can redistribute it and/or
-//  modify it under the terms of the GNU Lesser General Public
-//  License as published by the Free Software Foundation; either
-//  version 2.1 of the License.
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
 //
-//  This library is distributed in the hope that it will be useful,
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-//  Lesser General Public License for more details.
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
 //
-//  You should have received a copy of the GNU Lesser General Public
-//  License along with this library; if not, write to the Free Software
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
-//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
-//  SMESH SMESH : implementaion of SMESH idl descriptions
+
 //  File   : StdMeshers_Regular_1D.cxx
 //           Moved here from SMESH_Regular_1D.cxx
 //  Author : Paul RASCLE, EDF
 //  Module : SMESH
-
+//
 #include "StdMeshers_Regular_1D.hxx"
 #include "StdMeshers_Distribution.hxx"
 
@@ -61,6 +61,7 @@
 #include <TopExp_Explorer.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Edge.hxx>
+#include <TopoDS_Vertex.hxx>
 
 #include <string>
 #include <limits>
@@ -118,9 +119,8 @@ bool StdMeshers_Regular_1D::CheckHypothesis
   _hypType = NONE;
   _quadraticMesh = false;
 
-  const bool ignoreAuxiliaryHyps = false;
   const list <const SMESHDS_Hypothesis * > & hyps =
-    GetUsedHypothesis(aMesh, aShape, ignoreAuxiliaryHyps);
+    GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliaryHyps=*/false);
 
   // find non-auxiliary hypothesis
   const SMESHDS_Hypothesis *theHyp = 0;
@@ -325,6 +325,8 @@ static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last,
       return false;
     prevU = U;
   }
+  if ( theReverse )
+    theParams.reverse();
   return true;
 }
 
@@ -352,7 +354,7 @@ static void compensateError(double a1, double an,
                             bool              adjustNeighbors2an = false)
 {
   int i, nPar = theParams.size();
-  if ( a1 + an < length && nPar > 1 )
+  if ( a1 + an <= length && nPar > 1 )
   {
     bool reverse = ( U1 > Un );
     GCPnts_AbscissaPoint Discret(C3d, reverse ? an : -an, Un);
@@ -373,14 +375,23 @@ static void compensateError(double a1, double an,
       dUn = Utgt - theParams.back();
     }
 
-    double q  = dUn / ( nPar - 1 );
-    if ( !adjustNeighbors2an ) {
-      for ( itU = theParams.rbegin(), i = 1; i < nPar; itU++, i++ ) {
+    if ( !adjustNeighbors2an )
+    {
+      double q = dUn / ( Utgt - Un ); // (signed) factor of segment length change
+      for ( itU = theParams.rbegin(), i = 1; i < nPar; i++ ) {
+        double prevU = *itU;
         (*itU) += dUn;
-        dUn -= q;
+        ++itU;
+        dUn = q * (*itU - prevU) * (prevU-U1)/(Un-U1);
       }
     }
-    else {
+    else if ( nPar == 1 )
+    {
+      theParams.back() += dUn;
+    }
+    else
+    {
+      double q  = dUn / ( nPar - 1 );
       theParams.back() += dUn;
       double sign = reverse ? -1 : 1;
       double prevU = theParams.back();
@@ -602,12 +613,14 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
   case NB_SEGMENTS: {
 
     double eltSize = 1;
+    int nbSegments;
     if ( _hypType == MAX_LENGTH )
     {
       double nbseg = ceil(theLength / _value[ BEG_LENGTH_IND ]); // integer sup
       if (nbseg <= 0)
         nbseg = 1;                        // degenerated edge
       eltSize = theLength / nbseg;
+      nbSegments = (int) nbseg;
     }
     else if ( _hypType == LOCAL_LENGTH )
     {
@@ -624,7 +637,7 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
           bool computed = sm->IsMeshComputed();
           if (!computed) {
             if (sm->GetComputeState() == SMESH_subMesh::READY_TO_COMPUTE) {
-              sm->ComputeStateEngine(SMESH_subMesh::COMPUTE);
+              _gen->Compute( theMesh, _mainEdge, /*anUpward=*/true);
               computed = sm->IsMeshComputed();
             }
           }
@@ -648,13 +661,14 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
       if (nbseg <= 0)
         nbseg = 1;                        // degenerated edge
       eltSize = theLength / nbseg;
+      nbSegments = (int) nbseg;
     }
     else
     {
       // Number Of Segments hypothesis
-      int NbSegm = _ivalue[ NB_SEGMENTS_IND ];
-      if ( NbSegm < 1 )  return false;
-      if ( NbSegm == 1 ) return true;
+      nbSegments = _ivalue[ NB_SEGMENTS_IND ];
+      if ( nbSegments < 1 )  return false;
+      if ( nbSegments == 1 ) return true;
 
       switch (_ivalue[ DISTR_TYPE_IND ])
       {
@@ -664,8 +678,8 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
 
           if (fabs(scale - 1.0) < Precision::Confusion()) {
             // special case to avoid division by zero
-            for (int i = 1; i < NbSegm; i++) {
-              double param = f + (l - f) * i / NbSegm;
+            for (int i = 1; i < nbSegments; i++) {
+              double param = f + (l - f) * i / nbSegments;
               theParams.push_back( param );
             }
           } else {
@@ -673,14 +687,22 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
             if ( theReverse )
               scale = 1.0 / scale;
 
-            double alpha = pow(scale, 1.0 / (NbSegm - 1));
-            double factor = (l - f) / (1.0 - pow(alpha, NbSegm));
+            double alpha = pow(scale, 1.0 / (nbSegments - 1));
+            double factor = (l - f) / (1.0 - pow(alpha, nbSegments));
 
-            for (int i = 1; i < NbSegm; i++) {
+            for (int i = 1; i < nbSegments; i++) {
               double param = f + factor * (1.0 - pow(alpha, i));
               theParams.push_back( param );
             }
           }
+          const double lenFactor = theLength/(l-f);
+          list<double>::iterator u = theParams.begin(), uEnd = theParams.end();
+          for ( ; u != uEnd; ++u )
+          {
+            GCPnts_AbscissaPoint Discret( theC3d, ((*u)-f) * lenFactor, f );
+            if ( Discret.IsDone() )
+              *u = Discret.Parameter();
+          }
           return true;
         }
         break;
@@ -701,7 +723,7 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
         }
         break;
       case StdMeshers_NumberOfSegments::DT_Regular:
-        eltSize = theLength / _ivalue[ NB_SEGMENTS_IND ];
+        eltSize = theLength / nbSegments;
         break;
       default:
         return false;
@@ -711,13 +733,13 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
     if ( !Discret.IsDone() )
       return error( "GCPnts_UniformAbscissa failed");
 
-    int NbPoints = Discret.NbPoints();
-    for ( int i = 2; i < NbPoints; i++ )
+    int NbPoints = Min( Discret.NbPoints(), nbSegments + 1 );
+    for ( int i = 2; i < NbPoints; i++ ) // skip 1st and last points
     {
       double param = Discret.Parameter(i);
       theParams.push_back( param );
     }
-    compensateError( eltSize, eltSize, f, l, theLength, theC3d, theParams ); // for PAL9899
+    compensateError( eltSize, eltSize, f, l, theLength, theC3d, theParams, true ); // for PAL9899
     return true;
   }
 
@@ -728,6 +750,9 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
     double a1 = _value[ BEG_LENGTH_IND ];
     double an = _value[ END_LENGTH_IND ];
     double q  = ( theLength - a1 ) / ( theLength - an );
+    if ( q < theLength/1e6 || 1.01*theLength < a1 + an)
+      return error ( SMESH_Comment("Invalid segment lengths (")<<a1<<" and "<<an<<") "<<
+                     "for an edge of length "<<theLength);
 
     double U1 = theReverse ? l : f;
     double Un = theReverse ? f : l;
@@ -756,6 +781,9 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
 
     double a1 = _value[ BEG_LENGTH_IND ];
     double an = _value[ END_LENGTH_IND ];
+    if ( 1.01*theLength < a1 + an)
+      return error ( SMESH_Comment("Invalid segment lengths (")<<a1<<" and "<<an<<") "<<
+                     "for an edge of length "<<theLength);
 
     double  q = ( an - a1 ) / ( 2 *theLength/( a1 + an ) - 1 );
     int n = int(fabs(q) > numeric_limits<double>::min() ? ( 1+( an-a1 )/q ) : ( 1+theLength/a1 ));
@@ -938,23 +966,43 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & theMesh, const TopoDS_Shape & t
   if (!idFirst || !idLast)
     return error( COMPERR_BAD_INPUT_MESH, "No node on vertex");
 
+  // remove elements created by e.g. patern mapping (PAL21999)
+  // CLEAN event is incorrectly ptopagated seemingly due to Propagation hyp
+  // so TEMPORARY solution is to clean the submesh manually
+  //theMesh.GetSubMesh(theShape)->ComputeStateEngine( SMESH_subMesh::CLEAN );
+  if (SMESHDS_SubMesh * subMeshDS = meshDS->MeshElements(theShape))
+  {
+    SMDS_ElemIteratorPtr ite = subMeshDS->GetElements();
+    while (ite->more())
+      meshDS->RemoveFreeElement(ite->next(), subMeshDS);
+    SMDS_NodeIteratorPtr itn = subMeshDS->GetNodes();
+    while (itn->more()) {
+      const SMDS_MeshNode * node = itn->next();
+      if ( node->NbInverseElements() == 0 )
+        meshDS->RemoveFreeNode(node, subMeshDS);
+      else
+        meshDS->RemoveNode(node);
+    }
+  }
+
   if (!Curve.IsNull())
   {
     list< double > params;
     bool reversed = false;
     if ( theMesh.GetShapeToMesh().ShapeType() >= TopAbs_WIRE ) {
+      // if the shape to mesh is WIRE or EDGE
       reversed = ( EE.Orientation() == TopAbs_REVERSED );
     }
     if ( !_mainEdge.IsNull() ) {
+      // take into account reversing the edge the hypothesis is propagated from
       reversed = ( _mainEdge.Orientation() == TopAbs_REVERSED );
+      int mainID = meshDS->ShapeToIndex(_mainEdge);
+      if ( std::find( _revEdgesIDs.begin(), _revEdgesIDs.end(), mainID) != _revEdgesIDs.end())
+        reversed = !reversed;
     }
-    else if ( _revEdgesIDs.size() > 0 ) {
-      for ( int i = 0; i < _revEdgesIDs.size(); i++) {
-        if ( _revEdgesIDs[i] == shapeID ) {
-          reversed = !reversed;
-        }
-      }
-    }
+    // take into account this edge reversing
+    if ( std::find( _revEdgesIDs.begin(), _revEdgesIDs.end(), shapeID) != _revEdgesIDs.end())
+      reversed = !reversed;
 
     BRepAdaptor_Curve C3d( E );
     double length = EdgeLength( E );
@@ -979,7 +1027,6 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & theMesh, const TopoDS_Shape & t
       parLast = f;
     }
     */
-
     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++) {
       double param = *itU;
       gp_Pnt P = Curve->Value(param);
@@ -1155,10 +1202,9 @@ StdMeshers_Regular_1D::GetUsedHypothesis(SMESH_Mesh &         aMesh,
 
   SMESH_HypoFilter auxiliaryFilter, compatibleFilter;
   auxiliaryFilter.Init( SMESH_HypoFilter::IsAuxiliary() );
-  const bool ignoreAux = true;
-  InitCompatibleHypoFilter( compatibleFilter, ignoreAux );
+  InitCompatibleHypoFilter( compatibleFilter, /*ignoreAux=*/true );
 
-  // get non-auxiliary assigned to aShape
+  // get non-auxiliary assigned directly to aShape
   int nbHyp = aMesh.GetHypotheses( aShape, compatibleFilter, _usedHypList, false );
 
   if (nbHyp == 0 && aShape.ShapeType() == TopAbs_EDGE)