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

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?).
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):

pip install -e .[full]
parse_log_pytest
requirements.txt
pytest
apt-get update apt-get install -y gcc gmsh libsuitesparse-dev
3.7
requirements/base.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
certifi @ file:///croot/certifi_1671487769961/work/certifi
cycler==0.11.0
exceptiongroup==1.2.2
fonttools==4.38.0
importlib-metadata==6.7.0
iniconfig==2.0.0
kiwisolver==1.4.5
markdown-it-py==2.2.0
matplotlib==3.5.3
mdurl==0.1.2
meshio==5.3.4
numpy==1.21.6
packaging==24.0
Pillow==9.5.0
pluggy==1.2.0
Pygments==2.17.2
pyparsing==3.1.4
pytest==7.4.4
python-dateutil==2.9.0.post0
rich==13.8.1
-e git+https://github.com/kinnala/scikit-fem.git@ff9abd7af7abb5cb0707c76507da2895482244f0#egg=scikit_fem
scipy==1.7.3
six==1.17.0
tomli==2.0.1
typing_extensions==4.7.1
zipp==3.15.0

name: scikit-fem
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2022.12.7=py37h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.4.4=h6a678d5_1
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - pip=22.3.1=py37h06a4308_0
  - python=3.7.16=h7a1cb2a_0
  - readline=8.2=h5eee18b_0
  - setuptools=65.6.3=py37h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.38.4=py37h06a4308_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - cycler==0.11.0
      - exceptiongroup==1.2.2
      - fonttools==4.38.0
      - importlib-metadata==6.7.0
      - iniconfig==2.0.0
      - kiwisolver==1.4.5
      - markdown-it-py==2.2.0
      - matplotlib==3.5.3
      - mdurl==0.1.2
      - meshio==5.3.4
      - numpy==1.21.6
      - packaging==24.0
      - pillow==9.5.0
      - pluggy==1.2.0
      - pygments==2.17.2
      - pyparsing==3.1.4
      - pytest==7.4.4
      - python-dateutil==2.9.0.post0
      - rich==13.8.1
      - scipy==1.7.3
      - six==1.17.0
      - tomli==2.0.1
      - typing-extensions==4.7.1
      - zipp==3.15.0
prefix: /opt/conda/envs/scikit-fem

tests/test_basis.py::TestCompositeSplitting::runTest tests/test_mesh.py::MeshTests::runTest tests/test_mesh.py::RefinePreserveSubsets::runTest
tests/test_mesh.py::FaultyInputs::runTest tests/test_mesh.py::Loading::runTest tests/test_mesh.py::SaveLoadCycle::runTest tests/test_mesh.py::SaveLoadCycleHex::runTest tests/test_mesh.py::SerializeUnserializeCycle::runTest
skfem/io/meshio.py skfem/mesh/mesh.py
skfem/io/meshio.py skfem/mesh/mesh.py