remove default tags

README.md

@@ -219,6 +219,9 @@ with respect to documented and/or tested features.
 - Changed: Initializing `Basis` for `ElementTetP0` without specifying
   `intorder` or `quadrature` will now automatically fall back to a one
   point integration rule
+- Changed: Default tags ('left', 'right', 'top', ...) are no more
+  added automatically during mesh initialization, as a workaround you
+  can add them explicitly by calling `mesh = mesh.with_default_tags()`
 
 ### [9.1.1] - 2024-04-23
 

docs/examples/ex02.py

@@ -51,7 +51,10 @@
 from skfem.helpers import dd, ddot, trace, eye
 import numpy as np
 
-m = MeshTri.init_symmetric().refined(3)
+m = (MeshTri
+     .init_symmetric()
+     .refined(3)
+     .with_default_tags())
 basis = Basis(m, ElementTriMorley())
 
 d = 0.1

docs/examples/ex03.py

@@ -6,7 +6,7 @@
 
 m1 = MeshLine(np.linspace(0, 5, 50))
 m2 = MeshLine(np.linspace(0, 1, 10))
-m = m1 * m2
+m = (m1 * m2).with_default_tags()
 
 e1 = ElementQuad1()
 e = ElementVector(e1)

docs/examples/ex11.py

@@ -9,7 +9,7 @@
 from skfem.models.elasticity import lame_parameters
 
 
-m = MeshHex().refined(3)
+m = MeshHex().refined(3).with_default_tags()
 e = ElementVector(ElementHex1())
 basis = Basis(m, e, intorder=3)
 

skfem/mesh/mesh.py

@@ -216,6 +216,10 @@ def boundary_edges(self) -> ndarray:
         ))[0]
         return edge_candidates[ix]
 
+    def with_default_tags(self):
+        """Return a copy with the default tags ('left', 'right', ...)."""
+        return self.with_boundaries(self._build_default_tags())
+
     def with_boundaries(self,
                         boundaries: Dict[str, Union[Callable[[ndarray],
                                                              ndarray],
@@ -245,12 +249,14 @@ def with_boundaries(self,
         )
 
     def with_subdomains(self,
-                        subdomains: Dict[str, Callable[[ndarray], ndarray]]):
+                        subdomains: Dict[str, Union[Callable[[ndarray],
+                                                             ndarray],
+                                                    ndarray]]):
         """Return a copy of the mesh with named subdomains.
 
         Parameters
         ----------
-        boundaries
+        subdomains
             A dictionary of lambda functions with the names of the subdomains
             as keys.  The midpoint of the element should return ``True`` for
             the corresponding lambda function if the element belongs to the
@@ -267,6 +273,25 @@ def with_subdomains(self,
             },
         )
 
+    def _build_default_tags(self):
+
+        boundaries = {}
+        # default boundary names along the dimensions
+        minnames = ['left', 'bottom', 'front']
+        maxnames = ['right', 'top', 'back']
+        for d in range(self.doflocs.shape[0]):
+            dmin = np.min(self.doflocs[d])
+            ix = self.facets_satisfying(lambda x: x[d] == dmin)
+            if len(ix) >= 1:
+                boundaries[minnames[d]] = ix
+        for d in range(self.doflocs.shape[0]):
+            dmax = np.max(self.doflocs[d])
+            ix = self.facets_satisfying(lambda x: x[d] == dmax)
+            if len(ix) >= 1:
+                boundaries[maxnames[d]] = ix
+
+        return boundaries
+
     def _encode_point_data(self) -> Dict[str, List[ndarray]]:
 
         subdomains = {} if self._subdomains is None else self._subdomains
@@ -529,6 +554,7 @@ def __post_init__(self):
         M = self.elem.refdom.nnodes
 
         if self.nnodes > M and self.elem is not Element:
+            # this is for high-order meshes, input in a different format
             # reorder DOFs to the expected format: vertex DOFs are first
             # note: not run if elem is not set
             p, t = self.doflocs, self.t
@@ -557,28 +583,6 @@ def __post_init__(self):
                                "to C_CONTIGUOUS.")
             self.t = np.ascontiguousarray(self.t)
 
-        # try adding default boundary tags
-        if ((self._boundaries is None
-             and self.doflocs.shape[1] < 1e3
-             and self.elem is not Element)):
-            # note: not run if elem is not set
-            boundaries = {}
-            # default boundary names along the dimensions
-            minnames = ['left', 'bottom', 'front']
-            maxnames = ['right', 'top', 'back']
-            for d in range(self.doflocs.shape[0]):
-                dmin = np.min(self.doflocs[d])
-                ix = self.facets_satisfying(lambda x: x[d] == dmin)
-                if len(ix) >= 1:
-                    boundaries[minnames[d]] = ix
-            for d in range(self.doflocs.shape[0]):
-                dmax = np.max(self.doflocs[d])
-                ix = self.facets_satisfying(lambda x: x[d] == dmax)
-                if len(ix) >= 1:
-                    boundaries[maxnames[d]] = ix
-            if len(boundaries) > 0:
-                self._boundaries = boundaries
-
         # run validation
         if self.validate and logger.getEffectiveLevel() <= logging.DEBUG:
             self.is_valid()

tests/test_assembly.py

@@ -598,10 +598,14 @@ def proj(v, _):
 @pytest.mark.parametrize(
     "basis",
     [
-        Basis(MeshTri().refined(6), ElementTriRT1()),
-        Basis(MeshTri().refined(4), ElementTriRT2()),
-        Basis(MeshTri().refined(5), ElementTriBDM1()),
-        Basis(MeshQuad().refined(4), ElementQuadRT1()),
+        Basis(MeshTri().refined(6).with_default_tags(),
+              ElementTriRT1()),
+        Basis(MeshTri().refined(4).with_default_tags(),
+              ElementTriRT2()),
+        Basis(MeshTri().refined(5).with_default_tags(),
+              ElementTriBDM1()),
+        Basis(MeshQuad().refined(4).with_default_tags(),
+              ElementQuadRT1()),
     ]
 )
 def test_hdiv_boundary_integration(basis):
@@ -626,12 +630,10 @@ def test2(w):
 @pytest.mark.parametrize(
     "basis",
     [
-        Basis(MeshTet().refined(4).with_boundaries({
-            'right': lambda x: x[0] == 1
-        }), ElementTetRT1()),
-        Basis(MeshHex().refined(4).with_boundaries({
-            'right': lambda x: x[0] == 1
-        }), ElementHexRT1()),
+        Basis(MeshTet().refined(4).with_default_tags(),
+              ElementTetRT1()),
+        Basis(MeshHex().refined(4).with_default_tags(),
+              ElementHexRT1()),
     ]
 )
 def test_hdiv_boundary_integration_3d(basis):
@@ -655,9 +657,12 @@ def test2(w):
 @pytest.mark.parametrize(
     "basis",
     [
-        Basis(MeshTri().refined(4), ElementTriN1()),
-        Basis(MeshTri().refined(4), ElementTriN2()),
-        Basis(MeshQuad().refined(3), ElementQuadN1()),
+        Basis(MeshTri().refined(4).with_default_tags(),
+              ElementTriN1()),
+        Basis(MeshTri().refined(4).with_default_tags(),
+              ElementTriN2()),
+        Basis(MeshQuad().refined(3).with_default_tags(),
+              ElementQuadN1()),
     ]
 )
 def test_hcurl_boundary_integration(basis):
@@ -681,7 +686,8 @@ def test2(w):
 @pytest.mark.parametrize(
     "basis",
     [
-        Basis(MeshTet().refined(2), ElementTetN1()),
+        Basis(MeshTet().refined(2).with_default_tags(),
+              ElementTetN1()),
     ]
 )
 def test_hcurl_boundary_integration(basis):

tests/test_autodiff.py

@@ -68,7 +68,7 @@ def poisson(u, v, w):
 
 def test_navier_stokes():
 
-    m = MeshTri.init_sqsymmetric().refined(2)
+    m = MeshTri.init_sqsymmetric().refined(2).with_default_tags()
     basis = Basis(m, ElementVector(ElementTriP2()) * ElementTriP1())
     x = basis.zeros()
 
@@ -97,8 +97,11 @@ def navierstokes(u, p, v, q, w):
 
 def test_nonlin_elast():
 
-    m = MeshQuad.init_tensor(np.linspace(0, 5, 20),
-                         np.linspace(0, 0.5, 5)).to_meshtri(style='x')
+    m = (MeshQuad
+         .init_tensor(np.linspace(0, 5, 20),
+                      np.linspace(0, 0.5, 5))
+         .to_meshtri(style='x')
+         .with_default_tags())
     e = ElementVector(ElementTriP1())
     basis = Basis(m, e)
     x = basis.zeros()

tests/test_basis.py

@@ -92,7 +92,7 @@ class TestCompositeFacetAssembly(TestCase):
 
     def runTest(self):
 
-        m = MeshTri()
+        m = MeshTri().with_default_tags()
 
         fbasis1 = FacetBasis(m, ElementTriP1() * ElementTriP1(),
                              facets=m.facets_satisfying(lambda x: x[0] == 0))

tests/test_mesh.py

@@ -808,7 +808,7 @@ def test_incidence(mesh):
 
 def test_restrict_tags_boundary():
 
-    m = MeshTri().refined(3)
+    m = MeshTri().refined(3).with_default_tags()
     m = m.with_subdomains({
         'left': lambda x: x[0] <= 0.5,
         'bottom': lambda x: x[1] <= 0.5,

tests/test_utils.py

@@ -76,7 +76,7 @@ def test_simple_cg_solver():
 
 def test_mpc_periodic():
 
-    m = MeshQuad().refined(3)
+    m = MeshQuad().refined(3).with_default_tags()
     basis = Basis(m, ElementQuad1())
     A = laplace.assemble(basis)
     b = unit_load.assemble(basis)
@@ -93,7 +93,7 @@ def test_mpc_periodic():
 
 def test_mpc_2x_periodic():
 
-    m = MeshTri.init_sqsymmetric().refined(4)
+    m = MeshTri.init_sqsymmetric().refined(4).with_default_tags()
     e = ElementTriP2()
 
     basis = Basis(m, e)
@@ -118,7 +118,7 @@ def test_mpc_2x_periodic():
 
 def test_mpc_doubly_periodic():
 
-    m = MeshTri.init_sqsymmetric().refined(5)
+    m = MeshTri.init_sqsymmetric().refined(5).with_default_tags()
     e = ElementTriP1()
 
     basis = Basis(m, e)