python-adaptive
diff --git a/‎adaptive/domain.py
Copy file name to clipboardExpand all lines: adaptive/domain.py
+30Lines changed: 30 additions & 0 deletions b/‎adaptive/domain.py
Copy file name to clipboardExpand all lines: adaptive/domain.py
+30Lines changed: 30 additions & 0 deletions
diff --git a/‎adaptive/tests/test_domain.py
Copy file name to clipboardExpand all lines: adaptive/tests/test_domain.py
+21-1Lines changed: 21 additions & 1 deletion b/‎adaptive/tests/test_domain.py
Copy file name to clipboardExpand all lines: adaptive/tests/test_domain.py
+21-1Lines changed: 21 additions & 1 deletion
@@ -87,6 +87,21 @@ def which_subdomains(self, x):
         ValueError : if x is outside of the domain
         """
 
+    @abc.abstractmethod
+    def encloses(self, points):
+        """Return whether the domain encloses the points
+
+        Parameters
+        ----------
+        points : a point or sequence of points
+
+        Returns
+        -------
+        Boolean (if a single point was provided) or an array of booleans
+        if a sequence of points was provided) that is True when
+        the domain encloses the point.
+        """
+
     @abc.abstractmethod
     def vertices(self):
         """Returns the vertices of the domain."""
@@ -276,6 +291,14 @@ def contains_subdomain(self, subdomain):
     def vertices(self):
         return self.points
 
+    def encloses(self, points):
+        a, b = self.bounds
+        points = np.asarray(points)
+        if points.shape == ():  # single point
+            return a <= points <= b
+        else:
+            return np.logical_and(a <= points, points <= b)
+
     def neighbors(self, subdomain, n=1):
         a, b = subdomain
         p = self.points
@@ -679,6 +702,13 @@ def neighbors(self, subdomain, n=1):
     def subdomains(self):
         return self.triangulation.simplices
 
+    def encloses(self, points):
+        points = np.asarray(points).T
+        A, b = self.bounds.equations[:, :-1], self.bounds.equations[:, -1:]
+        if len(points.shape) == 1:
+            points = points[:, None]
+        return np.all(A @ points + b <= 0, axis=0)
+
     def vertices(self):
         return self.triangulation.vertices
 
 
@@ -6,8 +6,10 @@
     a_few_points_inside,
     make_hypercube_domain,
     point_inside,
-    point_outside,
     point_on_shared_face,
+    point_outside,
+    points_inside,
+    points_outside,
 )
 from hypothesis import given, settings
 
@@ -66,6 +68,24 @@ def test_insert_points_outside_domain_raises(data, ndim):
         domain.insert(x)
 
 
+@pytest.mark.parametrize("ndim", [1, 2, 3])
+@given(data=st.data())
+def test_encloses(data, ndim):
+    domain = data.draw(make_hypercube_domain(ndim))
+
+    xin = data.draw(point_inside(domain))
+    assert domain.encloses(xin)
+
+    xout = data.draw(point_outside(domain))
+    assert not domain.encloses(xout)
+
+    xins = data.draw(points_inside(domain, 20))
+    assert np.all(domain.encloses(xins))
+
+    xouts = data.draw(points_outside(domain, 20))
+    assert not np.any(domain.encloses(xouts))
+
+
 @pytest.mark.parametrize("ndim", [1, 2, 3])
 @given(data=st.data())
 def test_split_at_point_outside_domain_raises(data, ndim):