How to do Nonzero Neumann boundary condition for Laplace and heat equations

[aquillen]

I am having trouble understanding how to implement a non-zero Neumann boundary condition. For both Laplace's equation and the heat equation. The examples seem to not quite be simple enough for me to understand.

#Starting with this: (for Laplace's equation)

from skfem.helpers import dot, grad
m = skfem.MeshTri().refined(5)
e = skfem.ElementTriP1()
basis = skfem.Basis(m, e)
fbasis = skfem.FacetBasis(m,e)

@skfem.BilinearForm
def laplace(u, v, _):
    return dot(grad(u), grad(v))

@skfem.LinearForm
def rhs(v,_):  
    return 1.0 * v    # I think this will give me a Neumann boundary of normal derivative = 1 

Dlr =basis.get_dofs({'left', 'right'})    # trying a Neumann BC on left right
Dtb = basis.get_dofs({'top', 'bottom'})   # trying a Dirichlet BC on top bottom

rbasis = skfem.FacetBasis(m, e, facets=m.boundaries['right'])   # gives me facet basis for
# right boundary only?  Can I instead pass the dofs I made above somehow to this? No probably not. 
lbasis = skfem.FacetBasis(m, e, facets=m.boundaries['left'])   # gives me facet basis for
# left boundary only
B = skfem.asm(rhs,rbasis) + skfem.asm(rhs,lbasis)    # both left and right boundaries giving du/dn  = 1 
A = skfem.asm(laplace,basis)  # 

u = basis.zeros()
u[Dtb] = 1.0  # dirichlet BC on top and bottom?

x = skfem.solve(*skfem.condense(A, B, x=u, D=Dtb)) # Dirichlet boundary of 1 on top bottom ?

skfem.visuals.matplotlib.plot(basis, x)  # show solution

# I got something with this, but is it correct?

Thanks for any insight, suggestions or pointing me to more examples to study. I know how to write the boundary term in the weak form of the PDE, but I am not sure how to evaluate the boundary term in such a way to include it in a solve.

For the heat equation with similar boundaries:

from scipy.sparse.linalg import splu  # for sparse matrices

m = skfem.MeshTri().refined(5)
e = skfem.ElementTriP1()
basis = skfem.Basis(m, e)
#fbasis = skfem.FacetBasis(m,e)

#Dlr = basis.get_dofs({'left', 'right'})  # we don't use this for the Neumann boundary conditions (but why not?)
Dtb = basis.get_dofs({'top', 'bottom'})  # needed for Dirichlet boundary conditions 

@skfem.BilinearForm
def laplace(u, v, _):
    return skfem.helpers.dot(grad(u), grad(v))

@skfem.LinearForm
def rhs(v,_):
    return 1.0 * v

#Dlr = basis.get_dofs({'left', 'right'})  # we don't use this for the Neumann boundary conditions (but why not?)
Dtb = basis.get_dofs({'top', 'bottom'})  # needed for Dirichlet boundary conditions 

Diffusivity = 1.0
dt = 0.01
u_nb = 1.0  # value of normal deriv on right/left boundary 

L = Diffusivity * skfem.asm(laplace, basis)   # Laplacian 
M = skfem.asm(mass, basis)  # mass matrix (sizes of triangles) 

L0, M0 = skfem.penalize(L, M, D=Dtb)  #  Dirichlet BCs on top+bottom of the boundary 

theta = 0.5       # Crank–Nicolson algorithm for a diffusion PDE, rest has Neumman BC of zero so far 
A = M0 + theta * L0 * dt
B = M0 - (1 - theta) * L0 * dt  
backsolve = splu(A.T).solve

rbasis = skfem.FacetBasis(m, e, facets=m.boundaries['right'])   # gives me facet basis for
# right boundary only
lbasis = skfem.FacetBasis(m, e, facets=m.boundaries['left'])   # gives me facet basis for
# left boundary only

Bint = (skfem.asm(rhs,rbasis) + skfem.asm(rhs,lbasis) )*Diffusivity * u_nb *dt

def one_step(t,u):
    t += dt
    u = backsolve(B @ u + Bint)    # is this okay?
    # u[Dlr] += dt* u_nb * Diffusivity # definitely wrong 
    return t,u

#

it's this last evolve step that I think can't be right. Any suggestions on how to do this would be appreciated!

[kinnala]

I'm not exactly sure what you are looking for, but I will give you a simple example.

I will solve $-\Delta u = f$ in a square $[0, 1] \times [0, 1]$ with inhomogeneous $\partial u/\partial n = g$ on the right hand side boundary and $u=0$ on the remaining boundaries.

[kinnala]

I will denote the right hand side boundary as $\Gamma$. The weak formulation of this problem is: find $u \in V = \{v \in H^1([0,1]^2) : v|_{\partial \Omega \setminus \Gamma}=0 \}$ such that

$$(\nabla u, \nabla v) = (f, v) + (g, v)_\Gamma$$

for every $v \in V$.

[kinnala]

This will clearly become a matrix system of the form $Ax = b$ where $A$ consists of the Laplacian and $b$ consists of the two terms $(f, v)$ and $(g,v)_\Gamma$.

[kinnala]

Here is the code (scikit-fem==10.0.1):

from skfem import *
from skfem.helpers import *

m = MeshTri().refined(5).with_defaults()  # left, right, etc. tagged
e = ElementTriP1()
basis = Basis(m, e)

fbasis_gamma = basis.boundary('right')

@BilinearForm
def laplace(u, v, _):
    return dot(grad(u), grad(v))

@LinearForm
def bodyload(v, w):
    x, y = w.x
    f = 1 + 0 * x  # this can depend on x and y
    return f * v

@LinearForm
def neumann(v, w):
    x, y = w.x
    g = y ** 2  # this can depend on x and y
    return g * v


A = asm(laplace, basis)
b = asm(bodyload, basis) + asm(neumann, fbasis_gamma)

u = basis.zeros()
x = solve(*condense(A, b, D=basis.get_dofs({'left', 'top', 'bottom'})))

basis.plot(x, shading='gouraud', colorbar=True).show()

[kinnala]

In your code you had a comment:

Can I instead pass the dofs I made above somehow to this? No probably not.

DOFs cannot be passed when creating Basis objects. In fact, DOFs have not been defined before Basis object is created because the numbering of the DOFs happens while creating Basis objects.

[aquillen]

Thank you! Happily I seem to have been pretty close!
When I first started writing the question, I was very frustrated -- and giving up --- but then the format for using the facet basis started to make sense! I can follow your example! So nice!