34. Mixed Methods for second order equations

Again, we consider the second order equation

\[ -\Delta u = f, \qquad u = 0 \text{ on } \partial \Omega \]

We introduce now a new variable \(\sigma\) for \(\nabla u\), and rewrite the second order equation as first order system:

\[\begin{split} \begin{array}{ccccl} \sigma & - & \nabla u & = & 0 \\ \operatorname{div} \sigma & & & = & -f \end{array} \end{split}\]

We multiply the first equation by a test-function \(\tau\), and the secon d by a test-function \(v\). We use integration by parts \(\int_\Omega \nabla u \tau = -\int_\Omega u \operatorname{div} \tau + \int_{\partial \Omega} u \tau_n\).

The weak formulation is: find \(\sigma \in \Sigma := H(\operatorname{div})\) and \(u \in V := L_2\) such that

\[\begin{split} \begin{array}{ccccll} \int \sigma \tau & + & \int u \operatorname{div} \tau & = & 0 & \forall \, \tau \in \Sigma \\ \int v \operatorname{div} \sigma &&& = & 0 & \forall \, v \in V \end{array} \end{split}\]

The function space is

\[ H(\operatorname{div}) = \{ \tau \in [L_2(\Omega)]^d : \operatorname{div} \tau \in L_2(\Omega) \} \]

We don’t need derivatives of \(u\) and \(v\), thus these fields are chosen in \(L_2\).

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

Sigma = HDiv(mesh, order=2)
V = L2(mesh, order=1)
X = Sigma*V

sigma,u = X.TrialFunction()
tau,v = X.TestFunction()

a = BilinearForm( (sigma*tau+div(sigma)*v+div(tau)*u)*dx).Assemble()
f = LinearForm(-1*v*dx).Assemble()

gfu = GridFunction(X)
gfu.vec.data = a.mat.Inverse() * f.vec

Draw (gfu.components[0], mesh, "flux")
Draw (gfu.components[1], mesh, "u");

The total outflow is in balance with the source:

\[ \int_{\partial \Omega} \sigma \cdot n = \int_\Omega \operatorname{div} \sigma \, 1 = -\int_\Omega f \, 1 \]

n = specialcf.normal(mesh.dim)
Integrate (gfu.components[0]*n, mesh, BND)

-0.9999999999999974

34.1. Boundary conditions:

From integration by parts, Dirichlet boundary conditions \(u = u_D\) on \(\Gamma_D\) enter now in the linear form:

\[ f(v) = \int_{\Omega} f v \, dx + \int_{\Gamma_D} u_D \tau_n \, ds \]

However, Neumann boundary conditions are now built-in into the space for \(\sigma\):

\[ \sigma_n = g \qquad \text{on } \Gamma_N \]

The role of essential and natural boundary conditions is now the opposite as for a standard, primal formulation.

In NGSolve we use Sigma = HDiv(mesh, order=2, dirichlet="GammaN") for setting essential boundary conditions for \(\sigma\).