Salome HOME
0020557: EDF 1151 SMESH: Netgen 2D fail to mesh a ring
authoreap <eap@opencascade.com>
Thu, 29 Oct 2009 08:36:01 +0000 (08:36 +0000)
committereap <eap@opencascade.com>
Thu, 29 Oct 2009 08:36:01 +0000 (08:36 +0000)
   and merge with V5

src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx

index a5090cfa8992d47524e8120cee19682c1fc3fde2..e259be1b07e4d013080087de68c4d8ca809f6141 100644 (file)
@@ -206,10 +206,6 @@ static TError AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
     StdMeshers_FaceSidePtr wire = wires[ iW ];
     const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
 
-    bool reverse = // 20526: [CEA] Disk meshing fails
-      ( wire->NbEdges() == 1 && 
-        geom.emap(geom.emap.FindIndex(wire->Edge(0))).Orientation() == TopAbs_REVERSED );
-
     int firstPointID = ngMesh.GetNP() + 1;
     int edgeID = 1, posID = -2;
     for ( int i = 0; i < wire->NbSegments(); ++i ) // loop on segments
@@ -260,14 +256,6 @@ static TError AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
         }
         seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
       }
-      // 20526: [CEA] Disk meshing fails
-      if (reverse)
-      {
-        swap (seg.p1, seg.p2);
-        swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
-        swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
-        swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
-      }
 
       ngMesh.AddSegment (seg);
 
@@ -456,3 +444,90 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
   return !err;
 }
 
+
+//=============================================================================
+/*!
+ *
+ */
+//=============================================================================
+
+bool NETGENPlugin_NETGEN_2D_ONLY::Evaluate(SMESH_Mesh& aMesh,
+                                           const TopoDS_Shape& aShape,
+                                           MapShapeNbElems& aResMap)
+{
+  TopoDS_Face F = TopoDS::Face(aShape);
+  if(F.IsNull())
+    return false;
+
+  // collect info from edges
+  int nb0d = 0, nb1d = 0;
+  bool IsQuadratic = false;
+  bool IsFirst = true;
+  double fullLen = 0.0;
+  TopTools_MapOfShape tmpMap;
+  for (TopExp_Explorer exp(F, TopAbs_EDGE); exp.More(); exp.Next()) {
+    TopoDS_Edge E = TopoDS::Edge(exp.Current());
+    if( tmpMap.Contains(E) )
+      continue;
+    tmpMap.Add(E);
+    SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
+    MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
+    if( anIt==aResMap.end() ) {
+      SMESH_subMesh *sm = aMesh.GetSubMesh(F);
+      SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+      smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
+      return false;
+    }
+    std::vector<int> aVec = (*anIt).second;
+    nb0d += aVec[SMDSEntity_Node];
+    nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
+    double aLen = SMESH_Algo::EdgeLength(E);
+    fullLen += aLen;
+    if(IsFirst) {
+      IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
+      IsFirst = false;
+    }
+  }
+  tmpMap.Clear();
+
+  // compute edge length
+  double ELen = 0;
+  if (_hypLengthFromEdges || !_hypLengthFromEdges && !_hypMaxElementArea) {
+    if ( nb1d > 0 )
+      ELen = fullLen / nb1d;
+  }
+  if ( _hypMaxElementArea ) {
+    double maxArea = _hypMaxElementArea->GetMaxArea();
+    ELen = sqrt(2. * maxArea/sqrt(3.0));
+  }
+  if ( ELen < Precision::Confusion() ) {
+    SMESH_subMesh *sm = aMesh.GetSubMesh(F);
+    if ( sm ) {
+      SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+      smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated.\nToo small element length",this));
+    }
+    return false;
+  }
+
+  GProp_GProps G;
+  BRepGProp::SurfaceProperties(F,G);
+  double anArea = G.Mass();
+  int nbFaces = 0;
+  if ( ELen > Precision::Confusion() )
+    nbFaces = (int) ( anArea / ( ELen*ELen*sqrt(3.) / 4 ) );
+  int nbNodes = (int) ( ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
+  std::vector<int> aVec(SMDSEntity_Last);
+  for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
+  if( IsQuadratic ) {
+    aVec[SMDSEntity_Node] = nbNodes;
+    aVec[SMDSEntity_Quad_Triangle] = nbFaces;
+  }
+  else {
+    aVec[SMDSEntity_Node] = nbNodes;
+    aVec[SMDSEntity_Triangle] = nbFaces;
+  }
+  SMESH_subMesh *sm = aMesh.GetSubMesh(F);
+  aResMap.insert(std::make_pair(sm,aVec));
+
+  return true;
+}