how to simply assign a value on a function on a unique vertex with its coordinate

Question

Hi everyone,
Very simple question here, but i didnt find a simple way to do it.

Let's say f_ is a function defined on a non-mixed space.

If i want to evaluate the function on a point, i do:

current_value = f_([x_coor,y_coor,z_coor])

How can i do simply the opposite operation (assign a new value on this coordinnate). Ideally, something like this:

f_([x_coor,y_coor,z_coor]) = new_value # Does not work

I only want to change the function on this particular point, that i know exist in the discrete space.

Regards,

Answer 1

Hi!

If your coordinate lies exactly on a vertex of your mesh this can be done easily like this:

f_.vector().array()[index] = new_value

Comments

I agree,
My issue is to get index from x_coor,y_coor,z_coor such as:

index = some function i do not know (x_coor,y_coor,z_coor)

Answer 2

Hi, the following minimal working example shows how to do this (actually, the answer is largely based on this answer, so credits to MiroK ;) ):

from dolfin import * 
import numpy as np

mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, 'CG', 2)
v = Function(V)

dofmap = V.dofmap()
nvertices = mesh.ufl_cell().num_vertices()

# Set up a vertex_2_dof list
indices = [dofmap.tabulate_entity_dofs(0, i)[0] for i in range(nvertices)]

vertex_2_dof = dict()
[vertex_2_dof.update(dict(vd for vd in zip(cell.entities(0),
                                        dofmap.cell_dofs(cell.index())[indices])))
                        for cell in cells(mesh)]

# Get the vertex coordinates
X = mesh.coordinates()

# Set the vertex coordinate you want to modify
xcoord, ycoord = 0.5, 0.5

# Find the matching vertex (if it exists)
vertex_idx = np.where((X == (xcoord,ycoord)).all(axis = 1))[0] 
if not vertex_idx:
    print 'No matching vertex!'
else:
    vertex_idx = vertex_idx[0]
    dof_idx = vertex_2_dof[vertex_idx]
    print dof_idx
    v.vector()[dof_idx] = 1.

plot(v)
interactive()

Comments

Hi, it's great, thanks!
However, I am wandering if such a function is time efficient.
If i understand correctly,

vertex_idx = np.where((X == (xcoord,ycoord)).all(axis = 1))[0]

will compare all elements of X until it find the corresponding xcoord,ycoord. It can takes some times (from a CPU time perspective) to get it. My issue is that this operation (in my program) is called a lot of time, in a large mesh.

To avoid any np.where or equivalent find operation, i have finally used an equivalent method ( for my particular problem): i have basically defined an numpy array so that array[X][Y][Z] gives the value of the function f_ (that is no more used at all) at the coordinates x,y,z (the only operation here consists in transform the coordinates x in the corresponding index integer through a simple X=int(round(x/mesh_step)).
Since i am working with an array, and that i know (i) the coordinates where i want to change the value and (ii) the new value, it is know both simple and fast:

# Convert coordinates in index integer for the array
X=int(round(coor_x/mesh_step_along_x))
Y=int(round(coor_y/mesh_step_along_y))
Z=int(round(coor_z/mesh_step_along_z))
# Current value, that corresponds at the information i store at the node of coordinates coor_x coor_y coor_z
current_value = array[X][Y][Z]
# New value
array[X][Y][Z] = new_value

In my problem, i only want to be able to save/modify an information locally at the vertex locations, then this solution is enough.
Regards,

---

Well, I see which direction you want to go. And although your approach might be faster, you have to keep in mind that it is very non-generic (only works for regular meshes). Moreover, you need to extract the vertex index from X,Y,Z (for this, you have to know how the vertices are ordered). Once you know the vertex index, you cannot avoid using something like the vert_2_dof dictionary as proposed in my answer in order to get the dof index and set the corresponding entry in the function object.

---

I agree, it's not generic. I am working on a regular mesh made of cubic element (electrode microstructrure obtained through X-ray tomography), then it works well.

Basically, i check a field value on a particulate coordinates coor_x, coor_y ,coor_z (coordinates i already know).

current_value = f_([x_coor,y_coor,z_coor])

Then, i convert this coordinate in an integer that indicated the position (not a classic fenics index). For example, if coor_x=9.3 and each voxel (cubic element) has a lenght of 0.3, then X=9.3/0.3=31. I use round to get ridd off numerical truncation error and then int.

Depending on the current_value, i then assign array[X][Y][Z]=new_value

That's simply a way to store some additional information without dealing with fenics index. And as you said, it requires a regular mesh.

---

Allright, then just find the way the vertices are ordered in the (regular) mesh, get the vertex index based on the 'approximate' indices X,Y,Z, and use the vert_2_dof dictionary to find the dof index.
That should do the trick, I guess?