ForkMaster/netgenplugin
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

bos #18711 [CEA 18704] COTECH action 117.4: NETGEN 6 Integration in SALOME

Browse Source 
cleanup + PATCH
This commit is contained in:
eap 2021-09-06 15:52:51 +03:00
parent 74525024bf
commit a8606239fd
4 changed files with 368 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

293
src/NETGEN/netgen62ForSalome.patch Normal file
View File

@ -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 <edge length>/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<OCCGeometry&>(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

68
src/NETGENPlugin/NETGENPlugin_Mesher.cxx
View File

@ -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<const SMDS_MeshNode* >& 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;
}

12
src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx
View File

@ -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();

7
src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx
View File

@ -76,18 +76,13 @@
#include <ngexception.hpp>
#endif
#ifdef NETGEN_V6
#include <exception.hpp>
#include <core/exception.hpp>
#endif
namespace nglib {
#include <nglib.h>
}
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;