WIP: New vector interpolations and mappings

docs/Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.11.1"
 manifest_format = "2.0"
-project_hash = "4e52f4aa4cee9f66ec4f633f0ae538fbd227ac5e"
+project_hash = "c8cd4f9010376e633064839ca68de8e4b622021f"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "eea5d80188827b35333801ef97a40c2ed653b081"
@@ -69,6 +69,18 @@ git-tree-sha1 = "d57bd3762d308bded22c3b82d033bff85f6195c6"
 uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
 version = "0.4.0"
 
+[[deps.Arpack]]
+deps = ["Arpack_jll", "Libdl", "LinearAlgebra", "Logging"]
+git-tree-sha1 = "9b9b347613394885fd1c8c7729bfc60528faa436"
+uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+version = "0.5.4"
+
+[[deps.Arpack_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "OpenBLAS_jll", "Pkg"]
+git-tree-sha1 = "5ba6c757e8feccf03a1554dfaf3e26b3cfc7fd5e"
+uuid = "68821587-b530-5797-8361-c406ea357684"
+version = "3.5.1+1"
+
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra"]
 git-tree-sha1 = "3640d077b6dafd64ceb8fd5c1ec76f7ca53bcf76"
@@ -931,6 +943,16 @@ git-tree-sha1 = "267dad6b4b7b5d529c76d40ff48d33f7e94cb834"
 uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 version = "0.9.6"
 
+[[deps.KrylovKit]]
+deps = ["GPUArraysCore", "LinearAlgebra", "PackageExtensionCompat", "Printf", "Random", "VectorInterface"]
+git-tree-sha1 = "3ba86a12066c19a6ed0bc963984bcf89d495133e"
+uuid = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
+version = "0.8.2"
+weakdeps = ["ChainRulesCore"]
+
+    [deps.KrylovKit.extensions]
+    KrylovKitChainRulesCoreExt = "ChainRulesCore"
+
 [[deps.LAME_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "170b660facf5df5de098d866564877e119141cbd"
@@ -2325,6 +2347,12 @@ git-tree-sha1 = "c2d0db3ef09f1942d08ea455a9e252594be5f3b6"
 uuid = "4004b06d-e244-455f-a6ce-a5f9919cc534"
 version = "1.0.1"
 
+[[deps.VectorInterface]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "cea8abaa6e43f72f97a09cf95b80c9eb53ff75cf"
+uuid = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
+version = "0.4.9"
+
 [[deps.VectorizationBase]]
 deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static", "StaticArrayInterface"]
 git-tree-sha1 = "e7f5b81c65eb858bed630fe006837b935518aca5"

docs/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 Changelog = "5217a498-cd5d-4ec6-b8c2-9b85a09b6e3e"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
@@ -9,6 +10,7 @@ FerriteMeshParser = "0f8c756f-80dd-4a75-85c6-b0a5ab9d4620"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"

docs/liveserver.jl

@@ -14,7 +14,7 @@ push!(ARGS, "liveserver")
 # Run LiveServer.servedocs(...)
 import LiveServer
 LiveServer.servedocs(;
-    host = "0.0.0.0",
+    # host = "0.0.0.0",
     # Documentation root where make.jl and src/ are located
     foldername = joinpath(repo_root, "docs"),
     # Extra source folder to watch for changes

docs/make.jl

@@ -55,6 +55,7 @@ bibtex_plugin = CitationBibliography(
         "Tutorials" => [
             "Tutorials overview" => "tutorials/index.md",
             "tutorials/heat_equation.md",
+            "tutorials/heat_equation_hdiv.md",
             "tutorials/linear_elasticity.md",
             "tutorials/incompressible_elasticity.md",
             "tutorials/hyperelasticity.md",
@@ -67,6 +68,7 @@ bibtex_plugin = CitationBibliography(
             "tutorials/reactive_surface.md",
             "tutorials/linear_shell.md",
             "tutorials/dg_heat_equation.md",
+            "tutorials/maxwell.md",
         ],
         "Topic guides" => [
             "Topic guide overview" => "topics/index.md",

docs/src/assets/references.bib

@@ -190,3 +190,25 @@ @article{CroRav:1973:cnf
   year={1973},
   publisher={EDP Sciences}
 }
+@book{Gatica2014,
+title = {A Simple Introduction to the Mixed Finite Element Method: Theory and Applications},
+ISBN = {9783319036953},
+ISSN = {2191-8201},
+url = {http://dx.doi.org/10.1007/978-3-319-03695-3},
+DOI = {10.1007/978-3-319-03695-3},
+journal = {SpringerBriefs in Mathematics},
+publisher = {Springer International Publishing},
+author = {Gatica,  Gabriel N.},
+year = {2014}
+}
+@book{Boffi2013,
+  title = {Mixed Finite Element Methods and Applications},
+  ISBN = {9783642365195},
+  ISSN = {0179-3632},
+  url = {http://dx.doi.org/10.1007/978-3-642-36519-5},
+  DOI = {10.1007/978-3-642-36519-5},
+  journal = {Springer Series in Computational Mathematics},
+  publisher = {Springer Berlin Heidelberg},
+  author = {Boffi,  Daniele and Brezzi,  Franco and Fortin,  Michel},
+  year = {2013}
+}

docs/src/literate-tutorials/heat_equation_hdiv.jl

@@ -0,0 +1,234 @@
+# # [Heat equation - Mixed H(div) conforming formulation)](@id tutorial-heat-equation-hdiv)
+# As an alternative to the standard formulation for solving the heat equation used in
+# the [heat equation tutorial](@ref tutorial-heat-equation), we can used a mixed formulation
+# where both the temperature, $u(\mathbf{x})$, and the heat flux, $\boldsymbol{q}(\boldsymbol{x})$,
+# are primary variables. From a theoretical standpoint, there are many details on e.g. which combinations
+# of interpolations that are stable. See e.g. [Gatica2014](@cite) and [Boffi2013](@cite) for further reading.
+# This tutorial is based on the theory in
+# [Fenics' mixed poisson example](https://fenicsproject.org/olddocs/dolfin/1.4.0/python/demo/documented/mixed-poisson/python/documentation.html).
+#
+# ![Temperature solution](https://raw.githubusercontent.com/Ferrite-FEM/Ferrite.jl/refs/heads/gh-pages/assets/heat_equation_hdiv.png)
+# **Figure:** Temperature distribution considering a central part with lower heat conductivity.
+#
+# The advantage with the mixed formulation is that the heat flux is approximated better. However, the
+# temperature becomes discontinuous where the conductivity is discontinuous.
+#
+# ## Theory
+# We start with the strong form of the heat equation: Find the temperature, $u(\boldsymbol{x})$, and heat flux, $\boldsymbol{q}(x)$,
+# such that
+# ```math
+# \begin{align*}
+# \boldsymbol{\nabla}\cdot \boldsymbol{q} &= h(\boldsymbol{x}), \quad \forall \boldsymbol{x} \in \Omega \\
+# \boldsymbol{q}(\boldsymbol{x}) &= - k\ \boldsymbol{\nabla} u(\boldsymbol{x}), \quad \forall \boldsymbol{x} \in \Omega \\
+# \boldsymbol{q}(\boldsymbol{x})\cdot \boldsymbol{n}(\boldsymbol{x}) &= q_n, \quad \forall \boldsymbol{x} \in \Gamma_\mathrm{N}\\
+# u(\boldsymbol{x}) &= u_\mathrm{D}, \quad \forall \boldsymbol{x} \in \Gamma_\mathrm{D}
+# \end{align*}
+# ```
+#
+# From this strong form, we can formulate the weak form as a mixed formulation.
+# Find $u \in \mathbb{U}$ and $\boldsymbol{q}\in\mathbb{Q}$ such that
+# ```math
+# \begin{align*}
+# \int_{\Omega} \delta u [\boldsymbol{\nabla} \cdot \boldsymbol{q}]\ \mathrm{d}\Omega &= \int_\Omega \delta u h\ \mathrm{d}\Omega, \quad \forall\ \delta u \in \delta\mathbb{U} \\
+# \int_{\Omega} \boldsymbol{\delta q} \cdot \boldsymbol{q}\ \mathrm{d}\Omega &= -\int_\Omega \boldsymbol{\delta q} \cdot [k\ \boldsymbol{\nabla} u]\ \mathrm{d}\Omega \\
+# \int_{\Omega} \boldsymbol{\delta q} \cdot \boldsymbol{q}\ \mathrm{d}\Omega - \int_{\Omega} [\boldsymbol{\nabla} \cdot \boldsymbol{\delta q}] k u \ \mathrm{d}\Omega &=
+# -\int_\Gamma \boldsymbol{\delta q} \cdot \boldsymbol{n} k\ u\ \mathrm{d}\Omega, \quad \forall\ \boldsymbol{\delta q} \in \delta\mathbb{Q}
+# \end{align*}
+# ```
+# where we have the function spaces,
+# ```math
+# \begin{align*}
+# \mathbb{U} &= \delta\mathbb{U} = L^2 \\
+# \mathbb{Q} &= \lbrace \boldsymbol{q} \in H(\mathrm{div}) \text{ such that } \boldsymbol{q}\cdot\boldsymbol{n} = q_\mathrm{n} \text{ on } \Gamma_\mathrm{D}\rbrace \\
+# \delta\mathbb{Q} &= \lbrace \boldsymbol{q} \in H(\mathrm{div}) \text{ such that } \boldsymbol{q}\cdot\boldsymbol{n} = 0 \text{ on } \Gamma_\mathrm{D}\rbrace
+# \end{align*}
+# ```
+# A stable choice of finite element spaces for this problem on grid with triangles is using
+# * `DiscontinuousLagrange{RefTriangle, k-1}` for approximating $L^2$
+# * `BrezziDouglasMarini{RefTriangle, k}` for approximating $H(\mathrm{div})$
+#
+# We will also investigate the consequences of using $H^1$ `Lagrange` instead of $H(\mathrm{div})$ interpolations.
+#
+# ## Commented Program
+#
+# Now we solve the problem in Ferrite. What follows is a program spliced with comments.
+#
+# First we load Ferrite,
+using Ferrite
+# And define our grid, representing a channel with a central part having a lower
+# conductivity, $k$, than the surrounding.
+function create_grid(ny::Int)
+    width = 10.0
+    length = 40.0
+    center_width = 5.0
+    center_length = 20.0
+    upper_right = Vec((length / 2, width / 2))
+    grid = generate_grid(Triangle, (round(Int, ny * length / width), ny), -upper_right, upper_right)
+    addcellset!(grid, "center", x -> abs(x[1]) < center_length / 2 && abs(x[2]) < center_width / 2)
+    addcellset!(grid, "around", setdiff(1:getncells(grid), getcellset(grid, "center")))
+    return grid
+end
+
+grid = create_grid(10)
+
+# ### Setup
+# We define one `CellValues` for each field which share the same quadrature rule.
+ip_geo = geometric_interpolation(getcelltype(grid))
+ipu = DiscontinuousLagrange{RefTriangle, 0}()
+ipq = Ferrite.BrezziDouglasMarini{2, RefTriangle, 1}()
+qr = QuadratureRule{RefTriangle}(2)
+cellvalues = (u = CellValues(qr, ipu, ip_geo), q = CellValues(qr, ipq, ip_geo))
+
+# Distribute the degrees of freedom
+dh = DofHandler(grid)
+add!(dh, :u, ipu)
+add!(dh, :q, ipq)
+close!(dh);
+
+# In this problem, we have a zero temperature on the boundary, Γ, which is enforced
+# via zero Neumann boundary conditions. Hence, we don't need a `Constrainthandler`.
+Γ = union((getfacetset(grid, name) for name in ("left", "right", "bottom", "top"))...)
+
+# ### Element implementation
+
+function assemble_element!(Ke::Matrix, fe::Vector, cv::NamedTuple, dr::NamedTuple, k::Number)
+    cvu = cv[:u]
+    cvq = cv[:q]
+    dru = dr[:u]
+    drq = dr[:q]
+    h = 1.0 # Heat source
+    ## Loop over quadrature points
+    for q_point in 1:getnquadpoints(cvu)
+        ## Get the quadrature weight
+        dΩ = getdetJdV(cvu, q_point)
+        ## Loop over test shape functions
+        for (iu, Iu) in pairs(dru)
+            δNu = shape_value(cvu, q_point, iu)
+            ## Add contribution to fe
+            fe[Iu] += δNu * h * dΩ
+            ## Loop over trial shape functions
+            for (jq, Jq) in pairs(drq)
+                div_Nq = shape_divergence(cvq, q_point, jq)
+                ## Add contribution to Ke
+                Ke[Iu, Jq] += (δNu * div_Nq) * dΩ
+            end
+        end
+        for (iq, Iq) in pairs(drq)
+            δNq = shape_value(cvq, q_point, iq)
+            div_δNq = shape_divergence(cvq, q_point, iq)
+            for (ju, Ju) in pairs(dru)
+                Nu = shape_value(cvu, q_point, ju)
+                Ke[Iq, Ju] -= div_δNq * k * Nu * dΩ
+            end
+            for (jq, Jq) in pairs(drq)
+                Nq = shape_value(cvq, q_point, jq)
+                Ke[Iq, Jq] += (δNq ⋅ Nq) * dΩ
+            end
+        end
+    end
+    return Ke, fe
+end
+#md nothing # hide
+
+# ### Global assembly
+
+function assemble_global(cellvalues, dh::DofHandler)
+    grid = dh.grid
+    ## Allocate the element stiffness matrix and element force vector
+    dofranges = (u = dof_range(dh, :u), q = dof_range(dh, :q))
+    ncelldofs = ndofs_per_cell(dh)
+    Ke = zeros(ncelldofs, ncelldofs)
+    fe = zeros(ncelldofs)
+    ## Allocate global system matrix and vector
+    K = allocate_matrix(dh)
+    f = zeros(ndofs(dh))
+    ## Create an assembler
+    assembler = start_assemble(K, f)
+    x = copy(getcoordinates(grid, 1))
+    dofs = copy(celldofs(dh, 1))
+    ## Loop over all cells
+    for (cells, k) in (
+            (getcellset(grid, "center"), 0.1),
+            (getcellset(grid, "around"), 1.0),
+        )
+        for cellnr in cells
+            ## Reinitialize cellvalues for this cell
+            cell = getcells(grid, cellnr)
+            getcoordinates!(x, grid, cell)
+            celldofs!(dofs, dh, cellnr)
+            reinit!(cellvalues[:u], cell, x)
+            reinit!(cellvalues[:q], cell, x)
+            ## Reset to 0
+            fill!(Ke, 0)
+            fill!(fe, 0)
+            ## Compute element contribution
+            assemble_element!(Ke, fe, cellvalues, dofranges, k)
+            ## Assemble Ke and fe into K and f
+            assemble!(assembler, dofs, Ke, fe)
+        end
+    end
+    return K, f
+end
+#md nothing # hide
+
+# ### Solution of the system
+K, f = assemble_global(cellvalues, dh);
+u = K \ f;
+
+# ### Exporting to VTK
+# Currently, exporting discontinuous interpolations is not supported.
+# Since in this case, we have a single temperature degree of freedom
+# per cell, we work around this by calculating the per-cell temperature.
+temperature_dof = first(dof_range(dh, :u))
+u_cells = map(1:getncells(grid)) do i
+    u[celldofs(dh, i)[temperature_dof]]
+end
+VTKGridFile("heat_equation_hdiv", dh) do vtk
+    write_cell_data(vtk, u_cells, "temperature")
+end
+
+# ## Postprocess the total flux
+# We applied a constant unit heat source to the body, and the
+# total heat flux exiting across the boundary should therefore
+# match the area for the considered stationary case.
+function calculate_flux(dh, boundary_facets, ip, a)
+    grid = dh.grid
+    qr = FacetQuadratureRule{RefTriangle}(4)
+    ip_geo = geometric_interpolation(getcelltype(grid))
+    fv = FacetValues(qr, ip, ip_geo)
+
+    dofrange = dof_range(dh, :q)
+    flux = 0.0
+    dofs = celldofs(dh, 1)
+    ae = zeros(length(dofs))
+    x = getcoordinates(grid, 1)
+    for (cellnr, facetnr) in boundary_facets
+        getcoordinates!(x, grid, cellnr)
+        cell = getcells(grid, cellnr)
+        celldofs!(dofs, dh, cellnr)
+        map!(i -> a[i], ae, dofs)
+        reinit!(fv, cell, x, facetnr)
+        for q_point in 1:getnquadpoints(fv)
+            dΓ = getdetJdV(fv, q_point)
+            n = getnormal(fv, q_point)
+            q = function_value(fv, q_point, ae, dofrange)
+            flux += (q ⋅ n) * dΓ
+        end
+    end
+    return flux
+end
+
+println("Outward flux: ", calculate_flux(dh, Γ, ipq, u))
+
+# Note that this is not the case for the standard [Heat equation](@ref tutorial-heat-equation),
+# as the flux terms are less accurately approximated. A fine mesh is required to converge in that case.
+# However, the present example gives a worse approximation of the temperature field.
+
+#md # ## [Plain program](@id tutorial-heat-equation-hdiv-plain)
+#md #
+#md # Here follows a version of the program without any comments.
+#md # The file is also available here: [`heat_equation_hdiv.jl`](heat_equation_hdiv.jl).
+#md #
+#md # ```julia
+#md # @__CODE__
+#md # ```