From de760692837abc5191f3653b07a51bcc52bc32f8 Mon Sep 17 00:00:00 2001
From: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date: Wed, 11 Nov 2020 15:56:56 +0100
Subject: [PATCH 1/2] CSG2d bugfix

Fixes bug in IsInside(p) test for splines if p lies on an edge of the surrounding triangle.
Do a fast check using (new) function IsCloseToTrig() instead of IsInsideTrig().
---
 libsrc/geom2d/csg2d.cpp    |  8 +++++++-
 tests/pytest/test_csg2d.py | 29 +++++++++++++++++++++++++++++
 2 files changed, 36 insertions(+), 1 deletion(-)

diff --git a/libsrc/geom2d/csg2d.cpp b/libsrc/geom2d/csg2d.cpp
index fe0c8bd0..f97a0633 100644
--- a/libsrc/geom2d/csg2d.cpp
+++ b/libsrc/geom2d/csg2d.cpp
@@ -402,6 +402,12 @@ bool IsInsideTrig( const array<Point<2>,3> & t, Point<2> r )
   return ( (w % 2) != 0 );
 }
 
+bool IsCloseToTrig( const array<Point<2>,3> & t, Point<2> r, double eps=1e-4 )
+{
+  r += eps * (Center(t[0], t[1], t[2])-r); // move point a bit to center of trig
+  return IsInsideTrig( t, r );
+}
+
 
 IntersectionType IntersectTrig( Point<2> p0, Point<2> p1, const array<Point<2>,3> & trig)
 {
@@ -1612,7 +1618,7 @@ bool Loop :: IsInside( Point<2> r ) const
       auto s0 = s.StartPI();
       auto s1 = s.TangentPoint();
       auto s2 = s.EndPI();
-      if(!IsInsideTrig( {s0, s1, s2} , r ))
+      if(!IsCloseToTrig( {s0, s1, s2} , r ))
         w += w_simple;
       else
       {
diff --git a/tests/pytest/test_csg2d.py b/tests/pytest/test_csg2d.py
index 06e875a9..204e8ce1 100644
--- a/tests/pytest/test_csg2d.py
+++ b/tests/pytest/test_csg2d.py
@@ -54,6 +54,35 @@ def test_trig_and_circle():
     ngs.Draw(mesh)
 
 
+def test_circle_plus_rect():
+    circle = Circle( center=(0,0), radius=1 )
+    rect = Rectangle( pmin=(-0.5,0.0), pmax=(0.5,0.5) )
+
+    geo = CSG2d()
+    geo.Add(circle+rect)
+    m = geo.GenerateMesh(maxh=0.2)
+
+
+    ngs = pytest.importorskip("ngsolve")
+    mesh = ngs.Mesh(m)
+    mesh.Curve(5)
+    assert ngs.Integrate(1.0, mesh) == approx(math.pi)
+
+def test_circle_plus_rect1():
+    circle = Circle( center=(0,0), radius=1 )
+    rect = Rectangle( pmin=(-0.5,-0.5), pmax=(0.5,0.5) )
+
+    geo = CSG2d()
+    geo.Add(circle+rect)
+    m = geo.GenerateMesh(maxh=0.2)
+
+
+    ngs = pytest.importorskip("ngsolve")
+    mesh = ngs.Mesh(m)
+    mesh.Curve(5)
+    assert ngs.Integrate(1.0, mesh) == approx(math.pi)
+
+
 if __name__ == "__main__":
     test_two_circles()
     test_two_edge()

From 870b1479264afa9abd455dbd323ba075790be9ff Mon Sep 17 00:00:00 2001
From: Michael Neunteufel <michael.neunteufel@live.de>
Date: Wed, 11 Nov 2020 16:40:06 +0000
Subject: [PATCH 2/2] add ellipse as csg2d solid object

---
 python/geom2d.py | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/python/geom2d.py b/python/geom2d.py
index 886bfb60..3bd69d14 100644
--- a/python/geom2d.py
+++ b/python/geom2d.py
@@ -1,5 +1,6 @@
 from .libngpy._geom2d import SplineGeometry, Solid2d, CSG2d, Rectangle, Circle, EdgeInfo, PointInfo
 from .meshing import meshsize
+import math as math
 
 unit_square = SplineGeometry()
 _pnts = [ (0,0), (1,0), (1,1), (0,1) ]
@@ -144,3 +145,34 @@ def cp(p_or_px, py_or_none = None):
         return EdgeInfo(control_point=p)
     else:
         return EdgeInfo(control_point=(p_or_px,py_or_none))
+
+
+def Ellipse(center, a, b, bc="ellipse", mat="ellipse"):
+    """Creates ellipse centered at point center with principle axis a and b.
+
+    Parameters
+    ---------
+    center : Vec2
+      center of ellipse
+    a : Vec2
+      first principle axis, needs to be perpendicular to b
+    b : Vec2
+      second principle axis, needs to be perpendicular to a
+    bc : string
+      boundary name
+    mat : string
+      material name
+    """
+    if abs(a[0]*b[0] + a[1]*b[1]) > 1e-12:
+        raise Exception("In Ellipse: principle axis a and b are not perpendicular")
+    
+    ellipse = Circle( center=(0,0), radius=1.0, mat=mat, bc=bc )
+    
+    alpha = math.pi/2-math.atan2(a[0],a[1])
+    r_a = math.sqrt(a[0]**2+a[1]**2)
+    r_b = math.sqrt(b[0]**2+b[1]**2)
+    ellipse.Scale( (r_a,r_b) )
+    ellipse.Rotate( alpha/math.pi*180, center=(0,0) )
+    ellipse.Move( center )
+    
+    return ellipse