Solve Surface PDE on stl geometry with boundary

Main Forum Old-Forum - do not post here
1

Hi all :slight_smile:

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 :slight_smile:

2

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 …

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")
3

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”)[/code]
https://ngsolve.org/media/kunena/attachments/889/surface_pde.py

Best,
Michael

Attachment: 70741d8fdcc071269c78e9fa9db83d65