Ctrl+K

Application of the abstract theory

46. Application of the abstract theory#

\(\DeclareMathOperator{\opdiv}{div}\)

We consider now second order equations with mixed boundary conditions, and variable coefficients:

Find \(u\) such that

\[\begin{eqnarray*} -\opdiv \lambda(x) \nabla u & = & f \quad \; \; \; \text{ in } \Omega \\ u & = & u_D \quad \text{ on } \Gamma_D \\ \lambda \frac{\partial u}{\partial n} & = & g \quad \; \; \; \text{ on } \Gamma_N \end{eqnarray*}\]

Introduce \(\sigma = \lambda \nabla u\) to obtain the system

\[ \lambda^{-1} \sigma - \nabla u = 0 \quad \quad \opdiv \sigma = -f \]

and boundary conditions

\[ u = u_D \; \; \text{ on } \Gamma_D \qquad \sigma \cdot n = g \; \; \text{ on } \Gamma_N \]

We test the first equation by \(\tau\), the second by \(v\) and get

\[\begin{split} \begin{array}{ccccll} \int \lambda^{-1} \sigma \tau & - & \int \nabla u \tau & = & 0 & \forall \, \tau \\ \int \opdiv \sigma \, v & & & = & - \int f v & \forall \, v \end{array} \end{split}\]

To get a mixed system with the same bilinear-form at the lower left and upper right blocks, we have to integrate by parts either one of the terms.

46.1. Dual mixed formulation#

Integration by parts of the upper right term we get

\[ - \int_\Omega \nabla u \tau = \int_\Omega u \opdiv \tau - \int_{\Gamma_D} u \, \tau_n - \int_{\Gamma_N} u \, \tau_n \]

Now, we can plug in the boundary condition \(u = u_D\) for the Dirichlet boundary, and move the term to the right hand side. On the Neumann boundary, we have to restrict the test-function to \(\tau_n = 0\). Neumann boundary conditions become now essential ones.

The problem is now formulated as:

Find \(\sigma \in H(\opdiv)\) and \(u \in L_2\) such that \(\sigma_n = g\) on \(\Gamma_N\) and

\[\begin{split} \begin{array}{ccccll} \int \lambda^{-1} \sigma \tau & + & \int \opdiv \tau \, u & = & \int_{\Gamma_D} u_D \tau_n & \forall \, \tau, \; \tau_n = 0 \text{ on } \Gamma_N \\ \int \opdiv \sigma \, v & & & = & - \int f v & \forall \, v \end{array} \end{split}\]

Note that the role of essential and natural boundary conditions is now swapped.

There are some issues left to verify:

  • Is it allowed to take boundary values of \(H(\opdiv)\) functions ?

  • What kind of essential Neumann boundary data can we prescribe ?

  • Is \(\int u_D \tau_n\) a continuous linear form ?

from ngsolve import *
from ngsolve.webgui import Draw

from netgen.occ import *
outer = Rectangle(1,1).Face()
outer.edges.Max(X).name='r'
outer.edges.Min(X).name='l'
outer.edges.Max(Y).name='t'
outer.edges.Min(Y).name='b'
outer.faces.name="air"

rect = MoveTo(0.2,0.2).Rectangle(0.7,0.1).Face()
rect.faces.name="bar"
outer -= rect

circ = Circle( (0.3,0.7), 0.1).Face()
circ.faces.name="source"
outer -= circ
shape = Glue([outer,rect,circ])

mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.03)).Curve(3)

fes = H1(mesh, order=3, dirichlet="b|r")
u = fes.TrialFunction()
v = fes.TestFunction()

lam = mesh.MaterialCF( { "bar" : 100 }, default=1 )
a = BilinearForm(fes)
a += lam * grad(u) * grad(v)*dx
f = LinearForm(fes)
f += 1 * v *dx("source")
a.Assemble()
f.Assemble()
gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Draw (gfu, mesh, "stdtemp")
Draw (-lam*grad(gfu), mesh, "stdflux", min=0, max=0.1);

Sigma = HDiv(mesh, order=2, dirichlet="t|l")
V = L2(mesh, order=1)
X = Sigma*V

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

a = BilinearForm(X)
a += (1/lam*sigma*tau + div(sigma)*v + div(tau)*u) * dx
a.Assemble()

f = LinearForm(X)
f += 1 * v *dx("source")
f.Assemble()

gfu = GridFunction(X)
gfu.vec.data = a.mat.Inverse(X.FreeDofs()) * f.vec
Draw (-gfu.components[1], mesh, "mixedtemp")
Draw (gfu.components[0], mesh, "mixedflux", min=0, max=0.1);

We observe that the normal component of the flux is continuous arcross discontinuous coefficients. As we will see, this is a property of the function space \(H(\operatorname{div})\).

46.2. Primal mixed formulation#

We can also integrate the lower left term by part. This leads to the formulation

Find \(u \in H^1\) and \(\sigma \in [L_2]^d\) such that \(u = u_D\) on \(\Gamma_D\) and

\[\begin{split} \begin{array}{ccccll} \int \lambda^{-1} \sigma \tau & - & \int \tau \nabla u & = & 0 & \forall \, \tau \\ -\int \sigma \nabla v & & & = & - \int f v - \int_{\Gamma_N} g v & \forall \, v, \; \; v = 0 \text{ on } \Gamma_D \end{array} \end{split}\]

This formulation is very close to the primal method.

Contents