How to setup the Neumann Boundary Conditions?

More
4 years 11 months ago #2213 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?
More
4 years 11 months ago #2214 by liubenyuan
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()
More
4 years 11 months ago #2217 by christopher
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
More
4 years 11 months ago - 4 years 11 months ago #2218 by liubenyuan
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 11 months ago by liubenyuan. Reason: update ground truth in code
More
4 years 11 months ago - 4 years 11 months ago #2219 by christopher
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 11 months ago by christopher.
The following user(s) said Thank You: liubenyuan
Time to create page: 0.123 seconds