No typemap for adjacent_cells in Facet class

[Mikael Mortensen]

There function adjacent_cells does not seem to have a working typemap for its argument

adjacent_cells(const std::vector<std::size_t>* facet_orientation)

Any reason why this function uses a pointer to a const std::vector, and not the more common reference by the way?

[Johan Hake]

It would be easy to add. However I would suggest changing the signature to a reference. It looks like the pointer is used because the face_orientations is optionally provided by the mesh. When it is not provided a plain NULL pointer is passed to the Facet::adjacent_cells (see in Assembler.cpp). I think this logic should be handled inside Assembler.cpp instead of passed to adjacent_cells.

[Mikael Mortensen]

Agree. I can fix it if I get a green light for that approach.

[Prof Garth Wells]

It's ok, so long the std::vector is not created within an inner loop.

[Mikael Mortensen]

It´s not.

[Mikael Mortensen]

Next problem - there's no typemap for the return type

std::pair<const Cell, const Cell>

Could you please tell me how to implement this Johan:-) Alternatively it could return just the pair of cell-indices? It's only used by the Assembler/OpenMpAssembler.

[Johan Hake]

I would suggest return a tuple of two int. Return the passed by value Cell might end up with havoc due to garbage collection. The typemap is a simple one. Just add the following to dolfin/swig/typemaps/std_pair.i:

%typemap(out) std::pair<const dolfin::Cell, const dolfin::Cell>
{
  $result = Py_BuildValue("ii", $1.first.index(), $1.second.index());
}

[Mikael Mortensen]

Thanks Johan, that works.

For the original issue, now at least I think I understand why it was the way it was. A vector pointer must be used for facet_orientation because mesh.data().array() returns a reference to a vector. And facet_orientation must be declared before mesh.data().array() is called because we don't know if it exists (see initialization below).

Currently, my best solution to enable reference in adjacent_cells is to initialize like this (could also use a new std::vector<std::size_t> instead of &null)

  // Old code
  //const std::vector<std::size_t>* facet_orientation = NULL;
  //if (mesh.data().exists("facet_orientation", D - 1))
  //  facet_orientation = &(mesh.data().array("facet_orientation", D - 1));

  // New code
  std::vector<std::size_t> null;
  const std::vector<std::size_t>* facet_orientation = &null;
  if (mesh.data().exists("facet_orientation", D - 1))
    facet_orientation = &(mesh.data().array("facet_orientation", D - 1));

and then call like

//adjacent_cells(facet_orientation);
adjacent_cells(*facet_orientation);

Alternatively we could use the old signature and create a new typemap for the vector pointer. Let me know what you think.

[Johan Hake]

Could this work:

std::vector<std::size_t> null;
const std::vector<std::size_t>& facet_orientation = (mesh.data().exists("facet_orientation", D - 1) ? mesh.data().array("facet_orientation", D - 1) : null;
....
adjacent_cells(facet_orientation);

[Mikael Mortensen]

Sure does:-) But line is 157 long, is it more readable?

[Johan Hake]

Probably not but it avoids pointers. Your call.

Btw, I thought we had abandoned using Mesh::data. If facet_orientations was a speedy method returning a vector instead of a not so speedy map lookup, we could just postpone the access of the facet_orientation to the Facet::adjacent_cell method, which would be the natural way to do it as Facet keeps a reference to the Mesh.

[Mikael Mortensen]

Apparently there are still two MeshData, facet_orientation and parent_vertex_indices. I haven't really looked at the reasoning behind MeshData, but it's true, Facet keeps a reference to the mesh and a fast lookup of a vector item through that path would seem more natural.

[Prof Garth Wells]

Winding this back a step, why do you want to expose Facet::adjacent_cells? A reason to ask is that it will almost certainly be removed or changed once integration over interior facets is supported in parallel. This is planned for DOLFIN 1.4.

If it must be wrapped, I would prefer to not have it in DOLFIN 1.3 so that interface changes can be more easily made to support interior facet integration in v1.4.

On MeshData, I have said before that it should only be used for user data. Using it for library data is a confusing and bad design. I'd support any work to move data out of it.

[Mikael Mortensen]

Well, it was exposed, but not possible to use. So I thought it should be fixed. I don't really need it, I'm perfectly happy with using facet.entities from the Python side. Then again, I think the method has been there for a while and it makes perfect sense.

How about I create one overloaded adjacent_cells that take no arguments? This could be exposed to Python and the other with a facet_orientation pointer could be hidden.

[Johan Hake]

Sounds good to me.

[Prof Garth Wells]

Leaving this for now as it will probably redundant following some upcoming changes for ghost cells in parallel.