Neumann condition essential for Mixed Poisson equation

[mhtber]

Hi,
How to impose a non-homogenous Neumann "essential" condition $\sigma \cdot n = g$ for the mixed Poisson equation.
Thanks

[kinnala]

What do we know about $g$? In particular, is it some explicit function in x and y or some cellwise constant data array?

[mhtber]

g is an explicit function in x and y

[kinnala]

Let's say you use the lowest order Raviart-Thomas element. The DOFs are
$$\int_E \sigma \cdot n$$
over the element edge $E$.

I think that, given $g$, you can find the value
$$\int_E g$$
for each edge $E$ on the boundary and set $\sigma$ to the value using skfem.utils.condense.

To calculate the integral, you can use FacetBasis. Something like

fbasis = FacetBasis(mesh, ElementTriRT0())

@Functional
def flux(w):
    return g(w.x[0], w.x[1])

gE = flux.elemental(fbasis)

Then pass gE to condense. I don't have my laptop with me so cannot test right now.

[mhtber]

Sorry, I'm new to scikit-fem. For instance to solve the Poisson equation with $u_{ex} = 1 + x^2 + 2y^2$ as the exact solution and Neumann B.C. only on the right. Is this the way to proceed?

import numpy as np
from skfem import *
from skfem.helpers import dot, div

mesh = MeshTri().refined(6)
element = {'sigma': ElementTriRT0(), 'u': ElementDG(ElementTriP0())}

basis = {variable: Basis(mesh, e, intorder=4) for variable, e in element.items()}

@BilinearForm
def a(sigma, tau, _):
    return dot(sigma, tau)
A = asm(a, basis['sigma'])

def b(sigma, v, _):
    return div(sigma) * v
B = asm(b, basis['sigma'], basis['u'])

def f(v, _):
    return 6 * v
Fu = asm(f, basis['u'])

@LinearForm
def g_D(tau, w):
    x = w.x
    return (dot(tau, w.n) * (1 + x[0]*x[0] + 2*x[1]*x[1]) )
Fsigma = asm(g_D, basis['sigma'].boundary({"left","bottom","top"}))
   
from scipy.sparse import bmat
K = bmat([[A, B.T], [B, None]]).tocsr()
F = np.concatenate([Fsigma, Fu])

fbasis = FacetBasis(mesh, element["sigma"], facets=mesh.boundaries["right"])
@Functional
def flux(w):
    return 2*w.x[0]

gE = flux.elemental(fbasis)
D = fbasis.get_dofs("right")
xD = np.concatenate([basis["sigma"].zeros(), basis["u"].zeros()])
xD[D] = gE

x = solve(*condense(K, F, x=xD, D=D))

[kinnala]

I made some small changes, hopefully it is correct. You can compare against your analytical solution:

import numpy as np
from skfem import *
from skfem.helpers import dot, div

mesh = MeshTri().refined(6)
element = {'sigma': ElementTriRT0(), 'u': ElementTriP0()}

basis = {variable: Basis(mesh, e, intorder=4) for variable, e in element.items()}

@BilinearForm
def a(sigma, tau, _):
    return dot(sigma, tau)
A = asm(a, basis['sigma'])

@BilinearForm
def b(sigma, v, _):
    return div(sigma) * v
B = asm(b, basis['sigma'], basis['u'])

@LinearForm
def f(v, _):
    return 6 * v
Fu = asm(f, basis['u'])

@LinearForm
def g_D(tau, w):
    x = w.x
    return (dot(tau, w.n) * (1 + x[0]*x[0] + 2*x[1]*x[1]) )
Fsigma = asm(g_D, basis['sigma'].boundary({"left","bottom","top"}))
   
from scipy.sparse import bmat
K = bmat([[A, B.T], [B, None]]).tocsr()
F = np.concatenate([Fsigma, Fu])

fbasis = FacetBasis(mesh, element["sigma"], facets=mesh.boundaries["right"])
@Functional
def flux(w):
    return 2*w.x[0]

gE = flux.elemental(fbasis)
D = fbasis.get_dofs("right")
xD = np.concatenate([basis["sigma"].zeros(), basis["u"].zeros()])
xD[D] = gE

x = solve(*condense(K, F, x=xD, D=D))

basis["u"].plot(x[basis["sigma"].N:], colorbar=True).show()

[mhtber]

Thanks a lot, it works perfectly.

import numpy as np
from skfem import *
from skfem.helpers import dot, div

mesh = MeshTri().refined(6)
element = {'sigma': ElementTriRT0(), 'u': ElementTriP0()}

basis = {variable: Basis(mesh, e, intorder=4) for variable, e in element.items()}

@BilinearForm
def a(sigma, tau, _):
    return dot(sigma, tau)
A = asm(a, basis['sigma'])

@BilinearForm
def b(sigma, v, _):
    return div(sigma) * v
B = asm(b, basis['sigma'], basis['u'])

@LinearForm
def f(v, _):
    return 6 * v
Fu = asm(f, basis['u'])

@LinearForm
def g_D(tau, w):
    x = w.x
    return (dot(tau, w.n) * (1 + x[0]*x[0] + 2*x[1]*x[1]) )
Fsigma = asm(g_D, basis['sigma'].boundary({"left","bottom","top"}))
   
from scipy.sparse import bmat
K = bmat([[A, B.T], [B, None]]).tocsr()
F = np.concatenate([Fsigma, Fu])

fbasis = FacetBasis(mesh, element["sigma"], facets=mesh.boundaries["right"])
@Functional
def flux(w):
    return 2*w.x[0]

gE = flux.elemental(fbasis)
D = fbasis.get_dofs("right")
xD = np.concatenate([basis["sigma"].zeros(), basis["u"].zeros()])
xD[D] = gE

x = solve(*condense(K, F, x=xD, D=D))
sigma, u = np.split(x, [A.shape[0]])

basis['u'].plot(u, shading='gouraud', colorbar=True).show()

@Functional
def error_u(w):
    x = w.x
    u = w['u']
    return (u -(1 + x[0]*x[0] + 2*x[1]*x[1])) ** 2
err_u_L2 = np.sqrt(error_u.assemble(basis['u'], u=u))
print(err_u_L2)

@Functional
def error_sigma(w):
    x = w.x
    sigma = w['sigma']
    return (sigma[0] - 2*x[0]) ** 2 + (sigma[1] - 4*x[1]) ** 2
err_sigma_L2 = np.sqrt(error_sigma.assemble(basis['sigma'], sigma=sigma))
print(err_sigma_L2)

@Functional
def div_error(w):
    x = w.x
    sigma = w['sigma']    
    return (div(sigma) - 6)**2
err_div = div_error.assemble(basis['sigma'], sigma=sigma)
print(err_div)