Find point where solution equals 0

General
1

Hi all,
I have a working programme that allows me to solve a non-lineal problem on a 2D mesh. Now I would like to find the point where the solution is equal to 0. Is there a simple way to do this? In other words, is there a kind of inverse for func.eval(point, cells)?
Thanks!

2

This is not straightforward, as you would need to check this for every cell.

To get a good solution for this, you would need to be alot more specific about what you want to use these points for.

As I’m not sure if you are using DOLFIN or DOLFINx, here is a solution for identifying some cells:

from dolfin import *
import ufl
import numpy as np
mesh = UnitSquareMesh(10, 10)

expr = Expression("x[0]<0.2? 0: 0.5", degree=1)
V = FunctionSpace(mesh, "Lagrange", 1)
u = Function(V)
u.interpolate(expr)


W = FunctionSpace(mesh, "DG", 0)
w = Function(W)

check = ufl.conditional(ufl.eq(u, 0), 1, 0)
v = project(check, W)
cell_indices = np.argwhere(v.vector().get_local()==1)
for cell in cell_indices:
   print(Cell(mesh, cell).midpoint().array() )

For DOLFINx, I would do something like

import dolfinx
import ufl
import numpy as np
from mpi4py import MPI
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

x = ufl.SpatialCoordinate(mesh)
V =  dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
expr = dolfinx.fem.Expression(ufl.conditional(ufl.lt(x[0],0.2), 0, 0.5), V.element.interpolation_points())
u = dolfinx.fem.Function(V)
u.interpolate(expr)


W = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))
w = dolfinx.fem.Function(W)
check = dolfinx.fem.Expression(ufl.conditional(ufl.eq(u, 0), 1, 0), W.element.interpolation_points())
w.interpolate(check)

cell_indices = np.argwhere(w.x.array==1).astype(np.int32)
midpoints=dolfinx.mesh.compute_midpoints(mesh, mesh.topology.dim, cell_indices.reshape(-1))
print(midpoints)
1 Like
3

Hi dokken, thanks for your reply. I’m using dolfinx, I’m sorry I forgot to say that.
I tried to implement your suggestion but the ufl.conditional(ufl.eq(u, 0), 1, 0) complains (understandably) because my function is actually a vector field. Could you assist me in writing a component-wise auxiliary function? One that looks at the norm would also work.
Just in case it helps, in this problem I would expect the solution to vanish a) nowhere, b) at a point (maybe a handful of them), or c) on a curve, depending on the parameters.

4

For a component use u[i] where i is 0, 1 or 2 depending on the component.
For the norm replace u with dot(u, u).

This is still a fairly generic problem, as you would need to search through every cell and find the set of points where this happens.
What do you want this information for, post-processing or further computation?