diff --git a/docs/api.rst b/docs/api.rst index 3159b707..90167428 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -79,7 +79,7 @@ Subclasses of :class:`~skfem.assembly.basis.AbstractBasis` represent a global finite element basis evaluated at quadrature points. .. autoclass:: skfem.assembly.basis.AbstractBasis - :members: get_dofs + :members: get_dofs, interpolate, project Class: CellBasis **************** @@ -251,3 +251,8 @@ Module: skfem.helpers .. autofunction:: skfem.helpers.prod .. autofunction:: skfem.helpers.inv + +Module: skfem.visuals +===================== + +.. autofunction:: skfem.visuals.matplotlib.plot diff --git a/docs/conf.py b/docs/conf.py index 34509654..9bf49b1e 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -19,7 +19,7 @@ # -- Project information ----------------------------------------------------- project = 'scikit-fem' -copyright = '2018-2023, scikit-fem developers' +copyright = '2018-2024, scikit-fem developers' author = 'scikit-fem developers' @@ -92,7 +92,7 @@ # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # -html_theme = 'alabaster' +html_theme = 'classic' # Theme options are theme-specific and customize the look and feel of a theme # further. For a list of options available for each theme, see the diff --git a/docs/gettingstarted.rst b/docs/gettingstarted.rst index 009c81aa..4219d88b 100644 --- a/docs/gettingstarted.rst +++ b/docs/gettingstarted.rst @@ -13,8 +13,8 @@ install the package via Specifying ``[all]`` includes ``meshio`` for mesh input/output, and ``matplotlib`` for simple visualizations. The minimal dependencies are -``numpy`` and ``scipy``. You can also install scikit-fem in `Google Colab -`_ by executing +``numpy`` and ``scipy``. You can also install scikit-fem in Jupyter Notebook +or in `Google Colab `_ by executing .. code-block:: bash @@ -52,7 +52,7 @@ Next we write the forms .. math:: - a(u, v) = \int_\Omega \nabla u \cdot \nabla v \,\mathrm{d}x \quad \text{and} \quad L(v) = \int_\Omega f v \,\mathrm{d}x + a(u, v) = \int_\Omega \nabla u \cdot \nabla v \,\mathrm{d}x \quad \text{and} \quad l(v) = \int_\Omega f v \,\mathrm{d}x as source code. Each form is written as a function and decorated as follows: @@ -69,7 +69,7 @@ decorated as follows: >>> import numpy as np >>> @fem.LinearForm - ... def L(v, w): + ... def l(v, w): ... x, y = w.x # global coordinates ... f = np.sin(np.pi * x) * np.sin(np.pi * y) ... return f * v @@ -124,10 +124,10 @@ and the load vector has the type ``ndarray``. .. doctest:: >>> A = a.assemble(Vh) - >>> l = L.assemble(Vh) + >>> b = l.assemble(Vh) >>> A.shape (81, 81) - >>> l.shape + >>> b.shape (81,) Step 6: Find boundary DOFs @@ -144,18 +144,21 @@ the boundary. Empty call to Number of nodal DOFs: 32 ['u'] +:ref:`finddofs` explains how to match other subsets of DOFs. + Step 7: Eliminate boundary DOFs and solve ========================================= -The boundary DOFs must be eliminated from the linear system :math:`Ax=l` +The boundary DOFs must be eliminated from the linear system :math:`Ax=b` to set :math:`u=0` on the boundary. -This can be done using :func:`~skfem.utils.condense`. +This can be done using :func:`~skfem.utils.condense` +which can be useful also for inhomogeneous Dirichlet conditions. The output can be passed to :func:`~skfem.utils.solve` which is a simple wrapper to ``scipy`` sparse solver: .. doctest:: - >>> x = fem.solve(*fem.condense(A, l, D=D)) + >>> x = fem.solve(*fem.condense(A, b, D=D)) >>> x.shape (81,) @@ -167,12 +170,14 @@ which is a simple wrapper to ``scipy`` sparse solver: import numpy as np basis = Basis(MeshTri().refined(3), ElementTriP1()) a = BilinearForm(lambda u, v, _: dot(grad(u), grad(v))) - L = LinearForm(lambda v, w: np.sin(np.pi * w.x[0]) * np.sin(np.pi * w.x[1]) * v) - y = solve(*condense(a.assemble(basis), L.assemble(basis), D=basis.get_dofs())) + b = LinearForm(lambda v, w: np.sin(np.pi * w.x[0]) * np.sin(np.pi * w.x[1]) * v) + y = solve(*condense(a.assemble(basis), b.assemble(basis), D=basis.get_dofs())) ax = draw(basis) plot(basis, y, ax=ax, nrefs=2, colorbar=True, shading='gouraud') - +:ref:`visualizing` has some guidelines for visualization +and various other examples can be found in :ref:`gallery`. + Step 8: Calculate error ======================= @@ -182,7 +187,8 @@ The exact solution is known to be u(x, y) = \frac{1}{2 \pi^2} \sin \pi x \sin \pi y. -Thus, it makes sense to verify that the error is small. +Thus, it makes sense to verify that the error is small +by calculating the error in the :math:`L^2` norm: .. doctest:: @@ -194,3 +200,5 @@ Thus, it makes sense to verify that the error is small. ... return (uh - u) ** 2 >>> str(round(error.assemble(Vh, uh=Vh.interpolate(x)), 9)) '1.069e-06' + +:ref:`postprocessing` covers some ideas behind the use of :class:`~skfem.assembly.form.functional.Functional` wrapper. diff --git a/docs/howto.rst b/docs/howto.rst index 440561af..6d36d866 100644 --- a/docs/howto.rst +++ b/docs/howto.rst @@ -10,8 +10,17 @@ Finding degrees-of-freedom ========================== Often the goal is to constrain DOFs on a specific part of -the boundary. Currently the main tool for finding DOFs is +the boundary. The DOF numbers have one-to-one correspondence with +the rows and the columns of the finite element system matrix, +and the DOFs are typically passed to functions +such as :func:`~skfem.utils.condense` or :func:`~skfem.utils.enforce` +that modify the linear system to implement the boundary conditions. + +Currently the main tool for finding DOFs is :meth:`~skfem.assembly.basis.AbstractBasis.get_dofs`. +Let us look at the use of this method through an example +with a triangular mesh +and the quadratic Lagrange element: .. doctest:: @@ -19,19 +28,13 @@ the boundary. Currently the main tool for finding DOFs is >>> m = MeshTri().refined(2).with_defaults() >>> basis = Basis(m, ElementTriP2()) -.. plot:: - - from skfem import * - from skfem.visuals.matplotlib import * - m = MeshTri().refined(2) - basis = Basis(m, ElementTriP2()) - ax = draw(m) - for dof in basis.nodal_dofs.flatten(): - ax.text(*basis.doflocs[:, dof], str(dof)) - -We can provide an indicator function to +Normally, a list of facet indices is provided +for :meth:`~skfem.assembly.basis.AbstractBasis.get_dofs` +and it will find the matching DOF indices. +However, for convenience, +we can provide an indicator function to :meth:`~skfem.assembly.basis.AbstractBasis.get_dofs` and it will find the -DOFs on the matching facets: +DOFs on the matching facets (via their midpoints): .. doctest:: @@ -41,15 +44,30 @@ DOFs on the matching facets: >>> dofs.facet {'u': array([26, 30, 39, 40], dtype=int32)} -This element has one DOF per node and one DOF per facet. The facets have their -own numbering scheme starting from zero, however, the corresponding DOFs are -offset by the total number of nodal DOFs: +Here is a visualization of the nodal DOFs: + +.. plot:: + + from skfem import * + from skfem.visuals.matplotlib import * + m = MeshTri().refined(2) + basis = Basis(m, ElementTriP2()) + ax = draw(m) + for dof in basis.nodal_dofs.flatten(): + ax.text(*basis.doflocs[:, dof], str(dof)) + +This quadratic Lagrange element has one DOF per node and one DOF per +facet. Globally, the facets have their own numbering scheme starting +from zero, while the corresponding DOFs are offset by the total +number of nodal DOFs: .. doctest:: >>> dofs.facet['u'] array([26, 30, 39, 40], dtype=int32) +Here is a visualization of the facet DOFs: + .. plot:: from skfem import * @@ -84,10 +102,11 @@ the following table: +-----------+---------------------------------------------------------------+ | ``u^1^1`` | First component of the first component in a composite field | +-----------+---------------------------------------------------------------+ -| ``NA`` | Description not available (e.g., hierarchical or bubble DOF's)| +| ``NA`` | Description not available (e.g., hierarchical or bubble DOFs) | +-----------+---------------------------------------------------------------+ -An array of all DOFs with the key ``u`` can be obtained as follows: +Most of the time we just want an array of all DOFs with the key ``u``. +This can be obtained as follows: .. doctest:: @@ -97,7 +116,7 @@ An array of all DOFs with the key ``u`` can be obtained as follows: array([ 0, 2, 5, 10, 14, 26, 30, 39, 40], dtype=int32) If a set of facets is tagged, the name of the tag can be passed -to :meth:`~skfem.assembly.basis.AbstractBasis.get_dofs`: +also to :meth:`~skfem.assembly.basis.AbstractBasis.get_dofs`: .. doctest:: @@ -133,19 +152,25 @@ See :ref:`dofindexing` for more details. Performing projections ====================== +A common issue in finite element analysis is that you have either +an analytical function with a given expression or a finite element +function defined on one basis, while what you would like to have +instead is the same function defined on another finite element basis. + We can use :math:`L^2` projection to find discrete counterparts of functions or -transform from one finite element basis to another. Suppose we have +transform from one finite element basis to another. For example, +suppose we have :math:`u_0(x,y) = x^3 y^3` defined on the boundary of the domain and want to find the corresponding discrete function which is extended by zero in the -interior of the domain. You could explicitly assemble and solve the linear +interior of the domain. Technically, you could explicitly assemble and solve the linear system corresponding to: find :math:`\widetilde{u_0} \in V_h` satisfying .. math:: \int_{\partial \Omega} \widetilde{u_0} v\,\mathrm{d}s = \int_{\partial \Omega} u_0 v\,\mathrm{d}s\quad \forall v \in V_h. -However, this is so common that we have a shortcut -:meth:`~skfem.assembly.AbstractBasis.project`: +However, this is so common that we have a shortcut command +:meth:`~skfem.assembly.basis.AbstractBasis.project`: .. doctest:: @@ -174,7 +199,7 @@ However, this is so common that we have a shortcut ax = draw(ibasis) plot(ibasis, u0t, nrefs=3, ax=ax, colorbar=True, shading='gouraud') -We can also project over the entire domain: +As another example, we can also project over the entire domain: .. doctest:: @@ -198,7 +223,8 @@ We can also project over the entire domain: ax = draw(basis) plot(basis, fh, nrefs=3, ax=ax, colorbar=True, shading='gouraud') -We can project from one finite element basis to another: +Or alternatively, we can use the same command to +project from one finite element basis to another: .. doctest:: @@ -222,7 +248,7 @@ We can project from one finite element basis to another: ax = draw(basis) plot(basis0, fh, nrefs=3, ax=ax, colorbar=True, shading='gouraud') -We can interpolate the gradient at quadrature points and project: +We can also interpolate the gradient at quadrature points and then project: .. doctest:: @@ -250,11 +276,16 @@ We can interpolate the gradient at quadrature points and project: .. _predefined: -Discrete functions in forms -=========================== +Discrete functions in the forms +=============================== + +It is a common pattern to reuse an existing finite element function in the forms. +Everything within the form is expressed at the quadrature points and +the finite element functions must be interpolated +from nodes to the +quadrature points through :meth:`~skfem.assembly.basis.AbstractBasis.interpolate`. -We can use finite element functions inside the form by interpolating them at -quadrature points. For example, consider a fixed-point iteration for the +For example, consider a fixed-point iteration for the nonlinear problem .. math:: @@ -264,14 +295,14 @@ nonlinear problem u &= 0 \quad \text{on $\partial \Omega$}. \end{aligned} -We repeatedly find :math:`u_{k+1} \in H^1_0(\Omega)` which satisfies +We can repeatedly find :math:`u_{k+1} \in H^1_0(\Omega)` which satisfies .. math:: - \int_\Omega (u_{k} + \tfrac{1}{10}) \nabla u_{k+1} \cdot \nabla v \,\mathrm{d}x = \int_\Omega v\,\mathrm{d}x + \int_\Omega (u_{k} + \tfrac{1}{10}) \nabla u_{k+1} \cdot \nabla v \,\mathrm{d}x = \int_\Omega v\,\mathrm{d}x \quad \forall v \in H^1_0(\Omega). -for every :math:`v \in H^1_0(\Omega)`. -The bilinear form depends on the previous solution :math:`u_k`. +The bilinear form depends on the previous solution :math:`u_k` +which can be defined as follows: .. doctest:: @@ -342,38 +373,152 @@ The previous solution :math:`u_k` is interpolated at quadrature points using also as ``w.x``) corresponds to the global coordinates and ``w['h']`` (accessible also as ``w.h``) corresponds to the local mesh parameter. -Assembling jump terms -===================== +.. _postprocessing: -The shorthand :func:`~skfem.assembly.asm` -supports special syntax for assembling the same form over lists of -bases and summing the result. The form +Postprocessing the solution +=========================== -.. math:: +After solving the finite element system :math:`Ax=b`, it is common to +calculate derived quantities. +The most common techniques are: + +* Calculating gradient fields using the technique described at the end of :ref:`l2proj`. +* Using :class:`~skfem.assembly.form.functional.Functional` wrapper to calculate integrals + over the finite element solution. + +The latter consists of writing the integrand as a function +and decorating it using :class:`~skfem.assembly.form.functional.Functional`. +This is similar to the use of :class:`~skfem.assembly.form.bilinear_form.BilinearForm` +and :class:`~skfem.assembly.form.linear_form.LinearForm` wrappers +expect that the function wrapped by :class:`~skfem.assembly.form.functional.Functional` +should accept only a single argument ``w``. + +The parameter ``w`` is a dictionary containing all the default keys +(e.g., ``w['h']`` for mesh parameter and ``w['x']`` for global +coordinates) and any user provided arguments that can +be arbitrary finite element functions interpolated at the quadrature points +using :meth:`~skfem.assembly.basis.AbstractBasis.interpolate`. +As a simple example, we calculate the integral of the finite element +solution to the Poisson problem with a unit load: + +.. doctest:: + + >>> from skfem import * + >>> from skfem.models.poisson import laplace, unit_load + >>> mesh = MeshTri().refined(2).with_defaults() + >>> basis = Basis(mesh, ElementTriP2()) + >>> A = laplace.assemble(basis) + >>> b = unit_load.assemble(basis) + >>> x = solve(*condense(A, b, D=basis.get_dofs('left'))) + >>> @Functional + ... def integral(w): + ... return w['uh'] # grad, dot, etc. can be used here + >>> float(round(integral.assemble(basis, uh=basis.interpolate(x)), 5)) + 0.33333 + +.. plot:: - b(u,v) = \sum_{E \in \mathcal{E}_h} \int_{E} [u][v]\,\mathrm{d}s + from skfem import * + from skfem.models.poisson import laplace, unit_load + basis = Basis(MeshTri().refined(2).with_defaults(), ElementTriP2()) + A = laplace.assemble(basis) + b = unit_load.assemble(basis) + x = solve(*condense(A, b, D=basis.get_dofs('left'))) + from skfem.visuals.matplotlib import plot + plot(basis, x, nrefs=3, shading='gouraud', colorbar=True) -with jumps -:math:`[u] = u_1 - u_2` and :math:`[v] = v_1 - v_2` -over the interior edges can be split as +Similarly we can calculate the integral of its derivative: -.. math:: + >>> @Functional + ... def diffintegral(w): + ... return w['uh'].grad[0] # derivative wrt x + >>> float(round(diffintegral.assemble(basis, uh=basis.interpolate(x)), 5)) + 0.5 + +We can also calculate integrals over the boundary +using :class:`~skfem.assembly.basis.facet_basis.FacetBasis`: - b(u,v) = \sum_{E \in \mathcal{E}_h} \left(\int_{E} u_1 v_1\,\mathrm{d}s - \int_{E} u_1 v_2\,\mathrm{d}s - \int_{E} u_2 v_1\,\mathrm{d}s + \int_{E} u_2 v_2\,\mathrm{d}s\right) + >>> fbasis = basis.boundary('left') + >>> @Functional + ... def bndintegral(w): + ... return w['uh'].grad[1] # derivative wrt y + >>> float(round(bndintegral.assemble(fbasis, uh=fbasis.interpolate(x)), 5)) + 0.0 -and normally we would assemble all of the four forms separately. +.. _visualizing: -We can instead provide a list of bases during a call to :func:`skfem.assembly.asm`: +Visualizing the solution +======================== + +After solving the finite element system :math:`Ax=b`, it is common to +visualize the solution or other +derived fields. The library includes some basic 2D visualization +routines implemented using matplotlib. For more complex +visualizations, we suggest that the solution is saved to VTK and a +visualization is created using Paraview. + +The main routine for creating 2D visualizations using matplotlib is +:func:`skfem.visuals.matplotlib.plot` and it accepts the basis +object and the solution vector as its arguments: .. doctest:: - >>> import skfem as fem - >>> m = fem.MeshTri() - >>> e = fem.ElementTriP0() - >>> bases = [fem.InteriorFacetBasis(m, e, side=k) for k in [0, 1]] - >>> jumpform = fem.BilinearForm(lambda u, v, p: (-1) ** sum(p.idx) * u * v) - >>> fem.asm(jumpform, bases, bases).toarray() - array([[ 1.41421356, -1.41421356], - [-1.41421356, 1.41421356]]) - -For an example of practical usage, see :ref:`ex07`. + >>> from skfem import * + >>> from skfem.models.poisson import laplace, unit_load + >>> mesh = MeshTri().refined(2) + >>> basis = Basis(mesh, ElementTriP2()) + >>> A = laplace.assemble(basis) + >>> b = unit_load.assemble(basis) + >>> x = solve(*condense(A, b, D=basis.get_dofs())) + >>> from skfem.visuals.matplotlib import plot + >>> plot(basis, x) + + +.. plot:: + + from skfem import * + from skfem.models.poisson import laplace, unit_load + basis = Basis(MeshTri().refined(2), ElementTriP2()) + A = laplace.assemble(basis) + b = unit_load.assemble(basis) + x = solve(*condense(A, b, D=basis.get_dofs())) + from skfem.visuals.matplotlib import plot + plot(basis, x) + +It accepts various optional arguments to make the plots +nicer, some of which have been demonstrated in :ref:`gallery`. +For example, here is the same solution as above with different settings: + +.. doctest:: + + >>> plot(basis, x, shading='gouraud', colorbar={'orientation': 'horizontal'}, nrefs=3) + + +.. plot:: + + from skfem import * + from skfem.models.poisson import laplace, unit_load + basis = Basis(MeshTri().refined(2), ElementTriP2()) + A = laplace.assemble(basis) + b = unit_load.assemble(basis) + x = solve(*condense(A, b, D=basis.get_dofs())) + from skfem.visuals.matplotlib import plot + plot(basis, x, shading='gouraud', colorbar={'orientation': 'horizontal'}, nrefs=3) + +The routine is based on `matplotlib.pyplot.tripcolor `_ command and shares +its limitations. +For more control we suggest that the solution is saved to a VTK file for +visualization in Paraview. Saving of the solution is done through the +mesh object and it requires giving one number per node of the mesh. +Thus, for other than piecewise linear finite element bases, +the mesh must be refined and the solution interpolated for visualization: + +.. doctest:: + + >>> rmesh, rx = basis.refinterp(x, nrefs=3) # refine and interpolate + >>> rmesh.save('sol.vtk', point_data={'x': rx}) + +Another option would be to first project the solution +onto a piecewise linear finite element basis +as described in :ref:`l2proj`. +Please see :ref:`gallery` for more examples of visualization. diff --git a/docs/index.rst b/docs/index.rst index 78dc050a..4d964f22 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -3,7 +3,7 @@ ============================= `scikit-fem `_ is a pure -Python 3.7+ library for performing `finite element assembly +Python 3.8+ library for performing `finite element assembly `_. Its main purpose is the transformation of bilinear forms into sparse matrices and linear forms into vectors. The library supports triangular, quadrilateral, tetrahedral and diff --git a/docs/listofexamples.rst b/docs/listofexamples.rst index 1debffc3..b3dace55 100644 --- a/docs/listofexamples.rst +++ b/docs/listofexamples.rst @@ -1,3 +1,5 @@ +.. _gallery: + ===================== Gallery of examples ===================== diff --git a/skfem/assembly/basis/abstract_basis.py b/skfem/assembly/basis/abstract_basis.py index 3b31bb3e..81a99cb9 100644 --- a/skfem/assembly/basis/abstract_basis.py +++ b/skfem/assembly/basis/abstract_basis.py @@ -431,6 +431,11 @@ def _projection(self, interp, dtype=None): ) def project(self, interp, **kwargs): + """Perform :math:`L^2` projection onto the basis. + + See :ref:`l2proj` for more information. + + """ raise NotImplementedError def plot(self, x, visuals='matplotlib', **kwargs):