ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-01-12 22:20:35 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'test_refactoring_meshing_design' into 'master'

Browse Source 
Refactoring of surface meshing classes

See merge request jschoeberl/netgen!279
This commit is contained in:
Joachim Schöberl 2019-10-21 08:21:12 +00:00
commit 124ee905b2
48 changed files with 1721 additions and 2350 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

87
libsrc/csg/csgeom.cpp
View File

@ -72,6 +72,84 @@ namespace netgen
Clean();
}
void CSGeometry :: ProjectPoint(int surfind, Point<3> & p) const
{
Point<3> hp = p;
GetSurface(surfind)->Project (hp);
p = hp;
}
void CSGeometry :: ProjectPointEdge(int surfind, INDEX surfind2,
Point<3> & p) const
{
Point<3> hp = p;
ProjectToEdge (GetSurface(surfind),
GetSurface(surfind2), hp);
p = hp;
}
Vec<3> CSGeometry :: GetNormal(int surfind, const Point<3> & p) const
{
Vec<3> hn;
GetSurface(surfind)->CalcGradient(p, hn);
hn.Normalize();
return hn;
}
void CSGeometry ::
PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
Point<3> hnewp;
hnewp = p1+secpoint*(p2-p1);
if (surfi != -1)
{
GetSurface (surfi) -> Project (hnewp);
newgi.trignum = 1;
}
newp = hnewp;
}
void CSGeometry :: PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
Point<3> hnewp = p1+secpoint*(p2-p1);
//(*testout) << "hnewp " << hnewp << " s1 " << surfi1 << " s2 " << surfi2 << endl;
if (surfi1 != -1 && surfi2 != -1 && surfi1 != surfi2)
{
netgen::ProjectToEdge (GetSurface(surfi1),
GetSurface(surfi2),
hnewp);
// (*testout) << "Pointbetween, newp = " << hnewp << endl
// << ", err = " << sqrt (sqr (hnewp(0))+ sqr(hnewp(1)) + sqr (hnewp(2))) - 1 << endl;
newgi.edgenr = 1;
//(*testout) << "hnewp (a1) " << hnewp << endl;
}
else if (surfi1 != -1)
{
GetSurface (surfi1) -> Project (hnewp);
//(*testout) << "hnewp (a2) " << hnewp << endl;
}
newp = hnewp;
};
Vec<3> CSGeometry :: GetTangent(const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const
{
Vec<3> n1 = GetSurface (surfi1)->GetNormalVector (p);
Vec<3> n2 = GetSurface (surfi2)->GetNormalVector (p);
Vec<3> tau = Cross (n1, n2).Normalize();
return tau;
}
void CSGeometry :: Clean ()
{
@ -137,15 +215,6 @@ namespace netgen
return CSGGenerateMesh (*this, mesh, mparam);
}
const Refinement & CSGeometry :: GetRefinement () const
{
// cout << "get CSGeometry - Refinement" << endl;
// should become class variables
RefinementSurfaces * ref = new RefinementSurfaces(*this);
ref -> Set2dOptimizer(new MeshOptimize2dSurfaces(*this));
return *ref;
}
class WritePrimitivesIt : public SolidIterator
{
ostream & ost;

20
libsrc/csg/csgeom.hpp
View File

@ -188,6 +188,24 @@ namespace netgen
virtual void SaveToMeshFile (ostream & ost) const override;
void ProjectPoint(INDEX surfind, Point<3> & p) const override;
void ProjectPointEdge(INDEX surfind, INDEX surfind2, Point<3> & p) const override;
Vec<3> GetNormal(int surfind, const Point<3> & p) const override;
void PointBetween(const Point<3> & p1, const Point<3> & p2,
double secpoint, int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
Vec<3> GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const override;
int GetChangeVal() { return changeval; }
void Change() { changeval++; }
@ -348,8 +366,6 @@ namespace netgen
virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam) override;
virtual const Refinement & GetRefinement () const override;
void AddSplineSurface (shared_ptr<SplineSurface> ss) { spline_surfaces.Append(ss); }
};

16
libsrc/csg/genmesh.cpp
View File

@ -523,48 +523,48 @@ namespace netgen
if (multithread.terminate) return;
{
MeshOptimize2dSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.EdgeSwapping (mesh, (i > mparam.optsteps2d/2));
meshopt.EdgeSwapping (i > mparam.optsteps2d/2);
}
if (multithread.terminate) return;
{
// mesh.CalcSurfacesOfNode();
MeshOptimize2dSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.ImproveMesh (mesh, mparam);
meshopt.ImproveMesh(mparam);
}
{
MeshOptimize2dSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.CombineImprove (mesh);
meshopt.CombineImprove();
// mesh.CalcSurfacesOfNode();
}
if (multithread.terminate) return;
{
MeshOptimize2dSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.ImproveMesh (mesh, mparam);
meshopt.ImproveMesh(mparam);
}
}
}

143
libsrc/csg/meshsurf.cpp
View File

@ -59,147 +59,4 @@ double Meshing2Surfaces :: CalcLocalH (const Point<3> & p, double gh) const
return loch;
*/
}
MeshOptimize2dSurfaces :: MeshOptimize2dSurfaces (const CSGeometry & ageometry)
: MeshOptimize2d(), geometry(ageometry)
{
;
}
void MeshOptimize2dSurfaces :: ProjectPoint (INDEX surfind, Point<3> & p) const
{
Point<3> hp = p;
geometry.GetSurface(surfind)->Project (hp);
p = hp;
}
void MeshOptimize2dSurfaces :: ProjectPoint2 (INDEX surfind, INDEX surfind2,
Point<3> & p) const
{
Point<3> hp = p;
ProjectToEdge ( geometry.GetSurface(surfind),
geometry.GetSurface(surfind2), hp);
p = hp;
}
void MeshOptimize2dSurfaces ::
GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const
{
Vec<3> hn = n;
geometry.GetSurface(surfind)->CalcGradient (p, hn);
hn.Normalize();
n = hn;
/*
if (geometry.GetSurface(surfind)->Inverse())
n *= -1;
*/
}
RefinementSurfaces :: RefinementSurfaces (const CSGeometry & ageometry)
: Refinement(), geometry(ageometry)
{
if(geometry.GetNSurf() == 0)
*testout << endl
<< "WARNING: Initializing 2D refinement with 0-surface geometry" << endl
<< "==========================================================" << endl
<< endl << endl;
}
RefinementSurfaces :: ~RefinementSurfaces ()
{
;
}
void RefinementSurfaces ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
Point<3> hnewp;
hnewp = p1+secpoint*(p2-p1);
if (surfi != -1)
{
geometry.GetSurface (surfi) -> Project (hnewp);
newgi.trignum = 1;
}
newp = hnewp;
}
void RefinementSurfaces ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
Point<3> hnewp = p1+secpoint*(p2-p1);
//(*testout) << "hnewp " << hnewp << " s1 " << surfi1 << " s2 " << surfi2 << endl;
if (surfi1 != -1 && surfi2 != -1 && surfi1 != surfi2)
{
netgen::ProjectToEdge (geometry.GetSurface(surfi1),
geometry.GetSurface(surfi2),
hnewp);
// (*testout) << "Pointbetween, newp = " << hnewp << endl
// << ", err = " << sqrt (sqr (hnewp(0))+ sqr(hnewp(1)) + sqr (hnewp(2))) - 1 << endl;
newgi.edgenr = 1;
//(*testout) << "hnewp (a1) " << hnewp << endl;
}
else if (surfi1 != -1)
{
geometry.GetSurface (surfi1) -> Project (hnewp);
//(*testout) << "hnewp (a2) " << hnewp << endl;
}
newp = hnewp;
};
Vec<3> RefinementSurfaces :: GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const
{
Vec<3> n1 = geometry.GetSurface (surfi1)->GetNormalVector (p);
Vec<3> n2 = geometry.GetSurface (surfi2)->GetNormalVector (p);
Vec<3> tau = Cross (n1, n2).Normalize();
return tau;
}
Vec<3> RefinementSurfaces :: GetNormal (const Point<3> & p, int surfi1,
const PointGeomInfo & gi) const
{
return geometry.GetSurface (surfi1)->GetNormalVector (p);
}
void RefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi) const
{
if (surfi != -1)
geometry.GetSurface (surfi) -> Project (p);
};
void RefinementSurfaces :: ProjectToEdge (Point<3> & p, int surfi1, int surfi2, const EdgePointGeomInfo & egi) const
{
netgen::ProjectToEdge (geometry.GetSurface(surfi1),
geometry.GetSurface(surfi2),
p);
}
}

56
libsrc/csg/meshsurf.hpp
View File

@ -38,62 +38,6 @@ namespace netgen
///
double CalcLocalH(const Point<3> & p, double gh) const override;
};
///
class MeshOptimize2dSurfaces : public MeshOptimize2d
{
///
const CSGeometry & geometry;
public:
///
MeshOptimize2dSurfaces (const CSGeometry & ageometry);
///
virtual void ProjectPoint (INDEX surfind, Point<3> & p) const override;
///
virtual void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point<3> & p) const override;
///
virtual void GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const override;
};
class RefinementSurfaces : public Refinement
{
const CSGeometry & geometry;
public:
RefinementSurfaces (const CSGeometry & ageometry);
virtual ~RefinementSurfaces ();
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
virtual Vec<3> GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const override;
virtual Vec<3> GetNormal (const Point<3> & p, int surfi1,
const PointGeomInfo & gi) const override;
virtual void ProjectToSurface (Point<3> & p, int surfi) const override;
virtual void ProjectToEdge (Point<3> & p, int surfi1, int surfi2, const EdgePointGeomInfo & egi) const override;
};
}
#endif

4
libsrc/geom2d/CMakeLists.txt
View File

@ -1,5 +1,5 @@
add_definitions(-DNGLIB_EXPORTS)
add_library(geom2d ${NG_LIB_TYPE} genmesh2d.cpp geom2dmesh.cpp geometry2d.cpp python_geom2d.cpp )
add_library(geom2d ${NG_LIB_TYPE} genmesh2d.cpp geometry2d.cpp python_geom2d.cpp )
if(APPLE)
set_target_properties( geom2d PROPERTIES SUFFIX ".so")
endif(APPLE)
@ -18,7 +18,7 @@ if(USE_GUI)
endif(USE_GUI)
install(FILES
geom2dmesh.hpp geometry2d.hpp spline2d.hpp
geometry2d.hpp spline2d.hpp
vsgeom2d.hpp
DESTINATION ${NG_INSTALL_DIR_INCLUDE}/geom2d COMPONENT netgen_devel
)

119
libsrc/geom2d/geom2dmesh.cpp
View File

@ -1,119 +0,0 @@
#include <meshing.hpp>
#include <geometry2d.hpp>
namespace netgen
{
Refinement2d :: Refinement2d (const SplineGeometry2d & ageometry)
: Refinement(), geometry(ageometry)
{
;
}
Refinement2d :: ~Refinement2d ()
{
;
}
void Refinement2d ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
newp = p1+secpoint*(p2-p1);
newgi.trignum = 1;
}
void Refinement2d ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
Point<2> p2d;
double newdist;
auto spline = geometry.GetSplines().Get(ap1.edgenr);
if( (ap1.dist == 0.0) && (ap2.dist == 0.0) )
{
// used for manually generated meshes
const SplineSeg3<2> * ss3;
const LineSeg<2> * ls;
auto ext = dynamic_cast<const SplineSegExt *>(spline);
if(ext)
{
ss3 = dynamic_cast<const SplineSeg3<2> *>(ext->seg);
ls = dynamic_cast<const LineSeg<2> *>(ext->seg);
}
else
{
ss3 = dynamic_cast<const SplineSeg3<2> *>(spline);
ls = dynamic_cast<const LineSeg<2> *>(spline);
}
Point<2> p12d(p1(0),p1(1)), p22d(p2(0),p2(1));
Point<2> p1_proj(0.0,0.0), p2_proj(0.0,0.0);
double t1_proj = 0.0;
double t2_proj = 0.0;
if(ss3)
{
ss3->Project(p12d,p1_proj,t1_proj);
ss3->Project(p22d,p2_proj,t2_proj);
}
else if(ls)
{
ls->Project(p12d,p1_proj,t1_proj);
ls->Project(p22d,p2_proj,t2_proj);
}
p2d = spline->GetPoint (((1-secpoint)*t1_proj+secpoint*t2_proj));
newdist = (1-secpoint)*t1_proj+secpoint*t2_proj;
}
else
{
p2d = spline->GetPoint (((1-secpoint)*ap1.dist+secpoint*ap2.dist));
newdist = (1-secpoint)*ap1.dist+secpoint*ap2.dist;
}
// (*testout) << "refine 2d line, ap1.dist, ap2.dist = " << ap1.dist << ", " << ap2.dist << endl;
// (*testout) << "p1, p2 = " << p1 << p2 << ", newp = " << p2d << endl;
newp = Point3d (p2d(0), p2d(1), 0);
newgi.edgenr = ap1.edgenr;
newgi.dist = newdist;
};
Vec<3> Refinement2d :: GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const
{
Vec<2> t2d = geometry.GetSplines().Get(ap1.edgenr) -> GetTangent(ap1.dist);
return Vec<3> (t2d(0), t2d(1), 0);
}
Vec<3> Refinement2d :: GetNormal (const Point<3> & p, int surfi1,
const PointGeomInfo & gi) const
{
return Vec<3> (0,0,1);
}
void Refinement2d :: ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & /* gi */) const
{
p(2) = 0;
}
void Refinement2d :: ProjectToEdge (Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & egi) const
{
Point<2> p2d (p(0), p(1)), pp;
double t;
geometry.GetSplines().Get(egi.edgenr) -> Project (p2d, pp, t);
p = Point<3> (pp(0), pp(1), 0);
}
}

52
libsrc/geom2d/geom2dmesh.hpp
View File

@ -1,52 +0,0 @@
#ifndef FILE_GEOM2DMESH
#define FILE_GEOM2DMESH
/**************************************************************************/
/* File: geom2dmesh.hh */
/* Author: Joachim Schoeberl */
/* Date: 22. Jan. 01 */
/**************************************************************************/
namespace netgen
{
class Refinement2d : public Refinement
{
const class SplineGeometry2d & geometry;
public:
Refinement2d (const class SplineGeometry2d & ageometry);
virtual ~Refinement2d ();
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
virtual Vec<3> GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const override;
virtual Vec<3> GetNormal (const Point<3> & p, int surfi1,
const PointGeomInfo & gi) const override;
virtual void ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & /* gi */) const override;
virtual void ProjectToEdge (Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & egi) const override;
};
}
#endif

78
libsrc/geom2d/geometry2d.cpp
View File

@ -20,6 +20,76 @@ namespace netgen
delete [] materials[i];
}
void SplineGeometry2d :: PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
Point<2> p2d;
double newdist;
auto spline = GetSplines().Get(ap1.edgenr);
if( (ap1.dist == 0.0) && (ap2.dist == 0.0) )
{
// used for manually generated meshes
const SplineSeg3<2> * ss3;
const LineSeg<2> * ls;
auto ext = dynamic_cast<const SplineSegExt *>(spline);
if(ext)
{
ss3 = dynamic_cast<const SplineSeg3<2> *>(ext->seg);
ls = dynamic_cast<const LineSeg<2> *>(ext->seg);
}
else
{
ss3 = dynamic_cast<const SplineSeg3<2> *>(spline);
ls = dynamic_cast<const LineSeg<2> *>(spline);
}
Point<2> p12d(p1(0),p1(1)), p22d(p2(0),p2(1));
Point<2> p1_proj(0.0,0.0), p2_proj(0.0,0.0);
double t1_proj = 0.0;
double t2_proj = 0.0;
if(ss3)
{
ss3->Project(p12d,p1_proj,t1_proj);
ss3->Project(p22d,p2_proj,t2_proj);
}
else if(ls)
{
ls->Project(p12d,p1_proj,t1_proj);
ls->Project(p22d,p2_proj,t2_proj);
}
p2d = spline->GetPoint (((1-secpoint)*t1_proj+secpoint*t2_proj));
newdist = (1-secpoint)*t1_proj+secpoint*t2_proj;
}
else
{
p2d = spline->GetPoint (((1-secpoint)*ap1.dist+secpoint*ap2.dist));
newdist = (1-secpoint)*ap1.dist+secpoint*ap2.dist;
}
// (*testout) << "refine 2d line, ap1.dist, ap2.dist = " << ap1.dist << ", " << ap2.dist << endl;
// (*testout) << "p1, p2 = " << p1 << p2 << ", newp = " << p2d << endl;
newp = Point3d (p2d(0), p2d(1), 0);
newgi.edgenr = ap1.edgenr;
newgi.dist = newdist;
};
Vec<3> SplineGeometry2d :: GetTangent(const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const
{
Vec<2> t2d = GetSplines().Get(ap1.edgenr) -> GetTangent(ap1.dist);
return Vec<3> (t2d(0), t2d(1), 0);
}
Vec<3> SplineGeometry2d :: GetNormal(int surfi1, const Point<3> & p,
const PointGeomInfo & gi) const
{
return Vec<3> (0,0,1);
}
void SplineGeometry2d :: Load (const char * filename)
{
@ -992,14 +1062,6 @@ namespace netgen
return 0;
}
Refinement & SplineGeometry2d :: GetRefinement () const
{
return * new Refinement2d (*this);
}
class SplineGeometryRegister : public GeometryRegister
{
public:

30
libsrc/geom2d/geometry2d.hpp
View File

@ -13,7 +13,6 @@
// #include "../gprim/spline.hpp"
// #include "../gprim/splinegeometry.hpp"
#include "geom2dmesh.hpp"
namespace netgen
{
@ -151,12 +150,34 @@ namespace netgen
void TestComment ( ifstream & infile ) ;
void DoArchive(Archive& ar)
void DoArchive(Archive& ar) override
{
SplineGeometry<2>::DoArchive(ar);
ar & materials & maxh & quadmeshing & tensormeshing & layer & bcnames & elto0;
}
void PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override
{
newp = p1+secpoint*(p2-p1);
newgi.trignum = 1;
}
void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
Vec<3> GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const override;
Vec<3> GetNormal(int surfi1, const Point<3> & p,
const PointGeomInfo & gi) const override;
const SplineSegExt & GetSpline (const int i) const
{
return dynamic_cast<const SplineSegExt&> (*splines[i]);
@ -168,7 +189,7 @@ namespace netgen
}
DLL_HEADER virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
DLL_HEADER int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam) override;
void PartitionBoundary (MeshingParameters & mp, double h, Mesh & mesh2d);
@ -205,9 +226,6 @@ namespace netgen
int AddBCName (string name);
string * BCNamePtr ( const int bcnr );
DLL_HEADER virtual Refinement & GetRefinement () const;
};
}

107
libsrc/meshing/basegeom.cpp
View File

@ -10,6 +10,40 @@ namespace netgen
GeometryRegister :: ~GeometryRegister()
{ ; }
void NetgenGeometry :: OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam)
{
const auto savetask = multithread.task;
multithread.task = "Optimizing surface";
static Timer timer_opt2d("Optimization 2D");
RegionTimer reg(timer_opt2d);
auto meshopt = MeshOptimize2d(mesh);
for(auto i : Range(mparam.optsteps2d))
{
PrintMessage(2, "Optimization step ", i);
for(auto optstep : mparam.optimize2d)
{
switch(optstep)
{
case 's':
meshopt.EdgeSwapping(0);
break;
case 'S':
meshopt.EdgeSwapping(1);
break;
case 'm':
meshopt.ImproveMesh(mparam);
break;
case 'c':
meshopt.CombineImprove();
break;
}
}
}
mesh.CalcSurfacesOfNode();
mesh.Compress();
multithread.task = savetask;
}
shared_ptr<NetgenGeometry> GeometryRegisterArray :: LoadFromMeshFile (istream & ist) const
{
@ -27,26 +61,56 @@ 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_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;
@ -54,21 +118,14 @@ namespace netgen
if (mparam.perfstepsstart <= MESHCONST_OPTVOLUME)
{
multithread.task = "Volume optimization";
OptimizeVolume (mparam, *mesh);
if (multithread.terminate) return 0;
}
FinalizeMesh(*mesh);
return 0;
}
const Refinement & NetgenGeometry :: GetRefinement () const
{
return *new Refinement;;
}
void NetgenGeometry :: Save (string filename) const
{
throw NgException("Cannot save geometry - no geometry available");

69
libsrc/meshing/basegeom.hpp
View File

@ -15,16 +15,81 @@ namespace netgen
class DLL_HEADER NetgenGeometry
{
unique_ptr<Refinement> ref;
public:
NetgenGeometry()
{
ref = make_unique<Refinement>(*this);
}
virtual ~NetgenGeometry () { ; }
virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
virtual const Refinement & GetRefinement () const;
virtual const Refinement & GetRefinement () const
{
return *ref;
}
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 ProjectPoint (int surfind, Point<3> & p) const
{ }
virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p) const { }
virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo& gi) const
{ ProjectPointEdge(surfind, surfind2, p); }
virtual bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const {return false;}
virtual bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const
{
ProjectPoint(surfind, p);
return CalcPointGeomInfo(surfind, gi, p);
}
virtual Vec<3> GetNormal(int surfind, const Point<3> & p) const
{ return {0.,0.,1.}; }
virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const
{ return GetNormal(surfind, p); }
[[deprecated]]
void GetNormal(int surfind, const Point<3> & p, Vec<3> & n) const
{
n = GetNormal(surfind, p);
}
virtual void PointBetween (const Point<3> & p1,
const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp,
PointGeomInfo & newgi) const
{
newp = p1 + secpoint * (p2-p1);
}
virtual void PointBetweenEdge(const Point<3> & p1,
const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp,
EdgePointGeomInfo & newgi) const
{
newp = p1+secpoint*(p2-p1);
}
virtual Vec<3> GetTangent(const Point<3> & p, int surfi1,
int surfi2,
const EdgePointGeomInfo & egi) const
{ throw Exception("Call GetTangent of " + Demangle(typeid(*this).name())); }
virtual void Save (string filename) const;
virtual void SaveToMeshFile (ostream & /* ost */) const { ; }
};

114
libsrc/meshing/bisect.cpp
View File

@ -3430,11 +3430,11 @@ namespace netgen
PointGeomInfo npgi;
if (mesh[newp].Type() != EDGEPOINT)
PointBetween (mesh.Point (oldpi1), mesh.Point (oldpi2),
0.5, si,
oldtri.pgeominfo[(oldtri.markededge+1)%3],
oldtri.pgeominfo[(oldtri.markededge+2)%3],
mesh.Point (newp), npgi);
geo.PointBetween (mesh.Point (oldpi1), mesh.Point (oldpi2),
0.5, si,
oldtri.pgeominfo[(oldtri.markededge+1)%3],
oldtri.pgeominfo[(oldtri.markededge+2)%3],
mesh.Point (newp), npgi);
BTBisectTri (oldtri, newp, npgi, newtri1, newtri2);
@ -3508,28 +3508,16 @@ namespace netgen
PointGeomInfo npgi1, npgi2;
int si = mesh.GetFaceDescriptor (oldquad.surfid).SurfNr();
// geom->GetSurface(si)->Project (mesh.Point(newp1));
// geom->GetSurface(si)->Project (mesh.Point(newp2));
// (*testout)
// cerr << "project point 1 " << newp1 << " old: " << mesh.Point(newp1);
PointBetween (mesh.Point (edge1.I1()), mesh.Point (edge1.I2()),
0.5, si,
pgi11,
pgi12,
mesh.Point (newp1), npgi1);
// (*testout)
// cerr << " new: " << mesh.Point(newp1) << endl;
// cerr << "project point 2 " << newp2 << " old: " << mesh.Point(newp2);
PointBetween (mesh.Point (edge2.I1()), mesh.Point (edge2.I2()),
0.5, si,
pgi21,
pgi22,
mesh.Point (newp2), npgi2);
// cerr << " new: " << mesh.Point(newp2) << endl;
geo.PointBetween(mesh.Point (edge1.I1()), mesh.Point (edge1.I2()),
0.5, si,
pgi11,
pgi12,
mesh.Point (newp1), npgi1);
geo.PointBetween (mesh.Point (edge2.I1()), mesh.Point (edge2.I2()),
0.5, si,
pgi21,
pgi22,
mesh.Point (newp2), npgi2);
BTBisectQuad (oldquad, newp1, npgi1, newp2, npgi2,
newquad1, newquad2);
@ -3565,16 +3553,10 @@ namespace netgen
EdgePointGeomInfo newepgi;
//
// cerr << "move edgepoint " << newpi << " from " << mesh.Point(newpi);
PointBetween (mesh.Point (seg[0]), mesh.Point (seg[1]),
0.5, seg.surfnr1, seg.surfnr2,
seg.epgeominfo[0], seg.epgeominfo[1],
mesh.Point (newpi), newepgi);
// cerr << " to " << mesh.Point (newpi) << endl;
geo.PointBetweenEdge(mesh.Point (seg[0]), mesh.Point (seg[1]),
0.5, seg.surfnr1, seg.surfnr2,
seg.epgeominfo[0], seg.epgeominfo[1],
mesh.Point (newpi), newepgi);
nseg1.epgeominfo[1] = newepgi;
nseg2.epgeominfo[0] = newepgi;
@ -4141,62 +4123,4 @@ namespace netgen
refine_hp = 0;
refine_p = 0;
}
Refinement :: Refinement ()
{
optimizer2d = NULL;
}
Refinement :: ~Refinement ()
{
;
}
void Refinement :: PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
newp = p1+secpoint*(p2-p1);
}
void Refinement :: PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
//cout << "base class edge point between" << endl;
newp = p1+secpoint*(p2-p1);
}
Vec<3> Refinement :: GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const
{
cerr << "Refinement::GetTangent not overloaded" << endl;
return Vec<3> (0,0,0);
}
Vec<3> Refinement :: GetNormal (const Point<3> & p, int surfi1,
const PointGeomInfo & gi) const
{
cerr << "Refinement::GetNormal not overloaded" << endl;
return Vec<3> (0,0,0);
}
void Refinement :: ProjectToSurface (Point<3> & p, int surfi) const
{
if (printmessage_importance>0)
cerr << "Refinement :: ProjectToSurface ERROR: no geometry set" << endl;
};
void Refinement :: ProjectToEdge (Point<3> & p, int surfi1, int surfi2, const EdgePointGeomInfo & egi) const
{
cerr << "Refinement::ProjectToEdge not overloaded" << endl;
}
}

45
libsrc/meshing/bisect.hpp
View File

@ -38,11 +38,11 @@ DLL_HEADER extern void ZRefinement (Mesh &, const class NetgenGeometry *,
class DLL_HEADER Refinement
{
MeshOptimize2d * optimizer2d;
const NetgenGeometry& geo;
public:
Refinement ();
virtual ~Refinement ();
Refinement (const NetgenGeometry& ageo) : geo(ageo) {}
virtual ~Refinement () {}
void Refine (Mesh & mesh) const;
void Refine (Mesh & mesh);
@ -51,48 +51,9 @@ public:
void MakeSecondOrder (Mesh & mesh) const;
void MakeSecondOrder (Mesh & mesh);
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const;
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const;
virtual Vec<3> GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & egi) const;
virtual Vec<3> GetNormal (const Point<3> & p, int surfi1,
const PointGeomInfo & gi) const;
virtual void ProjectToSurface (Point<3> & p, int surfi) const;
virtual void ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & /* gi */) const
{
ProjectToSurface (p, surfi);
}
virtual void ProjectToEdge (Point<3> & p, int surfi1, int surfi2, const EdgePointGeomInfo & egi) const;
void ValidateSecondOrder (Mesh & mesh);
void ValidateRefinedMesh (Mesh & mesh,
NgArray<INDEX_2> & parents);
MeshOptimize2d * Get2dOptimizer(void) const
{
return optimizer2d;
}
void Set2dOptimizer(MeshOptimize2d * opti)
{
optimizer2d = opti;
}
virtual void LocalizeEdgePoints(Mesh & /* mesh */) const {;}
};

48
libsrc/meshing/curvedelems.cpp
View File

@ -539,7 +539,7 @@ namespace netgen
CurvedElements :: CurvedElements (const Mesh & amesh)
: mesh (amesh)
: mesh(amesh), geo(*mesh.GetGeometry())
{
order = 1;
rational = 0;
@ -838,8 +838,8 @@ namespace netgen
{
Point<3> pm = Center (p1, p2);
Vec<3> n1 = ref -> GetNormal (p1, surfnr[e], gi0[e]);
Vec<3> n2 = ref -> GetNormal (p2, surfnr[e], gi1[e]);
Vec<3> n1 = geo.GetNormal (surfnr[e], p1, gi0[e]);
Vec<3> n2 = geo.GetNormal (surfnr[e], p2, gi1[e]);
// p3 = pm + alpha1 n1 + alpha2 n2
@ -876,7 +876,7 @@ namespace netgen
Vec<3> v05 = 0.25 * Vec<3> (p1) + 0.5*w* Vec<3>(p3) + 0.25 * Vec<3> (p2);
v05 /= 1 + (w-1) * 0.5;
Point<3> p05 (v05), pp05(v05);
ref -> ProjectToSurface (pp05, surfnr[e], gi0[e]);
geo.ProjectPointGI(surfnr[e], pp05, gi0[e]);
double d = Dist (pp05, p05);
if (d < dold)
@ -911,16 +911,16 @@ namespace netgen
if (swap)
{
p = p1 + xi[j] * (p2-p1);
ref -> PointBetween (p1, p2, xi[j],
surfnr[e], gi0[e], gi1[e],
pp, ppgi);
geo.PointBetween (p1, p2, xi[j],
surfnr[e], gi0[e], gi1[e],
pp, ppgi);
}
else
{
p = p2 + xi[j] * (p1-p2);
ref -> PointBetween (p2, p1, xi[j],
surfnr[e], gi1[e], gi0[e],
pp, ppgi);
geo.PointBetween (p2, p1, xi[j],
surfnr[e], gi1[e], gi0[e],
pp, ppgi);
}
Vec<3> dist = pp - p;
@ -1053,10 +1053,10 @@ namespace netgen
if (rational)
{
Vec<3> tau1 = ref -> GetTangent (p1, edge_surfnr2[edgenr], edge_surfnr1[edgenr],
edge_gi0[edgenr]);
Vec<3> tau2 = ref -> GetTangent (p2, edge_surfnr2[edgenr], edge_surfnr1[edgenr],
edge_gi1[edgenr]);
Vec<3> tau1 = geo.GetTangent(p1, edge_surfnr2[edgenr], edge_surfnr1[edgenr],
edge_gi0[edgenr]);
Vec<3> tau2 = geo.GetTangent(p2, edge_surfnr2[edgenr], edge_surfnr1[edgenr],
edge_gi1[edgenr]);
// p1 + alpha1 tau1 = p2 + alpha2 tau2;
Mat<3,2> mat;
@ -1082,8 +1082,8 @@ namespace netgen
Vec<3> v05 = 0.25 * Vec<3> (p1) + 0.5*w* Vec<3>(p3) + 0.25 * Vec<3> (p2);
v05 /= 1 + (w-1) * 0.5;
Point<3> p05 (v05), pp05(v05);
ref -> ProjectToEdge (pp05, edge_surfnr1[edgenr], edge_surfnr2[edgenr],
edge_gi0[edgenr]);
geo.ProjectPointEdge(edge_surfnr1[edgenr], edge_surfnr2[edgenr], pp05,
edge_gi0[edgenr]);
double d = Dist (pp05, p05);
if (d < dold)
@ -1127,16 +1127,16 @@ namespace netgen
if (swap)
{
p = p1 + xi[j] * (p2-p1);
ref -> PointBetween (p1, p2, xi[j],
edge_surfnr2[edgenr], edge_surfnr1[edgenr],
edge_gi0[edgenr], edge_gi1[edgenr],
pp, ppgi);
geo.PointBetweenEdge(p1, p2, xi[j],
edge_surfnr2[edgenr], edge_surfnr1[edgenr],
edge_gi0[edgenr], edge_gi1[edgenr],
pp, ppgi);
}
else
{
p = p2 + xi[j] * (p1-p2);
ref -> PointBetween (p2, p1, xi[j],
edge_surfnr2[edgenr], edge_surfnr1[edgenr],
geo.PointBetweenEdge(p2, p1, xi[j],
edge_surfnr2[edgenr], edge_surfnr1[edgenr],
edge_gi1[edgenr], edge_gi0[edgenr],
pp, ppgi);
}
@ -1302,10 +1302,10 @@ namespace netgen
SurfaceElementIndex sei = top.GetFace2SurfaceElement (f+1)-1;
if (sei != SurfaceElementIndex(-1)) {
PointGeomInfo gi = mesh[sei].GeomInfoPi(1);
ref -> ProjectToSurface (pp, surfnr[facenr], gi);
geo.ProjectPointGI(surfnr[facenr], pp, gi);
}
else
{ ref -> ProjectToSurface (pp, surfnr[facenr]); }
{ geo.ProjectPoint(surfnr[facenr], pp); }
Vec<3> dist = pp-xa[jj];
CalcTrigShape (order1, lami[fnums[1]]-lami[fnums[0]],

1
libsrc/meshing/curvedelems.hpp
View File

@ -17,6 +17,7 @@ class Refinement;
class CurvedElements
{
const Mesh & mesh;
const NetgenGeometry& geo;
NgArray<int> edgeorder;
NgArray<int> faceorder;

162
libsrc/meshing/improve2.cpp
View File

@ -19,7 +19,7 @@ namespace netgen
};
bool MeshOptimize2d :: EdgeSwapping (Mesh & mesh, const int usemetric,
bool MeshOptimize2d :: EdgeSwapping (const int usemetric,
Array<Neighbour> &neighbors,
Array<bool> &swapped,
const SurfaceElementIndex t1, const int o1,
@ -85,12 +85,11 @@ namespace netgen
nv1.Normalize();
nv2.Normalize();
Vec<3> nvp3, nvp4;
GetNormalVector (surfnr, mesh.Point(pi3), gi3, nvp3);
auto nvp3 = geo.GetNormal (surfnr, mesh.Point(pi3), gi3);
nvp3.Normalize();
GetNormalVector (surfnr, mesh.Point(pi4), gi4, nvp4);
auto nvp4 = geo.GetNormal (surfnr, mesh.Point(pi4), gi4);
nvp4.Normalize();
@ -168,7 +167,7 @@ namespace netgen
}
void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
void MeshOptimize2d :: EdgeSwapping (int usemetric)
{
static Timer timer("EdgeSwapping (2D)"); RegionTimer reg(timer);
static Timer timer_nb("EdgeSwapping-Find neighbors");
@ -203,7 +202,7 @@ namespace netgen
}
if(mixed)
return GenericImprove(mesh);
return GenericImprove();
Array<Neighbour> neighbors(mesh.GetNSE());
auto elements_on_node = mesh.CreatePoint2SurfaceElementTable(faceindex);
@ -262,7 +261,7 @@ namespace netgen
}
});
ParallelFor( Range(seia), [&pdef, &neighbors, &mesh, &seia, &elements_on_node] (auto i) NETGEN_LAMBDA_INLINE
ParallelFor( Range(seia), [this, &pdef, &neighbors, &seia, &elements_on_node] (auto i) NETGEN_LAMBDA_INLINE
{
auto sei = seia[i];
for (PointIndex pi : mesh[sei].template PNums<3>())
@ -336,7 +335,7 @@ namespace netgen
throw NgException ("Meshing stopped");
for (int o1 = 0; o1 < 3; o1++)
if(EdgeSwapping(mesh, usemetric, neighbors, swapped, t1, o1, t, pdef, true))
if(EdgeSwapping(usemetric, neighbors, swapped, t1, o1, t, pdef, true))
improvement_candidates[cnt++]= std::make_pair(t1,o1);
});
@ -344,7 +343,7 @@ namespace netgen
QuickSort(elements_with_improvement);
for (auto [t1,o1] : elements_with_improvement)
done |= EdgeSwapping(mesh, usemetric, neighbors, swapped, t1, o1, t, pdef, false);
done |= EdgeSwapping(usemetric, neighbors, swapped, t1, o1, t, pdef, false);
t--;
}
@ -566,9 +565,9 @@ namespace netgen
return d_badness;
}
void MeshOptimize2d :: CombineImprove (Mesh & mesh)
void MeshOptimize2d :: CombineImprove ()
{
SplitImprove(mesh);
SplitImprove();
PrintMessage (3, "Combine improve");
if (multithread.terminate)
@ -649,7 +648,7 @@ namespace netgen
{
const int faceindex = hel.GetIndex();
const int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
GetNormalVector (surfnr, mesh[pi], hel.GeomInfoPi(k+1), normals[pi]);
normals[pi] = geo.GetNormal (surfnr, mesh[pi], hel.GeomInfoPi(k+1));
break;
}
}
@ -682,7 +681,7 @@ namespace netgen
mesh.SetNextTimeStamp();
}
void MeshOptimize2d :: SplitImprove (Mesh & mesh)
void MeshOptimize2d :: SplitImprove()
{
if (!faceindex)
{
@ -691,7 +690,7 @@ namespace netgen
mesh.CalcSurfacesOfNode(); // TODO: needed?
for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
{
SplitImprove (mesh);
SplitImprove();
if (multithread.terminate)
throw NgException ("Meshing stopped");
@ -750,70 +749,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;
geo.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();

49
libsrc/meshing/improve2.hpp
View File

@ -75,29 +75,31 @@ public:
///
class MeshOptimize2d
{
int faceindex;
int improveedges;
double metricweight;
int writestatus;
int faceindex = 0;
int improveedges = 0;
double metricweight = 0.;
int writestatus = 1;
Mesh& mesh;
const NetgenGeometry& geo;
public:
///
MeshOptimize2d ();
MeshOptimize2d(Mesh& amesh) : mesh(amesh), geo(*mesh.GetGeometry())
{}
virtual ~MeshOptimize2d() { ; }
///
void ImproveMesh (Mesh & mesh2d, const MeshingParameters & mp);
void ImproveMeshJacobian (Mesh & mesh2d, const MeshingParameters & mp);
void ImproveVolumeMesh (Mesh & mesh);
void ImproveMesh (const MeshingParameters & mp);
void ImproveMeshJacobian (const MeshingParameters & mp);
void ImproveVolumeMesh ();
void ProjectBoundaryPoints(NgArray<int> & surfaceindex,
const NgArray<Point<3>* > & from, NgArray<Point<3>* > & dest);
bool EdgeSwapping (Mesh & mesh, const int usemetric, Array<Neighbour> &neighbors, Array<bool> &swapped,
bool EdgeSwapping (const int usemetric, Array<Neighbour> &neighbors, Array<bool> &swapped,
const SurfaceElementIndex t1, const int edge, const int t, Array<int,PointIndex> &pdef, const bool check_only=false);
void EdgeSwapping (Mesh & mesh, int usemetric);
void CombineImprove (Mesh & mesh);
void SplitImprove (Mesh & mesh);
void EdgeSwapping (int usemetric);
void CombineImprove ();
void SplitImprove ();
void GenericImprove (Mesh & mesh);
void GenericImprove ();
void SetFaceIndex (int fi) { faceindex = fi; }
@ -106,28 +108,9 @@ public:
void SetWriteStatus (int ws) { writestatus = ws; }
///
virtual void ProjectPoint (INDEX /* surfind */, Point<3> & /* p */) const { };
/// project point, use gi as initial value, and compute new gi
virtual int ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const
{ ProjectPoint (surfind, p); return CalcPointGeomInfo (surfind, gi, p); }
///
virtual void ProjectPoint2 (INDEX /* surfind */, INDEX /* surfind2 */, Point<3> & /* p */) const { };
/// liefert zu einem 3d-Punkt die geominfo (Dreieck) und liefert 1, wenn erfolgreich,
/// 0, wenn nicht (Punkt ausserhalb von chart)
virtual int CalcPointGeomInfo(PointGeomInfo& gi, const Point<3> & /*p3*/) const
{ gi.trignum = 1; return 1;};
virtual int CalcPointGeomInfo(int /* surfind */, PointGeomInfo& gi, const Point<3> & p3) const
{ return CalcPointGeomInfo (gi, p3); }
///
virtual void GetNormalVector(INDEX surfind, const Point<3> & p, PointGeomInfo & gi, Vec<3> & n) const;
virtual void GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const;
void CheckMeshApproximation (Mesh & mesh);

9
libsrc/meshing/improve2gen.cpp
View File

@ -19,7 +19,7 @@ namespace netgen
};
void MeshOptimize2d :: GenericImprove (Mesh & mesh)
void MeshOptimize2d :: GenericImprove ()
{
static Timer timer("MeshOptimize2d::GenericImprove"); RegionTimer reg(timer);
if (!faceindex)
@ -28,7 +28,7 @@ namespace netgen
PrintMessage (3, "Generic Improve");
for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
GenericImprove (mesh);
GenericImprove ();
faceindex = 0;
}
@ -396,9 +396,8 @@ namespace netgen
// calc metric badness
double bad1 = 0, bad2 = 0;
Vec<3> n;
GetNormalVector (surfnr, mesh.Point(pmap.Get(1)), pgi.Elem(1), n);
// SelectSurfaceOfPoint (mesh.Point(pmap.Get(1)), pgi.Get(1));
auto n = geo.GetNormal(surfnr, mesh.Point(pmap.Get(1)), pgi.Elem(1));
for (int j = 0; j < rule.oldels.Size(); j++)
bad1 += mesh[elmap[j]].CalcJacobianBadness (mesh.Points(), n);

2
libsrc/meshing/meshclass.hpp
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;

16
libsrc/meshing/meshfunc2d.cpp
View File

@ -31,30 +31,30 @@ namespace netgen
{
case 's':
{ // topological swap
MeshOptimize2d meshopt;
MeshOptimize2d meshopt(mesh);
meshopt.SetMetricWeight (mp.elsizeweight);
meshopt.EdgeSwapping (mesh, 0);
meshopt.EdgeSwapping (0);
break;
}
case 'S':
{ // metric swap
MeshOptimize2d meshopt;
MeshOptimize2d meshopt(mesh);
meshopt.SetMetricWeight (mp.elsizeweight);
meshopt.EdgeSwapping (mesh, 1);
meshopt.EdgeSwapping (1);
break;
}
case 'm':
{
MeshOptimize2d meshopt;
MeshOptimize2d meshopt(mesh);
meshopt.SetMetricWeight (mp.elsizeweight);
meshopt.ImproveMesh(mesh, mp);
meshopt.ImproveMesh(mp);
break;
}
case 'c':
{
MeshOptimize2d meshopt;
MeshOptimize2d meshopt(mesh);
meshopt.SetMetricWeight (mp.elsizeweight);
meshopt.CombineImprove(mesh);
meshopt.CombineImprove();
break;
}
default:

2
libsrc/meshing/meshtype.hpp
View File

@ -1234,7 +1234,7 @@ namespace netgen
// P .. plot, pause
// c .. combine
**/
string optimize2d = "smsmsmSmSmSm";
string optimize2d = "smcmSmcmSmcm";
/// number of 2d optimization steps
int optsteps2d = 3;
/// power of error (to approximate max err optimization)

2
libsrc/meshing/python_mesh.hpp
View File

@ -66,7 +66,7 @@ optimize3d: str = "cmdmustm"
optsteps3d: int = 3
Number of 3d optimization steps.
optimize2d: str = "smsmsmSmSmSm"
optimize2d: str = "smcmSmcmSmcm"
2d optimization strategy:
s .. swap, opt 6 lines/node
S .. swap, optimal elements

34
libsrc/meshing/refine.cpp
View File

@ -147,11 +147,11 @@ namespace netgen
{
pointset[pinew] = true;
Point<3> pnew;
PointBetween (mesh.Point (el[0]),
mesh.Point (el[1]), 0.5,
el.surfnr1, el.surfnr2,
el.epgeominfo[0], el.epgeominfo[1],
pnew, ngi);
geo.PointBetweenEdge(mesh.Point (el[0]),
mesh.Point (el[1]), 0.5,
el.surfnr1, el.surfnr2,
el.epgeominfo[0], el.epgeominfo[1],
pnew, ngi);
// pinew = mesh.AddPoint (pnew);
mesh.Point(pinew) = pnew;
@ -216,12 +216,12 @@ namespace netgen
Point<3> pb;
PointGeomInfo pgi;
PointBetween (mesh.Point (pi1),
mesh.Point (pi2), 0.5,
mesh.GetFaceDescriptor(el.GetIndex ()).SurfNr(),
el.GeomInfoPi (betw[j][0]),
el.GeomInfoPi (betw[j][1]),
pb, pgi);
geo.PointBetween(mesh.Point (pi1),
mesh.Point (pi2), 0.5,
mesh.GetFaceDescriptor(el.GetIndex ()).SurfNr(),
el.GeomInfoPi (betw[j][0]),
el.GeomInfoPi (betw[j][1]),
pb, pgi);
pgis.Elem(4+j) = pgi;
@ -307,12 +307,12 @@ namespace netgen
else
{
Point<3> pb;
PointBetween (mesh.Point (pi1),
mesh.Point (pi2), 0.5,
mesh.GetFaceDescriptor(el.GetIndex ()).SurfNr(),
el.GeomInfoPi (betw[j][0]),
el.GeomInfoPi (betw[j][1]),
pb, pgis.Elem(5+j));
geo.PointBetween(mesh.Point (pi1),
mesh.Point (pi2), 0.5,
mesh.GetFaceDescriptor(el.GetIndex ()).SurfNr(),
el.GeomInfoPi (betw[j][0]),
el.GeomInfoPi (betw[j][1]),
pb, pgis.Elem(5+j));
pnums.Elem(5+j) = mesh.AddPoint (pb);

22
libsrc/meshing/secondorder.cpp
View File

@ -100,11 +100,11 @@ namespace netgen
{
Point<3> pb;
EdgePointGeomInfo ngi;
PointBetween (mesh.Point (el[0]),
mesh.Point (el[1]), 0.5,
el.surfnr1, el.surfnr2,
el.epgeominfo[0], el.epgeominfo[1],
pb, ngi);
geo.PointBetweenEdge(mesh.Point (el[0]),
mesh.Point (el[1]), 0.5,
el.surfnr1, el.surfnr2,
el.epgeominfo[0], el.epgeominfo[1],
pb, ngi);
el[2] = mesh.AddPoint (pb, mesh.Point(el[0]).GetLayer(),
EDGEPOINT);
@ -184,12 +184,12 @@ namespace netgen
{
Point<3> pb;
PointGeomInfo newgi;
PointBetween (mesh.Point (pi1),
mesh.Point (pi2), 0.5,
mesh.GetFaceDescriptor(el.GetIndex ()).SurfNr(),
el.GeomInfoPi (betw[j][0]+1),
el.GeomInfoPi (betw[j][1]+1),
pb, newgi);
geo.PointBetween(mesh.Point (pi1),
mesh.Point (pi2), 0.5,
mesh.GetFaceDescriptor(el.GetIndex ()).SurfNr(),
el.GeomInfoPi (betw[j][0]+1),
el.GeomInfoPi (betw[j][1]+1),
pb, newgi);
newel[onp+j] = mesh.AddPoint (pb, mesh.Point(pi1).GetLayer(),
SURFACEPOINT);

8
libsrc/meshing/smoothing2.5.cpp
View File

@ -15,14 +15,14 @@ namespace netgen
if(surfaceindex[i] >= 0)
{
*dest[i] = *from[i];
ProjectPoint(surfaceindex[i],*dest[i]);
geo.ProjectPoint(surfaceindex[i],*dest[i]);
}
}
}
void MeshOptimize2d :: ImproveVolumeMesh (Mesh & mesh)
void MeshOptimize2d :: ImproveVolumeMesh ()
{
if (!faceindex)
@ -31,7 +31,7 @@ namespace netgen
for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
{
ImproveVolumeMesh (mesh);
ImproveVolumeMesh ();
if (multithread.terminate)
throw NgException ("Meshing stopped");
}
@ -229,7 +229,7 @@ namespace netgen
//cout << "origp " << origp << " newp " << mesh[pi];
ngi = gi1;
moveisok = (ProjectPointGI (surfi, mesh[pi], ngi) != 0);
moveisok = (geo.ProjectPointGI(surfi, mesh[pi], ngi) != 0);
//cout << " projected " << mesh[pi] << endl;

72
libsrc/meshing/smoothing2.cpp
View File

@ -205,22 +205,20 @@ namespace netgen
class Opti2SurfaceMinFunction : public MinFunction
{
const Mesh & mesh;
Opti2dLocalData & ld;
const NetgenGeometry& geo;
public:
Opti2SurfaceMinFunction (const Mesh & amesh,
Opti2dLocalData & ald)
: mesh(amesh), ld(ald)
: ld(ald), geo(*amesh.GetGeometry())
{ } ;
virtual double Func (const Vector & x) const
{
Vec<3> n;
double badness = 0;
ld.meshthis -> GetNormalVector (ld.surfi, ld.sp1, ld.gi1, n);
auto n = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1);
Point<3> pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
for (int j = 0; j < ld.locelements.Size(); j++)
@ -355,13 +353,13 @@ namespace netgen
// static int timer = NgProfiler::CreateTimer ("opti2surface - funcgrad");
// NgProfiler::RegionTimer reg (timer);
Vec<3> n, vgrad;
Vec<3> vgrad;
Point<3> pp1;
vgrad = 0;
double badness = 0;
ld.meshthis -> GetNormalVector (ld.surfi, ld.sp1, ld.gi1, n);
auto n = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1);
pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
// meshthis -> ProjectPoint (surfi, pp1);
@ -410,13 +408,13 @@ namespace netgen
// static int timer = NgProfiler::CreateTimer ("opti2surface - funcderiv");
// NgProfiler::RegionTimer reg (timer);
Vec<3> n, vgrad;
Vec<3> vgrad;
Point<3> pp1;
vgrad = 0;
double badness = 0;
ld.meshthis -> GetNormalVector (ld.surfi, ld.sp1, ld.gi1, n);
auto n = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1);
pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
for (int j = 0; j < ld.locelements.Size(); j++)
@ -474,11 +472,12 @@ namespace netgen
{
const Mesh & mesh;
Opti2dLocalData & ld;
const NetgenGeometry& geo;
public:
Opti2EdgeMinFunction (const Mesh & amesh,
Opti2dLocalData & ald)
: mesh(amesh), ld(ald) { } ;
: mesh(amesh), ld(ald), geo(*amesh.GetGeometry()) { } ;
virtual double FuncGrad (const Vector & x, Vector & g) const;
virtual double Func (const Vector & x) const;
@ -493,7 +492,7 @@ namespace netgen
double Opti2EdgeMinFunction :: FuncGrad (const Vector & x, Vector & grad) const
{
int j, rot;
Vec<3> n1, n2, v1, v2, e1, e2, vgrad;
Vec<3> v1, v2, e1, e2, vgrad;
Point<3> pp1;
Vec<2> g1;
double badness, hbadness;
@ -502,7 +501,7 @@ namespace netgen
badness = 0;
pp1 = ld.sp1 + x(0) * ld.t1;
ld.meshthis -> ProjectPoint2 (ld.surfi, ld.surfi2, pp1);
geo.ProjectPointEdge(ld.surfi, ld.surfi2, pp1);
for (j = 0; j < ld.locelements.Size(); j++)
{
@ -526,8 +525,8 @@ namespace netgen
vgrad += g1(0) * e1 + g1(1) * e2;
}
ld.meshthis -> GetNormalVector (ld.surfi, pp1, n1);
ld.meshthis -> GetNormalVector (ld.surfi2, pp1, n2);
auto n1 = geo.GetNormal(ld.surfi, pp1);
auto n2 = geo.GetNormal(ld.surfi2, pp1);
v1 = Cross (n1, n2);
v1.Normalize();
@ -544,11 +543,12 @@ namespace netgen
{
const Mesh & mesh;
Opti2dLocalData & ld;
const NetgenGeometry& geo;
public:
Opti2SurfaceMinFunctionJacobian (const Mesh & amesh,
Opti2dLocalData & ald)
: mesh(amesh), ld(ald)
: mesh(amesh), ld(ald), geo(*amesh.GetGeometry())
{ } ;
virtual double FuncGrad (const Vector & x, Vector & g) const;
virtual double FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const;
@ -569,7 +569,7 @@ namespace netgen
// from 2d:
int lpi, gpi;
Vec<3> n, vgrad;
Vec<3> vgrad;
Point<3> pp1;
Vec<2> g1, vdir;
double badness, hbad, hderiv;
@ -577,7 +577,7 @@ namespace netgen
vgrad = 0;
badness = 0;
ld.meshthis -> GetNormalVector (ld.surfi, ld.sp1, ld.gi1, n);
auto n = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1);
pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
@ -641,7 +641,7 @@ namespace netgen
// from 2d:
int j, k, lpi, gpi;
Vec<3> n, vgrad;
Vec<3> vgrad;
Point<3> pp1;
Vec<2> g1, vdir;
double badness, hbad, hderiv;
@ -649,8 +649,6 @@ namespace netgen
vgrad = 0;
badness = 0;
ld.meshthis -> GetNormalVector (ld.surfi, ld.sp1, ld.gi1, n);
// pp1 = sp1;
// pp1.Add2 (x.Get(1), t1, x.Get(2), t2);
pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
@ -690,24 +688,7 @@ namespace netgen
return badness;
}
MeshOptimize2d dummy;
MeshOptimize2d :: MeshOptimize2d ()
{
SetFaceIndex (0);
SetImproveEdges (0);
SetMetricWeight (0);
SetWriteStatus (1);
}
void MeshOptimize2d :: ImproveMesh (Mesh & mesh, const MeshingParameters & mp)
void MeshOptimize2d :: ImproveMesh (const MeshingParameters & mp)
{
static Timer timer("MeshSmoothing 2D"); RegionTimer reg (timer);
@ -919,6 +900,7 @@ namespace netgen
}
ld.gi1 = hel.GeomInfoPi(hpi);
// SelectSurfaceOfPoint (ld.sp1, ld.gi1);
ld.locelements.SetSize(0);
ld.locrots.SetSize (0);
@ -951,7 +933,7 @@ namespace netgen
}
GetNormalVector (ld.surfi, ld.sp1, ld.gi1, ld.normal);
ld.normal = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1);
ld.t1 = ld.normal.GetNormal ();
ld.t2 = Cross (ld.normal, ld.t1);
@ -1029,7 +1011,7 @@ namespace netgen
PointGeomInfo ngi;
ngi = ld.gi1;
moveisok = ProjectPointGI (ld.surfi, mesh[pi], ngi);
moveisok = geo.ProjectPointGI(ld.surfi, mesh[pi], ngi);
// point lies on same chart in stlsurface
if (moveisok)
@ -1052,14 +1034,4 @@ namespace netgen
CheckMeshApproximation (mesh);
mesh.SetNextTimeStamp();
}
void MeshOptimize2d :: GetNormalVector(INDEX /* surfind */, const Point<3> & p, Vec<3> & nv) const
{
nv = Vec<3> (0, 0, 1);
}
void MeshOptimize2d :: GetNormalVector(INDEX surfind, const Point<3> & p, PointGeomInfo & gi, Vec<3> & n) const
{
GetNormalVector (surfind, p, n);
}
}

15
libsrc/meshing/validate.cpp
View File

@ -276,8 +276,8 @@ namespace netgen
double oldlamedge,oldlamface;
MeshOptimize2d * optimizer2d = refinement.Get2dOptimizer();
if(!optimizer2d)
auto geo = mesh.GetGeometry();
if(!geo)
{
cerr << "No 2D Optimizer!" << endl;
return;
@ -382,8 +382,15 @@ namespace netgen
for (int i = 1; i <= np; i++)
*can.Elem(i) = mesh.Point(i);
if(optimizer2d)
optimizer2d->ProjectBoundaryPoints(surfaceindex,can,should);
if(geo)
for(int i=0; i<surfaceindex.Size(); i++)
{
if(surfaceindex[i] >= 0)
{
*should[i] = *can[i];
geo->ProjectPoint(surfaceindex[i],*should[i]);
}
}
}

6
libsrc/occ/Partition_Inter3d.cxx
View File

@ -206,7 +206,7 @@ static void PutInBounds (const TopoDS_Face& F,
Handle (Geom_Surface) S = BRep_Tool::Surface(F,L);
if (S->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
S = (*(Handle_Geom_RectangularTrimmedSurface*)&S)->BasisSurface();
S = Handle(Geom_RectangularTrimmedSurface)::DownCast(S)->BasisSurface();
}
if (!S->IsUPeriodic() && !S->IsVPeriodic())
return;
@ -702,7 +702,7 @@ TopTools_MapOfShape& Partition_Inter3d::TouchedFaces()
//purpose :
//=======================================================================
Handle_BRepAlgo_AsDes Partition_Inter3d::AsDes() const
Handle(BRepAlgo_AsDes) Partition_Inter3d::AsDes() const
{
return myAsDes;
}
@ -829,7 +829,7 @@ TopoDS_Vertex Partition_Inter3d::ReplaceSameDomainV(const TopoDS_Vertex& V,
//purpose :
//=======================================================================
Handle_BRepAlgo_AsDes Partition_Inter3d::SectionEdgesAD() const
Handle(BRepAlgo_AsDes) Partition_Inter3d::SectionEdgesAD() const
{
return mySectionEdgesAD;
}

8
libsrc/occ/Partition_Inter3d.hxx
View File

@ -96,13 +96,13 @@ public:
void FacesPartition(const TopoDS_Face& F1,const TopoDS_Face& F2) ;
Standard_Boolean IsDone(const TopoDS_Face& F1,const TopoDS_Face& F2) const;
TopTools_MapOfShape& TouchedFaces() ;
Handle_BRepAlgo_AsDes AsDes() const;
Handle(BRepAlgo_AsDes) AsDes() const;
TopTools_MapOfShape& NewEdges() ;
Standard_Boolean HasSameDomainF(const TopoDS_Shape& F) const;
Standard_Boolean IsSameDomainF(const TopoDS_Shape& F1,const TopoDS_Shape& F2) const;
const TopTools_ListOfShape& SameDomain(const TopoDS_Face& F) const;
TopoDS_Vertex ReplaceSameDomainV(const TopoDS_Vertex& V,const TopoDS_Edge& E) const;
Handle_BRepAlgo_AsDes SectionEdgesAD() const;
Handle(BRepAlgo_AsDes) SectionEdgesAD() const;
Standard_Boolean IsSectionEdge(const TopoDS_Edge& E) const;
Standard_Boolean HasSectionEdge(const TopoDS_Face& F) const;
Standard_Boolean IsSplitOn(const TopoDS_Edge& NewE,const TopoDS_Edge& OldE,const TopoDS_Face& F) const;
@ -134,11 +134,11 @@ private:
// Fields PRIVATE
//
Handle_BRepAlgo_AsDes myAsDes;
Handle(BRepAlgo_AsDes) myAsDes;
TopTools_DataMapOfShapeListOfShape myDone;
TopTools_MapOfShape myTouched;
TopTools_MapOfShape myNewEdges;
Handle_BRepAlgo_AsDes mySectionEdgesAD;
Handle(BRepAlgo_AsDes) mySectionEdgesAD;
TopTools_DataMapOfShapeListOfShape mySameDomainFM;
TopTools_DataMapOfShapeShape mySameDomainVM;

2
libsrc/occ/Partition_Spliter.hxx
View File

@ -143,7 +143,7 @@ private:
TopTools_DataMapOfShapeShape myFaceShapeMap;
TopTools_DataMapOfShapeShape myInternalFaces;
TopTools_DataMapOfShapeShape myIntNotClFaces;
Handle_BRepAlgo_AsDes myAsDes;
Handle(BRepAlgo_AsDes) myAsDes;
BRepAlgo_Image myImagesFaces;
BRepAlgo_Image myImagesEdges;
BRepAlgo_Image myImageShape;

220
libsrc/occ/occgenmesh.cpp
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,6 @@ namespace netgen
// Philippose - 15/01/2009
double maxh = geom.face_maxh[k-1];
//double maxh = mparam.maxh;
mparam.checkoverlap = 0;
// int noldpoints = mesh->GetNP();
int noldsurfel = mesh.GetNSE();
@ -809,9 +808,13 @@ namespace netgen
MESHING2_RESULT res;
// TODO: check overlap not correctly working here
MeshingParameters mparam_without_overlap = mparam;
mparam_without_overlap.checkoverlap = false;
try {
static Timer t("GenerateMesh"); RegionTimer reg(t);
res = meshing.GenerateMesh (mesh, mparam, maxh, k);
res = meshing.GenerateMesh (mesh, mparam_without_overlap, maxh, k);
}
catch (SingularMatrixException)
@ -916,7 +919,7 @@ namespace netgen
}
void OCCOptimizeSurface(OCCGeometry & geom, Mesh & mesh,
MeshingParameters & mparam)
const MeshingParameters & mparam)
{
const char * savetask = multithread.task;
multithread.task = "Optimizing surface";
@ -941,41 +944,41 @@ namespace netgen
if (multithread.terminate) return;
{
MeshOptimize2dOCCSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.EdgeSwapping (mesh, (i > mparam.optsteps2d/2));
meshopt.EdgeSwapping (i > mparam.optsteps2d/2);
}
if (multithread.terminate) return;
{
MeshOptimize2dOCCSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.ImproveMesh (mesh, mparam);
meshopt.ImproveMesh (mparam);
}
{
MeshOptimize2dOCCSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.CombineImprove (mesh);
meshopt.CombineImprove ();
}
if (multithread.terminate) return;
{
MeshOptimize2dOCCSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.ImproveMesh (mesh, mparam);
meshopt.ImproveMesh (mparam);
}
}
@ -991,7 +994,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 +1282,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

229
libsrc/occ/occgeom.cpp
View File

@ -74,6 +74,30 @@ 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 :: FinalizeMesh(Mesh& mesh) const
{
for (int i = 0; i < mesh.GetNDomains(); i++)
if (snames.Size())
mesh.SetMaterial (i+1, snames[i]);
}
void OCCGeometry :: PrintNrShapes ()
{
@ -1010,10 +1034,7 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
SetCenter();
}
void OCCGeometry :: Project (int surfi, Point<3> & p) const
void OCCGeometry :: ProjectPoint(int surfi, Point<3> & p) const
{
static int cnt = 0;
if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl;
@ -1032,8 +1053,47 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
}
bool OCCGeometry :: ProjectPointGI(int surfind, Point<3>& p, PointGeomInfo& gi) const
{
double u = gi.u;
double v = gi.v;
Point<3> hp = p;
if (FastProject (surfind, hp, u, v))
{
p = hp;
return 1;
}
ProjectPoint (surfind, p);
return CalcPointGeomInfo (surfind, gi, p);
}
void OCCGeometry :: ProjectPointEdge(int surfind, INDEX surfind2,
Point<3> & p) const
{
TopExp_Explorer exp0, exp1;
bool done = false;
Handle(Geom_Curve) c;
for (exp0.Init(fmap(surfind), TopAbs_EDGE); !done && exp0.More(); exp0.Next())
for (exp1.Init(fmap(surfind2), TopAbs_EDGE); !done && exp1.More(); exp1.Next())
{
if (TopoDS::Edge(exp0.Current()).IsSame(TopoDS::Edge(exp1.Current())))
{
done = true;
double s0, s1;
c = BRep_Tool::Curve(TopoDS::Edge(exp0.Current()), s0, s1);
}
}
gp_Pnt pnt(p(0), p(1), p(2));
GeomAPI_ProjectPointOnCurve proj(pnt, c);
pnt = proj.NearestPoint();
p(0) = pnt.X();
p(1) = pnt.Y();
p(2) = pnt.Z();
}
bool OCCGeometry :: FastProject (int surfi, Point<3> & ap, double& u, double& v) const
{
@ -1091,7 +1151,148 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
return true;
}
Vec<3> OCCGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & geominfo) const
{
gp_Pnt pnt;
gp_Vec du, dv;
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
occface->D1(geominfo.u,geominfo.v,pnt,du,dv);
auto n = Cross (Vec<3>(du.X(), du.Y(), du.Z()),
Vec<3>(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1;
return n;
}
Vec<3> OCCGeometry :: GetNormal(int surfind, const Point<3> & p) const
{
Standard_Real u,v;
gp_Pnt pnt(p(0), p(1), p(2));
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface);
if (proj.NbPoints() < 1)
{
cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
<< endl;
cout << p << endl;
return;
}
proj.LowerDistanceParameters (u, v);
*/
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfind)) ) );
suval.Coord( u, v);
pnt = occface->Value( u, v );
gp_Vec du, dv;
occface->D1(u,v,pnt,du,dv);
/*
if (!occface->IsCNu (1) || !occface->IsCNv (1))
(*testout) << "SurfOpt: Differentiation FAIL" << endl;
*/
auto n = Cross (Vec3d(du.X(), du.Y(), du.Z()),
Vec3d(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1;
return n;
}
bool OCCGeometry :: CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p) const
{
Standard_Real u,v;
gp_Pnt pnt(p(0), p(1), p(2));
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface);
if (proj.NbPoints() < 1)
{
cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
<< endl;
cout << p << endl;
return 0;
}
proj.LowerDistanceParameters (u, v);
*/
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfind)) ) );
suval.Coord( u, v);
//pnt = occface->Value( u, v );
gi.u = u;
gi.v = v;
return true;
}
void OCCGeometry :: PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
Point<3> hnewp;
hnewp = p1+secpoint*(p2-p1);
if (surfi > 0)
{
double u = gi1.u+secpoint*(gi2.u-gi1.u);
double v = gi1.v+secpoint*(gi2.v-gi1.v);
auto savept = hnewp;
if (!FastProject(surfi, hnewp, u, v) || Dist(hnewp, savept) > Dist(p1,p2))
{
// cout << "Fast projection to surface fails! Using OCC projection" << endl;
hnewp = savept;
ProjectPoint(surfi, hnewp);
}
newgi.trignum = 1;
newgi.u = u;
newgi.v = v;
}
newp = hnewp;
}
void OCCGeometry :: PointBetweenEdge(const Point<3> & p1,
const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
double s0, s1;
Point<3> hnewp = p1+secpoint*(p2-p1);
gp_Pnt pnt(hnewp(0), hnewp(1), hnewp(2));
GeomAPI_ProjectPointOnCurve proj(pnt, BRep_Tool::Curve(TopoDS::Edge(emap(ap1.edgenr)), s0, s1));
pnt = proj.NearestPoint();
hnewp = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
newp = hnewp;
newgi = ap1;
};
// void OCCGeometry :: WriteOCC_STL(char * filename)
@ -1681,17 +1882,7 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
return false;
}
const Refinement & OCCGeometry :: GetRefinement () const
{
return * new OCCRefinementSurfaces (*this);
}
void OCCParameters :: Print(ostream & ost) const
void OCCParameters :: Print(ostream & ost) const
{
ost << "OCC Parameters:" << endl
<< "close edges: " << resthcloseedgeenable
@ -1703,10 +1894,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);
// }
}

107
libsrc/occ/occgeom.hpp
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,12 +260,44 @@ namespace netgen
emap.Clear();
vmap.Clear();
}
Mesh::GEOM_TYPE GetGeomType() const override
{ return Mesh::GEOM_OCC; }
void SetOCCParameters(const OCCParameters& par)
{ 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 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;
void ProjectPoint(int surfind, Point<3> & p) const override;
void ProjectPointEdge (int surfind, int surfind2, Point<3> & p) const override;
bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override;
Vec<3> GetNormal(int surfind, const Point<3> & p) const override;
Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const override;
bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override;
void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
void PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
DLL_HEADER void BuildFMap();
Box<3> GetBoundingBox() const
@ -264,9 +317,6 @@ namespace netgen
Point<3> Center() const
{ return center; }
void Project (int surfi, Point<3> & p) const;
bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
OCCSurface GetSurface (int surfi)
{
cout << "OCCGeometry::GetSurface using PLANESPACE" << endl;
@ -391,37 +441,9 @@ namespace netgen
// void WriteOCC_STL(char * filename);
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 virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
private:
bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
};
@ -434,15 +456,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);
}

252
libsrc/occ/occmeshsurf.cpp
View File

@ -532,188 +532,6 @@ namespace netgen
return gh;
}
MeshOptimize2dOCCSurfaces :: MeshOptimize2dOCCSurfaces (const OCCGeometry & ageometry)
: MeshOptimize2d(), geometry(ageometry)
{
;
}
void MeshOptimize2dOCCSurfaces :: ProjectPoint (INDEX surfind, Point<3> & p) const
{
geometry.Project (surfind, p);
}
int MeshOptimize2dOCCSurfaces :: ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const
{
double u = gi.u;
double v = gi.v;
Point<3> hp = p;
if (geometry.FastProject (surfind, hp, u, v))
{
p = hp;
return 1;
}
ProjectPoint (surfind, p);
return CalcPointGeomInfo (surfind, gi, p);
}
void MeshOptimize2dOCCSurfaces :: ProjectPoint2 (INDEX surfind, INDEX surfind2,
Point<3> & p) const
{
TopExp_Explorer exp0, exp1;
bool done = false;
Handle(Geom_Curve) c;
for (exp0.Init(geometry.fmap(surfind), TopAbs_EDGE); !done && exp0.More(); exp0.Next())
for (exp1.Init(geometry.fmap(surfind2), TopAbs_EDGE); !done && exp1.More(); exp1.Next())
{
if (TopoDS::Edge(exp0.Current()).IsSame(TopoDS::Edge(exp1.Current())))
{
done = true;
double s0, s1;
c = BRep_Tool::Curve(TopoDS::Edge(exp0.Current()), s0, s1);
}
}
gp_Pnt pnt(p(0), p(1), p(2));
GeomAPI_ProjectPointOnCurve proj(pnt, c);
pnt = proj.NearestPoint();
p(0) = pnt.X();
p(1) = pnt.Y();
p(2) = pnt.Z();
}
void MeshOptimize2dOCCSurfaces ::
GetNormalVector(INDEX surfind, const Point<3> & p, PointGeomInfo & geominfo, Vec<3> & n) const
{
gp_Pnt pnt;
gp_Vec du, dv;
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind)));
occface->D1(geominfo.u,geominfo.v,pnt,du,dv);
n = Cross (Vec<3>(du.X(), du.Y(), du.Z()),
Vec<3>(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
if (geometry.fmap(surfind).Orientation() == TopAbs_REVERSED) n = -1*n;
// GetNormalVector (surfind, p, n);
}
void MeshOptimize2dOCCSurfaces ::
GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const
{
// static int cnt = 0;
// if (cnt++ % 1000 == 0) cout << "GetNV cnt = " << cnt << endl;
Standard_Real u,v;
gp_Pnt pnt(p(0), p(1), p(2));
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind)));
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface);
if (proj.NbPoints() < 1)
{
cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
<< endl;
cout << p << endl;
return;
}
proj.LowerDistanceParameters (u, v);
*/
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(geometry.fmap(surfind)) ) );
suval.Coord( u, v);
pnt = occface->Value( u, v );
gp_Vec du, dv;
occface->D1(u,v,pnt,du,dv);
/*
if (!occface->IsCNu (1) || !occface->IsCNv (1))
(*testout) << "SurfOpt: Differentiation FAIL" << endl;
*/
n = Cross (Vec3d(du.X(), du.Y(), du.Z()),
Vec3d(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
if (geometry.fmap(surfind).Orientation() == TopAbs_REVERSED) n = -1*n;
}
int MeshOptimize2dOCCSurfaces ::
CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p) const
{
Standard_Real u,v;
gp_Pnt pnt(p(0), p(1), p(2));
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind)));
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface);
if (proj.NbPoints() < 1)
{
cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
<< endl;
cout << p << endl;
return 0;
}
proj.LowerDistanceParameters (u, v);
*/
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(geometry.fmap(surfind)) ) );
suval.Coord( u, v);
//pnt = occface->Value( u, v );
gi.u = u;
gi.v = v;
return 1;
}
OCCRefinementSurfaces :: OCCRefinementSurfaces (const OCCGeometry & ageometry)
: Refinement(), geometry(ageometry)
{
;
}
OCCRefinementSurfaces :: ~OCCRefinementSurfaces ()
{
;
}
/*
inline double Det3 (double a00, double a01, double a02,
double a10, double a11, double a12,
@ -772,76 +590,6 @@ namespace netgen
return true;
}
*/
void OCCRefinementSurfaces ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
Point<3> hnewp;
hnewp = p1+secpoint*(p2-p1);
if (surfi > 0)
{
double u = gi1.u+secpoint*(gi2.u-gi1.u);
double v = gi1.v+secpoint*(gi2.v-gi1.v);
auto savept = hnewp;
if (!geometry.FastProject (surfi, hnewp, u, v) || Dist(hnewp, savept) > Dist(p1,p2))
{
// cout << "Fast projection to surface fails! Using OCC projection" << endl;
hnewp = savept;
geometry.Project (surfi, hnewp);
}
newgi.trignum = 1;
newgi.u = u;
newgi.v = v;
}
newp = hnewp;
}
void OCCRefinementSurfaces ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
double s0, s1;
Point<3> hnewp = p1+secpoint*(p2-p1);
gp_Pnt pnt(hnewp(0), hnewp(1), hnewp(2));
GeomAPI_ProjectPointOnCurve proj(pnt, BRep_Tool::Curve(TopoDS::Edge(geometry.emap(ap1.edgenr)), s0, s1));
pnt = proj.NearestPoint();
hnewp = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
newp = hnewp;
newgi = ap1;
};
void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi) const
{
if (surfi > 0)
geometry.Project (surfi, p);
};
void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi) const
{
if (surfi > 0)
if (!geometry.FastProject (surfi, p, gi.u, gi.v))
{
cout << "Fast projection to surface fails! Using OCC projection" << endl;
geometry.Project (surfi, p);
}
};
}

56
libsrc/occ/occmeshsurf.hpp
View File

@ -141,64 +141,8 @@ protected:
};
///
class MeshOptimize2dOCCSurfaces : public MeshOptimize2d
{
///
const OCCGeometry & geometry;
public:
///
MeshOptimize2dOCCSurfaces (const OCCGeometry & ageometry);
///
virtual void ProjectPoint (INDEX surfind, Point<3> & p) const;
///
virtual void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point<3> & p) const;
///
virtual int ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const;
///
virtual void GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const;
///
virtual void GetNormalVector(INDEX surfind, const Point<3> & p, PointGeomInfo & gi, Vec<3> & n) const;
virtual int CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const;
};
class OCCGeometry;
class DLL_HEADER OCCRefinementSurfaces : public Refinement
{
const OCCGeometry & geometry;
public:
OCCRefinementSurfaces (const OCCGeometry & ageometry);
virtual ~OCCRefinementSurfaces ();
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
virtual void ProjectToSurface (Point<3> & p, int surfi) const override;
virtual void ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi) const override;
};
#endif

2
libsrc/occ/occpkg.cpp
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);
}
};

5
libsrc/occ/python_occ.cpp
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>(),

236
libsrc/stlgeom/meshstlsurface.cpp
View File

@ -299,14 +299,14 @@ int STLSurfaceMeshing (STLGeometry & geom, class Mesh & mesh, const MeshingParam
geom.SetMarkedTrig(seg.geominfo[1].trignum,1);
}
MeshOptimizeSTLSurface optmesh(geom);
MeshOptimize2d optmesh(mesh);
optmesh.SetFaceIndex (0);
optmesh.SetImproveEdges (0);
optmesh.SetMetricWeight (0);
mesh.CalcSurfacesOfNode();
optmesh.EdgeSwapping (mesh, 0);
optmesh.ImproveMesh (mesh, mparam);
optmesh.EdgeSwapping (0);
optmesh.ImproveMesh (mparam);
}
mesh.Compress();
@ -826,7 +826,7 @@ void STLSurfaceOptimization (STLGeometry & geom,
{
PrintFnStart("optimize STL Surface");
MeshOptimizeSTLSurface optmesh(geom);
MeshOptimize2d optmesh(mesh);
optmesh.SetFaceIndex (0);
optmesh.SetImproveEdges (0);
@ -847,25 +847,41 @@ void STLSurfaceOptimization (STLGeometry & geom,
{
case 's':
{
optmesh.EdgeSwapping (mesh, 0);
optmesh.EdgeSwapping(0);
break;
}
case 'S':
{
optmesh.EdgeSwapping (mesh, 1);
optmesh.EdgeSwapping(1);
break;
}
case 'm':
{
optmesh.ImproveMesh(mesh, mparam);
optmesh.ImproveMesh(mparam);
break;
}
case 'c':
{
optmesh.CombineImprove (mesh);
optmesh.CombineImprove();
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();
// }
//(*testout) << "optimize, after, step = " << meshparam.optimize2d[j-1] << mesh.Point (3679) << endl;
}
@ -1051,208 +1067,4 @@ double MeshingSTLSurface :: Area () const
return geom.Area();
}
MeshOptimizeSTLSurface :: MeshOptimizeSTLSurface (STLGeometry & ageom)
: MeshOptimize2d(), geom(ageom)
{
;
}
void MeshOptimizeSTLSurface :: ProjectPoint (INDEX surfind, Point<3> & p) const
{
if (!geom.Project (p))
{
PrintMessage(7,"project failed");
if (!geom.ProjectOnWholeSurface(p))
{
PrintMessage(7, "project on whole surface failed");
}
}
// geometry.GetSurface(surfind)->Project (p);
}
int MeshOptimizeSTLSurface :: ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const
{
int meshchart = geom.GetChartNr(gi.trignum);
const STLChart& chart = geom.GetChart(meshchart);
int trignum = chart.ProjectNormal(p);
if(trignum==0)
{
PrintMessage(7,"project failed");
geom.SelectChartOfTriangle (gi.trignum); // needed because ProjectOnWholeSurface uses meshchartnv (the normal vector of selected chart)
trignum = geom.ProjectOnWholeSurface(p);
if(trignum==0)
PrintMessage(7, "project on whole surface failed");
}
return trignum;
}
void MeshOptimizeSTLSurface :: ProjectPoint2 (INDEX surfind, INDEX surfind2, Point<3> & p) const
{
/*
ProjectToEdge ( geometry.GetSurface(surfind),
geometry.GetSurface(surfind2), p);
*/
}
int MeshOptimizeSTLSurface :: CalcPointGeomInfo(PointGeomInfo& gi, const Point<3> & p3) const
{
Point<3> hp = p3;
gi.trignum = geom.Project (hp);
if (gi.trignum) return 1;
return 0;
}
void MeshOptimizeSTLSurface :: GetNormalVector(INDEX surfind, const Point<3> & p, PointGeomInfo & gi, Vec<3> & n) const
{
n = geom.GetTriangle(gi.trignum).Normal();
}
void MeshOptimizeSTLSurface :: GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const
{
throw Exception("MeshOptimizeSTLSurface :: GetNormalVector called without PointGeomInfo");
}
RefinementSTLGeometry :: RefinementSTLGeometry (const STLGeometry & ageom)
: Refinement(), geom(ageom)
{
;
}
RefinementSTLGeometry :: ~RefinementSTLGeometry ()
{
;
}
void RefinementSTLGeometry ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
newp = p1+secpoint*(p2-p1);
/*
(*testout) << "surf-between: p1 = " << p1 << ", p2 = " << p2
<< ", gi = " << gi1 << " - " << gi2 << endl;
*/
if (gi1.trignum > 0)
{
// ((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
Point<3> np1 = newp;
Point<3> np2 = newp;
((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
int tn1 = geom.Project (np1);
((STLGeometry&)geom).SelectChartOfTriangle (gi2.trignum);
int tn2 = geom.Project (np2);
newgi.trignum = tn1; //urspruengliche version
newp = np1; //urspruengliche version
if (!newgi.trignum)
{ newgi.trignum = tn2; newp = np2; }
if (!newgi.trignum) newgi.trignum = gi1.trignum;
/*
if (tn1 != 0 && tn2 != 0 && ((STLGeometry&)geom).GetAngle(tn1,tn2) < M_PI*0.05) {
newgi.trignum = tn1;
newp = np1;
}
else
{
newp = ((STLGeometry&)geom).PointBetween(p1, gi1.trignum, p2, gi2.trignum);
tn1 = ((STLGeometry&)geom).Project(newp);
newgi.trignum = tn1;
if (!tn1)
{
newp = Center (p1, p2);
newgi.trignum = 0;
}
}
*/
}
else
{
// (*testout) << "WARNING: PointBetween got geominfo = 0" << endl;
newp = p1+secpoint*(p2-p1);
newgi.trignum = 0;
}
// (*testout) << "newp = " << newp << ", ngi = " << newgi << endl;
}
void RefinementSTLGeometry ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & gi1,
const EdgePointGeomInfo & gi2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
/*
(*testout) << "edge-between: p1 = " << p1 << ", p2 = " << p2
<< ", gi1,2 = " << gi1 << ", " << gi2 << endl;
*/
/*
newp = Center (p1, p2);
((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
newgi.trignum = geom.Project (newp);
*/
int hi;
newgi.dist = (1.0-secpoint) * gi1.dist + secpoint*gi2.dist;
newgi.edgenr = gi1.edgenr;
/*
(*testout) << "p1 = " << p1 << ", p2 = " << p2 << endl;
(*testout) << "refedge: " << gi1.edgenr
<< " d1 = " << gi1.dist << ", d2 = " << gi2.dist << endl;
*/
newp = geom.GetLine (gi1.edgenr)->GetPointInDist (geom.GetPoints(), newgi.dist, hi);
// (*testout) << "newp = " << newp << endl;
}
void RefinementSTLGeometry :: ProjectToSurface (Point<3> & p, int surfi) const
{
cout << "RefinementSTLGeometry :: ProjectToSurface not implemented!" << endl;
};
void RefinementSTLGeometry :: ProjectToSurface (Point<3> & p, int surfi,
PointGeomInfo & gi) const
{
((STLGeometry&)geom).SelectChartOfTriangle (gi.trignum);
gi.trignum = geom.Project (p);
// if (!gi.trignum)
// cout << "projectSTL failed" << endl;
};
}

54
libsrc/stlgeom/meshstlsurface.hpp
View File

@ -63,59 +63,5 @@ protected:
double Area () const override;
};
///
class MeshOptimizeSTLSurface : public MeshOptimize2d
{
///
STLGeometry & geom;
public:
///
MeshOptimizeSTLSurface (STLGeometry & ageom);
///
void ProjectPoint (INDEX surfind, Point<3> & p) const override;
///
int ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const override;
///
void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point<3> & p) const override;
///
int CalcPointGeomInfo(PointGeomInfo& gi, const Point<3> & p3) const override;
///
void GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const override;
void GetNormalVector(INDEX surfind, const Point<3> & p, PointGeomInfo & gi, Vec<3> & n) const override;
};
class RefinementSTLGeometry : public Refinement
{
const STLGeometry & geom;
public:
RefinementSTLGeometry (const STLGeometry & ageom);
virtual ~RefinementSTLGeometry ();
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
virtual void ProjectToSurface (Point<3> & p, int surfi) const override;
virtual void ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi) const override;
};
#endif

126
libsrc/stlgeom/stlgeom.cpp
View File

@ -44,7 +44,6 @@ void STLMeshing (STLGeometry & geom,
lineendpoints(), spiralpoints(), selectedmultiedge()
*/
{
ref = NULL;
edgedata = make_unique<STLEdgeDataList>(*this);
externaledges.SetSize(0);
Clear();
@ -66,7 +65,6 @@ STLGeometry :: ~STLGeometry()
{
// for (auto p : atlas) delete p;
// delete edgedata;
delete ref;
}
void STLGeometry :: Save (string filename) const
@ -102,17 +100,127 @@ int STLGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mp
return STLMeshingDummy (this, mesh, mparam, stlpar);
}
const Refinement & STLGeometry :: GetRefinement () const
Vec<3> STLGeometry :: GetNormal(INDEX surfind, const Point<3> & p) const
{
delete ref;
ref = new RefinementSTLGeometry(*this);
// ref -> Set2dOptimizer(new MeshOptimizeSTLSurface(*this)); ??? copied from CSG
return *ref;
throw Exception("STLGeometry::GetNormal without PointGeomInfo called");
}
Vec<3> STLGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const
{
return GetChart(GetChartNr(gi.trignum)).GetNormal();
}
bool STLGeometry :: CalcPointGeomInfo(int /*surfind*/, PointGeomInfo& gi, const Point<3> & p3) const
{
Point<3> hp = p3;
SelectChartOfTriangle(gi.trignum);
gi.trignum = Project (hp);
if (gi.trignum) return true;
return false;
}
bool STLGeometry :: ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const
{
static std::mutex mutex_project_whole_surface;
int meshchart = GetChartNr(gi.trignum);
const STLChart& chart = GetChart(meshchart);
int trignum = chart.ProjectNormal(p);
if(trignum==0)
{
// non-thread-safe implementation
std::lock_guard<std::mutex> guard(mutex_project_whole_surface);
PrintMessage(7,"project failed");
SelectChartOfTriangle (gi.trignum); // needed because ProjectOnWholeSurface uses meshchartnv (the normal vector of selected chart)
trignum = ProjectOnWholeSurface(p);
if(trignum==0)
{
PrintMessage(7, "project on whole surface failed");
return false;
}
}
return true;
}
void STLGeometry :: ProjectPoint (INDEX surfind, Point<3> & p) const
{
throw Exception("ProjectPoint without PointGeomInfo not implemented");
}
void STLGeometry ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
newp = p1+secpoint*(p2-p1);
/*
(*testout) << "surf-between: p1 = " << p1 << ", p2 = " << p2
<< ", gi = " << gi1 << " - " << gi2 << endl;
*/
if (gi1.trignum > 0)
{
// ((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
Point<3> np1 = newp;
Point<3> np2 = newp;
auto ngi1 = gi1;
auto ngi2 = gi2;
// SelectChartOfTriangle (gi1.trignum);
int tn1 = ProjectPointGI (surfi, np1, ngi1);
// SelectChartOfTriangle (gi2.trignum);
int tn2 = ProjectPointGI (surfi, np2, ngi2);
newgi.trignum = tn1; //urspruengliche version
newp = np1; //urspruengliche version
if (!newgi.trignum)
{ newgi.trignum = tn2; newp = np2; }
if (!newgi.trignum) newgi.trignum = gi1.trignum;
}
else
{
// (*testout) << "WARNING: PointBetween got geominfo = 0" << endl;
newp = p1+secpoint*(p2-p1);
newgi.trignum = 0;
}
}
void STLGeometry ::
PointBetweenEdge (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & gi1,
const EdgePointGeomInfo & gi2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
/*
(*testout) << "edge-between: p1 = " << p1 << ", p2 = " << p2
<< ", gi1,2 = " << gi1 << ", " << gi2 << endl;
*/
/*
newp = Center (p1, p2);
((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
newgi.trignum = geom.Project (newp);
*/
int hi;
newgi.dist = (1.0-secpoint) * gi1.dist + secpoint*gi2.dist;
newgi.edgenr = gi1.edgenr;
/*
(*testout) << "p1 = " << p1 << ", p2 = " << p2 << endl;
(*testout) << "refedge: " << gi1.edgenr
<< " d1 = " << gi1.dist << ", d2 = " << gi2.dist << endl;
*/
newp = GetLine (gi1.edgenr)->GetPointInDist (GetPoints(), newgi.dist, hi);
// (*testout) << "newp = " << newp << endl;
}
void STLGeometry :: STLInfo(double* data)
{

28
libsrc/stlgeom/stlgeom.hpp
View File

@ -148,7 +148,7 @@ namespace netgen
//for meshing and project:
NgArray<int> meshcharttrigs; //per trig: 1=belong to chart, 0 not
int meshchart;
mutable int meshchart;
NgArray<int> ha_points; // help array, np long, filled with 0
@ -159,12 +159,10 @@ namespace netgen
//transformation:
Vec<3> meshtrignv;
mutable Vec<3> meshtrignv;
Vec<3> ex, ey, ez;
Point<3> p1;
mutable class RefinementSTLGeometry * ref;
public:
int edgesfound;
int surfacemeshed;
@ -194,6 +192,24 @@ namespace netgen
virtual void Save (string filename) const override;
bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override;
void ProjectPoint(INDEX surfind, Point<3> & p) const override;
bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override;
Vec<3> GetNormal(int surfind, const Point<3> & p) const override;
Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const override;
void PointBetween(const Point<3> & p1, const Point<3> & p2,
double secpoint, int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
DLL_HEADER void STLInfo(double* data);
//stldoctor:
@ -419,7 +435,7 @@ namespace netgen
//
void DefineTangentialPlane(const Point<3> & ap1, const Point<3> & ap2, int trig);
//
void SelectChartOfTriangle (int trignum);
void SelectChartOfTriangle (int trignum) const;
//
void SelectChartOfPoint (const Point<3> & p);
//
@ -459,8 +475,6 @@ namespace netgen
int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam) override;
virtual const Refinement & GetRefinement () const override;
// Add additional Point to chart to close the surface and write the resulting stl to a file
DLL_HEADER void WriteChartToFile( ChartId chartnumber, string filename="chart.slb" );
};

2
libsrc/stlgeom/stlgeommesh.cpp
View File

@ -392,7 +392,7 @@ void STLGeometry :: DefineTangentialPlane (const Point<3> & ap1, const Point<3>
}
void STLGeometry :: SelectChartOfTriangle (int trignum)
void STLGeometry :: SelectChartOfTriangle (int trignum) const
{
meshchart = GetChartNr(trignum);
meshtrignv = GetTriangle(trignum).Normal();

1
ng/ngpkg.cpp
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;

3
nglib/CMakeLists.txt
View File

@ -7,7 +7,6 @@ if(WIN32)
$<TARGET_OBJECTS:interface>
$<TARGET_OBJECTS:geom2d>
$<TARGET_OBJECTS:csg>
$<TARGET_OBJECTS:stl>
$<TARGET_OBJECTS:visual>
$<TARGET_OBJECTS:occ>
@ -23,7 +22,7 @@ endif(WIN32)
add_library(nglib SHARED nglib.cpp ${nglib_objects})
if(NOT WIN32)
target_link_libraries( nglib PUBLIC mesh stl interface geom2d csg stl visual)
target_link_libraries( nglib PUBLIC mesh interface geom2d csg stl visual)
if(USE_GUI)
target_link_libraries( nglib PUBLIC stlvis geom2dvis csgvis )
endif(USE_GUI)

10
nglib/nglib.cpp
View File

@ -536,7 +536,7 @@ namespace nglib
Ng_Mesh * mesh,
int levels)
{
Refinement2d ref(*(SplineGeometry2d*)geom);
Refinement ref(*(SplineGeometry2d*)geom);
HPRefinement (*(Mesh*)mesh, &ref, levels);
}
@ -547,7 +547,7 @@ namespace nglib
Ng_Mesh * mesh,
int levels, double parameter)
{
Refinement2d ref(*(SplineGeometry2d*)geom);
Refinement ref(*(SplineGeometry2d*)geom);
HPRefinement (*(Mesh*)mesh, &ref, levels, parameter);
}
@ -1090,7 +1090,7 @@ namespace nglib
// ------------------ Begin - Second Order Mesh generation functions ----------------
DLL_HEADER void Ng_Generate_SecondOrder(Ng_Mesh * mesh)
{
Refinement ref;
Refinement ref(*((Mesh*) mesh)->GetGeometry());
ref.MakeSecondOrder(*(Mesh*) mesh);
}
@ -1139,8 +1139,8 @@ namespace nglib
// ------------------ Begin - Uniform Mesh Refinement functions ---------------------
DLL_HEADER void Ng_Uniform_Refinement (Ng_Mesh * mesh)
{
Refinement ref;
ref.Refine ( * (Mesh*) mesh );
Refinement ref(*((Mesh*)mesh)->GetGeometry());
ref.Refine ( * (Mesh*) mesh );
}

1340
tests/pytest/results.json
View File

File diff suppressed because it is too large Load Diff

2
tests/pytest/test_tutorials.py
View File

@ -28,7 +28,7 @@ def checkData(mesh, mp, ref):
assert ref['ne1d'] == data['ne1d']
assert ref['ne2d'] == data['ne2d']
assert ref['ne3d'] == data['ne3d']
assert ref['quality_histogram'] == data['quality_histogram']
assert json.loads(ref['quality_histogram']) == pytest.approx(json.loads(data['quality_histogram']), abs=1, rel=0.4)
assert ref['total_badness'] == pytest.approx(data['total_badness'], rel=1e-5)
# get tutorials