docs: def get_edge_mask(self, subdomain=None):
          """Get faces which are fully in subdomain.
          """
          if subdomain is None:
              # https://stackoverflow.com/a/42392791/353337
              return numpy.s_[:]
  
          if subdomain not in self.subdomains:
              self._mark_vertices(subdomain)
  
          # A face is inside if all its edges are in.
          # An edge is inside if all its nodes are in.
          is_in = self.subdomains[subdomain]["vertices"][self.idx_hierarchy]
          # Take `all()` over the first index
          is_inside = numpy.all(is_in, axis=tuple(range(1)))
  
          if subdomain.is_boundary_only:
              # Filter for boundary
              is_inside = is_inside & self.is_boundary_edge
  
          return is_inside