Merge branch 'fix_csg2d' into 'master'

Fix CSG2d bugs

See merge request jschoeberl/netgen!362

.gitlab-ci.yml

@@ -70,6 +70,7 @@ test_win:
   <<: *win
   stage: test
   script:
+    - pip install pytest-check
     - cd tests\pytest
     - cd %NETGEN_BUILD_DIR%\netgen
     - ctest -C Release -V --output-on-failure

libsrc/geom2d/csg2d.cpp

@@ -293,8 +293,11 @@ IntersectionType IntersectSplineSegment1( const Spline & s, const Point<2> & r0,
   double c_ = a0;
 
   double det =  b_*b_ - 4*a_*c_;
-  if(det<0.0)
-    return NO_INTERSECTION;
+  if(det<-EPSILON)
+      return NO_INTERSECTION;
+
+  if(det<EPSILON)
+      det = 0;
 
   double sqrt_det =  sqrt(det);
   double vbeta[2];
@@ -422,6 +425,64 @@ bool IsCloseToTrig( const array<Point<2>,3> & t, Point<2> r, double eps=1e-4 )
   return IsInsideTrig( t, r );
 }
 
+bool IsLeft( const Spline & s, Point<2> p )
+{
+  Point<2> a = s.StartPI();
+  Point<2> b = s.TangentPoint();
+  Point<2> c = s.EndPI();
+
+  // simple check by approximating spline with segment
+  bool is_left = Area(p, a, c) > 0.0;
+
+  // not close to spline -> simple check valid
+  if(!IsCloseToTrig( {a, b, c} , p ))
+      return is_left;
+
+  // p is control point -> simple check valid
+  auto bp = p-b;
+  if(bp.Length2() < EPSILON)
+      return is_left;
+
+  double sab = Area(p, a, b);
+  double sbc = Area(p, b, c);
+  if(fabs(sab)<EPSILON)
+      return is_left;
+  if(fabs(sbc)<EPSILON)
+      return is_left;
+
+  // 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();
+  auto s_tmp = s;
+  ComputeWeight( s_tmp, p );
+  double new_weight = s_tmp.GetWeight();
+
+  if(new_weight>old_weight)
+      return is_left;
+
+  double sabc = Area(a, b, c);
+
+  if (sabc > 0)
+  {
+    // chain makes a left turn
+    if (sab > 0 && sbc > 0)
+      return true;
+    else
+      return false;
+  }
+  else
+  {
+    // chain makes a right turn (or is straight)
+    if (sab < 0 && sbc < 0)
+      return false;
+    else
+      return true;
+  }
+}
+
+
 
 IntersectionType IntersectTrig( Point<2> p0, Point<2> p1, const array<Point<2>,3> & trig)
 {
@@ -743,8 +804,8 @@ void ComputeIntersections(Edge edgeP , Loop & l2)
 {
   for (Edge edgeQ : l2.Edges(SOURCE))
   {
-    double alpha = 0.0;
-    double beta = 0.0;
+    double alpha = -1;
+    double beta = -1;
     IntersectionType i = intersect(edgeP, edgeQ, alpha, beta);
     AddIntersectionPoint(edgeP, edgeQ, i, alpha, beta);
     if(i==X_INTERSECTION && (edgeP.v0->spline || edgeQ.v0->spline))
@@ -754,7 +815,6 @@ void ComputeIntersections(Edge edgeP , Loop & l2)
 
       // search for possible second intersection
       i = intersect(edgeP, edgeQ, alpha1, beta1);
-      // cout << "second intersection " << i << ',' << alpha1 << ',' << beta1 << ',' << alpha1-alpha << ',' << beta1-beta << endl;
       if(i!=NO_INTERSECTION && alpha+EPSILON<alpha1)
       {
         // Add midpoint of two intersection points to avoid false overlap detection of splines
@@ -819,47 +879,8 @@ enum RelativePositionType
   IS_P_p
 };
 
-RelativePositionType oracle(bool prev, Vertex* P1, Vertex* P2, Vertex* P3)
+inline RelativePositionType oracle_decide( double s1, double s2, double s3 )
 {
-  Vertex* Q;
-  Point<2> q;
-  if(prev)
-  {
-    Q = P2->neighbour->prev;
-    q = *Q;
-    if(Q->spline)
-      q = Q->spline->TangentPoint();
-  }
-  else
-  {
-    Q = P2->neighbour->next;
-    q = *Q;
-    if(P2->neighbour->spline)
-      q = P2->neighbour->spline->TangentPoint();
-  }
-
-  // is Q linked to P1 ?
-  if ( P1->is_intersection && (P1->neighbour == Q) )
-    return(IS_P_m);
-
-  // is Q linked to P2 ?
-  if ( P3->is_intersection && (P3->neighbour == Q) )
-    return(IS_P_p);
-
-  Point<2> p1 = *P1;
-  Point<2> p2 = *P2;
-  Point<2> p3 = *P3;
-
-  if(P1->spline)
-    p1 = P1->spline->TangentPoint();
-  if(P2->spline)
-    p3 = P2->spline->TangentPoint();
-
-  // check relative position of Q with respect to chain (P1,P2,P3)
-  double s1 = Area(  q, p1, p2);
-  double s2 = Area(  q, p2, p3);
-  double s3 = Area( p1, p2, p3);
-
   if (s3 > 0)
   {
     // chain makes a left turn
@@ -878,6 +899,139 @@ RelativePositionType oracle(bool prev, Vertex* P1, Vertex* P2, Vertex* P3)
   }
 }
 
+// no splines involved here
+// decides if Point q is left or right of chain (p1,p2,p3)
+RelativePositionType oracle_simple(Point<2> q, Point<2> p1, Point<2> p2, Point<2> p3)
+{
+  double s1 = Area(  q, p1, p2);
+  double s2 = Area(  q, p2, p3);
+  double s3 = Area( p1, p2, p3);
+
+  // check relative position of q with respect to chain (p1,p2,p3)
+  return oracle_decide(s1, s2, s3);
+}
+
+// (p1,p2) or (p2,p3) is a spline segment, compare with tangent (p1t,p2) instead of Segment (p1,p2)
+// BUT take care if tangent is collinear with (q,p2) (then use the segment (p1,p2) again)
+RelativePositionType oracle_spline_p(Point<2> q, Point<2> p1, Point<2> p1t, Point<2> p2, Point<2> p3, Point<2> p3t)
+{
+  double s1 = Area(  q, p1t, p2);
+  double s2 = Area(  q, p2, p3t);
+
+  if(fabs(s1) < EPSILON)
+    {
+      p1t = p1;
+      s1 = Area( q, p1t, p2 );
+    }
+
+  if(fabs(s2) < EPSILON)
+    {
+      p3t = p3;
+      s2 = Area( q, p2, p3t );
+    }
+
+  double s3 = Area( p1t, p2, p3t);
+
+  return oracle_decide(s1, s2, s3);
+}
+
+// (q,p2) is a spline segment, compare with tangent (qt,p2) instead of Segment (q,p2)
+// BUT take care if tangent at p2 is collinear with eiter (p1,p2) or (p2,p3) (then use the segment (q,p2) again)
+RelativePositionType oracle_spline_q(Point<2> q, Point<2> qt, Point<2> p1, Point<2> p2, Point<2> p3)
+{
+  double s1 = Area( qt, p1, p2);
+  double s2 = Area( qt, p2, p3);
+  double s3 = Area( p1, p2, p3);
+
+  if(fabs(s1) < EPSILON)
+      s1 = Area( q, p1, p2 );
+
+  if(fabs(s2) < EPSILON)
+      s2 = Area( q, p2, p3 );
+
+  return oracle_decide(s1, s2, s3);
+}
+
+// splines at (Q,P2) and either (P1,P2) or (P2,P3)
+// first use tangents to decide local orientation
+// if tangents of two splines match, use IsLeft(spline, other end point)
+// if tangent of spline and segment match, use simple methond (just end points)
+RelativePositionType oracle_spline(bool prev, Vertex *Q, Vertex *P1, Vertex *P2, Vertex *P3)
+{
+  Point<2> p1t =  *P1;
+  Point<2> p3t =  *P3;
+
+  auto sq = prev ? Q->spline : Q->prev->spline;
+  auto qt = sq->TangentPoint();
+  if(P1->spline) p1t = P1->spline->TangentPoint();
+  if(P2->spline) p3t = P2->spline->TangentPoint();
+
+  // Check using tangent directions first
+  double s1 = Area( qt, p1t, *P2 );
+  double s2 = Area( qt, *P2 , p3t);
+  double s3 = Area( p1t, *P2, p3t);
+
+  // tangents are facing in same direction
+  if(fabs(s1) < EPSILON)
+    {
+      if(P1->spline)
+          s1 = IsLeft(*P1->spline, *Q) ? 1 : -1;
+      else
+          s1 = Area( *Q, *P1, *P2 );
+    }
+
+  // tangents are facing in same direction
+  if(fabs(s2) < EPSILON)
+    {
+      if(P2->spline)
+          s2 = IsLeft(*P2->spline, *Q) ? 1 : -1;
+      else
+          s2 = Area( *Q, *P2, *P3 );
+    }
+
+  return oracle_decide(s1, s2, s3);
+}
+
+
+RelativePositionType oracle(bool prev, Vertex* P2)
+{
+  auto Q = prev ? P2->neighbour->prev : P2->neighbour->next;
+  auto sq = prev ? Q->spline : Q->prev->spline;
+  Vertex* P1 = P2->prev;
+  Vertex* P3 = P2->next;
+
+  // is Q linked to P1 ?
+  if ( P1->is_intersection && (P1->neighbour == Q) )
+    return(IS_P_m);
+
+  // is Q linked to P2 ?
+  if ( P3->is_intersection && (P3->neighbour == Q) )
+    return(IS_P_p);
+
+  // no splines -> simple variant
+  if(!P1->spline && !P2->spline && !Q->spline)
+      return oracle_simple(*Q, *P1, *P2, *P3);
+
+  Point<2> qt=*Q, p1t=*P1, p3t=*P3;
+
+  // splines -> also consider tangent points
+  if( sq) qt  =  Q->spline->TangentPoint();
+  if(P1->spline) p1t = P1->spline->TangentPoint();
+  if(P2->spline) p3t = P2->spline->TangentPoint();
+
+  // only spline at Q
+  if(!P1->spline && !P2->spline && Q->spline)
+      return oracle_spline_q(*Q, qt, *P1, *P2, *P3);
+
+  // only spline at P
+  if((P1->spline || !P2->spline) && !Q->spline)
+      return oracle_spline_p(*Q, *P1, p1t, *P2, *P3, p3t);
+
+  // spline at Q and P1 or P2
+  return oracle_spline(prev, Q, P1, P2, P3);
+}
+
+
 void LabelIntersections(Solid2d & sp, Solid2d & sq, Solid2d & sr, bool UNION)
 {
   auto & PP = sp.polys;
@@ -890,12 +1044,9 @@ void LabelIntersections(Solid2d & sp, Solid2d & sq, Solid2d & sr, bool UNION)
     {
 
       // determine local configuration at this intersection vertex
-      Vertex* P_m = I->prev;
-      Vertex* P_p = I->next;
-
       // check positions of Q- and Q+ relative to (P-, I, P+)
-      RelativePositionType Q_m_type = oracle(true,  P_m, I, P_p);
-      RelativePositionType Q_p_type = oracle(false, P_m, I, P_p);
+      RelativePositionType Q_m_type = oracle(true,  I);
+      RelativePositionType Q_p_type = oracle(false, I);
 
       // check non-overlapping cases
       if ((Q_m_type == LEFT  && Q_p_type == RIGHT) ||

tests/dockerfile

@@ -1,5 +1,6 @@
-FROM ubuntu:19.10
+FROM ubuntu:20.10
 ENV DEBIAN_FRONTEND=noninteractive
 MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>
-RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy python3-tk clang-tidy python3-distutils clang libocct-data-exchange-dev libcgns-dev libhdf5-dev
+RUN apt-get update && apt-get -y install python3 python3-pip libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy python3-tk clang-tidy python3-distutils clang libocct-data-exchange-dev libcgns-dev libhdf5-dev
+RUN python3 -m pip install pytest-check
 ADD . /root/src/netgen

tests/dockerfile_mpi

@@ -1,6 +1,6 @@
 FROM ubuntu:20.04
 ENV DEBIAN_FRONTEND=noninteractive
 MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>
-RUN apt-get update && apt-get -y install python3 libpython3-dev python3-pip libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy python3-tk python3-mpi4py clang-tidy python3-distutils clang libopenmpi-dev openmpi-bin gfortran
-RUN python3 -m pip install pytest-mpi
+RUN apt-get update && apt-get -y install python3 libpython3-dev python3-pip libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-numpy python3-tk python3-mpi4py clang-tidy python3-distutils clang libopenmpi-dev openmpi-bin gfortran
+RUN python3 -m pip install pytest-mpi pytest-check pytest
 ADD . /root/src/netgen

tests/pytest/test_csg2d.py

@@ -2,6 +2,20 @@ from netgen.geom2d import *
 import pytest
 import math
 from pytest import approx
+from pytest_check import check
+
+
+def check_area(geo, area):
+    if isinstance(geo, Solid2d):
+        g = CSG2d()
+        g.Add(geo)
+        geo = g
+
+    m = geo.GenerateMesh(maxh=0.2)
+    ngs = pytest.importorskip("ngsolve")
+    mesh = ngs.Mesh(m)
+    mesh.Curve(5)
+    with check: assert ngs.Integrate(1.0, mesh) == approx(area)
 
 def test_two_circles():
     c1 = Circle(center=(0,0), radius=1)
@@ -82,6 +96,18 @@ def test_circle_plus_rect1():
     mesh.Curve(5)
     assert ngs.Integrate(1.0, mesh) == approx(math.pi)
 
+def test_circle_and_rect():
+    c = Circle(center=(0,0),radius=1)
+    r = Rectangle((0,0),(1,1))
+
+    pi = math.pi
+    check_area(c-r, 3/4*pi)
+    check_area(c*r, 1/4*pi)
+    check_area(c+r, 3/4*pi+1)
+    check_area(r*c, 1/4*pi)
+    check_area(r+c, 3/4*pi+1)
+    check_area(r-c, 1-1/4*pi)
+
 
 if __name__ == "__main__":
     test_two_circles()