Change some values of a function over a mesh from a function defined on a BoundaryMesh

Question

Dear all,

I am trying to change the values of a function on some vertices of the mesh, but I am finding errors. To be more clear,

1. I first define a mesh over a domain and then I have create a submesh using BoundaryMesh.
2. I define a function over the main mesh and I set it to zero.
3. I try to change only its values over the boundary, based on another function, defined on the boundarymesh, using boundarymesh.entity_map(0) to have a map between vertices on the boundary and vertices on the domain.

However this does not work. I am copying part of the code below.
What am I doing wrong?

Thanks!

---

from dolfin import *
from fenics import *
from mshr import *

Omega = Circle(Point(0,0), 2)
mesh = generate_mesh(Omega,10)
bmesh = BoundaryMesh(mesh, 'exterior')
vb_2_v = bmesh.entity_map(0)

V = FunctionSpace(mesh, 'P', 1)
V_Boundary = FunctionSpace(bmesh,'P', 1)

f = Expression('x[1] < Eps + tol && x[1] > - Eps - tol ? 1 : 0', degree=1, tol=1E-14, Eps = 0.3)
fb = interpolate(f, V_Boundary)

F_domain = interpolate(Constant(0.0), V)
for i in range(bmesh.num_vertices()):
j = vb_2_v.array()[i]
F_domain.vector()[j] = fb.vector().array()[i]
interpolate(F_domain, V)
plot(F_domain, interactive = True)

Comments

The ordering of the vertices of the mesh does not necessary coincide with the ordering of the degree of freedom in the function space.

---

To read more about finding/extracting the dof indices at the boundary: see here

Answer 1

Hi, consider

from dolfin import *
import numpy as np

mesh = UnitCubeMesh(10, 10, 10)
bmesh = BoundaryMesh(mesh, 'exterior')

V = FunctionSpace(mesh, 'CG', 1)
Vb = FunctionSpace(bmesh, 'CG', 1)

F = Expression('2*(x[0]+x[1]) + 4*(x[1]-x[2]) - x[2]', degree=1)
fb = interpolate(F, Vb)

dofb_vb = np.array(dof_to_vertex_map(Vb), dtype=int)
vb_v = np.array(bmesh.entity_map(0), dtype=int)
v_dof = np.array(vertex_to_dof_map(V), dtype=int)

f = Function(V)
array = f.vector().array()
# Compose the maps dof of Vb --> vertices of Vb --> vertices of V --> dof of V
array[v_dof[vb_v[dofb_vb]]] = fb.vector().array()
f.vector()[:] = array

print assemble(inner(f-F, f-F)*ds(domain=mesh)) 

plot(f)
interactive() 

Comments

Thank you very much!

Answer 2

I agree with umberto: This is my way to do something like you want to do:

import numpy as np
dim = 1 # in my case 3 for 'CG' VectorFunctionSpace
dof_coords_V = V.tabulate_dof_coordinates().reshape(-1, dim)
dof_coords_V_Boundary = V_Boundary.tabulate_dof_coordinates().reshape(-1, dim)

Now you can do something like this in the 1D case, of course there are faster ways to get this information:

a1 = []
a2 = []
for j in range(len(dof_coords_V )):
    for k in range(len(dof_coords_V_Boundary)):
        if np.isclose(dof_coords_V[j], dof_coords_V_Boundary[k]) is True:
            a1.extend([j])
            a2.extend([k])

a1 and a2 contain the dofs which belong to the same coordinate of both meshes. Finally:

f1 = Function(V) # your "zero" function on the whole domain
f2 = Function(V_Boundary) # your interpolated boundary function containg the values
f1.vector()[a1] = f2.vector()[a2]

I extracted this code snippet from my solution, it should work fine in serial!

Comments

Thank you! It works! It can be solved faster in the way presented by MiroK. Thank you for your help!