Merge branch 'master' into chmerdon/diffinterface

chmerdon · web-flow · commit f4f4a38a369f · 2025-08-04T23:47:29.000+09:00
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "ExtendableFEM"
 uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
+version = "1.4"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
-version = "1.3.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -34,7 +34,7 @@ DifferentiationInterface = "0.7.4"
 DocStringExtensions = "0.8,0.9"
 ExampleJuggler = "2.2.1"
 ExplicitImports = "1"
-ExtendableFEMBase = "1"
+ExtendableFEMBase = "1.3.0"
 ExtendableGrids = "1.10.3"
 ExtendableSparse = "1.5.3"
 ForwardDiff = "0.10.35,1"
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,13 @@ 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...)
+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/.JuliaFormatter.toml b/src/.JuliaFormatter.toml
diff --git a/src/ExtendableFEM.jl b/src/ExtendableFEM.jl
@@ -67,7 +67,7 @@ using LinearSolve: LinearSolve, LinearProblem, UMFPACKFactorization, deleteat!,
     init, solve
 using Printf: Printf, @printf, @sprintf
 using SparseArrays: SparseArrays, AbstractSparseArray, SparseMatrixCSC, findnz, nnz,
-    nzrange, rowvals, sparse
+    nzrange, rowvals, sparse, SparseVector
 using ADTypes: ADTypes, KnownJacobianSparsityDetector
 using SparseConnectivityTracer: SparseConnectivityTracer, TracerSparsityDetector
 using DifferentiationInterface: DifferentiationInterface,
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]
@@ -790,7 +797,7 @@ function build_assembler!(A, O::BilinearOperatorDG{Tv}, FE_test, FE_ansatz; time
                     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[id], j + op_offsets_ansatz[idt]] > 0
+                            if sparsity_pattern[k + op_offsets_test[idt], j + op_offsets_ansatz[id]] > 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())
diff --git a/src/helper_functions.jl b/src/helper_functions.jl
@@ -244,8 +244,9 @@ function _get_periodic_coupling_matrix(
     # to be sure
     fill!(fe_vector.entries, 0.0)
 
-    # FE vector for interpolation
-    fe_vector_target = FEVector(FES)
+    # FE target vector for interpolation with sparse entries
+    fe_vector_target = FEVector(FES; entries = SparseVector{Float64, Int64}(FES.ndofs, Int64[], Float64[]))
+
 
     # resulting sparse matrix
     n = length(fe_vector.entries)
@@ -358,6 +359,13 @@ function _get_periodic_coupling_matrix(
         append!(searchareas, view(faces_to, 1:nfaces_to))
     end
 
+    # throw error if no search area had been found for a bface
+    for source in 1:num_sources(searchareas)
+        if num_targets(searchareas, source) == 0
+            throw("bface $source has no valid search area on the opposite side of the grid. Are from/to regions and give_opposite! function correct?")
+        end
+    end
+
     # loop over boundary face indices: we need this index for dofs_on_boundary
     for i_boundary_face in 1:n_boundary_faces
 
@@ -372,7 +380,8 @@ function _get_periodic_coupling_matrix(
                 end
 
                 # reset
-                fill!(fe_vector_target.entries, 0.0)
+                empty!(fe_vector_target.entries.nzind)
+                empty!(fe_vector_target.entries.nzval)
 
                 # activate one entry
                 fe_vector.entries[local_dof] = 1.0
@@ -394,7 +403,7 @@ function _get_periodic_coupling_matrix(
                 fe_vector.entries[local_dof] = 0.0
 
                 # set entries
-                for (i, target_entry) in enumerate(fe_vector_target.entries)
+                for (i, target_entry) in zip(findnz(fe_vector_target.entries)...)
                     if abs(target_entry) > sparsity_tol
                         result[local_dof, i] = target_entry
                     end
diff --git a/src/solvers.jl b/src/solvers.jl
