Access GridFunction with Coordinates

NGSolve
1

Hello all,

I am relatively new to NGSolve/NETGEN and came across the following problem:

A short description of what I am doing:

I am trying to solve Poisson’s Equation for a 3D model of a simple spine model I created with NETGEN. This should give me, for a given electric dipole in the spine, potentials on the surface of the spine. This is part of the forward problem so I then later can solve the inverse problem (e.g. with sLORETA)

After doing some basic definitions for my problem:

fes = H1(mesh, order = '1')     # define Finite Element Space
u,v = fes.TnT()                 # define Trial and Test Space

sigma_coeff_ = {"spine" : 0.22, "vertebre" : 0.0014, "disc" : 0.008, "neck" : 0.01 , "trachea" : 0.015 , "oesophagus" : 0.015}
sigma_coeff = CoefficientFunction([sigma_coeff_[mat] for mat in mesh.GetMaterials()])

# create bounded X-elliptic bilinear form
a = BilinearForm(fes)
a += grad(v)*sigma_coeff*grad(u) * dx + 1e-8*sigma_coeff*u*v*dx
c = Preconditioner(a, "bddc")

# create bounded linear form
f = LinearForm(fes)

→ the changes of the linear form f are not in the code since they are not relevant ( I think) ←

I solved the equation with:

# Assemble equations
with TaskManager():
     a.Assemble()

# Solve PDE
gfu = GridFunction(fes)
inv = CGSolver(a.mat, c.mat, printrates=False, precision=1e-8)
gfu.vec.data = inv * f.vec

So far this should do the trick but if there are more elegant ways please let me know!
Now to my first problem/question:

I only want the access some points of the GridFunction (those points where my electrodes will be placed). So when I write:

pnts_electrode = mesh(0,0,1) 
val_gfu = gfu(pnts_electrode)

Is this the correct way to access the value of the GridFunction gfu at some point in the mesh that is closest to (0,0,1) ?
Somehow I get an error that my point is not in the mesh, although it should be.

Many thanks in advance!
Regards,
Markus

2

So I just discovered one mistake and came right to the next problem:

I defined my mesh point as

pnts_electrode = mesh(np.array([0,0,1]))

and not as

pnts_electrode = mesh(0,0,1)

After correcting this mistake, I got the following error message:

Exception has occurred: TypeError
call(): incompatible function arguments. The following argument types are supported:
1. (self: ngsolve.fem.CoefficientFunction, mip: ngsolve.fem.BaseMappedIntegrationPoint) → object
2. (self: ngsolve.fem.CoefficientFunction, x: float, y: Optional[float] = None, z: Optional[float] = None) → ngcomp::PointEvaluationFunctional
3. (self: ngsolve.fem.CoefficientFunction, arg0: numpy.ndarray[ngsolve.fem.MeshPoint]) → numpy.ndarray

Invoked with: <ngsolve.comp.GridFunctionD object at 0x1176d6840>, <ngsolve.fem.MeshPoint object at 0x10088c570>

So it seems to me I just have to pass the actual coordinates like:

val_gfu = gfu(0,0,1)

Unfortunately I then get a PointEvaluationFunctional and still not the actual value for these coordinates. Is there something like

val_gfu.Get()

I can do to access the value of the PointEvaluationFunctional ?

3

Hi Markus,

your first attempt with

pnts_electrode = mesh(0,0,1) 
val_gfu = gfu(pnts_electrode)

is the way to go. Are you sure that the point (0,0,1) is in the mesh ?

If you post the whole example (including the mesh) we can have a look.

Joachim

1 Like
4

So I made a small example and while doing it I think I came across the mistake:

from ngsolve import Mesh, H1, CoefficientFunction, BilinearForm, LinearForm, grad, dx, Parameter, Preconditioner, Draw
from ngsolve import ElementId, VOL, TaskManager, GridFunction, CGSolver
from netgen.csg import Cylinder, OrthoBrick, CSGeometry, Pnt
from math import pi, cos, sin

##### CREATE MESH #####

maxh = 0.1

# create cylinder along z-axis with center at (x = 0.2 ,y = 0) and height going from z = 0 to z = 1
cyl = Cylinder(Pnt(0.2,0,0), Pnt(0.2,0,1) , 1) * OrthoBrick(Pnt(-1, -1, 0), Pnt(1.5, 1.5, 1))
geo = CSGeometry()
geo.Add (cyl.mat('test')) 
mesh = geo.GenerateMesh(maxh = maxh)
mesh = Mesh(mesh)

###### FEM SETTINGS #####

# define Finite Element Space
fes = H1(mesh, order = '1')     
# define Trial and Test Space
u,v = fes.TnT()                 

# define conductivity on body
sigma_coeff_ = {"test" : 0.22}
sigma_coeff = CoefficientFunction([sigma_coeff_[mat] for mat in mesh.GetMaterials()])

# create bounded X-elliptic bilinear form
a = BilinearForm(fes)
a += grad(v)*sigma_coeff*grad(u) * dx + 1e-8*sigma_coeff*u*v*dx
c = Preconditioner(a, "bddc")

# create bounded linear form
f = LinearForm(fes)

##### ADD DIPOLE #####

# -- adding dipole by creating two sources with different sign close to each other -- 
# distance between the sources greater then maximal mesh size
max_mesh= maxh *1.1 
# dipole location
x = 0
y = 0
z = 0.5
# dipole direction - in this case along z-axis
beta = pi / 2
alpha = 0
# dipole magnitude
p = Parameter(1)

# -- netgen implementation --
# add first source term
spc = f.space
mp1 = spc.mesh(x,y,z)
ei1 = ElementId(VOL, mp1.nr)    
fel1 = spc.GetFE(ei1)
dnums1 = spc.GetDofNrs(ei1)
shape1 = fel1.CalcShape(*mp1.pnt)

for d1,s1 in zip(dnums1, shape1):
    f.vec[d1] += ( p.Get()*s1 )

# change coordinates of second source term
x1 = x + max_mesh * cos(beta) * cos(alpha)
y1 = y + max_mesh * cos(beta) * sin(alpha)
z1 = z + max_mesh * sin(beta)

# add second source term
mp2 = spc.mesh(x1, y1, z1)
ei2 = ElementId(VOL, mp2.nr)
fel2 = spc.GetFE(ei2)
dnums2 = spc.GetDofNrs(ei2)
shape2 = fel2.CalcShape(*mp2.pnt)

for d2, s2 in zip(dnums2, shape2):
    f.vec[d2] += ( - p.Get()*s2 )

##### SOLVE PDE #####

# Assemble equations
with TaskManager():
    a.Assemble()

# Solve PDE
gfu = GridFunction(fes)
inv = CGSolver(a.mat, c.mat, precision=1e-8)
gfu.vec.data = inv * f.vec

# Visualize Solution
Draw(gfu, mesh, 'potential')

##### DEFINE ELECTRODE POSITION AND CALCULATE SOLUTION #####

#mesh_elec = mesh(-0.8, 0, 0.5)      # does not work for this point
mesh_elec = mesh(-0.799, 0, 0.5)      # works for this point

gfu_elec = gfu(mesh_elec)
print(gfu_elec)

So it seems like although the point ( -0.8, 0, 0.5) is theoretically in the mesh, it is only at the very edge… what would be a proper way to work around this problem?

Thanks!
Markus

5

Hi Markus,

Since it is at a curved boundary the point is actually not in the mesh. I think the best option is to curve the mesh to the geometry then the point will be found:

mesh.Curve(3)

Best
Christopher

1 Like
6

Works perfectly, thanks a lot!

Best,
Markus