Solve Surface PDE on stl geometry with boundary

Kai
Topic Author
Offline
Junior Member

3 years 3 weeks ago #4060 by Kai

Solve Surface PDE on stl geometry with boundary was created by Kai

Hi all

i like to do the surface pde example but with a geometry that comes from some stl file.

is this possible ?

i would be happy if some one could tell me

Kai
Topic Author
Offline
Junior Member

3 years 2 weeks ago - 3 years 2 weeks ago #4061 by Kai

Replied by Kai on topic Solve Surface PDE on stl geometry with boundary

I have now combined the surface pde example with the, manual mesh generation example, i think thats a good startingpoint but i do not know how to set dirichlet bnd conditions, would be great if some one can help me

.

I get in the console "used dof inconsistency" ... if i run the following code ...

Code:

from ngsolve import *
from netgen.csg import *
from netgen.geom2d import unit_square
from ngsolve import Draw

ngmesh = netgen.meshing.Mesh(dim=2)

N=5
pnums = []
for i in range(N + 1):
    for j in range(N + 1):
        pnums.append(ngmesh.Add(netgen.meshing.MeshPoint(Pnt(i / N, j / N, 0))))


idx_dom = ngmesh.AddRegion("mat", dim=2)
for j in range(N):
    for i in range(N):
        ngmesh.Add(netgen.meshing.Element2D(idx_dom, [pnums[i + j * (N + 1)],
                               pnums[i + (j + 1) * (N + 1)],
                               pnums[i + 1 + (j + 1) * (N + 1)],
                               pnums[i + 1 + j * (N + 1)]]))

# horizontal boundaries
for i in range(N):
   ngmesh.Add(netgen.meshing.Element1D([pnums[N + i * (N + 1)],
                       pnums[N + (i + 1) * (N + 1)]], index=1))
   ngmesh.Add(netgen.meshing.Element1D([pnums[0 + i * (N + 1)], pnums[0 + (i + 1) * (N + 1)]], index=1))

# vertical boundaries
for i in range(N):
   ngmesh.Add(netgen.meshing.Element1D([pnums[i], pnums[i + 1]], index=2))
   ngmesh.Add(netgen.meshing.Element1D([pnums[i + N * (N + 1)], pnums[i + 1 + N * (N + 1)]], index=2))

mesh = Mesh(ngmesh)
Draw(mesh)

# geo          = CSGeometry()
# sphere       = Sphere(Pnt(0,0,0), 1)
# bot          = Plane(Pnt(0,0,0), Vec(0,0,-1))
# finitesphere = sphere * bot
#
# geo.AddSurface(sphere, finitesphere.bc("surface"))
# geo.NameEdge(sphere,bot, "bottom")
#
# mesh = Mesh(geo.GenerateMesh(maxh=0.3))
# mesh.Curve(2)
# Draw(mesh)
#
fes = H1(mesh, order=2, dirichlet=[1, 2])
u, v = fes.TnT()
print(fes.FreeDofs())

a = BilinearForm(fes, symmetric=True)
a += grad(u).Trace()*grad(v).Trace()*ds
a.Assemble()

force = 1.0

f = LinearForm(fes)
f += force*v*ds
f.Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs())*f.vec
Draw(gfu, mesh, "u")

Last edit: 3 years 2 weeks ago by Kai.

mneunteufel
Offline
Premium Member

3 years 2 weeks ago #4062 by mneunteufel

Replied by mneunteufel on topic Solve Surface PDE on stl geometry with boundary

Hi Kai,

with ngmesh = netgen.meshing.Mesh(dim=2) you are creating a 2D mesh, not a surface mesh in 3D.
Thus, with a += grad(u).Trace()*grad(v).Trace()*ds you are integrating only over the four edges. This gives the used_dof_inconsistence warning.

I changed/added the following lines, compare attached file, to get the desired result. As stated in the surface pde tutorial you need the keyword dirichlet_bbnd instead of dirichlet to set Dirichlet boundary conditions for surface PDEs.

Code:

ngmesh = netgen.meshing.Mesh(dim=3)

# to access edges easier
ngmesh.SetCD2Name(1, "top_bot")
ngmesh.SetCD2Name(2, "left_right")

fes = H1(mesh, order=2, dirichlet_bbnd="top_bot|left_right")

File Attachment:

File Name: surface_pde.py
File Size:2 KB

Best,
Michael

Attachments:

surface_pde.py

Time to create page: 0.107 seconds