refactor a lot of the old code, stl still needs to be done

This commit is contained in:
Christopher Lackner 2019-10-02 20:29:18 +02:00 committed by Matthias Hochsteger
parent 43cc5e68b1
commit 05881c0eb5
36 changed files with 644 additions and 1173 deletions

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;

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++; }

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);
}
}
}

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);
}
}

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

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
)

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);
}
}

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

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:

View File

@ -13,7 +13,6 @@
// #include "../gprim/spline.hpp"
// #include "../gprim/splinegeometry.hpp"
#include "geom2dmesh.hpp"
namespace netgen
{
@ -151,12 +150,35 @@ 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 +190,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 +227,6 @@ namespace netgen
int AddBCName (string name);
string * BCNamePtr ( const int bcnr );
DLL_HEADER virtual Refinement & GetRefinement () const;
};
}

View File

@ -7,7 +7,7 @@ add_library(interface ${NG_LIB_TYPE}
wuchemnitz.cpp writegmsh2.cpp writeOpenFOAM15x.cpp
)
target_link_libraries(interface mesh csg geom2d stl visual)
target_link_libraries(interface mesh csg geom2d visual)
if(NOT WIN32)
install( TARGETS interface ${NG_INSTALL_DIR})

View File

@ -17,7 +17,7 @@ namespace netgen
static Timer timer_opt2d("Optimization 2D");
RegionTimer reg(timer_opt2d);
auto meshopt = GetMeshOptimizer();
auto meshopt = MeshOptimize2d(mesh);
for(auto i : Range(mparam.optsteps2d))
{
PrintMessage(2, "Optimization step ", i);
@ -26,16 +26,16 @@ namespace netgen
switch(optstep)
{
case 's':
meshopt->EdgeSwapping(mesh, 0);
meshopt.EdgeSwapping(0);
break;
case 'S':
meshopt->EdgeSwapping(mesh, 1);
meshopt.EdgeSwapping(1);
break;
case 'm':
meshopt->ImproveMesh(mesh, mparam);
meshopt.ImproveMesh(mparam);
break;
case 'c':
meshopt->CombineImprove (mesh);
meshopt.CombineImprove();
break;
}
}
@ -129,7 +129,7 @@ namespace netgen
const Refinement & NetgenGeometry :: GetRefinement () const
{
return *new Refinement;;
return *new Refinement(*this);
}

View File

@ -31,9 +31,57 @@ namespace netgen
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 unique_ptr<MeshOptimize2d> GetMeshOptimizer() const
{ return make_unique<MeshOptimize2d>(); }
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 { ; }
};

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;
}
}

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,49 +51,10 @@ 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 {;}
};

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,15 +1127,15 @@ 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],
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]],

View File

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

View File

@ -39,7 +39,7 @@ namespace netgen
void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
void MeshOptimize2d :: EdgeSwapping (int usemetric)
{
static Timer timer("EdgeSwapping (2D)"); RegionTimer reg(timer);
if (!faceindex)
@ -51,7 +51,7 @@ namespace netgen
for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
{
EdgeSwapping (mesh, usemetric);
EdgeSwapping (usemetric);
if (multithread.terminate)
throw NgException ("Meshing stopped");
@ -81,7 +81,7 @@ namespace netgen
for (SurfaceElementIndex sei : seia)
if (mesh[sei].GetNP() != 3)
{
GenericImprove (mesh);
GenericImprove();
return;
}
@ -317,14 +317,13 @@ namespace netgen
nv1.Normalize();
nv2.Normalize();
Vec<3> nvp3, nvp4;
SelectSurfaceOfPoint (mesh.Point(pi3), gi3);
GetNormalVector (surfnr, mesh.Point(pi3), gi3, nvp3);
// SelectSurfaceOfPoint (mesh.Point(pi3), gi3);
auto nvp3 = geo.GetNormal(surfnr, mesh.Point(pi3), gi3);
nvp3.Normalize();
SelectSurfaceOfPoint (mesh.Point(pi4), gi4);
GetNormalVector (surfnr, mesh.Point(pi4), gi4, nvp4);
// SelectSurfaceOfPoint (mesh.Point(pi4), gi4);
auto nvp4 = geo.GetNormal(surfnr, mesh.Point(pi4), gi4);
nvp4.Normalize();
@ -426,16 +425,16 @@ namespace netgen
void MeshOptimize2d :: CombineImprove (Mesh & mesh)
void MeshOptimize2d :: CombineImprove()
{
if (!faceindex)
{
SplitImprove(mesh);
SplitImprove();
PrintMessage (3, "Combine improve");
for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
{
CombineImprove (mesh);
CombineImprove();
if (multithread.terminate)
throw NgException ("Meshing stopped");
@ -530,8 +529,8 @@ namespace netgen
for (int k = 0; k < 3; k++)
if (hel[k] == pi)
{
SelectSurfaceOfPoint (mesh[pi], hel.GeomInfoPi(k+1));
GetNormalVector (surfnr, mesh[pi], hel.GeomInfoPi(k+1), normals[pi]);
// SelectSurfaceOfPoint (mesh[pi], hel.GeomInfoPi(k+1));
normals[pi] = geo.GetNormal(surfnr, mesh[pi], hel.GeomInfoPi(k+1));
break;
}
}
@ -624,9 +623,9 @@ namespace netgen
for (int k = 0; k < 3; k++)
if (hel[k] == pi1)
{
SelectSurfaceOfPoint (mesh[pi1],
hel.GeomInfoPi(k+1));
GetNormalVector (surfnr, mesh[pi1], hel.GeomInfoPi(k+1), nv);
// SelectSurfaceOfPoint (mesh[pi1],
// hel.GeomInfoPi(k+1));
nv = geo.GetNormal(surfnr, mesh[pi1], hel.GeomInfoPi(k+1));
break;
}
@ -794,7 +793,7 @@ namespace netgen
mesh.SetNextTimeStamp();
}
void MeshOptimize2d :: SplitImprove (Mesh & mesh)
void MeshOptimize2d :: SplitImprove()
{
if (!faceindex)
{
@ -803,7 +802,7 @@ namespace netgen
mesh.CalcSurfacesOfNode(); // TODO: needed?
for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
{
SplitImprove (mesh);
SplitImprove();
if (multithread.terminate)
throw NgException ("Meshing stopped");
@ -909,9 +908,9 @@ namespace netgen
PointIndex pi5;
PointGeomInfo gi5;
mesh.GetGeometry()->GetRefinement().PointBetween (mesh[pi1], mesh[pi2], 0.5,
faceindex,
gi1, gi2, p5, gi5);
geo.PointBetween(mesh[pi1], mesh[pi2], 0.5,
faceindex,
gi1, gi2, p5, gi5);
pi5 = mesh.AddPoint(p5);

View File

@ -6,27 +6,29 @@
///
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);
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; }
@ -35,31 +37,9 @@ public:
void SetWriteStatus (int ws) { writestatus = ws; }
///
virtual void SelectSurfaceOfPoint (const Point<3> & p,
const PointGeomInfo & gi);
///
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);

View File

@ -19,7 +19,7 @@ namespace netgen
};
void MeshOptimize2d :: GenericImprove (Mesh & mesh)
void MeshOptimize2d :: GenericImprove ()
{
if (!faceindex)
{
@ -27,7 +27,7 @@ namespace netgen
PrintMessage (3, "Generic Improve");
for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
GenericImprove (mesh);
GenericImprove ();
faceindex = 0;
}
@ -395,10 +395,8 @@ namespace netgen
// calc metric badness
double bad1 = 0, bad2 = 0;
Vec<3> n;
SelectSurfaceOfPoint (mesh.Point(pmap.Get(1)), pgi.Get(1));
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);

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:

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

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

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;

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,30 +688,7 @@ namespace netgen
return badness;
}
MeshOptimize2d dummy;
MeshOptimize2d :: MeshOptimize2d ()
{
SetFaceIndex (0);
SetImproveEdges (0);
SetMetricWeight (0);
SetWriteStatus (1);
}
void MeshOptimize2d :: SelectSurfaceOfPoint (const Point<3> & p,
const PointGeomInfo & gi)
{
;
}
void MeshOptimize2d :: ImproveMesh (Mesh & mesh, const MeshingParameters & mp)
void MeshOptimize2d :: ImproveMesh (const MeshingParameters & mp)
{
if (!faceindex)
{
@ -721,7 +696,7 @@ namespace netgen
for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
{
ImproveMesh (mesh, mp);
ImproveMesh (mp);
if (multithread.terminate)
throw NgException ("Meshing stopped");
}
@ -959,7 +934,7 @@ namespace netgen
}
ld.gi1 = hel.GeomInfoPi(hpi);
SelectSurfaceOfPoint (ld.sp1, ld.gi1);
// SelectSurfaceOfPoint (ld.sp1, ld.gi1);
ld.locelements.SetSize(0);
ld.locrots.SetSize (0);
@ -992,7 +967,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);
@ -1070,7 +1045,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)
@ -1094,14 +1069,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);
}
}

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]);
}
}
}

View File

@ -941,41 +941,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);
}
}

View File

@ -1034,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;
@ -1056,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
{
@ -1115,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)
@ -1705,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

View File

@ -273,8 +273,6 @@ namespace netgen
const MeshingParameters& mparam) override;
void MeshSurface(Mesh& mesh,
const MeshingParameters& mparam) override;
unique_ptr<MeshOptimize2d> GetMeshOptimizer() const override
{ return make_unique<MeshOptimize2dOCCSurfaces>(*this); }
void FinalizeMesh(Mesh& mesh) const override;
@ -282,6 +280,24 @@ namespace netgen
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
@ -301,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;
@ -429,8 +442,8 @@ namespace netgen
// void WriteOCC_STL(char * filename);
// DLL_HEADER virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
DLL_HEADER const Refinement & GetRefinement () const override;
private:
bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
};

View File

@ -463,188 +463,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,
@ -703,76 +521,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);
}
};
}

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

View File

@ -299,15 +299,15 @@ 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.EdgeSwapping(0);
mesh.CalcSurfacesOfNode();
optmesh.ImproveMesh (mesh, mparam);
optmesh.ImproveMesh(mparam);
}
mesh.Compress();
@ -827,7 +827,7 @@ void STLSurfaceOptimization (STLGeometry & geom,
{
PrintFnStart("optimize STL Surface");
MeshOptimizeSTLSurface optmesh(geom);
MeshOptimize2d optmesh(mesh);
optmesh.SetFaceIndex (0);
optmesh.SetImproveEdges (0);
@ -848,22 +848,22 @@ 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;
}
}
@ -881,7 +881,7 @@ void STLSurfaceOptimization (STLGeometry & geom,
}
}
}
optmesh.SplitImprove(mesh);
optmesh.SplitImprove();
}
//(*testout) << "optimize, after, step = " << meshparam.optimize2d[j-1] << mesh.Point (3679) << endl;
}
@ -1068,18 +1068,6 @@ double MeshingSTLSurface :: Area () const
return geom.Area();
}
MeshOptimizeSTLSurface :: MeshOptimizeSTLSurface (STLGeometry & ageom)
: MeshOptimize2d(), geom(ageom)
{
;
}
void MeshOptimizeSTLSurface :: SelectSurfaceOfPoint (const Point<3> & p,
const PointGeomInfo & gi)
{

View File

@ -63,59 +63,5 @@ protected:
double Area () const override;
};
///
class MeshOptimizeSTLSurface : public MeshOptimize2d
{
///
STLGeometry & geom;
public:
///
MeshOptimizeSTLSurface (STLGeometry & ageom);
///
virtual void SelectSurfaceOfPoint (const Point<3> & p,
const PointGeomInfo & gi);
///
virtual void ProjectPoint (INDEX surfind, Point<3> & p) const;
///
virtual void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point<3> & p) const;
///
virtual int CalcPointGeomInfo(PointGeomInfo& gi, const Point<3> & p3) const;
///
virtual void GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const;
};
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

View File

@ -66,7 +66,6 @@ STLGeometry :: ~STLGeometry()
{
// for (auto p : atlas) delete p;
// delete edgedata;
delete ref;
}
void STLGeometry :: Save (string filename) const
@ -102,18 +101,6 @@ int STLGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mp
return STLMeshingDummy (this, mesh, mparam, stlpar);
}
const Refinement & STLGeometry :: GetRefinement () const
{
delete ref;
ref = new RefinementSTLGeometry(*this);
// ref -> Set2dOptimizer(new MeshOptimizeSTLSurface(*this)); ??? copied from CSG
return *ref;
}
void STLGeometry :: STLInfo(double* data)
{
data[0] = GetNT();

View File

@ -3,18 +3,15 @@ add_definitions(-DNGLIB_EXPORTS)
if(WIN32)
set(nglib_objects
$<TARGET_OBJECTS:mesh>
$<TARGET_OBJECTS:stl>
$<TARGET_OBJECTS:interface>
$<TARGET_OBJECTS:geom2d>
$<TARGET_OBJECTS:csg>
$<TARGET_OBJECTS:stl>
$<TARGET_OBJECTS:visual>
$<TARGET_OBJECTS:occ>
)
if(USE_GUI)
set(nglib_objects ${nglib_objects}
$<TARGET_OBJECTS:stlvis>
$<TARGET_OBJECTS:geom2dvis>
$<TARGET_OBJECTS:csgvis>
)
@ -23,9 +20,9 @@ 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 visual)
if(USE_GUI)
target_link_libraries( nglib PUBLIC stlvis geom2dvis csgvis )
target_link_libraries( nglib PUBLIC geom2dvis csgvis )
endif(USE_GUI)
endif(NOT WIN32)

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 );
}