WIAS-PDELib · chmerdon · Jul 23, 2025 · Jul 14, 2025 · Jul 16, 2025 · Jul 23, 2025
diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "ExtendableFEM"
 uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
-version = "1.3.0"
+version = "1.4"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
 
 [deps]

diff --git a/docs/src/itemintegrators.md b/docs/src/itemintegrators.md
@@ -20,4 +20,4 @@ Pages = ["common_operators/item_integrator_dg.jl"]
 Order   = [:type, :function]
 ```
 
-See the [examples](https://wias-pdelib.github.io/ExtendableFEM.jl/stable/examples/) for some practical use cases.
+See, e.g., Example207, Example210 or Example245 for some practical use cases.
diff --git a/docs/src/nonlinearoperator.md b/docs/src/nonlinearoperator.md
@@ -1,7 +1,7 @@
 
 # NonlinearOperator
 
-`NonlinearOperator` provides a high-level interface for defining nonlinear forms in finite element problems. It automatically assembles all terms required for Newton's method, including the residual and Jacobian, using user-supplied kernel functions. For alternative linearizations, use `BilinearOperator` or `LinearOperator` with appropriate kernels.
+`NonlinearOperator` provides a high-level interface for defining nonlinear forms in finite element problems. It automatically assembles all terms required for Newton's method, including the residual and Jacobian, using user-supplied kernel functions. For alternative linearizations, use [`BilinearOperator`](@ref "BilinearOperator") or [`LinearOperator`](@ref "LinearOperator") with appropriate kernels.
 
 
 ## Constructor

diff --git a/docs/src/pdesolvers.md b/docs/src/pdesolvers.md
@@ -13,7 +13,7 @@ To solve a `ProblemDescription`, you must provide discretization information:
     ```julia
     FESpace{FEType}(grid::ExtendableGrid)
     ```
-    where `FEType` is the finite element type. See the [list of available FETypes](https://wias-pdelib.github.io/ExtendableFEMBase.jl/dev/fems/) in the [ExtendableFEMBase.jl documentation](https://wias-pdelib.github.io/ExtendableFEMBase.jl/dev/).
+    where `FEType` is the finite element type. See the [list of available FETypes](https://wias-pdelib.github.io/ExtendableFEMBase.jl/stable/fems/) in the [ExtendableFEMBase.jl documentation](https://wias-pdelib.github.io/ExtendableFEMBase.jl/stable/).
 
 ## Monolithic Solve
 

diff --git a/docs/src/pdesolvers_dt.md b/docs/src/pdesolvers_dt.md
@@ -55,9 +55,3 @@ Given a finite element space `FES` and an initial `FEVector` `sol` for the unkno
 prob = generate_ODEProblem(PD, FES, (0, T); init = sol)
 DifferentialEquations.solve(prob, ImplicitEuler(autodiff = false), dt = 1e-3, dtmin = 1e-6, adaptive = true)
 ```
-
-## Tips
-
-- For more advanced time-stepping schemes, manage previous solutions and time-derivative terms manually in the `ProblemDescription`.
-- See the [examples](https://wias-pdelib.github.io/ExtendableFEM.jl/stable/examples/) for practical implementations of time-dependent problems.
-- For further details on the ODE interface, see the [ExtendableFEMDiffEQExt.jl documentation](https://github.com/WIAS-PDELib/ExtendableFEMDiffEQExt.jl).
diff --git a/docs/src/problemdescription.md b/docs/src/problemdescription.md
@@ -34,15 +34,15 @@ Operators define how terms contribute to the system matrix or right-hand side. T
 ### Types of Operators
 
 The main operator classes are:
-- `NonlinearOperator` (e.g., nonlinear convection in Navier–Stokes)
-- `BilinearOperator` (e.g., Laplacian in Poisson)
-- `LinearOperator` (e.g., right-hand side in Poisson or Navier–Stokes)
+- [NonlinearOperator](@ref) (e.g., nonlinear convection in Navier–Stokes)
+- [BilinearOperator](@ref) (e.g., Laplacian in Poisson)
+- [LinearOperator](@ref) (e.g., right-hand side in Poisson or Navier–Stokes)
 
 For boundary conditions or global constraints, use:
-- `InterpolateBoundaryData`
-- `HomogeneousBoundaryData`
-- `FixDofs`
-- `CombineDofs`
+- [InterpolateBoundaryData](@ref)
+- [HomogeneousData](@ref)
+- [FixDofs](@ref)
+- [CombineDofs](@ref)
 
 ### Entities and Regions
 
@@ -55,8 +55,8 @@ Each operator assembles on certain mesh entities. The default is cell-wise assem
 | ON_FACES         | Assemble/interpolate on all mesh faces                           |
 | ON_IFACES        | Assemble/interpolate on interior mesh faces                      |
 | ON_BFACES        | Assemble/interpolate on boundary mesh faces                      |
-| ON_EDGES (*)     | Assemble/interpolate on all mesh edges (3D only, experimental)   |
-| ON_BEDGES (*)    | Assemble/interpolate on boundary mesh edges (3D only, experimental) |
+| ON_EDGES (*)     | Assemble/interpolate on all mesh edges   |
+| ON_BEDGES (*)    | Assemble/interpolate on boundary mesh edges |
 
 !!! note
     (*) = Only reasonable in 3D and still experimental; may have some issues.

diff --git a/examples/Example201_PoissonProblem.jl b/examples/Example201_PoissonProblem.jl
@@ -29,13 +29,24 @@ using Test #hide
 ## define variables
 u = Unknown("u"; name = "potential")
 
-## define data functions
-function f!(fval, qpinfo)
+## default right-hand side f(x,y) = xy
+function default_f!(fval, qpinfo)
     fval[1] = qpinfo.x[1] * qpinfo.x[2]
     return nothing
 end
 
-function main(; μ = 1.0, nrefs = 4, order = 2, Plotter = nothing, parallel = false, npart = parallel ? 8 : 1, kwargs...)
+@info "This example solve the 2D Poisson problem -μΔu = f.
+    To run it call Example201_PoissonProblem.main(; kwargs...) from the Julia REPL.
+
+    Most important kwargs:
+    μ       - diffusion coefficient (default: 1.0)
+    f!      - right-hand side function (default: (fval, qpinfo) -> (fval[1] = qpinfo.x[1] * qpinfo.x[2]))
+    nrefs   - number of uniform grid refinements (default: 4)
+    order   - polynomial order of the finite element space (default: 2)
+    Plotter - plotting backend (e.g. GLMakie, PyPlot, default: nothing)
+    parallel- enable parallel assembly (default: false)
+    npart   - number of partitions for parallel run (default: 8 if parallel, else 1)"
+function main(; μ = 1.0, nrefs = 4, order = 2, Plotter = nothing, parallel = false, f! = default_f!, npart = parallel ? 8 : 1, kwargs...)
     ## problem description
     PD = ProblemDescription()
     assign_unknown!(PD, u)

diff --git a/src/common_operators/bilinear_operator.jl b/src/common_operators/bilinear_operator.jl
@@ -108,9 +108,7 @@ Constructs a bilinear operator from a user-supplied matrix `A`, which can be a s
 - `A::AbstractMatrix`: The matrix representing the bilinear form, e.g., a sparse matrix.
 - `u_test::Vector{<:Union{Unknown, Int}}`: Identifiers or indices for the test functions.
 - `u_ansatz::Vector{<:Union{Unknown, Int}}` (optional): Identifiers or indices for the ansatz (trial) functions. Defaults to `u_test`.
-
-# Keyword Arguments
-$(ExtendableFEMBase._myprint(default_blfop_kwargs()))
+- `u_args::Vector{<:Union{Unknown, Int}}` (optional): Identifiers or indices for unknowns on which the matrix depends in a nonlinear way (this tells the solver which solution blocks trigger reassembly).
 
 # Returns
 A `BilinearOperatorFromMatrix` object that can be used in the assembly process.

diff --git a/src/common_operators/bilinear_operator_dg.jl b/src/common_operators/bilinear_operator_dg.jl
@@ -159,23 +159,30 @@ function BilinearOperatorDG(
 Generates a bilinear form that evaluates the vector product of the (discontinuous)
 operator evaluation(s) of the test function(s) with the (discontinuous) operator evaluation(s)
 of the ansatz function(s). If a function is provided in the first argument,
-the ansatz function evaluations can be customized by the kernel function
+the ansatz function evaluations can be customized by the kernel function,
 and its result vector is then used in a dot product with the test function evaluations.
-In this case the header of the kernel functions needs to be conform
-to the interface
+In this case, the header of the kernel function must conform to the interface:
 
-	kernel!(result, eval_ansatz, qpinfo)
+    kernel!(result, eval_ansatz, qpinfo)
 
-where qpinfo allows to access information at the current quadrature point.
+where `qpinfo` allows access to information at the current quadrature point.
 
 Operator evaluations are tuples that pair an unknown identifier or integer
-with a Function operator.
+with a function operator (such as those generated by `jump(grad(u))`, `id(u)`, etc).
 
-Example: BilinearOperatorDG([jump(grad(1))], [jump(grad(1))]; kwargs...) generates an interior penalty stabilisation.
+# Arguments
+- `oa_test`: Array of tuples `(unknown, operator)` for test functions.
+- `oa_ansatz`: Array of tuples `(unknown, operator)` for ansatz (trial) functions. Defaults to `oa_test`.
+- `kwargs...`: Additional keyword arguments for operator assembly (see below).
 
-Keyword arguments:
+# Keyword arguments
 $(_myprint(default_blfopdg_kwargs()))
 
+# Example
+```julia
+BilinearOperatorDG([jump(grad(1))], [jump(grad(1))]; entities=ON_FACES)
+```
+
 """
 function BilinearOperatorDG(oa_test::Array{<:Tuple{Union{Unknown, Int}, DataType}, 1}, oa_ansatz::Array{<:Tuple{Union{Unknown, Int}, DataType}, 1} = oa_test; kwargs...)
     u_test = [oa[1] for oa in oa_test]
@@ -375,7 +382,7 @@ function build_assembler!(A, O::BilinearOperatorDG{Tv}, FE_test, FE_ansatz, FE_a
                     couple = false
                     for j in 1:op_lengths_ansatz[id]
                         for k in 1:op_lengths_test[idt]
-                            if sparsity_pattern[k + op_offsets_test[idt], j + op_offsets_ansatz[id]] > 0
+                            if sparsity_pattern[k + op_offsets_test[id], j + op_offsets_ansatz[idt]] > 0
                                 couple = true
                             end
                         end

diff --git a/src/common_operators/linear_operator.jl b/src/common_operators/linear_operator.jl
@@ -171,9 +171,6 @@ This wraps an explicit vector `b` (which can be an `AbstractVector` or a block-s
 - `b`: The vector to use as the right-hand side. Can be an `AbstractVector` or an `FEVector` with multiple blocks.
 - `u_test`: Array specifying which unknowns or blocks the vector should be assigned to.
 
-# Keyword Arguments
-$(_myprint(default_linop_kwargs()))
-
 """
 function LinearOperator(b, u_test; kwargs...)
     parameters = Dict{Symbol, Any}(k => v[1] for (k, v) in default_linop_kwargs())
@@ -193,9 +190,6 @@ This wraps an explicit matrix `A` (which can be an `AbstractMatrix` or a block-s
 - `u_test`: Array specifying which unknowns or blocks the result should be assigned to.
 - `u_args`: Array specifying which unknowns or blocks are used as arguments (i.e., which solution coefficients are multiplied by `A`).
 
-# Keyword Arguments
-$(_myprint(default_linop_kwargs()))
-
 """
 function LinearOperator(A::AbstractMatrix, u_test::Array{<:Union{Unknown, Int}, 1}, u_args::Array{<:Union{Unknown, Int}, 1}; kwargs...)
     parameters = Dict{Symbol, Any}(k => v[1] for (k, v) in default_linop_kwargs())