Help with Example Electrostatics Model

mneunteufel
Offline
Premium Member

4 years 7 months ago #2555 by mneunteufel

Replied by mneunteufel on topic Help with Example Electrostatics Model

Hi jleman,

if I load the saved mesh in python via

Code:

mesh = Mesh('dome.vol')

all boundaries seem to be saved and loaded properly.

With

Code:

for v in mesh.vertices:
            print(v.point)

you get all mesh points.

If you solve the problem in NGSolve and assume your solution is saved in the GridFunction gfu then you can access data at a point p=(p_x, p_y, p_z) by

Code:

value = gfu(mesh(p_x,p_y,p_z))

If you have a CoefficientFunction Q representing the charge density you can get values on the surface by

Code:

value = Q(mesh(p_x,p_y,p_z, BND))

Of course the point must be a surface point on the mesh.

If you have the equations by hand I would recommend doing everything in NGSolve via Python.

Best
Michael

The following user(s) said Thank You: jleman

Yash.sooriya
New Member

2 years 1 month ago #4501 by Yash.sooriya

Replied by Yash.sooriya on topic Help with Example Electrostatics Model

Hello mneunteufel,

I was wondering if you still have access to this dome.py file or you remember how to do this, as it would be very helpful for my current project. Thanks!

christopher
Offline
Administrator

2 years 1 month ago #4504 by christopher

Replied by christopher on topic Help with Example Electrostatics Model

Hi,
This code snippet should solve this problem:

Code:

from netgen.occ import *
from ngsolve import *
air = Sphere((0,0,0), 1)
air.faces.name = "air"

ground = HalfSpace((0,0,0), (0,0,-1))
ground.faces.name = "ground"

box = Box((-0.1,-0.1,-0.1), (0.1,0.1,0.1))
box.faces.name = "volt"
box.faces.maxh = 0.02

domain = (air * ground) - box

geo = OCCGeometry(domain)

mesh = Mesh(geo.GenerateMesh(maxh=0.2))
mesh.Curve(3)
Draw(mesh)

fes = H1(mesh, order=3, dirichlet="ground|volt")
u,v = fes.TnT()

eps = 1

a = BilinearForm(fes)
a += eps * grad(u) * grad(v) * dx

u = GridFunction(fes)

with TaskManager():
    r = u.vec.CreateVector()
    a.Assemble()
    u.Set(1, definedon=mesh.Boundaries("volt"))
    r.data = - a.mat * u.vec
    u.vec.data += a.mat.Inverse(fes.FreeDofs()) * r
    Draw(u, mesh, "potential")
    Draw(-grad(u), mesh, "E")

Best
Christopher

Time to create page: 0.113 seconds