From a8606239fd767f5f572f411ddbefb24439d26ff7 Mon Sep 17 00:00:00 2001 From: eap Date: Mon, 6 Sep 2021 15:52:51 +0300 Subject: [PATCH] bos #18711 [CEA 18704] COTECH action 117.4: NETGEN 6 Integration in SALOME cleanup + PATCH --- src/NETGEN/netgen62ForSalome.patch | 293 ++++++++++++++++++ src/NETGENPlugin/NETGENPlugin_Mesher.cxx | 68 +++- .../NETGENPlugin_NETGEN_2D_ONLY.cxx | 12 +- src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx | 7 +- 4 files changed, 368 insertions(+), 12 deletions(-) create mode 100644 src/NETGEN/netgen62ForSalome.patch diff --git a/src/NETGEN/netgen62ForSalome.patch b/src/NETGEN/netgen62ForSalome.patch new file mode 100644 index 0000000..20c187b --- /dev/null +++ b/src/NETGEN/netgen62ForSalome.patch @@ -0,0 +1,293 @@ +diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp +index 314d405a..7c912bcc 100644 +--- a/libsrc/occ/occgenmesh.cpp ++++ b/libsrc/occ/occgenmesh.cpp +@@ -49,7 +49,7 @@ namespace netgen + double ComputeH (double kappa, const MeshingParameters & mparam) + { + kappa *= mparam.curvaturesafety; +- /* ++ /**/ + double hret; + + if (mparam.maxh * kappa < 1) +@@ -61,7 +61,7 @@ namespace netgen + hret = mparam.maxh; + + return hret; +- */ ++ /**/ + // return min(mparam.maxh, 1/kappa); + return (mparam.maxh*kappa < 1) ? mparam.maxh : 1/kappa; + } +@@ -253,10 +253,12 @@ namespace netgen + { + oldpnt = pnt; + pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0)); ++ // !!! no more than 1 segment per /DIVIDEEDGESECTIONS + hvalue[i] = hvalue[i-1] + +- 1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))* +- pnt.Distance(oldpnt); ++ min( 1.0, ++ 1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))* ++ pnt.Distance(oldpnt)); + + //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) + // << " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl; + +@@ -299,7 +301,7 @@ namespace netgen + cout << "CORRECTED" << endl; + ps.SetSize (nsubedges-2); + params.SetSize (nsubedges); +- params[nsubedges] = s1; ++ params[nsubedges-1] = s1; + } + } + +@@ -323,6 +325,8 @@ namespace netgen + + double eps = 1e-6 * geom.GetBoundingBox().Diam(); + ++ int first_vp = mesh.GetNP()+1; ++ + tsearch.Start(); + for (int i = 1; i <= nvertices; i++) + { +@@ -481,22 +485,107 @@ namespace netgen + + if (!merge_solids) + { +- pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1; +- pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1; ++ //pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1; ++ //pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1; ++ const int dpi = PointIndex::BASE-1; ++ MeshPoint dfltP ( Point<3> ( 0, 0, 0 ) ); ++ PointIndex *ipp[] = { &pnums[0], &pnums.Last() }; ++ TopoDS_Iterator vIt( edge, false ); ++ TopoDS_Vertex v[2]; ++ v[0] = TopoDS::Vertex( vIt.Value() ); vIt.Next(); ++ v[1] = TopoDS::Vertex( vIt.Value() ); ++ if ( v[0].Orientation() == TopAbs_REVERSED ) ++ std::swap( v[0], v[1] ); ++ for ( int i = 0; i < 2; ++i) ++ { ++ PointIndex &ip = *ipp[i]; ++ ip = geom.vmap.FindIndex ( v[i] ) + dpi; ++ if ( ip == dpi || ip > nvertices + dpi ) ++ { ++ PointIndex iv = ip; ++ if ( ip == dpi ) ++ ip = iv = const_cast(geom).vmap.Add( v[i] ); ++ gp_Pnt pnt = BRep_Tool::Pnt( v[i] ); ++ MeshPoint mp( Point<3>(pnt.X(), pnt.Y(), pnt.Z()) ); ++ for (PointIndex pi = 1; pi < first_vp; pi++) ++ if ( Dist2 (mesh.Point(pi), Point<3>(mp)) < 1e-100 ) ++ { ++ ip = pi; ++ if ( mesh.Point(ip).GetLayer() != dfltP.GetLayer() && mesh.Point(ip).GetLayer() != iv ) ++ continue; ++ if ( mesh.Point(ip).GetLayer() == dfltP.GetLayer()) ++ mesh.Point(ip) = MeshPoint( mesh.Point(ip), iv ); ++ break; ++ } ++ } ++ else ++ { ++ ip += first_vp - 1; ++ } ++ } + } + else + { +- Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge))); +- Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge))); ++ TopoDS_Iterator vIt( edge, false ); ++ TopoDS_Vertex v1 = TopoDS::Vertex( vIt.Value() ); vIt.Next(); ++ TopoDS_Vertex v2 = TopoDS::Vertex( vIt.Value() ); ++ if ( v1.Orientation() == TopAbs_REVERSED ) ++ std::swap( v1, v2 ); ++ const bool isClosedEdge = v1.IsSame( v2 ); ++ ++ Point<3> fp = occ2ng (BRep_Tool::Pnt (v1)); ++ Point<3> lp = occ2ng (BRep_Tool::Pnt (v2)); ++ double tol2 = std::min( eps*eps, 1e-6 * Dist2( fp, lp )); ++ if ( isClosedEdge ) ++ tol2 = BRep_Tool::Tolerance( v1 ) * BRep_Tool::Tolerance( v1 ); + + pnums[0] = PointIndex::INVALID; + pnums.Last() = PointIndex::INVALID; + for (PointIndex pi : vertexrange) + { +- if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi; +- if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi; ++ if (Dist2 (mesh[pi], fp) < tol2) pnums[0] = pi; ++ if (Dist2 (mesh[pi], lp) < tol2) pnums.Last() = pi; + } +- } ++ if (( isClosedEdge && pnums[0] != pnums.Last() ) || ++ ( !isClosedEdge && pnums[0] == pnums.Last() )) ++ pnums[0] = pnums.Last() = PointIndex::INVALID; ++ if ( pnums[0] < PointIndex::BASE || pnums.Last() < PointIndex::BASE ) ++ { ++ // take into account a possible large gap between a vertex and an edge curve ++ // end and a large vertex tolerance covering the whole edge ++ if ( pnums[0] < PointIndex::BASE ) ++ { ++ double tol = BRep_Tool::Tolerance( v1 ); ++ for (PointIndex pi = 1; pi < first_ep; pi++) ++ if (pi != pnums.Last() && Dist2 (mesh[pi], fp) < 2*tol*tol) ++ pnums[0] = pi; ++ ++ if ( pnums[0] < PointIndex::BASE ) ++ pnums[0] = first_ep-1 - nvertices + geom.vmap.FindIndex ( v1 ); ++ } ++ if ( isClosedEdge ) ++ { ++ pnums.Last() = pnums[0]; ++ } ++ else ++ { ++ if ( pnums.Last() < PointIndex::BASE ) ++ { ++ double tol = BRep_Tool::Tolerance( v2 ); ++ for (PointIndex pi = 1; pi < first_ep; pi++) ++ if (pi != pnums[0] && Dist2 (mesh[pi], lp) < 2*tol*tol) ++ pnums.Last() = pi; ++ ++ if ( pnums.Last() < PointIndex::BASE ) ++ pnums.Last() = first_ep-1-nvertices + geom.vmap.FindIndex ( v2 ); ++ } ++ ++ if ( Dist2( fp, mesh[PointIndex(pnums[0])]) > ++ Dist2( lp, mesh[PointIndex(pnums.Last())])) ++ std::swap( pnums[0], pnums.Last() ); ++ } ++ } ++ } + + for (size_t i = 1; i <= mp.Size(); i++) + { +@@ -505,17 +594,19 @@ namespace netgen + + // for (PointIndex j = first_ep; j < mesh.Points().End(); j++) + for (PointIndex j = first_ep; j < *mesh.Points().Range().end(); j++) ++ { ++ if (!merge_solids && mesh.Point(j).GetLayer() != geomedgenr ) continue; + if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps) + { + exists = true; + pnums[i] = j; + break; + } +- ++ } + tsearch.Stop(); + + if (!exists) +- pnums[i] = mesh.AddPoint (mp[i-1]); ++ pnums[i] = mesh.AddPoint (mp[i-1], geomedgenr); + } + if(geom.enames.Size() && geom.enames[curr-1] != "") + mesh.SetCD2Name(geomedgenr, geom.enames[curr-1]); +@@ -599,6 +690,9 @@ namespace netgen + // << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl; + // exit(10); + ++ for (int j = 1; j <= mesh.GetNP(); j++) // support SALOME fuse edges: set level to zero ++ mesh.Point(j) = MeshPoint( (Point<3>&) mesh.Point(j) ); ++ + mesh.CalcSurfacesOfNode(); + multithread.task = savetask; + } +@@ -652,7 +746,7 @@ namespace netgen + + Box<3> bb = geom.GetBoundingBox(); + +- int projecttype = PLANESPACE; ++ //int projecttype = PLANESPACE; + + static Timer tinit("init"); + tinit.Start(); +diff --git a/libsrc/occ/occmeshsurf.cpp b/libsrc/occ/occmeshsurf.cpp +index 7fa5d237..9e05fd95 100644 +--- a/libsrc/occ/occmeshsurf.cpp ++++ b/libsrc/occ/occmeshsurf.cpp +@@ -173,50 +173,6 @@ namespace netgen + gp_Vec du, dv; + occface->D1 (geominfo1.u, geominfo1.v, pnt, du, dv); + +- // static Timer t("occ-defintangplane calculations"); +- // RegionTimer reg(t); +- +- Mat<3,2> D1_; +- D1_(0,0) = du.X(); D1_(1,0) = du.Y(); D1_(2,0) = du.Z(); +- D1_(0,1) = dv.X(); D1_(1,1) = dv.Y(); D1_(2,1) = dv.Z(); +- auto D1T_ = Trans(D1_); +- auto D1TD1_ = D1T_*D1_; +- if (Det (D1TD1_) == 0) throw SingularMatrixException(); +- Mat<2,2> DDTinv_; +- CalcInverse (D1TD1_, DDTinv_); +- +- Mat<3,2> Y_; +- Vec<3> y1_ = (ap2-ap1).Normalize(); +- Vec<3> y2_ = Cross(n, y1_).Normalize(); +- for (int i = 0; i < 3; i++) +- { +- Y_(i,0) = y1_(i); +- Y_(i,1) = y2_(i); +- } +- +- auto A_ = DDTinv_ * D1T_ * Y_; +- Mat<2,2> Ainv_; +- if (Det(A_) == 0) throw SingularMatrixException(); +- CalcInverse (A_, Ainv_); +- +- Vec<2> temp_ = Ainv_ * (psp2-psp1); +- double r_ = temp_.Length(); +- Mat<2,2> R_; +- R_(0,0) = temp_(0)/r_; +- R_(1,0) = temp_(1)/r_; +- R_(0,1) = -R_(1,0); +- R_(1,1) = R_(0,0); +- +- A_ = A_ * R_; +- Ainv_ = Trans(R_) * Ainv_; +- +- Amat = A_; +- Amatinv = Ainv_; +- +- // temp = Amatinv * (psp2-psp1); +- +- +-#ifdef OLD + DenseMatrix D1(3,2), D1T(2,3), DDTinv(2,2); + D1(0,0) = du.X(); D1(1,0) = du.Y(); D1(2,0) = du.Z(); + D1(0,1) = dv.X(); D1(1,1) = dv.Y(); D1(2,1) = dv.Z(); +@@ -289,8 +245,7 @@ namespace netgen + } + // cout << "=?= Ainv = " << endl << Ainv << endl; + temp = Amatinv * (psp2-psp1); +- cout << " =?= Amatinv = " << Amatinv << endl; +-#endif ++ // cout << " =?= Amatinv = " << Amatinv << endl; + }; + + } +@@ -380,6 +335,7 @@ namespace netgen + + double u = gi.u; + double v = gi.v; ++ if ( 0 ) { // exclude Newton's method + gp_Pnt x = occface->Value (u,v); + + if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return; +@@ -429,7 +385,7 @@ namespace netgen + } + } + while (count < 20); +- ++ } + + // Newton did not converge, use OCC projection + diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx index ed60210..5853a7e 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx @@ -669,6 +669,12 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp) } } } + +#ifdef NETGEN_V6 + + netgen::mparam.closeedgefac = 2; + +#endif } //============================================================================= @@ -732,7 +738,7 @@ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo, if ( faceNgID >= 1 ) { #ifdef NETGEN_V6 - occgeo.SetFaceMaxH(faceNgID-1, val); + occgeo.SetFaceMaxH(faceNgID-1, val, netgen::mparam); #else occgeo.SetFaceMaxH(faceNgID, val); #endif @@ -801,7 +807,7 @@ void NETGENPlugin_Mesher::SetLocalSizeForChordalError( netgen::OCCGeometry& occg double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() )); double size = elemSizeForChordalError( _chordalError, 1 / maxCurv ); #ifdef NETGEN_V6 - occgeo.SetFaceMaxH( i-1, size * sizeCoef ); + occgeo.SetFaceMaxH( i-1, size * sizeCoef, netgen::mparam ); #else occgeo.SetFaceMaxH( i, size * sizeCoef ); #endif @@ -3814,7 +3820,7 @@ NETGENPlugin_Mesher::ReadErrors(const vector& nodeVec) ok = ok && file.getInts( three2 ); for ( int i = 0; ok && i < 3; ++i ) ok = ( three1[i] < nbNodes && nodeVec[ three1[i]]); - for ( int i = 0; ok && i < 3; ++i ) + for ( int i = 0; ok && i < 3; ++i ) ok = ( three2[i] < nbNodes && nodeVec[ three2[i]]); if ( ok ) { @@ -3865,6 +3871,59 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh ) << "mesh = smesh.Mesh()" << std::endl << std::endl; using namespace netgen; + +#ifdef NETGEN_V6 + + for ( int i = 1; i <= ngMesh->GetNP(); i++) + { + const Point3d & p = ngMesh->Point(i); + outfile << "mesh.AddNode( "; + outfile << p.X() << ", "; + outfile << p.Y() << ", "; + outfile << p.Z() << ") ## "<< i << std::endl; + } + + int nbDom = ngMesh->GetNDomains(); + for ( int i = 0; i < nbDom; ++i ) + outfile<< "grp" << i+1 << " = mesh.CreateEmptyGroup( SMESH.FACE, 'domain"<< i+1 << "')"<< std::endl; + + for (int i = 1; i <= ngMesh->GetNSE(); i++) + { + outfile << "mesh.AddFace([ "; + Element2d sel = ngMesh->SurfaceElement(i); + for (int j = 1; j <= sel.GetNP(); j++) + outfile << sel.PNum(j) << ( j < sel.GetNP() ? ", " : " ])"); + if ( sel.IsDeleted() ) outfile << " ## IsDeleted "; + outfile << std::endl; + + if (sel.GetIndex()) + { + if ( int dom1 = ngMesh->GetFaceDescriptor(sel.GetIndex ()).DomainIn()) + outfile << "grp"<< dom1 <<".Add([ " << i << " ])" << std::endl; + if ( int dom2 = ngMesh->GetFaceDescriptor(sel.GetIndex ()).DomainOut()) + outfile << "grp"<< dom2 <<".Add([ " << i << " ])" << std::endl; + } + } + + for (int i = 1; i <= ngMesh->GetNE(); i++) + { + Element el = ngMesh->VolumeElement(i); + outfile << "mesh.AddVolume([ "; + for (int j = 1; j <= el.GetNP(); j++) + outfile << el.PNum(j) << ( j < el.GetNP() ? ", " : " ])"); + outfile << std::endl; + } + + for (int i = 1; i <= ngMesh->GetNSeg(); i++) + { + const Segment & seg = ngMesh->LineSegment (i); + outfile << "mesh.AddEdge([ " + << seg[0]+1 << ", " + << seg[1]+1 << " ])" << std::endl; + } + +#else //////// V 5 + PointIndex pi; for (pi = PointIndex::BASE; pi < ngMesh->GetNP()+PointIndex::BASE; pi++) @@ -3914,6 +3973,9 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh ) << seg[0] << ", " << seg[1] << " ])" << std::endl; } + +#endif + std::cout << "Write " << pyFile << std::endl; } diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx index 84d4363..cef4819 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx @@ -65,7 +65,9 @@ namespace nglib { namespace netgen { NETGENPLUGIN_DLL_HEADER extern MeshingParameters mparam; +#ifdef NETGEN_V5 extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh); +#endif } using namespace std; @@ -280,15 +282,19 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, } // set local size depending on curvature and NOT closeness of EDGEs #ifdef NETGEN_V6 - const double factor = *netgen::mparam.closeedgefac; - netgen::mparam.closeedgefac = std::nullopt; + const double factor = 2; //netgen::occparam.resthcloseedgefac; #else const double factor = netgen::occparam.resthcloseedgefac; netgen::occparam.resthcloseedgeenable = false; - //netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading; + netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading; #endif occgeoComm.face_maxh = netgen::mparam.maxh; +#ifdef NETGEN_V6 + netgen::OCCParameters occparam; + netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0], netgen::mparam, occparam ); +#else netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0] ); +#endif occgeoComm.emap.Clear(); occgeoComm.vmap.Clear(); diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx index de45d54..88abe10 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx @@ -76,18 +76,13 @@ #include #endif #ifdef NETGEN_V6 -#include +#include #endif namespace nglib { #include } namespace netgen { -#ifdef NETGEN_V5 - extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, MeshingParameters&, int, int); -#else - extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*); -#endif NETGENPLUGIN_DLL_HEADER extern MeshingParameters mparam; -- 2.30.2