From c2f3a3996574f077ed5012ddd7ffd82e8ec57b2f Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Tue, 24 Jun 2025 22:04:02 +0200 Subject: [PATCH 1/8] documentation overhaul, supported by github copilot --- CHANGELOG.md | 3 + Project.toml | 2 +- docs/src/examples_intro.md | 32 +- docs/src/feevaluator.md | 3 +- docs/src/fematrix.md | 2 +- docs/src/fems.md | 34 +- docs/src/fespace.md | 2 +- docs/src/fevector.md | 2 +- docs/src/functionoperators.md | 32 +- docs/src/index.md | 4 +- docs/src/interpolations.md | 8 +- docs/src/meshing.md | 8 +- docs/src/plots.md | 12 +- docs/src/pointevaluators.md | 2 +- docs/src/quadrature.md | 13 +- docs/src/segmentintegrators.md | 2 +- examples/Example200_LowLevelPoisson.jl | 91 ++--- src/ExtendableFEMBase.jl | 56 ++- src/dofmaps.jl | 18 +- src/feevaluator.jl | 40 ++- src/fematrix.jl | 100 +++--- src/fevector.jl | 85 +++-- src/finiteelements.jl | 98 ++++-- src/interpolations.jl | 450 ++++++++++++++++--------- src/interpolators.jl | 93 ++++- src/lazy_interpolate.jl | 46 ++- src/point_evaluator.jl | 42 ++- src/qpinfos.jl | 25 +- src/quadrature.jl | 71 ++-- src/segment_integrator.jl | 36 +- 30 files changed, 918 insertions(+), 494 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9d39de0..9dca4c5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,8 @@ # CHANGES +## v1.2.0 June 24, 2025 + - major documentation overhaul + ## v1.1.1 June 16, 2025 - docstring improvements diff --git a/Project.toml b/Project.toml index 6729757..9beed4e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ExtendableFEMBase" uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c" authors = ["Christian Merdon ", "Patrick Jaap "] -version = "1.1.1" +version = "1.2.0" [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" diff --git a/docs/src/examples_intro.md b/docs/src/examples_intro.md index b8f27ca..3756002 100644 --- a/docs/src/examples_intro.md +++ b/docs/src/examples_intro.md @@ -1,28 +1,24 @@ -# About the examples +# 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 provided examples are designed with the following goals in mind: +- They can be run directly from the Julia REPL. +- Each example is implemented as a Julia module, named similarly to the file's basename. +- Examples can serve as starting points for your own projects. +- Some examples include test cases that are integrated into the test suite. +- Assembly and solve times (especially for the first run) may be significantly higher due to Julia's just-in-time compilation. ## 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")` +- Download the example file (e.g., using the source code link at the top of the page). +- Ensure all required packages are installed in your Julia environment. +- In the Julia REPL, run: +```julia +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) + +- Some examples provide visual output via the optional argument `Plotter = PyPlot` or `Plotter = GLMakie` (these require the corresponding plotting package to be installed and loaded). diff --git a/docs/src/feevaluator.md b/docs/src/feevaluator.md index 87d5f01..cd60eb8 100644 --- a/docs/src/feevaluator.md +++ b/docs/src/feevaluator.md @@ -1,7 +1,8 @@ # FEEvaluator -FEEvaluators provide a structure that handles the evaluation of finite element basis functions for a given function operator, quadrature rule and item geometry. It stores the evaluations on the reference geometry (where derivatives are computed by automatic differentiation) and on the current mesh item. The current mesh item can be changed via the update! call. +The `FEEvaluator` provides a unified interface for evaluating finite element basis functions, their derivatives, and related quantities for a given function operator, quadrature rule, and mesh entity. It manages the storage and reuse of basis evaluations both on the reference element (where derivatives are computed via automatic differentiation) and on the current mesh item. The mesh item context can be updated dynamically using the `update!` function. + ```@autodocs Modules = [ExtendableFEMBase] diff --git a/docs/src/fematrix.md b/docs/src/fematrix.md index d3edd8e..834fcfc 100644 --- a/docs/src/fematrix.md +++ b/docs/src/fematrix.md @@ -1,6 +1,6 @@ # FEMatrix -A FEMatrix consists of FEMatrixBlocks that share a common ExtendableSparseMatrix. Each block is associated to two FESpaces and can only write into a submatrix of the common sparse matrix specified by offsets. +A `FEMatrix` represents a block-structured finite element matrix, where each block (a `FEMatrixBlock`) corresponds to a pair of finite element spaces and operates on a submatrix of a shared `ExtendableSparseMatrix`. ```@autodocs Modules = [ExtendableFEMBase] diff --git a/docs/src/fems.md b/docs/src/fems.md index 2eb5312..f13e7bc 100644 --- a/docs/src/fems.md +++ b/docs/src/fems.md @@ -1,13 +1,10 @@ - # Implemented Finite Elements -This page describes the finite element type-tree and lists all implemented finite elements. - - +This page provides an overview of the finite element type hierarchy and lists all finite elements currently implemented in ExtendableFEMBase. -## The Finite Element Type-Tree +## The Finite Element Type Hierarchy -Finite elements are abstract type leaves in a type-tree. The complete tree looks like this: +Finite elements in ExtendableFEMBase are organized as leaves in an abstract type hierarchy. The complete type tree is as follows: ``` AbstractFiniteElement @@ -40,21 +37,20 @@ AbstractFiniteElement #### Remarks -- each type depends on one/two or three parameters, the first one is always the number of components (ncomponents) that determines if the - finite element is scalar- or veector-valued; some elements additionally require the parameter edim <: Int if they are structurally different in different space dimensions; arbitrary order elements require a third parameter that determines the order -- each finite elements mainly comes with a set of basis functions in reference coordinates for each applicable AbstractElementGeometry and degrees of freedom maps for each mesh entity -- broken finite elements are possible via the broken switch in the [FESpace](@ref) constructor -- the type steers how the basis functions are transformed from local to global coordinates and how FunctionOperators are evaluated -- depending on additional continuity properties of the element types more basis function sets are defined: - - AbstractH1FiniteElements additionally have evaluations of nonzero basisfunctions on faces/bfaces - - AbstractHdivFiniteElements additionally have evaluations of nonzero normalfluxes of basisfunctions on faces/bfaces - - AbstractHcurlFiniteElements additionally have evaluations of nonzero tangentfluxes of basisfunctions on edges/bedges -- each finite element has its own implemented standard interpolation interpolate! (see [Finite Element Interpolations](@ref)) that can be applied to a function with header function(result, qpinfo), below it is shortly described what this means for each finite element +- Each finite element type depends on one, two, or three parameters. The first parameter is always the number of components (`ncomponents`), which determines whether the element is scalar- or vector-valued. Some elements also require the parameter `edim <: Int` if their structure differs by spatial dimension. Arbitrary order elements require a third parameter specifying the polynomial order. +- Each finite element provides a set of basis functions in reference coordinates for each applicable `AbstractElementGeometry`, as well as degree-of-freedom (dof) maps for each mesh entity. +- Discontinuous (broken) finite elements can be created using the `broken` switch in the [`FESpace`](@ref) constructor. +- The element type determines how basis functions are transformed from local to global coordinates and how `FunctionOperators` are evaluated. +- Additional continuity properties of element types lead to more specialized basis function sets: + - `AbstractH1FiniteElement` types provide evaluations of nonzero basis functions on faces/boundary faces. + - `AbstractHdivFiniteElement` types provide evaluations of nonzero normal fluxes of basis functions on faces/boundary faces. + - `AbstractHcurlFiniteElement` types provide evaluations of nonzero tangential fluxes of basis functions on edges/boundary edges. +- Each finite element has its own standard interpolation routine `interpolate!` (see [Finite Element Interpolations](@ref)), which can be applied to a function with the signature `function(result, qpinfo)`. The specific interpolation behavior is described for each element below. -## List of implemented Finite Elements +## List of Implemented Finite Elements -The following table lists all currently implemented finite elements and on which geometries they are available (in brackets a dofmap pattern for CellDofs is shown and the number of local degrees of freedom for a vector-valued realisation). Click on the FEType to find out more details. +The following table summarizes all finite elements currently implemented in ExtendableFEMBase and indicates the reference geometries on which they are available. For each entry, the dofmap pattern for cell degrees of freedom is shown in brackets, along with the number of local degrees of freedom for a vector-valued realization. Click on an FEType to view more details. | FEType | Triangle2D | Parallelogram2D | Tetrahedron3D | Parallelepiped3D | | :----------------: | :----------------: | :----------------: | :----------------: | :----------------: | @@ -86,7 +82,7 @@ The following table lists all currently implemented finite elements and on which | [`HDIVRTkENRICH`](@ref) | ✓ (order-dep) | | ✓ (order-dep) | | -Note: the dofmap pattern describes the connection of the local degrees of freedom to entities of the grid and also hints to the continuity. Here, "N" or "n" means nodes, "F" or "f" means faces, "E" or "e" means edges and "I" means interior (dofs without any continuity across elements). Capital letters cause that every component has its own degree of freedom, while small letters signalize that only one dof is associated to the entity. As an example "N1f1" (for the Bernardi-Raugel element) means that at each node sits one dof per component and at each face sits a single dof. Usually finite elements that involve small letters are only defined vector-valued (i.e. the number of components has to match the element dimension), while finite elements that only involve capital letters are available for any number of components. +Note: The dofmap pattern describes how local degrees of freedom are associated with grid entities and provides insight into the continuity properties of the element. Here, "N" or "n" denotes nodes, "F" or "f" denotes faces, "E" or "e" denotes edges, and "I" denotes interior degrees of freedom (i.e., those without continuity across elements). Capital letters indicate that each component has its own degree of freedom, while lowercase letters mean only one degree of freedom is associated with the entity. For example, "N1f1" (as in the Bernardi-Raugel element) means that each node has one dof per component and each face has a single dof. Typically, finite elements involving lowercase letters are only defined for vector-valued cases (i.e., the number of components must match the element dimension), while those with only capital letters are available for any number of components. ## H1-conforming finite elements diff --git a/docs/src/fespace.md b/docs/src/fespace.md index d87db7c..a83cb33 100644 --- a/docs/src/fespace.md +++ b/docs/src/fespace.md @@ -1,7 +1,7 @@ # FESpace -To generate a finite element space only a finite element type and a grid is needed, dofmaps are generated automatically on their first demand. +A **finite element space** (FESpace) represents the set of all functions that can be expressed using a particular finite element type on a given computational grid. In ExtendableFEMBase, constructing a finite element space requires only specifying the finite element type and the grid; all necessary degree-of-freedom (dof) mappings are generated automatically on first access. ```@autodocs Modules = [ExtendableFEMBase] diff --git a/docs/src/fevector.md b/docs/src/fevector.md index 3d9f4b7..2b2631f 100644 --- a/docs/src/fevector.md +++ b/docs/src/fevector.md @@ -1,6 +1,6 @@ # FEVector -A FEVector consists of FEVectorBlocks that share a common one-dimensional array. Each block is associated to a FESpace and can only write into a region of the common array specified by offsets that stores the degrees of freedom of that FEspace. +A `FEVector` represents a block-structured vector used in finite element computations. It consists of one or more `FEVectorBlock`s, each associated with a specific `FESpace`. All blocks share a common one-dimensional array that stores the global degrees of freedom (DoFs). Each block can only write to a region of the array specified by offsets, ensuring that each `FESpace` manages its own DoFs within the global vector. ```@autodocs diff --git a/docs/src/functionoperators.md b/docs/src/functionoperators.md index ea74a84..be289bf 100644 --- a/docs/src/functionoperators.md +++ b/docs/src/functionoperators.md @@ -1,12 +1,10 @@ - # Function Operators ## StandardFunctionOperators -StandardFunctionOperators are abstract types that encode primitive (linear) operators (like Identity, Gradient etc.) -used to dispatch different evaluations of finite element basis functions. +`StandardFunctionOperators` are abstract types that represent fundamental (linear) operators, such as the identity, gradient, divergence, and others. These operators provide a unified interface for evaluating finite element basis functions in various ways. -### List of primitive operators +### List of Primitive Operators | StandardFunctionOperator | Description | Mathematically | | :--------------------------------------------------- | :------------------------------------------------------- | :------------------------------------------------------------------------ | @@ -28,12 +26,9 @@ used to dispatch different evaluations of finite element basis functions. !!! note - As each finite element type is transformed differently from the reference domain to the general domain, - the evaluation of each function operator has to be implemented for each finite element class. - Currently, not every function operator works in any dimension and for any finite element. More evaluations - are added as soon as they are needed (and possibly upon request). - Also, the function operators can be combined with user-defined actions to evaluate other operators that - can be build from the ones available (e.g. the deviator). + The transformation from the reference domain to the physical domain differs for each finite element class. As a result, the evaluation of each function operator must be implemented specifically for every finite element class. Not all function operators are currently available for every dimension or element type, but new implementations are added as needed or upon request. + + Additionally, function operators can be combined with user-defined kernels to postprocess/construct more advanced operators from the available primitives (for example, the deviatoric part of a tensor). ```@autodocs @@ -45,8 +40,7 @@ Order = [:type, :function] ## ReconstructionOperators -There are special operators that allow to evaluate a primitive operator of some discrete -reconstructed version of a testfunction. +Special operators are provided to evaluate a primitive operator on a reconstructed version of a test function. These are useful for advanced discretizations and post-processing. ```@autodocs Modules = [ExtendableFEMBase] @@ -54,16 +48,14 @@ Pages = ["reconstructionoperators.jl"] Order = [:type, :function] ``` -### Divergence-free reconstruction operators -For gradient-robust discretisations of certain classical non divergence-conforming ansatz spaces, -reconstruction operators are available that map a discretely divergence-free H1 function to a pointwise divergence-free -Hdiv function. So far such operators are available for the vector-valued Crouzeix-Raviart (H1CR) and Bernardi--Raugel (H1BR) finite element types, -as well as for the P2-bubble (H1P2B) finite element type in two dimensions. +### Divergence-Free Reconstruction Operators + +For gradient-robust discretizations of certain classical non-divergence-conforming ansatz spaces, reconstruction operators are available that map a discretely divergence-free H1 function to a pointwise divergence-free H(div) function. Currently, such operators are implemented for the vector-valued Crouzeix-Raviart (H1CR) and Bernardi–Raugel (H1BR) finite element types, as well as for the P2-bubble (H1P2B) element in two dimensions. -**Example:** Reconst{HDIVRT0{d}, Identity} gives the reconstruction of the Identity operator into HDIVRT0 (and is available for H1BR{d} and H1CR{d} for d = 1,2) +**Example:** `Reconst{HDIVRT0{d}, Identity}` reconstructs the identity operator into HDIVRT0, and is available for `H1BR{d}` and `H1CR{d}` for `d = 1, 2`. -## Operator Pairs (experimental) +## Operator Pairs (Experimental) -Two function operators can be put into an OperatorPair so that one can provide effectively two operators in each argument of an assembly pattern. However, the user should make sure that both operators can be evaluated together reasonably (meaning both should be well-defined on the element geometries and the finite element space where the argument will be evaluated, and the action of the operator has to operate with coressponding input and result fields). This feature is still experimental and might have issues in some cases. OperatorTriple for a combination of three operators is also available. +Two function operators can be combined into an `OperatorPair`, allowing you to provide two operators in each argument of an assembly pattern. Both operators must be well-defined on the relevant element geometries and finite element spaces, and their actions must be compatible with the input and result fields. This feature is experimental and may have limitations in some cases. An `OperatorTriple` is also available for combining three operators. diff --git a/docs/src/index.md b/docs/src/index.md index 4b0a72b..b9ecd85 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -5,7 +5,9 @@ # ExtendableFEMBase.jl -This package provides some low level structures like finite element spaces, interpolors, matrices and vectors to assemble custom finite element solvers based on [ExtendableGrids.jl](https://github.com/j-fu/ExtendableGrids.jl) infrastructure. +ExtendableFEMBase.jl provides foundational data structures and tools for building custom finite element solvers in Julia. The package includes flexible representations for finite element spaces, interpolators, matrices, and vectors, all designed to work seamlessly with the [ExtendableGrids.jl](https://github.com/j-fu/ExtendableGrids.jl) infrastructure. + +With ExtendableFEMBase.jl, users can efficiently assemble and manipulate finite element systems, enabling the development of advanced simulation workflows and research codes. ### Dependencies on other Julia packages diff --git a/docs/src/interpolations.md b/docs/src/interpolations.md index a8f7b53..e073c0d 100644 --- a/docs/src/interpolations.md +++ b/docs/src/interpolations.md @@ -1,4 +1,3 @@ - # Finite Element Interpolations ## Source functions and QPInfo @@ -22,13 +21,14 @@ The qpinfo argument communicates vast information of the current quadrature poin ## Standard Interpolations -Each finite element has its standard interpolator that can be applied to some user-defined source Function. Instead of interpolating on the full cells, the interpolation can be restricted to faces or edges via an AssemblyType. +Each finite element type provides a standard interpolation routine that can be applied to user-defined source functions. By default, interpolation is performed over all cells, but it can also be restricted to faces or edges using an appropriate `AssemblyType`. ```@docs interpolate! ``` -It is also possible to interpolate finite element functions on one grid onto a finite element function on another grid via the lazy_interpolate routine. +Additionally, you can transfer finite element functions from one grid to another using the `lazy_interpolate!` routine, which interpolates between different meshes. + ```@docs lazy_interpolate! @@ -43,7 +43,7 @@ continuify ## Nodal Evaluations -Usually, Plotters need nodal values, so there is a generic function that evaluates any finite element function at the nodes of the grids (possibly by averaging if discontinuous). In case of Identity evaluations of an H1-conforming finite element, the function nodevalues_view can generate a view into the coefficient field that avoids further allocations. +Plotting routines require nodal values, i.e., the values of a finite element function at the mesh nodes. The generic `nodevalues!` function evaluates any finite element function at the grid nodes, averaging values if the function is discontinuous. For H1-conforming finite elements and identity evaluations, the `nodevalues_view` function can provide a direct view into the coefficient field, avoiding unnecessary memory allocations. ```@docs diff --git a/docs/src/meshing.md b/docs/src/meshing.md index 06da415..c3c52e8 100644 --- a/docs/src/meshing.md +++ b/docs/src/meshing.md @@ -1,10 +1,10 @@ - # Meshing -Meshes are stored as an ExtendableGrid, see [ExtendableGrids.jl](https://github.com/j-fu/ExtendableGrids.jl) for details and constructors. -Grid generators for simplex grids can be found e.g. in the external module [SimplexGridFactory.jl](https://github.com/j-fu/SimplexGridFactory.jl) +Meshes in ExtendableFEMBase are represented using the `ExtendableGrid` type. For detailed information on grid structures and constructors, refer to the [ExtendableGrids.jl documentation](https://github.com/j-fu/ExtendableGrids.jl). + +To generate simplex grids, you can use external tools such as [SimplexGridFactory.jl](https://github.com/j-fu/SimplexGridFactory.jl). -Cells, faces and edges of the mesh are associated to AbstractElementGeometries (defined by [ExtendableGrids.jl](https://github.com/j-fu/ExtendableGrids.jl)) that are used to dispatch functionality (local/global transformation, enumeration rules, set of basis functions, volume calculation, refinements etc.). See further below for a list of recognized element geometries. +Each cell, face, and edge in a mesh is associated with an `AbstractElementGeometry` (as defined in [ExtendableGrids.jl](https://github.com/j-fu/ExtendableGrids.jl)). These geometries are used to dispatch key functionality, including local-to-global transformations, enumeration rules, basis function definitions, volume calculations, and mesh refinements. See below for a list of supported element geometries. ## Recognized Geometries and Reference Domains diff --git a/docs/src/plots.md b/docs/src/plots.md index 44055ff..97aefab 100644 --- a/docs/src/plots.md +++ b/docs/src/plots.md @@ -1,17 +1,15 @@ # Plots -## GridVisualize/PlutoVista +## GridVisualize and PlutoVista -Plotting is possible e.g. via [Nodal Evaluations](@ref) and the plot routines from +Visualization of finite element solutions is possible, for example, via [Nodal Evaluations](@ref) and the plotting routines from [ExtendableGrids.jl](https://github.com/WIAS-PDELib/ExtendableGrids.jl). -In Pluto notebooks it is recommended to use [PlutoVista.jl](https://github.com/j-fu/PlutoVista.jl) as an backend. - +For interactive Pluto notebooks, it is recommended to use [PlutoVista.jl](https://github.com/j-fu/PlutoVista.jl) as the backend for high-quality plots. ## UnicodePlots -For a fast and rough peak several UnicodePlots plotters are available via an extension (ExtendableFEMBaseUnicodePlotsExt) -that is loaded when UnicodePlots is available. - +For quick, in-terminal visualization, several UnicodePlots-based plotters are available via the `ExtendableFEMBaseUnicodePlotsExt` extension. +This extension is loaded automatically when UnicodePlots is available. ```@docs unicode_gridplot diff --git a/docs/src/pointevaluators.md b/docs/src/pointevaluators.md index 83efc0e..902c764 100644 --- a/docs/src/pointevaluators.md +++ b/docs/src/pointevaluators.md @@ -1,6 +1,6 @@ # PointEvaluator -Point evaluators allow to evaluate a finite element function (FEVector) at arbitrary points. +Point evaluators provide a convenient interface to evaluate finite element functions (`FEVector`) at arbitrary spatial points, not restricted to mesh nodes or quadrature points. This is useful for post-processing, visualization, or extracting solution values at specific locations. ```@autodocs Modules = [ExtendableFEMBase] diff --git a/docs/src/quadrature.md b/docs/src/quadrature.md index 30a54c8..751644c 100644 --- a/docs/src/quadrature.md +++ b/docs/src/quadrature.md @@ -1,12 +1,6 @@ - # Quadrature -Usually quadrature is a hidden layer as quadrature rules are chosen automatically based on the polynomial degree of the ansatz functions and the specified quadorder of the user data. - -Hence, quadrature rules are only needed if the user wants write his own low-level assembly. - - -Quadrature rules consist of points (coordinates of evaluation points with respect to reference geometry) and weights. There are constructors for several AbstractElementGeometries (from ExtendableGrids) and different order (some have generic formulas for arbitrary order), see below for a detailed list. +A quadrature rule consists of a set of points (coordinates of evaluation points in the reference geometry) and associated weights. Constructors are provided for various `AbstractElementGeometries` (from ExtendableGrids) and for different orders; some element types support generic formulas for arbitrary order. See below for a detailed list. ```@autodocs Modules = [ExtendableFEMBase] @@ -14,10 +8,9 @@ Pages = ["quadrature.jl"] Order = [:type, :function] ``` +#### Accumulating Vector (internal, for completeness) -#### Accumulating Vector (not relevant for users, but for completeness) - -Internally a global integration uses an accumulating vector and calls the cell-wise integration. +Internally, global integration uses an accumulating vector and calls cell-wise integration routines. ```@docs AccumulatingVector diff --git a/docs/src/segmentintegrators.md b/docs/src/segmentintegrators.md index 7298cf0..54bdb98 100644 --- a/docs/src/segmentintegrators.md +++ b/docs/src/segmentintegrators.md @@ -1,6 +1,6 @@ # SegmentIntegrator -Segment integrators allow to integrate a finite element function (FEVector) along arbitrary lines through mesh cells. +Segment integrators provide tools to integrate finite element functions (`FEVector`) along arbitrary line segments that pass through mesh cells. ```@autodocs Modules = [ExtendableFEMBase] diff --git a/examples/Example200_LowLevelPoisson.jl b/examples/Example200_LowLevelPoisson.jl index 18d08db..fb72582 100644 --- a/examples/Example200_LowLevelPoisson.jl +++ b/examples/Example200_LowLevelPoisson.jl @@ -30,20 +30,21 @@ using GridVisualize using UnicodePlots using Test # -## data for Poisson problem -const μ = 1.0 -const f = x -> x[1] - x[2] - -function main(; maxnref = 8, order = 2, Plotter = nothing) - +function main(; + maxnref = 8, + order = 2, + Plotter = nothing, + mu = 1.0, + rhs = x -> x[1] - x[2] + ) ## Finite element type FEType = H1Pk{1, 2, order} ## run once on a tiny mesh for compiling X = LinRange(0, 1, 4) xgrid = simplexgrid(X, X) - FES = FESpace{FEType}(xgrid) - sol, time_assembly, time_solve = solve_poisson_lowlevel(FES, μ, f) + fe_space = FESpace{FEType}(xgrid) + solution, time_assembly, time_solve = solve_poisson_lowlevel(fe_space, mu, rhs) ## loop over uniform refinements + timings plt = GridVisualizer(; Plotter = Plotter, layout = (1, 1), clear = true, resolution = (500, 500)) @@ -52,61 +53,60 @@ function main(; maxnref = 8, order = 2, Plotter = nothing) X = LinRange(0, 1, 2^level + 1) time_grid = @elapsed xgrid = simplexgrid(X, X) time_facenodes = @elapsed xgrid[FaceNodes] - FES = FESpace{FEType}(xgrid) - println("\nLEVEL = $level, ndofs = $(FES.ndofs)\n") + fe_space = FESpace{FEType}(xgrid) + println("\nLEVEL = $level, ndofs = $(fe_space.ndofs)\n") if level < 4 println(stdout, unicode_gridplot(xgrid)) end - time_dofmap = @elapsed FES[CellDofs] - sol, time_assembly, time_solve = solve_poisson_lowlevel(FES, μ, f) + time_dofmap = @elapsed fe_space[CellDofs] + solution, time_assembly, time_solve = solve_poisson_lowlevel(fe_space, mu, rhs) ## plot statistics println(stdout, barplot(["Grid", "FaceNodes", "celldofs", "Assembly", "Solve"], [time_grid, time_facenodes, time_dofmap, time_assembly, time_solve], title = "Runtimes")) ## plot if Plotter !== nothing - scalarplot!(plt[1, 1], xgrid, view(sol.entries, 1:num_nodes(xgrid)), limits = (-0.0125, 0.0125)) + scalarplot!(plt[1, 1], xgrid, view(solution.entries, 1:num_nodes(xgrid)), limits = (-0.0125, 0.0125)) else - sol_grad = continuify(sol[1], Gradient) - println(stdout, unicode_scalarplot(sol[1])) + solution_grad = continuify(solution[1], Gradient) + println(stdout, unicode_scalarplot(solution[1])) end end - return sol, plt + return solution, plt end -function solve_poisson_lowlevel(FES, μ, f) - Solution = FEVector(FES) - FES = Solution[1].FES - A = FEMatrix(FES, FES) - b = FEVector(FES) +function solve_poisson_lowlevel(fe_space, mu, rhs) + solution = FEVector(fe_space) + stiffness_matrix = FEMatrix(fe_space, fe_space) + rhs_vector = FEVector(fe_space) println("Assembling...") time_assembly = @elapsed @time begin - loop_allocations = assemble!(A.entries, b.entries, FES, f, μ) + loop_allocations = assemble!(stiffness_matrix.entries, rhs_vector.entries, fe_space, rhs, mu) ## fix boundary dofs begin - bdofs = boundarydofs(FES) + bdofs = boundarydofs(fe_space) for dof in bdofs - A.entries[dof, dof] = 1.0e60 - b.entries[dof] = 0 + stiffness_matrix.entries[dof, dof] = 1.0e60 + rhs_vector.entries[dof] = 0 end end - ExtendableSparse.flush!(A.entries) + ExtendableSparse.flush!(stiffness_matrix.entries) end ## solve println("Solving linear system...") - time_solve = @elapsed @time copyto!(Solution.entries, A.entries \ b.entries) + time_solve = @elapsed @time copyto!(solution.entries, stiffness_matrix.entries \ rhs_vector.entries) - return Solution, time_assembly, time_solve + return solution, time_assembly, time_solve end -function assemble!(A::ExtendableSparseMatrix, b::Vector, FES, f, μ = 1) - xgrid = FES.xgrid +function assemble!(A::ExtendableSparseMatrix, b::Vector, fe_space, rhs, mu = 1) + xgrid = fe_space.xgrid EG = xgrid[UniqueCellGeometries][1] - FEType = eltype(FES) + FEType = eltype(fe_space) L2G = L2GTransformer(EG, xgrid, ON_CELLS) ## quadrature formula @@ -117,11 +117,11 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES, f, μ = 1) cellvolumes = xgrid[CellVolumes] ## FE basis evaluator and dofmap - FEBasis_∇ = FEEvaluator(FES, Gradient, qf) + FEBasis_∇ = FEEvaluator(fe_space, Gradient, qf) ∇vals = FEBasis_∇.cvals - FEBasis_id = FEEvaluator(FES, Identity, qf) + FEBasis_id = FEEvaluator(fe_space, Identity, qf) idvals = FEBasis_id.cvals - celldofs = FES[CellDofs] + celldofs = fe_space[CellDofs] ## ASSEMBLY LOOP loop_allocations = 0 @@ -146,7 +146,7 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES, f, μ = 1) end Aloc[j, k] = temp end - Aloc .*= μ * cellvolumes[cell] + Aloc .*= mu * cellvolumes[cell] ## add local matrix to global matrix for j in 1:ndofs4cell @@ -173,7 +173,7 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES, f, μ = 1) ## get global x for quadrature point eval_trafo!(x, L2G, xref[qp]) ## evaluate (f(x), v_j(x)) - temp += weights[qp] * idvals[1, j, qp] * f(x) + temp += weights[qp] * idvals[1, j, qp] * rhs(x) end ## write into global vector dof_j = celldofs[j, cell] @@ -192,19 +192,24 @@ function generateplots(dir = pwd(); Plotter = nothing, kwargs...) return GridVisualize.save(joinpath(dir, "example200.png"), scene; Plotter = Plotter) end -function runtests(; order = 2, kwargs...) #hide +function runtests(; + order = 2, + mu = 1.0, + rhs = x -> x[1] - x[2], + kwargs... + ) FEType = H1Pk{1, 2, order} X = LinRange(0, 1, 64) xgrid = simplexgrid(X, X) - FES = FESpace{FEType}(xgrid) - A = FEMatrix(FES, FES) - b = FEVector(FES) - @info "ndofs = $(FES.ndofs)" + fe_space = FESpace{FEType}(xgrid) + stiffness_matrix = FEMatrix(fe_space, fe_space) + rhs_vector = FEVector(fe_space) + @info "ndofs = $(fe_space.ndofs)" ## first assembly causes allocations when filling sparse matrix - loop_allocations = assemble!(A.entries, b.entries, FES, f, μ) + loop_allocations = assemble!(stiffness_matrix.entries, rhs_vector.entries, fe_space, rhs, mu) @info "allocations in 1st assembly: $loop_allocations" ## second assembly in same matrix should have allocation-free inner loop - loop_allocations = assemble!(A.entries, b.entries, FES, f, μ) + loop_allocations = assemble!(stiffness_matrix.entries, rhs_vector.entries, fe_space, rhs, mu) @info "allocations in 2nd assembly: $loop_allocations" return @test loop_allocations == 0 end #hide diff --git a/src/ExtendableFEMBase.jl b/src/ExtendableFEMBase.jl index df3983c..f2f4e14 100644 --- a/src/ExtendableFEMBase.jl +++ b/src/ExtendableFEMBase.jl @@ -188,9 +188,36 @@ function unicode_gridplot( kwargs... ```` -(via extension that requires UnicodePlots) -Plots the grid on a UnicodePlots canvas (default: BrailleCanvas) by drawing all edges in the triangulation. +Render a quick terminal plot of a 2D grid using UnicodePlots, suitable for quick visualization in the terminal. + +This function draws the edges and boundary faces of the given `ExtendableGrid` using Unicode characters, with customizable resolution, colors, and plotting style. The plot can be based on cells, faces, or edges, and boundary faces are highlighted with a separate color. + +# Arguments +- `xgrid`: The `ExtendableGrid` to visualize. + +# Keyword Arguments +- `title`: Title of the plot (default: `"gridplot"`). +- `resolution`: Tuple specifying the number of columns and rows (characters) in the plot (default: `(40, 20)`). +- `autoscale`: Whether to automatically adjust the aspect ratio for the grid (default: `true`). +- `color`: RGB tuple for the main grid lines (default: `(200, 200, 200)`). +- `bface_color`: RGB tuple for boundary face lines (default: `(255, 0, 0)`). +- `CanvasType`: UnicodePlots canvas type to use (default: `BrailleCanvas`). +- `plot_based`: Plot based on `ON_CELLS`, `ON_FACES`, or `ON_EDGES` (default: `ON_CELLS`). +- `kwargs...`: Additional keyword arguments passed to the canvas or plot. + +# Returns +- A `UnicodePlots.Plot` object displaying the grid in the terminal. + +# Notes +- Requires the UnicodePlots extension to be loaded. + +# Example +```julia +using ExtendableFEMBase, ExtendableGrids, UnicodePlots +grid = simplexgrid(LinRange(0, 1, 5), LinRange(0, 1, 5)) +unicode_gridplot(grid; title = "My Grid", resolution = (60, 30)) +``` """ function unicode_gridplot end @@ -206,16 +233,27 @@ function unicode_scalarplot( kwargs...) ```` -(via extension that requires UnicodePlots) +Render a quick terminal plot of one or more components of a finite element function using UnicodePlots. + +This function visualizes the nodal values of the FE function stored in `u` by interpolating onto a coarse uniform mesh and plotting the result using UnicodePlots.jl. In 1D, all selected components are shown in a single line plot; in 2D, each component is shown as a separate heatmap. If `abs = true`, the Euclidean norm over all components is plotted instead. + +# Arguments +- `u`: The `FEVectorBlock` containing the finite element function to plot. -Plots components of the finite element function in the FEVectorBlock u by -using a lazy_interpolate! onto a coarse uniform mesh and UnicodePlots.jl -lineplot or heatmap for 1D or 2D, respectively. +# Keyword Arguments +- `components`: Indices of the components to plot (default: all components). +- `abs`: If `true`, plot the Euclidean norm over all components (default: `false`). +- `nrows`: Number of rows in the grid layout for 2D plots (default: `1`). +- `resolution`: Tuple specifying the plot resolution (characters) in each direction (default: `(30, 30)`). +- `colormap`: Colormap symbol for 2D heatmaps (default: `:viridis`). +- `title`: Title of the plot (default: `u.name`). +- `kwargs...`: Additional keyword arguments passed to the plotting routines. -In 1D all components all plotted in the same lineplot, while -in 2D all components are plotted in a separate heatmap. +# Returns +- A `UnicodePlots.Plot` object (or a vector of plots for multiple components in 2D). -If abs = true, only the absolute value over the components is plotted. +# Notes +- Requires the UnicodePlots extension to be loaded. """ function unicode_scalarplot end export unicode_gridplot, unicode_scalarplot diff --git a/src/dofmaps.jl b/src/dofmaps.jl index 4704bd2..a7f9417 100644 --- a/src/dofmaps.jl +++ b/src/dofmaps.jl @@ -1,9 +1,21 @@ """ $(TYPEDEF) -Dofmaps are stored as an ExtendableGrids.AbstractGridAdjacency in the finite element space and collect -information with respect to different AssemblyTypes. They are generated automatically on demand and the dofmaps -associated to each subtype can be accessed via FESpace[DofMap]. +A `DofMap` encodes the association of global DoF indices to mesh entities of a given type (e.g., cells, faces, edges) for a particular finite element space. This mapping is essential for assembling system matrices and vectors, applying boundary conditions, and extracting solution values. +DofMaps are typically not constructed directly by users. Instead, they are generated/managed internally by the finite element space and accessed as needed for assembly or evaluation tasks. + +# Implementation + +- All concrete `DofMap` subtypes (e.g., `CellDofs`, `FaceDofs`, `EdgeDofs`, etc.) specify the mesh entity type to which DoFs are attached. +- DofMaps are stored as `ExtendableGrids.AbstractGridAdjacency` (usually VariableTargetAdjacency, SerialVariableTargetAdjacency or AbstractGridIntegerArray2D) objects within the finite element space and are generated automatically as needed. +- The appropriate DofMap for a given assembly type can be accessed via `FESpace[DofMapSubtype]`. + +# Example + +```julia +cell_dofs = FESpace[CellDofs] # Get the cell-based DoF map +face_dofs = FESpace[FaceDofs] # Get the face-based DoF map +``` """ abstract type DofMap <: AbstractGridAdjacency end diff --git a/src/feevaluator.jl b/src/feevaluator.jl index a92a8fc..febbf89 100644 --- a/src/feevaluator.jl +++ b/src/feevaluator.jl @@ -33,12 +33,42 @@ Base.getindex(FEB::FEEvaluator, c, dof, qp) = FEB.cvals[c, dof, qp] function FEEvaluator(FE::FESpace, operator::AbstractFunctionOperator, qrule::QuadratureRule; T = Float64, AT = ON_CELLS, L2G = nothing) ```` -Constructs a FEEvaluator that handles evaluations of finite element basis function evaluation for the given FESpace, operator -at the quadrature points of the given QuadratureRule. It has an update! function to update the evaluation upon entry to a new cell. -Evaluations can be accessed via FEEvaluator.cvals[j,k,i] where i is the quadrature point id, k is the local dof number and j -is the component. +Construct a finite element basis function evaluator for a given finite element space, operator, and quadrature rule. -Note that matrix-valued operators evaluations, e.g. for Gradient, are given as a long vector (in component-wise order). +# Arguments +- `FE::FESpace`: The finite element space whose basis functions are to be evaluated. +- `operator::AbstractFunctionOperator`: The operator to apply to the basis functions (e.g., `Identity`, `Gradient`, etc.). +- `qrule::QuadratureRule`: The quadrature rule specifying the evaluation points (quadrature points) on the reference element. +- `xgrid`: (optional, defaults to `FE.dofgrid`) The grid on which the FE space is defined. Used for geometric information. +- `L2G`: (optional) A local-to-global transformation object. If not provided, it is constructed automatically. +- `T`: (optional, default `Float64`) The floating-point type for computations. +- `AT`: (optional, default `ON_CELLS`) The assembly type, specifying the geometric entity (cells, faces, etc.) for evaluation. + +# Returns +A `FEEvaluator` object that can efficiently evaluate (and cache) the values of the specified operator applied to the FE basis functions at the quadrature points of each cell (or other entity). The evaluator supports fast updates to new cells via `update_basis!`. + +# Usage + +- After construction, call `update_basis!(FEB, cellid)` to update the evaluator to a new cell. +- Access the evaluated values via `FEB.cvals[component, dof, qp]`, where: + - `component`: the output component (e.g., vector or matrix entry) + - `dof`: the local basis function index + - `qp`: the quadrature point index + +# Notes + +- For matrix-valued operators (e.g., `Gradient`), the result is stored as a long vector in component-wise order. +- The evaluator is designed for high performance in assembly loops and supports both scalar and vector-valued elements and operators. +- For advanced use, the evaluator exposes internal fields for basis values, derivatives, and transformation matrices. + +# Example + +```julia +FEB = FEEvaluator(FE, Gradient, qrule) +for cell in 1:ncells + update_basis!(FEB, cell) + # Access FEB.cvals for basis gradients at quadrature points +end """ function FEEvaluator( FE::FESpace{TvG, TiG, FEType, FEAPT}, diff --git a/src/fematrix.jl b/src/fematrix.jl index ac9993f..a986ced 100644 --- a/src/fematrix.jl +++ b/src/fematrix.jl @@ -140,76 +140,70 @@ $(TYPEDSIGNATURES) Custom `show` function for `FEMatrix` that prints some information on its blocks. """ function Base.show(io::IO, FEM::FEMatrix{TvM, TiM, TvG, TiG, nbrow, nbcol, nbtotal}) where {TvM, TiM, TvG, TiG, nbrow, nbcol, nbtotal} + if length(FEM) == 0 + println(io, "FEMatrix is empty.") + return + end println(io, "\nFEMatrix information") println(io, "====================") - println(io, " block | starts | ends | size | name ") + @printf(io, " block | starts (row,col) | ends (row,col) | size (row,col) | name\n") for j in 1:length(FEM) n = mod(j - 1, nbrow) + 1 m = Int(ceil(j / nbrow)) - @printf(io, " [%2d,%2d] |", m, n) - @printf(io, " (%6d,%6d) |", FEM[j].offset + 1, FEM[j].offsetY + 1) - @printf(io, " (%6d,%6d) |", FEM[j].last_index, FEM[j].last_indexY) - @printf(io, " (%6d,%6d) |", FEM[j].FES.ndofs, FEM[j].FESY.ndofs) - @printf(io, " %s \n", FEM[j].name) + @printf( + io, " [%2d,%2d] | (%6d,%6d) | (%6d,%6d) | (%6d,%6d) | %s\n", + m, n, + FEM[j].offset + 1, FEM[j].offsetY + 1, + FEM[j].last_index, FEM[j].last_indexY, + FEM[j].FES.ndofs, FEM[j].FESY.ndofs, + FEM[j].name + ) end - return println(io, "\n nnzvals = $(length(FEM.entries.cscmatrix.nzval))") + println(io, "\n total size = $(size(FEM.entries, 1)) x $(size(FEM.entries, 2))") + println(io, " nnzvals = $(length(FEM.entries.cscmatrix.nzval))") + return end """ ```` -FEMatrix{TvM,TiM}(name::String, FES::FESpace{TvG,TiG,FETypeX,APTX}) where {TvG,TiG,FETypeX,APTX} +FEMatrix(FESX::Union{FESpace, Vector{<:FESpace}}, FESY::Union{FESpace, Vector{<:FESpace}}=FESX; TvM=Float64, TiM=Int64, ...) ```` -Creates FEMatrix with one square block (FES,FES). -""" -function FEMatrix(FES::FESpace; name = :automatic) - return FEMatrix{Float64, Int64}(FES; name = name) -end -function FEMatrix{TvM}(FES::FESpace; name = :automatic) where {TvM} - return FEMatrix{TvM, Int64}(FES; name = name) -end -function FEMatrix{TvM, TiM}(FES::FESpace{TvG, TiG, FETypeX, APTX}; name = :automatic) where {TvM, TiM, TvG, TiG, FETypeX, APTX} - return FEMatrix{TvM, TiM}([FES], [FES]; name = name) -end +Constructs an `FEMatrix` for storing (sparse) matrix representations associated with one or more pairs of finite element spaces (`FESpace`). -""" -```` -FEMatrix{TvM,TiM}(FESX, FESY; name = :automatic) -```` +-- If `FESX` and `FESY` are single `FESpace` objects, the resulting `FEMatrix` contains one rectangular block. +- If `FESX` and/or `FESY` are vectors of `FESpace` objects, the resulting `FEMatrix` is block-structured, with one block for each pair. -Creates FEMatrix with one rectangular block (FESX,FESY) if FESX and FESY are single FESpaces, or -a rectangular block matrix with blocks corresponding to the entries of the FESpace vectors FESX and FESY. -Optionally a name for the matrix can be given. -""" -function FEMatrix(FESX::FESpace, FESY::FESpace; kwargs...) - return FEMatrix{Float64, Int64}(FESX, FESY; kwargs...) -end -function FEMatrix(FESX::Vector{<:FESpace}, FESY::Vector{<:FESpace}; kwargs...) - return FEMatrix{Float64, Int64}(FESX, FESY; kwargs...) -end -function FEMatrix{TvM}(FESX::FESpace, FESY::FESpace; kwargs...) where {TvM} - return FEMatrix{TvM, Int64}(FESX, FESY; kwargs...) -end -function FEMatrix{TvM, TiM}(FESX::FESpace, FESY::FESpace; kwargs...) where {TvM, TiM} - return FEMatrix{TvM, TiM}([FESX], [FESY]; kwargs...) -end -function FEMatrix(FES::Array{<:FESpace{TvG, TiG}, 1}; kwargs...) where {TvG, TiG} - return FEMatrix{Float64, Int64}(FES; kwargs...) -end -function FEMatrix{TvM}(FES::Array{<:FESpace{TvG, TiG}, 1}; kwargs...) where {TvM, TvG, TiG} - return FEMatrix{TvM, Int64}(FES, FES; kwargs...) -end -function FEMatrix{TvM, TiM}(FES::Array{<:FESpace{TvG, TiG}, 1}; kwargs...) where {TvM, TiM, TvG, TiG} - return FEMatrix{TvM, TiM}(FES, FES; kwargs...) -end +Optionally, you can assign a name (as a `String` for all blocks) and/or tags (as arrays for rows and columns) to the blocks for identification and access. -""" -```` -FEMatrix{TvM,TiM}(FESX, FESY; name = :automatic) -```` +# Arguments +- `FESX::FESpace` or `FESX::Vector{<:FESpace}`: Row finite element space(s). +- `FESY::FESpace` or `FESY::Vector{<:FESpace}`: Column finite element space(s) (default: same as FESX). + +# Keyword Arguments +- `entries`: Optional sparse matrix of coefficients. If not provided, a new sparse matrix of appropriate size is created. +- `name`: Name for the matrix or for each block (default: `:automatic`). +- `tags`: Array of tags for both rows and columns (default: `nothing`). +- `tagsX`: Array of tags for the row blocks (default: `tags`). +- `tagsY`: Array of tags for the column blocks (default: `tagsX`). +- `npartitions`: Number of partitions for the underlying sparse matrix (default: `1`). +- Additional keyword arguments are passed to the underlying block constructors. + +# Returns +- An `FEMatrix` object with one or more `FEMatrixBlock`s, each corresponding to a given pair of `FESpace` objects. -Creates an FEMatrix with blocks coressponding to the ndofs of FESX (rows) and FESY (columns). """ +function FEMatrix( + FESX::Union{FESpace, Vector{<:FESpace}}, + FESY::Union{FESpace, Vector{<:FESpace}} = FESX; + TvM = Float64, TiM = Int64, kwargs... + ) + # Convert single FESpace to vector + FESXv = isa(FESX, FESpace) ? [FESX] : FESX + FESYv = isa(FESY, FESpace) ? [FESY] : FESY + return FEMatrix{TvM, TiM}(FESXv, FESYv; kwargs...) +end + function FEMatrix{TvM, TiM}(FESX::Array{<:FESpace{TvG, TiG}, 1}, FESY::Array{<:FESpace{TvG, TiG}, 1}; entries = nothing, name = :automatic, tags = nothing, tagsX = tags, tagsY = tagsX, npartitions = 1, kwargs...) where {TvM, TiM, TvG, TiG} ndofsX, ndofsY = 0, 0 for j in 1:length(FESX) diff --git a/src/fevector.jl b/src/fevector.jl index bb51b30..c53e436 100644 --- a/src/fevector.jl +++ b/src/fevector.jl @@ -8,7 +8,26 @@ """ $(TYPEDEF) -block of an FEVector that carries coefficients for an associated FESpace and can be assigned as an AbstractArray (getindex, setindex, size, length) +A block within an `FEVector` representing a contiguous segment of coefficients associated with a specific finite element space (`FESpace`). + +Each `FEVectorBlock` provides array-like access to the degrees of freedom (DOFs) for its associated `FESpace`, mapping local indices to a shared global coefficient array. This enables efficient block-wise operations, assembly, and extraction of sub-vectors corresponding to different FE spaces. + +# Type Parameters +- `T`: Value type of the vector entries (e.g., `Float64`). +- `Tv`: Value type for the associated `FESpace`. +- `Ti`: Integer type for the associated `FESpace`. +- `FEType`: Type of the finite element. +- `APT`: Assembly type for the finite element. + +# Fields +- `name::String`: Name or label for this block (for identification or debugging). +- `FES::FESpace{Tv, Ti, FEType, APT}`: The finite element space associated with this block. +- `offset::Int`: Global offset (start index in the global vector). +- `last_index::Int`: Global end index (inclusive). +- `entries::Array{T, 1}`: Reference to the global coefficient array (shared with the parent `FEVector`). + +# Usage +`FEVectorBlock` is typically created internally by `FEVector` constructors and provides efficient access to the coefficients for a particular FE space. Supports standard array operations (`getindex`, `setindex!`, `size`, `length`, etc.) and can be used for block-wise assembly, extraction, and manipulation. """ struct FEVectorBlock{T, Tv, Ti, FEType, APT} <: AbstractArray{T, 1} name::String @@ -42,9 +61,20 @@ end """ $(TYPEDEF) -a plain array but with an additional layer of several FEVectorBlock subdivisions each carrying coefficients for their associated FESpace. -The j-th block can be accessed by getindex(::FEVector, j) or by getindex(::FEVector, tag) if tags are associated. -The full vector can be accessed via FEVector.entries +A block-structured vector for storing coefficients associated with one or more finite element spaces (`FESpace`). + +An `FEVector` consists of a global coefficient array subdivided into multiple `FEVectorBlock`s, each corresponding to a specific `FESpace`. This structure enables efficient block-wise access, assembly, and manipulation of solution vectors in finite element computations, especially for multi-field or mixed problems. + +# Type Parameters +- `T`: Value type of the vector entries (e.g., `Float64`). +- `Tv`: Value type for the associated `FESpace`. +- `Ti`: Integer type for the associated `FESpace`. + +# Fields +- `FEVectorBlocks::Array{FEVectorBlock{T, Tv, Ti}, 1}`: Array of blocks, each representing a segment of the global vector for a specific `FESpace`. +- `entries::Array{T, 1}`: The global coefficient array, shared by all blocks. +- `tags::Vector{Any}`: Optional tags for identifying or accessing blocks (e.g., by name or symbol). + """ struct FEVector{T, Tv, Ti} #<: AbstractVector{T} FEVectorBlocks::Array{FEVectorBlock{T, Tv, Ti}, 1} @@ -135,9 +165,24 @@ Base.length(FEB::FEVectorBlock) = FEB.last_index - FEB.offset FEVector{T}(FES; name = nothing, tags = nothing, kwargs...) where T <: Real ```` -Creates FEVector that has one block if FES is a single FESpace, and a blockwise FEVector if FES is a vector of FESpaces. -Optionally a name for the vector (as a String) or each of the blocks (as a vector of Strings), or tags (as an Array{Any}) -for the blocks can be specified. +Constructs an `FEVector` for storing coefficients associated with one or more finite element spaces (`FESpace`). + +- If `FES` is a single `FESpace`, the resulting `FEVector` contains one block. +- If `FES` is a vector of `FESpace` objects, the resulting `FEVector` is block-structured, with one block per space. + +Optionally, you can assign a name (as a `String` for all blocks, or a vector of `String` for each block) and/or tags (as an array of any type) to the blocks for identification and access. + +# Arguments +- `FES::FESpace` or `FES::Vector{<:FESpace}`: The finite element space(s) for which to create the vector. + +# Keyword Arguments +- `entries`: Optional array of coefficients. If not provided, a zero vector of appropriate length is created. +- `name`: Name for the vector or for each block (default: `nothing` causes auto naming by index or tag). +- `tags`: Array of tags for the blocks (default: `[]`, i.e. block access only by index). + +# Returns +- An `FEVector` object with one or more `FEVectorBlock`s, each corresponding to a given `FESpace`. + """ function FEVector(FES::FESpace{Tv, Ti, FEType, APT}; kwargs...) where {Tv, Ti, FEType, APT} return FEVector{Float64}([FES]; kwargs...) @@ -192,23 +237,25 @@ $(TYPEDSIGNATURES) Custom `show` function for `FEVector` that prints some information on its blocks. """ function Base.show(io::IO, FEF::FEVector) + if length(FEF) == 0 + println(io, "FEVector is empty.") + return + end println(io, "\nFEVector information") println(io, "====================") - print(io, " block | starts | ends | length | min / max \t| FEType \t\t (name/tag)") + @printf(io, " block | starts | ends | length | min / max | FEType | (name/tag)\n") for j in 1:length(FEF) - @printf(io, "\n [%5d] |", j) - @printf(io, " %6d |", FEF[j].offset + 1) - @printf(io, " %6d |", FEF[j].last_index) - @printf(io, " %6d |", FEF[j].FES.ndofs) ext = extrema(view(FEF[j])) - @printf(io, " %.2e/%.2e \t|", ext[1], ext[2]) - if length(FEF.tags) >= j - @printf(io, " %s \t (%s)", FEF[j].FES.name, FEF.tags[j]) - else - @printf(io, " %s \t (%s)", FEF[j].FES.name, FEF[j].name) - end + name = FEF[j].FES.name + tag = length(FEF.tags) >= j ? FEF.tags[j] : FEF[j].name + @printf( + io, " [%5d] | %6d | %6d | %6d | %12.4e / %12.4e | %-16s | %s\n", + j, FEF[j].offset + 1, FEF[j].last_index, FEF[j].FES.ndofs, ext[1], ext[2], name, tag + ) end - return + total_dofs = sum(FEF[j].FES.ndofs for j in 1:length(FEF)) + println(io, "\n total size = $total_dofs") + return nothing end diff --git a/src/finiteelements.jl b/src/finiteelements.jl index b616ac7..74ef9cd 100644 --- a/src/finiteelements.jl +++ b/src/finiteelements.jl @@ -57,18 +57,36 @@ abstract type AbstractInterpolationOperator end """ ```` struct FESpace{Tv, Ti, FEType<:AbstractFiniteElement,AT<:AssemblyType} - name::String # full name of finite element space (used in messages) - broken::Bool # if true, broken dofmaps are generated - ndofs::Int # total number of dofs - coffset::Int # offset for component dofs - xgrid::ExtendableGrid[Tv,Ti} # link to xgrid - dofgrid::ExtendableGrid{Tv,Ti} # link to (sub) grid used for dof numbering (expected to be equal to or child grid of xgrid) - dofmaps::Dict{Type{<:AbstractGridComponent},Any} # backpack with dofmaps - interpolators::Dict{Type{<:AssemblyType}, Any} # backpack with interpolators + name::String + broken::Bool + ndofs::Int + coffset::Int + xgrid::ExtendableGrid[Tv,Ti} + dofgrid::ExtendableGrid{Tv,Ti} + dofmaps::Dict{Type{<:AbstractGridComponent},Any} + interpolators::Dict{Type{<:AssemblyType}, Any} end ```` +A finite element space representing the global collection of degrees of freedom (DOFs) for a given finite element type and assembly type on a computational grid. + +`FESpace` encapsulates the mapping between mesh entities (cells, faces, edges, etc.) and DOFs, as well as metadata and auxiliary structures needed for assembly, interpolation, and evaluation of finite element functions. + +# Type Parameters +- `Tv`: Value type for grid coordinates (e.g., `Float64`). +- `Ti`: Integer type for grid indices (e.g., `Int64`). +- `FEType`: The finite element type (e.g., `H1P1`). +- `AT`: The assembly type (e.g., `ON_CELLS`). + +# Fields +- `name::String`: Name of the finite element space. +- `broken::Bool`: Whether the space is "broken" (discontinuous across elements). +- `ndofs::Int64`: Total number of global degrees of freedom. +- `coffset::Int`: Offset for component DOFs (for vector-valued or mixed spaces). +- `xgrid::ExtendableGrid{Tv, Ti}`: Reference to the computational grid. +- `dofgrid::ExtendableGrid{Tv, Ti}`: Grid used for DOF numbering (may be a subgrid of `xgrid`). +- `dofmaps::Dict{Type{<:AbstractGridComponent}, Any}`: Dictionary of DOF maps for different grid components (cells, faces, edges, etc.). +- `interpolators::Dict{Type{<:AssemblyType}, Any}`: Dictionary of interpolation operators for different assembly types. -A struct that has a finite element type as parameter and carries dofmaps (CellDofs, FaceDofs, BFaceDofs) plus additional grid information and access to arrays holding coefficients if needed. """ struct FESpace{Tv, Ti, FEType <: AbstractFiniteElement, AT <: AssemblyType} name::String # full name of finite element space (used in messages) @@ -88,6 +106,27 @@ end """ $(TYPEDSIGNATURES) +returns the name of the finite element space. +""" +name(FES::FESpace) = FES.name + +""" +$(TYPEDSIGNATURES) + +returns the computational grid of the finite element space. +""" +xgrid(FES::FESpace) = FES.xgrid + +""" +$(TYPEDSIGNATURES) + +returns the dofgrid of the finite element space. +""" +dofgrid(FES::FESpace) = FES.dofgrid + +""" +$(TYPEDSIGNATURES) + returns the total number of degrees of freedom of the finite element space. """ ndofs(FES::FESpace) = FES.ndofs @@ -132,15 +171,24 @@ Base.setindex!(FES::FESpace, v, DM::Type{<:DofMap}) = FES.dofmaps[DM] = v """ ```` -function FESpace{FEType<:AbstractFiniteElement,AT<:AssemblyType}( - xgrid::ExtendableGrid{Tv,Ti}; - name = "", - broken::Bool = false) +FESpace{FEType, AT}(xgrid::ExtendableGrid{Tv, Ti}; name = "", regions = nothing, broken::Bool = false) ```` -Constructor for FESpace of the given FEType, AT = ON_CELLS/ON_FACES/ON_EDGES generates a finite elements space on the cells/faces/edges of the provided xgrid (if omitted ON_CELLS is used as default). -The broken switch allows to generate a broken finite element space (that is piecewise H1/Hdiv/HCurl). If no name is provided it is generated automatically from FEType. -If no AT is provided, the space is generated ON_CELLS. +Constructs a finite element space (`FESpace`) of the given finite element type (`FEType`) and assembly type (`AT`) on the provided computational grid. + +- The `FESpace` represents the global collection of degrees of freedom (DOFs) for the specified finite element type and assembly type, and manages the mapping between mesh entities (cells, faces, edges, etc.) and DOFs. +- The `broken` switch allows the creation of a "broken" (discontinuous) finite element space, where basis functions are not required to be continuous across element boundaries. +- The `regions` argument can be used to restrict the space to a subset of the grid. +- If no `AT` is provided, the space is generated on cells (`ON_CELLS`). + +# Arguments +- `xgrid::ExtendableGrid{Tv, Ti}`: The computational grid on which the finite element space is defined. + +# Keyword Arguments +- `name::String: Name for the finite element space (for identification and debugging) (default: "", triggers automatic generation from FEType and broken). +- `regions`: Optional subset of the grid to restrict the space (default: all regions). +- `broken`: Whether to create a broken (discontinuous) space (default: false). + """ function FESpace{FEType, AT}( xgrid::ExtendableGrid{Tv, Ti}; @@ -226,14 +274,20 @@ function Base.show(io::IO, FES::FESpace{Tv, Ti, FEType, APT}) where {Tv, Ti, FET println(io, " name = $(FES.name)") println(io, " FEType = $FEType ($(FES.broken ? "$APT, broken" : "$APT"))") println(io, " FEClass = $(supertype(FEType))") - println(io, " ndofs = $(FES.ndofs) $(FES.coffset !== FES.ndofs ? "(coffset = $(FES.coffset))" : "")\n") + println(io, " ndofs = $(FES.ndofs) $(FES.coffset !== FES.ndofs ? "(coffset = $(FES.coffset))" : "")") println(io, " xgrid = $(FES.xgrid)") println(io, " dofgrid = $(FES.dofgrid !== FES.xgrid ? FES.dofgrid : "xgrid")") - println(io, "") - println(io, "DofMaps") - println(io, "=======") - for tuple in FES.dofmaps - println(io, "> $(tuple[1])") + if !isempty(FES.dofmaps) + println(io, "\nDofMaps:") + for (k, v) in FES.dofmaps + println(io, " > $(k): $(typeof(v))") + end + end + if !isempty(FES.interpolators) + println(io, "\nInterpolators:") + for (k, v) in FES.interpolators + println(io, " > $(k): $(typeof(v))") + end end return end diff --git a/src/interpolations.jl b/src/interpolations.jl index d76b980..71cac3c 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -14,26 +14,31 @@ end """ ```` -function ExtendableGrids.interpolate!(target::FEVectorBlock, - AT::Type{<:AssemblyType}, - source!::Function; - items = [], - bonus_quadorder = 0, - time = 0, - kwargs...) +function ExtendableGrids.interpolate!( + target::FEVectorBlock{T, Tv, Ti}, + AT::Type{<:AssemblyType}, + source; + items = [], + kwargs... +) where {T, Tv, Ti} ```` -Interpolates the given source into the finite elements space assigned to the target FEVectorBlock with the specified AssemblyType -(usually ON_CELLS). +Interpolate a function or data into the finite element space (i.e. computes the coefficients) associated with `target`, using the specified assembly type. + +# Arguments +- `target::FEVectorBlock`: The block of the FE vector to store the interpolated coefficients. +- `AT::Type{<:AssemblyType}`: The assembly type specifying where interpolation is performed (e.g., `ON_CELLS`). +- `source`: The function or callable to interpolate. Should have the signature `source!(result, qpinfo)`. -The source functions should adhere to the interface -```julia - source!(result, qpinfo) -``` -The qpinfo argument communicates vast information of the current quadrature/evaluation point. +# Keyword Arguments +- `items`: List of mesh entities (cells, faces, etc.) to interpolate on. If empty, all entities of the specified type are used. +- `kwargs...`: Additional keyword arguments passed to lower-level routines (e.g., `bonus_quadorder`, `time`). + +# Notes +- For "broken" FE spaces, interpolation is performed in a continuous auxiliary space and then mapped to the broken space. +- The `source!` function is called at each quadrature point and should fill `result` with the function values at that point. +- The `qpinfo` argument provides information about the current quadrature point, including coordinates, weights, and possibly time. -The bonus_quadorder argument can be used to steer the quadrature order of integrals that needs to be computed -for the interpolation (the default quadrature order corresponds to the polynomial order of the finite element). """ function ExtendableGrids.interpolate!( target::FEVectorBlock{T, Tv, Ti}, @@ -77,22 +82,10 @@ end ```` function ExtendableGrids.interpolate!(target::FEVectorBlock, source::Function; - items = [], - bonus_quadorder = 0, - time = 0, kwargs...) ```` -Interpolates the given source function into the finite element space assigned to the target FEVectorBlock. - -The source functions should adhere to the interface -```julia - source!(result, qpinfo) -``` -The qpinfo argument communicates vast information of the current quadrature/evaluation point. - -The bonus_quadorder argument can be used to steer the quadrature order of integrals that needs to be computed -for the interpolation (the default quadrature order corresponds to the polynomial order of the finite element). +see interpolate!(target, ON_CELLS, source; kwargs...) """ function ExtendableGrids.interpolate!(target::FEVectorBlock, source; kwargs...) return interpolate!(target, ON_CELLS, source; kwargs...) @@ -115,11 +108,26 @@ function nodevalues_subset!( continuous::Bool = false) ```` -Evaluates the finite element function with the coefficient vector source (interpreted as a coefficient vector for the FESpace FE) -and the specified FunctionOperator at the specified list of nodes of the grid and writes the values in that order into target. -Node values for nodes that are not part of the specified regions (default = all regions) are set to zero. -Discontinuous (continuous = false) quantities are evaluated in all neighbouring cells (in the specified regions) -of each node and then averaged. Continuous (continuous = true) quantities are only evaluated once at each node. +Evaluate (an operator of) a finite element function (given by the coefficient vector `source` and FE space `FE`) at a specified subset of nodes, and write the results into `target`. + +For each node in `nodes`, the function is evaluated (optionally with `operator`) and the result is written to the corresponding column of `target`. If `continuous` is `false`, values are averaged over all neighboring cells; if `true`, only one cell is used per node. If `abs` is `true`, the Euclidean norm is computed instead of the raw values. + +# Arguments +- `target`: Output array to store the evaluated values (size: result dimension × number of nodes). +- `source`: Coefficient vector for the FE function. +- `FE`: The finite element space. +- `operator`: Function operator to apply (default: `Identity`). +- `abs`: If `true`, store the Euclidean norm of the result at each node (default: `false`). +- `factor`: Scaling factor applied to the result (default: `1`). +- `nodes`: List of node indices to evaluate (default: all nodes). +- `regions`: List of region indices to restrict evaluation (default: all regions). +- `target_offset`: Offset for writing into `target` (default: `0`). +- `source_offset`: Offset for reading from `source` (default: `0`). +- `zero_target`: If `true`, zero out `target` before writing (default: `true`). +- `continuous`: If `true`, evaluate only once per node; otherwise, average over all neighboring cells (default: `false`). + +# Notes +- The result dimension is determined by the FE space and operator and the `abs` argument. """ function nodevalues_subset!( target::AbstractArray{T, 2}, @@ -129,7 +137,7 @@ function nodevalues_subset!( abs::Bool = false, factor = 1, nodes = [], - regions::Array{Int, 1} = [0], # are ignored at the moment + regions::Array{Int, 1} = [0], target_offset::Int = 0, source_offset::Int = 0, zero_target::Bool = true, @@ -256,20 +264,27 @@ function nodevalues!( source::AbstractArray{T,1}, FE::FESpace{Tv,Ti,FEType,AT}, operator::Type{<:AbstractFunctionOperator} = Identity; - regions::Array{Int,1} = [0], - abs::Bool = false, - factor = 1, - target_offset::Int = 0, # start to write into target after offset - zero_target::Bool = true, # target vector is zeroed - continuous::Bool = false) + kwargs...) ```` -Evaluates the finite element function with the coefficient vector source (interpreted as a coefficient vector for the FESpace FE) -and the specified FunctionOperator at all the nodes of the (specified regions of the) grid and writes the values into target. -Discontinuous (continuous = false) quantities are evaluated in all neighbouring cells of each node and then averaged. Continuous -(continuous = true) quantities are only evaluated once at each node. +calls the nodevalues_subset! function with all nodes of the dofgrid, see nodevalues_subset!(target, source, FE, operator; nodes = 1:num_nodes(FE.dofgrid), kwargs...) """ function nodevalues!( + target::AbstractArray{T, 2}, + source::AbstractArray{T, 1}, + FE::FESpace{Tv, Ti, FEType, AT}, + operator::Type{<:AbstractFunctionOperator} = Identity; + kwargs... + ) where {T, Tv, Ti, FEType, AT} + + nodevalues_subset!(target, source, FE, operator; nodes = 1:num_nodes(FE.dofgrid), kwargs...) + + return nothing +end + +""" +```` +piecewise_nodevalues!( target::AbstractArray{T, 2}, source::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, AT}, @@ -281,6 +296,38 @@ function nodevalues!( source_offset::Int = 0, zero_target::Bool = true, continuous::Bool = false + ) +```` + +Evaluate a finite element function (given by the coefficient vector `source` and FE space `FE`) at all nodes, but store the results in a piecewise (cellwise) fashion, i.e., for each cell, the values at its local nodes are written into `target`. The result is organized so that each column of `target` corresponds to a cell, and each row corresponds to the values at the cell's nodes. + +# Arguments +- `target`: Output array to store the evaluated values (size: result dimension × number of nodes per cell, number of cells). +- `source`: Coefficient vector for the FE function. +- `FE`: The finite element space. +- `operator`: Function operator to apply (default: `Identity`). +- `abs`: If `true`, store the Euclidean norm of the result at each node (default: `false`). +- `factor`: Scaling factor applied to the result (default: `1`). +- `regions`: List of region indices to restrict evaluation (default: all regions). +- `target_offset`: Offset for writing into `target` (default: `0`). +- `source_offset`: Offset for reading from `source` (default: `0`). +- `zero_target`: If `true`, zero out `target` before writing (default: `true`). +- `continuous`: If `true`, evaluate only once per node; otherwise, average over all neighboring cells (default: `false`). + +# Notes +- The result dimension is determined by the FE space, the operator, and the `abs` argument. +""" +function piecewise_nodevalues!( + target::AbstractArray{T, 2}, + source::AbstractArray{T, 1}, + FE::FESpace{Tv, Ti, FEType, AT}, + operator::Type{<:AbstractFunctionOperator} = Identity; + abs::Bool = false, + factor = 1, + regions::Array{Int, 1} = [0], + target_offset::Int = 0, + source_offset::Int = 0, + zero_target::Bool = true ) where {T, Tv, Ti, FEType, AT} xgrid = FE.dofgrid @@ -307,14 +354,11 @@ function nodevalues!( fill!(target, 0) end - nnodes::Int = num_sources(xgrid[Coordinates]) - nneighbours::Array{Int, 1} = zeros(Int, nnodes) - flag4node::Array{Bool, 1} = zeros(Bool, nnodes) - function barrier(EG, qf, BE) node::Int = 0 dof::Ti = 0 weights::Array{T, 1} = qf.w + nweights = length(weights) basisvals = BE.cvals ndofs = size(basisvals, 2) cvals_resultdim::Int = size(basisvals, 1) @@ -334,27 +378,30 @@ function nodevalues!( for i in eachindex(weights) # vertices node = xItemNodes[i, item] fill!(localT, 0) - if continuous == false || flag4node[node] == false - nneighbours[node] += 1 - flag4node[node] = true - for dof_i in 1:ndofs - dof = xItemDofs[dof_i, item] - eval_febe!(temp, BE, dof_i, i) - for k in 1:cvals_resultdim - localT[k] += source[source_offset + dof] * temp[k] - #target[k+target_offset,node] += temp[k] * source[source_offset + dof] - end + + for dof_i in 1:ndofs + dof = xItemDofs[dof_i, item] + eval_febe!(temp, BE, dof_i, i) + for k in 1:cvals_resultdim + localT[k] += source[source_offset + dof] * temp[k] + #target[k+target_offset,node] += temp[k] * source[source_offset + dof] end - localT .*= factor - if abs - for k in 1:cvals_resultdim - target[1 + target_offset, node] += localT[k]^2 - end - else - for k in 1:cvals_resultdim - target[k + target_offset, node] += localT[k] - end + end + localT .*= factor + if abs + for k in 1:cvals_resultdim + target[target_offset + i, item] += localT[k]^2 end + else + for k in 1:cvals_resultdim + target[target_offset + i + (k - 1) * nweights, item] += localT[k] + end + end + end + + if abs + for i in 1:nweights + target[i + target_offset, item] = sqrt(target[i + target_offset, item]) end end break # region for loop @@ -370,26 +417,53 @@ function nodevalues!( barrier(EG[j], qf, BE) end + return nothing +end + - if continuous == false - for n in 1:nnodes, k in 1:target_resultdim - if nneighbours[n] > 0 - target[k + target_offset, n] /= nneighbours[n] - end - end - end +""" +```` +function nodevalues!( + target::AbstractArray{<:Real,2}, + source::AbstractArray{T,1}, + FE::FESpace{Tv,Ti,FEType,AT}, + operator::Type{<:AbstractFunctionOperator} = Identity; + abs::Bool = false, + factor = 1, + regions::Array{Int,1} = [0], + target_offset::Int = 0, + source_offset::Int = 0, + zero_target::Bool = true, + continuous::Bool = false + ) +```` - if abs - for n in 1:nnodes - target[1 + target_offset, n] = sqrt(target[1 + target_offset, n]) - end - end +Evaluate a finite element function (given by the coefficient vector `source` and FE space `FE`) at all nodes of the grid, applying an optional function operator, and write the results into `target`. - return nothing -end +By default, the function is evaluated at every node in the mesh. If `continuous` is `false`, the value at each node is averaged over all neighboring cells (suitable for discontinuous quantities). If `continuous` is `true`, the value is taken from a single cell per node (suitable for continuous quantities). The result can optionally be scaled, offset, or restricted to specific regions. +# Arguments +- `target`: Output array to store the evaluated values (size: result dimension × number of nodes). +- `source`: Coefficient vector for the FE function. +- `FE`: The finite element space. +- `operator`: Function operator to apply at each node (default: `Identity`). -function piecewise_nodevalues!( +# Keyword Arguments +- `abs`: If `true`, store the Euclidean norm of the result at each node (default: `false`). +- `factor`: Scaling factor applied to the result (default: `1`). +- `regions`: List of region indices to restrict evaluation (default: all regions). +- `target_offset`: Offset for writing into `target` (default: `0`). +- `source_offset`: Offset for reading from `source` (default: `0`). +- `zero_target`: If `true`, zero out `target` before writing (default: `true`). +- `continuous`: If `true`, evaluate only once per node; otherwise, average over all neighboring cells (default: `false`). + +# Notes +- The result dimension is determined by the FE space, the operator, and the `abs` argument. +- The function modifies `target` in-place. +- For vector-valued or higher-dimensional results, the first dimension of `target` corresponds to the result dimension. + +""" +function nodevalues!( target::AbstractArray{T, 2}, source::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, AT}, @@ -435,7 +509,6 @@ function piecewise_nodevalues!( node::Int = 0 dof::Ti = 0 weights::Array{T, 1} = qf.w - nweights = length(weights) basisvals = BE.cvals ndofs = size(basisvals, 2) cvals_resultdim::Int = size(basisvals, 1) @@ -455,32 +528,29 @@ function piecewise_nodevalues!( for i in eachindex(weights) # vertices node = xItemNodes[i, item] fill!(localT, 0) - - for dof_i in 1:ndofs - dof = xItemDofs[dof_i, item] - eval_febe!(temp, BE, dof_i, i) - for k in 1:cvals_resultdim - localT[k] += source[source_offset + dof] * temp[k] - #target[k+target_offset,node] += temp[k] * source[source_offset + dof] - end - end - localT .*= factor - if abs - for k in 1:cvals_resultdim - target[target_offset + i, item] += localT[k]^2 + if continuous == false || flag4node[node] == false + nneighbours[node] += 1 + flag4node[node] = true + for dof_i in 1:ndofs + dof = xItemDofs[dof_i, item] + eval_febe!(temp, BE, dof_i, i) + for k in 1:cvals_resultdim + localT[k] += source[source_offset + dof] * temp[k] + #target[k+target_offset,node] += temp[k] * source[source_offset + dof] + end end - else - for k in 1:cvals_resultdim - target[target_offset + i + (k - 1) * nweights, item] += localT[k] + localT .*= factor + if abs + for k in 1:cvals_resultdim + target[1 + target_offset, node] += localT[k]^2 + end + else + for k in 1:cvals_resultdim + target[k + target_offset, node] += localT[k] + end end end end - - if abs - for i in 1:nweights - target[i + target_offset, item] = sqrt(target[i + target_offset, item]) - end - end break # region for loop end # if in region end # region for loop @@ -494,62 +564,69 @@ function piecewise_nodevalues!( barrier(EG[j], qf, BE) end + + if continuous == false + for n in 1:nnodes, k in 1:target_resultdim + if nneighbours[n] > 0 + target[k + target_offset, n] /= nneighbours[n] + end + end + end + + if abs + for n in 1:nnodes + target[1 + target_offset, n] = sqrt(target[1 + target_offset, n]) + end + end + return nothing end - """ ```` -function nodevalues!( - target::AbstractArray{<:Real,2}, - source::FEVectorBlock, - operator::Type{<:AbstractFunctionOperator} = Identity; - regions::Array{Int,1} = [0], - abs::Bool = false, - factor = 1, - cellwise = false, # return cellwise nodevalues ncells x nnodes_on_cell - target_offset::Int = 0, # start to write into target after offset - zero_target::Bool = true, # target vector is zeroed - continuous::Bool = false) + nodevalues( + source::FEVectorBlock, + operator::Type{<:AbstractFunctionOperator} = Identity; + continuous = "auto", + nodes = [], + cellwise = false, + abs = false, + kwargs... + ) ```` -Evaluates the finite element function with the coefficient vector source -and the specified FunctionOperator at all the nodes of the (specified regions of the) grid and writes the values into target. -Discontinuous (continuous = false) quantities are evaluated in all neighbouring cells of each node and then averaged. Continuous -(continuous = true) quantities are only evaluated once at each node. -""" -function nodevalues!(target, source::FEVectorBlock, operator::Type{<:AbstractFunctionOperator} = Identity; cellwise = false, kwargs...) - return if cellwise - piecewise_nodevalues!(target, source.entries, source.FES, operator; source_offset = source.offset, kwargs...) - else - nodevalues!(target, source.entries, source.FES, operator; source_offset = source.offset, kwargs...) - end -end +Evaluate a finite element function (given by the coefficient vector `source`) at nodes of the grid, applying an optional function operator, and return the result as a newly allocated array of the appropriate size. +This function provides a flexible interface for extracting nodal or cellwise values from a finite element solution. By default, it evaluates at all nodes, but a subset of nodes can be specified. The result can be returned in a cellwise (piecewise) layout if desired. -""" -```` -function nodevalues( - source::FEVectorBlock, - operator::Type{<:AbstractFunctionOperator} = Identity; - regions::Array{Int,1} = [0], - abs::Bool = false, - factor = 1, - nodes = [], - cellwise = false, # return cellwise nodevalues ncells x nnodes_on_cell (only if nodes == []) - target_offset::Int = 0, # start to write into target after offset - zero_target::Bool = true, # target vector is zeroed - continuous::Bool = false) -```` +# Arguments +- `source`: The `FEVectorBlock` containing the coefficients of the FE function. +- `operator`: The function operator to apply at each node (default: `Identity`). + +# Keyword Arguments +- `continuous`: If `"auto"`, automatically choose continuous/discontinuous evaluation based on FE type and operator. If `true`, evaluate only once per node; if `false`, average over all neighboring cells. Is ignored when `cellwise` is `true`. +- `nodes`: List of node indices to evaluate (default: all nodes). +- `cellwise`: If `true`, return values in a cellwise (piecewise) layout (default: `false`). +- `abs`: If `true`, return the Euclidean norm at each node (default: `false`). +- `kwargs...`: Additional keyword arguments passed to lower-level routines (e.g., `regions`, `factor`, `target_offset`, `zero_target`, etc.). -Evaluates the finite element function with the coefficient vector source -and the specified FunctionOperator at the specified list of nodes of the grid (default = all nodes) -and writes the values in that order into target. Nodes that are not part of the specified regions (default = all regions) -are set to zero. -Discontinuous (continuous = false) quantities are evaluated in all neighbouring cells of each node and then averaged. Continuous -(continuous = true) quantities are only evaluated once at each node. +# Returns +- A newly allocated array containing the evaluated values, with shape depending on the options chosen. + +# Notes +- If `nodes` is empty, all nodes are evaluated. +- If `cellwise` is `true`, the result is organized per cell (suitable for discontinuous or element-wise quantities). +- The result dimension is determined by the FE space, the operator, and the `abs` argument. """ -function nodevalues(source::FEVectorBlock{T, Tv, Ti, FEType, APT}, operator::Type{<:AbstractFunctionOperator} = Identity; continuous = "auto", nodes = [], cellwise = false, abs = false, kwargs...) where {T, Tv, Ti, APT, FEType} +function nodevalues( + source::FEVectorBlock{T, Tv, Ti, FEType, APT}, + operator::Type{<:AbstractFunctionOperator} = Identity; + continuous = "auto", + nodes = [], + cellwise = false, + abs = false, + kwargs... + ) where {T, Tv, Ti, APT, FEType} if continuous == "auto" if FEType <: AbstractH1FiniteElement && operator == Identity && !source.FES.broken && !(FEType <: H1CR) continuous = true @@ -567,15 +644,15 @@ function nodevalues(source::FEVectorBlock{T, Tv, Ti, FEType, APT}, operator::Typ if nodes == [] if cellwise target = zeros(T, nvals * max_num_targets_per_source(source.FES.dofgrid[CellNodes]), num_cells(source.FES.dofgrid)) - piecewise_nodevalues!(target, source.entries, source.FES, operator; continuous = continuous, source_offset = source.offset, abs = abs, kwargs...) + piecewise_nodevalues!(target, source.entries, source.FES, operator; abs = abs, kwargs...) else target = zeros(T, nvals, num_nodes(source.FES.dofgrid)) - nodevalues!(target, source.entries, source.FES, operator; continuous = continuous, source_offset = source.offset, abs = abs, kwargs...) + nodevalues!(target, source.entries, source.FES, operator; continuous = continuous, abs = abs, kwargs...) end else @assert cellwise == false target = zeros(T, nvals, length(nodes)) - nodevalues_subset!(target, source.entries, source.FES, operator; continuous = continuous, source_offset = source.offset, abs = abs, nodes = nodes, kwargs...) + nodevalues_subset!(target, source.entries, source.FES, operator; continuous = continuous, abs = abs, nodes = nodes, kwargs...) end return target end @@ -587,7 +664,23 @@ function nodevalues_view( operator::Type{<:AbstractFunctionOperator} = Identity) ```` -Returns a vector of views of the nodal values of the source block (currently works for unbroken H1-conforming elements) that directly accesses the coefficients. +Return a vector of views into the nodal values of the given finite element function, allowing direct access to the underlying coefficient storage for each component. + +This function provides efficient, zero-copy access to the nodal values of an `FEVectorBlock` for unbroken H1-conforming finite element spaces with the identity operator. Each entry in the returned vector is a view into the coefficients corresponding to one component of the FE function, optionally restricted to a subset of nodes. + +# Arguments +- `source`: The `FEVectorBlock` containing the coefficients of the FE function. +- `operator`: The function operator to apply (must be `Identity` for direct views; default: `Identity`). + +# Keyword Arguments +- `nodes`: List of node indices to view (default: all nodes). + +# Returns +- A vector of `SubArray` views, one for each component, directly referencing the coefficients for the specified nodes. + +# Notes +- Only available for unbroken H1-conforming elements and the Identity operator. + """ function nodevalues_view(source::FEVectorBlock{T, Tv, Ti, FEType, APT}, operator::Type{<:AbstractFunctionOperator} = Identity; nodes = [0]) where {T, Tv, Ti, APT, FEType} @@ -635,8 +728,24 @@ function continuify( regions::Array{Int,1} = [0]) where {T,Tv,Ti,FEType,APT} ```` -interpolates operator evaluation of source into a FE function of FEType H1Pk, i.e., Lagrange interpolation of arbitrary -operator evaluations of the source finite element type, broken = true generates a piecewise interpolation +Interpolate the evaluation of an operator applied to a finite element function onto a continuous Lagrange finite element space (`H1Pk`), returning a new `FEVector` with the interpolated values. + +This function performs nodal interpolation of the (possibly vector-valued) result of applying `operator` to the FE function represented by `source`. The result is a new FE function in a continuous Lagrange space of the specified order and dimension. If `broken = true`, the interpolation is performed in a piecewise (discontinuous) fashion. + +# Arguments +- `source`: The `FEVectorBlock` containing the coefficients of the original FE function. +- `operator`: The function operator to apply before interpolation (default: `Identity`). + +# Keyword Arguments +- `abs`: If `true`, interpolate the Euclidean norm of the result (default: `false`). +- `broken`: If `true`, generate a piecewise (discontinuous) interpolation (default: `false`). +- `order`: Polynomial order of the target Lagrange space (default: `"auto"`, which chooses an appropriate order). +- `factor`: Scaling factor applied to the result before interpolation (default: `1`). +- `regions`: List of region indices to restrict interpolation (default: all regions). + +# Returns +- A new `FEVector` in a continuous (or broken) Lagrange space, containing the interpolated values. + """ function continuify( source::FEVectorBlock{T, Tv, Ti, FEType, APT}, @@ -776,8 +885,17 @@ end ```` function displace_mesh!(xgrid::ExtendableGrid, source::FEVectorBlock; magnify = 1) ```` -Moves all nodes of the grid by adding the displacement field in source (expects a vector-valued finite element) -times a magnify value. +Displace all nodes of the given grid by adding a vector-valued finite element field as a displacement, scaled by an optional magnification factor. + +This function modifies the coordinates of `xgrid` in-place by adding the nodal values of the FE function represented by `source` (typically a displacement field) to each node. The displacement can be scaled by the `magnify` parameter. After the update, all cached geometric quantities in the grid are invalidated and recomputed. + +# Arguments +- `xgrid`: The `ExtendableGrid` whose node coordinates will be updated. +- `source`: An `FEVectorBlock` representing the displacement field (must be vector-valued and defined on the grid). + +# Keyword Arguments +- `magnify`: Scaling factor for the displacement field (default: `1`). + """ function displace_mesh!(xgrid::ExtendableGrid, source::FEVectorBlock; magnify = 1) nnodes = size(xgrid[Coordinates], 2) @@ -800,8 +918,18 @@ end ```` function displace_mesh(xgrid::ExtendableGrid, source::FEVectorBlock; magnify = 1) ```` -Returns a new grid by adding the displacement field in source (expects a vector-valued finite element) -to the coordinates of the provided xgrid times a magnify value. +Return a new grid with node coordinates displaced by a vector-valued finite element field, optionally scaled by a magnification factor. + +# Arguments +- `xgrid`: The `ExtendableGrid` to be copied and displaced. +- `source`: An `FEVectorBlock` representing the displacement field (must be vector-valued and defined on the grid). + +# Keyword Arguments +- `magnify`: Scaling factor for the displacement field (default: `1`). + +# Returns +- A new `ExtendableGrid` with displaced node coordinates. + """ function displace_mesh(xgrid::ExtendableGrid, source::FEVectorBlock; kwargs...) xgrid_displaced = deepcopy(xgrid) diff --git a/src/interpolators.jl b/src/interpolators.jl index 211d6b9..1a56e83 100644 --- a/src/interpolators.jl +++ b/src/interpolators.jl @@ -22,9 +22,27 @@ end function NodalInterpolator(FES::FESpace{T}, grid = FES.dofgrid; broken = FES.broken, component_offset = FES.coffset, kwargs...) ```` -Constructs an interpolator structure that has an evaluate! function that evaluates some -function at the requested nodes (items) of the grid. In broken mode only ON_CELL the items -refer to cells and all nodes of these cells are evaluated. +Constructs a nodal interpolation operator for a given finite element space. The resulting object provides an `evaluate!` function that interpolates a user-supplied function at the nodal points (degrees of freedom) of the provided grid. + +# Arguments +- `FES::FESpace{T}`: The finite element space for which the interpolator is constructed. +- `grid`: The grid or mesh on which interpolation is performed. Defaults to `FES.dofgrid`. + +# Keywords +- `broken`: If `true`, interpolation is performed in "broken" mode, i.e., independently on each cell (typically for discontinuous spaces). Defaults to `FES.broken`. +- `component_offset`: Offset for vector-valued spaces, specifying the stride between components. Defaults to `FES.coffset`. +- `kwargs...`: Additional keyword arguments passed to internal structures (e.g., `QPInfos`). + +# Returns +- A `NodalInterpolator` struct containing an `evaluate!` function with the signature: + `evaluate!(target, exact_function!, items; time=0, params=[], kwargs...)` + which fills `target` with the values of `exact_function!` at the nodal points specified by `items`. + +# Notes +- In "broken" mode, `items` refers to cell indices, and all nodes of each cell are evaluated. +- In continuous mode, `items` refers to global node indices. +- The `exact_function!` should have the signature `exact_function!(result, QP)` where `QP` is a `QPInfos` object. + """ function NodalInterpolator( FES::FESpace{T}, @@ -113,11 +131,38 @@ end function MomentInterpolator(FE::FESpace{Tv, Ti, FEType, APT}, AT::Type{<:AssemblyType}, xgrid = FE.dofgrid; operator = Identity, FEType_ref = :auto, FEType_moments = :auto, moments_operator = operator, moments_dofs = Int[], bestapprox = false, order = 0, coffset::Int = -1, componentwise = true, kwargs...) where {Tv, Ti, FEType <: AbstractFiniteElement, APT} ```` -Constructs an interpolation structure that has an evaluate! function that sets the interior degrees of freedom -such that the moments of the given function are preserved up to the given order. For this small local problems -are solved with a mass matrix of interior basis functions (evaluated with operator) and (by moments_dofs selected) -basis functions of the moments of type FEType_ref and with moments_operator. In the bestapprox mode -the mass matrix instead is formed from the scalar product of the interior basis functions. +Constructs a moment-based interpolation operator for a given finite element space. The resulting object provides an `evaluate!` function that sets the interior degrees of freedom (DOFs) so that the moments of a user-supplied function are preserved up to a specified order. This is achieved by solving small local problems involving a mass matrix of interior basis functions and selected moment basis functions. + +# Arguments +- `FE::FESpace{Tv, Ti, FEType, APT}`: The finite element space for which the interpolator is constructed. +- `AT::Type{<:AssemblyType}`: The assembly type (e.g., `ON_CELLS`, `ON_FACES`) specifying the geometric entity for interpolation. +- `xgrid`: The grid or mesh on which interpolation is performed. Defaults to `FE.dofgrid`. + +# Keywords +- `operator`: Operator used to evaluate the basis functions (default: `Identity`). +- `FEType_ref`: Reference finite element type for the moments (default: `:auto`). +- `FEType_moments`: Finite element type for the moment basis (default: `:auto`). +- `moments_operator`: Operator for evaluating the moment basis (default: `operator`). +- `moments_dofs`: Indices of moment DOFs to use (default: all). +- `bestapprox`: If `true`, uses best-approximation (L2 projection) for interior DOFs (default: `false`). +- `order`: Order of moments to preserve (default: `0`). +- `coffset`: Component offset for vector-valued spaces (default: `-1`, auto-detected). +- `componentwise`: If `true`, moments are enforced componentwise (default: `true`). +- `kwargs...`: Additional keyword arguments passed to internal structures (e.g., `QPInfos`). + +# Returns +- A `MomentInterpolator` struct containing an `evaluate!` function with the signature: + `evaluate!(target, exact_function!, items; time=0, quadorder=..., params=[], bonus_quadorder=0, kwargs...)` + which fills `target` with DOFs such that the prescribed moments of `exact_function!` are matched on the specified entities. + +# Notes +- The `exact_function!` should have the signature `exact_function!(result, QP)` where `QP` is a `QPInfos` object. +- In `bestapprox` mode, the mass matrix is formed from the scalar product of the interior basis functions (L2 projection). +- The interpolator currently supports grids with a single element geometry. +- The `items` argument specifies which entities (cells/faces) to interpolate; if empty, all are used. +- The moment basis and reference basis types are auto-selected for common cases, but can be overridden. +- The resulting interpolator is useful for constructing higher-order or non-nodal interpolations. + """ function MomentInterpolator( FE::FESpace{Tv, Ti, FEType, APT}, @@ -453,11 +498,33 @@ function FunctionalInterpolator( operator = NormalFlux, nfluxes = 0, dofs = [], kwargs...) where {Tv, Ti, FEType <: AbstractFiniteElement, APT} ```` -Constructs an interpolation structure that has an evaluate! function that determines the interior degrees of freedom -(or the specified local dofs) by evaluating functionals! (the result dimension nfluxes should correspond to the number of dofs, -both are set to the number of interior dofs as defaults). -The functionals are corrected by the operator evaluations of the fixed dofs. -The mean parameter decides if the functional is divided by the area in the end (if mean = true). +Constructs a FunctionalInterpolator for a given finite element space. The resulting object provides an `evaluate!` function that sets the interior degrees of freedom (DOFs) or the specified local DOFs by evaluating the supplied functionals. The number of functionals (`nfluxes`) should match the number of DOFs to be set (by default, all interior DOFs). The functionals are corrected by subtracting the operator evaluations of the fixed DOFs. Optionally, the result can be averaged by the entity volume if `mean = true`. + +# Arguments +- `functionals!::Function`: A function with signature `functionals!(result, values, QP)` that computes the functionals to be interpolated, where `result` is filled with the functional values, `values` is the evaluation of the target function, and `QP` is a `QPInfos` object. +- `FE::FESpace{Tv, Ti, FEType, APT}`: The finite element space for which the interpolator is constructed. +- `AT::Type{<:AssemblyType}`: The assembly type (e.g., `ON_FACES`, `ON_CELLS`) specifying the geometric entity for interpolation. Defaults to `ON_FACES`. +- `xgrid`: The grid or mesh on which interpolation is performed. Defaults to `FE.dofgrid`. + +# Keywords +- `operator`: Operator used to evaluate the basis functions (default: `NormalFlux`). +- `nfluxes`: Number of functionals/DOFs to interpolate (default: number of interior DOFs). +- `dofs`: Indices of DOFs to set (default: all interior DOFs). +- `mean`: If `true`, divides the functional by the entity volume (default: `false`). +- `bonus_quadorder`: Additional quadrature order for integration (default: `0`). +- `kwargs...`: Additional keyword arguments passed to internal structures (e.g., `QPInfos`). + +# Returns +- A `FunctionalInterpolator` struct containing an `evaluate!` function with the signature: + `evaluate!(target, exact_function!, items; time=0, quadorder=..., params=[], bonus_quadorder=0, kwargs...)` + which fills `target` with DOFs such that the prescribed functionals of `exact_function!` are matched on the specified entities. + +# Notes +- The `exact_function!` should have the signature `exact_function!(result, QP)` where `QP` is a `QPInfos` object. +- The `items` argument specifies which entities (cells/faces) to interpolate; if empty, all are used. +- The interpolator is useful for constructing DOFs associated with functionals (e.g., fluxes, averages) rather than nodal values. +- The number of functionals and DOFs must match; both default to the number of interior DOFs if not specified. + """ function FunctionalInterpolator( functionals!::Function, diff --git a/src/lazy_interpolate.jl b/src/lazy_interpolate.jl index ba1702c..d85b35a 100644 --- a/src/lazy_interpolate.jl +++ b/src/lazy_interpolate.jl @@ -17,28 +17,44 @@ function lazy_interpolate!( ```` Interpolates (operator-evaluations of) the given FEVector source (or an array of FEVectorBlocks) -into the (finite element space assigned to) the `target` FEVectorBlock. -By the given `postprocess!` function that is conforming to the interface +into the finite element space assigned to the `target` FEVectorBlock. - postprocess!(result, input, qpinfo) +The interpolation is performed using a point evaluation pattern and cell search. If `CellParents` information +is available in the target grid, enabling `use_cellparents=true` can improve the efficiency of the search. -the operator evaluations (=input) can be further manipulated (default is unmodified input). The qpinfo argument -allows to access information at the current quadrature point. The `xtrafo!` function with the interface +A custom postprocessing function can be provided via the `postprocess!` argument, which should have the interface: + postprocess!(result, input, qpinfo) +where `result` is the output buffer, `input` is the operator evaluation, and `qpinfo` provides quadrature point information. +If the source and target grids have different coordinate dimensions, a coordinate transformation function +`xtrafo!` must be provided, with the interface: xtrafo!(x_source, x) +which maps coordinates `x` from the target grid to coordinates in the source grid. -maps coordinates x from the target grid to coordinates in the source grid in case the grids -(the default is the identity). If `x_source` cannot be found in the source_grid the value -`not_in_domain_value` is used as a function value. With the `items` arguments the -target cells for the interpolation can be restricted. +If a point cannot be found in the source grid, the value `not_in_domain_value` is used as the function value. +The `items` argument can be used to restrict the interpolation to specific target cells. -Note: discontinuous quantities at vertices of the target grid will be evaluated in the first found cell of the -source grid. No averaging is performed. With eps the tolerances of the cell search via ExtendableGrids.CellFinder can be steered. +# Arguments +- `target::FEVectorBlock`: The target finite element vector block. +- `source`: The source array of FEVectorBlocks. +- `operators`: Array of operator argument tuples (source block tag, operator type) (default: `[(1, Identity)]`). + +# Keyword Arguments +- `postprocess!`: Function to postprocess operator evaluations (default: `standard_kernel`). +- `xtrafo!`: Optional coordinate transformation function (default: `nothing`). +- `items`: List of target cells to interpolate (default: `[]`). +- `resultdim`: Result dimension (default: `get_ncomponents(eltype(target.FES))`). +- `not_in_domain_value`: Value assigned if a point is not found in the source domain (default: `1.0e30`). +- `start_cell`: Starting cell index for cell search (default: `1`). +- `only_localsearch`: Restrict cell search to local neighborhood (default: `false`). +- `use_cellparents`: Use parent cell information for search (default: `false`). +- `eps`: Tolerance for cell search (default: `1.0e-13`). +- `kwargs...`: Additional keyword arguments passed to `interpolate!`. + +# Notes +- Discontinuous quantities at target grid vertices are evaluated in the first found cell of the source grid; no averaging is performed. +- The function is not the most efficient for large-scale problems due to its reliance on pointwise cell search. -Note 2: This is not the most efficient way (therefore lazy) as it is based on the PointEvaluation pattern and cell search (with -tolerance `eps`). -If CellParents are available in the grid components of the target grid, parent cell information can be used to (slightly) improve the -search. To activate this put `use_cellparents` = true. """ function lazy_interpolate!( target::FEVectorBlock{T1, Tv, Ti}, diff --git a/src/point_evaluator.jl b/src/point_evaluator.jl index 1bc8095..8e8758b 100644 --- a/src/point_evaluator.jl +++ b/src/point_evaluator.jl @@ -25,26 +25,29 @@ default_peval_kwargs() = Dict{Symbol, Tuple{Any, String}}( ```` function Pointevaluator( [kernel!::Function], - oa_args::Array{<:Tuple{<:Any, DataType},1}; + oa_args::Array{<:Tuple{<:Any, DataType},1} + sol = nothing; kwargs...) ```` -Generates a PointEvaluator that can evaluate the specified operator evaluations -at arbitrary points. If no kernel function is given, the arguments -are given directly. If a kernel is provided, the arguments are postprocessed -accordingly and the kernel has to be conform to the interface +Constructs a `PointEvaluator` object for evaluating operator expressions at arbitrary points in a finite element space. - kernel!(result, eval_args, qpinfo) +A `PointEvaluator` can be used to evaluate one or more operator evaluations (e.g., function values, gradients) at arbitrary points, optionally postprocessed by a user-supplied kernel function. The evaluation is performed with respect to a given solution and finite element basis. -where qpinfo allows to access information at the current evaluation point. -Additionally the length of the result needs to be specified via the kwargs. +If a kernel function is provided, it must have the interface: + kernel!(result, eval_args, qpinfo) +where `result` is the output buffer, `eval_args` contains the operator evaluations, and `qpinfo` provides information about the current evaluation point. -Evaluation can be triggered via the evaluate function after an initialize! call. +Operator evaluations are specified as tuples pairing a tag (to identify the unknown or vector position) with a `FunctionOperator`. -Operator evaluations are tuples that pair a tag (to identify an unknown or the position in the vector) -with a FunctionOperator. +After construction, the `PointEvaluator` must be initialized with a solution using `initialize!`. Evaluation at a point is then performed using `evaluate!` or `evaluate_bary!`. -Keyword arguments: +# Arguments +- `kernel!` (optional): Postprocessing function for operator evaluations. +- `oa_args: Array of operator argument tuples (source block tag, operator type). +- `sol` (optional): Solution object for immediate initialization. + +# Keyword arguments: $(_myprint(default_peval_kwargs())) """ function PointEvaluator(kernel, u_args, ops_args, sol = nothing; Tv = Float64, kwargs...) @@ -78,7 +81,20 @@ function initialize!( kwargs...) ```` -Initializes the PointEvaluator for the specified solution. +Initializes the given `PointEvaluator` for a specified solution (FEVector or vector of FEVectorBlocks). + +This function prepares the `PointEvaluator` for evaluation by associating it with the provided solution vector. It sets up the necessary finite element basis evaluators, local-to-global transformations, and cell finder structures for the underlying grid. + +# Arguments +- `O::PointEvaluator`: The `PointEvaluator` instance to initialize. +- `sol`: The solution object (e.g., array of FEVectorBlocks) to be used for evaluations. + +# Keyword Arguments +- `time`: (default: `0`) Time value to be passed to the quadrature point info structure. +$(_myprint(default_peval_kwargs())) + +# Notes +- This function must be called before using `evaluate!` or `evaluate_bary!` with the `PointEvaluator`. """ function initialize!(O::PointEvaluator{T, UT}, sol; time = 0, kwargs...) where {T, UT} _update_params!(O.parameters, kwargs) diff --git a/src/qpinfos.jl b/src/qpinfos.jl index f86c38c..5f57711 100644 --- a/src/qpinfos.jl +++ b/src/qpinfos.jl @@ -1,11 +1,26 @@ """ $(TYPEDEF) -object that shares information about the current quadrature point -(item number w.r.t. current AssemblyType, cell number when reasonable, -region number w.r.t. current AssemblyType, volume if the item, -normal when reasonable, current time, current world coordinates, -current reference coordinates, set of parameters) +A mutable struct that encapsulates information about the current quadrature point during finite element assembly or evaluation. + +# Fields +- `item::Ti`: Index of the current item (with respect to the current assembly type). +- `cell::Ti`: Index of the current cell, if applicable. +- `region::Ti`: Index of the current region (with respect to the current assembly type). +- `volume::TvG`: Volume associated with the current item. +- `normal::Vector{TvG}`: Outward normal vector at the quadrature point, if applicable (e.g., for face assembly/integration). +- `time::Ttime`: Current time value, useful for time-dependent problems. +- `x::Vector{Tx}`: World (physical) coordinates of the quadrature point. +- `xref::Vector{Txref}`: Reference (local) coordinates of the quadrature point within the reference element. +- `grid::ExtendableGrid{TvG, TiG}`: Reference to the underlying grid or mesh structure. +- `params::PT`: Additional user- or problem-specific parameters. + +# Description +`QPInfos` is used to pass contextual information about the current quadrature point to finite element kernels, integrators, and evaluators. It provides access to geometric, topological, and problem-specific data required for local assembly, evaluation, or postprocessing routines. + +# Notes +- Not all fields are always meaningful for every assembly type (e.g., `normal` may be unused for volume integrals). + """ mutable struct QPInfos{Ti, Tv, Ttime, Tx, Txref, TvG, TiG, PT} item::Ti diff --git a/src/quadrature.jl b/src/quadrature.jl index 22e771e..9bf58e3 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -645,22 +645,29 @@ end """ $(TYPEDSIGNATURES) -Integration of an integrand of the signature - - integrand!(result, qpinfo) - -over the entities AT of the grid. The result on every item is written into integral4items -(at least of size nitems x length of result). -As usual, qpinfo allows to access information -at the current quadrature point. - -Keyword arguments: -- quadorder = specifies the quadrature order (default = 0) -- regions = restricts integration to these regions (default = [] meaning all regions) -- items = restricts integration to these item numbers (default = [] meaning all items) -- time = sets the time to be passed to qpinfo (default = 0) -- params = sets the params array to be passed to qpinfo (default = []) -- offset = offset for writing into result array integral4items (default = [0]), +Compute cellwise (or per-entity) integrals of a user-supplied integrand over entities of type `AT` in the given `grid`, writing the result for each entity into `integral4items`. + +# Arguments +- `integral4items::AbstractArray{T}`: Preallocated array to store the integral for each entity (e.g., cell, edge, or face). The shape should be compatible with the number of entities and the integrand's output dimension. +- `grid::ExtendableGrid{Tv, Ti}`: The grid or mesh over which to integrate. +- `AT::Type{<:AssemblyType}`: The entity type to integrate over (e.g., `ON_CELLS`, `ON_EDGES`, `ON_FACES`). +- `integrand!`: A function with signature `integrand!(result, qpinfo)` that computes the integrand at a quadrature point. The function should write its output into `result` (a preallocated vector) and use `qpinfo` to access quadrature point data. + +# Keyword Arguments +- `offset`: Offset(s) for writing into `integral4items` (default: `[0]`). +- `bonus_quadorder`: Additional quadrature order to add to `quadorder` (default: `0`). +- `quadorder`: Quadrature order (default: `0`). +- `regions`: Restrict integration to these region indices (default: `[]`, meaning all regions). +- `items`: Restrict integration to these item numbers (default: `[]`, meaning all items). +- `time`: Time value to be passed to `qpinfo` (default: `0`). +- `force_quadrature_rule`: Use this quadrature rule instead of the default (default: `nothing`). +- Additional keyword arguments (`kwargs...`) are forwarded to the quadrature point info constructor. + +# Notes +- The function loops over all specified entities (cells, edges, or faces), applies the quadrature rule, and accumulates the result for each entity in `integral4items`. +- The integrand function is called at each quadrature point and should write its output in-place to the provided result vector. +- For total (global) integrals, use [`integrate`](@ref) instead, which is more memory-efficient. +- The shape of `integral4items` determines whether the result is stored as a vector per entity or as a matrix (e.g., for multiple components). """ function integrate!( integral4items::AbstractArray{T}, @@ -823,20 +830,30 @@ end """ $(TYPEDSIGNATURES) -Integration that returns total integral of an integrand of the signature +Compute the total integral of a user-supplied integrand over entities of type `AT` in the given `grid`. - integrand!(result, qpinfo) +# Arguments +- `grid::ExtendableGrid`: The grid or mesh over which to integrate. +- `AT::Type{<:AssemblyType}`: The entity type to integrate over (e.g., `ON_CELLS`, `ON_EDGES`). +- `integrand!`: A function with signature `integrand!(result, qpinfo)` that computes the integrand at a quadrature point. The function should write its output into `result` (a preallocated vector) and use `qpinfo` to access quadrature point data. +- `resultdim::Int`: The length of the result vector expected from `integrand!` (i.e., the number of components to integrate). -over the entities AT of the grid. Here, qpinfo allows to access information -at the current quadrature point. -The length of the result needs to be specified via resultdim. +# Keyword Arguments +- `T`: The number type for accumulation (default: `Float64`). +- `quadorder`: Quadrature order (default: `0`). +- `regions`: Restrict integration to these region indices (default: `[]`, meaning all regions). +- `items`: Restrict integration to these item numbers (default: `[]`, meaning all items). +- `time`: Time value to be passed to `qpinfo` (default: `0`). +- `params`: Parameter array to be passed to `qpinfo` (default: `[]`). +- Additional keyword arguments are forwarded to the quadrature point info constructor. + +# Returns +- If `resultdim == 1`, returns a scalar value (the total integral). +- If `resultdim > 1`, returns a vector of length `resultdim` (componentwise integrals). + +# Notes +- This function is memory-efficient and accumulates the total integral directly, without storing per-entity results. For cellwise or per-entity integration, use [`integrate!`](@ref) instead. -Keyword arguments: -- quadorder = specifies the quadrature order (default = 0) -- regions = restricts integration to these regions (default = [] meaning all regions) -- items = restricts integration to these item numbers (default = [] meaning all items) -- time = sets the time to be passed to qpinfo (default = 0) -- params = sets the params array to be passed to qpinfo (default = []) """ function integrate( grid::ExtendableGrid, diff --git a/src/segment_integrator.jl b/src/segment_integrator.jl index 9e5d9c6..53b5335 100644 --- a/src/segment_integrator.jl +++ b/src/segment_integrator.jl @@ -44,27 +44,31 @@ function SegmentIntegrator( kwargs...) ```` -Generates an SegmentIntegrator that can integrate over segments -of the specified geometry EG. -To do so, it evaluates, at each quadrature point, the specified operator evaluations, -postprocesses them with the kernel function (if provided) -and accumulates the results with the quadrature weights. -If no kernel is given, the arguments -are integrated directly. If a kernel is provided it has be conform -to the interface +Construct a `SegmentIntegrator` for integrating over segments of the given element geometry `EG`. - kernel!(result, eval_args, qpinfo) +At each quadrature point, the specified operator evaluations are performed, optionally postprocessed by the provided `kernel!` function, and accumulated using the quadrature weights. If no kernel is given, the operator arguments are integrated directly. -where qpinfo allows to access information at the current quadrature point. -Additionally the length of the result needs to be specified via the kwargs. +# Arguments +- `EG::ElementGeometry`: The segment geometry over which to integrate. +- `kernel!::Function` (optional): A function of the form `kernel!(result, eval_args, qpinfo)` that postprocesses operator evaluations at each quadrature point. If omitted, a default kernel is used. +- `oa_args::Array{<:Tuple{Any, DataType},1}`: Array of tuples `(tag, FunctionOperator)`, where `tag` identifies the unknown or vector position, and `FunctionOperator` specifies the operator to evaluate. -Evaluation can be triggered via the integrate_segment! function after an initialize! +# Keyword arguments: +$(_myprint(default_segint_kwargs())) + +# Kernel Function Interface + + kernel!(result, eval_args, qpinfo) + +- `result`: Preallocated array to store the kernel output. +- `eval_args`: Array of operator evaluations at the current quadrature point. +- `qpinfo`: `QPInfos` struct with information about the current quadrature point. + +# Usage + +After construction, call `initialize!` to prepare the integrator for a given solution, then use `integrate_segment!` to perform integration over a segment. -Operator evaluations are tuples that pair a tag (to identify an unknown or the position in the vector) -with a FunctionOperator. -Keyword arguments: -$(_myprint(default_segint_kwargs())) """ function SegmentIntegrator(EG, kernel, oa_args::Array{<:Tuple{Union{<:Any, Int}, DataType}, 1}; kwargs...) u_args = [oa[1] for oa in oa_args] From 845d7c5e803e3bc63c9d64e442a0b7df9bbc026b Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Tue, 24 Jun 2025 22:25:03 +0200 Subject: [PATCH 2/8] removed duplicate nodevalues! --- src/interpolations.jl | 202 ++++++++++++++++++------------------------ 1 file changed, 88 insertions(+), 114 deletions(-) diff --git a/src/interpolations.jl b/src/interpolations.jl index 71cac3c..e6b193b 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -260,38 +260,13 @@ end """ ```` function nodevalues!( - target::AbstractArray{<:Real,2}, - source::AbstractArray{T,1}, - FE::FESpace{Tv,Ti,FEType,AT}, - operator::Type{<:AbstractFunctionOperator} = Identity; - kwargs...) -```` - -calls the nodevalues_subset! function with all nodes of the dofgrid, see nodevalues_subset!(target, source, FE, operator; nodes = 1:num_nodes(FE.dofgrid), kwargs...) -""" -function nodevalues!( - target::AbstractArray{T, 2}, - source::AbstractArray{T, 1}, - FE::FESpace{Tv, Ti, FEType, AT}, - operator::Type{<:AbstractFunctionOperator} = Identity; - kwargs... - ) where {T, Tv, Ti, FEType, AT} - - nodevalues_subset!(target, source, FE, operator; nodes = 1:num_nodes(FE.dofgrid), kwargs...) - - return nothing -end - -""" -```` -piecewise_nodevalues!( - target::AbstractArray{T, 2}, - source::AbstractArray{T, 1}, - FE::FESpace{Tv, Ti, FEType, AT}, + target::AbstractArray{<:Real,2}, + source::AbstractArray{T,1}, + FE::FESpace{Tv,Ti,FEType,AT}, operator::Type{<:AbstractFunctionOperator} = Identity; abs::Bool = false, factor = 1, - regions::Array{Int, 1} = [0], + regions::Array{Int,1} = [0], target_offset::Int = 0, source_offset::Int = 0, zero_target::Bool = true, @@ -299,13 +274,17 @@ piecewise_nodevalues!( ) ```` -Evaluate a finite element function (given by the coefficient vector `source` and FE space `FE`) at all nodes, but store the results in a piecewise (cellwise) fashion, i.e., for each cell, the values at its local nodes are written into `target`. The result is organized so that each column of `target` corresponds to a cell, and each row corresponds to the values at the cell's nodes. +Evaluate a finite element function (given by the coefficient vector `source` and FE space `FE`) at all nodes of the grid, applying an optional function operator, and write the results into `target`. + +By default, the function is evaluated at every node in the mesh. If `continuous` is `false`, the value at each node is averaged over all neighboring cells (suitable for discontinuous quantities). If `continuous` is `true`, the value is taken from a single cell per node (suitable for continuous quantities). The result can optionally be scaled, offset, or restricted to specific regions. # Arguments -- `target`: Output array to store the evaluated values (size: result dimension × number of nodes per cell, number of cells). +- `target`: Output array to store the evaluated values (size: result dimension × number of nodes). - `source`: Coefficient vector for the FE function. - `FE`: The finite element space. -- `operator`: Function operator to apply (default: `Identity`). +- `operator`: Function operator to apply at each node (default: `Identity`). + +# Keyword Arguments - `abs`: If `true`, store the Euclidean norm of the result at each node (default: `false`). - `factor`: Scaling factor applied to the result (default: `1`). - `regions`: List of region indices to restrict evaluation (default: all regions). @@ -316,8 +295,11 @@ Evaluate a finite element function (given by the coefficient vector `source` and # Notes - The result dimension is determined by the FE space, the operator, and the `abs` argument. +- The function modifies `target` in-place. +- For vector-valued or higher-dimensional results, the first dimension of `target` corresponds to the result dimension. + """ -function piecewise_nodevalues!( +function nodevalues!( target::AbstractArray{T, 2}, source::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, AT}, @@ -327,7 +309,8 @@ function piecewise_nodevalues!( regions::Array{Int, 1} = [0], target_offset::Int = 0, source_offset::Int = 0, - zero_target::Bool = true + zero_target::Bool = true, + continuous::Bool = false ) where {T, Tv, Ti, FEType, AT} xgrid = FE.dofgrid @@ -354,11 +337,14 @@ function piecewise_nodevalues!( fill!(target, 0) end + nnodes::Int = num_sources(xgrid[Coordinates]) + nneighbours::Array{Int, 1} = zeros(Int, nnodes) + flag4node::Array{Bool, 1} = zeros(Bool, nnodes) + function barrier(EG, qf, BE) node::Int = 0 dof::Ti = 0 weights::Array{T, 1} = qf.w - nweights = length(weights) basisvals = BE.cvals ndofs = size(basisvals, 2) cvals_resultdim::Int = size(basisvals, 1) @@ -378,32 +364,29 @@ function piecewise_nodevalues!( for i in eachindex(weights) # vertices node = xItemNodes[i, item] fill!(localT, 0) - - for dof_i in 1:ndofs - dof = xItemDofs[dof_i, item] - eval_febe!(temp, BE, dof_i, i) - for k in 1:cvals_resultdim - localT[k] += source[source_offset + dof] * temp[k] - #target[k+target_offset,node] += temp[k] * source[source_offset + dof] - end - end - localT .*= factor - if abs - for k in 1:cvals_resultdim - target[target_offset + i, item] += localT[k]^2 + if continuous == false || flag4node[node] == false + nneighbours[node] += 1 + flag4node[node] = true + for dof_i in 1:ndofs + dof = xItemDofs[dof_i, item] + eval_febe!(temp, BE, dof_i, i) + for k in 1:cvals_resultdim + localT[k] += source[source_offset + dof] * temp[k] + #target[k+target_offset,node] += temp[k] * source[source_offset + dof] + end end - else - for k in 1:cvals_resultdim - target[target_offset + i + (k - 1) * nweights, item] += localT[k] + localT .*= factor + if abs + for k in 1:cvals_resultdim + target[1 + target_offset, node] += localT[k]^2 + end + else + for k in 1:cvals_resultdim + target[k + target_offset, node] += localT[k] + end end end end - - if abs - for i in 1:nweights - target[i + target_offset, item] = sqrt(target[i + target_offset, item]) - end - end break # region for loop end # if in region end # region for loop @@ -417,20 +400,34 @@ function piecewise_nodevalues!( barrier(EG[j], qf, BE) end + + if continuous == false + for n in 1:nnodes, k in 1:target_resultdim + if nneighbours[n] > 0 + target[k + target_offset, n] /= nneighbours[n] + end + end + end + + if abs + for n in 1:nnodes + target[1 + target_offset, n] = sqrt(target[1 + target_offset, n]) + end + end + return nothing end - """ ```` -function nodevalues!( - target::AbstractArray{<:Real,2}, - source::AbstractArray{T,1}, - FE::FESpace{Tv,Ti,FEType,AT}, +piecewise_nodevalues!( + target::AbstractArray{T, 2}, + source::AbstractArray{T, 1}, + FE::FESpace{Tv, Ti, FEType, AT}, operator::Type{<:AbstractFunctionOperator} = Identity; abs::Bool = false, factor = 1, - regions::Array{Int,1} = [0], + regions::Array{Int, 1} = [0], target_offset::Int = 0, source_offset::Int = 0, zero_target::Bool = true, @@ -438,17 +435,13 @@ function nodevalues!( ) ```` -Evaluate a finite element function (given by the coefficient vector `source` and FE space `FE`) at all nodes of the grid, applying an optional function operator, and write the results into `target`. - -By default, the function is evaluated at every node in the mesh. If `continuous` is `false`, the value at each node is averaged over all neighboring cells (suitable for discontinuous quantities). If `continuous` is `true`, the value is taken from a single cell per node (suitable for continuous quantities). The result can optionally be scaled, offset, or restricted to specific regions. +Evaluate a finite element function (given by the coefficient vector `source` and FE space `FE`) at all nodes, but store the results in a piecewise (cellwise) fashion, i.e., for each cell, the values at its local nodes are written into `target`. The result is organized so that each column of `target` corresponds to a cell, and each row corresponds to the values at the cell's nodes. # Arguments -- `target`: Output array to store the evaluated values (size: result dimension × number of nodes). +- `target`: Output array to store the evaluated values (size: result dimension × number of nodes per cell, number of cells). - `source`: Coefficient vector for the FE function. - `FE`: The finite element space. -- `operator`: Function operator to apply at each node (default: `Identity`). - -# Keyword Arguments +- `operator`: Function operator to apply (default: `Identity`). - `abs`: If `true`, store the Euclidean norm of the result at each node (default: `false`). - `factor`: Scaling factor applied to the result (default: `1`). - `regions`: List of region indices to restrict evaluation (default: all regions). @@ -459,11 +452,8 @@ By default, the function is evaluated at every node in the mesh. If `continuous` # Notes - The result dimension is determined by the FE space, the operator, and the `abs` argument. -- The function modifies `target` in-place. -- For vector-valued or higher-dimensional results, the first dimension of `target` corresponds to the result dimension. - """ -function nodevalues!( +function piecewise_nodevalues!( target::AbstractArray{T, 2}, source::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, AT}, @@ -473,8 +463,7 @@ function nodevalues!( regions::Array{Int, 1} = [0], target_offset::Int = 0, source_offset::Int = 0, - zero_target::Bool = true, - continuous::Bool = false + zero_target::Bool = true ) where {T, Tv, Ti, FEType, AT} xgrid = FE.dofgrid @@ -501,14 +490,11 @@ function nodevalues!( fill!(target, 0) end - nnodes::Int = num_sources(xgrid[Coordinates]) - nneighbours::Array{Int, 1} = zeros(Int, nnodes) - flag4node::Array{Bool, 1} = zeros(Bool, nnodes) - function barrier(EG, qf, BE) node::Int = 0 dof::Ti = 0 weights::Array{T, 1} = qf.w + nweights = length(weights) basisvals = BE.cvals ndofs = size(basisvals, 2) cvals_resultdim::Int = size(basisvals, 1) @@ -528,27 +514,30 @@ function nodevalues!( for i in eachindex(weights) # vertices node = xItemNodes[i, item] fill!(localT, 0) - if continuous == false || flag4node[node] == false - nneighbours[node] += 1 - flag4node[node] = true - for dof_i in 1:ndofs - dof = xItemDofs[dof_i, item] - eval_febe!(temp, BE, dof_i, i) - for k in 1:cvals_resultdim - localT[k] += source[source_offset + dof] * temp[k] - #target[k+target_offset,node] += temp[k] * source[source_offset + dof] - end + + for dof_i in 1:ndofs + dof = xItemDofs[dof_i, item] + eval_febe!(temp, BE, dof_i, i) + for k in 1:cvals_resultdim + localT[k] += source[source_offset + dof] * temp[k] + #target[k+target_offset,node] += temp[k] * source[source_offset + dof] end - localT .*= factor - if abs - for k in 1:cvals_resultdim - target[1 + target_offset, node] += localT[k]^2 - end - else - for k in 1:cvals_resultdim - target[k + target_offset, node] += localT[k] - end + end + localT .*= factor + if abs + for k in 1:cvals_resultdim + target[target_offset + i, item] += localT[k]^2 end + else + for k in 1:cvals_resultdim + target[target_offset + i + (k - 1) * nweights, item] += localT[k] + end + end + end + + if abs + for i in 1:nweights + target[i + target_offset, item] = sqrt(target[i + target_offset, item]) end end break # region for loop @@ -564,21 +553,6 @@ function nodevalues!( barrier(EG[j], qf, BE) end - - if continuous == false - for n in 1:nnodes, k in 1:target_resultdim - if nneighbours[n] > 0 - target[k + target_offset, n] /= nneighbours[n] - end - end - end - - if abs - for n in 1:nnodes - target[1 + target_offset, n] = sqrt(target[1 + target_offset, n]) - end - end - return nothing end From b6e9d3499f583cb7a687f9b8b3e30ddf0832bcda Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Wed, 25 Jun 2025 08:30:55 +0200 Subject: [PATCH 3/8] fixed FEMatrix constructor issues leading to errors in downstream tests, some more documentation improvements --- README.md | 22 ++++++++-------------- docs/src/index.md | 2 -- src/fematrix.jl | 36 +++++++++++++++++++++++++++++++++--- 3 files changed, 41 insertions(+), 19 deletions(-) diff --git a/README.md b/README.md index efc80e3..09e5b36 100644 --- a/README.md +++ b/README.md @@ -6,19 +6,13 @@ # ExtendableFEMBase -This package provides basic finite element structures to setup finite element schemes on ExtendableGrids. For a full high-level API -see [ExtendableFEM.jl](https://github.com/WIAS-PDELib/ExtendableFEM.jl). +ExtendableFEMBase.jl provides foundational data structures and tools for assembling custom finite element solvers in Julia. It is designed for flexibility, efficiency, and extensibility, and serves as the low-level backend for [ExtendableFEM.jl](https://github.com/WIAS-PDELib/ExtendableFEM.jl). All functionality is built on top of [ExtendableGrids.jl](https://github.com/WIAS-PDELib/ExtendableGrids.jl). -This low level structures in the package incorporate: +## Features -- Finite element types (Basis functions on reference geometries and dof management for several H1, Hdiv and Hcurl elements) -- FESpace (Discrete finite element space with respect to a mesh from ExtendableGrids, knows the Dofmaps) -- FEMatrix (block overlay for an ExtendableSparse matrix, where each block corresponds to a coupling between two FESpaces in a system) -- FEVector (block overlay for an array, where each block corresponds to a FESpace) -- FunctionOperators (primitive linear operators like Identity, Gradient, Divergence) and rules how to evaluate them for for different finite element types -- FEEvaluator (finite element basis evaluators for different FunctionOperators and entities of the grid) -- QuadratureRule (basic quadrature rules for different ElementGeometries from ExtendableGrids) -- interpolations (standard interpolations into the provided finite element spaces, averaging routines and interpolations between meshes/FESpaces) -- reconstruction operators (special FunctionOperators that involve an interpolation into a different finite element type) - - +- Wide range of finite element types (H1, Hdiv, Hcurl, etc.) +- Flexible finite element spaces (`FESpace`) +- Block-structured matrices and vectors (`FEMatrix`, `FEVector`) +- Primitive and composite function operators (e.g., gradient, divergence) +- Efficient basis evaluation and quadrature routines +- Interpolation and reconstruction operators diff --git a/docs/src/index.md b/docs/src/index.md index b9ecd85..9f39ae6 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -7,8 +7,6 @@ ExtendableFEMBase.jl provides foundational data structures and tools for building custom finite element solvers in Julia. The package includes flexible representations for finite element spaces, interpolators, matrices, and vectors, all designed to work seamlessly with the [ExtendableGrids.jl](https://github.com/j-fu/ExtendableGrids.jl) infrastructure. -With ExtendableFEMBase.jl, users can efficiently assemble and manipulate finite element systems, enabling the development of advanced simulation workflows and research codes. - ### Dependencies on other Julia packages diff --git a/src/fematrix.jl b/src/fematrix.jl index a986ced..224a3cf 100644 --- a/src/fematrix.jl +++ b/src/fematrix.jl @@ -9,7 +9,22 @@ """ $(TYPEDEF) -block of an FEMatrix that carries coefficients for an associated pair of FESpaces and can be assigned as an two-dimensional AbstractArray (getindex, setindex, size) +A block of an `FEMatrix` representing the coupling between two finite element spaces. + +`FEMatrixBlock` acts as a two-dimensional array (subclassing `AbstractArray{TvM, 2}`) and stores the coefficients for a specific pair of row and column finite element spaces (`FESpace`). +Each block is mapped to a submatrix of the global sparse matrix, with offsets and sizes corresponding to the degrees of freedom of the associated spaces. + +# Fields +- `name::String`: Name of the block (for identification and display). +- `FES::FESpace{TvG, TiG, FETypeX, APTX}`: Row finite element space. +- `FESY::FESpace{TvG, TiG, FETypeY, APTY}`: Column finite element space. +- `offset::Int64`: Row offset in the global matrix. +- `offsetY::Int64`: Column offset in the global matrix. +- `last_index::Int64`: Last row index for this block. +- `last_indexY::Int64`: Last column index for this block. +- `entries::AbstractSparseMatrix{TvM, TiM}`: Reference to the underlying global sparse matrix (shared with the parent `FEMatrix`). + +See also: [`FEMatrix`], [`FESpace`] """ struct FEMatrixBlock{TvM, TiM, TvG, TiG, FETypeX, FETypeY, APTX, APTY} <: AbstractArray{TvM, 2} name::String @@ -40,7 +55,22 @@ end """ $(TYPEDEF) -an AbstractMatrix (e.g. an ExtendableSparseMatrix) with an additional layer of several FEMatrixBlock subdivisions each carrying coefficients for their associated pair of FESpaces +A block-structured sparse matrix type for finite element assembly. + +`FEMatrix` represents a (typically sparse) matrix with an additional layer of structure: it is subdivided into multiple `FEMatrixBlock`s, each associated with a specific pair of finite element spaces (`FESpace`). + +# Fields +- `FEMatrixBlocks::Array{FEMatrixBlock{TvM, TiM, TvG, TiG}, 1}`: The array of matrix blocks, each corresponding to a pair of row and column finite element spaces. +- `entries::AbstractSparseMatrix{TvM, TiM}`: The underlying sparse matrix storing all coefficients. +- `tags::Matrix{Any}`: Optional tags for identifying or accessing blocks. + +# Type Parameters +- `TvM`: Value type for matrix entries (e.g., `Float64`). +- `TiM`: Integer type for matrix indices (e.g., `Int64`). +- `TvG`, `TiG`: Value and index types for the associated finite element spaces. +- `nbrow`: Number of block rows. +- `nbcol`: Number of block columns. +- `nbtotal`: Total number of blocks. """ struct FEMatrix{TvM, TiM, TvG, TiG, nbrow, nbcol, nbtotal} <: AbstractSparseMatrix{TvM, TiM} FEMatrixBlocks::Array{FEMatrixBlock{TvM, TiM, TvG, TiG}, 1} @@ -204,7 +234,7 @@ function FEMatrix( return FEMatrix{TvM, TiM}(FESXv, FESYv; kwargs...) end -function FEMatrix{TvM, TiM}(FESX::Array{<:FESpace{TvG, TiG}, 1}, FESY::Array{<:FESpace{TvG, TiG}, 1}; entries = nothing, name = :automatic, tags = nothing, tagsX = tags, tagsY = tagsX, npartitions = 1, kwargs...) where {TvM, TiM, TvG, TiG} +function FEMatrix{TvM, TiM}(FESX::Array{<:FESpace{TvG, TiG}, 1}, FESY::Array{<:FESpace{TvG, TiG}, 1} = FESX; entries = nothing, name = :automatic, tags = nothing, tagsX = tags, tagsY = tagsX, npartitions = 1, kwargs...) where {TvM, TiM, TvG, TiG} ndofsX, ndofsY = 0, 0 for j in 1:length(FESX) ndofsX += FESX[j].ndofs From 79484cbd5613870009f6b565fd95d508a8f0e756 Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Wed, 25 Jun 2025 09:00:13 +0200 Subject: [PATCH 4/8] nodevalues! for FEVectorBlocks should work again, added range checks and Array conversions for FEVector --- src/fevector.jl | 45 +++++++++++++++++++++++++++++++++++++++++-- src/interpolations.jl | 15 ++++++++++++++- 2 files changed, 57 insertions(+), 3 deletions(-) diff --git a/src/fevector.jl b/src/fevector.jl index c53e436..ab8f7b0 100644 --- a/src/fevector.jl +++ b/src/fevector.jl @@ -88,8 +88,19 @@ function Base.copy(FEV::FEVector{T, Tv, Ti}) where {T, Tv, Ti} end # overload stuff for AbstractArray{T,1} behaviour -Base.getindex(FEF::FEVector{T, Tv, Ti}, tag) where {T, Tv, Ti} = FEF.FEVectorBlocks[findfirst(==(tag), FEF.tags)] -Base.getindex(FEF::FEVector, i::Int) = FEF.FEVectorBlocks[i] +Base.getindex(FEF::FEVector{T, Tv, Ti}, tag) where {T, Tv, Ti} = begin + idx = findfirst(==(tag), FEF.tags) + if idx === nothing + error("Tag '$(tag)' not found in FEVector. Available tags: $(FEF.tags)") + end + FEF.FEVectorBlocks[idx] +end +Base.getindex(FEF::FEVector, i::Int) = begin + if i < 1 || i > length(FEF.FEVectorBlocks) + error("Index $(i) out of bounds for FEVector with $(length(FEF.FEVectorBlocks)) blocks.") + end + FEF.FEVectorBlocks[i] +end Base.getindex(FEB::FEVectorBlock, i::Int) = FEB.entries[FEB.offset + i] Base.getindex(FEB::FEVectorBlock, i::AbstractArray) = FEB.entries[FEB.offset .+ i] Base.getindex(FEB::FEVectorBlock, ::Colon) = FEB.entries[(FEB.offset + 1):FEB.last_index] @@ -110,6 +121,21 @@ Returns a view of the part of the full FEVector that coressponds to the block. """ Base.view(FEB::FEVectorBlock) = view(FEB.entries, (FEB.offset + 1):FEB.last_index) + +""" +Returns a view of a slice of the FEVectorBlock, specified by local indices `inds` (which can be an integer, a range, or an array of indices). +The indices are relative to the block (i.e., `1` corresponds to the first entry of the block). + +# Arguments +- `FEB::FEVectorBlock`: The FEVectorBlock to view. +- `inds`: Indices relative to the block (e.g., `1:30`, `[2,4,6]`). + +# Returns +- A view into the underlying entries array for the specified slice. +""" +Base.view(FEB::FEVectorBlock, inds) = view(FEB.entries, FEB.offset .+ inds) + + function LinearAlgebra.norm(FEV::FEVector, p::Real = 2) return norm(FEV.entries, p) end @@ -317,3 +343,18 @@ Scalar product between two FEVEctorBlocks. function LinearAlgebra.dot(a::FEVectorBlock{T}, b::FEVectorBlock{T}) where {T} return dot(view(a), view(b)) end + +""" +$(TYPEDSIGNATURES) + +Convert an `FEVector` to a standard Julia array containing all coefficients. +""" +Base.Array(FEV::FEVector) = copy(FEV.entries) + +""" +$(TYPEDSIGNATURES) + +Convert an `FEVectorBlock` to a standard Julia array containing the coefficients for that block. +""" +Base.Array(FEB::FEVectorBlock) = copy(FEB.entries[(FEB.offset + 1):FEB.last_index]) +`````` diff --git a/src/interpolations.jl b/src/interpolations.jl index e6b193b..3b025a1 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -432,7 +432,6 @@ piecewise_nodevalues!( source_offset::Int = 0, zero_target::Bool = true, continuous::Bool = false - ) ```` Evaluate a finite element function (given by the coefficient vector `source` and FE space `FE`) at all nodes, but store the results in a piecewise (cellwise) fashion, i.e., for each cell, the values at its local nodes are written into `target`. The result is organized so that each column of `target` corresponds to a cell, and each row corresponds to the values at the cell's nodes. @@ -910,3 +909,17 @@ function displace_mesh(xgrid::ExtendableGrid, source::FEVectorBlock; kwargs...) displace_mesh!(xgrid_displaced, source; kwargs...) return xgrid_displaced end + +""" +Convenience method: Evaluate a finite element function (given by the FEVectorBlock), applying an optional function operator, and write the results into `target`. + +This forwards to the main nodevalues! method using a view of the block's entries and its FESpace. +""" +function nodevalues!( + target::AbstractArray{T, 2}, + block::FEVectorBlock{T, Tv, Ti, FEType, AT}, + operator::Type{<:AbstractFunctionOperator} = Identity; + kwargs... + ) where {T, Tv, Ti, FEType, AT} + return nodevalues!(target, view(block), block.FES, operator; kwargs...) +end From b701db71999f5391f163297f82336f38bbbbb6e3 Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Wed, 25 Jun 2025 09:12:11 +0200 Subject: [PATCH 5/8] remove ambuigity with view function from ExtendableGrids.subgrid.jl:331 --- src/fevector.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/fevector.jl b/src/fevector.jl index ab8f7b0..6469037 100644 --- a/src/fevector.jl +++ b/src/fevector.jl @@ -115,7 +115,7 @@ Base.last(FEB::FEVectorBlock) = FEB.last_index """ -$(TYPEDEF) +$(TYPEDSIGNATURES) Returns a view of the part of the full FEVector that coressponds to the block. """ @@ -123,6 +123,8 @@ Base.view(FEB::FEVectorBlock) = view(FEB.entries, (FEB.offset + 1):FEB.last_inde """ +$(TYPEDSIGNATURES) + Returns a view of a slice of the FEVectorBlock, specified by local indices `inds` (which can be an integer, a range, or an array of indices). The indices are relative to the block (i.e., `1` corresponds to the first entry of the block). @@ -133,7 +135,7 @@ The indices are relative to the block (i.e., `1` corresponds to the first entry # Returns - A view into the underlying entries array for the specified slice. """ -Base.view(FEB::FEVectorBlock, inds) = view(FEB.entries, FEB.offset .+ inds) +Base.view(FEB::FEVectorBlock, inds::Union{Integer, AbstractArray{<:Integer}, UnitRange{<:Integer}}) = view(FEB.entries, FEB.offset .+ inds) function LinearAlgebra.norm(FEV::FEVector, p::Real = 2) @@ -145,7 +147,7 @@ function LinearAlgebra.norm(FEV::FEVectorBlock, p::Real = 2) end """ -$(TYPEDEF) +$(TYPEDSIGNATURES) Returns a vector with the individual norms of all blocks. """ From dfdeaff81b4e947d6208f786b4c5d7b99eff8d6b Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Wed, 25 Jun 2025 14:21:17 +0200 Subject: [PATCH 6/8] improved docstrings in quadrature, updated changelog --- CHANGELOG.md | 8 ++++++-- src/quadrature.jl | 31 +++++++++++++++++++++++++++---- 2 files changed, 33 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9dca4c5..229d348 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,7 +1,11 @@ # CHANGES -## v1.2.0 June 24, 2025 - - major documentation overhaul +## v1.2.0 June 25, 2025 + - major documentation and docstring overhaul + - improved show functions and constructors for FEMatrix, FEVector + - added range check in getindex functions for FEVector + - added slice views for FEVectorBlocks + - added Array conversions for FEVector and FEVectorBlock ## v1.1.1 June 16, 2025 - docstring improvements diff --git a/src/quadrature.jl b/src/quadrature.jl index 9bf58e3..f64c720 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -19,7 +19,21 @@ abstract type QuadratureRule{T <: Real, ET <: AbstractElementGeometry} end """ $(TYPEDEF) -A struct that contains the name of the quadrature rule, the reference points and the weights for the parameter-determined element geometry. +A concrete quadrature rule for a given element geometry and number type. + +It represents a set of quadrature (integration) points and weights for a specific reference element geometry (such as an interval, triangle, quadrilateral, tetrahedron, or hexahedron) and number type. + +# Fields +- `name::String`: A descriptive name for the quadrature rule (e.g., "midpoint rule", "Gauss rule order 3"). +- `xref::Vector{Vector{T}}`: Reference coordinates of the quadrature points, given as a vector of coordinate vectors (one per point, each of length `dim`). +- `w::Vector{T}`: Weights associated with each quadrature point, typically summing to the measure of the reference element. + +# Type Parameters +- `T <: Real`: Number type for coordinates and weights (e.g., `Float64`, `Rational{Int}`). +- `ET <: AbstractElementGeometry`: The reference element geometry type (e.g., `Edge1D`, `Triangle2D`). +- `dim`: The topological dimension of the element geometry. +- `npoints`: The number of quadrature points. + """ struct SQuadratureRule{T <: Real, ET <: AbstractElementGeometry, dim, npoints} <: QuadratureRule{T, ET} name::String @@ -61,9 +75,18 @@ end """ $(TYPEDSIGNATURES) -sets up a quadrature rule that evaluates at vertices of element geometry; -not optimal from quadrature point of view, but helpful when interpolating. -Note, that order of xref matches dof order of H1Pk element +Constructs a quadrature rule that evaluates at the vertices of the reference element geometry `ET`. + +This rule is not optimal for numerical integration, but is especially useful for nodal interpolation, visualization, and extracting nodal values in finite element computations. The order parameter determines the inclusion of higher-order nodes (e.g., edge, face or cell nodes for higher-order Lagrange elements). + +# Arguments +- `ET::Type{<:AbstractElementGeometry}`: The reference element geometry (e.g., `Edge1D`, `Triangle2D`, `Parallelogram2D`, `Tetrahedron3D`, `Parallelepiped3D`). +- `order::Integer`: Polynomial order of the finite element (default: `1`). Higher orders include additional points corresponding to edge, face, or cell dofs. +- `T`: Number type for the coordinates and weights (default: `Float64`). + +# Returns +- A quadrature rule containing the nodal points (`xref`) and equal weights (`w`), matching the dof structure of the corresponding Lagrange element. + """ function VertexRule(ET::Type{Edge1D}, order = 1; T = Float64) if order == 0 From d9b5053e04d7322294be9ef9221b0e7802c3b55e Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Wed, 25 Jun 2025 18:58:41 +0200 Subject: [PATCH 7/8] allow iterate through blocks of FEVector --- src/fevector.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/fevector.jl b/src/fevector.jl index 6469037..33fead0 100644 --- a/src/fevector.jl +++ b/src/fevector.jl @@ -112,6 +112,8 @@ Base.size(FEF::FEVector) = size(FEF.FEVectorBlocks) Base.size(FEB::FEVectorBlock) = FEB.last_index - FEB.offset Base.first(FEB::FEVectorBlock) = FEB.offset + 1 Base.last(FEB::FEVectorBlock) = FEB.last_index +Base.iterate(FEV::FEVector) = iterate(FEV.FEVectorBlocks) +Base.iterate(FEV::FEVector, state) = iterate(FEV.FEVectorBlocks, state) """ From ae47f2a4d3383822ce1f120160914a814ed7acca Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Thu, 26 Jun 2025 12:48:45 +0200 Subject: [PATCH 8/8] fixed source_offset issue --- src/interpolations.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/interpolations.jl b/src/interpolations.jl index 3b025a1..007c978 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -617,15 +617,15 @@ function nodevalues( if nodes == [] if cellwise target = zeros(T, nvals * max_num_targets_per_source(source.FES.dofgrid[CellNodes]), num_cells(source.FES.dofgrid)) - piecewise_nodevalues!(target, source.entries, source.FES, operator; abs = abs, kwargs...) + piecewise_nodevalues!(target, source.entries, source.FES, operator; abs = abs, source_offset = source.offset, kwargs...) else target = zeros(T, nvals, num_nodes(source.FES.dofgrid)) - nodevalues!(target, source.entries, source.FES, operator; continuous = continuous, abs = abs, kwargs...) + nodevalues!(target, source.entries, source.FES, operator; continuous = continuous, source_offset = source.offset, abs = abs, kwargs...) end else @assert cellwise == false target = zeros(T, nvals, length(nodes)) - nodevalues_subset!(target, source.entries, source.FES, operator; continuous = continuous, abs = abs, nodes = nodes, kwargs...) + nodevalues_subset!(target, source.entries, source.FES, operator; continuous = continuous, source_offset = source.offset, abs = abs, nodes = nodes, kwargs...) end return target end