- Thank you received: 0
How to setup the Neumann Boundary Conditions?
- liubenyuan
- Topic Author
- Offline
- New Member
Less
More
4 years 11 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
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?
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
Less
More
- Thank you received: 0
4 years 11 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
Less
More
- Thank you received: 101
4 years 11 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
Then you need to add neumann boundary term to rhs
and for the inhomogeneous dirichlet you need to homogenize the problem as shown here:
ngsolve.org/docu/latest/i-tutorials/unit...chlet/dirichlet.html
Best
Christopher
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
Less
More
- Thank you received: 0
4 years 11 months ago - 4 years 11 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)
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 11 months ago by liubenyuan. Reason: update ground truth in code
- christopher
- Offline
- Administrator
Less
More
- Thank you received: 101
4 years 11 months ago - 4 years 11 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:
(inhomogeneous solution = homogeneous solution + inhomogeneous bc)
Best
Christopher
Code:
# solve
gfu.vec.data += a.mat.Inverse(fes.FreeDofs()) * r
(inhomogeneous solution = homogeneous solution + inhomogeneous bc)
Best
Christopher
Last edit: 4 years 11 months ago by christopher.
The following user(s) said Thank You: liubenyuan
Time to create page: 0.123 seconds