From a8d5b09b87b162c25dec8e539e3597ff7fbb748f Mon Sep 17 00:00:00 2001
From: kinnala <kinnala@users.noreply.github.com>
Date: Thu, 13 Jun 2024 10:40:58 +0300
Subject: [PATCH] Fix expansion of facets to edges for interior facets (#1132)

---
 skfem/element/__init__.py |  3 ++-
 skfem/mesh/mesh.py        | 12 ++----------
 2 files changed, 4 insertions(+), 11 deletions(-)

diff --git a/skfem/element/__init__.py b/skfem/element/__init__.py
index 1ec3e674..282c102b 100644
--- a/skfem/element/__init__.py
+++ b/skfem/element/__init__.py
@@ -54,7 +54,8 @@
     ElementTetP2: ElementTriP2,
     ElementHex0: ElementQuad0,
     ElementHex1: ElementQuad1,
-    ElementHex2: ElementQuad2
+    ElementHex2: ElementQuad2,
+    ElementHex1DG: ElementQuad1,
 }
 
 # for backwards compatibility
diff --git a/skfem/mesh/mesh.py b/skfem/mesh/mesh.py
index 326c0068..f85a10b1 100644
--- a/skfem/mesh/mesh.py
+++ b/skfem/mesh/mesh.py
@@ -458,16 +458,8 @@ def _expand_facets(self, ix: ndarray) -> Tuple[ndarray, ndarray]:
         """
         vertices = np.unique(self.facets[:, ix].flatten())
 
-        if self.dim() == 3:
-            edge_candidates = self.t2e[:, self.f2t[0, ix]].flatten()
-            # subset of edges that share all points with the given facets
-            subset = np.nonzero(
-                np.prod(np.isin(self.edges[:, edge_candidates],
-                                self.facets[:, ix].flatten()),
-                        axis=0)
-            )[0]
-            edges = np.intersect1d(self.boundary_edges(),
-                                   edge_candidates[subset])
+        if self.dim() == 3 and self.bndelem is not None:
+            edges = np.unique(self.f2e[:, ix])
         else:
             edges = np.array([], dtype=np.int64)