Get/Plot data from function over mesh

Question

Dear all.
I want to solve a time dependent problem where I have a function f1 which depends on a second function f2. Both functions are evaluated in the mesh coordinates.

def f2(x, y, f1):
  value = numpy.empty(1, dtype=float)     
  if f1(x,y) <= 0 + DOLFIN_EPS:
     value[0] = 1           
  elif f1(x,y) >= 0 - DOLFIN_EPS:
     value[0] = 2 
  return value

[...]

coor = mesh.coordinates()
number = len(coor)
f2_vector = numpy.empty(number, dtype=float)
for i in range(number):
   f2_vector[i] = f2(coor[i][0], coor[i][1], f1)

My problem is now the following: I want to plot the data given by function f2 against the mesh. Obviously, the data of f2 is not of the "correct" type for doing this (because f2 is not an expression etc but a vector). So here is my question: How can I plot the obtained data against the mesh? If I use an expression for f2 instead of a function, I cannot pass the (time-dependent function) f1.

I tried a quick workaround similar to the example in the "FEniCS Tutorial" by using a help function w and override the entries w.vector() but the order of the values is not the same (in general).

Cheers,
Mischa

Answer 1

Hi, as you say

the order of the values is not the same

The reason for this is that w.vector() entries are ordered by dofmap of function space of w, whereas you ordered f2_vector by coordinates. These orderings are not the same and also the number of mesh vertices in general does not equal the number of degrees of freedom. Finally, you might want to vectorize your code.

from dolfin import *
import numpy as np                                                               

mesh = UnitSquareMesh(50, 50)                                                    
V = FunctionSpace(mesh, "CG", 1)                                                 

u = Function(V)                                                                  
u_vector = u.vector()                                                            
dim = V.dim()                                                                    
N = mesh.geometry().dim()                                                        
dof_coordinates = V.dofmap().tabulate_all_coordinates(mesh).reshape(dim, N)      

x = dof_coordinates[:, 0]                                                        
y = dof_coordinates[:, 1]                                                        

# serial approach                                                                
def f2_s(x, y, f1):                                                              
  return 1 if f1(x, y) <= DOLFIN_EPS else 2                                      

def f1_s(x, y):                                                                  
  return sin(4*pi*x)*cos(4*pi*y)                                                 

tic()                                                                            
for i in xrange(dim):                                                            
  u_vector[i] = f2_s(x[i], y[i], f1_s)                                           
ser = toc(); print "Serial", ser                                                 

# vectorized approach                                                            
def f2(x, y, f1):                                                                
  return np.where(f1(x, y) <= DOLFIN_EPS, 1, 2)                                  

def f1(x, y):                                                                    
  return np.sin(4*pi*x)*np.cos(4*pi*y)                                           

tic()                                                                            
u_vector[:] = f2(x, y, f1)                                                       
vec = toc(); print "Vectorized", vec, "Speed up", int(ser/vec)  

Comments

Thank you very much, this solved my problem. I tried something quite similar using the dofs but could not get it working. So thank you again!

Cheers,
Mischa