X-Git-Url: http://git.salome-platform.org/gitweb/?p=plugins%2Fnetgenplugin.git;a=blobdiff_plain;f=src%2FNETGENPlugin%2FNETGENPlugin_NETGEN_2D_ONLY.cxx;fp=src%2FNETGENPlugin%2FNETGENPlugin_NETGEN_2D_ONLY.cxx;h=4c4e1e6ab20cf51c17c301cee5163f2de2aa04ce;hp=16213e86aa9a684bfe210312eae33dc5754e5d2f;hb=ccecac5e15a24f8da73b0ed44cfbbbb20b15b0db;hpb=3af617de2150a1f2aace89182033d20e6fc0b237 diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx index 16213e8..4c4e1e6 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx @@ -22,8 +22,6 @@ // Project : SALOME // #include "NETGENPlugin_NETGEN_2D_ONLY.hxx" - -#include "NETGENPlugin_Mesher.hxx" #include "NETGENPlugin_Hypothesis_2D.hxx" #include @@ -46,6 +44,8 @@ #include +#include + #include #include #include @@ -61,6 +61,23 @@ namespace nglib { #endif #include #include + +#include "SMESH_Octree.hxx" + +//// Used for node projection in curve +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + //#include namespace netgen { NETGENPLUGIN_DLL_HEADER @@ -74,6 +91,264 @@ using namespace std; using namespace netgen; using namespace nglib; +namespace // copied class from StdMesher_Importe_1D +{ + /*! + * \brief Compute point position on a curve. Use octree to fast reject far points + */ + class CurveProjector : public SMESH_Octree + { + public: + CurveProjector( const TopoDS_Edge& edge, double enlarge ); + + bool IsOnCurve( const gp_XYZ& point, double & distance2, double & u ); + + bool IsOut( const gp_XYZ& point ) const { return getBox()->IsOut( point ); } + + protected: + CurveProjector() {} + SMESH_Octree* newChild() const { return new CurveProjector; } + void buildChildrenData(); + Bnd_B3d* buildRootBox(); + + private: + struct CurveSegment : public Bnd_B3d + { + double _chord, _chord2, _length2; + gp_Pnt _pFirst, _pLast; + gp_Lin _line; + Handle(Geom_Curve) _curve; + + CurveSegment() {} + void Init( const gp_Pnt& pf, const gp_Pnt& pl, + double uf, double ul, double tol, Handle(Geom_Curve)& curve ); + bool IsOn( const gp_XYZ& point, double & distance2, double & u ); + bool IsInContact( const Bnd_B3d& bb ); + }; + std::vector< CurveSegment > _segments; + }; + + //=============================================================================== + /*! + * \brief Create an octree of curve segments + */ + //================================================================================ + + CurveProjector::CurveProjector( const TopoDS_Edge& edge, double enlarge ) + :SMESH_Octree( 0 ) + { + double f,l; + Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, f, l ); + double curDeflect = 0.3; // Curvature deflection + double angDeflect = 1e+100; // Angular deflection - don't control chordal error + GCPnts_TangentialDeflection div( BRepAdaptor_Curve( edge ), angDeflect, curDeflect ); + _segments.resize( div.NbPoints() - 1 ); + for ( int i = 1; i < div.NbPoints(); ++i ) + try { + _segments[ i - 1 ].Init( div.Value( i ), div.Value( i+1 ), + div.Parameter( i ), div.Parameter( i+1 ), + enlarge, curve ); + } + catch ( Standard_Failure ) { + _segments.resize( _segments.size() - 1 ); + --i; + } + if ( _segments.size() < 3 ) + myIsLeaf = true; + + compute(); + + if ( _segments.size() == 1 ) + myBox->Enlarge( enlarge ); + } + + //================================================================================ + /*! + * \brief Return the maximal box + */ + //================================================================================ + + Bnd_B3d* CurveProjector::buildRootBox() + { + Bnd_B3d* box = new Bnd_B3d; + for ( size_t i = 0; i < _segments.size(); ++i ) + box->Add( _segments[i] ); + return box; + } + + //================================================================================ + /*! + * \brief Redistribute segments among children + */ + //================================================================================ + + void CurveProjector::buildChildrenData() + { + bool allIn = true; + for ( size_t i = 0; i < _segments.size(); ++i ) + { + for (int j = 0; j < 8; j++) + { + if ( _segments[i].IsInContact( *myChildren[j]->getBox() )) + ((CurveProjector*)myChildren[j])->_segments.push_back( _segments[i]); + else + allIn = false; + } + } + if ( allIn && _segments.size() < 3 ) + { + myIsLeaf = true; + for (int j = 0; j < 8; j++) + static_cast( myChildren[j])->myIsLeaf = true; + } + else + { + SMESHUtils::FreeVector( _segments ); // = _segments.clear() + free memory + + for (int j = 0; j < 8; j++) + { + CurveProjector* child = static_cast( myChildren[j]); + if ( child->_segments.size() < 3 ) + child->myIsLeaf = true; + } + } + } + + //================================================================================ + /*! + * \brief Return true if a point is close to the curve + * \param [in] point - the point + * \param [out] distance2 - distance to the curve + * \param [out] u - parameter on the curve + * \return bool - is the point is close to the curve + */ + //================================================================================ + + bool CurveProjector::IsOnCurve( const gp_XYZ& point, double & distance2, double & u ) + { + if ( getBox()->IsOut( point )) + return false; + + bool ok = false; + double dist2, param; + distance2 = Precision::Infinite(); + + if ( isLeaf() ) + { + for ( size_t i = 0; i < _segments.size(); ++i ) + if ( !_segments[i].IsOut( point ) && + _segments[i].IsOn( point, dist2, param ) && + dist2 < distance2 ) + { + distance2 = dist2; + u = param; + ok = true; + } + return ok; + } + else + { + for (int i = 0; i < 8; i++) + if (((CurveProjector*) myChildren[i])->IsOnCurve( point, dist2, param ) && + dist2 < distance2 ) + { + distance2 = dist2; + u = param; + ok = true; + } + } + return ok; + } + + //================================================================================ + /*! + * \brief Initialize + */ + //================================================================================ + + void CurveProjector::CurveSegment::Init(const gp_Pnt& pf, + const gp_Pnt& pl, + const double uf, + const double ul, + const double tol, + Handle(Geom_Curve)& curve ) + { + _pFirst = pf; + _pLast = pl; + _curve = curve; + _length2 = pf.SquareDistance( pl ); + _line.SetLocation( pf ); + _line.SetDirection( gp_Vec( pf, pl )); + _chord2 = Max( _line. SquareDistance( curve->Value( uf + 0.25 * ( ul - uf ))), + Max( _line.SquareDistance( curve->Value( uf + 0.5 * ( ul - uf ))), + _line.SquareDistance( curve->Value( uf + 0.75 * ( ul - uf ))))); + _chord2 *= ( 1.05 * 1.05 ); // +5% + _chord2 = Max( tol, _chord2 ); + _chord = Sqrt( _chord2 ); + + Bnd_Box bb; + BndLib_Add3dCurve::Add( GeomAdaptor_Curve( curve, uf, ul ), tol, bb ); + Add( bb.CornerMin() ); + Add( bb.CornerMax() ); + } + + //================================================================================ + /*! + * \brief Return true if a point is close to the curve segment + * \param [in] point - the point + * \param [out] distance2 - distance to the curve + * \param [out] u - parameter on the curve + * \return bool - is the point is close to the curve segment + */ + //================================================================================ + + bool CurveProjector::CurveSegment::IsOn( const gp_XYZ& point, double & distance2, double & u ) + { + distance2 = _line.SquareDistance( point ); + if ( distance2 > _chord2 ) + return false; + + // check if the point projection falls into the segment range + { + gp_Vec edge( _pFirst, _pLast ); + gp_Vec n1p ( _pFirst, point ); + u = ( edge * n1p ) / _length2; // param [0,1] on the edge + if ( u < 0. ) + { + if ( _pFirst.SquareDistance( point ) > _chord2 ) + return false; + } + else if ( u > 1. ) + { + if ( _pLast.SquareDistance( point ) > _chord2 ) + return false; + } + } + gp_Pnt proj; + distance2 = ShapeAnalysis_Curve().Project( _curve, point, Precision::Confusion(), + proj, u, false ); + distance2 *= distance2; + return true; + } + + //================================================================================ + /*! + * \brief Check if the segment is in contact with a box + */ + //================================================================================ + + bool CurveProjector::CurveSegment::IsInContact( const Bnd_B3d& bb ) + { + if ( bb.IsOut( _line.Position(), /*isRay=*/true, _chord )) + return false; + + gp_Ax1 axRev = _line.Position().Reversed(); + axRev.SetLocation( _pLast ); + return !bb.IsOut( axRev, /*isRay=*/true, _chord ); + } + +} + //============================================================================= /*! * @@ -217,43 +492,216 @@ bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh& aMesh, // } // } -//============================================================================= -/*! - *Here we are going to use the NETGEN mesher - */ -//============================================================================= -bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape) -{ - netgen::multithread.terminate = 0; - //netgen::multithread.task = "Surface meshing"; - SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); +/** + * @brief MapSegmentsToEdges. + * @remark To feed 1D segments not associated to any geometry we need: + * 1) For each face, check all segments that are in the face (use ShapeAnalysis_Surface class) + * 2) Check to which edge the segment below to, use the copied [from StdMesher_Importe_1D] CurveProjector class + * 3) Create new netgen segments with the (u,v) parameters obtained from the ShapeAnalysis_Surface projector + * 4) also define the 'param' value of the nodes relative to the edges obtained from CurveProjector + * 5) Add the new netgen segments IN ORDER into the netgen::mesh data structure to form a closed chain + * 6) Beware with the occ::edge orientation + * + * @param aMesh Mesh file (containing 1D elements) + * @param aShape Shape file (BREP or STEP format) + * @param ngLib netgenlib library wrapper + * @param nodeVec vector of nodes used internally to feed smesh aMesh after computation + * @param premeshedNodes map of prmeshed nodes and the smesh nodeID associate to it + * @param newNetgenCoordinates map of 3D coordinate of new points created by netgen + * @param newNetgenElements map of triangular or quadrangular elements ID and the nodes defining the 2D element + * @return false + */ +bool NETGENPlugin_NETGEN_2D_ONLY::MapSegmentsToEdges(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape, NETGENPlugin_NetgenLibWrapper &ngLib, + vector< const SMDS_MeshNode* >& nodeVec, std::map& premeshedNodes, + std::map>& newNetgenCoordinates, std::map>& newNetgenElements ) +{ + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); SMESH_MesherHelper helper(aMesh); helper.SetElementsOnShape( true ); + const int numberOfPremeshedNodes = aMesh.NbNodes(); + TopTools_IndexedMapOfShape faces; + TopExp::MapShapes( aShape, TopAbs_FACE, faces ); - NETGENPlugin_NetgenLibWrapper ngLib; - ngLib._isComputeOk = false; + for ( int i = 1; i <= faces.Size(); ++i ) + { + int numOfEdges = 0; + int totalEdgeLenght = 0; - netgen::Mesh ngMeshNoLocSize; - netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh, & ngMeshNoLocSize }; - netgen::OCCGeometry occgeoComm; + TopoDS_Face face = TopoDS::Face(faces( i ) ); + Bnd_Box FaceBox; + GEOMUtils::PreciseBoundingBox( face, FaceBox ); + if ( face.Orientation() != TopAbs_FORWARD && face.Orientation() != TopAbs_REVERSED ) + face.Orientation( TopAbs_FORWARD ); - // min / max sizes are set as follows: - // if ( _hypParameters ) - // min and max are defined by the user - // else if ( _hypLengthFromEdges ) - // min = aMesher.GetDefaultMinSize() - // max = average segment len of a FACE - // else if ( _hypMaxElementArea ) - // min = aMesher.GetDefaultMinSize() - // max = f( _hypMaxElementArea ) - // else - // min = aMesher.GetDefaultMinSize() - // max = max segment len of a FACE + // Set the occgeom to be meshed and some ngMesh parameteres + netgen::OCCGeometry occgeom; + netgen::Mesh * ngMesh = (netgen::Mesh*) ngLib._ngMesh; + ngMesh->DeleteMesh(); - NETGENPlugin_Mesher aMesher( &aMesh, aShape, /*isVolume=*/false); + occgeom.shape = face; + occgeom.fmap.Add( face ); + occgeom.CalcBoundingBox(); + occgeom.facemeshstatus.SetSize(1); + occgeom.facemeshstatus = 0; + occgeom.face_maxh_modified.SetSize(1); + occgeom.face_maxh_modified = 0; + occgeom.face_maxh.SetSize(1); + occgeom.face_maxh = netgen::mparam.maxh; + + // Set the face descriptor + const int solidID = 0, faceID = 1; /*always 1 because faces are meshed one by one*/ + if ( ngMesh->GetNFD() < 1 ) + ngMesh->AddFaceDescriptor( netgen::FaceDescriptor( faceID, solidID, solidID, 0 )); + + ngMesh->SetGlobalH ( netgen::mparam.maxh ); + ngMesh->SetMinimalH( netgen::mparam.minh ); + Box<3> bb = occgeom.GetBoundingBox(); + bb.Increase (bb.Diam()/10); + ngMesh->SetLocalH (bb.PMin(), bb.PMax(), netgen::mparam.grading); + // end set the occgeom to be meshed and some ngMesh parameteres + + Handle(ShapeAnalysis_Surface) sprojector = new ShapeAnalysis_Surface( BRep_Tool::Surface( face )); + double tol = BRep_Tool::MaxTolerance( face, TopAbs_FACE ); + gp_Pnt surfPnt(0,0,0); + + map node2ngID; + map ng2smesh; + // Figure out which edge is this onde in! + TopTools_IndexedMapOfShape edges; + TopExp::MapShapes( face, TopAbs_EDGE, edges ); + TopoDS_Edge meshingEdge; + // Check wich nodes are in this face! + SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge); + std::unique_ptr myCurveProjector; + while ( iteratorElem->more() ) // loop on elements on a geom face + { + // check mesh face + const SMDS_MeshElement* elem = iteratorElem->next(); + const SMDS_MeshNode* node0 = elem->GetNode( 0 ); + const SMDS_MeshNode* node1 = elem->GetNode( 1 ); + SMESH_NodeXYZ nXYZ0( node0 ); + SMESH_NodeXYZ nXYZ1( node1 ); + double segmentLength = ( nXYZ0 - nXYZ1 ).Modulus(); + nXYZ0 += nXYZ1; + nXYZ0 /= 2.0; + if( FaceBox.IsOut( nXYZ0 ) ) + continue; + + gp_XY uv = sprojector->ValueOfUV( nXYZ0, tol ).XY(); + surfPnt = sprojector->Value( uv ); + double dist = surfPnt.Distance( nXYZ0 ); + + // used for the curve projector of edges + double geomTol = Precision::Confusion(); + if ( dist < tol /*element is on face*/ ) + { + numOfEdges++; + totalEdgeLenght += segmentLength; + int occEdgeIdFound = -1; + for ( int edgeId = 1; edgeId <= edges.Size(); ++edgeId ) /*find in which edge the node is placed*/ + { + meshingEdge = TopoDS::Edge(edges( edgeId )); + myCurveProjector = std::unique_ptr( new CurveProjector(meshingEdge, geomTol) ); + if ( myCurveProjector->IsOut( nXYZ0 ) /*keep searching*/) + continue; + else + { + occEdgeIdFound = edgeId; + break; + } + } + + netgen::Segment seg; + for ( size_t nodeId = 0; nodeId < 2; nodeId++) + { + const SMDS_MeshNode* node = elem->GetNode( nodeId ); + int ngId = ngMesh->GetNP() + 1; + ngId = node2ngID.insert( make_pair( node, ngId )).first->second; + if ( ngId > ngMesh->GetNP() /* mean it is a new node to be add to the mesh*/) + { + // Restric size of mesh based on the edge dimension + { + SMESH_NodeXYZ nXYZ( node ); + netgen::Point3d pi(nXYZ.X(), nXYZ.Y(), nXYZ.Z()); + ngMesh->RestrictLocalH( pi, segmentLength ); + } + + netgen::MeshPoint mp( netgen::Point<3> (node->X(), node->Y(), node->Z()) ); + ngMesh->AddPoint ( mp, 1, netgen::EDGEPOINT ); + ng2smesh.insert( make_pair(ngId, node) ); + premeshedNodes.insert( make_pair( (int)node->GetID(), node ) ); + } + seg[nodeId] = ngId; // ng node id + SMESH_NodeXYZ nXYZ( node ); + gp_XY uv = sprojector->ValueOfUV( nXYZ, tol ).XY(); + + // Compute the param (relative distance) of the node in relation the edge + // fundamental for netgen working properly + double dist2, param; + if ( myCurveProjector->IsOnCurve( nXYZ, dist2, param ) ) + { + seg.epgeominfo[ nodeId ].dist = param; + seg.epgeominfo[ nodeId ].u = uv.X(); + seg.epgeominfo[ nodeId ].v = uv.Y(); + seg.epgeominfo[ nodeId ].edgenr = occEdgeIdFound; + } + } + + if ( meshingEdge.Orientation() != TopAbs_FORWARD ) + { + swap(seg[0], seg[1]); + swap(seg.epgeominfo[0], seg.epgeominfo[1] ); + } + + seg.edgenr = ngMesh->GetNSeg() + 1; // netgen segment id + seg.si = faceID; + ngMesh->AddSegment( seg ); + } // end if edge is on the face + } // end iteration on elements + + // set parameters from _hypLengthFromEdges if needed + if ( !_hypParameters && _hypLengthFromEdges && numOfEdges ) + { + netgen::mparam.maxh = totalEdgeLenght / numOfEdges; + if ( netgen::mparam.maxh < DBL_MIN ) + netgen::mparam.maxh = occgeom.GetBoundingBox().Diam(); + + occgeom.face_maxh = netgen::mparam.maxh; + ngMesh->SetGlobalH ( netgen::mparam.maxh ); + } + // end set parameters + + ngMesh->CalcSurfacesOfNode(); + const int startWith = MESHCONST_MESHSURFACE; + const int endWith = MESHCONST_OPTSURFACE; + int err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh); + + // Ng_Mesh * ngMeshptr = (Ng_Mesh*) ngLib._ngMesh; + // int NetgenNodes = Ng_GetNP(ngMeshptr); + // int NetgenSeg_2D = Ng_GetNSeg_2D(ngMeshptr); + // int NetgenFaces = Ng_GetNSE(ngMeshptr); + // std::cout << "\n"; + // std::cout << "Number of nodes, seg, faces, vols: " << NetgenNodes << ", " << NetgenSeg_2D << ", " << NetgenFaces << "\n"; + // std::cout << "err from netgen computation: " << err << "\n"; + + // ---------------------------------------------------- + // Fill the SMESHDS with the generated nodes and faces + // ---------------------------------------------------- + nodeVec.clear(); + FillNodesAndElements( aMesh, helper, ngMesh, nodeVec, ng2smesh, newNetgenCoordinates, newNetgenElements, numberOfPremeshedNodes ); + } // Face iteration + + return false; +} + +std::tuple NETGENPlugin_NETGEN_2D_ONLY::SetParameteres( SMESH_Mesh& aMesh, const TopoDS_Shape& aShape, + NETGENPlugin_Mesher& aMesher, netgen::Mesh * ngMeshes, + netgen::OCCGeometry& occgeoComm, bool isSubMeshSupported ) +{ + SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); + aMesher.SetParameters( _hypParameters ); // _hypParameters -> netgen::mparam const bool toOptimize = _hypParameters ? _hypParameters->GetOptimize() : true; if ( _hypMaxElementArea ) @@ -269,7 +717,6 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, if ( isCommonLocalSize ) // compute common local size in ngMeshes[0] { - //list< SMESH_subMesh* > meshedSM[4]; --> all sub-shapes are added to occgeoComm aMesher.PrepareOCCgeometry( occgeoComm, aShape, aMesh );//, meshedSM ); // local size set at MESHCONST_ANALYSE step depends on @@ -291,44 +738,268 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, occgeoComm.face_maxh = netgen::mparam.maxh; #ifdef NETGEN_V6 netgen::OCCParameters occparam; - netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0], netgen::mparam, occparam ); + netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes, netgen::mparam, occparam ); #else - netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0] ); + netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes ); #endif occgeoComm.emap.Clear(); occgeoComm.vmap.Clear(); - // set local size according to size of existing segments - TopTools_IndexedMapOfShape edgeMap; - TopExp::MapShapes( aMesh.GetShapeToMesh(), TopAbs_EDGE, edgeMap ); - for ( int iE = 1; iE <= edgeMap.Extent(); ++iE ) + if ( isSubMeshSupported ) { - const TopoDS_Shape& edge = edgeMap( iE ); - if ( SMESH_Algo::isDegenerated( TopoDS::Edge( edge ))) - continue; - SMESHDS_SubMesh* smDS = meshDS->MeshElements( edge ); - if ( !smDS ) continue; - SMDS_ElemIteratorPtr segIt = smDS->GetElements(); - while ( segIt->more() ) + // set local size according to size of existing segments + TopTools_IndexedMapOfShape edgeMap; + TopExp::MapShapes( aMesh.GetShapeToMesh(), TopAbs_EDGE, edgeMap ); + for ( int iE = 1; iE <= edgeMap.Extent(); ++iE ) { - const SMDS_MeshElement* seg = segIt->next(); - SMESH_TNodeXYZ n1 = seg->GetNode(0); - SMESH_TNodeXYZ n2 = seg->GetNode(1); - gp_XYZ p = 0.5 * ( n1 + n2 ); + const TopoDS_Shape& edge = edgeMap( iE ); + if ( SMESH_Algo::isDegenerated( TopoDS::Edge( edge ))) + continue; + SMESHDS_SubMesh* smDS = meshDS->MeshElements( edge ); + if ( !smDS ) continue; + SMDS_ElemIteratorPtr segIt = smDS->GetElements(); + while ( segIt->more() ) + { + const SMDS_MeshElement* seg = segIt->next(); + SMESH_TNodeXYZ n1 = seg->GetNode(0); + SMESH_TNodeXYZ n2 = seg->GetNode(1); + gp_XYZ p = 0.5 * ( n1 + n2 ); + netgen::Point3d pi(p.X(), p.Y(), p.Z()); + ngMeshes->RestrictLocalH( pi, factor * ( n1 - n2 ).Modulus() ); + } + } + } + else + { + SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge); + while ( iteratorElem->more() ) // loop on elements on a geom face + { + const SMDS_MeshElement* elem = iteratorElem->next(); + const SMDS_MeshNode* node0 = elem->GetNode( 0 ); + const SMDS_MeshNode* node1 = elem->GetNode( 1 ); + SMESH_NodeXYZ nXYZ0( node0 ); + SMESH_NodeXYZ nXYZ1( node1 ); + double segmentLength = ( nXYZ0 - nXYZ1 ).Modulus(); + gp_XYZ p = 0.5 * ( nXYZ0 + nXYZ1 ); netgen::Point3d pi(p.X(), p.Y(), p.Z()); - ngMeshes[0]->RestrictLocalH( pi, factor * ( n1 - n2 ).Modulus() ); + ngMeshes->RestrictLocalH( pi, factor * segmentLength ); } } // set local size defined on shapes - aMesher.SetLocalSize( occgeoComm, *ngMeshes[0] ); - aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes[0] ); + aMesher.SetLocalSize( occgeoComm, *ngMeshes ); + aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes ); try { - ngMeshes[0]->LoadLocalMeshSize( mparam.meshsizefilename ); + ngMeshes->LoadLocalMeshSize( mparam.meshsizefilename ); } catch (NgException & ex) { - return error( COMPERR_BAD_PARMETERS, ex.What() ); + throw error( COMPERR_BAD_PARMETERS, ex.What() ); + } + } + + return std::make_tuple( isCommonLocalSize, isDefaultHyp ); +} + +bool NETGENPlugin_NETGEN_2D_ONLY::ComputeMaxhOfFace( TopoDS_Face& Face, NETGENPlugin_Mesher& aMesher, TSideVector& wires, + netgen::OCCGeometry& occgeoComm, bool isDefaultHyp, bool isCommonLocalSize ) +{ + size_t nbWires = wires.size(); + if ( !_hypParameters ) + { + double edgeLength = 0; + if (_hypLengthFromEdges ) + { + // compute edgeLength as an average segment length + smIdType nbSegments = 0; + for ( size_t iW = 0; iW < nbWires; ++iW ) + { + edgeLength += wires[ iW ]->Length(); + nbSegments += wires[ iW ]->NbSegments(); + } + if ( nbSegments ) + edgeLength /= double( nbSegments ); + netgen::mparam.maxh = edgeLength; + } + else if ( isDefaultHyp ) + { + // set edgeLength by a longest segment + double maxSeg2 = 0; + for ( size_t iW = 0; iW < nbWires; ++iW ) + { + const UVPtStructVec& points = wires[ iW ]->GetUVPtStruct(); + if ( points.empty() ) + return error( COMPERR_BAD_INPUT_MESH ); + gp_Pnt pPrev = SMESH_TNodeXYZ( points[0].node ); + for ( size_t i = 1; i < points.size(); ++i ) + { + gp_Pnt p = SMESH_TNodeXYZ( points[i].node ); + maxSeg2 = Max( maxSeg2, p.SquareDistance( pPrev )); + pPrev = p; + } + } + edgeLength = sqrt( maxSeg2 ) * 1.05; + netgen::mparam.maxh = edgeLength; + } + if ( netgen::mparam.maxh < DBL_MIN ) + netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam(); + + if ( !isCommonLocalSize ) + { + netgen::mparam.minh = aMesher.GetDefaultMinSize( Face, netgen::mparam.maxh ); + } + } + return true; +} + +void NETGENPlugin_NETGEN_2D_ONLY::FillNodesAndElements( SMESH_Mesh& aMesh, SMESH_MesherHelper& helper, netgen::Mesh * ngMesh, + vector< const SMDS_MeshNode* >& nodeVec, map& ng2smesh, + std::map>& newNetgenCoordinates, + std::map>& newNetgenElements, const int numberOfPremeshedNodes ) +{ + SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); + + int nbNodes = ngMesh->GetNP(); + int nbFaces = ngMesh->GetNSE(); + nodeVec.resize( nbNodes + 1, 0 ); + // map to index local node numeration to global numeration used by Remote mesher to write the final global resulting mesh + std::map local2globalMap; + + smIdType myNewCoordinateCounter = newNetgenCoordinates.size() > 0 ? newNetgenCoordinates.rbegin()->first + 1: numberOfPremeshedNodes+1; + int myNewFaceCounter = newNetgenElements.size() > 0 ? newNetgenElements.rbegin()->first + 1 : 1; + // add nodes + for ( int ngID = 1; ngID <= nbNodes; ++ngID ) + { + const MeshPoint& ngPoint = ngMesh->Point( ngID ); + // Check if ngPoint is not already present because was in the premeshed mesh boundary + if ( ng2smesh.count( ngID ) == 0 ) + { + std::vector netgenCoordinates = {ngPoint(0), ngPoint(1), ngPoint(2)}; + SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2)); + nodeVec[ ngID ] = node; + newNetgenCoordinates.insert( make_pair( myNewCoordinateCounter, std::move(netgenCoordinates)) ); + local2globalMap.insert( std::make_pair( node->GetID(), myNewCoordinateCounter ) ); + myNewCoordinateCounter++; + } + else + { + nodeVec[ ngID ] = ng2smesh[ ngID ]; + local2globalMap.insert( std::make_pair( nodeVec[ ngID ]->GetID(), nodeVec[ ngID ]->GetID() ) ); + } + } + // create faces + int i,j; + vector nodes; + for ( i = 1; i <= nbFaces ; ++i ) + { + const Element2d& elem = ngMesh->SurfaceElement(i); + nodes.resize( elem.GetNP() ); + for (j=1; j <= elem.GetNP(); ++j) + { + int pind = elem.PNum(j); + if ( pind < 1 ) + break; + nodes[ j-1 ] = nodeVec[ pind ]; + } + if ( j > elem.GetNP() ) + { + std::vector netgenCoordinates = { local2globalMap[nodes[0]->GetID()], local2globalMap[nodes[1]->GetID()], local2globalMap[nodes[2]->GetID()] }; + newNetgenElements.insert( std::make_pair( myNewFaceCounter, std::move( netgenCoordinates ) ) ); + helper.AddFace(nodes[0],nodes[1],nodes[2]); + myNewFaceCounter++; } } +} + +void NETGENPlugin_NETGEN_2D_ONLY::FillNodesAndElements( SMESH_Mesh& aMesh, SMESH_MesherHelper& helper, netgen::Mesh * ngMesh, vector< const SMDS_MeshNode* >& nodeVec, int faceId ) +{ + SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); + + int nbNodes = ngMesh->GetNP(); + int nbFaces = ngMesh->GetNSE(); + + int nbInputNodes = (int) nodeVec.size()-1; + nodeVec.resize( nbNodes+1, 0 ); + + // add nodes + for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID ) + { + const MeshPoint& ngPoint = ngMesh->Point( ngID ); + SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2)); + nodeVec[ ngID ] = node; + } + + // create faces + int i,j; + vector nodes; + for ( i = 1; i <= nbFaces ; ++i ) + { + const Element2d& elem = ngMesh->SurfaceElement(i); + nodes.resize( elem.GetNP() ); + for (j=1; j <= elem.GetNP(); ++j) + { + int pind = elem.PNum(j); + if ( pind < 1 ) + break; + nodes[ j-1 ] = nodeVec[ pind ]; + if ( nodes[ j-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) + { + const PointGeomInfo& pgi = elem.GeomInfoPi(j); + meshDS->SetNodeOnFace( nodes[ j-1 ], faceId, pgi.u, pgi.v); + } + } + if ( j > elem.GetNP() ) + { + if ( elem.GetType() == TRIG ) + helper.AddFace(nodes[0],nodes[1],nodes[2]); + else + helper.AddFace(nodes[0],nodes[1],nodes[2],nodes[3]); + } + } +} + + +//============================================================================= +/*! + *Here we are going to use the NETGEN mesher + */ +//============================================================================= + +bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape) +{ + netgen::multithread.terminate = 0; + //netgen::multithread.task = "Surface meshing"; + + SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); + SMESH_MesherHelper helper(aMesh); + helper.SetElementsOnShape( true ); + + NETGENPlugin_NetgenLibWrapper ngLib; + ngLib._isComputeOk = false; + + netgen::Mesh ngMeshNoLocSize; + netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh, & ngMeshNoLocSize }; + + netgen::OCCGeometry occgeoComm; + + // min / max sizes are set as follows: + // if ( _hypParameters ) + // min and max are defined by the user + // else if ( _hypLengthFromEdges ) + // min = aMesher.GetDefaultMinSize() + // max = average segment len of a FACE + // else if ( _hypMaxElementArea ) + // min = aMesher.GetDefaultMinSize() + // max = f( _hypMaxElementArea ) + // else + // min = aMesher.GetDefaultMinSize() + // max = max segment len of a FACE + + NETGENPlugin_Mesher aMesher( &aMesh, aShape, /*isVolume=*/false); + auto options = SetParameteres( aMesh, aShape, aMesher, ngMeshes[0], occgeoComm ); + const bool isCommonLocalSize = std::get<0>( options ); + const bool isDefaultHyp = std::get<1>( options ); + const bool toOptimize = _hypParameters ? _hypParameters->GetOptimize() : true; + netgen::mparam.uselocalh = toOptimize; // restore as it is used at surface optimization // ================== @@ -382,51 +1053,10 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, // compute maxh of a FACE // ---------------------- - if ( !_hypParameters ) - { - double edgeLength = 0; - if (_hypLengthFromEdges ) - { - // compute edgeLength as an average segment length - smIdType nbSegments = 0; - for ( size_t iW = 0; iW < nbWires; ++iW ) - { - edgeLength += wires[ iW ]->Length(); - nbSegments += wires[ iW ]->NbSegments(); - } - if ( nbSegments ) - edgeLength /= double( nbSegments ); - netgen::mparam.maxh = edgeLength; - } - else if ( isDefaultHyp ) - { - // set edgeLength by a longest segment - double maxSeg2 = 0; - for ( size_t iW = 0; iW < nbWires; ++iW ) - { - const UVPtStructVec& points = wires[ iW ]->GetUVPtStruct(); - if ( points.empty() ) - return error( COMPERR_BAD_INPUT_MESH ); - gp_Pnt pPrev = SMESH_TNodeXYZ( points[0].node ); - for ( size_t i = 1; i < points.size(); ++i ) - { - gp_Pnt p = SMESH_TNodeXYZ( points[i].node ); - maxSeg2 = Max( maxSeg2, p.SquareDistance( pPrev )); - pPrev = p; - } - } - edgeLength = sqrt( maxSeg2 ) * 1.05; - netgen::mparam.maxh = edgeLength; - } - if ( netgen::mparam.maxh < DBL_MIN ) - netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam(); - - if ( !isCommonLocalSize ) - { - netgen::mparam.minh = aMesher.GetDefaultMinSize( F, netgen::mparam.maxh ); - } - } - + bool setMaxh = ComputeMaxhOfFace( F, aMesher, wires, occgeoComm, isDefaultHyp, isCommonLocalSize ); + if (!setMaxh) + return setMaxh; + // prepare occgeom netgen::OCCGeometry occgeom; occgeom.shape = F; @@ -546,52 +1176,10 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, } } - // ---------------------------------------------------- // Fill the SMESHDS with the generated nodes and faces // ---------------------------------------------------- - - int nbNodes = ngMesh->GetNP(); - int nbFaces = ngMesh->GetNSE(); - - int nbInputNodes = (int) nodeVec.size()-1; - nodeVec.resize( nbNodes+1, 0 ); - - // add nodes - for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID ) - { - const MeshPoint& ngPoint = ngMesh->Point( ngID ); - SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2)); - nodeVec[ ngID ] = node; - } - - // create faces - int i,j; - vector nodes; - for ( i = 1; i <= nbFaces ; ++i ) - { - const Element2d& elem = ngMesh->SurfaceElement(i); - nodes.resize( elem.GetNP() ); - for (j=1; j <= elem.GetNP(); ++j) - { - int pind = elem.PNum(j); - if ( pind < 1 ) - break; - nodes[ j-1 ] = nodeVec[ pind ]; - if ( nodes[ j-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) - { - const PointGeomInfo& pgi = elem.GeomInfoPi(j); - meshDS->SetNodeOnFace( nodes[ j-1 ], faceID, pgi.u, pgi.v); - } - } - if ( j > elem.GetNP() ) - { - if ( elem.GetType() == TRIG ) - helper.AddFace(nodes[0],nodes[1],nodes[2]); - else - helper.AddFace(nodes[0],nodes[1],nodes[2],nodes[3]); - } - } + FillNodesAndElements( aMesh, helper, ngMesh, nodeVec, faceID ); break; } // two attempts