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?
Agree. I can fix it if I get a green light for that approach.
It's ok, so long the std::vector is not created within an inner loop.
It´s not.
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.
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());
}
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.
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);
Sure does:-) But line is 157 long, is it more readable?
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.
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.
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.
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.
Sounds good to me.
Leaving this for now as it will probably redundant following some upcoming changes for ghost cells in parallel.
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_orientationsis optionally provided by the mesh. When it is not provided a plainNULLpointer is passed to theFacet::adjacent_cells(see inAssembler.cpp). I think this logic should be handled insideAssembler.cppinstead of passed toadjacent_cells.