Subdomains in parallel

Question

Dear all,

I have the following reaction term in my bilinear form that uses subdomains:

... # some things

psi = TrialFunction(V)  # solution
v   = TestFunction(V)
Sigma_a  = Function(V0)

sigma_a_values = [sigma1,sigma2]

for cell_no in range(len(subdomains.array())):
    dofs_0 = V0.dofmap().cell_dofs(cell_no)
    subdomain_no = subdomains.array()[cell_no]
    Sigma_a.vector()[cell_no] = sigma_a_values[subdomain_no]

Sigma_a0 = interpolate(Sigma_a,V)

...  # some other things

a += Sigma_a0*psi*v*dx(mesh)

...  # some other things

If I run this code, I get different results in serial or in parallel...

If I replace the bilinear form term mentioned above with:

a += sigma1*psi*v*dx(mesh)

it gives the same results in serial and in parallel.

As my subdomains have been working fine in serial, I am assuming that the problem is caused by the use of subdomains in parallel.

I would think that this is related to this question but using:

Sigma_a.vector().apply('')

did not solve the problem...

Does anyone have an idea what the problem might be?

Thanks a lot!

Answer 1

This should work just fine if using the development version (see here):

dofmap = V0.dofmap()
dof_x = dofmap.tabulate_all_coordinates(mesh).reshape((-1, 2))
first_dof, last_dof = dofmap.ownership_range()  # U.local_size()

new_values = np.zeros(last_dof - first_dof)
for local_cell_no in range(len(new_values)):
  subdomain_no = subdomains.array()[local_cell_no]
  new_values[local_cell_no] = sigma_a_values[subdomain_no]

Sigma_a.vector().set_local(new_values)
Sigma_a.vector().apply('insert')