ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-04-10 23:27:28 +05:00
Code Issues Packages Projects Releases Wiki Activity

little OCC-meshing cleanup

Browse Source
This commit is contained in:
Joachim Schöberl 2019-07-28 20:22:48 +02:00
parent 954cae2686
commit 321bee9b02
4 changed files with 1235 additions and 1231 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
libsrc/meshing/meshing2.cpp
View File

@ -21,6 +21,7 @@ namespace netgen
Meshing2 :: Meshing2 (const MeshingParameters & mp, const Box<3> & aboundingbox) Meshing2 :: Meshing2 (const MeshingParameters & mp, const Box<3> & aboundingbox)
{ {
static Timer t("Mesing2::Meshing2"); RegionTimer r(t);
boundingbox = aboundingbox; boundingbox = aboundingbox;
LoadRules (NULL, mp.quad); LoadRules (NULL, mp.quad);

82
libsrc/occ/occgenmesh.cpp
View File

@ -41,24 +41,14 @@ namespace netgen
return 1e99; return 1e99;
} }
double Line :: Length ()
{
return (p1-p0).Length();
}
inline Point<3> occ2ng (const gp_Pnt & p) inline Point<3> occ2ng (const gp_Pnt & p)
{ {
return Point<3> (p.X(), p.Y(), p.Z()); return Point<3> (p.X(), p.Y(), p.Z());
} }
double ComputeH (double kappa, const MeshingParameters & mparam) double ComputeH (double kappa, const MeshingParameters & mparam)
{ {
/*
double hret; double hret;
kappa *= mparam.curvaturesafety; kappa *= mparam.curvaturesafety;
@ -70,12 +60,14 @@ namespace netgen
if (mparam.maxh < hret) if (mparam.maxh < hret)
hret = mparam.maxh; hret = mparam.maxh;
return (hret); return hret;
*/
// return min(mparam.maxh, 1/kappa);
return (mparam.maxh*kappa < 1) ? mparam.maxh : 1/kappa;
} }
void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2, void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h, BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h,
const MeshingParameters & mparam) const MeshingParameters & mparam)
@ -316,6 +308,8 @@ namespace netgen
void OCCFindEdges (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam) void OCCFindEdges (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam)
{ {
static Timer t("OCCFindEdges"); RegionTimer r(t);
static Timer tsearch("OCCFindEdges - search point");
const char * savetask = multithread.task; const char * savetask = multithread.task;
multithread.task = "Edge meshing"; multithread.task = "Edge meshing";
@ -329,6 +323,7 @@ namespace netgen
double eps = 1e-6 * geom.GetBoundingBox().Diam(); double eps = 1e-6 * geom.GetBoundingBox().Diam();
tsearch.Start();
for (int i = 1; i <= nvertices; i++) for (int i = 1; i <= nvertices; i++)
{ {
gp_Pnt pnt = BRep_Tool::Pnt (TopoDS::Vertex(geom.vmap(i))); gp_Pnt pnt = BRep_Tool::Pnt (TopoDS::Vertex(geom.vmap(i)));
@ -346,6 +341,7 @@ namespace netgen
if (!exists) if (!exists)
mesh.AddPoint (mp); mesh.AddPoint (mp);
} }
tsearch.Stop();
(*testout) << "different vertices = " << mesh.GetNP() << endl; (*testout) << "different vertices = " << mesh.GetNP() << endl;
@ -504,6 +500,7 @@ namespace netgen
for (int i = 1; i <= mp.Size(); i++) for (int i = 1; i <= mp.Size(); i++)
{ {
bool exists = 0; bool exists = 0;
tsearch.Start();
int j; int j;
for (j = first_ep; j <= mesh.GetNP(); j++) for (j = first_ep; j <= mesh.GetNP(); j++)
if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps) if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps)
@ -511,6 +508,7 @@ namespace netgen
exists = 1; exists = 1;
break; break;
} }
tsearch.Stop();
if (exists) if (exists)
pnums[i] = j; pnums[i] = j;
@ -611,8 +609,10 @@ namespace netgen
void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend, void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend,
MeshingParameters & mparam) MeshingParameters & mparam)
{ {
int i, j, k; static Timer t("OCCMeshSurface"); RegionTimer r(t);
int changed;
// int i, j, k;
// int changed;
const char * savetask = multithread.task; const char * savetask = multithread.task;
multithread.task = "Surface meshing"; multithread.task = "Surface meshing";
@ -633,7 +633,7 @@ namespace netgen
int surfmesherror = 0; int surfmesherror = 0;
for (k = 1; k <= mesh.GetNFD(); k++) for (int k = 1; k <= mesh.GetNFD(); k++)
{ {
if(1==0 && !geom.fvispar[k-1].IsDrawable()) if(1==0 && !geom.fvispar[k-1].IsDrawable())
{ {
@ -655,7 +655,6 @@ namespace netgen
} }
*/ */
FaceDescriptor & fd = mesh.GetFaceDescriptor(k); FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
int oldnf = mesh.GetNSE(); int oldnf = mesh.GetNSE();
@ -664,7 +663,10 @@ namespace netgen
// int projecttype = PLANESPACE; // int projecttype = PLANESPACE;
static Timer tinit("init");
tinit.Start();
Meshing2OCCSurfaces meshing(TopoDS::Face(geom.fmap(k)), bb, projecttype, mparam); Meshing2OCCSurfaces meshing(TopoDS::Face(geom.fmap(k)), bb, projecttype, mparam);
tinit.Stop();
if (meshing.GetProjectionType() == PLANESPACE) if (meshing.GetProjectionType() == PLANESPACE)
PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (plane space projection)"); PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (plane space projection)");
@ -676,20 +678,20 @@ namespace netgen
// Meshing2OCCSurfaces meshing(f2, bb); // Meshing2OCCSurfaces meshing(f2, bb);
meshing.SetStartTime (starttime); meshing.SetStartTime (starttime);
//(*testout) << "Face " << k << endl << endl; //(*testout) << "Face " << k << endl << endl;
if (meshing.GetProjectionType() == PLANESPACE) if (meshing.GetProjectionType() == PLANESPACE)
{ {
static Timer t("MeshSurface: Find edges and points - Physical"); RegionTimer r(t);
int cntp = 0; int cntp = 0;
glob2loc = 0; glob2loc = 0;
for (i = 1; i <= mesh.GetNSeg(); i++) for (int i = 1; i <= mesh.GetNSeg(); i++)
{ {
Segment & seg = mesh.LineSegment(i); Segment & seg = mesh.LineSegment(i);
if (seg.si == k) if (seg.si == k)
{ {
for (j = 1; j <= 2; j++) for (int j = 1; j <= 2; j++)
{ {
int pi = (j == 1) ? seg[0] : seg[1]; int pi = (j == 1) ? seg[0] : seg[1];
if (!glob2loc.Get(pi)) if (!glob2loc.Get(pi))
@ -702,7 +704,7 @@ namespace netgen
} }
} }
for (i = 1; i <= mesh.GetNSeg(); i++) for (int i = 1; i <= mesh.GetNSeg(); i++)
{ {
Segment & seg = mesh.LineSegment(i); Segment & seg = mesh.LineSegment(i);
if (seg.si == k) if (seg.si == k)
@ -722,9 +724,11 @@ namespace netgen
} }
else else
{ {
static Timer t("MeshSurface: Find edges and points - Parameter"); RegionTimer r(t);
int cntp = 0; int cntp = 0;
for (i = 1; i <= mesh.GetNSeg(); i++) for (int i = 1; i <= mesh.GetNSeg(); i++)
if (mesh.LineSegment(i).si == k) if (mesh.LineSegment(i).si == k)
cntp+=2; cntp+=2;
@ -734,7 +738,7 @@ namespace netgen
gis.SetAllocSize (cntp); gis.SetAllocSize (cntp);
gis.SetSize (0); gis.SetSize (0);
for (i = 1; i <= mesh.GetNSeg(); i++) for (int i = 1; i <= mesh.GetNSeg(); i++)
{ {
Segment & seg = mesh.LineSegment(i); Segment & seg = mesh.LineSegment(i);
if (seg.si == k) if (seg.si == k)
@ -748,7 +752,7 @@ namespace netgen
int locpnum[2] = {0, 0}; int locpnum[2] = {0, 0};
for (j = 0; j < 2; j++) for (int j = 0; j < 2; j++)
{ {
PointGeomInfo gi = (j == 0) ? gi0 : gi1; PointGeomInfo gi = (j == 0) ? gi0 : gi1;
@ -798,6 +802,7 @@ namespace netgen
MESHING2_RESULT res; MESHING2_RESULT res;
try { try {
static Timer t("GenerateMesh"); RegionTimer reg(t);
res = meshing.GenerateMesh (mesh, mparam, maxh, k); res = meshing.GenerateMesh (mesh, mparam, maxh, k);
} }
@ -814,6 +819,7 @@ namespace netgen
} }
projecttype = PARAMETERSPACE; projecttype = PARAMETERSPACE;
static Timer t1("rest of loop"); RegionTimer reg1(t1);
if (res != MESHING2_OK) if (res != MESHING2_OK)
{ {
@ -846,14 +852,14 @@ namespace netgen
notrys = 1; notrys = 1;
for (i = oldnf+1; i <= mesh.GetNSE(); i++) for (int i = oldnf+1; i <= mesh.GetNSE(); i++)
mesh.SurfaceElement(i).SetIndex (k); mesh.SurfaceElement(i).SetIndex (k);
} }
// ofstream problemfile("occmesh.rep"); // ofstream problemfile("occmesh.rep");
// problemfile << "SURFACEMESHING" << endl << endl; // problemfile << "SURFACEMESHING" << endl << endl;
if (surfmesherror) if (surfmesherror)
{ {
@ -863,8 +869,8 @@ namespace netgen
if (geom.facemeshstatus[i-1] == -1) if (geom.facemeshstatus[i-1] == -1)
{ {
cout << "Face " << i << endl; cout << "Face " << i << endl;
// problemfile << "problem with face " << i << endl; // problemfile << "problem with face " << i << endl;
// problemfile << "vertices: " << endl; // problemfile << "vertices: " << endl;
TopExp_Explorer exp0,exp1,exp2; TopExp_Explorer exp0,exp1,exp2;
for ( exp0.Init(TopoDS::Face (geom.fmap(i)), TopAbs_WIRE); exp0.More(); exp0.Next() ) for ( exp0.Init(TopoDS::Face (geom.fmap(i)), TopAbs_WIRE); exp0.More(); exp0.Next() )
{ {
@ -876,22 +882,22 @@ namespace netgen
{ {
TopoDS_Vertex vertex = TopoDS::Vertex(exp2.Current()); TopoDS_Vertex vertex = TopoDS::Vertex(exp2.Current());
gp_Pnt point = BRep_Tool::Pnt(vertex); gp_Pnt point = BRep_Tool::Pnt(vertex);
// problemfile << point.X() << " " << point.Y() << " " << point.Z() << endl; // problemfile << point.X() << " " << point.Y() << " " << point.Z() << endl;
} }
} }
} }
// problemfile << endl; // problemfile << endl;
} }
cout << endl << endl; cout << endl << endl;
cout << "for more information open IGES/STEP Topology Explorer" << endl; cout << "for more information open IGES/STEP Topology Explorer" << endl;
// problemfile.close(); // problemfile.close();
throw NgException ("Problem in Surface mesh generation"); throw NgException ("Problem in Surface mesh generation");
} }
else else
{ {
// problemfile << "OK" << endl << endl; // problemfile << "OK" << endl << endl;
// problemfile.close(); // problemfile.close();
} }
@ -905,7 +911,7 @@ namespace netgen
static Timer timer_opt2d("Optimization 2D"); static Timer timer_opt2d("Optimization 2D");
timer_opt2d.Start(); timer_opt2d.Start();
for (k = 1; k <= mesh.GetNFD(); k++) for (int k = 1; k <= mesh.GetNFD(); k++)
{ {
// if (k != 42) continue; // if (k != 42) continue;
// if (k != 36) continue; // if (k != 36) continue;
@ -916,7 +922,7 @@ namespace netgen
FaceDescriptor & fd = mesh.GetFaceDescriptor(k); FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
PrintMessage (1, "Optimize Surface ", k); PrintMessage (1, "Optimize Surface ", k);
for (i = 1; i <= mparam.optsteps2d; i++) for (int i = 1; i <= mparam.optsteps2d; i++)
{ {
// (*testout) << "optstep " << i << endl; // (*testout) << "optstep " << i << endl;
if (multithread.terminate) return; if (multithread.terminate) return;
@ -1411,7 +1417,7 @@ namespace netgen
MESHING3_RESULT res = MeshVolume (mparam, *mesh); MESHING3_RESULT res = MeshVolume (mparam, *mesh);
/* /*
ofstream problemfile("occmesh.rep",ios_base::app); ofstream problemfile("occmesh.rep",ios_base::app);
problemfile << "VOLUMEMESHING" << endl << endl; problemfile << "VOLUMEMESHING" << endl << endl;
@ -1422,7 +1428,7 @@ namespace netgen
<< mesh->GetNE() << " elements" << endl << endl; << mesh->GetNE() << " elements" << endl << endl;
problemfile.close(); problemfile.close();
*/ */
if (res != MESHING3_OK) return TCL_ERROR; if (res != MESHING3_OK) return TCL_ERROR;

28
libsrc/occ/occgeom.hpp
View File

@ -171,10 +171,8 @@ namespace netgen
{ {
public: public:
Point<3> p0, p1; Point<3> p0, p1;
double Dist (Line l); double Dist (Line l);
double Length () { return (p1-p0).Length(); }
double Length ();
}; };
@ -250,22 +248,22 @@ namespace netgen
DLL_HEADER void BuildFMap(); DLL_HEADER void BuildFMap();
Box<3> GetBoundingBox() Box<3> GetBoundingBox() const
{ return boundingbox;} { return boundingbox; }
int NrSolids() int NrSolids() const
{ return somap.Extent();} { return somap.Extent(); }
// Philippose - 17/01/2009 // Philippose - 17/01/2009
// Total number of faces in the geometry // Total number of faces in the geometry
int NrFaces() int NrFaces() const
{ return fmap.Extent();} { return fmap.Extent(); }
void SetCenter() void SetCenter()
{ center = boundingbox.Center();} { center = boundingbox.Center(); }
Point<3> Center() Point<3> Center() const
{ return center;} { return center; }
void Project (int surfi, Point<3> & p) const; void Project (int surfi, Point<3> & p) const;
bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const; bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
@ -347,9 +345,7 @@ namespace netgen
// Returns the index of the currently selected face // Returns the index of the currently selected face
int SelectedFace() int SelectedFace()
{ {
int i; for(int i = 1; i <= fmap.Extent(); i++)
for(i = 1; i <= fmap.Extent(); i++)
{ {
if(face_sel_status[i-1]) if(face_sel_status[i-1])
{ {
@ -386,7 +382,7 @@ namespace netgen
DLL_HEADER void GetNotDrawableFaces (stringstream & str); DLL_HEADER void GetNotDrawableFaces (stringstream & str);
DLL_HEADER bool ErrorInSurfaceMeshing (); DLL_HEADER bool ErrorInSurfaceMeshing ();
// void WriteOCC_STL(char * filename); // void WriteOCC_STL(char * filename);
DLL_HEADER virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam); DLL_HEADER virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);

1
libsrc/occ/occmeshsurf.hpp
View File

@ -55,6 +55,7 @@ protected:
public: public:
OCCSurface (const TopoDS_Face & aface, int aprojecttype) OCCSurface (const TopoDS_Face & aface, int aprojecttype)
{ {
static Timer t("occurface ctor"); RegionTimer r(t);
topods_face = aface; topods_face = aface;
occface = BRep_Tool::Surface(topods_face); occface = BRep_Tool::Surface(topods_face);
orient = topods_face.Orientation(); orient = topods_face.Orientation();