Binding with PYBIND11 module

1

Hi FEniCS community,

I am facing an issue regarding the binding of the following python and c++ script using pybind. Both codes work well, except the “PYBIND11 module” which creates an error. See the following MWE for FEniCS 2019.1.0. The purpose is to store the output of “extract_dof_component_map” from the c++ code to python.

I only have limited knowledge at the root level of FEniCS. Can someone provide a solution (syntax) to this?

from dolfin import *
import cppimport

compiled_cpp_module = cppimport.imp(‘delta_interpolation’)

mesh = UnitCubeMesh(21, 21, 21)
V = VectorFunctionSpace(mesh, ‘CG’, 2)

a = -1
compiled_cpp_module.extract_dof_component_map(dof_component_map, V, a)

/*
<%
from dolfin.jit.jit import dolfin_pc
setup_pybind11(cfg)
cfg[‘include_dirs’] = dolfin_pc[‘include_dirs’]
cfg[‘library_dirs’] = dolfin_pc[‘library_dirs’]
cfg[‘compiler_args’] = [‘-std=c++11’, ‘-DHAS_MPI’]
%>
*/

#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
#include <pybind11/numpy.h>
#include <pybind11/stl.h>

#include <dolfin/la/GenericVector.h>
#include <dolfin/function/Function.h>
#include <dolfin/function/FunctionSpace.h>
#include <dolfin/fem/GenericDofMap.h>
#include <dolfin/fem/FiniteElement.h>
#include <dolfin/common/RangedIndexSet.h>
#include <dolfin/geometry/BoundingBoxTree.h>
#include <dolfin/mesh/Edge.h>
#include <dolfin/mesh/Mesh.h>
#include <dolfin/mesh/Cell.h>
#include <dolfin/mesh/CellType.h>
#include <dolfin/mesh/MeshEntityIterator.h>
#include <dolfin/mesh/Vertex.h>
#include <dolfin/geometry/Point.h>
#include <dolfin/mesh/MeshEntity.h>

using namespace dolfin;
namespace py = pybind11;

void extract_dof_component_map(std::unordered_map<std::size_t,
std::size_t>& dof_component_map,
const FunctionSpace& V,
int* component)
{
// Extract sub dofmaps recursively and store dof to component map
if (V.element()->num_sub_elements() == 0)
{
std::unordered_map<std::size_t, std::size_t> collapsed_map;
std::shared_ptr dummy
= V.dofmap()->collapse(collapsed_map, *V.mesh());
(*component)++;
for (const auto &map_it : collapsed_map)
dof_component_map[map_it.second] = (*component);
}
else
{
for (std::size_t i = 0; i < V.element()->num_sub_elements(); ++i)
{
const std::vectorstd::size_t comp = {i};
std::shared_ptr Vs = V.extract_sub_space(comp);
extract_dof_component_map(dof_component_map, *Vs, component);
}
}}

PYBIND11_MODULE(delta_interpolation, m)
{

m.def(“extract_dof_component_map”, (void ()(std::unordered_map<std::size_t,
std::size_t>& , const FunctionSpace& V, int ))
&extract_dof_component_map);
// m.def(“extract_dof_component_map”, (py::object u0, py::object v, py::str abc){
// auto _u = u0.attr(“_cpp_object”);
// auto _v = v.attr(“_cpp_object”).cast<Function*>();
// auto _u0 = _u.cast<const Function*>();
// interpolate_delta(*_u0, *_v, abc.caststd::string());
// });
}

Reason why I am doing this (though not relevant to the problem statement) :

I have a user-defined interpolation algorithm (which I use for my FSI problem) very similar to the interpolate_non_matching_mesh in fenicstools, both of which use extract_dof_component_map. For my problem the fluid mesh is stationary and has close to 10 million DOF’s as a result of which running extract_dof_component_map (which is by itself a recursive function) is the bottleneck (consumes 80% of wall-time) of the entire FSI algorithm (even with mpi).
The idea is (since the mesh is stationary), the dof_component_map can be evaluated once at the beginning of time directly from the python code and later on be fed as an argument to my interpolation algorithm for further time steps. A more detailed understanding can be developed by going through the flow of “interpolate_non_matching_mesh” module in fenicstools