Salome HOME
Join modifications from branch OCC_development_for_3_2_0a2
[modules/smesh.git] / src / StdMeshers / StdMeshers_Regular_1D.cxx
index 1723d6549b3627d60e9388ded200925cdddf6982..dc1ec9bb42064f558075805889cb3727060f847c 100644 (file)
 //  Module : SMESH
 //  $Header$
 
+using namespace std;
+
 #include "StdMeshers_Regular_1D.hxx"
 #include "StdMeshers_Distribution.hxx"
 #include "SMESH_Gen.hxx"
 #include "SMESH_Mesh.hxx"
+#include "SMESH_HypoFilter.hxx"
+#include "SMESH_subMesh.hxx"
 
 #include "StdMeshers_LocalLength.hxx"
 #include "StdMeshers_NumberOfSegments.hxx"
@@ -42,7 +46,6 @@
 #include "SMDS_MeshElement.hxx"
 #include "SMDS_MeshNode.hxx"
 #include "SMDS_EdgePosition.hxx"
-#include "SMESH_subMesh.hxx"
 
 #include "Utils_SALOME_Exception.hxx"
 #include "utilities.h"
@@ -88,6 +91,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 +112,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")
@@ -352,17 +372,27 @@ bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge
       {
       case StdMeshers_NumberOfSegments::DT_Scale:
         {
+          int NbSegm   = _ivalue[ NB_SEGMENTS_IND ];
           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 +420,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 +499,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 +556,25 @@ 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 ( ! computeInternalParameters( E, params, reversed )) {
+        //cout << "computeInternalParameters() failed" <<endl;
         return false;
+      }
     }
     catch ( Standard_Failure ) {
+      //cout << "computeInternalParameters() failed, Standard_Failure" <<endl;
       return false;
     }
 
@@ -552,9 +582,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 +597,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 +639,38 @@ 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.;
+        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);
+      }
       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.;
+      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);
+    }
   }
   return true;
 }
@@ -604,40 +681,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, false );
     }
   }
-  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;
 }