start unify meshing with occ meshing

This commit is contained in:
Christopher Lackner 2019-10-02 17:20:13 +02:00 committed by Matthias Hochsteger
parent 86342534c1
commit bee097b153
11 changed files with 235 additions and 323 deletions

View File

@ -27,29 +27,59 @@ namespace netgen
int NetgenGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam)
{
if (!mesh) return 1;
multithread.percent = 0;
if (mparam.perfstepsstart <= MESHCONST_MESHVOLUME)
if(mparam.perfstepsstart <= MESHCONST_ANALYSE)
{
multithread.task = "Volume meshing";
MESHING3_RESULT res =
MeshVolume (mparam, *mesh);
if (res != MESHING3_OK) return 1;
if (multithread.terminate) return 0;
RemoveIllegalElements (*mesh);
if (multithread.terminate) return 0;
MeshQuality3d (*mesh);
if(!mesh)
mesh = make_shared<Mesh>();
mesh->geomtype = GetGeomType();
Analyse(*mesh, mparam);
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHVOLUME)
if(multithread.terminate || mparam.perfstepsend <= MESHCONST_ANALYSE)
return 0;
if(mparam.perfstepsstart <= MESHCONST_MESHEDGES)
FindEdges(*mesh, mparam);
if(multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHEDGES)
return 0;
if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE)
{
MeshSurface(*mesh, mparam);
mesh->CalcSurfacesOfNode();
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHSURFACE)
return 0;
if (mparam.perfstepsstart <= MESHCONST_OPTSURFACE)
OptimizeSurface(*mesh, mparam);
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_OPTSURFACE)
return 0;
if(mparam.perfstepsstart <= MESHCONST_MESHVOLUME)
{
multithread.task = "Volume meshing";
MESHING3_RESULT res = MeshVolume (mparam, *mesh);
if (res != MESHING3_OK) return 1;
if (multithread.terminate) return 0;
RemoveIllegalElements (*mesh);
if (multithread.terminate) return 0;
MeshQuality3d (*mesh);
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHVOLUME)
return 0;
if (mparam.perfstepsstart <= MESHCONST_OPTVOLUME)
{
@ -58,9 +88,9 @@ namespace netgen
OptimizeVolume (mparam, *mesh);
if (multithread.terminate) return 0;
}
FinalizeMesh(*mesh);
return 0;
}
}
const Refinement & NetgenGeometry :: GetRefinement () const

View File

@ -22,9 +22,16 @@ namespace netgen
virtual const Refinement & GetRefinement () const;
virtual void DoArchive(Archive&)
virtual void DoArchive(Archive&)
{ throw NgException("DoArchive not implemented for " + Demangle(typeid(*this).name())); }
virtual Mesh::GEOM_TYPE GetGeomType() const { return Mesh::NO_GEOM; }
virtual void Analyse(Mesh& mesh,
const MeshingParameters& mparam) {}
virtual void FindEdges(Mesh& mesh, const MeshingParameters& mparam) {}
virtual void MeshSurface(Mesh& mesh, const MeshingParameters& mparam) {}
virtual void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) {}
virtual void FinalizeMesh(Mesh& mesh) const {}
virtual void Save (string filename) const;
virtual void SaveToMeshFile (ostream & /* ost */) const { ; }
};

View File

@ -862,70 +862,81 @@ namespace netgen
// TODO: split also bad trigs, nut just illegal ones
if (mesh.LegalTrig(sel)) continue;
for (int j = 0; j < 3; j++)
// find longest edge
INDEX_2 edge;
double edge_len = 0;
PointIndex pi1, pi2, pi3, pi4;
PointGeomInfo gi1, gi2, gi3, gi4;
for(auto j : Range(1,4))
{
PointIndex pi1 = sel.PNumMod(j+2);
PointIndex pi2 = sel.PNumMod(j+3);
PointIndex pi3 = sel.PNumMod(j+1);
PointIndex pi4;
PointGeomInfo gi1 = sel.GeomInfoPiMod(j+2);
PointGeomInfo gi2 = sel.GeomInfoPiMod(j+3);
PointGeomInfo gi3 = sel.GeomInfoPiMod(j+1);
PointGeomInfo gi4;
if (mesh.IsSegment (pi1, pi2)) continue;
// get neighbor element
INDEX_2 ii2 (pi1, pi2);
ii2.Sort();
auto els = els_on_edge.Get(ii2);
SurfaceElementIndex other_i = get<0>(els);
if(other_i==sei) other_i = get<1>(els);
auto & other = mesh[other_i];
// find opposite point of neighbor element
for (int j = 0; j < 3; j++)
if(other[j]!=pi1 && other[j]!=pi2)
{
pi4 = other[j];
gi4 = other.GeomInfoPi(j);
break;
}
// split edge pi1,pi2
Point<3> p5;
PointIndex pi5;
PointGeomInfo gi5;
mesh.GetGeometry()->GetRefinement().PointBetween (mesh[pi1], mesh[pi2], 0.5,
faceindex,
gi1, gi2, p5, gi5);
pi5 = mesh.AddPoint(p5);
Element2d e1(3);
e1.SetIndex(faceindex);
e1={ {pi1,gi1}, {pi5,gi5}, {pi3,gi3} };
mesh.AddSurfaceElement( e1 );
Element2d e2(3);
e2.SetIndex(faceindex);
e2 ={ {pi5,gi5}, {pi2,gi2}, {pi3,gi3} };
mesh.AddSurfaceElement( e2 );
Element2d e3(3);
e3.SetIndex(faceindex);
e3 ={ {pi1,gi1}, {pi4,gi4}, {pi5,gi5} };
mesh.AddSurfaceElement( e3 );
Element2d e4(3);
e4.SetIndex(faceindex);
e4 ={ {pi4,gi4}, {pi2,gi2}, {pi5,gi5} };
mesh.AddSurfaceElement( e4 );
sel.Delete();
other.Delete();
auto test_pi1 = sel.PNumMod(j);
auto test_pi2 = sel.PNumMod(j+1);
if (mesh.IsSegment(test_pi1, test_pi2))
continue;
auto len = (mesh[test_pi2]-mesh[test_pi1]).Length();
if(len > edge_len)
{
edge = {test_pi1, test_pi2};
edge.Sort();
edge_len = len;
pi1 = test_pi1;
pi2 = test_pi2;
pi3 = sel.PNumMod(j+2);
gi1 = sel.GeomInfoPiMod(j);
gi2 = sel.GeomInfoPiMod(j+1);
gi3 = sel.GeomInfoPiMod(j+2);
}
}
if(!edge_len)
throw Exception("Couldn't find edge to split, something is wrong");
// get neighbor element
auto els = els_on_edge.Get(edge);
SurfaceElementIndex other_i = get<0>(els);
if(other_i==sei) other_i = get<1>(els);
auto & other = mesh[other_i];
// find opposite point of neighbor element
for (int j = 0; j < 3; j++)
if(other[j]!=pi1 && other[j]!=pi2)
{
pi4 = other[j];
gi4 = other.GeomInfoPi(j);
break;
}
// split edge pi1,pi2
Point<3> p5;
PointIndex pi5;
PointGeomInfo gi5;
mesh.GetGeometry()->GetRefinement().PointBetween (mesh[pi1], mesh[pi2], 0.5,
faceindex,
gi1, gi2, p5, gi5);
pi5 = mesh.AddPoint(p5);
Element2d e1(3);
e1.SetIndex(faceindex);
e1={ {pi1,gi1}, {pi5,gi5}, {pi3,gi3} };
mesh.AddSurfaceElement( e1 );
Element2d e2(3);
e2.SetIndex(faceindex);
e2 ={ {pi5,gi5}, {pi2,gi2}, {pi3,gi3} };
mesh.AddSurfaceElement( e2 );
Element2d e3(3);
e3.SetIndex(faceindex);
e3 ={ {pi1,gi1}, {pi4,gi4}, {pi5,gi5} };
mesh.AddSurfaceElement( e3 );
Element2d e4(3);
e4.SetIndex(faceindex);
e4 ={ {pi4,gi4}, {pi2,gi2}, {pi5,gi5} };
mesh.AddSurfaceElement( e4 );
sel.Delete();
other.Delete();
}
mesh.SetNextTimeStamp();

View File

@ -864,7 +864,7 @@ namespace netgen
///
friend class Meshing3;
// only for saving the geometry
enum GEOM_TYPE { NO_GEOM = 0, GEOM_2D = 1, GEOM_CSG = 10, GEOM_STL = 11, GEOM_OCC = 12, GEOM_ACIS = 13 };
GEOM_TYPE geomtype;

View File

@ -602,7 +602,7 @@ namespace netgen
void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh,
MeshingParameters & mparam)
const MeshingParameters & mparam)
{
static Timer t("OCCMeshSurface"); RegionTimer r(t);
@ -796,7 +796,7 @@ namespace netgen
// Philippose - 15/01/2009
double maxh = geom.face_maxh[k-1];
//double maxh = mparam.maxh;
mparam.checkoverlap = 0;
// mparam.checkoverlap = 0;
// int noldpoints = mesh->GetNP();
int noldsurfel = mesh.GetNSE();
@ -916,7 +916,7 @@ namespace netgen
}
void OCCOptimizeSurface(OCCGeometry & geom, Mesh & mesh,
MeshingParameters & mparam)
const MeshingParameters & mparam)
{
const char * savetask = multithread.task;
multithread.task = "Optimizing surface";
@ -991,7 +991,7 @@ namespace netgen
void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh,
void OCCSetLocalMeshSize(const OCCGeometry & geom, Mesh & mesh,
const MeshingParameters & mparam, const OCCParameters& occparam)
{
static Timer t1("OCCSetLocalMeshSize");
@ -1279,197 +1279,6 @@ namespace netgen
mesh.LoadLocalMeshSize (mparam.meshsizefilename);
}
int OCCGenerateMesh (OCCGeometry & geom, shared_ptr<Mesh> & mesh, MeshingParameters & mparam,
const OCCParameters& occparam)
{
multithread.percent = 0;
if (mparam.perfstepsstart <= MESHCONST_ANALYSE)
{
if(mesh.get() == nullptr)
mesh = make_shared<Mesh>();
mesh->geomtype = Mesh::GEOM_OCC;
OCCSetLocalMeshSize(geom,*mesh, mparam, occparam);
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_ANALYSE)
return TCL_OK;
if (mparam.perfstepsstart <= MESHCONST_MESHEDGES)
{
OCCFindEdges (geom, *mesh, mparam);
/*
cout << "Removing redundant points" << endl;
int i, j;
int np = mesh->GetNP();
NgArray<int> equalto;
equalto.SetSize (np);
equalto = 0;
for (i = 1; i <= np; i++)
{
for (j = i+1; j <= np; j++)
{
if (!equalto[j-1] && (Dist2 (mesh->Point(i), mesh->Point(j)) < 1e-12))
equalto[j-1] = i;
}
}
for (i = 1; i <= np; i++)
if (equalto[i-1])
{
cout << "Point " << i << " is equal to Point " << equalto[i-1] << endl;
for (j = 1; j <= mesh->GetNSeg(); j++)
{
Segment & seg = mesh->LineSegment(j);
if (seg[0] == i) seg[0] = equalto[i-1];
if (seg[1] == i) seg[1] = equalto[i-1];
}
}
cout << "Removing degenerated segments" << endl;
for (j = 1; j <= mesh->GetNSeg(); j++)
{
Segment & seg = mesh->LineSegment(j);
if (seg[0] == seg[1])
{
mesh->DeleteSegment(j);
cout << "Deleting Segment " << j << endl;
}
}
mesh->Compress();
*/
/*
for (int i = 1; i <= geom.fmap.Extent(); i++)
{
Handle(Geom_Surface) hf1 =
BRep_Tool::Surface(TopoDS::Face(geom.fmap(i)));
for (int j = i+1; j <= geom.fmap.Extent(); j++)
{
Handle(Geom_Surface) hf2 =
BRep_Tool::Surface(TopoDS::Face(geom.fmap(j)));
if (hf1 == hf2) cout << "face " << i << " and face " << j << " lie on same surface" << endl;
}
}
*/
#ifdef LOG_STREAM
(*logout) << "Edges meshed" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
#endif
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHEDGES)
return TCL_OK;
if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE)
{
OCCMeshSurface (geom, *mesh, mparam);
if (multithread.terminate) return TCL_OK;
#ifdef LOG_STREAM
(*logout) << "Surfaces meshed" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
#endif
#ifdef STAT_STREAM
(*statout) << mesh->GetNSeg() << " & "
<< mesh->GetNSE() << " & - &"
<< GetTime() << " & " << endl;
#endif
// MeshQuality2d (*mesh);
mesh->CalcSurfacesOfNode();
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHSURFACE)
return TCL_OK;
if (mparam.perfstepsstart <= MESHCONST_OPTSURFACE)
{
OCCOptimizeSurface(geom, *mesh, mparam);
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_OPTSURFACE)
return TCL_OK;
if (mparam.perfstepsstart <= MESHCONST_MESHVOLUME)
{
multithread.task = "Volume meshing";
MESHING3_RESULT res = MeshVolume (mparam, *mesh);
if (res != MESHING3_OK) return TCL_ERROR;
if (multithread.terminate) return TCL_OK;
RemoveIllegalElements (*mesh);
if (multithread.terminate) return TCL_OK;
MeshQuality3d (*mesh);
#ifdef STAT_STREAM
(*statout) << GetTime() << " & ";
#endif
#ifdef LOG_STREAM
(*logout) << "Volume meshed" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
#endif
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHVOLUME)
return TCL_OK;
if (mparam.perfstepsstart <= MESHCONST_OPTVOLUME)
{
multithread.task = "Volume optimization";
OptimizeVolume (mparam, *mesh);
if (multithread.terminate) return TCL_OK;
#ifdef STAT_STREAM
(*statout) << GetTime() << " & "
<< mesh->GetNE() << " & "
<< mesh->GetNP() << " " << '\\' << '\\' << " \\" << "hline" << endl;
#endif
#ifdef LOG_STREAM
(*logout) << "Volume optimized" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
#endif
// cout << "Optimization complete" << endl;
}
/*
(*testout) << "NP: " << mesh->GetNP() << endl;
for (int i = 1; i <= mesh->GetNP(); i++)
(*testout) << mesh->Point(i) << endl;
(*testout) << endl << "NSegments: " << mesh->GetNSeg() << endl;
for (int i = 1; i <= mesh->GetNSeg(); i++)
(*testout) << mesh->LineSegment(i) << endl;
*/
for (int i = 0; i < mesh->GetNDomains(); i++)
if (geom.snames.Size())
mesh->SetMaterial (i+1, geom.snames[i]);
return TCL_OK;
}
}
#endif

View File

@ -74,6 +74,36 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
}
}
void OCCGeometry :: Analyse(Mesh& mesh,
const MeshingParameters& mparam)
{
OCCSetLocalMeshSize(*this, mesh, mparam, occparam);
}
void OCCGeometry :: FindEdges(Mesh& mesh,
const MeshingParameters& mparam)
{
OCCFindEdges(*this, mesh, mparam);
}
void OCCGeometry :: MeshSurface(Mesh& mesh,
const MeshingParameters& mparam)
{
OCCMeshSurface(*this, mesh, mparam);
}
void OCCGeometry :: OptimizeSurface(Mesh& mesh,
const MeshingParameters& mparam)
{
OCCOptimizeSurface(*this, mesh, mparam);
}
void OCCGeometry :: FinalizeMesh(Mesh& mesh) const
{
for (int i = 0; i < mesh.GetNDomains(); i++)
if (snames.Size())
mesh.SetMaterial (i+1, snames[i]);
}
void OCCGeometry :: PrintNrShapes ()
{
@ -1703,10 +1733,10 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
DLL_HEADER extern OCCParameters occparam;
OCCParameters occparam;
int OCCGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam)
{
return OCCGenerateMesh (*this, mesh, mparam, occparam);
}
// int OCCGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam)
// {
// return OCCGenerateMesh (*this, mesh, mparam, occparam);
// }
}

View File

@ -183,12 +183,33 @@ namespace netgen
return a00*a11*a22 + a01*a12*a20 + a10*a21*a02 - a20*a11*a02 - a10*a01*a22 - a21*a12*a00;
}
class DLL_HEADER OCCParameters
{
public:
/// Factor for meshing close edges
double resthcloseedgefac = 2.;
/// Enable / Disable detection of close edges
int resthcloseedgeenable = true;
/// Minimum edge length to be used for dividing edges to mesh points
double resthminedgelen = 0.001;
/// Enable / Disable use of the minimum edge length (by default use 1e-4)
int resthminedgelenenable = true;
/*!
Dump all the OpenCascade specific meshing parameters
to console
*/
void Print (ostream & ost) const;
};
class OCCGeometry : public NetgenGeometry
{
Point<3> center;
OCCParameters occparam;
public:
TopoDS_Shape shape;
@ -239,11 +260,26 @@ namespace netgen
emap.Clear();
vmap.Clear();
}
Mesh::GEOM_TYPE GetGeomType() const override
{ return Mesh::GEOM_OCC; }
void SetOCCParameters(const OCCParameters& par)
{ cout << "set occ pars to = " << par.resthcloseedgefac << endl; occparam = par; }
void Analyse(Mesh& mesh,
const MeshingParameters& mparam) override;
void FindEdges(Mesh& mesh,
const MeshingParameters& mparam) override;
void MeshSurface(Mesh& mesh,
const MeshingParameters& mparam) override;
void OptimizeSurface(Mesh& mesh,
const MeshingParameters& mparam) override;
void FinalizeMesh(Mesh& mesh) const override;
DLL_HEADER void Save (string filename) const override;
DLL_HEADER virtual void Save (string filename) const;
void DoArchive(Archive& ar);
void DoArchive(Archive& ar) override;
DLL_HEADER void BuildFMap();
@ -391,37 +427,9 @@ namespace netgen
// 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);
DLL_HEADER virtual const Refinement & GetRefinement () const;
};
class DLL_HEADER OCCParameters
{
public:
/// Factor for meshing close edges
double resthcloseedgefac = 2.;
/// Enable / Disable detection of close edges
int resthcloseedgeenable = true;
/// Minimum edge length to be used for dividing edges to mesh points
double resthminedgelen = 0.001;
/// Enable / Disable use of the minimum edge length (by default use 1e-4)
int resthminedgelenenable = true;
/*!
Dump all the OpenCascade specific meshing parameters
to console
*/
void Print (ostream & ost) const;
DLL_HEADER const Refinement & GetRefinement () const override;
};
@ -434,15 +442,12 @@ namespace netgen
// Philippose - 31.09.2009
// External access to the mesh generation functions within the OCC
// subsystem (Not sure if this is the best way to implement this....!!)
DLL_HEADER extern int OCCGenerateMesh (OCCGeometry & occgeometry, shared_ptr<Mesh> & mesh,
MeshingParameters & mparam, const OCCParameters& occparam);
DLL_HEADER extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam,
DLL_HEADER extern void OCCSetLocalMeshSize(const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam,
const OCCParameters& occparam);
DLL_HEADER extern void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, MeshingParameters & mparam);
DLL_HEADER extern void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
DLL_HEADER extern void OCCOptimizeSurface (OCCGeometry & geom, Mesh & mesh, MeshingParameters & mparam);
DLL_HEADER extern void OCCOptimizeSurface (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
DLL_HEADER extern void OCCFindEdges (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
}

View File

@ -53,6 +53,8 @@ namespace netgen
atof (Tcl_GetVar (interp, "::stloptions.resthminedgelen", 0));
occparam.resthminedgelenenable =
atoi (Tcl_GetVar (interp, "::stloptions.resthminedgelenenable", 0));
if(auto geo = dynamic_pointer_cast<OCCGeometry>(ng_geometry); geo)
geo->SetOCCParameters(occparam);
}
};

View File

@ -183,11 +183,12 @@ DLL_HEADER void ExportNgOCC(py::module &m)
CreateOCCParametersFromKwargs(occparam, kwargs);
CreateMPfromKwargs(mp, kwargs);
}
geo->SetOCCParameters(occparam);
auto mesh = make_shared<Mesh>();
SetGlobalMesh(mesh);
geo->GenerateMesh(mesh, mp);
mesh->SetGeometry(geo);
SetGlobalMesh(mesh);
ng_geometry = geo;
OCCGenerateMesh(*geo, mesh, mp, occparam);
return mesh;
}, py::arg("mp") = nullptr,
py::call_guard<py::gil_scoped_release>(),

View File

@ -867,6 +867,22 @@ void STLSurfaceOptimization (STLGeometry & geom,
break;
}
}
while(mesh.CheckOverlappingBoundary())
{
for(const auto & el : mesh.SurfaceElements())
{
if(el.BadElement())
{
cout << "Restrict localh at el nr " << el << endl;
for(const auto& p : el.PNums())
{
const auto& pnt = mesh[p];
mesh.RestrictLocalH(pnt, 0.5*mesh.GetH(pnt));
}
}
}
optmesh.SplitImprove(mesh);
}
//(*testout) << "optimize, after, step = " << meshparam.optimize2d[j-1] << mesh.Point (3679) << endl;
}

View File

@ -537,6 +537,7 @@ namespace netgen
// delete ng_geometry;
// ng_geometry = hgeom;
ng_geometry = shared_ptr<NetgenGeometry> (hgeom);
geometryregister[i]->SetParameters(interp);
mesh.reset();
return TCL_OK;