From: cconopoima Date: Tue, 30 Apr 2024 11:20:03 +0000 (+0100) Subject: Using RestrictLocalH netgen function to simulate feed of 1D elements for 2D_Only... X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=e7e012257be9b644c826b144578649b77e5e16ab;p=plugins%2Fnetgenplugin.git Using RestrictLocalH netgen function to simulate feed of 1D elements for 2D_Only and submeshing of 2D! --- diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx index efa5e70..7a49e1f 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx @@ -52,6 +52,8 @@ #include +#include + #include #include #include @@ -957,7 +959,8 @@ void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry& occgeo, updateTriangulation( shape ); Bnd_Box bb; - BRepBndLib::Add (shape, bb); + // BRepBndLib::Add (shape, bb); + GEOMUtils::PreciseBoundingBox(shape,bb); double x1,y1,z1,x2,y2,z2; bb.Get (x1,y1,z1,x2,y2,z2); netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1); @@ -1026,6 +1029,9 @@ void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry& occgeo, occgeo.face_maxh_modified = 0; occgeo.face_maxh.SetSize(totNbFaces); occgeo.face_maxh = netgen::mparam.maxh; + + // + occgeo.BuildFMapFromPrefilled(); } //================================================================================ @@ -1153,10 +1159,11 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, const SMESHDS_SubMesh * smDS = sm->GetSubMeshDS(); if ( !smDS ) continue; + int countEdgeId = 0; switch ( sm->GetSubShape().ShapeType() ) { case TopAbs_EDGE: { // EDGE - // ---------------------- + TopoDS_Edge geomEdge = TopoDS::Edge( sm->GetSubShape() ); if ( geomEdge.Orientation() >= TopAbs_INTERNAL ) geomEdge.Orientation( TopAbs_FORWARD ); // issue 0020676 @@ -1180,7 +1187,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, // get all meshed EDGEs of the FACE connected to geomEdge (issue 0021140) helper.SetSubShape( face ); list< TopoDS_Edge > edges = getConnectedEdges( geomEdge, face, computedSM, helper, - visitedEdgeSM2Faces ); + visitedEdgeSM2Faces ); if ( edges.empty() ) continue; // wrong ancestor? @@ -1203,119 +1210,191 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, if ( points.empty() ) return false; // invalid node params? smIdType i, nbSeg = fSide.NbSegments(); + for ( i = 0; i < nbSeg; ++i ) + { + const SMDS_MeshNode * n0 = points[ i ].node; + const SMDS_MeshNode * n1 = points[ i+1 ].node; + double size = SMESH_TNodeXYZ( n0 ).Distance( n1 ); + netgen::Point3d p0(n0->X(), n0->Y(), n0->Z()); + netgen::Point3d p1(n1->X(), n1->Y(), n1->Z()); + ngMesh.RestrictLocalHLine( p0, p1, size, true ); + } - // remember EDGEs of fSide to treat only once - for ( int iE = 0; iE < fSide.NbEdges(); ++iE ) - visitedEdgeSM2Faces[ helper.GetMesh()->GetSubMesh( fSide.Edge(iE )) ].insert(faceSMDSId); + } +/////////////////////////////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////////////////////// + { + /* + // ---------------------- + TopoDS_Edge geomEdge = TopoDS::Edge( sm->GetSubShape() ); + if ( geomEdge.Orientation() >= TopAbs_INTERNAL ) + geomEdge.Orientation( TopAbs_FORWARD ); // issue 0020676 + + // Add ng segments for each not meshed FACE the EDGE bounds + PShapeIteratorPtr fIt = helper.GetAncestors( geomEdge, *sm->GetFather(), TopAbs_FACE ); + while ( const TopoDS_Shape * anc = fIt->next() ) + { + faceNgID = occgeom.fmap.FindIndex( *anc ); + if ( faceNgID < 1 ) + continue; // meshed face + + int faceSMDSId = meshDS->ShapeToIndex( *anc ); + if ( visitedEdgeSM2Faces[ sm ].count( faceSMDSId )) + continue; // already treated EDGE + + TopoDS_Face face = TopoDS::Face( occgeom.fmap( faceNgID )); + if ( face.Orientation() >= TopAbs_INTERNAL ) + face.Orientation( TopAbs_FORWARD ); // issue 0020676 + + // get all meshed EDGEs of the FACE connected to geomEdge (issue 0021140) + helper.SetSubShape( face ); + list< TopoDS_Edge > edges = getConnectedEdges( geomEdge, face, computedSM, helper, + visitedEdgeSM2Faces ); + if ( edges.empty() ) + continue; // wrong ancestor? + + // find out orientation of within + TopoDS_Edge eNotSeam = edges.front(); + if ( helper.HasSeam() ) + { + list< TopoDS_Edge >::iterator eIt = edges.begin(); + while ( helper.IsRealSeam( *eIt )) ++eIt; + if ( eIt != edges.end() ) + eNotSeam = *eIt; + } + TopAbs_Orientation fOri = helper.GetSubShapeOri( face, eNotSeam ); + bool isForwad = ( fOri == eNotSeam.Orientation() || fOri >= TopAbs_INTERNAL ); - double otherSeamParam = 0; - bool isSeam = false; + // get all nodes from connected + const bool skipMedium = netgen::mparam.secondorder;//smDS->IsQuadratic(); + StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, skipMedium, &helper ); + const vector& points = fSide.GetUVPtStruct(); + if ( points.empty() ) + return false; // invalid node params? + smIdType i, nbSeg = fSide.NbSegments(); - // add segments + // remember EDGEs of fSide to treat only once + for ( int iE = 0; iE < fSide.NbEdges(); ++iE ) + visitedEdgeSM2Faces[ helper.GetMesh()->GetSubMesh( fSide.Edge(iE )) ].insert(faceSMDSId); - int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap ); + double otherSeamParam = 0; + bool isSeam = false; - for ( i = 0; i < nbSeg; ++i ) - { - const UVPtStruct& p1 = points[ i ]; - const UVPtStruct& p2 = points[ i+1 ]; + // add segments - if ( p1.node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX ) //an EDGE begins - { - isSeam = false; - if ( helper.IsRealSeam( p1.node->GetShapeID() )) + int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap ); + + for ( i = 0; i < nbSeg; ++i ) { - TopoDS_Edge e = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam ))); - isSeam = helper.IsRealSeam( e ); + const UVPtStruct& p1 = points[ i ]; + const UVPtStruct& p2 = points[ i+1 ]; + + if ( p1.node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX ) //an EDGE begins + { + isSeam = false; + if ( helper.IsRealSeam( p1.node->GetShapeID() )) + { + TopoDS_Edge e = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam ))); + isSeam = helper.IsRealSeam( e ); + if ( isSeam ) + { + otherSeamParam = helper.GetOtherParam( helper.GetPeriodicIndex() & 1 ? p2.u : p2.v ); + } + } + } + netgen::Segment seg; + // ng node ids + seg[0] = prevNgId; + seg[1] = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap ); + seg.edgenr = countEdgeId+1; + seg.si = countEdgeId+1; + seg.epgeominfo[0].edgenr = countEdgeId; + seg.epgeominfo[1].edgenr = countEdgeId; + // node param on curve + seg.epgeominfo[ 0 ].dist = p1.param; + seg.epgeominfo[ 1 ].dist = p2.param; + // uv on face + seg.epgeominfo[ 0 ].u = p1.u; + seg.epgeominfo[ 0 ].v = p1.v; + seg.epgeominfo[ 1 ].u = p2.u; + seg.epgeominfo[ 1 ].v = p2.v; + + //geomEdge = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam ))); + //seg.epgeominfo[ 0 ].edgenr = seg.epgeominfo[ 1 ].edgenr = occgeom.emap.FindIndex( geomEdge ); + + //seg.epgeominfo[ iEnd ].edgenr = edgeID; // = geom.emap.FindIndex(edge); + // seg.si = faceNgID; // = geom.fmap.FindIndex (face); + // seg.edgenr = ngMesh.GetNSeg() + 1; // segment id + ngMesh.AddSegment (seg); + + SMESH_TNodeXYZ np1( p1.node ), np2( p2.node ); + RestrictLocalSize( ngMesh, 0.5*(np1+np2), (np1-np2).Modulus() ); + + #ifdef DUMP_SEGMENTS + cout << "Segment: " << seg.edgenr << " on SMESH face " << meshDS->ShapeToIndex( face ) << endl + << "\tface index: " << seg.si << endl + << "\tp1: " << seg[0] << endl + << "\tp2: " << seg[1] << endl + << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl + << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl + //<< "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl + << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl + << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl; + //<< "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl; + #endif if ( isSeam ) { - otherSeamParam = helper.GetOtherParam( helper.GetPeriodicIndex() & 1 ? p2.u : p2.v ); + if ( helper.GetPeriodicIndex() && 1 ) { + seg.epgeominfo[ 0 ].u = otherSeamParam; + seg.epgeominfo[ 1 ].u = otherSeamParam; + swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v); + } else { + seg.epgeominfo[ 0 ].v = otherSeamParam; + seg.epgeominfo[ 1 ].v = otherSeamParam; + swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u); + } + swap( seg[0], seg[1] ); + swap( seg.epgeominfo[0].dist, seg.epgeominfo[1].dist ); + seg.edgenr = ngMesh.GetNSeg() + 1; // segment id + ngMesh.AddSegment( seg ); + #ifdef DUMP_SEGMENTS + cout << "Segment: " << seg.edgenr << endl + << "\t is SEAM (reverse) of the previous. " + << " Other " << (helper.GetPeriodicIndex() && 1 ? "U" : "V") + << " = " << otherSeamParam << endl; + #endif + } + else if ( fOri == TopAbs_INTERNAL ) + { + swap( seg[0], seg[1] ); + swap( seg.epgeominfo[0], seg.epgeominfo[1] ); + seg.edgenr = ngMesh.GetNSeg() + 1; // segment id + ngMesh.AddSegment( seg ); + #ifdef DUMP_SEGMENTS + cout << "Segment: " << seg.edgenr << endl << "\t is REVERSE of the previous" << endl; + #endif } } - } - netgen::Segment seg; - // ng node ids - seg[0] = prevNgId; - seg[1] = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap ); - // node param on curve - seg.epgeominfo[ 0 ].dist = p1.param; - seg.epgeominfo[ 1 ].dist = p2.param; - // uv on face - seg.epgeominfo[ 0 ].u = p1.u; - seg.epgeominfo[ 0 ].v = p1.v; - seg.epgeominfo[ 1 ].u = p2.u; - seg.epgeominfo[ 1 ].v = p2.v; - - //geomEdge = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam ))); - //seg.epgeominfo[ 0 ].edgenr = seg.epgeominfo[ 1 ].edgenr = occgeom.emap.FindIndex( geomEdge ); - - //seg.epgeominfo[ iEnd ].edgenr = edgeID; // = geom.emap.FindIndex(edge); - seg.si = faceNgID; // = geom.fmap.FindIndex (face); - seg.edgenr = ngMesh.GetNSeg() + 1; // segment id - ngMesh.AddSegment (seg); - - SMESH_TNodeXYZ np1( p1.node ), np2( p2.node ); - RestrictLocalSize( ngMesh, 0.5*(np1+np2), (np1-np2).Modulus() ); + } // loop on geomEdge ancestors -#ifdef DUMP_SEGMENTS - cout << "Segment: " << seg.edgenr << " on SMESH face " << meshDS->ShapeToIndex( face ) << endl - << "\tface index: " << seg.si << endl - << "\tp1: " << seg[0] << endl - << "\tp2: " << seg[1] << endl - << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl - << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl - //<< "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl - << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl - << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl; - //<< "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl; -#endif - if ( isSeam ) + if ( quadHelper ) // remember medium nodes of sub-meshes { - if ( helper.GetPeriodicIndex() && 1 ) { - seg.epgeominfo[ 0 ].u = otherSeamParam; - seg.epgeominfo[ 1 ].u = otherSeamParam; - swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v); - } else { - seg.epgeominfo[ 0 ].v = otherSeamParam; - seg.epgeominfo[ 1 ].v = otherSeamParam; - swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u); + SMDS_ElemIteratorPtr edges = smDS->GetElements(); + while ( edges->more() ) + { + const SMDS_MeshElement* e = edges->next(); + if ( !quadHelper->AddTLinks( static_cast< const SMDS_MeshEdge*>( e ))) + break; } - swap( seg[0], seg[1] ); - swap( seg.epgeominfo[0].dist, seg.epgeominfo[1].dist ); - seg.edgenr = ngMesh.GetNSeg() + 1; // segment id - ngMesh.AddSegment( seg ); -#ifdef DUMP_SEGMENTS - cout << "Segment: " << seg.edgenr << endl - << "\t is SEAM (reverse) of the previous. " - << " Other " << (helper.GetPeriodicIndex() && 1 ? "U" : "V") - << " = " << otherSeamParam << endl; -#endif } - else if ( fOri == TopAbs_INTERNAL ) - { - swap( seg[0], seg[1] ); - swap( seg.epgeominfo[0], seg.epgeominfo[1] ); - seg.edgenr = ngMesh.GetNSeg() + 1; // segment id - ngMesh.AddSegment( seg ); -#ifdef DUMP_SEGMENTS - cout << "Segment: " << seg.edgenr << endl << "\t is REVERSE of the previous" << endl; -#endif - } - } - } // loop on geomEdge ancestors - - if ( quadHelper ) // remember medium nodes of sub-meshes - { - SMDS_ElemIteratorPtr edges = smDS->GetElements(); - while ( edges->more() ) - { - const SMDS_MeshElement* e = edges->next(); - if ( !quadHelper->AddTLinks( static_cast< const SMDS_MeshEdge*>( e ))) - break; - } - } + + countEdgeId++; // Assume edge is iterate incrementally as done by the emap! + break; + + */ - break; + } /* scope of original code*/ } // case TopAbs_EDGE case TopAbs_FACE: { // FACE @@ -1439,6 +1518,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, int fID = occgeom.fmap.Add( geomFace ); if ( occgeom.facemeshstatus.Size() < fID ) occgeom.facemeshstatus.SetSize( fID ); occgeom.facemeshstatus[ fID-1 ] = netgen::FACE_MESHED_OK; + while ( fID < faceNgID ) // geomFace is already in occgeom.fmap, add a copy { fID = occgeom.fmap.Add( BRepBuilderAPI_Copy( geomFace, /*copyGeom=*/false )); @@ -2219,6 +2299,61 @@ void NETGENPlugin_Mesher::AddIntVerticesInSolids(const netgen::OCCGeometry& } // loop on solids with internal vertices } + +// Iterate in all edges of the geometry, and use RestrictLocalHLine function to set the value of allow edges in 1D mesh +SMESH_ComputeErrorPtr +NETGENPlugin_Mesher::AddSegmentsToMesh2(netgen::Mesh& ngMesh, + netgen::OCCGeometry& geom, + const TSideVector& wires, + SMESH_MesherHelper& helper, + vector< const SMDS_MeshNode* > & nodeVec, + const bool overrideMinH) +{ + map node2ngID; + + for ( size_t iW = 0; iW < wires.size(); ++iW ) + { + StdMeshers_FaceSidePtr wire = wires[ iW ]; + const vector& uvPtVec = wire->GetUVPtStruct(); + const smIdType nbSegments = wire->NbPoints() - 1; + + for ( int i = 0; i < nbSegments; ++i ) // loop on segments + { + const SMDS_MeshNode * n0 = uvPtVec[ i ].node; + const SMDS_MeshNode * n1 = uvPtVec[ i + 1 ].node; + const int posShapeID = n0->GetShapeID(); + // skip nodes on degenerated edges + if ( helper.IsDegenShape( posShapeID ) && helper.IsDegenShape( n1->GetShapeID() )) + continue; + + double size = SMESH_TNodeXYZ( n0 ).Distance( n1 ); + + netgen::Point3d p0(n0->X(), n0->Y(), n0->Z()); + netgen::Point3d p1(n1->X(), n1->Y(), n1->Z()); + // ngMesh.RestrictLocalHLine( p0, p1, size, overrideMinH ); + ngMesh.RestrictLocalH( p0, size, overrideMinH ); + + // if ( node2ngID.count(n0) == 0 ) + // { + // netgen::MeshPoint mp0( p0 ); + // ngMesh.AddPoint ( mp0, 1, netgen::FIXEDPOINT ); + // node2ngID.insert( make_pair( n0, ngMesh.GetNP() + 1 ) ); + // nodeVec.push_back( n0 ); + // } + // if ( node2ngID.count(n1) == 0 ) + // { + // netgen::MeshPoint mp1( p1 ); + // ngMesh.AddPoint ( mp1, 1, netgen::FIXEDPOINT ); + // node2ngID.insert( make_pair( n1, ngMesh.GetNP() + 1 ) ); + // nodeVec.push_back( n1 ); + // } + + } + } + + return TError(); +} + //================================================================================ /*! * \brief Fill netgen mesh with segments of a FACE @@ -2291,9 +2426,14 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh& ngMesh, } const int solidID = 0, faceID = geom.fmap.FindIndex( helper.GetSubShape() ); - if ( ngMesh.GetNFD() < 1 ) - ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceID, solidID, solidID, 0 )); - + // if ( ngMesh.GetNFD() < 1 ) + // ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceID, solidID, solidID, 0 )); + auto numEdges = geom.GetNEdges(); + auto numVertex = geom.GetNVertices(); + auto numFaces = geom.GetNFaces(); + printf("Geom elements: %d, %d, %d\n", numEdges, numVertex, numFaces ); + + // Index meshed edges that were already add to the mesh to avoid add then twice as required by netge! for ( size_t iW = 0; iW < wires.size(); ++iW ) { StdMeshers_FaceSidePtr wire = wires[ iW ]; @@ -2303,6 +2443,7 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh& ngMesh, // assure the 1st node to be in node2ngID, which is needed to correctly // "close chain of segments" (see below) in case if the 1st node is not // onVertex because it is on a Viscous layer + // printf("First reference node is: %0.1lf, %0.1lf, %0.1lf\n", uvPtVec[ 0 ].node->X(), uvPtVec[ 0 ].node->Y(), uvPtVec[ 0 ].node->Z() ); node2ngID.insert( make_pair( uvPtVec[ 0 ].node, ngMesh.GetNP() + 1 )); // compute length of every segment @@ -2310,10 +2451,13 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh& ngMesh, for ( int i = 0; i < nbSegments; ++i ) segLen[i] = SMESH_TNodeXYZ( uvPtVec[ i ].node ).Distance( uvPtVec[ i+1 ].node ); - int edgeID = 1, posID = -2; + int edgeID = -1, posID = -2; bool isInternalWire = false; double vertexNormPar = 0; const int prevNbNGSeg = ngMesh.GetNSeg(); + int previusEdgeId = -1; + const SMDS_MeshNode * lastIndexedNode; + for ( int i = 0; i < nbSegments; ++i ) // loop on segments { // Add the first point of a segment @@ -2328,11 +2472,29 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh& ngMesh, helper.IsDegenShape( uvPtVec[ i+1 ].node->GetShapeID() )) continue; + double normParam = 0.5 * ( uvPtVec[ i ].normParam + uvPtVec[ i+1 ].normParam ); + int edgeIndexInWire = wire->EdgeIndex( normParam ); + const TopoDS_Edge& edge = wire->Edge( edgeIndexInWire ); + edgeID = geom.GetEdge(edge).nr; + int ngID1 = ngMesh.GetNP() + 1, ngID2 = ngID1+1; - if ( onVertex || ( !wasNgMeshEmpty && onEdge ) || helper.IsRealSeam( posShapeID )) - ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second; - + // Treat seam edge differently + // if ( helper.IsRealSeam( posShapeID ) && node2ngID.count( n ) != 0 ) + // continue; + + if ( onVertex || ( !wasNgMeshEmpty && onEdge ) || helper.IsRealSeam( posShapeID ) ) + { + ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second; + + // printf("OnVertex: %d \n", ngID1 ); + // if ( node2ngID.count( n ) == 0 ) + // node2ngID.insert( make_pair( n, ngID1 )); + // ngID1 = node2ngID[ n ]; + + // printf("OnVertex: %d, %0.1lf, %0.1lf, %0.1lf\n", ngID1, n->X(), n->Y(), n->Z() ); + } + if ( ngID1 > ngMesh.GetNP() ) { netgen::MeshPoint mp( netgen::Point<3> (n->X(), n->Y(), n->Z()) ); @@ -2349,70 +2511,64 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh& ngMesh, } } - // Add the segment - netgen::Segment seg; seg[0] = ngID1; // ng node id seg[1] = ngID2; // ng node id - seg.edgenr = ngMesh.GetNSeg() + 1; // ng segment id - seg.si = faceID; // = geom.fmap.FindIndex (face); - for ( int iEnd = 0; iEnd < 2; ++iEnd) - { - const UVPtStruct& pnt = uvPtVec[ i + iEnd ]; - seg.epgeominfo[ iEnd ].dist = pnt.param; // param on curve - seg.epgeominfo[ iEnd ].u = pnt.u; - seg.epgeominfo[ iEnd ].v = pnt.v; - // find out edge id and node parameter on edge - onVertex = ( pnt.normParam + 1e-10 > vertexNormPar ); - if ( onVertex || posShapeID != posID ) + // netgen6 + seg.edgenr = edgeID+1; // ng segment id + seg.si = edgeID+1; // = geom.fmap.FindIndex (face); + + const UVPtStruct& pnt0 = uvPtVec[ i ]; + const UVPtStruct& pnt1 = uvPtVec[ i + 1 ]; + + seg.epgeominfo[ 0 ].dist = pnt0.param; // param on curve + seg.epgeominfo[ 0 ].u = pnt0.u; + seg.epgeominfo[ 0 ].v = pnt0.v; + seg.epgeominfo[ 0 ].edgenr = edgeID; + + seg.epgeominfo[ 1 ].dist = pnt1.param; // param on curve + seg.epgeominfo[ 1 ].u = pnt1.u; + seg.epgeominfo[ 1 ].v = pnt1.v; + seg.epgeominfo[ 1 ].edgenr = edgeID; // + + // seg.singedge_left = geom.GetEdge(edge).properties.hpref; + // seg.singedge_right = geom.GetEdge(edge).properties.hpref; + // seg.domin = geom.GetEdge(edge).domin+1; + // seg.domout = geom.GetEdge(edge).domout+1; + + ngMesh.AddSegment (seg); + { - // get edge id - double normParam = pnt.normParam; - if ( onVertex ) - normParam = 0.5 * ( uvPtVec[ i ].normParam + uvPtVec[ i+1 ].normParam ); - int edgeIndexInWire = wire->EdgeIndex( normParam ); - vertexNormPar = wire->LastParameter( edgeIndexInWire ); - const TopoDS_Edge& edge = wire->Edge( edgeIndexInWire ); - edgeID = geom.emap.FindIndex( edge ); - posID = posShapeID; - isInternalWire = ( edge.Orientation() == TopAbs_INTERNAL ); - // if ( onVertex ) // param on curve is different on each of two edges - // seg.epgeominfo[ iEnd ].dist = helper.GetNodeU( edge, pnt.node ); + // restrict size of elements near the segment + SMESH_TNodeXYZ np1( n ), np2( uvPtVec[ i+1 ].node ); + // get an average size of adjacent segments to avoid sharp change of + // element size (regression on issue 0020452, note 0010898) + int iPrev = SMESH_MesherHelper::WrapIndex( i-1, (int) nbSegments ); + int iNext = SMESH_MesherHelper::WrapIndex( i+1, (int) nbSegments ); + double sumH = segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ]; + int nbSeg = ( int( segLen[ iPrev ] > sumH / 100.) + + int( segLen[ i ] > sumH / 100.) + + int( segLen[ iNext ] > sumH / 100.)); + if ( nbSeg > 0 ) + RestrictLocalSize( ngMesh, 0.5*(np1+np2), sumH / nbSeg, overrideMinH ); } - seg.epgeominfo[ iEnd ].edgenr = edgeID; // = geom.emap.FindIndex(edge); - } - ngMesh.AddSegment (seg); - { - // restrict size of elements near the segment - SMESH_TNodeXYZ np1( n ), np2( uvPtVec[ i+1 ].node ); - // get an average size of adjacent segments to avoid sharp change of - // element size (regression on issue 0020452, note 0010898) - int iPrev = SMESH_MesherHelper::WrapIndex( i-1, (int) nbSegments ); - int iNext = SMESH_MesherHelper::WrapIndex( i+1, (int) nbSegments ); - double sumH = segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ]; - int nbSeg = ( int( segLen[ iPrev ] > sumH / 100.) + - int( segLen[ i ] > sumH / 100.) + - int( segLen[ iNext ] > sumH / 100.)); - if ( nbSeg > 0 ) - RestrictLocalSize( ngMesh, 0.5*(np1+np2), sumH / nbSeg, overrideMinH ); - } - if ( isInternalWire ) - { - swap (seg[0], seg[1]); - swap( seg.epgeominfo[0], seg.epgeominfo[1] ); - seg.edgenr = ngMesh.GetNSeg() + 1; // segment id - ngMesh.AddSegment (seg); - } + // lastIndexedNode = pnt1.node; + + // printf("Segments init %d %d %0.1lf, %0.1lf, %0.1lf\n", edgeID, ngID1, pnt0.node->X(), pnt0.node->Y(), pnt0.node->Z() ); + // printf("Segments end %d %d %0.1lf, %0.1lf, %0.1lf\n", edgeID, ngID2, pnt1.node->X(), pnt1.node->Y(), pnt1.node->Z() ); + } // loop on segments on a wire + // printf("Out of the segment loop \n"); // close chain of segments if ( nbSegments > 0 ) { netgen::Segment& lastSeg = ngMesh.LineSegment( ngMesh.GetNSeg() - int( isInternalWire )); const SMDS_MeshNode * lastNode = uvPtVec.back().node; + lastSeg[1] = node2ngID.insert( make_pair( lastNode, lastSeg[1] )).first->second; if ( lastSeg[1] > ngMesh.GetNP() ) { @@ -2458,16 +2614,34 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh& ngMesh, } // loop on WIREs of a FACE // add a segment instead of an internal vertex - if ( wasNgMeshEmpty ) - { - NETGENPlugin_Internals intShapes( *helper.GetMesh(), helper.GetSubShape(), /*is3D=*/false ); - AddIntVerticesInFaces( geom, ngMesh, nodeVec, intShapes ); - } + // if ( wasNgMeshEmpty ) + // { + // NETGENPlugin_Internals intShapes( *helper.GetMesh(), helper.GetSubShape(), /*is3D=*/false ); + // AddIntVerticesInFaces( geom, ngMesh, nodeVec, intShapes ); + // } + + + // Get the add segment and check + + printf("Going to call surface\n"); ngMesh.CalcSurfacesOfNode(); + printf("End to call surface, %d,%d\n", ngMesh.GetNP(), ngMesh.GetNSeg() ); return TError(); } + +// int NETGENPlugin_Mesher::CheckSegmentChain(const netgen::OCCGeometry& occgeo, +// netgen::Mesh& ngMesh ) +// { + +// auto segments = occgeo.GetFace(0).GetBoundary(); +// for (auto s : segments) +// { + +// } +// } + //================================================================================ /*! * \brief Fill SMESH mesh according to contents of netgen mesh @@ -2494,6 +2668,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, int nbSeg = ngMesh.GetNSeg(); int nbFac = ngMesh.GetNSE(); int nbVol = ngMesh.GetNE(); + printf( "Numbers: %d,%d,%d,%d\n", nbNod, nbSeg, nbFac, nbVol ); SMESHDS_Mesh* meshDS = sMesh.GetMeshDS(); // quadHelper is used for either @@ -3262,6 +3437,7 @@ void NETGENPlugin_Mesher::SetBasicMeshParameters( NETGENPlugin_NetgenLibWrapper& _ngMesh->ClearFaceDescriptors(); // we make descriptors our-self + /* This line breaks when using _simpleHyp for 1D2D version*/ if ( !mparams.uselocalh ) // mparams.grading is not taken into account yet _ngMesh->LocalHFunction().SetGrading( mparams.grading ); @@ -3273,7 +3449,7 @@ void NETGENPlugin_Mesher::SetBasicMeshParameters( NETGENPlugin_NetgenLibWrapper& double segSize = _simpleHyp->GetLocalLength(); for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE ) { - const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE)); + const TopoDS_Edge& e = TopoDS::Edge(occgeo.emap(iE)); if ( nbSeg ) segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 ); setLocalSize( e, segSize, *_ngMesh ); @@ -3420,7 +3596,8 @@ void NETGENPlugin_Mesher::InitialSetupSA( NETGENPlugin_NetgenLibWrapper& ngLib, updateTriangulation( _shape ); Bnd_Box bb; - BRepBndLib::Add (_shape, bb); + // BRepBndLib::Add (_shape, bb); + GEOMUtils::PreciseBoundingBox( _shape, bb ); double x1,y1,z1,x2,y2,z2; bb.Get (x1,y1,z1,x2,y2,z2); netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1); @@ -3665,6 +3842,7 @@ bool NETGENPlugin_Mesher::Compute() //int nbSeg = _ngMesh->GetNSeg(); int nbFac = _ngMesh->GetNSE(); int nbVol = _ngMesh->GetNE(); + bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) ); // Feed back the SMESHDS with the generated Nodes and Elements diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx index 173ca08..b56a988 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx @@ -293,6 +293,13 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher std::vector< const SMDS_MeshNode* > & nodeVec, const bool overrideMinH=true); + static SMESH_ComputeErrorPtr AddSegmentsToMesh2(netgen::Mesh& ngMesh, + netgen::OCCGeometry& geom, + const TSideVector& wires, + SMESH_MesherHelper& helper, + std::vector< const SMDS_MeshNode* > & nodeVec, + const bool overrideMinH=true); + void SetDefaultParameters(); static SMESH_ComputeErrorPtr ReadErrors(const std::vector< const SMDS_MeshNode* >& nodeVec); diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx index 6d37205..1ea1a83 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx @@ -549,6 +549,9 @@ bool NETGENPlugin_NETGEN_2D_ONLY::MapSegmentsToEdges(SMESH_Mesh& aMesh, const To occgeom.face_maxh_modified = 0; occgeom.face_maxh.SetSize(1); occgeom.face_maxh = netgen::mparam.maxh; + occgeom.BuildFMapFromPrefilled(); + + printf("Calling NETGENPlugin_NETGEN_2D_ONLY::MapSegmentsToEdges\n"); // Set the face descriptor const int solidID = 0, faceID = 1; /*always 1 because faces are meshed one by one*/ @@ -922,7 +925,7 @@ void NETGENPlugin_NETGEN_2D_ONLY::FillNodesAndElements( SMESH_Mesh& aMesh, SMESH 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)); + SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2)); nodeVec[ ngID ] = node; } @@ -1057,17 +1060,13 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, return setMaxh; // prepare occgeom - netgen::OCCGeometry occgeom; - occgeom.shape = F; - occgeom.fmap.Add( F ); - 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; + printf( "going to define occgeom\n"); + netgen::OCCGeometry occgeom; + + occgeom.shape = F; + occgeom.BuildFMapFace(); + // ------------------------- // Fill netgen mesh // ------------------------- @@ -1102,8 +1101,12 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, } nodeVec.clear(); - faceErr = aMesher.AddSegmentsToMesh( *ngMesh, occgeom, wires, helper, nodeVec, + // faceErr = aMesher.AddSegmentsToMesh( *ngMesh, occgeom, wires, helper, nodeVec, + // /*overrideMinH=*/!_hypParameters); + + faceErr = aMesher.AddSegmentsToMesh2( *ngMesh, occgeom, wires, helper, nodeVec, /*overrideMinH=*/!_hypParameters); + if ( faceErr && !faceErr->IsOK() ) break; @@ -1113,15 +1116,19 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, // ------------------------- // Generate surface mesh // ------------------------- - - const int startWith = MESHCONST_MESHSURFACE; + const int startWith = netgen::MESHCONST_MESHEDGES; const int endWith = toOptimize ? MESHCONST_OPTSURFACE : MESHCONST_MESHSURFACE; + // const int startWith = MESHCONST_MESHSURFACE; + // const int endWith = toOptimize ? MESHCONST_OPTSURFACE : MESHCONST_MESHSURFACE; + // const int endWith = MESHCONST_MESHSURFACE; + SMESH_Comment str; try { OCC_CATCH_SIGNALS; - + printf("Going to generate mesh\n"); err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh); + printf("Mesh Generated: %d,%d,%d \n", ngMesh->GetNP(), ngMesh->GetNSeg(), ngMesh->GetNSE() ); if ( netgen::multithread.terminate ) return false; diff --git a/src/NETGENPlugin/NETGENPlugin_Remesher_2D.cxx b/src/NETGENPlugin/NETGENPlugin_Remesher_2D.cxx index 2162e04..c046ed5 100644 --- a/src/NETGENPlugin/NETGENPlugin_Remesher_2D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Remesher_2D.cxx @@ -231,7 +231,7 @@ namespace void HoleFiller::AddHoleBordersAndEdges( Ng_STL_Geometry * ngStlGeo, bool toAddEdges ) { - nglib::readedges.SetSize(0); + // nglib::readedges.SetSize(0); for ( size_t i = 0; i < myHole.size(); ++i ) for ( size_t iP = 1; iP < myHole[i].size(); ++iP ) @@ -516,9 +516,9 @@ namespace n->setIsMarked( true ); SMESH_NodeXYZ p( n ); - int id = stlGeom->GetPointNum( netgen::Point<3>( p.X(),p.Y(),p.Z() )); - if ( id > 0 ) - stlGeom->SetLineEndPoint( id ); + // int id = stlGeom->GetPointNum( netgen::Point<3>( p.X(),p.Y(),p.Z() )); + // if ( id > 0 ) + // stlGeom->SetLineEndPoint( id ); } } }