How to setup the Neumann Boundary Conditions?

Main Forum Old-Forum - do not post here
1

Hi, I want to set up the Neumann boundary conditions. Firstly, I want to recreate the tutorial in

Section 4 3D Elliptic PDEs in NETGEN

But I do not know how to create the Neumann BCs in Python from these codes

## coefficient for source integral of linear form
define coefficient f
(−4*(x+1)*(z*z+y*y−2)) ,
## coefficient g/beta for Neumann integral of linear form
define coefficient g_beta
0 , (2*(y*y−1)*(z*z−1)) ,
## f i n i t e element space with order=3 and Dirichlet−boundary at −bc=1
define fespace v −order=3 −dirichlet =[1]
## grid function to store solution
define gridfunction u −fespace=v
## linear form
define linearform f −fespace=v
source f
neumann g_beta

Moreover, for the Neumann BCs, I need to set the Zero potential point in the mesh. I found a post about it here:

How could I do this? should I force a Dirichlet BCs on a Node ID?

Is there any tutorial in ngspy on the Neumann BCs?

2

For example, In the following 1D poisson example code, how could I set both the dirichlet and Neumann boundary conditions? For example, f(left) = 10, and f’(right) = -1.

# https://ngsolve.org/forum/ngspy-forum/167-1-d-poisson-with-constant-forcing-function
import matplotlib.pyplot as plt

import ngsolve as ngs
from ngsolve import grad, dx
import netgen.meshing as ng

"""1. generate a 1D mesh"""
m = ng.Mesh()
m.dim = 1
nel = 10
pnums = []

# add nodes and ID
for i in range(nel+1):
    node = ng.Pnt(i/nel, 0, 0)
    m.Add(ng.MeshPoint(node))
    pnums.append(i+1)  # node index starts 1

# add elements
for i in range(nel):
    m.Add(ng.Element1D([pnums[i], pnums[i+1]], index=1))

# set materials
m.SetMaterial(1, 'material')

# add BC points
m.Add(ng.Element0D(pnums[0], index=1))
m.Add(ng.Element0D(pnums[nel], index=2))

# set boundary condition names
m.SetBCName(0, 'left')
m.SetBCName(1, 'right')

print(pnums)

"""2. compute on this mesh"""
mesh = ngs.Mesh(m)
print(mesh.GetBoundaries())

# Specify Dirichlet boundary conditions
fes = ngs.H1(mesh, order=2, dirichlet='left|right')
# fes = ngs.H1(mesh, order=1, dirichlet=[1, 2])
print("ndof: ", fes.ndof)
print("freedofs: ", fes.FreeDofs())

u = fes.TrialFunction()  # symbolic object
v = fes.TestFunction()   # symbolic object

# coefficient/stiffness matrix
a = ngs.BilinearForm(fes)
a += grad(u) * grad(v) * dx
a.Assemble()
print("mat = \n", a.mat)

# Forcing function
f = ngs.LinearForm(fes)    
f += 1*v*dx
f.Assemble()
print("rhs = \n", f.vec)

# boundary conditions
# how to set:
# 1. f(left) = a, dirichlet boundary
# 2. f'(right)*n = b, Neumann boundary

# solution
u = ngs.GridFunction(fes)
u.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
print ("sol = \n", u.vec)

"""3. plot"""
pnts = [i/100 for i in range(101)]
pnts_vals = [(x, u(x)) for x in pnts if mesh.Contains(x)]

fig, ax = plt.subplots()

pnts, vals = zip(*pnts_vals)
ax.plot(pnts, vals, "-b")
# ax.plot(u.vec, "-b")
plt.show()
3

Dear Christopher, Thank you!

I rewrite the code, and now it can handle both the dirichelt and the neumann BCs! But when I plot the results, it seems like that gfu interpolate the value between [0, 0.1], Is it the right behavior of ngsolve? The problem is now -f’‘=2, f’(1)=-1, f(0)=1, and the ground truth is -x^2 + x + 1.

For example, in the following code, gfu(0.05) is 0.5474999999999999, which should be a value larger than f(0)

# https://ngsolve.org/forum/ngspy-forum/167-1-d-poisson-with-constant-forcing-function
import matplotlib.pyplot as plt

import ngsolve as ngs
from ngsolve import grad, dx, ds, BND
import netgen.meshing as ng

"""1. generate a 1D mesh"""
m = ng.Mesh()
m.dim = 1
nel = 10
pnums = []

# add nodes and ID
for i in range(nel+1):
    node = ng.Pnt(i/nel, 0, 0)
    m.Add(ng.MeshPoint(node))
    pnums.append(i+1)  # node index starts 1

# add elements
for i in range(nel):
    m.Add(ng.Element1D([pnums[i], pnums[i+1]], index=1))

# set materials
m.SetMaterial(1, 'material')

# add BC points
m.Add(ng.Element0D(pnums[0], index=1))
m.Add(ng.Element0D(pnums[nel], index=2))

# set boundary condition names
m.SetBCName(0, 'left')
m.SetBCName(1, 'right')

print(pnums)

"""2. compute on this mesh"""
mesh = ngs.Mesh(m)
print(mesh.GetBoundaries())

# Specify Dirichlet boundary conditions
fes = ngs.H1(mesh, order=2, dirichlet='left')
# fes = ngs.H1(mesh, order=1, dirichlet=[1, 2])
print("ndof: ", fes.ndof)
print("freedofs: ", fes.FreeDofs())

u = fes.TrialFunction()  # symbolic object
v = fes.TestFunction()   # symbolic object

# coefficient/stiffness matrix
a = ngs.BilinearForm(fes)
a += grad(u) * grad(v) * dx
a.Assemble()
print("mat = \n", a.mat)

# theoretic solution
# source: -f''= 2
# neumann: f'(1) = -1
# dirichlet: f(0) = 1
# f = -x^2 + x + 1

# Source
f = ngs.LinearForm(fes)
source_b = 2
f += source_b * v * dx

# boundary conditions
# 1. f'(right)*n = b, Neumann boundary
neumann_b = -1
f += neumann_b * v * ds(definedon='right')

# form
f.Assemble()
print("rhs = \n", f.vec)

# solution
gfu = ngs.GridFunction(fes)

# 1. f(left) = a, dirichlet boundary
dirichlet_b = 1
# gfu.Set(dirichlet_b, BND)
gfu.Set(dirichlet_b, definedon=mesh.Boundaries("left"))
r = f.vec.CreateVector()
r.data = f.vec - a.mat * gfu.vec

# solve
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * r
print ("sol = \n", gfu.vec)

"""3. plot"""
pnts = [i/100 for i in range(101)]
pnts_vals = [(x, gfu(x)) for x in pnts if mesh.Contains(x)]

fig, ax = plt.subplots()

pnts, vals = zip(*pnts_vals)
ax.plot(pnts, vals, "-b")
# ax.plot(gfu.vec, "-b")
plt.show()

# fig.savefig('ngsolve-2217.png')
4

You have a += for the inhomogeneous dirichlet problem:

# solve
gfu.vec.data += a.mat.Inverse(fes.FreeDofs()) * r

(inhomogeneous solution = homogeneous solution + inhomogeneous bc)

Best
Christopher