Started adding BrezziDouglasMarini

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/src/literate-tutorials/heat_equation_rt.jl

@@ -1,5 +1,6 @@
 # # [Heat equation (Mixed, RaviartThomas)](@id tutorial-heat-equation-rt)
-# Note, there are a lot to consider here it seems like. Good ref,
+# Note, there are a lot to consider here it seems like. Good refs,
+# ```
 # @book{Gatica2014,
 # title = {A Simple Introduction to the Mixed Finite Element Method: Theory and Applications},
 # ISBN = {9783319036953},
@@ -24,54 +25,63 @@
 #   year = {2013}
 # }
 # for a(n even) more comprehensive book.
-#
-# ## Strong form
-# ```math
-# \nabla \cdot \boldsymbol{q} = h \in \Omega \\
-# \boldsymbol{q} = - k\ \nabla u \in \Omega \\
-# \boldsymbol{q}\cdot \boldsymbol{n} = q_n \in \Gamma_\mathrm{N}\\
-# u = u_\mathrm{D} \in \Gamma_\mathrm{D}
-# ```
-#
-# ## Weak form
-# ### Part 1
-# ```math
-# \int_{\Omega} \delta u \nabla \cdot \boldsymbol{q}\ \mathrm{d}\Omega = \int_{\Omega} \delta u\ h\ \mathrm{d}\Omega \\
-# \int_{\Gamma} \delta u \boldsymbol{n} \cdot \boldsymbol{q}\ \mathrm{d}\Gamma -
-# \int_{\Omega} \nabla (\delta u) \cdot \boldsymbol{q}\ \mathrm{d}\Omega = \int_{\Omega} \delta u\ h\ \mathrm{d}\Omega \\
 # ```
 #
-# ### Part 2
+# ## 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
-# \int_{\Omega} \boldsymbol{\delta q} \cdot \boldsymbol{q}\ \mathrm{d}\Omega = - \int_{\Omega} \boldsymbol{\delta q} \cdot \left[k\ \nabla u\right]\ \mathrm{d}\Omega
+# \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*}
 # ```
-# where no Green-Gauss theorem is applied.
 #
-# ### Summary
-# The weak form becomes, find $u\in H^1$ and $\boldsymbol{q} \in H\mathrm{(div)}$, such that
+# 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} \nabla (\delta u) \cdot \boldsymbol{q}\ \mathrm{d}\Omega &= \int_{\Omega} \delta u\ h\ \mathrm{d}\Omega -
-# \int_{\Gamma} \delta u\ q_\mathrm{n}\ \mathrm{d}\Gamma
-# \quad
-# \forall\ \delta u \in \delta H^1 \\
-# \int_{\Omega} \boldsymbol{\delta q} \cdot \boldsymbol{q}\ \mathrm{d}\Omega &= - \int_{\Omega} \boldsymbol{\delta q} \cdot \left[k\ \nabla u\right]\ \mathrm{d}\Omega
-#  \quad \forall\ \boldsymbol{\delta q} \in \delta H\mathrm{(div)}
+# \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
+# * $\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$
+# * $\mathbb{Q} = \lbrace \boldsymbol{q} \in H(\mathrm{div}) \text{such that} \boldsymbol{q}\cdot\boldsymbol{n} = 0 \text{ on } \Gamma_\mathrm{D}\rbrace$
+#
+# 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})$
+# following [fenics](https://fenicsproject.org/olddocs/dolfin/1.4.0/python/demo/documented/mixed-poisson/python/documentation.html).
+# For further details, see Boffi2013.
+# We will also see what happens if we instead use `Lagrange` elements which are a subspace of $H^1$ instead of $H(\mathrm{div})$ elements.
 #
 # ## Commented Program
 #
 # Now we solve the problem in Ferrite. What follows is a program spliced with comments.
-#md # The full program, without comments, can be found in the next [section](@ref heat_equation-plain-program).
 #
-# First we load Ferrite, and some other packages we need
-using Ferrite, SparseArrays
-# We start by generating a simple grid with 20x20 quadrilateral elements
-# using `generate_grid`. The generator defaults to the unit square,
-# so we don't need to specify the corners of the domain.
-#grid = generate_grid(QuadraticTriangle, (20, 20));
-grid = generate_grid(Triangle, (20, 20));
+# 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_width/2 && abs(x[2]) < center_length / 2)
+    addcellset!(grid, "around", setdiff(1:getncells(grid), getcellset(grid, "center")))
+    return grid
+end
+
+grid = create_grid(10)
 
 # ### Trial and test functions
 # A `CellValues` facilitates the process of evaluating values and gradients of
@@ -82,7 +92,7 @@ grid = generate_grid(Triangle, (20, 20));
 # the same reference element. We combine the interpolation and the quadrature rule
 # to a `CellValues` object.
 ip_geo = geometric_interpolation(getcelltype(grid))
-ipu = Lagrange{RefTriangle, 1}() # Why does it "explode" for 2nd order ipu?
+ipu = DiscontinuousLagrange{RefTriangle, 1}() # Why does it "explode" for 2nd order ipu?
 ipq = RaviartThomas{2,RefTriangle, 1}()
 qr = QuadratureRule{RefTriangle}(2)
 cellvalues = (u=CellValues(qr, ipu, ip_geo), q=CellValues(qr, ipq, ip_geo))

docs/src/literate-tutorials/heat_equation_triangle.jl

@@ -32,7 +32,6 @@
 # ## Commented Program
 #
 # Now we solve the problem in Ferrite. What follows is a program spliced with comments.
-#md # The full program, without comments, can be found in the next [section](@ref heat_equation-plain-program).
 #
 # First we load Ferrite, and some other packages we need
 using Ferrite, SparseArrays

src/interpolations.jl

@@ -1855,13 +1855,64 @@ function get_direction(::RaviartThomas{2, RefTriangle, 2}, j, cell)
     return ifelse(edge[2] > edge[1], 1, -1)
 end
 
+#####################################
+# Brezzi-Douglas–Marini, H(div)     #
+#####################################
+struct BrezziDouglasMarini{vdim, shape, order} <: VectorInterpolation{vdim, shape, order} end
+mapping_type(::BrezziDouglasMarini) = ContravariantPiolaMapping()
+reference_coordinates(ip::BrezziDouglasMarini{vdim}) where vdim = fill(NaN*zero(Vec{vdim}), getnbasefunctions(ip))
+dirichlet_facedof_indices(ip::BrezziDouglasMarini) = facetdof_interior_indices(ip)
+n_dbc_components(::BrezziDouglasMarini) = 1
+#=
+----------------+--------------------
+Vertex numbers: | Vertex coordinates:
+    2           |
+    | \         | v1: 𝛏 = (1.0, 0.0)
+    |   \       | v2: 𝛏 = (0.0, 1.0)
+ξ₂^ |     \     | v3: 𝛏 = (0.0, 0.0)
+  | 3-------1   |
+  +--> ξ₁       |
+----------------+--------------------
+Edge numbers:   | Edge identifiers:
+    +           |
+    | \         | e1: (v1, v2)
+    2   1       | e2: (v2, v3)
+    |     \     | e3: (v3, v1)
+    +---3---+   |
+----------------+--------------------
+=#
+
+# RefTriangle, 1st order Lagrange
+function reference_shape_value(ip::BrezziDouglasMarini{2, RefTriangle, 1}, ξ::Vec{2}, i::Int)
+    x, y = ξ
+    # Edge 1
+    i == 1 && return Vec(4x, -2y) # Changed sign to make positive outwards
+    i == 2 && return Vec(-2x, 4y) # Changed sign to make positive outwards
+    # Edge 2 (reverse order to follow Ferrite convention)
+    i == 3 && return Vec(-2x - 6y + 2, 4y)
+    i == 4 && return Vec(4x + 6y - 4, -2y)
+    # Edge 3
+    i == 5 && return Vec(-2x, 6x + 4y - 4) # Changed sign to make positive outwards
+    i == 6 && return Vec(4x, -6x - 2y + 2) # Changed sign to make positive outwards
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::BrezziDouglasMarini{2, RefTriangle, 1}) = 6
+edgedof_interior_indices(::BrezziDouglasMarini{2, RefTriangle, 1}) = ((1, 2), (3, 4), (5, 6))
+adjust_dofs_during_distribution(::BrezziDouglasMarini{2, RefTriangle, 1}) = false
+
+function get_direction(::BrezziDouglasMarini{2, RefTriangle, 1}, j, cell)
+    edge = edges(cell)[(j + 1) ÷ 2]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end
+
 #####################################
 # Nedelec (1st kind), H(curl)       #
 #####################################
 struct Nedelec{vdim, shape, order} <: VectorInterpolation{vdim, shape, order} end
 mapping_type(::Nedelec) = CovariantPiolaMapping()
 reference_coordinates(ip::Nedelec{vdim}) where vdim = fill(NaN*zero(Vec{vdim}), getnbasefunctions(ip))
-dirichlet_facedof_indices(ip::Nedelec) = facedof_interior_indices(ip)
+dirichlet_facedof_indices(ip::Nedelec) = facetdof_interior_indices(ip)
 n_dbc_components(::Nedelec) = 1
 # RefTriangle, 1st order Lagrange
 # https://defelement.com/elements/examples/triangle-nedelec1-lagrange-1.html

test/test_interpolations.jl

@@ -270,7 +270,7 @@ end
 #    Where f(s) = 1 for linear interpolation and f(s)=1-s and f(s)=s for 2nd order interpolation (first and second shape function)
 #    And s is the path parameter ∈[0,1] along the positive direction of the path.
 # 2) Zero along other edges: Nⱼ ⋅ v = 0 if j∉𝔇
-@testset "Nedelec" begin
+@testset "H(curl) on RefCell" begin
     lineqr = QuadratureRule{RefLine}(20)
     for ip in (Nedelec{2,RefTriangle,1}(), Nedelec{2,RefTriangle,2}(), Nedelec{3,RefTetrahedron,1}(), Nedelec{3,RefHexahedron,1}())
         cell = reference_cell(getrefshape(ip))
@@ -302,6 +302,46 @@ end
     end
 end
 
+# Required properties of shape value Nⱼ of an edge-elements (Hdiv) on an edge with normal n, length L, and dofs ∈ 𝔇
+# 1) Unit property: ∫(Nⱼ ⋅ n f(s) dS) = 1
+#    Where f(s) = 1 for single shape function on edge, and f(s)=1-s and f(s)=s for two shape functions on edge
+#    s is the path parameter ∈[0,1] along the positive direction of the path.
+# 2) Zero normal component on other edges: Nⱼ ⋅ n = 0 if j∉𝔇
+@testset "H(div) on RefCell" begin
+    lineqr = QuadratureRule{RefLine}(20)
+    for ip in (RaviartThomas{2, RefTriangle, 1}(), Ferrite.BrezziDouglasMarini{2, RefTriangle, 1}(), )
+        cell = reference_cell(getrefshape(ip))
+        cell_facets = Ferrite.facets(cell)
+        dofs = Ferrite.facetdof_interior_indices(ip)
+        x = Ferrite.reference_coordinates(geometric_interpolation(typeof(cell)))
+        normals = reference_normals(geometric_interpolation(typeof(cell)))
+        @testset "$ip" begin
+            for (facet_nr, (i1, i2)) in enumerate(cell_facets)
+                @testset "Facet $facet_nr" begin
+                    Δx = x[i2]-x[i1]
+                    x0 = (x[i1]+x[i2])/2
+                    L = norm(Δx)
+                    n = normals[facet_nr]
+                    for (idof, shape_nr) in enumerate(dofs[facet_nr])
+                        nfacetdofs = length(dofs[facet_nr])
+                        f(x) = nfacetdofs == 1 ? 1.0 : (idof == 1 ? 1-x : x)
+                        s = line_integral(lineqr, ip, shape_nr, x0, Δx, L, n, f)
+                        @test s ≈ one(s)
+                    end
+                    for (j_facet, shape_nrs) in enumerate(dofs)
+                        j_facet == facet_nr && continue
+                        for shape_nr in shape_nrs
+                            for ξ in (x[i1] + r*Δx for r in [0.0, rand(3)..., 1.0])
+                                @test abs(reference_shape_value(ip, ξ, shape_nr) ⋅ n) < eps()*100
+                            end
+                        end
+                    end
+                end
+            end
+        end
+    end
+end
+
 tupleshift(t::NTuple{N}, shift::Int) where N = ntuple(i -> t[mod(i - 1 - shift, N) + 1], N)
 #tupleshift(t::NTuple, shift::Int) = tuple(circshift(SVector(t), shift)...)
 cell_permutations(cell::Quadrilateral)      = (Quadrilateral(tupleshift(cell.nodes, shift)) for shift in 0:3)
@@ -336,9 +376,10 @@ end
     include(joinpath(@__DIR__, "InterpolationTestUtils.jl"))
     import .InterpolationTestUtils as ITU
     nel = 3
-    transformation_functions = Dict(
-        Nedelec=>(v,n)-> v - n*(v⋅n),   # Hcurl (tangent continuity)
-        RaviartThomas=>(v,n) -> v ⋅ n)  # Hdiv (normal continuity)
+    hdiv_check(v, n) = v ⋅ n        # Hdiv (normal continuity)
+    hcurl_check(v, n) = v - n*(v⋅n) # Hcurl (tangent continuity)
+    transformation_functions = ((Nedelec, hcurl_check), (RaviartThomas, hdiv_check), (Ferrite.BrezziDouglasMarini, hdiv_check))
+
     for CT in (Triangle, QuadraticTriangle, Tetrahedron, Hexahedron)
         dim = Ferrite.getrefdim(CT) # dim = sdim = rdim
         p1, p2 = (rand(Vec{dim}), ones(Vec{dim})+rand(Vec{dim}))
@@ -352,11 +393,11 @@ end
         basecell = getcells(grid, cellnr)
         RefShape = Ferrite.getrefshape(basecell)
         for order in (1, 2)
-            for IPT in (Nedelec, RaviartThomas)
+            for (IPT, transformation_function) in transformation_functions
                 dim == 3 && order > 1 && continue
                 IPT == RaviartThomas && (dim == 3 || order > 1) && continue
                 IPT == RaviartThomas && (RefShape == RefHexahedron) && continue
-                transformation_function=transformation_functions[IPT]
+                IPT == Ferrite.BrezziDouglasMarini && !(RefShape == RefTriangle && order == 1) && continue
                 ip = IPT{dim, RefShape, order}()
                 @testset "$CT, $ip" begin
                     for testcell in cell_permutations(basecell)

visualization/2d_vector_interpolation_checkplots.jl

@@ -0,0 +1,23 @@
+using Ferrite
+import CairoMakie as Plt
+
+function plot_shape_function(ip::VectorInterpolation{2, RefShape}, qr::Int, i::Int) where {RefShape<:Ferrite.AbstractRefShape{2}}
+    return plot_shape_function(ip, QuadratureRule{RefShape}(qr), i)
+end
+
+function plot_shape_function(ip::VectorInterpolation{2, RefShape}, qr::QuadratureRule{RefShape}, i::Int) where {RefShape<:Ferrite.AbstractRefShape{2}}
+    points = Ferrite.getpoints(qr)
+    N = map(ξ -> Ferrite.reference_shape_value(ip, ξ, i), points)
+    vertices = Ferrite.reference_coordinates(Lagrange{RefShape, 1}())
+    push!(vertices, first(vertices))
+
+    fig = Plt.Figure()
+    ax = Plt.Axis(fig[1,1])
+    Plt.lines!(ax, first.(vertices), last.(vertices))
+    Plt.arrows!(ax, first.(points), last.(points), first.(N), last.(N); lengthscale = 0.3)
+    return fig
+end
+
+function plot_global_shape_function(ip, qr, nx, ny, i)
+    #TODO: Plot a single global shape function to investigate continuity
+end