Salome HOME
Merge from V5_1_main branch 24/11/2010
[modules/smesh.git] / src / StdMeshers / StdMeshers_Regular_1D.cxx
index 477b809f9b5628e2c4ab0918061bfba28900d661..951414c53b7673d1b0d19bfdde8f9ef3d0461053 100644 (file)
@@ -1,4 +1,4 @@
-//  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
+//  Copyright (C) 2007-2010  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
 //
 //  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"
 
@@ -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;
 }
 
@@ -374,10 +376,14 @@ static void compensateError(double a1, double an,
     }
 
     double q  = dUn / ( nPar - 1 );
-    if ( !adjustNeighbors2an ) {
-      for ( itU = theParams.rbegin(), i = 1; i < nPar; itU++, i++ ) {
+    if ( !adjustNeighbors2an )
+    {
+      q = Abs( dUn / ( Utgt - Un )); // 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 {
@@ -624,7 +630,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();
             }
           }
@@ -728,6 +734,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 +765,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 ));
@@ -804,38 +816,83 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
       }
       if(!IsExist) Params.InsertBefore(j,aPnts[i]);
     }
-    double pf, pl, par2, par1, psize;
-    if (theReverse) {
-      pf = l;
-      pl = f;
+    double par2, par1, lp;
+    par1 = f;
+    lp = l;
+    double sign = 1.0;
+    if(theReverse) {
+      par1 = l;
+      lp = f;
+      sign = -1.0;
     }
-    else {
-      pf = f;
-      pl = l;
-    }
-    psize = pl - pf;
-    par1 = pf;
-    //cout<<"aPnts.size() = "<<aPnts.size()<<"  Params.Length() = "
-    //    <<Params.Length()<<"   nbsegs.size() = "<<nbsegs.size()<<endl;
+    double eltSize, segmentSize = 0.;
+    double currAbscissa = 0;
     for(i=0; i<Params.Length(); i++) {
-      par2 = pf + Params.Value(i+1)*psize;
       int nbseg = ( i > nbsegs.size()-1 ) ? nbsegs[0] : nbsegs[i];
-      double dp = (par2-par1)/nbseg;
-      int j = 1;
-      for(; j<=nbseg; j++) {
-        double param = par1 + dp*j;
-        theParams.push_back( param );
+      segmentSize = Params.Value(i+1)*theLength - currAbscissa;
+      currAbscissa += segmentSize;
+      GCPnts_AbscissaPoint APnt(theC3d, sign*segmentSize, par1);
+      if( !APnt.IsDone() )
+        return error( "GCPnts_AbscissaPoint failed");
+      par2 = APnt.Parameter();
+      eltSize = segmentSize/nbseg;
+      GCPnts_UniformAbscissa Discret(theC3d, eltSize, par1, par2);
+      if(theReverse)
+        Discret.Initialize(theC3d, eltSize, par2, par1);
+      else
+        Discret.Initialize(theC3d, eltSize, par1, par2);
+      if ( !Discret.IsDone() )
+        return error( "GCPnts_UniformAbscissa failed");
+      int NbPoints = Discret.NbPoints();
+      list<double> tmpParams;
+      for(int i=2; i<NbPoints; i++) {
+        double param = Discret.Parameter(i);
+        tmpParams.push_back( param );
       }
+      if (theReverse) {
+        compensateError( eltSize, eltSize, par2, par1, segmentSize, theC3d, tmpParams );
+        tmpParams.reverse();
+      }
+      else {
+        compensateError( eltSize, eltSize, par1, par2, segmentSize, theC3d, tmpParams );
+      }
+      list<double>::iterator itP = tmpParams.begin();
+      for(; itP != tmpParams.end(); itP++) {
+        theParams.push_back( *(itP) );
+      }
+      theParams.push_back( par2 );
+
       par1 = par2;
     }
     // add for last
     int nbseg = ( nbsegs.size() > Params.Length() ) ? nbsegs[Params.Length()] : nbsegs[0];
-    double dp = (pl-par1)/nbseg;
-    int j = 1;
-    for(; j<nbseg; j++) {
-      double param = par1 + dp*j;
-      theParams.push_back( param );
+    segmentSize = theLength - currAbscissa;
+    eltSize = segmentSize/nbseg;
+    GCPnts_UniformAbscissa Discret;
+    if(theReverse)
+      Discret.Initialize(theC3d, eltSize, par1, lp);
+    else
+      Discret.Initialize(theC3d, eltSize, lp, par1);
+    if ( !Discret.IsDone() )
+      return error( "GCPnts_UniformAbscissa failed");
+    int NbPoints = Discret.NbPoints();
+    list<double> tmpParams;
+    for(int i=2; i<NbPoints; i++) {
+      double param = Discret.Parameter(i);
+      tmpParams.push_back( param );
     }
+    if (theReverse) {
+      compensateError( eltSize, eltSize, lp, par1, segmentSize, theC3d, tmpParams );
+      tmpParams.reverse();
+    }
+    else {
+      compensateError( eltSize, eltSize, par1, lp, segmentSize, theC3d, tmpParams );
+    }
+    list<double>::iterator itP = tmpParams.begin();
+    for(; itP != tmpParams.end(); itP++) {
+      theParams.push_back( *(itP) );
+    }
+
     if (theReverse) {
       theParams.reverse(); // NPAL18025
     }
@@ -893,23 +950,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 );
@@ -934,7 +1011,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);
@@ -1110,10 +1186,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)