documentation etc overhaul (#69)

* improved docstrings, show functions, documentation etc. with assistance of github copilot

Author: chmerdon

CHANGELOG.md

@@ -1,5 +1,12 @@
 # CHANGES
 
+
+## v1.4.0
+
+### Changed
+  - improved docstrings, documentation
+  - improved print_convergencehistory and show functions
+
 ## v1.3.0 June 17, 2025
 
 ### Added

README.md

@@ -4,54 +4,62 @@
 [![DOI](https://zenodo.org/badge/668345991.svg)](https://zenodo.org/doi/10.5281/zenodo.10563834)
 [![code style: runic](https://img.shields.io/badge/code_style-%E1%9A%B1%E1%9A%A2%E1%9A%BE%E1%9B%81%E1%9A%B2-black)](https://github.com/fredrikekre/Runic.jl)
 
-# ExtendableFEM
-High Level API Finite Element Methods based on [ExtendableGrids.jl](https://github.com/WIAS-PDELib/ExtendableGrids.jl) (for grid management)
-and [ExtendableFEMBase.jl](https://github.com/WIAS-PDELib/ExtendableFEMBase.jl) (for finite element basis functions and dof management). 
-It offers a ProblemDescription interface, that basically involves assigning Unknowns and Operators. Such operators usually stem from a weak formulation of the problem and mainly consist of three types that can be customized via kernel functions:
+**ExtendableFEM.jl** is a high-level, extensible finite element method (FEM) library for Julia, supporting flexible problem descriptions, custom operators, and advanced grid management.
 
-- BilinearOperator,
-- LinearOperator,
-- NonlinearOperator (that automatically assemble Newton's method by automatic differentiation)
+## Features
 
-### Quick Example
+- High-level, extensible API for solving PDE problems by finite element methods
+- Flexible `ProblemDescription` interface for assigning unknowns and operators
+- Supports custom kernel functions for bilinear, linear, and nonlinear forms
+- Automatic assembly and Newton's method for NonlinearOperators
+- Builds upon [ExtendableGrids.jl](https://github.com/WIAS-PDELib/ExtendableGrids.jl) and low level structures from [ExtendableFEMBase.jl](https://github.com/WIAS-PDELib/ExtendableFEMBase.jl)
 
-The following minimal example demonstrates how to setup a Poisson problem.
+## Quick Example
+
+The following minimal example demonstrates how to set up and solve a Poisson problem:
 
 ```julia
 using ExtendableFEM
 using ExtendableGrids
 
-# build/load any grid (here: a uniform-refined 2D unit square into triangles)
+# Build a uniform-refined 2D unit square grid with triangles
 xgrid = uniform_refine(grid_unitsquare(Triangle2D), 4)
 
-# create empty PDE description
+# Create a new PDE description
 PD = ProblemDescription()
 
-# create and assign unknown
+# Define and assign the unknown
 u = Unknown("u"; name = "potential")
 assign_unknown!(PD, u)
 
-# assign Laplace operator
+# Assign Laplace operator (diffusion term)
 assign_operator!(PD, BilinearOperator([grad(u)]; factor = 1e-3))
 
-# assign right-hand side data
+# Assign right-hand side data
 function f!(fval, qpinfo)
     x = qpinfo.x # global coordinates of quadrature point
     fval[1] = x[1] * x[2]
 end
 assign_operator!(PD, LinearOperator(f!, [id(u)]))
 
-# assign boundary data (here: u = 0)
+# Assign Dirichlet boundary data (u = 0)
 assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4))
 
-# discretise = choose FESpace
+# Discretize: choose finite element space
 FEType = H1Pk{1,2,3} # cubic H1-conforming element with 1 component in 2D
 FES = FESpace{FEType}(xgrid)
 
-# solve
-sol = solve!(Problem, [FES])
+# Solve the problem
+sol = solve!(PD, [FES])
 
-# plot
+# Plot the solution
 using PyPlot
 plot(id(u), sol; Plotter = PyPlot)
 ```
+
+## Citing
+
+If you use ExtendableFEM.jl in your research, please cite [this Zenodo record](https://zenodo.org/doi/10.5281/zenodo.10563834).
+
+## License
+ExtendableFEM.jl is licensed under the MIT License. See [LICENSE](https://github.com/WIAS-PDELib/ExtendableFEM.jl/blob/master/LICENSE) for details.

docs/src/bilinearoperator.md

@@ -1,16 +1,12 @@
-
 # BilinearOperator
 
-A bilinear operator allows to add matrices to the system matrix that usually refer to
-linearisations of the PDE operators or stabilisations. If the bilinear operator
-lives on face entities, also jumps of operators can be involved, if they are naturally
-continuous for the finite element space in operation (also jumps for broken spaces)
-and only involve degrees of freedom on the face, e.g.
-normal jumps for Hdiv spaces or jumps for H1-conforming spaces or tangential jumps
-of Hcurl spaces. For all other discontinuous operator evaluations
-(that needs to evaluat more than the degrees of freedom on the face)
-there is the possibility to use BilinearOperatorDG.
-It is also possible to assign a matrix assembled by the user as a BilinearOperator.
+`BilinearOperator` provides a flexible interface for defining bilinear forms (matrices) in finite element problems. These operators typically represent (linearizations of) PDE operators, stabilization terms, or custom user-defined couplings. The interface supports both standard and discontinuous Galerkin (DG) forms, and allows for custom assembly on cells, faces, or other grid entities.
+
+## Overview
+
+- **Standard Bilinear Operators:** Assemble contributions to the system matrix, including diffusion, mass, and convection terms, as well as custom couplings. Supports both cell- and face-based assembly, including jump terms for conforming and piecewise spaces that only require the dofs on that entity.
+- **Discontinuous Galerkin (DG) Operators:** Use `BilinearOperatorDG` for forms that require evaluation of jumps or averages involving all degrees of freedom on neighboring cells (e.g., interior penalty stabilization, gradient jumps for H1 elements).
+- **Custom Matrices:** You can also assign a user-assembled matrix as a `BilinearOperator`.
 
 ## Constructors
 
@@ -22,10 +18,7 @@ Order   = [:type, :function]
 
 ## BilinearOperatorDG
 
-BilinearOperatorDG is intended for bilinear forms that involves jumps of discontinuous quantities
-on faces whose assembly requires evaluation of all degrees of freedom on the neighbouring cells,
-e.g. gradient jumps for H1 conforming functions. In this case the assembly loop triggers
-integration along the boundary of the cells.
+`BilinearOperatorDG` is intended for bilinear forms that involve jumps or averages of discontinuous quantities on faces, requiring access to all degrees of freedom on neighboring cells. This is essential for DG methods and certain stabilization techniques.
 
 ```@autodocs
 Modules = [ExtendableFEM]
@@ -35,11 +28,10 @@ Order   = [:type, :function]
 
 ## Examples
 
-Below two examples illustrate some use cases.
+### Example: Stokes Operator
 
-### Example - Stokes operator
+For the linear operator of a Stokes problem, a kernel could look like:
 
-For the linear operator of a Stokes problem a kernel could look like
 ```julia
 μ = 0.1 # viscosity parameter
 function kernel!(result, input, qpinfo)
@@ -52,26 +44,29 @@ function kernel!(result, input, qpinfo)
     return nothing
 end
 ```
-and the coressponding BilinearOperator constructor call reads
+
+The corresponding `BilinearOperator` constructor call reads:
+
 ```julia
 u = Unknown("u"; name = "velocity")
 p = Unknown("p"; name = "pressure")
 BilinearOperator(kernel!, [grad(u), id(p)]; use_sparsity_pattern = true)
 ```
-The additional argument causes that the zero pressure-pressure block of the matrix is not (even tried to be) assembled,
-since ```input[5]``` does not couple with ```result[5]```.
 
+The `use_sparsity_pattern` argument ensures that the zero pressure-pressure block of the matrix is not assembled, since `input[5]` does not couple with `result[5]`.
+
+### Example: Interior Penalty Stabilization
 
-### Example - interior penalty stabilization
+A popular stabilization for convection-dominated problems is based on jumps of the gradient, which can be realized with the following kernel:
 
-A popular convection stabilization is based on the jumps of the gradient, which can be realised with
-the kernel
 ```julia
 function stab_kernel!(result, input, qpinfo)
     result .= input .* qpinfo.volume^2
 end
 ```
-and the BilinearOperatorDG constructor call
+
+and the `BilinearOperatorDG` constructor call:
+
 ```julia
 u = Unknown("u")
 assign_operator!(PD, BilinearOperatorDG(stab_kernel!, [jump(grad(u))]; entities = ON_IFACES, factor = 0.01))

docs/src/callbackoperator.md

@@ -1,13 +1,37 @@
-
 # CallbackOperator
 
-A callback operator passes the matrix and rhs to a user-defined function where they can be modified as desired.
-An example where this is used is Example265.
+`CallbackOperator` provides a flexible interface for injecting custom assembly logic into the finite element workflow. It delegates the assembly of the system matrix and right-hand side to a user-supplied callback function, allowing for advanced or nonstandard modifications that are difficult to express with standard operators.
 
-## Constructors
+## Constructor
 
 ```@autodocs
 Modules = [ExtendableFEM]
 Pages = ["common_operators/callback_operator.jl"]
 Order   = [:type, :function]
 ```
+
+## Example
+
+```julia
+function my_callback!(A, b, args; assemble_matrix=true, assemble_rhs=true, time=0, kwargs...)
+    # Example: add a constant to the diagonal
+    if assemble_matrix && A !== nothing
+        for i in 1:min(size(A)...)
+            A[i, i] += 1.0
+        end
+    end
+    if assemble_rhs && b !== nothing
+        b .+= 2.0
+    end
+end
+op = CallbackOperator(my_callback!; u_args=[1], name="CustomOp")
+```
+
+## When to Use
+
+- Custom or experimental PDE terms
+- Coupling to external codes or data
+- Advanced boundary or interface conditions
+- Prototyping new assembly strategies
+
+See also: [Example265](https://wias-pdelib.github.io/ExtendableFEM.jl/stable/examples/) for a practical use case.

docs/src/combinedofs.md

@@ -1,17 +1,39 @@
-
 # CombineDofs
 
+`CombineDofs` provides a mechanism to couple degrees of freedom (DOFs) in a finite element system. This is especially useful for enforcing periodic boundary conditions, multi-point constraints, or other situations where DOFs from different locations or unknowns should be treated as identical in the assembled system.
+
 ```@autodocs
 Modules = [ExtendableFEM]
 Pages = ["common_operators/combinedofs.jl"]
 Order   = [:type, :function]
 ```
 
-The following function might be useful to find out the dofs
-the need to be coupled for periodic
-boundary conditions:
+## Periodic Boundary Conditions
 
+To set up periodic boundary conditions, you often need to determine which DOFs should be coupled. The following utility functions can help:
 
 ```@docs
 get_periodic_coupling_info
+get_periodic_coupling_matrix
 ```
+
+## Example: Periodic Boundary Coupling (extract from Example212)
+
+Suppose you want to enforce periodicity between the left and right boundaries of a 2D domain. You can use `get_periodic_coupling_matrix` to find the corresponding DOFs, and then use `CombineDofs` to couple them in the system. The following code is adapted from [Example212_PeriodicBoundary2D.jl](https://github.com/WIAS-PDELib/ExtendableFEM.jl/blob/main/examples/Example212_PeriodicBoundary2D.jl):
+
+```julia
+# Define a function to map points on the left to the right boundary
+function give_opposite!(y, x)
+    y .= x
+    y[1] = width - x[1]
+    return nothing
+end
+
+# Compute the coupling matrix
+coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!)
+
+# Assign the CombineDofs operator to the problem description
+assign_operator!(PD, CombineDofs(u, u, coupling_matrix))
+```
+
+See also [Example212](https://wias-pdelib.github.io/ExtendableFEM.jl/stable/examples/) and [Example312](https://wias-pdelib.github.io/ExtendableFEM.jl/stable/examples/) for more advanced use cases and details on periodic boundary conditions.

docs/src/examples_intro.md

@@ -1,29 +1,29 @@
-# About the examples
-
-The examples have been designed with the following issues in mind:
-- they run from the Julia REPL
-- each example is a Julia module named similar to the basename of the example file.
-- an example can be used as the starting point for a project 
-- some examples define test cases for the test suite
-- ExampleXYZ with X = A can be considered advanced and uses low-level structures
-  and/or demonstrates customisation features or experimental features
-- the default output of the main function is printed on the website and can be
-  used to check if the code runs as expected (unfortunately REPL messages are not recorded)
-- printed assembly and solving times (especially in a first iteration) can be much larger due to first-run compilation times,
-  the printouts in the documentation are taken from a second run after compilations are done
-
-
-## Running the examples
-
-In order to run `ExampleXXX`, perform the following steps:
-
-- Download the example file (e.g. via the source code link at the top)
-- Make sure all used packages are installed in your Julia environment
-- In the REPL: 
-```
-julia> include("ExampleXXX.jl")`
-
-julia> ExampleXXX.main()
-```
-- Some examples offer visual output via the optional argument Plotter = PyPlot or Plotter = GLMakie
-(provided the package PyPlot/GLMakie is installed and loaded)
+# About the Examples
+
+The examples in this package are designed to be practical, reproducible, and educational. They demonstrate a wide range of finite element applications and PDE model problems.
+
+## Design Principles
+
+- All examples can be run directly from the Julia REPL.
+- Each example is a Julia module named after the file.
+- Examples can serve as templates for your own projects.
+- Many examples include test cases for automated verification.
+
+## Running the Examples
+
+To run an example (e.g., `Example212_PeriodicBoundary2D`):
+
+1. Download the example file (see the source code link at the top of the example page).
+2. Ensure all required packages are installed in your Julia environment.
+3. In the Julia REPL:
+
+    ```julia
+    julia> include("Example212_PeriodicBoundary2D.jl")
+    julia> Example212_PeriodicBoundary2D.main()
+    ```
+
+4. Some examples offer visual output via the optional argument `Plotter = PyPlot` or `Plotter = GLMakie` (provided the package is installed and loaded):
+
+    ```julia
+    julia> Example212_PeriodicBoundary2D.main(Plotter = PyPlot)
+    ```

docs/src/faceinterpolator.md

@@ -1,7 +1,26 @@
-# Face interpolator
+# Face Interpolator
+
+The face interpolator provides tools for evaluating and interpolating finite element functions on mesh faces. This is particularly useful for postprocessing tasks that require access to traces or jumps of functions across element boundaries.
+
+## API Reference
 
 ```@autodocs
 Modules = [ExtendableFEM]
 Pages = ["common_operators/discface_interpolator.jl"]
 Order   = [:type, :function]
 ```
+
+## Example Usage (extracted from Example210)
+
+Suppose you want to compute the jumps of the gradient of a scalar-valued
+Lagrange finite element function on the interior edges, e.g. to compute an a posteriori error estimator.
+
+```julia
+function gradnormalflux!(result, ∇u, qpinfo)
+    result[1] = dot(∇u, qpinfo.normal)
+end
+NormalJumpProjector = FaceInterpolator(gradnormalflux!, [jump(grad(u))]; resultdim = 1, only_interior = true)
+Jumps4Faces = evaluate!(NormalJumpProjector, sol)
+```
+
+See the [Example212](https://wias-pdelib.github.io/ExtendableFEM.jl/stable/examples/) for the complete example.

docs/src/fixdofs.md

@@ -1,6 +1,9 @@
-
 # FixDofs
 
+`FixDofs` provides a mechanism to fix (constrain) specific degrees of freedom (DOFs) in a finite element system to given values.
+
+## API Reference
+
 ```@autodocs
 Modules = [ExtendableFEM]
 Pages = ["common_operators/fixdofs_operator.jl"]

docs/src/homogeneousdata.md

@@ -1,8 +1,18 @@
-
 # HomogeneousData
 
+`HomogeneousData` provides a convenient way to enforce homogeneous (zero) Dirichlet boundary conditions or constraints in a finite element problem. It automatically sets the solution to zero on specified regions or entities.
+
+## API Reference
+
 ```@autodocs
 Modules = [ExtendableFEM]
 Pages = ["common_operators/homogeneousdata_operator.jl"]
 Order   = [:type, :function]
 ```
+
+## Example: Zero Dirichlet Boundary Condition
+
+```julia
+# Impose u = 0 on boundary region 1
+assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1]))
+```