+ return;
+}
+
+//================================================================================
+/*!
+ * \brief Restrict local size to achieve a required _chordalError
+ */
+//================================================================================
+
+void NETGENPlugin_Mesher::SetLocalSizeForChordalError( netgen::OCCGeometry& occgeo,
+ netgen::Mesh& ngMesh)
+{
+ if ( _chordalError <= 0. )
+ return;
+
+ TopLoc_Location loc;
+ BRepLProp_SLProps surfProp( 2, 1e-6 );
+ const double sizeCoef = 0.95;
+
+ // find non-planar FACEs with non-constant curvature
+ std::vector<int> fInd;
+ for ( int i = 1; i <= occgeo.fmap.Extent(); ++i )
+ {
+ const TopoDS_Face& face = TopoDS::Face( occgeo.fmap( i ));
+ BRepAdaptor_Surface surfAd( face, false );
+ switch ( surfAd.GetType() )
+ {
+ case GeomAbs_Plane:
+ continue;
+ case GeomAbs_Cylinder:
+ case GeomAbs_Sphere:
+ case GeomAbs_Torus: // constant curvature
+ {
+ surfProp.SetSurface( surfAd );
+ surfProp.SetParameters( 0, 0 );
+ double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() ));
+ double size = elemSizeForChordalError( _chordalError, 1 / maxCurv );
+ occgeo.SetFaceMaxH( i, size * sizeCoef );
+ // limit size one edges
+ TopTools_MapOfShape edgeMap;
+ for ( TopExp_Explorer eExp( face, TopAbs_EDGE ); eExp.More(); eExp.Next() )
+ if ( edgeMap.Add( eExp.Current() ))
+ setLocalSize( TopoDS::Edge( eExp.Current() ), size, ngMesh, /*overrideMinH=*/false );
+ break;
+ }
+ default:
+ Handle(Geom_Surface) surf = BRep_Tool::Surface( face, loc );
+ if ( GeomLib_IsPlanarSurface( surf ).IsPlanar() )
+ continue;
+ fInd.push_back( i );
+ }
+ }
+ // set local size
+ if ( !fInd.empty() )
+ {
+ BRep_Builder b;
+ TopoDS_Compound allFacesComp;
+ b.MakeCompound( allFacesComp );
+ for ( size_t i = 0; i < fInd.size(); ++i )
+ b.Add( allFacesComp, occgeo.fmap( fInd[i] ));
+
+ // copy the shape to avoid spoiling its triangulation
+ TopoDS_Shape allFacesCompCopy = BRepBuilderAPI_Copy( allFacesComp );
+
+ // create triangulation with desired chordal error
+ BRepMesh_IncrementalMesh( allFacesCompCopy,
+ _chordalError,
+ /*isRelative = */Standard_False,
+ /*theAngDeflection = */ 0.5,
+ /*isInParallel = */Standard_True);
+
+ // loop on FACEs
+ for ( TopExp_Explorer fExp( allFacesCompCopy, TopAbs_FACE ); fExp.More(); fExp.Next() )
+ {
+ const TopoDS_Face& face = TopoDS::Face( fExp.Current() );
+ Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation ( face, loc );
+ if ( triangulation.IsNull() ) continue;
+
+ BRepAdaptor_Surface surf( face, false );
+ surfProp.SetSurface( surf );
+
+ gp_XY uv[3];
+ gp_XYZ p[3];
+ double size[3];
+ for ( int i = 1; i <= triangulation->NbTriangles(); ++i )
+ {
+ Standard_Integer n1,n2,n3;
+ triangulation->Triangles()(i).Get( n1,n2,n3 );
+ p [0] = triangulation->Nodes()(n1).Transformed(loc).XYZ();
+ p [1] = triangulation->Nodes()(n2).Transformed(loc).XYZ();
+ p [2] = triangulation->Nodes()(n3).Transformed(loc).XYZ();
+ uv[0] = triangulation->UVNodes()(n1).XY();
+ uv[1] = triangulation->UVNodes()(n2).XY();
+ uv[2] = triangulation->UVNodes()(n3).XY();
+ surfProp.SetParameters( uv[0].X(), uv[0].Y() );
+ if ( !surfProp.IsCurvatureDefined() )
+ break;
+
+ for ( int n = 0; n < 3; ++n ) // get size at triangle nodes
+ {
+ surfProp.SetParameters( uv[n].X(), uv[n].Y() );
+ double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() ));
+ size[n] = elemSizeForChordalError( _chordalError, 1 / maxCurv );
+ }
+ for ( int n1 = 0; n1 < 3; ++n1 ) // limit size along each triangle edge
+ {
+ int n2 = ( n1 + 1 ) % 3;
+ double minSize = size[n1], maxSize = size[n2];
+ if ( size[n1] > size[n2] )
+ minSize = size[n2], maxSize = size[n1];
+
+ if ( maxSize / minSize < 1.2 ) // netgen ignores size difference < 1.2
+ {
+ ngMesh.RestrictLocalHLine ( netgen::Point3d( p[n1].X(), p[n1].Y(), p[n1].Z() ),
+ netgen::Point3d( p[n2].X(), p[n2].Y(), p[n2].Z() ),
+ sizeCoef * minSize );
+ }
+ else
+ {
+ gp_XY uvVec( uv[n2] - uv[n1] );
+ double len = ( p[n1] - p[n2] ).Modulus();
+ int nb = int( len / minSize ) + 1;
+ for ( int j = 0; j <= nb; ++j )
+ {
+ double r = double( j ) / nb;
+ gp_XY uvj = uv[n1] + r * uvVec;
+
+ surfProp.SetParameters( uvj.X(), uvj.Y() );
+ double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() ));
+ double h = elemSizeForChordalError( _chordalError, 1 / maxCurv );
+
+ const gp_Pnt& pj = surfProp.Value();
+ netgen::Point3d ngP( pj.X(), pj.Y(), pj.Z());
+ ngMesh.RestrictLocalH( ngP, h * sizeCoef );
+ }
+ }
+ }
+ }
+ }
+ }