mirror of
https://github.com/NGSolve/netgen.git
synced 2024-12-24 21:10:33 +05:00
Merge branch 'csg2d' into master
This commit is contained in:
commit
20b82ae7af
@ -24,7 +24,7 @@ void ComputeWeight( Spline & s, Point<2> p )
|
|||||||
double A = (p[1]-a[1])*(b[0]-p[0]) - (p[0]-a[0])*(b[1]-p[1]);
|
double A = (p[1]-a[1])*(b[0]-p[0]) - (p[0]-a[0])*(b[1]-p[1]);
|
||||||
double B = (p[1]-c[1])*(b[0]-p[0]) - (p[0]-c[0])*(b[1]-p[1]);
|
double B = (p[1]-c[1])*(b[0]-p[0]) - (p[0]-c[0])*(b[1]-p[1]);
|
||||||
double det = sqrt(-A*B);
|
double det = sqrt(-A*B);
|
||||||
double tt = (B-det)/(A+det);
|
double tt = fabs(A+det)<EPSILON ? 1 : (B-det)/(A+det);
|
||||||
auto v = b-p;
|
auto v = b-p;
|
||||||
int dim = fabs(v[0]) > fabs(v[1]) ? 0 : 1;
|
int dim = fabs(v[0]) > fabs(v[1]) ? 0 : 1;
|
||||||
double weight = fabs(tt*(p[dim]-a[dim])/v[dim] + 1.0/tt*(p[dim]-c[dim])/v[dim]);
|
double weight = fabs(tt*(p[dim]-a[dim])/v[dim] + 1.0/tt*(p[dim]-c[dim])/v[dim]);
|
||||||
@ -268,7 +268,7 @@ IntersectionType IntersectSplineSegment( const Spline & s, const Point<2> & r0,
|
|||||||
return ClassifyNonOverlappingIntersection(alpha, beta);
|
return ClassifyNonOverlappingIntersection(alpha, beta);
|
||||||
}
|
}
|
||||||
|
|
||||||
IntersectionType IntersectSplineSegment1( const Spline & s, const Point<2> & r0, const Point<2> & r1, double& alpha, double& beta )
|
IntersectionType IntersectSplineSegment1( const Spline & s, const Point<2> & r0, const Point<2> & r1, double& alpha, double& beta, bool first=false)
|
||||||
{
|
{
|
||||||
Point<2> p0 = s.StartPI();
|
Point<2> p0 = s.StartPI();
|
||||||
Point<2> p1 = s.TangentPoint();
|
Point<2> p1 = s.TangentPoint();
|
||||||
@ -310,11 +310,14 @@ IntersectionType IntersectSplineSegment1( const Spline & s, const Point<2> & r0,
|
|||||||
}
|
}
|
||||||
|
|
||||||
int choice = 0;
|
int choice = 0;
|
||||||
if(vtype[0]==NO_INTERSECTION && vtype[1]!=NO_INTERSECTION)
|
if(!first)
|
||||||
choice = 1;
|
{
|
||||||
|
if(vtype[0]==NO_INTERSECTION && vtype[1]!=NO_INTERSECTION)
|
||||||
|
choice = 1;
|
||||||
|
|
||||||
if(valpha[0] < alpha+EPSILON)
|
if(valpha[0] < alpha+EPSILON)
|
||||||
choice = 1;
|
choice = 1;
|
||||||
|
}
|
||||||
|
|
||||||
if(valpha[choice] < alpha+EPSILON)
|
if(valpha[choice] < alpha+EPSILON)
|
||||||
return NO_INTERSECTION;
|
return NO_INTERSECTION;
|
||||||
@ -332,26 +335,38 @@ bool IsOverlapping( Spline p, Spline s, double & alpha, double & beta, Intersect
|
|||||||
|
|
||||||
double lam0 = -1e3*EPSILON;
|
double lam0 = -1e3*EPSILON;
|
||||||
double lam1 = -1e3*EPSILON;
|
double lam1 = -1e3*EPSILON;
|
||||||
|
double lam2 = -1e3*EPSILON;
|
||||||
|
double lam3 = -1e3*EPSILON;
|
||||||
alpha=-1e8;
|
alpha=-1e8;
|
||||||
beta=-1e8;
|
beta=-1e8;
|
||||||
|
double alpha_mid=-1e8;
|
||||||
|
double beta_mid=-1e8;
|
||||||
|
|
||||||
// Check if s.p0 lies on p and vice versa, also check if tangents are in same direction (TODO: TEST)
|
// Check if s.p0 lies on p and vice versa, also check if tangents are in same direction (TODO: TEST)
|
||||||
// If so, assume overlapping splines
|
// If so, assume overlapping splines
|
||||||
// TODO: Better checks! False positives could happen here!
|
// TODO: Better checks! False positives could happen here!
|
||||||
IntersectSplineSegment1( p, s.StartPI(), p_mid, lam0, alpha );
|
IntersectSplineSegment1( p, s.StartPI(), p_mid, lam0, alpha, true );
|
||||||
IntersectSplineSegment1( s, p.StartPI(), s_mid, lam1, beta );
|
IntersectSplineSegment1( s, p.StartPI(), s_mid, lam1, beta, true );
|
||||||
|
|
||||||
|
// Also check if midpoints lie on other spline
|
||||||
|
IntersectSplineSegment1( p, s.GetPoint(0.5), p_mid, lam2, alpha_mid, true );
|
||||||
|
IntersectSplineSegment1( s, p.GetPoint(0.5), s_mid, lam3, beta_mid, true );
|
||||||
|
|
||||||
auto tang0 = s.GetTangent(0.);
|
auto tang0 = s.GetTangent(0.);
|
||||||
auto tang1 = p.GetTangent(alpha);
|
auto tang1 = p.GetTangent(alpha);
|
||||||
double err = tang0*tang1;
|
double err = tang0*tang1;
|
||||||
err*=err;
|
err*=err;
|
||||||
err *= 1.0/(tang0.Length2()*tang1.Length2());
|
err *= 1.0/(tang0.Length2()*tang1.Length2());
|
||||||
|
|
||||||
if(fabs(lam0) < 1e3*EPSILON && fabs(lam1) < 1e3*EPSILON /*&& err < EPSILON*/)
|
double constexpr eps = 1e3*EPSILON;
|
||||||
{
|
if(fabs(lam0)>eps) return false;
|
||||||
type = ClassifyOverlappingIntersection( alpha, beta );
|
if(fabs(lam1)>eps) return false;
|
||||||
return true;
|
if(fabs(lam2)>eps) return false;
|
||||||
}
|
if(fabs(lam3)>eps) return false;
|
||||||
return false;
|
if(fabs(1.0-err)>eps) return false;
|
||||||
|
|
||||||
|
type = ClassifyOverlappingIntersection( alpha, beta );
|
||||||
|
return true;
|
||||||
}
|
}
|
||||||
|
|
||||||
bool IsInsideTrig( const array<Point<2>,3> & t, Point<2> r )
|
bool IsInsideTrig( const array<Point<2>,3> & t, Point<2> r )
|
||||||
@ -1322,6 +1337,42 @@ Solid2d ClipSolids ( Solid2d s1, Solid2d s2, bool intersect)
|
|||||||
return res;
|
return res;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
bool Loop :: IsInside( Point<2> r ) const
|
||||||
|
{
|
||||||
|
int w = 0;
|
||||||
|
for(auto e : Edges(ALL))
|
||||||
|
{
|
||||||
|
int w_simple = CalcSide(*e.v0, *e.v1, r);
|
||||||
|
if(!e.v0->spline)
|
||||||
|
w += w_simple;
|
||||||
|
else
|
||||||
|
{
|
||||||
|
auto s = *e.v0->spline;
|
||||||
|
auto s0 = s.StartPI();
|
||||||
|
auto s1 = s.TangentPoint();
|
||||||
|
auto s2 = s.EndPI();
|
||||||
|
if(!IsInsideTrig( {s0, s1, s2} , r ))
|
||||||
|
w += w_simple;
|
||||||
|
else
|
||||||
|
{
|
||||||
|
// r close to spline, need exact test
|
||||||
|
// idea: compute weight, such that r lies on spline
|
||||||
|
// weight increases -> same side of spline as control point, simple test gives correct result
|
||||||
|
// weight decreases -> opposite side of spline as control point, adding control point to test polygon gives correct result
|
||||||
|
double old_weight = s.GetWeight();
|
||||||
|
ComputeWeight( s, r );
|
||||||
|
double new_weight = s.GetWeight();
|
||||||
|
|
||||||
|
if(new_weight >= old_weight)
|
||||||
|
w += w_simple;
|
||||||
|
else
|
||||||
|
w += CalcSide(s0, s1, r) + CalcSide(s1, s2, r);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return ( (w % 2) != 0 );
|
||||||
|
}
|
||||||
|
|
||||||
Solid2d :: Solid2d(const Array<std::variant<Point<2>, EdgeInfo>> & points, string name_, string bc)
|
Solid2d :: Solid2d(const Array<std::variant<Point<2>, EdgeInfo>> & points, string name_, string bc)
|
||||||
: name(name_)
|
: name(name_)
|
||||||
{
|
{
|
||||||
@ -1370,29 +1421,23 @@ Solid2d Solid2d :: operator-(const Solid2d & other_) const
|
|||||||
other.Append(RectanglePoly(-1e8, 1e8, -1e8, 1e8, "JUST_FOR_CLIPPING"));
|
other.Append(RectanglePoly(-1e8, 1e8, -1e8, 1e8, "JUST_FOR_CLIPPING"));
|
||||||
auto res = ClipSolids(*this, other);
|
auto res = ClipSolids(*this, other);
|
||||||
|
|
||||||
for (auto i : Range(other.polys))
|
|
||||||
{
|
|
||||||
auto & first = *other.polys[i].first;
|
|
||||||
if(first[0] == -1e8)
|
|
||||||
other.polys.DeleteElement(i);
|
|
||||||
}
|
|
||||||
res.name = name;
|
res.name = name;
|
||||||
return res;
|
return res;
|
||||||
}
|
}
|
||||||
|
|
||||||
Solid2d Solid2d :: operator+=(const Solid2d & other)
|
Solid2d & Solid2d :: operator+=(const Solid2d & other)
|
||||||
{
|
{
|
||||||
*this = *this + other;
|
*this = *this + other;
|
||||||
return *this;
|
return *this;
|
||||||
}
|
}
|
||||||
|
|
||||||
Solid2d Solid2d :: operator*=(const Solid2d & other)
|
Solid2d & Solid2d :: operator*=(const Solid2d & other)
|
||||||
{
|
{
|
||||||
*this = *this * other;
|
*this = *this * other;
|
||||||
return *this;
|
return *this;
|
||||||
}
|
}
|
||||||
|
|
||||||
Solid2d Solid2d :: operator-=(const Solid2d & other)
|
Solid2d & Solid2d :: operator-=(const Solid2d & other)
|
||||||
{
|
{
|
||||||
*this = *this - other;
|
*this = *this - other;
|
||||||
return *this;
|
return *this;
|
||||||
@ -1432,25 +1477,42 @@ bool Solid2d :: IsInside( Point<2> r ) const
|
|||||||
{
|
{
|
||||||
int w = 0;
|
int w = 0;
|
||||||
for(auto & poly : polys)
|
for(auto & poly : polys)
|
||||||
for(auto v : poly.Vertices(ALL))
|
w += poly.IsInside(r);
|
||||||
w += CalcSide(*v, *v->next, r);
|
|
||||||
return ( (w % 2) != 0 );
|
return ( (w % 2) != 0 );
|
||||||
}
|
}
|
||||||
|
|
||||||
bool Solid2d :: IsLeftInside( const Vertex & p0 )
|
bool Solid2d :: IsLeftInside( const Vertex & p0 )
|
||||||
{
|
{
|
||||||
auto & p1 = *p0.next;
|
auto & p1 = *p0.next;
|
||||||
|
if(p0.spline)
|
||||||
|
{
|
||||||
|
auto s = *p0.spline;
|
||||||
|
auto v = s.GetTangent(0.5);
|
||||||
|
auto n = Vec<2>{-v[1], v[0]};
|
||||||
|
auto q = s.GetPoint(0.5) + 1e-6*n;
|
||||||
|
return IsInside(q);
|
||||||
|
}
|
||||||
auto v = p1-p0;
|
auto v = p1-p0;
|
||||||
auto n = Vec<2>{v[1], -v[0]};
|
auto n = Vec<2>{-v[1], v[0]};
|
||||||
auto q = p0 + 0.5*v + 1e-6*n;
|
auto q = p0 + 0.5*v + 1e-6*n;
|
||||||
|
|
||||||
return IsInside(q);
|
return IsInside(q);
|
||||||
}
|
}
|
||||||
|
|
||||||
bool Solid2d :: IsRightInside( const Vertex & p0 )
|
bool Solid2d :: IsRightInside( const Vertex & p0 )
|
||||||
{
|
{
|
||||||
auto & p1 = *p0.next;
|
auto & p1 = *p0.next;
|
||||||
|
if(p0.spline)
|
||||||
|
{
|
||||||
|
auto s = *p0.spline;
|
||||||
|
auto v = s.GetTangent(0.5);
|
||||||
|
auto n = Vec<2>{v[1], -v[0]};
|
||||||
|
auto q = s.GetPoint(0.5) + 1e-6*n;
|
||||||
|
return IsInside(q);
|
||||||
|
}
|
||||||
|
|
||||||
auto v = p1-p0;
|
auto v = p1-p0;
|
||||||
auto n = Vec<2>{-v[1], v[0]};
|
auto n = Vec<2>{v[1], -v[0]};
|
||||||
auto q = p0 + 0.5*v + 1e-6*n;
|
auto q = p0 + 0.5*v + 1e-6*n;
|
||||||
return IsInside(q);
|
return IsInside(q);
|
||||||
}
|
}
|
||||||
@ -1584,7 +1646,7 @@ shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
|
|||||||
|
|
||||||
if(li!=ri)
|
if(li!=ri)
|
||||||
{
|
{
|
||||||
if(s.IsLeftInside(p0) == flip)
|
if(s.IsLeftInside(p0) != flip)
|
||||||
ls.left = dom;
|
ls.left = dom;
|
||||||
else
|
else
|
||||||
ls.right = dom;
|
ls.right = dom;
|
||||||
|
@ -447,13 +447,7 @@ struct Loop
|
|||||||
v->prev->pnext = std::move(v->pnext);
|
v->prev->pnext = std::move(v->pnext);
|
||||||
}
|
}
|
||||||
|
|
||||||
bool IsInside( Point<2> r ) const
|
bool IsInside( Point<2> r ) const;
|
||||||
{
|
|
||||||
int w = 0;
|
|
||||||
for(auto e : Edges(ALL))
|
|
||||||
w += CalcSide(*e.v0, *e.v1, r);
|
|
||||||
return ( (w % 2) != 0 );
|
|
||||||
}
|
|
||||||
|
|
||||||
EdgeIterator Edges(IteratorType iterType) const
|
EdgeIterator Edges(IteratorType iterType) const
|
||||||
{
|
{
|
||||||
@ -575,9 +569,9 @@ struct Solid2d
|
|||||||
Solid2d operator-(const Solid2d & other) const;
|
Solid2d operator-(const Solid2d & other) const;
|
||||||
|
|
||||||
Solid2d& operator=(const Solid2d & other) = default;
|
Solid2d& operator=(const Solid2d & other) = default;
|
||||||
Solid2d operator+=(const Solid2d & other);
|
Solid2d& operator+=(const Solid2d & other);
|
||||||
Solid2d operator*=(const Solid2d & other);
|
Solid2d& operator*=(const Solid2d & other);
|
||||||
Solid2d operator-=(const Solid2d & other);
|
Solid2d& operator-=(const Solid2d & other);
|
||||||
|
|
||||||
void Append( const Loop & poly )
|
void Append( const Loop & poly )
|
||||||
{
|
{
|
||||||
|
60
tests/pytest/test_csg2d.py
Normal file
60
tests/pytest/test_csg2d.py
Normal file
@ -0,0 +1,60 @@
|
|||||||
|
from netgen.geom2d import *
|
||||||
|
import pytest
|
||||||
|
import math
|
||||||
|
from pytest import approx
|
||||||
|
|
||||||
|
def test_two_circles():
|
||||||
|
c1 = Circle(center=(0,0), radius=1)
|
||||||
|
c2 = c1.Rotate(deg=45)
|
||||||
|
s = c1*c2
|
||||||
|
geo = CSG2d()
|
||||||
|
geo.Add(s)
|
||||||
|
m = geo.GenerateMesh()
|
||||||
|
assert len(m.Elements2D()) > 0
|
||||||
|
|
||||||
|
ngs = pytest.importorskip("ngsolve")
|
||||||
|
mesh = ngs.Mesh(m)
|
||||||
|
mesh.Curve(5)
|
||||||
|
assert ngs.Integrate(1.0, mesh) == approx(math.pi)
|
||||||
|
ngs.Draw(mesh)
|
||||||
|
|
||||||
|
def test_two_edge():
|
||||||
|
s = Solid2d( [(-1,0), cp(0,1), (1,0), cp(0,2)] )
|
||||||
|
geo = CSG2d()
|
||||||
|
geo.Add(s)
|
||||||
|
m = geo.GenerateMesh()
|
||||||
|
assert len(m.Elements2D()) > 0
|
||||||
|
|
||||||
|
ngs = pytest.importorskip("ngsolve")
|
||||||
|
g = geo.GenerateSplineGeometry()
|
||||||
|
ngs.Draw(g)
|
||||||
|
mesh = ngs.Mesh(m)
|
||||||
|
mesh.Curve(5)
|
||||||
|
ngs.Draw(mesh)
|
||||||
|
|
||||||
|
def test_trig_and_circle():
|
||||||
|
g = CSG2d()
|
||||||
|
|
||||||
|
trig = Solid2d( [(0,0), (1,1), (-1,1) ] ).BC("diamond")
|
||||||
|
circle = Circle( center=(0,0.101), radius=0.1).BC("circle") # TODO: Failing with center=(0,0.1)
|
||||||
|
|
||||||
|
d = trig-circle
|
||||||
|
g.Add(d)
|
||||||
|
g.Add(circle)
|
||||||
|
|
||||||
|
m = g.GenerateMesh(maxh=0.1)
|
||||||
|
assert len(m.Elements2D()) > 0
|
||||||
|
|
||||||
|
ngs = pytest.importorskip("ngsolve")
|
||||||
|
geo = g.GenerateSplineGeometry()
|
||||||
|
ngs.Draw(geo)
|
||||||
|
|
||||||
|
mesh = ngs.Mesh(m)
|
||||||
|
mesh.Curve(3)
|
||||||
|
ngs.Draw(mesh)
|
||||||
|
|
||||||
|
|
||||||
|
if __name__ == "__main__":
|
||||||
|
test_two_circles()
|
||||||
|
test_two_edge()
|
||||||
|
test_trig_and_circle()
|
Loading…
Reference in New Issue
Block a user