Salome HOME
Join modifications from branch BR_For_OCT_611: migration to OCCT6.1.1 with new except...
[modules/smesh.git] / src / StdMeshers / StdMeshers_Regular_1D.cxx
index 1723d6549b3627d60e9388ded200925cdddf6982..d933aebbfd6e1c967fe7bc5c9a6c1af29b5445cc 100644 (file)
@@ -17,7 +17,7 @@
 //  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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org 
+//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
 //
 //
 //  Module : SMESH
 //  $Header$
 
+using namespace std;
+
 #include "StdMeshers_Regular_1D.hxx"
 #include "StdMeshers_Distribution.hxx"
-#include "SMESH_Gen.hxx"
-#include "SMESH_Mesh.hxx"
 
 #include "StdMeshers_LocalLength.hxx"
 #include "StdMeshers_NumberOfSegments.hxx"
 #include "StdMeshers_Deflection1D.hxx"
 #include "StdMeshers_AutomaticLength.hxx"
 
+#include "SMESH_Gen.hxx"
+#include "SMESH_Mesh.hxx"
+#include "SMESH_HypoFilter.hxx"
+#include "SMESH_subMesh.hxx"
+
 #include "SMDS_MeshElement.hxx"
 #include "SMDS_MeshNode.hxx"
 #include "SMDS_EdgePosition.hxx"
-#include "SMESH_subMesh.hxx"
 
 #include "Utils_SALOME_Exception.hxx"
 #include "utilities.h"
 #include <GCPnts_AbscissaPoint.hxx>
 #include <GCPnts_UniformAbscissa.hxx>
 #include <GCPnts_UniformDeflection.hxx>
-#include <Standard_ErrorHandler.hxx>
 #include <Precision.hxx>
 #include <Expr_GeneralExpression.hxx>
 #include <Expr_NamedUnknown.hxx>
 #include <Expr_Array1OfNamedUnknown.hxx>
-#include <TColStd_Array1OfReal.hxx>
 #include <ExprIntrp_GenExp.hxx>
+#include <TColStd_Array1OfReal.hxx>
 #include <OSD.hxx>
 
+#include <Standard_ErrorHandler.hxx>
+#include <Standard_Failure.hxx>
+
 #include <string>
 #include <math.h>
 
@@ -88,6 +94,8 @@ StdMeshers_Regular_1D::StdMeshers_Regular_1D(int hypId, int studyId,
        _compatibleHypothesis.push_back("Deflection1D");
        _compatibleHypothesis.push_back("Arithmetic1D");
        _compatibleHypothesis.push_back("AutomaticLength");
+
+       _compatibleHypothesis.push_back("QuadraticMesh"); // auxiliary !!!
 }
 
 //=============================================================================
@@ -107,22 +115,37 @@ StdMeshers_Regular_1D::~StdMeshers_Regular_1D()
 //=============================================================================
 
 bool StdMeshers_Regular_1D::CheckHypothesis
-                         (SMESH_Mesh& aMesh,
-                          const TopoDS_Shape& aShape,
+                         (SMESH_Mesh&                          aMesh,
+                          const TopoDS_Shape&                  aShape,
                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 {
   _hypType = NONE;
+  _quadraticMesh = false;
+
+  const bool ignoreAuxiliaryHyps = false;
+  const list <const SMESHDS_Hypothesis * > & hyps =
+    GetUsedHypothesis(aMesh, aShape, ignoreAuxiliaryHyps);
+
+  // find non-auxiliary hypothesis
+  const SMESHDS_Hypothesis *theHyp = 0;
+  list <const SMESHDS_Hypothesis * >::const_iterator h = hyps.begin();
+  for ( ; h != hyps.end(); ++h ) {
+    if ( static_cast<const SMESH_Hypothesis*>(*h)->IsAuxiliary() ) {
+      if ( strcmp( "QuadraticMesh", (*h)->GetName() ) == 0 )
+        _quadraticMesh = true;
+    }
+    else {
+      if ( !theHyp )
+        theHyp = *h; // use only the first non-auxiliary hypothesis
+    }
+  }
 
-  const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
-  if (hyps.size() == 0)
+  if ( !theHyp )
   {
     aStatus = SMESH_Hypothesis::HYP_MISSING;
     return false;  // can't work without a hypothesis
   }
 
-  // use only the first hypothesis
-  const SMESHDS_Hypothesis *theHyp = hyps.front();
-
   string hypName = theHyp->GetName();
 
   if (hypName == "LocalLength")
@@ -267,7 +290,8 @@ static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last,
                                int nbSeg, Function& func,
                                list<double>& theParams)
 {
-  OSD::SetSignal( true );
+  // never do this way
+  //OSD::SetSignal( true );
 
   if( nbSeg<=0 )
     return false;
@@ -327,7 +351,7 @@ bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge
 
   double f, l;
   Handle(Geom_Curve) Curve = BRep_Tool::Curve(theEdge, f, l);
-  GeomAdaptor_Curve C3d(Curve);
+  GeomAdaptor_Curve C3d (Curve, f, l);
 
   double length = EdgeLength(theEdge);
 
@@ -348,21 +372,34 @@ bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge
     else
     {
       // Number Of Segments hypothesis
+      int NbSegm = _ivalue[ NB_SEGMENTS_IND ];
+      if ( NbSegm < 1 )  return false;
+      if ( NbSegm == 1 ) return true;
+
       switch (_ivalue[ DISTR_TYPE_IND ])
       {
       case StdMeshers_NumberOfSegments::DT_Scale:
         {
           double scale = _value[ SCALE_FACTOR_IND ];
-          if ( theReverse )
-            scale = 1. / scale;
-          double alpha = pow( scale , 1.0 / (_ivalue[ NB_SEGMENTS_IND ] - 1));
-          double factor = (l - f) / (1 - pow( alpha,_ivalue[ NB_SEGMENTS_IND ]));
-
-          int i, NbPoints = 1 + _ivalue[ NB_SEGMENTS_IND ];
-          for ( i = 2; i < NbPoints; i++ )
-          {
-            double param = f + factor * (1 - pow(alpha, i - 1));
-            theParams.push_back( param );
+
+          if (fabs(scale - 1.0) < Precision::Confusion()) {
+            // special case to avoid division on zero
+            for (int i = 1; i < NbSegm; i++) {
+              double param = f + (l - f) * i / NbSegm;
+              theParams.push_back( param );
+            }
+          } else {
+            // general case of scale distribution
+            if ( theReverse )
+              scale = 1.0 / scale;
+
+            double alpha = pow(scale, 1.0 / (NbSegm - 1));
+            double factor = (l - f) / (1.0 - pow(alpha, NbSegm));
+
+            for (int i = 1; i < NbSegm; i++) {
+              double param = f + factor * (1.0 - pow(alpha, i));
+              theParams.push_back( param );
+            }
           }
           return true;
         }
@@ -390,7 +427,6 @@ bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge
         return false;
       }
     }
-
     GCPnts_UniformAbscissa Discret(C3d, eltSize, f, l);
     if ( !Discret.IsDone() )
       return false;
@@ -470,7 +506,7 @@ bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge
 
   case DEFLECTION: {
 
-    GCPnts_UniformDeflection Discret(C3d, _value[ DEFLECTION_IND ], true);
+    GCPnts_UniformDeflection Discret(C3d, _value[ DEFLECTION_IND ], f, l, true);
     if ( !Discret.IsDone() )
       return false;
 
@@ -527,24 +563,28 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh
 
   ASSERT(!VLast.IsNull());
   lid=aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
-  if (!lid->more())
-  {
+  if (!lid->more()) {
     MESSAGE (" NO NODE BUILT ON VERTEX ");
     return false;
   }
   const SMDS_MeshNode * idLast = lid->next();
 
-  if (!Curve.IsNull())
-  {
+  if (!Curve.IsNull()) {
     list< double > params;
     bool reversed = false;
     if ( !_mainEdge.IsNull() )
       reversed = aMesh.IsReversedInChain( EE, _mainEdge );
     try {
-      if ( ! computeInternalParameters( E, params, reversed ))
+#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
+      OCC_CATCH_SIGNALS;
+#endif
+      if ( ! computeInternalParameters( E, params, reversed )) {
+        //cout << "computeInternalParameters() failed" <<endl;
         return false;
+      }
     }
     catch ( Standard_Failure ) {
+      //cout << "computeInternalParameters() failed, Standard_Failure" <<endl;
       return false;
     }
 
@@ -552,9 +592,14 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh
     // only internal nodes receive an edge position with param on curve
 
     const SMDS_MeshNode * idPrev = idFirst;
+    double parPrev = f;
+    double parLast = l;
+//     if(reversed) {
+//       parPrev = l;
+//       parLast = f;
+//     }
     
-    for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
-    {
+    for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++) {
       double param = *itU;
       gp_Pnt P = Curve->Value(param);
 
@@ -562,17 +607,39 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh
       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
       meshDS->SetNodeOnEdge(node, shapeID, param);
 
-      SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
-      meshDS->SetMeshElementOnShape(edge, shapeID);
+      if(_quadraticMesh) {
+        // create medium node
+        double prm = ( parPrev + param )/2;
+        gp_Pnt PM = Curve->Value(prm);
+        SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
+        meshDS->SetNodeOnEdge(NM, shapeID, prm);
+        SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node, NM);
+        meshDS->SetMeshElementOnShape(edge, shapeID);
+      }
+      else {
+        SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
+        meshDS->SetMeshElementOnShape(edge, shapeID);
+      }
+
       idPrev = node;
+      parPrev = param;
+    }
+    if(_quadraticMesh) {
+      double prm = ( parPrev + parLast )/2;
+      gp_Pnt PM = Curve->Value(prm);
+      SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
+      meshDS->SetNodeOnEdge(NM, shapeID, prm);
+      SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast, NM);
+      meshDS->SetMeshElementOnShape(edge, shapeID);
+    }
+    else {
+      SMDS_MeshEdge* edge = meshDS->AddEdge(idPrev, idLast);
+      meshDS->SetMeshElementOnShape(edge, shapeID);
     }
-    SMDS_MeshEdge* edge = meshDS->AddEdge(idPrev, idLast);
-    meshDS->SetMeshElementOnShape(edge, shapeID);
   }
-  else
-  {
+  else {
     // Edge is a degenerated Edge : We put n = 5 points on the edge.
-    int NbPoints = 5;
+    const int NbPoints = 5;
     BRep_Tool::Range(E, f, l);
     double du = (l - f) / (NbPoints - 1);
     //MESSAGE("************* Degenerated edge! *****************");
@@ -582,18 +649,36 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh
     gp_Pnt P = BRep_Tool::Pnt(V1);
 
     const SMDS_MeshNode * idPrev = idFirst;
-    for (int i = 2; i < NbPoints; i++)
-    {
+    for (int i = 2; i < NbPoints; i++) {
       double param = f + (i - 1) * du;
       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
+      if(_quadraticMesh) {
+        // create medium node
+        double prm = param - du/2.;
+        SMDS_MeshNode * NM = meshDS->AddNode(P.X(), P.Y(), P.Z());
+        meshDS->SetNodeOnEdge(NM, shapeID, prm);
+        SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node, NM);
+        meshDS->SetMeshElementOnShape(edge, shapeID);
+      }
+      else {
+        SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
+        meshDS->SetMeshElementOnShape(edge, shapeID);
+      }
       meshDS->SetNodeOnEdge(node, shapeID, param);
-
-      SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
-      meshDS->SetMeshElementOnShape(edge, shapeID);
       idPrev = node;
     }
-    SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast);
-    meshDS->SetMeshElementOnShape(edge, shapeID);
+    if(_quadraticMesh) {
+      // create medium node
+      double prm = l - du/2.;
+      SMDS_MeshNode * NM = meshDS->AddNode(P.X(), P.Y(), P.Z());
+      meshDS->SetNodeOnEdge(NM, shapeID, prm);
+      SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast, NM);
+      meshDS->SetMeshElementOnShape(edge, shapeID);
+    }
+    else {
+      SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast);
+      meshDS->SetMeshElementOnShape(edge, shapeID);
+    }
   }
   return true;
 }
@@ -604,40 +689,47 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh
  */
 //=============================================================================
 
-const list <const SMESHDS_Hypothesis *> & StdMeshers_Regular_1D::GetUsedHypothesis(
-       SMESH_Mesh & aMesh, const TopoDS_Shape & aShape)
+const list <const SMESHDS_Hypothesis *> &
+StdMeshers_Regular_1D::GetUsedHypothesis(SMESH_Mesh &         aMesh,
+                                         const TopoDS_Shape & aShape,
+                                         const bool           ignoreAuxiliary)
 {
   _usedHypList.clear();
-  _usedHypList = GetAppliedHypothesis(aMesh, aShape);  // copy
-  int nbHyp = _usedHypList.size();
   _mainEdge.Nullify();
+
+  SMESH_HypoFilter auxiliaryFilter, compatibleFilter;
+  auxiliaryFilter.Init( SMESH_HypoFilter::IsAuxiliary() );
+  const bool ignoreAux = true;
+  InitCompatibleHypoFilter( compatibleFilter, ignoreAux );
+
+  // get non-auxiliary assigned to aShape
+  int nbHyp = aMesh.GetHypotheses( aShape, compatibleFilter, _usedHypList, false );
+
   if (nbHyp == 0)
   {
     // Check, if propagated from some other edge
     if (aShape.ShapeType() == TopAbs_EDGE &&
         aMesh.IsPropagatedHypothesis(aShape, _mainEdge))
     {
-      // Propagation of 1D hypothesis from <aMainEdge> on this edge
-      //_usedHypList = GetAppliedHypothesis(aMesh, _mainEdge); // copy
-      // use a general method in order not to nullify _mainEdge
-      _usedHypList = SMESH_Algo::GetUsedHypothesis(aMesh, _mainEdge);  // copy
-      nbHyp = _usedHypList.size();
+      // Propagation of 1D hypothesis from <aMainEdge> on this edge;
+      // get non-auxiliary assigned to _mainEdge
+      nbHyp = aMesh.GetHypotheses( _mainEdge, compatibleFilter, _usedHypList, true );
     }
   }
-  if (nbHyp == 0)
+
+  if (nbHyp == 0) // nothing propagated nor assigned to aShape
   {
-    TopTools_ListIteratorOfListOfShape ancIt( aMesh.GetAncestors( aShape ));
-    for (; ancIt.More(); ancIt.Next())
-    {
-      const TopoDS_Shape& ancestor = ancIt.Value();
-      _usedHypList = GetAppliedHypothesis(aMesh, ancestor);    // copy
-      nbHyp = _usedHypList.size();
-      if (nbHyp == 1)
-        break;
-    }
+    SMESH_Algo::GetUsedHypothesis( aMesh, aShape, ignoreAuxiliary );
+    nbHyp = _usedHypList.size();
   }
-  if (nbHyp > 1)
-    _usedHypList.clear();      //only one compatible hypothesis allowed
+  else
+  {
+    // get auxiliary hyps from aShape
+    aMesh.GetHypotheses( aShape, auxiliaryFilter, _usedHypList, true );
+  }
+  if ( nbHyp > 1 && ignoreAuxiliary )
+    _usedHypList.clear(); //only one compatible non-auxiliary hypothesis allowed
+
   return _usedHypList;
 }