### nodal and vertex values in parallel fenics

545
views
0
8 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
return(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
return(k)
​
Community: FEniCS Project
I found the solution....now 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
return(k)
​
written 6 months ago by Nasibeh Hassanjanikhoshkroud

0
6 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:
idcs.append(i)
values.append(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
k.vector().apply('insert')

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