Forum Message

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

How to setup the Neumann Boundary Conditions?

liubenyuan
Topic Author
Offline
New Member

4 years 4 months ago #2213 by liubenyuan

How to setup the Neumann Boundary Conditions? was created by liubenyuan

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

www.math.uzh.ch/compmath/fileadmin/user/...ld_Noam/myManual.pdf

Section 4 3D Elliptic PDEs in NETGEN

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

Code:

## 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:
ngsolve.org/forum/ngspy-forum/672-neumann-boundary#1258
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?

liubenyuan
Topic Author
Offline
New Member

4 years 4 months ago #2214 by liubenyuan

Replied by liubenyuan on topic How to setup the Neumann Boundary Conditions?

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.

Code:

# 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()

christopher
Offline
Administrator

4 years 4 months ago #2217 by christopher

Replied by christopher on topic How to setup the Neumann Boundary Conditions?

You need to change the fes line to not have dirichlet on the right side

Code:

fes = ngs.H1(mesh, order=2, dirichlet='left')

Then you need to add neumann boundary term to rhs

Code:

# boundary term for neumann 
f += neuval*v*ds(definedon="right")

and for the inhomogeneous dirichlet you need to homogenize the problem as shown here:
ngsolve.org/docu/latest/i-tutorials/unit...chlet/dirichlet.html

Code:

# homogenization of inhomogeneous dirichlet problem
u.Set(dirval, definedon=mesh.Boundaries("left"))
res = f.vec.CreateVector()
res.data = f.vec - a.mat * u.vec
u.vec.data += a.mat.Inverse(fes.FreeDofs()) * res

Best
Christopher

The following user(s) said Thank You: liubenyuan

liubenyuan
Topic Author
Offline
New Member

4 years 4 months ago - 4 years 4 months ago #2218 by liubenyuan

Replied by liubenyuan on topic How to setup the Neumann Boundary Conditions?

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)

Code:

# 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')

Last edit: 4 years 4 months ago by liubenyuan. Reason: update ground truth in code

christopher
Offline
Administrator

4 years 4 months ago - 4 years 4 months ago #2219 by christopher

Replied by christopher on topic How to setup the Neumann Boundary Conditions?

You have a += for the inhomogeneous dirichlet problem:

Code:

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

(inhomogeneous solution = homogeneous solution + inhomogeneous bc)

Best
Christopher

Last edit: 4 years 4 months ago by christopher.

The following user(s) said Thank You: liubenyuan

Time to create page: 0.140 seconds

Login

NGSolve on Twitter