recording location of values above a threshold

Question

I would like to track the location of high values of an evolving scalar field, say v, in a given function space

Vh = FunctionSpace(mesh,"Lagrange",1)
v = Function(Vh)

so I create an index set related to the sought threshold (0.75)

indexA = list(set(v.vector()>0.75))

and create an array to store the values associated to these indexes

valueA = v.vector().array()[indexA]

Then, if I try to see the results (using Vh.dofmap().tabulate_all_coordinates(mesh)), apparently the listing does not work as expected. Perhaps this is due to the syntax and I'm missing something? Any help would be appreciated...

Comments

How does it not work as expected? How do you 'see the results'? Keep in mind that the array you've created is a copy of the one for v. If you want a Function w in Vh with the values you've saved, you need to use w.vector().set_local() (after declaring w of course). Please post more details

---

Thank you for your prompt answer!

A minimal working example is as follows.

from dolfin import *

# Create mesh and define function space
mesh = RectangleMesh(Point(0,0),Point(1,1),10,10)
Mh = FunctionSpace(mesh, "Lagrange", 1)

# Tabulating coordinates
dof_coord = Mh.dofmap().tabulate_all_coordinates(mesh)   
dof_coord.resize((Mh.dim(),mesh.geometry().dim()))
dof_x = dof_coord[:, 0]
dof_y = dof_coord[:, 1]     

# Create a function with maximum near the center
f = Expression("1.0/exp(pow(x[0] - 0.5,2)+pow(x[1]-0.5,2))")
v = Function(Mh); v=interpolate(f,Mh)
plot(v, interactive=True)

# Getting the indexes of the values above a certain threshold 
indexes = list(set(v.vector()>0.9))
values  = v.vector().array()[indexes]

# Printing the value and locations 
for j in range(len(values)):
     print 'v(%2g,%2g) = %g' % (dof_x[j], dof_y[j], values[j])

And as a result I get

v(0, 1) = 0.606531;   v(0,0.9) = 0.66365

which do not correspond with the expected values (above 0.9) and located near the center of the domain... I suspect the issue is the definition of the indexes, but haven't been able to resolve the problem.

Answer 1

Hi, consider the following

from dolfin import *
import numpy as np

# Create mesh and define function space
mesh = RectangleMesh(Point(0,0),Point(1,1),10,10)
Mh = FunctionSpace(mesh, "Lagrange", 1)

# Tabulating coordinates
# dof_coord = Mh.dofmap().tabulate_all_coordinates(mesh)   
dof_coord = Mh.tabulate_dof_coordinates().reshape((-1, 2))  # 2016.2.0

# Create a function with maximum near the center
f = Expression("1.0/exp(pow(x[0] - 0.5,2)+pow(x[1]-0.5,2))", degree=3)
v = Function(Mh); v=interpolate(f,Mh)
plot(v, interactive=True)

# Getting the indexes of the values above a certain threshold 
indexes = np.where(v.vector()>0.9)[0]
values  = v.vector().array()[indexes] 

# Printing the value and locations 
for j, dof_index in enumerate(indexes):
    x, y = dof_coord[dof_index]
    val = values[j]
    val0 = v(x, y)
    print 'v(%2g,%2g) = %g [%g]' % (x, y, val, val0)

Comments

Excellent! Thank you so much

np.where

was indeed the correct syntax.