Subdomains with different properties not working in parallel

Question

Hello,

I am trying to divide my mesh in subdomains (in c++) and assign different material properties to them. The problem that I am trying to run is a demo of the FEniCS Solid Mech. app.

I attempted this code in c++, which runs in serial but fails in parallel. This is minimal, but I can post the full code if someone is interested. Can anyone please advise on how to accomplish the purpose in question, while running in parallel.

  auto V0 = std::make_shared<DG0::FunctionSpace>(mesh);

  MeshFunction<std::size_t> subdomains(mesh, mesh.topology().dim(), true);

  Omega0 Omega_0;  // This subdomain is defined as (x[1] <= 0.5) for a Gmsh unitsquare
  Omega_0.mark(subdomains, 0);

  Omega1 Omega_1;  // This subdomain is defined as (x[1] >= 0.5) for a Gmsh unitsquare
  Omega_1.mark(subdomains, 1);

  // Material parameters -->  mass density
  Function rho(V0);
  Array<double> rho_values(2);
  rho_values[0] = 0.5;
  rho_values[1] = 2;

  std:size_t cell_no, subdomain_no;
  for (CellIterator c(mesh); !c.end(); ++c)
  {
   cell_no = c->index();                                 // Cell number
   subdomain_no = subdomains.values()[cell_no];         
   // Number subdomain identifier that owns the cell number " cell_no " 
   rho.vector()->setitem(cell_no, rho_values[subdomain_no]);  

   std::cout << "" << std::endl;
   std::cout << "- sub_domain_number: " << subdomains.values()[cell_no] << std::endl;
   std::cout << "- cell_number: " << cell_no << std::endl;
   std::cout << "- rho_in_cell: " << rho.vector()->getitem(cell_no) << std::endl;
  }

  plot(mesh);
  interactive();
  plot(subdomains);
  interactive();  
  plot(rho);
  interactive();  

I think the problem with this code in parallel may be the loop where it iterates over cells. This will not work in parallel since each process has only a subset of the cells.

Thanks for your help in advance!

Answer 1

What exactly happens when it fails? Does the code crash or just give an unexpected result?

You can't be sure that the DoFs of rho line up with the cell indices. It may just be coincidence that it works in serial.

The DofMap may be completely different to the cell ordering.

Why not just project or interpolate a piecewise expression into rho instead? Simple Python example:

from dolfin import *

mesh = UnitSquareMesh(32, 32)
V0 = FunctionSpace(mesh, 'DG', 0)
rho = interpolate(Expression('(x[1] < 0.5) ? 0.5 : 2.0'), V0)
plot(rho, interactive=True)

Comments

Hi nate,

Thank you very much for your reply. It runs even when run in parallel but does not produce correct result, as I observe from the plots at the end. The plots are correct when run in serial though. It does not crash in any case. I can post the full code if you want to run it. Can the command "interpolate(Expression('(x[1] < 0.5) ? 0.5 : 2.0'), V0)" be modified for 3 or 4 layers in c++?

Thanks!!

---

Hi Nate,

I tried your approach in C++ and it worked fr a simple problem, but I got stuck where I actually need it - in my problem. I have posted the Question here.

https://groups.google.com/forum/#!topic/fenics-support/FNu7ZTXJwEY

I can post it in this forum if you want.

Please tell me if you have a solution to this.

Thanks a lot!!