Modernize examples (#1133)

* Wrap forms in asm()

* make helpers.jump() do nothing if not assembled through asm()

* fix vedo to work with the latest release

* remove numpy 2.0 rc tests (numpy 2 released)

* remove default tags

* move to/from_dict to mesh.py

README.md

@@ -214,9 +214,14 @@ with respect to documented and/or tested features.
 - Fixed: `Mesh.p2e` returned incorrect incidence
 - Fixed: `InteriorFacetBasis.get_dofs` did not return all edge DOFs for 3D elements
 - Added: The lowest order, one point integration rule for tetrahedral elements
+- Added: `asm` will now wrap functions with three arguments using `BilinearForm`,
+  functions with two arguments using `LinearForm`, etc.
 - 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_defaults()`
 
 ### [9.1.1] - 2024-04-23
 

docs/advanced.rst

@@ -167,7 +167,6 @@ cube mesh:
      Number of elements: 1
      Number of vertices: 8
      Number of nodes: 8
-     Named boundaries [# facets]: left [1], bottom [1], front [1], right [1], top [1], back [1]
    >>> basis = Basis(m, ElementHex2())
    >>> basis
    <skfem CellBasis(MeshHex1, ElementHex2) object>

docs/examples/ex01.py

@@ -10,26 +10,25 @@
 m = MeshTri().refined(6)
 # or, with your own points and cells:
 # m = MeshTri(points, cells)
+# or, load from file
+# m = MeshTri.load("mesh.msh")
 
 e = ElementTriP1()
 basis = Basis(m, e)
 
-# this method could also be imported from skfem.models.laplace
+
 @BilinearForm
 def laplace(u, v, _):
     return dot(grad(u), grad(v))
 
 
-# this method could also be imported from skfem.models.unit_load
 @LinearForm
 def rhs(v, _):
     return 1.0 * v
 
-A = asm(laplace, basis)
-b = asm(rhs, basis)
-# or:
-# A = laplace.assemble(basis)
-# b = rhs.assemble(basis)
+
+A = laplace.assemble(basis)
+b = rhs.assemble(basis)
 
 # enforce Dirichlet boundary conditions
 A, b = enforce(A, b, D=m.boundary_nodes())

docs/examples/ex02.py

@@ -1,8 +1,6 @@
-r"""
+r"""Kirchhoff plate problem.
 
-This example demonstrates the solution of a slightly more complicated problem
-with multiple boundary conditions and a fourth-order differential operator. We
-consider the `Kirchhoff plate bending problem
+This example demonstrates the solution a fourth order `Kirchhoff plate bending problem
 <https://en.wikipedia.org/wiki/Kirchhoff%E2%80%93Love_plate_theory>`_ which
 finds its applications in solid mechanics. For a stationary plate of constant
 thickness :math:`d`, the governing equation reads: find the deflection :math:`u
@@ -16,20 +14,6 @@
 The Young's modulus of steel is :math:`E = 200 \cdot 10^9\,\text{Pa}` and Poisson
 ratio :math:`\nu = 0.3`.
 
-In reality, the operator
-
-.. math::
-    \frac{Ed^3}{12(1-\nu^2)} \Delta^2 
-is a combination of multiple first-order operators:
-
-.. math::
-    \boldsymbol{K}(u) = - \boldsymbol{\varepsilon}(\nabla u), \quad \boldsymbol{\varepsilon}(\boldsymbol{w}) = \frac12(\nabla \boldsymbol{w} + \nabla \boldsymbol{w}^T),
-.. math::
-    \boldsymbol{M}(u) = \frac{d^3}{12} \mathbb{C} \boldsymbol{K}(u), \quad \mathbb{C} \boldsymbol{T} = \frac{E}{1+\nu}\left( \boldsymbol{T} + \frac{\nu}{1-\nu}(\text{tr}\,\boldsymbol{T})\boldsymbol{I}\right),
-where :math:`\boldsymbol{I}` is the identity matrix. In particular,
-
-.. math::
-    \frac{Ed^3}{12(1-\nu^2)} \Delta^2 u = - \text{div}\,\textbf{div}\,\boldsymbol{M}(u).
 There are several boundary conditions that the problem can take.
 The *fully clamped* boundary condition reads
 
@@ -61,54 +45,51 @@
 <https://users.aalto.fi/~jakke74/WebFiles/Slides-Niiranen-ADMOS-09.pdf>`_ which
 is a piecewise quadratic :math:`C^0`-continuous element for biharmonic problems.
 
-The full source code of the example reads as follows:
-
-.. literalinclude:: examples/ex02.py
-    :start-after: EOF"""
+"""
 from skfem import *
 from skfem.models.poisson import unit_load
+from skfem.helpers import dd, ddot, trace, eye
 import numpy as np
 
-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,
-        }
-    )
-)
+m = (MeshTri
+     .init_symmetric()
+     .refined(3)
+     .with_defaults())
+basis = Basis(m, ElementTriMorley())
+
+d = 0.1
+E = 200e9
+nu = 0.3
+
 
-e = ElementTriMorley()
-ib = Basis(m, e)
+def C(T):
+    return E / (1 + nu) * (T + nu / (1 - nu) * eye(trace(T), 2))
 
 
 @BilinearForm
-def bilinf(u, v, w):
-    from skfem.helpers import dd, ddot, trace, eye
-    d = 0.1
-    E = 200e9
-    nu = 0.3
+def bilinf(u, v, _):
+    return d ** 3 / 12.0 * ddot(C(dd(u)), dd(v))
 
-    def C(T):
-        return E / (1 + nu) * (T + nu / (1 - nu) * eye(trace(T), 2))
 
-    return d**3 / 12.0 * ddot(C(dd(u)), dd(v))
+@LinearForm
+def load(v, _):
+    return 1e6 * v
 
 
-K = asm(bilinf, ib)
-f = 1e6 * asm(unit_load, ib)
+K = bilinf.assemble(basis)
+f = load.assemble(basis)
 
-D = np.hstack([ib.get_dofs("left"), ib.get_dofs({"right", "top"}).all("u")])
+D = np.hstack((
+    basis.get_dofs("left"),
+    basis.get_dofs({"right", "top"}).all("u"),
+))
 
 x = solve(*condense(K, f, D=D))
 
 def visualize():
     from skfem.visuals.matplotlib import draw, plot
     ax = draw(m)
-    return plot(ib,
+    return plot(basis,
                 x,
                 ax=ax,
                 shading='gouraud',

docs/examples/ex03.py

@@ -1,58 +1,56 @@
 """Linear elastic eigenvalue problem."""
 
 from skfem import *
-from skfem.helpers import dot, ddot, sym_grad
-from skfem.models.elasticity import linear_elasticity, linear_stress
+from skfem.helpers import dot, ddot, sym_grad, eye, trace
 import numpy as np
 
 m1 = MeshLine(np.linspace(0, 5, 50))
 m2 = MeshLine(np.linspace(0, 1, 10))
-m = (m1 * m2).with_boundaries(
-    {
-        "left": lambda x: x[0] == 0.0
-    }
-)
+m = (m1 * m2).with_defaults()
 
 e1 = ElementQuad1()
-
-mapping = MappingIsoparametric(m, e1)
-
 e = ElementVector(e1)
 
-gb = Basis(m, e, mapping, 2)
+basis = Basis(m, e, intorder=2)
 
 lam = 1.
 mu = 1.
-K = asm(linear_elasticity(lam, mu), gb)
+
+
+def C(T):
+    return 2. * mu * T + lam * eye(trace(T), T.shape[0])
+
+
+@BilinearForm
+def stiffness(u, v, w):
+    return ddot(C(sym_grad(u)), sym_grad(v))
+
 
 @BilinearForm
 def mass(u, v, w):
     return dot(u, v)
 
-M = asm(mass, gb)
 
-D = gb.get_dofs("left")
-y = gb.zeros()
+K = stiffness.assemble(basis)
+M = mass.assemble(basis)
 
-I = gb.complement_dofs(D)
+D = basis.get_dofs("left")
 
-L, x = solve(*condense(K, M, I=I),
+L, x = solve(*condense(K, M, D=D),
              solver=solver_eigen_scipy_sym(k=6, sigma=0.0))
 
-y = x[:, 4]
-
 # calculate stress
-sgb = gb.with_element(ElementVector(e))
-C = linear_stress(lam, mu)
-yi = gb.interpolate(y)
-sigma = sgb.project(C(sym_grad(yi)))
+y = x[:, 4]
+sbasis = basis.with_element(ElementVector(e))
+yi = basis.interpolate(y)
+sigma = sbasis.project(C(sym_grad(yi)))
 
 def visualize():
     from skfem.visuals.matplotlib import plot, draw
-    M = MeshQuad(np.array(m.p + .5 * y[gb.nodal_dofs]), m.t)
+    M = MeshQuad(np.array(m.p + .5 * y[basis.nodal_dofs]), m.t)
     ax = draw(M)
     return plot(M,
-                sigma[sgb.nodal_dofs[0]],
+                sigma[sbasis.nodal_dofs[0]],
                 ax=ax,
                 colorbar='$\sigma_{xx}$',
                 shading='gouraud')

docs/examples/ex05.py

@@ -15,17 +15,19 @@
  to a saddle point system.
 
 """
-
 from skfem import *
 from skfem.helpers import dot, grad
 from skfem.models.poisson import laplace
+import numpy as np
+import scipy.sparse
+
 
 m = MeshTri().refined(5).with_boundaries({"plate": lambda x: x[1] == 0.0})
 
 e = ElementTriP1()
 
-ib = Basis(m, e)
-fb = FacetBasis(m, e)
+basis = Basis(m, e)
+fbasis = FacetBasis(m, e)
 
 
 @BilinearForm
@@ -42,19 +44,18 @@ def facetlinf(v, w):
     return -dot(grad(v), n) * (x[0] == 1.0)
 
 
-A = asm(laplace, ib)
-B = asm(facetbilinf, fb)
+A = laplace.assemble(basis)
+B = facetbilinf.assemble(fbasis)
+b = facetlinf.assemble(fbasis)
 
-b = asm(facetlinf, fb)
+I = basis.complement_dofs(basis.get_dofs("plate"))
 
-I = ib.complement_dofs(ib.get_dofs("plate"))
 
-import scipy.sparse
 b = scipy.sparse.csr_matrix(b)
-K = scipy.sparse.bmat([[A+B, b.T], [b, None]], 'csr')
 
-import numpy as np
-f = np.concatenate((ib.zeros(), -1.0*np.ones(1)))
+# create a block system with an extra Lagrange multiplier
+K = scipy.sparse.bmat([[A + B, b.T], [b, None]], 'csr')
+f = np.concatenate((basis.zeros(), -1.0 * np.ones(1)))
 
 I = np.append(I, K.shape[0] - 1)
 

docs/examples/ex06.py

@@ -33,16 +33,15 @@
 
 e1 = ElementQuad1()
 e = ElementQuad2()
-mapping = MappingIsoparametric(m, e1)
-ib = Basis(m, e, mapping, 4)
+basis = Basis(m, e, intorder=4)
 
-K = asm(laplace, ib)
+K = laplace.assemble(basis)
 
-f = asm(unit_load, ib)
+f = unit_load.assemble(basis)
 
-x = solve(*condense(K, f, D=ib.get_dofs()))
+x = solve(*condense(K, f, D=basis.get_dofs()))
 
-M, X = ib.refinterp(x, 3)
+M, X = basis.refinterp(x, 3)
 
 if __name__ == "__main__":
     from os.path import splitext