fix marking of vertical interfaces

LDDsimulation/boundary_and_interface.py

@@ -171,17 +171,21 @@ class BoundaryPart(df.SubDomain):
         ymin = min(p1[1], p2[1])
         ymax = max(p1[1], p2[1])
         # print(f"test if point {p} is on line segment between {p1} or {p2}")
-        # equality check
-        if np.absolute((p[0] - xmax)*(p[0] - xmin)) < tol and np.absolute((p[1] - ymax)*(p[1] - ymin)) < tol:
+        # check if p == p1 or p == p2
+        if np.fabs((p[0] - xmax)*(p[0] - xmin)) < tol and np.fabs((p[1] - ymax)*(p[1] - ymin)) < tol:
             #print(f"point {p} is close to either {p1} or {p2}")
             return True
 
         # check there holds p1[0] < p[0] < p2[0]. If not, p  cannot be on the line segment
-
-        # the test ((p[0] - xmax)*(p[0] - xmin) < 0) might fail if p[0] and one of p1[0]
-        # or p2[0] are equal. In pp1 or pp2 is vertical. We need to still calculate
-        # the distance to the line segment in this case
-        if ((p[0] - xmax)*(p[0] - xmin) < tol) or np.absolute((p[0] - xmax)*(p[0] - xmin)) < tol:
+        # same needs to be done for p1[1] < p[1] < p2[1]
+
+        # the test ((p[0] - xmax)*(p[0] - xmin) < 0) might fail if p[0] and one
+        # of p1[0] or p2[0] are equal. In this case the line pp1 or pp2 is
+        # vertical. We need to still calculate the distance to the line segment
+        # in this case. Same is true for ((p[1] - ymax)*(p[1] - ymin) < 0)
+        x_is_in_range = ((p[0] - xmax)*(p[0] - xmin) < tol) or np.fabs((p[0] - xmax)*(p[0] - xmin)) < tol
+        y_is_in_range = ((p[1] - ymax)*(p[1] - ymin) < tol) or np.fabs((p[1] - ymax)*(p[1] - ymin)) < tol
+        if x_is_in_range and y_is_in_range:
             # here we have a chance to actually lie on the line segment.
             segment_direction = (p2 - p1)/np.linalg.norm(p2 - p1)
             distance = (p - p1) - np.dot(segment_direction, p - p1)*segment_direction