This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

recording location of values above a threshold

0 votes

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...

asked Apr 24, 2017 by manucomp FEniCS Novice (150 points)

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.

1 Answer

+1 vote
 
Best answer

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)
answered Apr 25, 2017 by MiroK FEniCS Expert (80,920 points)
selected Apr 25, 2017 by manucomp

Excellent! Thank you so much

np.where

was indeed the correct syntax.

...