Ferrite with Robin Bcs?

[xtalax]

Noticed someone having trouble here

https://www.reddit.com/r/Julia/comments/14wrcyq/ferrite_fem_package_with_robin_boundary_conditions/

[fredrikekre]

Thanks for the repost, I don't keep tabs on reddit at all.

I have a student who did this recently and suggested that he add a section in the docs about this, I will check what the status is, but can also post some notes about it here later.

[fredrikekre]

Consider the following strong form

$$ \begin{aligned} \boldsymbol{\nabla} \cdot \boldsymbol{q} &= f \quad &\forall \boldsymbol{x} \in \Omega \ u &= g_\mathrm{D} \quad &\forall \boldsymbol{x} \in \Gamma_\mathrm{D} \ q_\mathrm{n} := \boldsymbol{q} \cdot \boldsymbol{n} &= g_\mathrm{N} \quad &\forall \boldsymbol{x} \in \Gamma_\mathrm{N}\ a\ u + b\ q_\mathrm{n} = a\ u + b\ \boldsymbol{q} \cdot \boldsymbol{n} &= g_\mathrm{R} \quad &\forall \boldsymbol{x} \in \Gamma_\mathrm{R} \end{aligned} $$

where $u = u(\boldsymbol{x})$ is the temperature, $\boldsymbol{q} = -k_\Omega \boldsymbol{\nabla} u$ the flux (Fourier's law), $\Omega$ the domain, $\Gamma = \Gamma_\mathrm{D} \cup \Gamma_\mathrm{N} \cup \Gamma_\mathrm{R}$ the boundary of $\Omega$ (split into Dirichlet, Neumann, and Robin parts), and $\boldsymbol{n}$ the outwards normal vector. On the Dirichlet part the temperature is prescribed ($u = g_\mathrm{D}$), on the Neumann part the (normal) flux is prescribed $q_\mathrm{n} = g_\mathrm{N}$), and on the Robin part a linear combination of temperature and (normal) flux is prescribed ($a\ u + b\ q_\mathrm{n} = g_\mathrm{R}$).

The corresponding weak form, with test function $v$, can then be derived to: Find $u \in \mathbb{U}$ s.t.

$$ \int_\Omega k_\Omega\ \boldsymbol{\nabla} v \cdot \boldsymbol{\nabla} u\ \mathrm{d}\Omega - \int_{\Gamma_\mathrm{R}}v \frac{a}{b} u\ \mathrm{d}\Gamma = \int_\Omega v\ f\ \mathrm{d}\Omega - \int_{\Gamma_\mathrm{N}} v\ g_\mathrm{N}\ \mathrm{d}\Gamma - \int_{\Gamma_\mathrm{R}} v\ \frac{g_{\mathrm{R}}}{b}\ \mathrm{d}\Gamma \quad \forall v \in \mathbb{U}^0. $$

Of interest for the Robin boundary is the boundary integrals over $\Gamma_\mathrm{R}$ on the left and right hand sides, which will give contributions to the tangent matrix, and the right hand side, respectively.

With FE discretizations for $u$ and $v$:

$$ u \approx u_\mathrm{h} = \sum_{i = 1}^{N} \phi_i a_i, \quad u \approx v_\mathrm{h} = \sum_{i = 1}^{N} \phi_i b_i $$

we obtain the discrete system $K\ a = f$, where

$$ \begin{aligned} K_{ij} = \int_\Omega k_\Omega\ \boldsymbol{\nabla} \phi_i \cdot \boldsymbol{\nabla} \phi_j\ \mathrm{d}\Omega - \int_{\Gamma_\mathrm{R}}\phi_i \frac{a}{b} \phi_j\ \mathrm{d}\Gamma\ f_{i} = \int_\Omega \phi_i\ f\ \mathrm{d}\Omega - \int_{\Gamma_\mathrm{N}} \phi_i\ g_\mathrm{N}\ \mathrm{d}\Gamma - \int_{\Gamma_\mathrm{R}} \phi_i\ \frac{g_{\mathrm{R}}}{b}\ \mathrm{d}\Gamma \end{aligned} $$

In Ferrite this would be implemented something like (showing only the evaluation of the Robin integrals here):

using Ferrite

# FE basis and quadrature (2D example with linear quads)
ip = Lagrange{2, RefCube, 1}()
qr = QuadratureRule{2, RefCube}(1)
fv = FaceScalarValues(qr, ip)

# Global matrix/vector
K = create_sparsity_pattern(dh)
f = zeros(ndofs(dh))

# Local matrix/vector
Ke = zeros(ndofs_per_cell(dh), ndofs_per_cell(dh))
fe = zeros(ndofs_per_cell(dh))

# Loop for the Robin boundary
for (cellid, faceid) in getfaceset(grid, "RobinBC")
    # Reset the buffers
    fill!(Ke, 0)
    fill!(fe, 0)
    # Update FE basis
    reinit!(fv, getcoordinates(grid, cellid), faceid)
    # Compute the local contribution
    for qp in 1:getnquadpoints(fv)
        dΓ = getdetJdV(fv, qp)
        for i in 1:getnbasefunctions(fv)
            ϕᵢ = shape_value(fv, qp, i)
            fe[i] -= ( ϕᵢ * gR / b ) * dΓ
            for j in 1:getnbasefunction(fv)
                ϕⱼ = shape_value(fv, qp, j)
                Ke[i, j] -= ( ϕᵢ * a / b * ϕⱼ ) * dΓ
            end
        end
    end
    # Assemble local contribution
    dofs = celldofs(dh, cellid)
    for (i, I) in pairs(dofs)
        f[I] += fe[i]
        for (j, J) in pairs(dofs)
            K[I, J] += Ke[i, j]
        end
    end
end

Another common way to define the Robin boundary condition is

$$ q_\mathrm{n} = -k_\Gamma (u_\Gamma - u) $$

where $k_\Gamma$ is the interface conductivity, and $u_\Gamma$ the ambient temperature. With this parametrization we can identify e.g.

$$ \begin{aligned} b &= 1\ a &= - k_\Gamma\ g_\mathrm{R} &= - k_\Gamma u_\Gamma \end{aligned} $$

which would result in the (perhaps more recognizable?) system

$$ \int_\Omega k_\Omega\ \boldsymbol{\nabla} v \cdot \boldsymbol{\nabla} u\ \mathrm{d}\Omega + \int_{\Gamma_\mathrm{R}}v\ k_\Gamma\ u\ \mathrm{d}\Gamma = \int_\Omega v\ f\ \mathrm{d}\Omega - \int_{\Gamma_\mathrm{N}} v\ g_\mathrm{N}\ \mathrm{d}\Gamma + \int_{\Gamma_\mathrm{R}} v\ k_\Gamma\ u_\Gamma\ \mathrm{d}\Gamma. $$

(Someone please check the signs :))