{find->get}_dofs in examples (#776)

docs/examples/ex02.py

@@ -69,7 +69,17 @@
 from skfem.models.poisson import unit_load
 import numpy as np
 
-m = MeshTri.init_symmetric().refined(3)
+m = (
+    MeshTri.init_symmetric()
+    .refined(3)
+    .with_boundaries(
+        {
+            "left": lambda x: x[0] == 0,
+            "right": lambda x: x[0] == 1,
+            "top": lambda x: x[1] == 1,
+        }
+    )
+)
 
 e = ElementTriMorley()
 ib = Basis(m, e)
@@ -91,18 +101,7 @@ def C(T):
 K = asm(bilinf, ib)
 f = 1e6 * asm(unit_load, ib)
 
-dofs = ib.find_dofs({
-    'left':  m.facets_satisfying(lambda x: x[0] == 0),
-    'right': m.facets_satisfying(lambda x: x[0] == 1),
-    'top':   m.facets_satisfying(lambda x: x[1] == 1),
-})
-
-D = np.concatenate((
-    dofs['left'].nodal['u'],
-    dofs['left'].facet['u_n'],
-    dofs['right'].nodal['u'],
-    dofs['top'].nodal['u'],
-))
+D = np.hstack([ib.get_dofs("left"), ib.get_dofs({"right", "top"}).all("u")])
 
 x = solve(*condense(K, f, D=D))
 

docs/examples/ex03.py

@@ -7,7 +7,11 @@
 
 m1 = MeshLine(np.linspace(0, 5, 50))
 m2 = MeshLine(np.linspace(0, 1, 10))
-m = m1 * m2
+m = (m1 * m2).with_boundaries(
+    {
+        "left": lambda x: x[0] == 0.0
+    }
+)
 
 e1 = ElementQuad1()
 
@@ -27,7 +31,7 @@ def mass(u, v, w):
 
 M = asm(mass, gb)
 
-D = gb.find_dofs({'left': m.facets_satisfying(lambda x: x[0] == 0.)})
+D = gb.get_dofs("left")
 y = gb.zeros()
 
 I = gb.complement_dofs(D)

docs/examples/ex04.py

@@ -65,9 +65,8 @@
 m = from_file(mesh_file)
 
 M = (
-    (MeshLine(np.linspace(0, 1, 6)) * MeshLine(np.linspace(-1, 1, 10)))
-    .translated((1.0, 0.0))
-    .refined()
+    (MeshLine(np.linspace(1, 2, 6)) * MeshLine(np.linspace(-1, 1, 10)))
+    .refined().with_boundaries({"contact": lambda x: x[0] == 1.0})
 )
 
 # define elements and bases
@@ -86,8 +85,8 @@
 ]
 
 mapping = MappingMortar.init_2D(m, M,
-                                m.boundaries['contact'],
-                                M.facets_satisfying(lambda x: x[0] == 1.0),
+                                m.boundaries["contact"],
+                                M.boundaries["contact"],
                                 np.array([0.0, 1.0]))
 
 mb = [

docs/examples/ex05.py

@@ -20,7 +20,7 @@
 from skfem.helpers import dot, grad
 from skfem.models.poisson import laplace
 
-m = MeshTri().refined(5)
+m = MeshTri().refined(5).with_boundaries({"plate": lambda x: x[1] == 0.0})
 
 e = ElementTriP1()
 
@@ -47,8 +47,7 @@ def facetlinf(v, w):
 
 b = asm(facetlinf, fb)
 
-D = ib.find_dofs({'plate': m.facets_satisfying(lambda x: (x[1] == 0.0))})
-I = ib.complement_dofs(D)
+I = ib.complement_dofs(ib.get_dofs("plate"))
 
 import scipy.sparse
 b = scipy.sparse.csr_matrix(b)

docs/examples/ex06.py

@@ -40,7 +40,7 @@
 
 f = asm(unit_load, ib)
 
-x = solve(*condense(K, f, D=ib.find_dofs()))
+x = solve(*condense(K, f, D=ib.get_dofs()))
 
 M, X = ib.refinterp(x, 3)
 

docs/examples/ex12.py

@@ -32,7 +32,7 @@
 A = asm(laplace, basis)
 b = asm(unit_load, basis)
 
-x = solve(*condense(A, b, D=basis.find_dofs()))
+x = solve(*condense(A, b, D=basis.get_dofs()))
 
 area = sum(b)
 k = b @ x / area**2

docs/examples/ex13.py

@@ -40,11 +40,9 @@
 basis = Basis(mesh, elements)
 A = asm(laplace, basis)
 
-boundary_dofs = basis.find_dofs()
-
 u = basis.zeros()
-u[boundary_dofs['positive'].all()] = 1.
-u = solve(*condense(A, x=u, D=boundary_dofs))
+u[basis.get_dofs("positive")] = 1.
+u = solve(*condense(A, x=u, D=basis.get_dofs({"positive", "ground"})))
 
 M = asm(mass, basis)
 u_exact = 2 * np.arctan2(*basis.doflocs[::-1]) / np.pi

docs/examples/ex14.py

@@ -39,7 +39,7 @@ def dirichlet(x):
 
 
 boundary_basis = FacetBasis(m, e)
-boundary_dofs = boundary_basis.find_dofs()['all'].all()
+boundary_dofs = basis.get_dofs()
 
 u = basis.zeros()
 u[boundary_dofs] = projection(dirichlet, boundary_basis, I=boundary_dofs)

docs/examples/ex15.py

@@ -12,7 +12,7 @@
 A = asm(laplace, basis)
 b = asm(unit_load, basis)
 
-x = solve(*condense(A, b, D=basis.find_dofs()))
+x = solve(*condense(A, b, D=basis.get_dofs()))
 
 if __name__ == "__main__":
     from skfem.visuals.matplotlib import plot, show

docs/examples/ex16.py

@@ -39,7 +39,8 @@ def stiffness(u, v, w):
 M = asm(mass, basis)
 
 ks, u = eigsh(L, M=M, sigma=0.)
-u /= u[basis.find_dofs()['all'].nodal['u'][-1], :]
+u /= u[basis.get_dofs().nodal["u"][-1], :]
+
 
 if __name__ == "__main__":
     fig, ax = subplots()

docs/examples/ex18.py

@@ -74,14 +74,14 @@ def body_force(v, w):
 f = np.concatenate([asm(body_force, basis['u']),
                     basis['p'].zeros()])
 
-uvp = solve(*condense(K, f, D=basis['u'].find_dofs()))
+uvp = solve(*condense(K, f, D=basis['u'].get_dofs()))
 
 velocity, pressure = np.split(uvp, [A.shape[0]])
 
 basis['psi'] = basis['u'].with_element(ElementTriP2())
 A = asm(laplace, basis['psi'])
 vorticity = asm(rot, basis['psi'], w=basis['u'].interpolate(velocity))
-psi = solve(*condense(A, vorticity, D=basis['psi'].find_dofs()))
+psi = solve(*condense(A, vorticity, D=basis['psi'].get_dofs()))
 
 
 if __name__ == '__main__':

docs/examples/ex19.py

@@ -64,7 +64,7 @@
 dt = .01
 print('dt =', dt)
 theta = 0.5                     # Crank–Nicolson
-L0, M0 = penalize(L, M, D=basis.find_dofs())
+L0, M0 = penalize(L, M, D=basis.get_dofs())
 A = M0 + theta * L0 * dt
 B = M0 - (1 - theta) * L0 * dt
 

docs/examples/ex20.py

@@ -49,7 +49,7 @@ def biharmonic(u, v, w):
 stokes = asm(biharmonic, ib)
 rotf = asm(unit_load, ib)
 
-psi = solve(*condense(stokes, rotf, D=ib.find_dofs()))
+psi = solve(*condense(stokes, rotf, D=ib.get_dofs()))
 (psi0,) = ib.probes(np.zeros((2, 1))) @ psi
 
 velocity = asm(

docs/examples/ex21.py

@@ -86,7 +86,7 @@ def mass(u, v, w):
 M = asm(mass, ib)
 
 L, x = solve(
-    *condense(K, M, D=ib.find_dofs()["fixed"]), solver=solver_eigen_scipy_sym()
+    *condense(K, M, D=ib.get_dofs("fixed")), solver=solver_eigen_scipy_sym()
 )
 
 if __name__ == "__main__":

docs/examples/ex23.py

@@ -61,7 +61,7 @@ def __init__(self, n: int):
                            ElementLineP1())
         self.lap = asm(laplace, self.basis)
         self.mass = asm(mass, self.basis)
-        self.D = self.basis.find_dofs()['all'].nodal['u']
+        self.D = self.basis.get_dofs()
 
     def inner(self, a: np.ndarray, b: np.ndarray) -> float:
         """return the inner product of two solutions"""

docs/examples/ex24.py

@@ -29,7 +29,7 @@
 basis = {variable: Basis(mesh, e, intorder=3)
          for variable, e in element.items()}
 
-D = np.concatenate([b.all() for b in basis['u'].find_dofs().values()])
+D = basis['u'].get_dofs(mesh.boundaries)
 
 A = asm(vector_laplace, basis['u'])
 B = -asm(divergence, basis['u'], basis['p'])
@@ -39,7 +39,7 @@
 uvp = np.zeros(K.shape[0])
 
 inlet_basis = FacetBasis(mesh, element['u'], facets=mesh.boundaries['inlet'])
-inlet_dofs = inlet_basis.find_dofs()['inlet'].all()
+inlet_dofs = inlet_basis.get_dofs('inlet')
 
 
 def parabolic(x):
@@ -56,7 +56,7 @@ def parabolic(x):
 A = asm(laplace, basis['psi'])
 psi = basis['psi'].zeros()
 vorticity = asm(rot, basis['psi'], w=basis['u'].interpolate(velocity))
-psi = solve(*condense(A, vorticity, D=basis['psi'].find_dofs()['floor'].all()))
+psi = solve(*condense(A, vorticity, D=basis['psi'].get_dofs('floor')))
 
 
 if __name__ == '__main__':

docs/examples/ex25.py

@@ -29,7 +29,8 @@
 
 mesh = MeshQuad.init_tensor(
     np.linspace(0, length, ceil(mesh_inlet_n / height * length)),
-    np.linspace(0, height / 2, mesh_inlet_n))
+    np.linspace(0, height / 2, mesh_inlet_n),
+).with_boundaries({"inlet": lambda x: x[0] == 0.0, "floor": lambda x: x[1] == 0.0})
 basis = Basis(mesh, ElementQuad2())
 
 
@@ -41,13 +42,11 @@ def advection(u, v, w):
     return v * velocity_0 * grad(u)[0]
 
 
-dofs = basis.find_dofs({'inlet': mesh.facets_satisfying(lambda x: x[0] == 0.),
-                        'floor': mesh.facets_satisfying(lambda x: x[1] == 0.)})
-interior = basis.complement_dofs(dofs)
+interior = basis.complement_dofs(basis.get_dofs({"inlet", "floor"}))
 
 A = asm(laplace, basis) + peclet * asm(advection, basis)
 t = basis.zeros()
-t[dofs['floor'].all()] = 1.
+t[basis.get_dofs("floor")] = 1.0
 t = solve(*condense(A, x=t, I=interior))
 
 basis0 = basis.with_element(ElementQuad0())

docs/examples/ex27.py

@@ -112,7 +112,7 @@ def __init__(self,
         self.basis['inlet'] = FacetBasis(self.mesh, self.element['u'],
                                          facets=self.mesh.boundaries['inlet'])
         self.basis['psi'] = self.basis['u'].with_element(ElementTriP2())
-        self.D = np.concatenate([b.all() for b in self.basis['u'].find_dofs().values()])
+        self.D = self.basis['u'].get_dofs([b for b in self.mesh.boundaries])
 
         A = asm(vector_laplace, self.basis['u'])
         B = asm(divergence, self.basis['u'], self.basis['p'])
@@ -121,7 +121,7 @@ def __init__(self,
         self.I = np.setdiff1d(np.arange(self.S.shape[0]), self.D)
 
     def inlet_dofs(self):
-        return self.basis['inlet'].find_dofs()['inlet'].all()
+        return self.basis["inlet"].get_dofs("inlet")
 
     @staticmethod
     def parabolic(x):
@@ -145,7 +145,7 @@ def streamfunction(self, velocity: np.ndarray) -> np.ndarray:
         psi = self.basis['psi'].zeros()
         vorticity = asm(rot, self.basis['psi'],
                         w=self.basis['u'].interpolate(velocity))
-        psi = solve(*condense(A, vorticity, D=self.basis['psi'].find_dofs()['floor'].all()))
+        psi = solve(*condense(A, vorticity, D=self.basis['psi'].get_dofs('floor')))
         return psi
 
     def mesh_plot(self):

docs/examples/ex28.py

@@ -107,10 +107,7 @@ def advection(u, v, w):
      * (asm(unit_load, basis['fluid-outlet'])
         + kratio * asm(unit_load, basis['solid-outlet'])))
 
-D = basis['heat'].find_dofs(
-    {label: boundary for
-     label, boundary in mesh.boundaries.items()
-     if label.endswith('-inlet')})
+D = basis["heat"].get_dofs([b for b in mesh.boundaries if b.endswith("-inlet")])
 I = basis['heat'].complement_dofs(D)
 
 
@@ -127,14 +124,17 @@ def exact(x: np.ndarray, y: np.ndarray) -> np.ndarray:
 
 temperature = solve(*condense(A, b, temperature, I=I))
 
-dofs = basis['heat'].find_dofs(
-    {label: facets for label, facets in mesh.boundaries.items()
-    if label.endswith('let')})
-
 exit_interface_temperature = {
-    'skfem': temperature[np.intersect1d(dofs['fluid-outlet'].all(),
-                                        dofs['solid-outlet'].all())[0]],
-    'exact': exact(length, -1.)
+    "skfem": temperature[
+        np.intersect1d(
+            *(
+                basis["heat"].get_dofs(b)
+                for b in mesh.boundaries
+                if b.endswith("-outlet")
+            )
+        )
+    ][0],
+    "exact": exact(length, -1.0),
 }
 
 if __name__ == '__main__':

docs/examples/ex29.py

@@ -89,7 +89,12 @@ def base_shear(u, v, w):
     return v * U.deriv()(w.x[0]) * u
 
 
-mesh = MeshLine(np.linspace(0, 1, 2**6))
+mesh = MeshLine(np.linspace(0, 1, 2**6)).with_boundaries(
+    {
+        "centre": lambda x: x[0] == 0,
+        "wall": lambda x: x[0] == 1
+    }
+)
 element = {'u': getattr(element_line, f'ElementLine{u_element}')(),
            'p': ElementLineP1()}
 basis = {v: Basis(mesh, e, intorder=4) for v, e in element.items()}
@@ -115,9 +120,8 @@ def base_shear(u, v, w):
 # the perturbation to the velocity vanishes on the centre-line z = 0,
 # z here being the sole coordinate.
 
-u_boundaries = basis['u'].find_dofs()['all'].all()
-walls = np.concatenate([u_boundaries,
-                        u_boundaries[1:] + basis['u'].N])
+walls = np.hstack([basis["u"].get_dofs(),
+                   basis["u"].get_dofs("wall").all() + basis['u'].N])
 
 pencil = condense(stiffness, mass_matrix, D=walls, expand=False)
 c = {'Criminale et al': np.loadtxt(Path(__file__).with_suffix('.csv'),

docs/examples/ex30.py

@@ -77,7 +77,7 @@ def body_force(v, w):
 A = asm(vector_laplace, basis['u'])
 B = -asm(divergence, basis['u'], basis['p'])
 f = asm(body_force, basis['u'])
-D = basis['u'].find_dofs()['all'].all()
+D = basis['u'].get_dofs()
 Aint = condense(A, D=D, expand=False)
 solver = solver_iter_pcg(M=build_pc_ilu(Aint))
 I = basis['u'].complement_dofs(D)
@@ -112,7 +112,7 @@ def dilatation(pressure: np.ndarray) -> np.ndarray:
 
 basis['psi'] = basis['u'].with_element(ElementQuad2())
 vorticity = asm(rot, basis['psi'], w=basis['u'].interpolate(velocity))
-psi = solve(*condense(asm(laplace, basis['psi']), vorticity, D=basis['psi'].find_dofs()))
+psi = solve(*condense(asm(laplace, basis['psi']), vorticity, D=basis['psi'].get_dofs()))
 psi0 = (basis['psi'].probes(np.zeros((2, 1))) @ psi)[0]
 
 

docs/examples/ex31.py

@@ -40,7 +40,7 @@
 A = asm(laplace, basis)
 M = asm(mass, basis)
 
-L, x = solve(*condense(A, M, D=basis.find_dofs()), solver=solver_eigen_scipy_sym(k=8))
+L, x = solve(*condense(A, M, D=basis.get_dofs()), solver=solver_eigen_scipy_sym(k=8))
 
 if __name__ == '__main__':
 

docs/examples/ex32.py

@@ -128,7 +128,7 @@ def body_force(v, w):
 f = np.concatenate([asm(body_force, basis['u']),
                     basis['p'].zeros()])
 
-D = basis['u'].find_dofs()
+D = basis['u'].get_dofs()
 Kint, fint, u, I = condense(K, f, D=D)
 Aint = Kint[:-(basis['p'].N), :-(basis['p'].N)]