diff --git a/.github/workflows/numpy2.yml b/.github/workflows/numpy2.yml deleted file mode 100644 index b226d2ad..00000000 --- a/.github/workflows/numpy2.yml +++ /dev/null @@ -1,42 +0,0 @@ -name: tests with numpy 2 - -on: - pull_request: - push: - branches: master - -jobs: - build: - - runs-on: ubuntu-latest - strategy: - matrix: - python-version: [3.9, '3.10', 3.11, 3.12] - - steps: - - uses: actions/checkout@v3 - - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v4 - with: - python-version: ${{ matrix.python-version }} - - name: Install dependencies - run: | - python -m pip install --upgrade pip - pip install -r requirements.txt - pip install --pre numpy==2.0.0rc1 - pip install . - - name: Run flake8 - run: | - flake8 skfem - - name: Run sphinx build - run: | - sphinx-build -W -a -b html docs docs/_build - - name: Run sphinx doctests - run: | - sphinx-build -W -a -b doctest docs docs/_build - - name: Run mypy - run: | - mypy skfem - - name: Run pytest - run: | - pytest -k 'not TestEx09 and not TestEx32 and not TestEx49' diff --git a/README.md b/README.md index f24a9308..549245db 100644 --- a/README.md +++ b/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 diff --git a/docs/advanced.rst b/docs/advanced.rst index 83e3b42d..8731a040 100644 --- a/docs/advanced.rst +++ b/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 diff --git a/docs/examples/ex01.py b/docs/examples/ex01.py index dd250bb9..53681b77 100644 --- a/docs/examples/ex01.py +++ b/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()) diff --git a/docs/examples/ex02.py b/docs/examples/ex02.py index 0c17b084..722e7f0e 100644 --- a/docs/examples/ex02.py +++ b/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 `_ 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 @@ `_ 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', diff --git a/docs/examples/ex03.py b/docs/examples/ex03.py index 855c731a..919f51b3 100644 --- a/docs/examples/ex03.py +++ b/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') diff --git a/docs/examples/ex05.py b/docs/examples/ex05.py index e25c0415..6bea0209 100644 --- a/docs/examples/ex05.py +++ b/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) diff --git a/docs/examples/ex06.py b/docs/examples/ex06.py index 15d0421a..3bd07967 100644 --- a/docs/examples/ex06.py +++ b/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 diff --git a/docs/examples/ex07.py b/docs/examples/ex07.py index 3e2443a3..58e1368e 100644 --- a/docs/examples/ex07.py +++ b/docs/examples/ex07.py @@ -12,23 +12,23 @@ bb = FacetBasis(m, e) fb = [InteriorFacetBasis(m, e, side=i) for i in [0, 1]] + @BilinearForm -def dgform(u, v, p): - ju, jv = jump(p, u, v) - h = p.h - n = p.n +def dgform(u, v, w): + ju, jv = jump(w, u, v) + h = w.h + n = w.n return ju * jv / (alpha * h) - dot(grad(u), n) * jv - dot(grad(v), n) * ju -@BilinearForm -def nitscheform(u, v, p): - h = p.h - n = p.n - return u * v / (alpha * h) - dot(grad(u), n) * v - dot(grad(v), n) * u -A = asm(laplace, ib) +A = laplace.assemble(ib) +C = dgform.assemble(bb) +b = unit_load.assemble(ib) + +# calling asm(form, [...], [...]) will automatically +# assemble all combinations from the lists and sum +# the result B = asm(dgform, fb, fb) -C = asm(nitscheform, bb) -b = asm(unit_load, ib) x = solve(A + B + C, b) diff --git a/docs/examples/ex09.py b/docs/examples/ex09.py index 9c6df02d..aed714bf 100644 --- a/docs/examples/ex09.py +++ b/docs/examples/ex09.py @@ -1,9 +1,5 @@ r"""Preconditioned conjugate gradient for 3-D Poisson. -.. note:: - - This example will make use of the external packages `PyAMG `_ or `pyamgcl `_, if installed. - Whereas most of the examples thus far have used direct linear solvers, this is not appropriate for larger problems, which includes most of those posed in three dimensions. @@ -35,35 +31,42 @@ * Demidov, D. (2019). AMGCL: an efficient, flexible, and extensible algebraic multigrid implementation. `arXiv:1811.05704 `_ """ - from skfem import * -from skfem.models.poisson import * +from skfem.helpers import * import numpy as np from scipy.sparse import spmatrix from scipy.sparse.linalg import LinearOperator + p = np.linspace(0, 1, 16) m = MeshTet.init_tensor(*(p,) * 3) +basis = Basis(m, ElementTetP1()) -e = ElementTetP1() -basis = Basis(m, e) -A = asm(laplace, basis) -b = asm(unit_load, basis) +@BilinearForm +def laplace(u, v, _): + return dot(grad(u), grad(v)) -I = m.interior_nodes() +@LinearForm +def unit_load(v, _): + return 1. * v + + +A = laplace.assemble(basis) +b = unit_load.assemble(basis) + +I = m.interior_nodes() x = 0. * b -if __name__ == "__main__": - verbose = True -else: - verbose = False Aint, bint = condense(A, b, I=I, expand=False) preconditioners = [None, build_pc_ilu(Aint, drop_tol=1e-3)] + +# note: pyamg does not support numpy 2.0 +# import algebraic multigrid from external package try: from pyamg import smoothed_aggregation_solver @@ -73,32 +76,20 @@ def build_pc_amgsa(A: spmatrix, **kwargs) -> LinearOperator: preconditioners.append(build_pc_amgsa(Aint)) -except ImportError: +except: print('Skipping PyAMG') -try: - import pyamgcl - - def build_pc_amgcl(A: spmatrix, **kwargs) -> LinearOperator: - """AMG preconditioner""" - - if hasattr(pyamgcl, 'amgcl'): # v. 1.3.99+ - return pyamgcl.amgcl(A, **kwargs) - else: - return pyamgcl.amg(A, **kwargs) - - preconditioners.append(build_pc_amgcl(Aint)) - -except ImportError: - print('Skipping pyamgcl') +# solve for each preconditioner for pc in preconditioners: - x[I] = solve(Aint, bint, solver=solver_iter_pcg(verbose=verbose, M=pc)) + x[I] = solve(Aint, bint, solver=solver_iter_pcg(verbose=True, M=pc)) -if verbose: +if __name__ == "__main__": from os.path import splitext from sys import argv - m.draw('vedo', point_data={'potential': x}).show() - #m.save(splitext(argv[0])[0] + ".vtk", {'potential': x}) + # use vedo: press 5 to visualize potential, X for cutter tool + basis.plot(x, 'vedo').show() + # preferred: save to vtk for visualization in Paraview + #m.save(splitext(argv[0])[0] + ".vtk", point_data{'potential': x}) diff --git a/docs/examples/ex10.py b/docs/examples/ex10.py index 5e5ff1e9..168071ec 100644 --- a/docs/examples/ex10.py +++ b/docs/examples/ex10.py @@ -3,7 +3,6 @@ This example solves the nonlinear minimal surface problem using Newton's method. """ - from skfem import * from skfem.helpers import grad, dot import numpy as np @@ -32,6 +31,8 @@ def rhs(v, p): D = m.boundary_nodes() x[D] = np.sin(np.pi * m.p[0, D]) + +# Newton iterations for itr in range(100): w = basis.interpolate(x) J = asm(jacobian, basis, prev=w) @@ -43,6 +44,7 @@ def rhs(v, p): if __name__ == "__main__": print(np.linalg.norm(x - x_prev)) + if __name__ == "__main__": from skfem.visuals.matplotlib import plot3, show plot3(m, x) diff --git a/docs/examples/ex11.py b/docs/examples/ex11.py index 8bf53c7d..eb07e930 100644 --- a/docs/examples/ex11.py +++ b/docs/examples/ex11.py @@ -1,39 +1,44 @@ r"""Linear elasticity. -This example solves the linear elasticity problem using trilinear elements. The -weak form of the linear elasticity problem is defined in -:func:`skfem.models.elasticity.linear_elasticity`. +This example solves the linear elasticity problem using trilinear elements. """ - import numpy as np from skfem import * -from skfem.models.elasticity import linear_elasticity, lame_parameters +from skfem.helpers import ddot, sym_grad, eye, trace +from skfem.models.elasticity import lame_parameters + + +m = MeshHex().refined(3).with_defaults() +e = ElementVector(ElementHex1()) +basis = Basis(m, e, intorder=3) + +# calculate Lamé parameters from Young's modulus and Poisson ratio +lam, mu = lame_parameters(1e3, 0.3) -m = MeshHex().refined(3) -e1 = ElementHex1() -e = ElementVector(e1) -ib = Basis(m, e, MappingIsoparametric(m, e1), 3) -K = asm(linear_elasticity(*lame_parameters(1e3, 0.3)), ib) +def C(T): + return 2. * mu * T + lam * eye(trace(T), T.shape[0]) -dofs = { - 'left' : ib.get_dofs(lambda x: x[0] == 0.0), - 'right': ib.get_dofs(lambda x: x[0] == 1.0), -} -u = ib.zeros() -u[dofs['right'].nodal['u^1']] = 0.3 +@BilinearForm +def stiffness(u, v, w): + return ddot(C(sym_grad(u)), sym_grad(v)) -I = ib.complement_dofs(dofs) -u = solve(*condense(K, x=u, I=I)) +K = stiffness.assemble(basis) + +u = basis.zeros() +u[basis.get_dofs('right').nodal['u^1']] = 0.3 + +u = solve(*condense(K, x=u, D=basis.get_dofs({'left', 'right'}))) sf = 1.0 -m = m.translated(sf * u[ib.nodal_dofs]) +m = m.translated(sf * u[basis.nodal_dofs]) if __name__ == "__main__": from os.path import splitext from sys import argv - + + # save to VTK for visualization in Paraview m.save(splitext(argv[0])[0] + '.vtk') diff --git a/docs/examples/ex32.py b/docs/examples/ex32.py index a3fa935a..d04f08f7 100644 --- a/docs/examples/ex32.py +++ b/docs/examples/ex32.py @@ -61,22 +61,14 @@ try: - try: - from pyamgcl import amgcl # v. 1.3.99+ - except ImportError: - from pyamgcl import amg as amgcl - from scipy.sparse.linalg import aslinearoperator - - def build_pc_amg(A: spmatrix, **kwargs) -> LinearOperator: - """AMG preconditioner""" - return aslinearoperator(amgcl(A, **kwargs)) - -except ImportError: from pyamg import smoothed_aggregation_solver - def build_pc_amg(A: spmatrix, **kwargs) -> LinearOperator: + def build_pc(A: spmatrix, **kwargs) -> LinearOperator: return smoothed_aggregation_solver(A, **kwargs).aspreconditioner() +except: + build_pc = lambda A: build_pc_ilu(A, drop_tol=1e-3) + class Sphere(NamedTuple): @@ -132,7 +124,7 @@ def body_force(v, w): Kint, fint, u, I = condense(K, f, D=D) Aint = Kint[:-(basis['p'].N), :-(basis['p'].N)] -Apc = build_pc_amg(Aint) +Apc = build_pc(Aint) diagQ = Q.diagonal() diff --git a/docs/gettingstarted.rst b/docs/gettingstarted.rst index 31485af7..009c81aa 100644 --- a/docs/gettingstarted.rst +++ b/docs/gettingstarted.rst @@ -90,13 +90,12 @@ unit square: Number of elements: 128 Number of vertices: 81 Number of nodes: 81 - Named boundaries [# facets]: left [8], bottom [8], right [8], top [8] .. plot:: from skfem import * - MeshTri().refined(3).draw(boundaries=True) + MeshTri().refined(3).draw() Step 4: Define a basis diff --git a/docs/howto.rst b/docs/howto.rst index 8c6fa778..3780887e 100644 --- a/docs/howto.rst +++ b/docs/howto.rst @@ -16,7 +16,7 @@ the boundary. Currently the main tool for finding DOFs is .. doctest:: >>> from skfem import MeshTri, Basis, ElementTriP2 - >>> m = MeshTri().refined(2) + >>> m = MeshTri().refined(2).with_defaults() >>> basis = Basis(m, ElementTriP2()) .. plot:: @@ -118,7 +118,7 @@ as follows: from skfem import * from skfem.visuals.matplotlib import * - m = MeshTri().refined(2) + m = MeshTri().refined(2).with_defaults() basis = Basis(m, ElementTriP2()) dofs = basis.get_dofs('left') ax = draw(m) diff --git a/skfem/assembly/__init__.py b/skfem/assembly/__init__.py index c76b3f89..368f0f9a 100644 --- a/skfem/assembly/__init__.py +++ b/skfem/assembly/__init__.py @@ -11,7 +11,6 @@ Number of elements: 2 Number of vertices: 4 Number of nodes: 4 - Named boundaries [# facets]: left [1], bottom [1], right [1], top [1] 2. Create :class:`~skfem.assembly.CellBasis` or :class:`~skfem.assembly.FacetBasis` objects. @@ -77,6 +76,17 @@ def asm(form: Form, supports assembling multiple bases at once and summing the result. """ + if not isinstance(form, Form) and callable(form): + # wrap the function + nargs = form.__code__.co_argcount + wrapper = [ + Functional, + LinearForm, + BilinearForm, + TrilinearForm, + ][nargs - 1] + form = wrapper(form) + assert form.form is not None logger.info("Assembling '{}'.".format(form.form.__name__)) nargs = [[arg] if not isinstance(arg, list) else arg for arg in args] diff --git a/skfem/helpers.py b/skfem/helpers.py index 2779bde7..422068f0 100644 --- a/skfem/helpers.py +++ b/skfem/helpers.py @@ -12,8 +12,7 @@ def jump(w: FormExtraParams, *args): if not hasattr(w, 'idx'): - raise NotImplementedError("jump() can be used only if the form is " - "assembled through asm().") + return args out = [] for i, arg in enumerate(args): out.append((-1.) ** w.idx[i] * arg) diff --git a/skfem/io/json.py b/skfem/io/json.py index a762bf0f..4acbdd10 100644 --- a/skfem/io/json.py +++ b/skfem/io/json.py @@ -1,47 +1,11 @@ -"""Import mesh from JSON.""" +"""Import mesh from JSON file.""" import json from os import PathLike from typing import Type - -import numpy as np - from skfem.mesh import MeshLine1, MeshTri, MeshQuad, MeshTet, MeshHex, Mesh -def to_dict(m): - boundaries = None - subdomains = None - if m.boundaries is not None: - boundaries = {k: v.tolist() for k, v in m.boundaries.items()} - if m.subdomains is not None: - subdomains = {k: v.tolist() for k, v in m.subdomains.items()} - return { - 'p': m.p.T.tolist(), - 't': m.t.T.tolist(), - 'boundaries': boundaries, - 'subdomains': subdomains, - } - - -def from_dict(cls, data): - if 'p' not in data or 't' not in data: - raise ValueError("Dictionary must contain keys 'p' and 't'.") - else: - data['p'] = np.ascontiguousarray(np.array(data['p']).T) - data['t'] = np.ascontiguousarray(np.array(data['t']).T) - if 'boundaries' in data and data['boundaries'] is not None: - data['boundaries'] = {k: np.array(v) - for k, v in data['boundaries'].items()} - if 'subdomains' in data and data['subdomains'] is not None: - data['subdomains'] = {k: np.array(v) - for k, v in data['subdomains'].items()} - data['doflocs'] = data.pop('p') - data['_subdomains'] = data.pop('subdomains') - data['_boundaries'] = data.pop('boundaries') - return cls(**data) - - def from_file(filename: PathLike) -> Mesh: with open(filename, 'r') as handle: d = json.load(handle) @@ -65,7 +29,7 @@ def from_file(filename: PathLike) -> Mesh: else: raise NotImplementedError("The given mesh is not supported.") - return from_dict(mesh_type, d) + return mesh_type.from_dict(d) def to_file(mesh: Mesh, filename: str): diff --git a/skfem/mesh/__init__.py b/skfem/mesh/__init__.py index 40b1b1f0..dd6ec989 100644 --- a/skfem/mesh/__init__.py +++ b/skfem/mesh/__init__.py @@ -16,7 +16,6 @@ Number of elements: 2 Number of vertices: 4 Number of nodes: 4 - Named boundaries [# facets]: left [1], bottom [1], right [1], top [1] Each mesh type has several constructors; see the docstring, e.g., ``help(MeshTri)`` or click :class:`~skfem.mesh.MeshTri` in the online diff --git a/skfem/mesh/mesh.py b/skfem/mesh/mesh.py index 2f2a5b22..17293974 100644 --- a/skfem/mesh/mesh.py +++ b/skfem/mesh/mesh.py @@ -216,6 +216,10 @@ def boundary_edges(self) -> ndarray: ))[0] return edge_candidates[ix] + def with_defaults(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() @@ -756,13 +760,38 @@ def load(cls, **kwargs) @classmethod - def from_dict(cls, d): - from skfem.io.json import from_dict - return from_dict(cls, d) + def from_dict(cls, data): + + if 'p' not in data or 't' not in data: + raise ValueError("Dictionary must contain keys 'p' and 't'.") + else: + data['p'] = np.ascontiguousarray(np.array(data['p']).T) + data['t'] = np.ascontiguousarray(np.array(data['t']).T) + if 'boundaries' in data and data['boundaries'] is not None: + data['boundaries'] = {k: np.array(v) + for k, v in data['boundaries'].items()} + if 'subdomains' in data and data['subdomains'] is not None: + data['subdomains'] = {k: np.array(v) + for k, v in data['subdomains'].items()} + data['doflocs'] = data.pop('p') + data['_subdomains'] = data.pop('subdomains') + data['_boundaries'] = data.pop('boundaries') + return cls(**data) def to_dict(self): - from skfem.io.json import to_dict - return to_dict(self) + + boundaries = None + subdomains = None + if self.boundaries is not None: + boundaries = {k: v.tolist() for k, v in self.boundaries.items()} + if self.subdomains is not None: + subdomains = {k: v.tolist() for k, v in self.subdomains.items()} + return { + 'p': self.p.T.tolist(), + 't': self.t.T.tolist(), + 'boundaries': boundaries, + 'subdomains': subdomains, + } @classmethod def from_mesh(cls, mesh, t: Optional[ndarray] = None): @@ -803,7 +832,15 @@ def init_refdom(cls): return cls(cls.elem.refdom.p, cls.elem.refdom.t, validate=False) def morphed(self, *args): - """Morph the mesh using functions.""" + """Morph the mesh using functions. + + Parameters + ---------- + funcs + One function per dimension, input is `p` and output is the + new coordinate of `p[i]', with `i=1,..,ndim`. + + """ p = self.p.copy() for i, arg in enumerate(args): if arg is None: diff --git a/skfem/visuals/vedo.py b/skfem/visuals/vedo.py index 75630fcd..b7586715 100644 --- a/skfem/visuals/vedo.py +++ b/skfem/visuals/vedo.py @@ -3,7 +3,7 @@ def draw(m, backend=False, **kwargs): """Visualize meshes.""" - from vedo import Plotter, UGrid + from vedo import Plotter, UnstructuredGrid vp = Plotter() grid = None with tempfile.NamedTemporaryFile() as tmp: @@ -11,9 +11,10 @@ def draw(m, backend=False, **kwargs): encode_cell_data=False, encode_point_data=True, **kwargs) - grid = UGrid(tmp.name + '.vtk') + grid = UnstructuredGrid(tmp.name + '.vtk') + vp += grid.tomesh() # save these for further use - grid.show = lambda: vp.show([grid.tomesh()]).close() + grid.show = lambda: vp.show() grid.plotter = vp return grid diff --git a/tests/test_assembly.py b/tests/test_assembly.py index befa3297..6faa4bcd 100644 --- a/tests/test_assembly.py +++ b/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_defaults(), + ElementTriRT1()), + Basis(MeshTri().refined(4).with_defaults(), + ElementTriRT2()), + Basis(MeshTri().refined(5).with_defaults(), + ElementTriBDM1()), + Basis(MeshQuad().refined(4).with_defaults(), + 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_defaults(), + ElementTetRT1()), + Basis(MeshHex().refined(4).with_defaults(), + 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_defaults(), + ElementTriN1()), + Basis(MeshTri().refined(4).with_defaults(), + ElementTriN2()), + Basis(MeshQuad().refined(3).with_defaults(), + 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_defaults(), + ElementTetN1()), ] ) def test_hcurl_boundary_integration(basis): diff --git a/tests/test_autodiff.py b/tests/test_autodiff.py index b8fb299f..f39af4b5 100644 --- a/tests/test_autodiff.py +++ b/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_defaults() 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_defaults()) e = ElementVector(ElementTriP1()) basis = Basis(m, e) x = basis.zeros() diff --git a/tests/test_basis.py b/tests/test_basis.py index 338e0b88..8ac092c0 100644 --- a/tests/test_basis.py +++ b/tests/test_basis.py @@ -92,7 +92,7 @@ class TestCompositeFacetAssembly(TestCase): def runTest(self): - m = MeshTri() + m = MeshTri().with_defaults() fbasis1 = FacetBasis(m, ElementTriP1() * ElementTriP1(), facets=m.facets_satisfying(lambda x: x[0] == 0)) diff --git a/tests/test_examples.py b/tests/test_examples.py index fae1453c..777a36f6 100644 --- a/tests/test_examples.py +++ b/tests/test_examples.py @@ -16,7 +16,7 @@ class TestEx02(TestCase): def runTest(self): import docs.examples.ex02 as ex02 - self.assertAlmostEqual(np.max(ex02.x[ex02.ib.nodal_dofs[0]]), + self.assertAlmostEqual(np.max(ex02.x[ex02.basis.nodal_dofs[0]]), 0.00033840961095522285) @@ -82,7 +82,7 @@ class TestEx11(TestCase): def runTest(self): import docs.examples.ex11 as ex11 u = ex11.u - ib = ex11.ib + ib = ex11.basis # since the mesh is symmetric, the mean values should equal to zero self.assertAlmostEqual(np.mean(u[ib.nodal_dofs[2, :]]), 0.0) self.assertAlmostEqual(np.mean(u[ib.nodal_dofs[1, :]]), 0.0) diff --git a/tests/test_mesh.py b/tests/test_mesh.py index de9f74ad..87a355ef 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -14,7 +14,6 @@ ElementHex0) from skfem.utils import projection from skfem.io.meshio import to_meshio, from_meshio -from skfem.io.json import to_dict, from_dict MESH_PATH = Path(__file__).parents[1] / 'docs' / 'examples' / 'meshes' @@ -114,7 +113,7 @@ def runTest(self): .refined(2) .with_boundaries({'down': lambda x: x[0] == 0,}) .with_subdomains({'up': lambda x: x[0] > 0.5})) - M = from_dict(cls, to_dict(m)) + M = cls.from_dict(m.to_dict()) self.assertTrue(np.sum(m.p - M.p) < 1e-13) self.assertTrue(np.sum(m.t - M.t) < 1e-13) for k in m.boundaries: @@ -808,7 +807,7 @@ def test_incidence(mesh): def test_restrict_tags_boundary(): - m = MeshTri().refined(3) + m = MeshTri().refined(3).with_defaults() m = m.with_subdomains({ 'left': lambda x: x[0] <= 0.5, 'bottom': lambda x: x[1] <= 0.5, diff --git a/tests/test_utils.py b/tests/test_utils.py index 741d7463..01fb1731 100644 --- a/tests/test_utils.py +++ b/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_defaults() 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_defaults() 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_defaults() e = ElementTriP1() basis = Basis(m, e)