Assign MeshFunction/MeshValueCollection to MeshDomain

Question

Hi,
I'm using the Dolfin C++ interface.
I would like to store in a Mesh some information related to facet/cell markers, in such a way that when I do

File mesh_file ("my_mesh.xml");
mesh_file << mesh;

also the information about markers is saved.

I've seen that one of the suggested possibility is to use the mesh.domains().
Is this the right way of doing it?
If it is the case, how can I assign values to mesh.domains.markers(n) ??
At the moment I'm storing information in a MeshFunction/MeshValueCollection and I would like to assign it to mesh.domains().markers( dim);

A very simple example is:

int main()
{
  using namespace dolfin;
  UnitSquareMesh mesh(4,4);
  MeshFunction<uint> sub_domains(mesh, 2);
  sub_domains.set_all(0);
  for (CellIterator cell(mesh); !cell.end(); ++cell)
    {
      Point p = (*cell).midpoint();
      if (p.x() > 0.5)
    sub_domains[*cell] = 1;
    }
  MeshValueCollection<uint> my_sub_domains ;
  my_sub_domains =  sub_domains;

  //how can I do something like
      mesh.domains().markers(2) = my_sub_domains;

  File mesh_file ("my_mesh.xml");
  mesh_file << mesh;

}

Thanks

Marco

Answer 1

In DOLFIN version 1.2 you should be able to just do:

mesh.domains().markers(2) = sub_domains;

as what is returned by MeshDomains::markers is a MeshValueCollection.

In the development version you need to set them manually, at least at the moment. Something like:

for (CellIterator cell(mesh); !cell.end(); ++cell)
{
  Point p = (*cell).midpoint();
  if (p.x() > 0.5)
    mesh.domains().markers(2)[cell.index()] = 1;
}

Comments

Thanks for your answer.
I'm using version 1.2.0 and so the operator [ ] is not available.
As you suggested, code works with
*(mesh.domains().markers(2)) = my_sub_domains;
Thanks

Answer 2

You can probably avoid assigning markers to mesh.domains(). You can handle IO with these

void operator>> (Mesh& input);
void operator<< (const Mesh& output);

void operator>> (MeshFunction<std::size_t>& input)
void operator<< (const MeshFunction<std::size_t>& output)

void operator>> (MeshValueCollection<std::size_t>& input)
void operator<< (const MeshValueCollection<std::size_t>& output)

You simply need to save Mesh and MeshFunction to separate XML file.

Comments

I would like to keep them together.
I've seen that the meshes available on your website which are provided with boundary markers (like Calcium and Aneurysm) are just one single file, and I would like to do the same.
Thanks
Marco