nodal and vertex values in parallel fenics

9 months ago by
I want to parallel my fenics code. Everything is fine but the problem is I need to change a function value(k) which makes some problem.
As I found till now the mesh will be partitioned based on mesh vertex coordinate. However when I want to change function value I need to use nodal value which makes problem.
The mesh, function space and partition parameter:
mesh = BoxMesh(Point(x0, y0, z0), Point(x1, y1, z1), nx, ny, nz)
parameters["mesh_partitioner"] = "SCOTCH"#"ParMETIS"s
V = FunctionSpace(mesh, 'CG', 1)​
my k function:
tol = 1e-6
muz = 0.1
k  = Expression('x[2]>zf+tol ? k_a  : k_p' ,degree=0, zf=muz, k_a=20 , k_p=1, tol=tol)​
The k function:
def Kfunction(mesh,f,k,lw,u,V):
	k_p  = 20		
	v2d = vertex_to_dof_map(V)
	coordinates = mesh.coordinates()
	k = interpolate(k,V)
	nodal_values_k  = k.vector()
	for i, x in enumerate(coordinates):
		if( (x[2] <= f.muz) and (x[2] > f.muz-lw)):
			nodal_values_k[v2d[i]]  = k_p

	k.vector()[:]  = nodal_values_k

It gives me in line : "nodal_values_k[v2d[i]] = k_p"    following error:
IndexError: expected indices to be in [0..106]

which  is from v2d... means"vertex_to_dof_map" is not partitioned as desire....Is there any way to correct it?
In addition If I could write the Kfunction with vertex value as follow is there any way the vertex value_k to k without using dof to vertex map
def Kfunction(mesh,f,k,lw,u,V):
	k_p  = 20		
	coordinates = mesh.coordinates()
	k = interpolate(k,V)
	vertex_values_k = k.compute_vertex_values(mesh)
	for i, x in enumerate(coordinates):
		if( (x[2] <= f.muz) and (x[2] > f.muz-lw)):
			vertex_values_k[i] = k_p

	#??? how should I assign the vertex value_k to k witouth using dof to vertex map
Community: FEniCS Project
I found the it works perfectly :)

def Kfunction(mesh,f,k,lw,u,V):
    k_p  = 20
    k = interpolate(k,V)
    k_array  = k.vector().array()
    X = V.tabulate_dof_coordinates()
    X.resize((V.dim(), 3))	
    X.resize((V.dim(), 3))	
    for i, (x, v) in enumerate(zip(X, k_array)):
        if( (x[2] <= f.muz) and (x[2] > f.muz-lw)):
            k_array[i] = k_p
        k.vector()[:] = k_array
written 8 months ago by Nasibeh Hassanjanikhoshkroud  

1 Answer

8 months ago by
If I may add something to this monologue... ;)

Using the 'add_local' / 'set_local' methods are a more FEniCSian way of achieving what you are after. Consider the following MWE
from dolfin import *
import numpy as np

mesh = UnitSquareMesh(10,10)
V    = FunctionSpace(mesh, 'CG', 1)

k = interpolate(Expression('x[0]', degree = 1), V)
X = V.tabulate_dof_coordinates().reshape((-1,mesh.topology().dim()))

# Set value of k to 0.5 when x>=0.5
idcs = []; values = []
for i, dof_coords in enumerate(X):
    if dof_coords[0] >= 0.5:

## And set values in k vector
k.vector().set_local(np.asarray(values, dtype = np.float), np.asarray(idcs, dtype = np.intc))
# Apply changes

# Plot result for checking
outfile = File('test_set_local.pvd')
outfile << k​

Please login to add an answer/comment or follow this question.

Similar posts:
Search »