From ac87e9b62cc77bf373b00e1e7d547b26b1cb7bc5 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 28 Aug 2020 14:21:45 +0200 Subject: [PATCH 1/6] csg2d - proper +=/-=/*= operator --- libsrc/geom2d/csg2d.cpp | 6 +++--- libsrc/geom2d/csg2d.hpp | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/libsrc/geom2d/csg2d.cpp b/libsrc/geom2d/csg2d.cpp index d1871166..af05475a 100644 --- a/libsrc/geom2d/csg2d.cpp +++ b/libsrc/geom2d/csg2d.cpp @@ -1380,19 +1380,19 @@ Solid2d Solid2d :: operator-(const Solid2d & other_) const return res; } -Solid2d Solid2d :: operator+=(const Solid2d & other) +Solid2d & Solid2d :: operator+=(const Solid2d & other) { *this = *this + other; return *this; } -Solid2d Solid2d :: operator*=(const Solid2d & other) +Solid2d & Solid2d :: operator*=(const Solid2d & other) { *this = *this * other; return *this; } -Solid2d Solid2d :: operator-=(const Solid2d & other) +Solid2d & Solid2d :: operator-=(const Solid2d & other) { *this = *this - other; return *this; diff --git a/libsrc/geom2d/csg2d.hpp b/libsrc/geom2d/csg2d.hpp index a7a92b12..42054bd5 100644 --- a/libsrc/geom2d/csg2d.hpp +++ b/libsrc/geom2d/csg2d.hpp @@ -575,9 +575,9 @@ struct Solid2d Solid2d operator-(const Solid2d & other) const; 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 ) { From f559cdef16c94ef88be0cee244bfd9ae0d76e853 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 28 Aug 2020 14:22:39 +0200 Subject: [PATCH 2/6] csg2d - better IsInside() check for splines --- libsrc/geom2d/csg2d.cpp | 63 +++++++++++++++++++++++++++++++++++------ libsrc/geom2d/csg2d.hpp | 8 +----- 2 files changed, 56 insertions(+), 15 deletions(-) diff --git a/libsrc/geom2d/csg2d.cpp b/libsrc/geom2d/csg2d.cpp index af05475a..31f707e5 100644 --- a/libsrc/geom2d/csg2d.cpp +++ b/libsrc/geom2d/csg2d.cpp @@ -1322,6 +1322,42 @@ Solid2d ClipSolids ( Solid2d s1, Solid2d s2, bool intersect) 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, EdgeInfo>> & points, string name_, string bc) : name(name_) { @@ -1370,12 +1406,6 @@ Solid2d Solid2d :: operator-(const Solid2d & other_) const other.Append(RectanglePoly(-1e8, 1e8, -1e8, 1e8, "JUST_FOR_CLIPPING")); 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; return res; } @@ -1432,23 +1462,40 @@ bool Solid2d :: IsInside( Point<2> r ) const { int w = 0; for(auto & poly : polys) - for(auto v : poly.Vertices(ALL)) - w += CalcSide(*v, *v->next, r); + w += poly.IsInside(r); return ( (w % 2) != 0 ); } bool Solid2d :: IsLeftInside( const Vertex & p0 ) { 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 n = Vec<2>{v[1], -v[0]}; auto q = p0 + 0.5*v + 1e-6*n; + return IsInside(q); } bool Solid2d :: IsRightInside( const Vertex & p0 ) { 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 n = Vec<2>{-v[1], v[0]}; auto q = p0 + 0.5*v + 1e-6*n; diff --git a/libsrc/geom2d/csg2d.hpp b/libsrc/geom2d/csg2d.hpp index 42054bd5..4fe00fe0 100644 --- a/libsrc/geom2d/csg2d.hpp +++ b/libsrc/geom2d/csg2d.hpp @@ -447,13 +447,7 @@ struct Loop v->prev->pnext = std::move(v->pnext); } - 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 ); - } + bool IsInside( Point<2> r ) const; EdgeIterator Edges(IteratorType iterType) const { From 1c825ebddf08d7e6ce16a47719a13d1b94751d37 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 28 Aug 2020 14:26:57 +0200 Subject: [PATCH 3/6] csg2d - better check for spline overlapping --- libsrc/geom2d/csg2d.cpp | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/libsrc/geom2d/csg2d.cpp b/libsrc/geom2d/csg2d.cpp index 31f707e5..d3e92c57 100644 --- a/libsrc/geom2d/csg2d.cpp +++ b/libsrc/geom2d/csg2d.cpp @@ -332,26 +332,38 @@ bool IsOverlapping( Spline p, Spline s, double & alpha, double & beta, Intersect double lam0 = -1e3*EPSILON; double lam1 = -1e3*EPSILON; + double lam2 = -1e3*EPSILON; + double lam3 = -1e3*EPSILON; alpha=-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) // If so, assume overlapping splines // TODO: Better checks! False positives could happen here! IntersectSplineSegment1( p, s.StartPI(), p_mid, lam0, alpha ); IntersectSplineSegment1( s, p.StartPI(), s_mid, lam1, beta ); + + // Also check if midpoints lie on other spline + IntersectSplineSegment1( p, s.GetPoint(0.5), p_mid, lam2, alpha_mid ); + IntersectSplineSegment1( s, p.GetPoint(0.5), s_mid, lam3, beta_mid ); + auto tang0 = s.GetTangent(0.); auto tang1 = p.GetTangent(alpha); double err = tang0*tang1; err*=err; err *= 1.0/(tang0.Length2()*tang1.Length2()); - if(fabs(lam0) < 1e3*EPSILON && fabs(lam1) < 1e3*EPSILON /*&& err < EPSILON*/) - { - type = ClassifyOverlappingIntersection( alpha, beta ); - return true; - } - return false; + double constexpr eps = 1e3*EPSILON; + if(fabs(lam0)>eps) return false; + if(fabs(lam1)>eps) return false; + if(fabs(lam2)>eps) return false; + if(fabs(lam3)>eps) return false; + if(fabs(1.0-err)>eps) return false; + + type = ClassifyOverlappingIntersection( alpha, beta ); + return true; } bool IsInsideTrig( const array,3> & t, Point<2> r ) From 956b06f90786fca6a53f453fdb30e06765360ae7 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 28 Aug 2020 17:26:43 +0200 Subject: [PATCH 4/6] csg2d - fix inside tests --- libsrc/geom2d/csg2d.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/libsrc/geom2d/csg2d.cpp b/libsrc/geom2d/csg2d.cpp index d3e92c57..5f004349 100644 --- a/libsrc/geom2d/csg2d.cpp +++ b/libsrc/geom2d/csg2d.cpp @@ -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 B = (p[1]-c[1])*(b[0]-p[0]) - (p[0]-c[0])*(b[1]-p[1]); double det = sqrt(-A*B); - double tt = (B-det)/(A+det); + double tt = fabs(A+det) fabs(v[1]) ? 0 : 1; double weight = fabs(tt*(p[dim]-a[dim])/v[dim] + 1.0/tt*(p[dim]-c[dim])/v[dim]); @@ -1485,12 +1485,12 @@ bool Solid2d :: IsLeftInside( const Vertex & p0 ) { auto s = *p0.spline; auto v = s.GetTangent(0.5); - auto n = Vec<2>{v[1], -v[0]}; + 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 n = Vec<2>{v[1], -v[0]}; + auto n = Vec<2>{-v[1], v[0]}; auto q = p0 + 0.5*v + 1e-6*n; return IsInside(q); @@ -1503,13 +1503,13 @@ bool Solid2d :: IsRightInside( const Vertex & p0 ) { auto s = *p0.spline; auto v = s.GetTangent(0.5); - auto n = Vec<2>{-v[1], v[0]}; + 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 n = Vec<2>{-v[1], v[0]}; + auto n = Vec<2>{v[1], -v[0]}; auto q = p0 + 0.5*v + 1e-6*n; return IsInside(q); } @@ -1643,7 +1643,7 @@ shared_ptr CSG2d :: GenerateSplineGeometry() if(li!=ri) { - if(s.IsLeftInside(p0) == flip) + if(s.IsLeftInside(p0) != flip) ls.left = dom; else ls.right = dom; From f2b9251032360f9f1fc2b8e5250b6f7ae82cbb55 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 28 Aug 2020 17:28:35 +0200 Subject: [PATCH 5/6] csg2d - tests --- tests/pytest/test_csg2d.py | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 tests/pytest/test_csg2d.py diff --git a/tests/pytest/test_csg2d.py b/tests/pytest/test_csg2d.py new file mode 100644 index 00000000..3afc29ef --- /dev/null +++ b/tests/pytest/test_csg2d.py @@ -0,0 +1,33 @@ +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 m.ne > 0 + + ngs = pytest.importorskip("ngsolve") + mesh = ngs.Mesh(m) + mesh.Curve(5) + assert ngs.Integrate(1.0, mesh) == approx(math.pi) + +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 m.ne > 0 + + ngs = pytest.importorskip("ngsolve") + g = geo.GenerateSplineGeometry() + ngs.Draw(g) + +if __name__ == "__main__": + test_two_circles() + test_two_edge() From 2a7d6bb55edc11353c967930175380389359092b Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 28 Aug 2020 18:35:35 +0200 Subject: [PATCH 6/6] csg2d - fix overlap detection, test --- libsrc/geom2d/csg2d.cpp | 21 ++++++++++++--------- tests/pytest/test_csg2d.py | 31 +++++++++++++++++++++++++++++-- 2 files changed, 41 insertions(+), 11 deletions(-) diff --git a/libsrc/geom2d/csg2d.cpp b/libsrc/geom2d/csg2d.cpp index 5f004349..c3e0c950 100644 --- a/libsrc/geom2d/csg2d.cpp +++ b/libsrc/geom2d/csg2d.cpp @@ -268,7 +268,7 @@ IntersectionType IntersectSplineSegment( const Spline & s, const Point<2> & r0, 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> p1 = s.TangentPoint(); @@ -310,11 +310,14 @@ IntersectionType IntersectSplineSegment1( const Spline & s, const Point<2> & r0, } int choice = 0; - if(vtype[0]==NO_INTERSECTION && vtype[1]!=NO_INTERSECTION) - choice = 1; + if(!first) + { + if(vtype[0]==NO_INTERSECTION && vtype[1]!=NO_INTERSECTION) + choice = 1; - if(valpha[0] < alpha+EPSILON) - choice = 1; + if(valpha[0] < alpha+EPSILON) + choice = 1; + } if(valpha[choice] < alpha+EPSILON) return NO_INTERSECTION; @@ -342,12 +345,12 @@ bool IsOverlapping( Spline p, Spline s, double & alpha, double & beta, Intersect // 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 // TODO: Better checks! False positives could happen here! - IntersectSplineSegment1( p, s.StartPI(), p_mid, lam0, alpha ); - IntersectSplineSegment1( s, p.StartPI(), s_mid, lam1, beta ); + IntersectSplineSegment1( p, s.StartPI(), p_mid, lam0, alpha, true ); + 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 ); - IntersectSplineSegment1( s, p.GetPoint(0.5), s_mid, lam3, beta_mid ); + 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 tang1 = p.GetTangent(alpha); diff --git a/tests/pytest/test_csg2d.py b/tests/pytest/test_csg2d.py index 3afc29ef..d962085c 100644 --- a/tests/pytest/test_csg2d.py +++ b/tests/pytest/test_csg2d.py @@ -10,24 +10,51 @@ def test_two_circles(): geo = CSG2d() geo.Add(s) m = geo.GenerateMesh() - assert m.ne > 0 + 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 m.ne > 0 + 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()