Repository: kinnala/scikit-fem

Issue:
internal boundaries from meshio not in bndfacets
I was porting some old one-dimensional code to scikit-fem and found that whereas previously I had identified an internal vertex facet as a `Physical Point` in the Gmsh GEO, this wasn't picked up by `skfem.io.meshio.from_meshio`.

Roughly the code is like:
```python
from pygmsh import generate_mesh
from pygmsh.built_in import Geometry

import meshio
import skfem
from skfem.io import from_meshio

geom = Geometry()
points = []
lines = []
for x in range(3):
    points.append(geom.add_point([x, 0., 0.]))
    if len(points) > 1:
        lines.append(geom.add_line(*points[-2:]))                     
geom.add_physical(points[0], 'start')
geom.add_physical(points[1], 'middle')
geom.add_physical(points[2], 'end')
geom.add_physical(lines, 'domain')
mesh = from_meshio(generate_mesh(geom, dim=1))
``` 
Then `mesh.boundaries` is missing `'middle'`:
```python
{'start': array([0]), 'end': array([2])}
```

The reason is simple:

https://github.com/kinnala/scikit-fem/blob/61424c55508954e71302c78dc5e6637883f68a27/skfem/io/meshio.py#L101

`index` doesn't include `'middle'` because it's not in `bndfacets`.

I haven't decided whether the current behaviour is better or it should be changed to permit the definition of internal facets.  Here I was just going to use it as a probe point at which to simply sample the solution without having to interpolate.  I wonder whether the latter might be interesting for other applications (interfacial sources? domain decomposition into subdomains?).

Discussion:

Patch:
diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py
index 304ffe0c..8f75d157 100644
--- a/skfem/io/meshio.py
+++ b/skfem/io/meshio.py
@@ -76,7 +76,6 @@ def from_meshio(m, force_mesh_type=None):
                     return key
             return None
 
-        bndfacets = mtmp.boundary_facets()
         if m.cell_sets:         # MSH 4.1
             subdomains = {k: v[meshio_type]
                           for k, v in m.cell_sets_dict.items()
@@ -87,7 +86,7 @@ def from_meshio(m, force_mesh_type=None):
                       if bnd_type in v}
             boundaries = {k: np.array([i for i, f in
                                        enumerate(map(tuple, mtmp.facets.T))
-                                       if f in v and i in bndfacets])
+                                       if f in v])
                           for k, v in facets.items()}
         else: # MSH 2.2?
             elements_tag = m.cell_data_dict['gmsh:physical'][meshio_type]
@@ -110,7 +109,7 @@ def from_meshio(m, force_mesh_type=None):
             # get index of corresponding Mesh.facets for each meshio
             # facet found in the dict
             index = np.array([[dic[tuple(np.sort(mtmp.facets[:, i]))], i]
-                              for i in bndfacets
+                              for i in mtmp.boundary_facets()
                               if tuple(np.sort(mtmp.facets[:, i])) in dic])
 
             # read meshio tag numbers and names
diff --git a/skfem/mesh/mesh.py b/skfem/mesh/mesh.py
index b533f249..5e6deee5 100644
--- a/skfem/mesh/mesh.py
+++ b/skfem/mesh/mesh.py
@@ -402,7 +402,8 @@ class Mesh:
             'subdomains': subdomains,
         }
 
-    def define_boundary(self, name: str, test: Callable[[ndarray], ndarray]):
+    def define_boundary(self, name: str, test: Callable[[ndarray], ndarray],
+                        boundaries_only: bool = True):
         """Define a named boundary via function handle.
 
         Parameters
@@ -412,12 +413,13 @@ class Mesh:
         test
             A function which returns True for facet midpoints belonging to the
             boundary.
-
+        boundaries_only
+            If True, include only facets on the boundary of the mesh.
+            
         """
         if self.boundaries is None:
             self.boundaries = {}
-        self.boundaries[name] =\
-            self.facets_satisfying(test, boundaries_only=True)
+        self.boundaries[name] = self.facets_satisfying(test, boundaries_only)
 
     def copy(self):
         from copy import deepcopy


Test patch:
diff --git a/tests/test_basis.py b/tests/test_basis.py
index 6afaf95a..24bee71a 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -11,14 +11,24 @@ class TestCompositeSplitting(TestCase):
         """Solve Stokes problem, try splitting and other small things."""
 
         m = MeshTri()
-        m.refine(4)
+        m.refine()
+        m.define_boundary('centreline', lambda x: x[0] == .5,
+                          boundaries_only=False)
+        m.refine(3)
+
         e = ElementVectorH1(ElementTriP2()) * ElementTriP1()
 
         m.define_boundary('up', lambda x: x[1] == 1.)
         m.define_boundary('rest', lambda x: x[1] != 1.)
 
         basis = InteriorBasis(m, e)
+        self.assertEqual(
+            basis.get_dofs(m.boundaries['centreline']).all().size,
+            (2 + 1) * (2**(1 + 3) + 1) + 2 * 2**(1 + 3))
+        self.assertEqual(basis.boundary_dofs()['centreline'].all().size,
+                         (2 + 1) * (2**(1 + 3) + 1) + 2 * 2**(1 + 3))
 
+        
         @BilinearForm
         def bilinf(u, p, v, q, w):
             from skfem.helpers import grad, ddot, div
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 9f3753e1..030178d6 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -17,6 +17,14 @@ class MeshTests(unittest.TestCase):
         M = m.remove_elements(np.array([0]))
         self.assertEqual(M.t.shape[1], 7)
 
+        # Mesh.define_boundary
+        m.define_boundary('foo', lambda x: x[0] == 0.)
+        self.assertEqual(m.boundaries['foo'].size, 2)
+
+        # Mesh.define_boundary (internal)
+        m.define_boundary('bar', lambda x: x[0] == 1./2, boundaries_only=False)
+        self.assertEqual(m.boundaries['bar'].size, 2)
+
         # Mesh.scale, Mesh.translate
         m = MeshHex()
         m.scale(0.5)
@@ -75,17 +83,21 @@ class RefinePreserveSubsets(unittest.TestCase):
     """Check that uniform refinement preserves named boundaries."""
 
     def runTest(self):
-        for mtype in (MeshTri, MeshQuad):
+        for mtype in (MeshLine, MeshTri, MeshQuad):
             m = mtype()
             m.refine(2)
-            handle = lambda x: x[0] == 0.0
-            m.boundaries = {'test': m.facets_satisfying(handle)}
+            boundaries = [('external', lambda x: x[0] * (1. - x[0]) == 0.0),
+                          ('internal', lambda x: x[0] == 0.5)]
+            for name, handle in boundaries:
+                m.define_boundary(name, handle, boundaries_only=False)
             m.refine()
-            self.assertTrue((np.sort(m.boundaries['test'])
-                             == np.sort(m.facets_satisfying(handle))).all())
+            for name, handle in boundaries:
+                self.assertTrue((np.sort(m.boundaries[name])
+                                 == np.sort(m.facets_satisfying(handle))).all())
             m.refine(2)
-            self.assertTrue((np.sort(m.boundaries['test'])
-                             == np.sort(m.facets_satisfying(handle))).all())
+            for name, handle in boundaries:
+                self.assertTrue((np.sort(m.boundaries[name])
+                                 == np.sort(m.facets_satisfying(handle))).all())
 
 
 class SaveLoadCycle(unittest.TestCase):