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.

Solve Surface PDE on stl geometry with boundary

More
2 years 5 months 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 :)
More
2 years 5 months ago - 2 years 5 months 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: 2 years 5 months ago by Kai.
More
2 years 5 months 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:
Time to create page: 0.119 seconds

Login