ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-24 21:10:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' into boundarylayers

Browse Source
This commit is contained in:
Christopher Lackner 2020-11-11 12:09:59 +01:00
commit 7beb82af04
7 changed files with 206 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
libsrc/geom2d/csg2d.cpp
View File

@ -1292,20 +1292,38 @@ void CreateResult(Solid2d & sp, Solid2d & sr, bool UNION)
}
}
// Check if vertex v is not necessary (i.e. is on the line v->prev, v->next and has same info as v->prev and no pinfo
bool canRemoveVertex( Vertex * v )
{
return false;
if(v->spline)
return false;
if(v->pinfo.name != POINT_NAME_DEFAULT)
return false;
if(v->pinfo.maxh != MAXH_DEFAULT)
return false;
if(v->info.bc != v->prev->info.bc || v->info.maxh != v->prev->info.maxh )
return false;
if(fabs(Area(*v->prev,*v,*v->next)) >= EPSILON)
return false;
return true;
}
void CleanUpResult(Solid2d & sr)
{
auto & RR = sr.polys;
for (Loop& R : RR)
{
while ( (R.first.get() != NULL) && (fabs(Area(*R.first->prev,*R.first,*R.first->next)) < EPSILON) )
while ( (R.first.get() != NULL) && canRemoveVertex(R.first.get()))
R.Remove(R.first.get());
if (R.first.get() != NULL)
for (Vertex* V : R.Vertices(ALL))
if (!V->spline && !V->prev->spline && fabs(Area(*V->prev,*V,*V->next)) < EPSILON)
{
if (canRemoveVertex(V))
R.Remove(V);
}
}
for (int i = RR.Size()-1; i>=0; i--)
if(RR[i].Size()==0)
@ -1547,6 +1565,39 @@ Solid2d ClipSolids ( Solid2d && s1, Solid2d && s2, char op)
return std::move(res);
}
Vertex* Loop :: getNonIntersectionVertex()
{
for (Vertex* v : Vertices(ALL))
if (!v->is_intersection)
return(v);
// no non-intersection vertex found -> generate and return temporary vertex
for (Vertex* v : Vertices(ALL))
// make sure that edge from V to V->next is not collinear with other polygon
if ( (v->next->neighbour != v->neighbour->prev) && (v->next->neighbour != v->neighbour->next) )
{
// add edge midpoint as temporary vertex
if(v->spline)
{
auto p = v->spline->GetPoint(0.5);
auto s = *v->spline;
v->spline = Split(s, 0, 0.5);
auto vnew = v->Insert(p);
vnew->info = v->info;
vnew->spline = Split(s, 0.5, 1.0);
return vnew;
}
else
{
auto p = Center(*v, *v->next);
auto vnew = v->Insert(p);
vnew->info = v->info;
return vnew;
}
}
return(NULL);
}
bool Loop :: IsInside( Point<2> r ) const
{
int w = 0;
@ -1652,11 +1703,14 @@ Solid2d & Solid2d :: Move( Vec<2> v )
return Transform( [v](Point<2> p) -> Point<2> { return p+v; } );
}
Solid2d & Solid2d :: Scale( double sx, double sy )
Solid2d & Solid2d :: Scale( double s )
{
if(sy==0.0)
sy=sx;
return Transform( [sx,sy](Point<2> p) -> Point<2> { return{p[0]*sx, p[1]*sy}; } );
return Transform( [s](Point<2> p) -> Point<2> { return{p[0]*s, p[1]*s}; } );
}
Solid2d & Solid2d :: Scale( Vec<2> s )
{
return Transform( [s](Point<2> p) -> Point<2> { return{p[0]*s[0], p[1]*s[1]}; } );
}
Solid2d & Solid2d :: RotateRad( double ang, Point<2> center )
@ -1669,8 +1723,8 @@ Solid2d & Solid2d :: RotateRad( double ang, Point<2> center )
p -= c;
double x = p[0];
double y = p[1];
p[0] = cosa*x+sina*y;
p[1] = -sina*x+cosa*y;
p[0] = cosa*x-sina*y;
p[1] = sina*x+cosa*y;
p += c;
return p;
} );

31
libsrc/geom2d/csg2d.hpp
View File

@ -29,10 +29,10 @@ enum IntersectionType
T_INTERSECTION_Q,
T_INTERSECTION_P,
V_INTERSECTION,
X_OVERLAP,
T_OVERLAP_Q,
T_OVERLAP_P,
V_OVERLAP
X_OVERLAP, // Q0 -- P1 -- Q1 -- P0 (different direction)
T_OVERLAP_Q, // same direction or P inside Q
T_OVERLAP_P, // same direction or Q inside P
V_OVERLAP // one common point
};
enum IntersectionLabel
@ -588,23 +588,7 @@ struct Loop
//
// return and insert a non-intersection vertex
//
Vertex* getNonIntersectionVertex()
{
for (Vertex* v : Vertices(ALL))
if (!v->is_intersection)
return(v);
// no non-intersection vertex found -> generate and return temporary vertex
for (Vertex* v : Vertices(ALL))
// make sure that edge from V to V->next is not collinear with other polygon
if ( (v->next->neighbour != v->neighbour->prev) && (v->next->neighbour != v->neighbour->next) )
{
// add edge midpoint as temporary vertex
auto p = Center(*v, *v->next);
return v->Insert(p);
}
return(NULL);
}
Vertex* getNonIntersectionVertex();
void SetBC(string bc)
{
@ -694,11 +678,12 @@ struct Solid2d
}
Solid2d & Move( Vec<2> v );
Solid2d & Scale( double sx, double sy=0.0 );
Solid2d & Scale( double s );
Solid2d & Scale( Vec<2> s );
Solid2d & RotateRad( double ang, Point<2> center = {0,0} );
Solid2d & RotateDeg( double ang, Point<2> center = {0,0} )
{
return RotateRad( ang/180.*M_PI );
return RotateRad( ang/180.*M_PI, center );
}
Solid2d & BC(string bc)

13
libsrc/geom2d/python_geom2d.cpp
View File

@ -414,16 +414,9 @@ DLL_HEADER void ExportGeom2d(py::module &m)
.def("Copy", [](Solid2d & self) -> Solid2d { return self; })
.def("Move", &Solid2d::Move)
.def("Scale", &Solid2d::Scale)
.def("Rotate", [](Solid2d & self, optional<double> rad, optional<double> deg, Point<2> center )
{
if(rad)
self.RotateRad(*rad, center);
if(deg)
self.RotateDeg(*deg, center);
return self;
}, py::arg("rad")=nullopt, py::arg("deg")=nullopt, py::arg("center")=Point<2>{0,0})
.def("Scale", static_cast<Solid2d& (Solid2d::*)(double)>(&Solid2d::Scale))
.def("Scale", static_cast<Solid2d& (Solid2d::*)(Vec<2>)>(&Solid2d::Scale))
.def("Rotate", &Solid2d::RotateDeg, py::arg("angle"), py::arg("center")=Point<2>{0,0})
;

30
libsrc/meshing/python_mesh.cpp
View File

@ -211,6 +211,10 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
py::class_<Vec<2>> (m, "Vec2d")
.def(py::init<double,double>())
.def(py::init( [] (std::pair<double,double> xy)
{
return Vec<2>{xy.first, xy.second};
}))
.def ("__str__", &ToString<Vec<3>>)
.def(py::self==py::self)
.def(py::self+py::self)
@ -223,6 +227,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
.def("__len__", [](Vec<2>& /*unused*/) { return 2; })
;
py::implicitly_convertible<py::tuple, Vec<2>>();
py::class_<Vec<3>> (m, "Vec3d")
.def(py::init<double,double,double>())
.def(py::init([](py::tuple v)
@ -1261,27 +1267,45 @@ project_boundaries : Optional[str] = None
}), py::arg("mapping"))
.def(NGSPickle<SurfaceGeometry>())
.def("GenerateMesh", [](shared_ptr<SurfaceGeometry> geo,
bool quads, int nx, int ny, bool flip_triangles, py::list py_bbbpts, py::list py_bbbnames)
bool quads, int nx, int ny, bool flip_triangles, py::list py_bbbpts, py::list py_bbbnames, py::list py_hppts, py::list py_hpbnd)
{
if (py::len(py_bbbpts) != py::len(py_bbbnames))
throw Exception("In SurfaceGeometry::GenerateMesh bbbpts and bbbnames do not have same lengths.");
Array<Point<3>> bbbpts(py::len(py_bbbpts));
Array<string> bbbname(py::len(py_bbbpts));
Array<Point<3>> hppts(py::len(py_hppts));
Array<float> hpptsfac(py::len(py_hppts));
Array<string> hpbnd(py::len(py_hpbnd));
Array<float> hpbndfac(py::len(py_hpbnd));
for(int i = 0; i<py::len(py_bbbpts);i++)
{
py::tuple pnt = py::extract<py::tuple>(py_bbbpts[i])();
bbbpts[i] = Point<3>(py::extract<double>(pnt[0])(),py::extract<double>(pnt[1])(),py::extract<double>(pnt[2])());
bbbname[i] = py::extract<string>(py_bbbnames[i])();
}
for(int i = 0; i<py::len(py_hppts);i++)
{
py::tuple pnt = py::extract<py::tuple>(py_hppts[i])();
hppts[i] = Point<3>(py::extract<double>(pnt[0])(),py::extract<double>(pnt[1])(),py::extract<double>(pnt[2])());
//hpptsfac[i] = py::len(pnt) > 3 ? py::extract<double>(pnt[3])() : 0.0;
hpptsfac[i] = py::extract<double>(pnt[3])();
}
for(int i = 0; i<py::len(py_hpbnd);i++)
{
py::tuple bnd = py::extract<py::tuple>(py_hpbnd[i])();
hpbnd[i] = py::extract<string>(bnd[0])();
hpbndfac[i] = py::extract<double>(bnd[1])();
}
auto mesh = make_shared<Mesh>();
SetGlobalMesh (mesh);
mesh->SetGeometry(geo);
ng_geometry = geo;
auto result = geo->GenerateMesh (mesh, quads, nx, ny, flip_triangles, bbbpts, bbbname);
auto result = geo->GenerateMesh (mesh, quads, nx, ny, flip_triangles, bbbpts, bbbname, hppts, hpptsfac, hpbnd, hpbndfac);
if(result != 0)
throw Exception("SurfaceGeometry: Meshing failed!");
return mesh;
}, py::arg("quads")=true, py::arg("nx")=10, py::arg("ny")=10, py::arg("flip_triangles")=false, py::arg("bbbpts")=py::list(), py::arg("bbbnames")=py::list())
}, py::arg("quads")=true, py::arg("nx")=10, py::arg("ny")=10, py::arg("flip_triangles")=false, py::arg("bbbpts")=py::list(), py::arg("bbbnames")=py::list(), py::arg("hppts")=py::list(), py::arg("hpbnd")=py::list())
;
;

127
libsrc/meshing/surfacegeom.cpp
View File

@ -45,10 +45,10 @@ namespace netgen
prr = p[1]+2*eps;
pll = p[1]-2*eps;
dr = GetTangentVectors(u, pr);
dl = GetTangentVectors(u, pl);
drr = GetTangentVectors(u, prr);
dll = GetTangentVectors(u, pll);
GetTangentVectors(u, pr, dr);
GetTangentVectors(u, pl, dl);
GetTangentVectors(u, prr, drr);
GetTangentVectors(u, pll, dll);
f_vv = (1.0/(12.0*eps)) * (8.0*dr[1]-8.0*dl[1]-drr[1]+dll[1]);
}
@ -74,11 +74,31 @@ namespace netgen
return tang;
}
void SurfaceGeometry :: GetTangentVectors(double u, double v, Array<Vec<3>>& tang) const
{
Point<2> pru = Point<2>(u+eps,v);
Point<2> plu = Point<2>(u-eps,v);
Point<2> prru = Point<2>(u+2*eps,v);
Point<2> pllu = Point<2>(u-2*eps,v);
Point<2> prv = Point<2>(u,v+eps);
Point<2> plv = Point<2>(u,v-eps);
Point<2> prrv = Point<2>(u,v+2*eps);
Point<2> pllv = Point<2>(u,v-2*eps);
tang[0] = 1/(12.0*eps)*( 8.0*func(pru) - 8.0*func(plu) - func(prru) + func(pllu) );
tang[1] = 1/(12.0*eps)*( 8.0*func(prv) - 8.0*func(plv) - func(prrv) + func(pllv) );
}
Vec<3> SurfaceGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const
{
Array<Vec<3>> tang = GetTangentVectors(gi->u, gi->v);
auto normal = Cross(tang[0], tang[1]);
return Cross(tang[0], tang[1]);
normal.Normalize();
return normal;
}
@ -97,14 +117,13 @@ namespace netgen
bool SurfaceGeometry :: ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const
{
Array<Vec<3>> tangs;
Array<Vec<3>> tangs(2);
Vec<3> diff, f_uu, f_vv, f_uv;
Vec<2> r, dx;
double norm_r, det, energy=0.0, new_energy=0.0, alpha=2.0,u=0.0,v=0.0;
double norm_r, det, energy=0.0, new_energy=0.0, alpha=2.0,u=0.0,v=0.0,maxerr=1e-16;
Mat<2,2> mat, inv;
int num=0, maxit=20;
double damping=0.2;
int num=0, maxit=25;
double damping=0.5;
//Solve minimization problem
// argmin_(u,v) 0.5*\| f(u,v)-p\|^2
@ -115,7 +134,7 @@ namespace netgen
do
{
num++;
tangs = GetTangentVectors(gi.u, gi.v);
GetTangentVectors(gi.u, gi.v,tangs);
diff = func(Point<2>(gi.u, gi.v)) - Vec<3>(p);
energy = diff.Length2();
r = Vec<2>( diff*tangs[0], diff*tangs[1] );
@ -146,12 +165,13 @@ namespace netgen
while (alpha > 1e-10 && new_energy > energy+1e-14);
if (alpha <= 1e-10)
throw Exception("In SurfaceGeometry::ProjectPointGI: Linesearch min alpha reached!");
gi.u = u;
gi.v = v;
}
while ( norm_r > 1e-12 && num < maxit);
while ( norm_r > maxerr && num < maxit);
//Stay in reference domain [0,1]^2
if (gi.u < 0 || gi.u > 1 || gi.v < 0 || gi.v > 1)
@ -179,19 +199,13 @@ namespace netgen
void SurfaceGeometry :: 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);
PointGeomInfo pgi;
pgi.u = ap1.u+secpoint*(ap2.u-ap1.u);
pgi.v = ap1.v+secpoint*(ap2.v-ap1.v);
ProjectPointGI(surfi1, newp, pgi);
newgi.u = pgi.u;
newgi.v = pgi.v;
newgi.u = ap1.u+secpoint*(ap2.u-ap1.u);
newgi.v = ap1.v+secpoint*(ap2.v-ap1.v);
newgi.edgenr = ap1.edgenr;
newgi.body = -1;
newgi.dist = -1.0;
newp = Point<3>(func(Point<2>(newgi.u, newgi.v)));
}
void SurfaceGeometry :: PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
@ -200,16 +214,16 @@ namespace netgen
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
newp = p1+secpoint*(p2-p1);
newgi.u = gi1.u+secpoint*(gi2.u-gi1.u);
newgi.v = gi1.v+secpoint*(gi2.v-gi1.v);
newgi.trignum = -1;
ProjectPointGI(surfi, newp, newgi);
newp = Point<3>(func(Point<2>(newgi.u, newgi.v)));
//newp = p1+secpoint*(p2-p1);
//ProjectPointGI(surfi, newp, newgi);
}
int SurfaceGeometry :: GenerateMesh(shared_ptr<Mesh> & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array<Point<3>>& bbbpts, const Array<string>& bbbnames)
int SurfaceGeometry :: GenerateMesh(shared_ptr<Mesh> & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array<Point<3>>& bbbpts, const Array<string>& bbbnames, const Array<Point<3>>& hppoints, const Array<float>& hpptsfac, const Array<string>& hpbnd, const Array<float>& hpbndfac)
{
mesh->SetDimension(3);
@ -241,6 +255,16 @@ namespace netgen
indbbbpts[k] = pids[pids.Size()-1];
}
}
for (int k = 0; k < hppoints.Size(); k++)
{
auto diff = pnt - hppoints[k];
if(diff.Length2() < 1e-14)
{
(*mesh)[pids[pids.Size()-1]].Singularity(hpptsfac[k]);
}
}
}
for (bool f : found)
@ -318,8 +342,19 @@ namespace netgen
Segment seg;
seg.si = 1;
seg.edgenr = 1;
seg.epgeominfo[0].edgenr = 1;
seg.epgeominfo[1].edgenr = 1;
seg.epgeominfo[0].edgenr = 0;
seg.epgeominfo[1].edgenr = 0;
//for hp refinement
seg.singedge_left = 0;
seg.singedge_right = 0;
for (size_t i=0; i < hpbnd.Size(); i++)
{
if (hpbnd[i] == "bottom")
{
seg.singedge_left = hpbndfac[i];
seg.singedge_right = hpbndfac[i];
}
}
// needed for codim2 in 3d
seg.edgenr = 1;
for(int i=0; i < nx; i++)
@ -341,6 +376,18 @@ namespace netgen
seg.si = 2;
seg.edgenr = 2;
seg.singedge_left = 0;
seg.singedge_right = 0;
for (size_t i=0; i < hpbnd.Size(); i++)
{
if (hpbnd[i] == "right")
{
seg.singedge_left = hpbndfac[i];
seg.singedge_right = hpbndfac[i];
}
}
for(int i=0; i<ny; i++)
{
seg[0] = pids[i*(nx+1)+nx];
@ -360,6 +407,18 @@ namespace netgen
seg.si = 3;
seg.edgenr = 3;
seg.singedge_left = 0;
seg.singedge_right = 0;
for (size_t i=0; i < hpbnd.Size(); i++)
{
if (hpbnd[i] == "top")
{
seg.singedge_left = hpbndfac[i];
seg.singedge_right = hpbndfac[i];
}
}
for(int i=0; i<nx; i++)
{
seg[0] = pids[ny*(nx+1)+i+1];
@ -379,6 +438,18 @@ namespace netgen
seg.si = 4;
seg.edgenr = 4;
seg.singedge_left = 0;
seg.singedge_right = 0;
for (size_t i=0; i < hpbnd.Size(); i++)
{
if (hpbnd[i] == "left")
{
seg.singedge_left = hpbndfac[i];
seg.singedge_right = hpbndfac[i];
}
}
for(int i=0; i<ny; i++)
{
seg[0] = pids[(i+1)*(nx+1)];

5
libsrc/meshing/surfacegeom.hpp
View File

@ -36,6 +36,9 @@ namespace netgen
Array<Vec<3>> GetTangentVectors(double u, double v) const;
void GetTangentVectors(double u, double v, Array<Vec<3>>& tang) const;
virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const override;
virtual PointGeomInfo ProjectPoint(int surfind, Point<3> & p) const override;
@ -59,7 +62,7 @@ namespace netgen
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
int GenerateMesh(shared_ptr<Mesh> & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array<Point<3>>& bbbpts, const Array<string>& bbbnames);
int GenerateMesh(shared_ptr<Mesh> & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array<Point<3>>& bbbpts, const Array<string>& bbbnames, const Array<Point<3>>& hppoints, const Array<float>& hpptsfac, const Array<string>& hpbnd, const Array<float>& hpbndfac);
};

2
tests/pytest/test_csg2d.py
View File

@ -5,7 +5,7 @@ from pytest import approx
def test_two_circles():
c1 = Circle(center=(0,0), radius=1)
c2 = c1.Rotate(deg=45)
c2 = c1.Rotate(45)
s = c1*c2
geo = CSG2d()
geo.Add(s)