IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Mesh const& mesh,
                                                     std::string const& name)
{
    if (!mesh.getProperties().existsPropertyVector<char>(
            "IntegrationPointMetaData"))
    {
        OGS_FATAL(
            "Integration point data '%s' is present in the vtk field "
            "data but the required 'IntegrationPointMetaData' array "
            "is not available.",
            name.c_str());
    }
    auto const& mesh_property_ip_meta_data =
        *mesh.getProperties().template getPropertyVector<char>(
            "IntegrationPointMetaData");

    if (mesh_property_ip_meta_data.getMeshItemType() !=
        MeshLib::MeshItemType::IntegrationPoint)
    {
        OGS_FATAL("IntegrationPointMetaData array must be field data.");
    }

    // Find the current integration point data entry and extract the
    // meta data.
    auto const ip_meta_data = extractIntegrationPointMetaData(
        json::parse(mesh_property_ip_meta_data.begin(),
                    mesh_property_ip_meta_data.end()),
        name);

    return ip_meta_data;
}

void writeToFile(std::string const& id_area_fname, std::string const& csv_fname,
    std::vector<std::pair<std::size_t, double>> const& ids_and_areas,
    std::vector<MeshLib::Node*> const& mesh_nodes)
{
    std::ofstream ids_and_area_out(id_area_fname);
    if (!ids_and_area_out) {
        OGS_FATAL("Unable to open the file \"%s\" - aborting.", id_area_fname.c_str());
    }
    std::ofstream csv_out(csv_fname);
    if (!csv_out) {
        OGS_FATAL("Unable to open the file \"%s\" - aborting.", csv_fname.c_str());
    }

    ids_and_area_out << std::setprecision(20);
    csv_out << std::setprecision(20);

    ids_and_area_out << ids_and_areas[0].first << " " << ids_and_areas[0].second;
    csv_out << "ID x y z area node_id\n"; // CSV header
    csv_out << 0 << " " << *mesh_nodes[ids_and_areas[0].first]
            << ids_and_areas[0].second << " " << ids_and_areas[0].first;
    for (std::size_t k(1); k<ids_and_areas.size(); k++) {
        ids_and_area_out << "\n"
            << ids_and_areas[k].first << " " << ids_and_areas[k].second;
        csv_out << "\n" << k << " " << *mesh_nodes[ids_and_areas[k].first]
            << ids_and_areas[k].second << " " << ids_and_areas[k].first;
    }
    ids_and_area_out << "\n";
    csv_out << "\n";
}

static void checkJacobianDeterminant(const double detJ,
                                     MeshLib::Element const& element)
{
    if (detJ > 0)  // The usual case
        return;

    if (detJ < 0)
    {
        ERR("det J = %g is negative for element %d.", detJ, element.getID());
#ifndef NDEBUG
        std::cerr << element << "\n";
#endif  // NDEBUG
        OGS_FATAL("Please check the node numbering of the element.");
    }

    if (detJ == 0)
    {
        ERR("det J is zero for element %d.", element.getID());
#ifndef NDEBUG
        std::cerr << element << "\n";
#endif  // NDEBUG
        OGS_FATAL(
            "Please check whether:\n"
            "\t the element nodes may have the same coordinates,\n"
            "\t or the coordinates of all nodes are not given on the x-axis "
            "for a 1D problem or in the xy-plane in the 2D case.");
    }
}

std::unique_ptr<PythonBoundaryCondition> createPythonBoundaryCondition(
    BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
    NumLib::LocalToGlobalIndexMap const& dof_table, std::size_t bulk_mesh_id,
    int const variable_id, int const component_id,
    unsigned const integration_order, unsigned const shapefunction_order,
    unsigned const global_dim)
{
    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
    config.checkConfigParameter("type", "Python");

    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Python__bc_object}
    auto const bc_object = config.getConfigParameter<std::string>("bc_object");
    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Python__flush_stdout}
    auto const flush_stdout = config.getConfigParameter("flush_stdout", false);

    // Evaluate Python code in scope of main module
    pybind11::object scope =
        pybind11::module::import("__main__").attr("__dict__");

    if (!scope.contains(bc_object))
        OGS_FATAL(
            "Function `%s' is not defined in the python script file, or there "
            "was no python script file specified.",
            bc_object.c_str());

    auto* bc = scope[bc_object.c_str()]
                   .cast<PythonBoundaryConditionPythonSideInterface*>();

    if (variable_id >= static_cast<int>(dof_table.getNumberOfVariables()) ||
        component_id >= dof_table.getNumberOfVariableComponents(variable_id))
    {
        OGS_FATAL(
            "Variable id or component id too high. Actual values: (%d, %d), "
            "maximum values: (%d, %d).",
            variable_id, component_id, dof_table.getNumberOfVariables(),
            dof_table.getNumberOfVariableComponents(variable_id));
    }

    // In case of partitioned mesh the boundary could be empty, i.e. there is no
    // boundary condition.
#ifdef USE_PETSC
    // This can be extracted to createBoundaryCondition() but then the config
    // parameters are not read and will cause an error.
    // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
    // subtree and move the code up in createBoundaryCondition().
    if (boundary_mesh.getDimension() == 0 &&
        boundary_mesh.getNumberOfNodes() == 0 &&
        boundary_mesh.getNumberOfElements() == 0)
    {
        return nullptr;
    }
#endif  // USE_PETSC

    return std::make_unique<PythonBoundaryCondition>(
        PythonBoundaryConditionData{
            bc, dof_table, bulk_mesh_id,
            dof_table.getGlobalComponent(variable_id, component_id),
            boundary_mesh},
        integration_order, shapefunction_order, global_dim, flush_stdout);
}

std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
createNonlinearSolver(GlobalLinearSolver& linear_solver,
                      BaseLib::ConfigTree const& config)
{
    //! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__type}
    auto const type = config.getConfigParameter<std::string>("type");
    //! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__max_iter}
    auto const max_iter = config.getConfigParameter<unsigned>("max_iter");

    if (type == "Picard") {
        auto const tag = NonlinearSolverTag::Picard;
        using ConcreteNLS = NonlinearSolver<tag>;
        return std::make_pair(
            std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
    }
    if (type == "Newton")
    {
        //! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__damping}
        auto const damping = config.getConfigParameter<double>("damping", 1.0);
        if (damping <= 0)
        {
            OGS_FATAL(
                "The damping factor for the Newon method must be positive, got "
                "%g.",
                damping);
        }
        auto const tag = NonlinearSolverTag::Newton;
        using ConcreteNLS = NonlinearSolver<tag>;
        return std::make_pair(
            std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
            tag);
    }
    OGS_FATAL("Unsupported nonlinear solver type");
}

ConfigTreeTopLevel
makeConfigTree(const std::string& filepath, const bool be_ruthless,
               const std::string& toplevel_tag)
{
    ConfigTree::PTree ptree;

    // note: Trimming whitespace and ignoring comments is crucial in order
    //       for our configuration tree implementation to work!
    try {
        read_xml(filepath, ptree,
                 boost::property_tree::xml_parser::no_comments |
                 boost::property_tree::xml_parser::trim_whitespace);
    } catch (boost::property_tree::xml_parser_error e) {
        OGS_FATAL("Error while parsing XML file `%s' at line %lu: %s.",
                  e.filename().c_str(), e.line(), e.message().c_str());
    }

    DBUG("Project configuration from file \'%s\' read.", filepath.c_str());

    if (auto child = ptree.get_child_optional(toplevel_tag)) {
        return ConfigTreeTopLevel(filepath, be_ruthless, std::move(*child));
    } else {
        OGS_FATAL("Tag <%s> has not been found in file `%s'.",
                  toplevel_tag.c_str(), filepath.c_str());
    }
}

std::unique_ptr<ParameterBase> createMeshElementParameter(
    BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh)
{
    //! \ogs_file_param{parameter__type}
    config.checkConfigParameter("type", "MeshElement");
    //! \ogs_file_param{parameter__MeshElement__field_name}
    auto const field_name = config.getConfigParameter<std::string>("field_name");
    DBUG("Using field_name %s", field_name.c_str());

    if (!mesh.getProperties().hasPropertyVector(field_name)) {
        OGS_FATAL("The required property %s does not exists in the mesh.",
                  field_name.c_str());
    }

    // TODO other data types than only double
    auto const& property =
        mesh.getProperties().getPropertyVector<double>(field_name);
    if (!property) {
        OGS_FATAL("The mesh property `%s' is not of the requested type.",
                  field_name.c_str());
    }

    if (property->getMeshItemType() != MeshLib::MeshItemType::Cell) {
        OGS_FATAL("The mesh property `%s' is not an element property.",
                  field_name.c_str());
    }

    return std::unique_ptr<ParameterBase>(
        new MeshElementParameter<double>(*property));
}

void LocalLinearLeastSquaresExtrapolator::calculateResiduals(
    const unsigned num_components,
    ExtrapolatableElementCollection const& extrapolatables,
    const double t,
    GlobalVector const& current_solution,
    LocalToGlobalIndexMap const& dof_table)
{
    auto const num_element_dof_result = static_cast<GlobalIndexType>(
        _dof_table_single_component.size() * num_components);

    if (!_residuals || _residuals->size() != num_element_dof_result)
    {
#ifndef USE_PETSC
        _residuals.reset(new GlobalVector{num_element_dof_result});
#else
        _residuals.reset(new GlobalVector{num_element_dof_result, false});
#endif
    }

    if (static_cast<std::size_t>(num_element_dof_result) !=
        extrapolatables.size() * num_components)
    {
        OGS_FATAL("mismatch in number of D.o.F.");
    }

    auto const size = extrapolatables.size();
    for (std::size_t i = 0; i < size; ++i)
    {
        calculateResidualElement(i, num_components, extrapolatables, t,
                                 current_solution, dof_table);
    }
    MathLib::LinAlg::finalizeAssembly(*_residuals);
}

std::unique_ptr<ParameterBase> createCurveScaledParameter(
    std::string const& name,
    BaseLib::ConfigTree const& config,
    std::map<std::string,
             std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
        curves)
{
    //! \ogs_file_param{prj__parameters__parameter__type}
    config.checkConfigParameter("type", "CurveScaled");

    //! \ogs_file_param{prj__parameters__parameter__CurveScaled__curve}
    auto curve_name = config.getConfigParameter<std::string>("curve");
    DBUG("Using curve %s", curve_name.c_str());

    auto const curve_it = curves.find(curve_name);
    if (curve_it == curves.end())
        OGS_FATAL("Curve `%s' does not exists.", curve_name.c_str());

    auto referenced_parameter_name =
        //! \ogs_file_param{prj__parameters__parameter__CurveScaled__parameter}
        config.getConfigParameter<std::string>("parameter");
    DBUG("Using parameter %s", referenced_parameter_name.c_str());

    // TODO other data types than only double
    return std::make_unique<CurveScaledParameter<double>>(
        name, *curve_it->second, referenced_parameter_name);
}

MeshLib::Element* BoundaryElementsAlongPolyline::modifyEdgeNodeOrdering(
    const MeshLib::Element& edge, const GeoLib::Polyline& ply,
    const std::vector<std::size_t>& edge_node_distances_along_ply,
    const std::vector<std::size_t>& node_ids_on_poly) const
{
    // The first node of the edge should be always closer to the beginning of
    // the polyline than other nodes.
    if (edge_node_distances_along_ply.front() >
            edge_node_distances_along_ply.back() ||
        (ply.isClosed() &&
         edge_node_distances_along_ply.back() == node_ids_on_poly.size() - 1))
    {  // Create a new element with reversed local node index
        auto new_nodes = new MeshLib::Node*[edge.getNumberOfNodes()];
        if (auto const* e = dynamic_cast<MeshLib::Line const*>(&edge))
        {
            new_nodes[0] = const_cast<MeshLib::Node*>(e->getNode(1));
            new_nodes[1] = const_cast<MeshLib::Node*>(e->getNode(0));
        }
        else if (auto const* e = dynamic_cast<MeshLib::Line3 const*>(&edge))
        {
            new_nodes[0] = const_cast<MeshLib::Node*>(e->getNode(1));
            new_nodes[1] = const_cast<MeshLib::Node*>(e->getNode(0));
            new_nodes[2] = const_cast<MeshLib::Node*>(e->getNode(2));
        }
        else
            OGS_FATAL("Not implemented for element type %s", typeid(edge).name());

        return edge.clone(new_nodes, edge.getID());
    }

    // Return the original edge otherwise.
    return const_cast<MeshLib::Element*>(&edge);
}

ProcessVariable& findProcessVariable(
    std::vector<ProcessVariable> const& variables,
    BaseLib::ConfigTree const& pv_config, std::string const& tag)
{
    // Find process variable name in process config.
    //! \ogs_file_special
    std::string const name = pv_config.getConfigParameter<std::string>(tag);

        // Find corresponding variable by name.
    auto variable = std::find_if(variables.cbegin(), variables.cend(),
                                 [&name](ProcessVariable const& v)
                                 {
                                     return v.getName() == name;
                                 });

    if (variable == variables.end())
    {
        OGS_FATAL(
            "Could not find process variable '%s' in the provided variables "
            "list for config tag <%s>.",
            name.c_str(), tag.c_str());
    }
    DBUG("Found process variable \'%s\' for config tag <%s>.",
         variable->getName().c_str(), tag.c_str());

    // Const cast is needed because of variables argument constness.
    return const_cast<ProcessVariable&>(*variable);
}

std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
createNonlinearSolver(GlobalLinearSolver& linear_solver,
                      BaseLib::ConfigTree const& config)
{
    using AbstractNLS = NonlinearSolverBase;

    //! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__type}
    auto const type = config.getConfigParameter<std::string>("type");
    //! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__tol}
    auto const tol = config.getConfigParameter<double>("tol");
    //! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__max_iter}
    auto const max_iter = config.getConfigParameter<unsigned>("max_iter");

    if (type == "Picard")
    {
        auto const tag = NonlinearSolverTag::Picard;
        using ConcreteNLS = NonlinearSolver<tag>;
        return std::make_pair(std::unique_ptr<AbstractNLS>(new ConcreteNLS{
                                  linear_solver, tol, max_iter}),
                              tag);
    }
    else if (type == "Newton")
    {
        auto const tag = NonlinearSolverTag::Newton;
        using ConcreteNLS = NonlinearSolver<tag>;
        return std::make_pair(std::unique_ptr<AbstractNLS>(new ConcreteNLS{
                                  linear_solver, tol, max_iter}),
                              tag);
    }
    OGS_FATAL("Unsupported nonlinear solver type");
}

PhaseFieldProcess<DisplacementDim>::PhaseFieldProcess(
    MeshLib::Mesh& mesh,
    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
    unsigned const integration_order,
    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
        process_variables,
    PhaseFieldProcessData<DisplacementDim>&& process_data,
    SecondaryVariableCollection&& secondary_variables,
    NumLib::NamedFunctionCaller&& named_function_caller,
    bool const use_monolithic_scheme)
    : Process(mesh, std::move(jacobian_assembler), parameters,
              integration_order, std::move(process_variables),
              std::move(secondary_variables), std::move(named_function_caller),
              use_monolithic_scheme),
      _process_data(std::move(process_data))
{
    if (use_monolithic_scheme)
    {
        OGS_FATAL(
            "Monolithic scheme is not implemented for the PhaseField process.");
    }

    _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
        mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
}

std::unique_ptr<BHEBottomDirichletBoundaryCondition>
createBHEBottomDirichletBoundaryCondition(
    std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices)
{
    DBUG("Constructing BHEBottomDirichletBoundaryCondition.");

    // In case of partitioned mesh the boundary could be empty, i.e. there is no
    // boundary condition.
#ifdef USE_PETSC
    // For this special boundary condition the boundary condition is not empty
    // if the global indices are non-negative.
    if (in_out_global_indices.first < 0 && in_out_global_indices.second < 0)
    {
        return nullptr;
    }
    // If only one of the global indices (in or out) is negative the
    // implementation is not valid.
    if (in_out_global_indices.first < 0 || in_out_global_indices.second < 0)
    {
        OGS_FATAL(
            "The partition cuts the BHE into two independent parts. This "
            "behaviour is not implemented.");
    }
#endif  // USE_PETSC

    return std::make_unique<BHEBottomDirichletBoundaryCondition>(
        std::move(in_out_global_indices));
}

void ConvergenceCriterionPerComponentResidual::checkDeltaX(
    const GlobalVector& minus_delta_x, GlobalVector const& x)
{
    if ((!_dof_table) || (!_mesh))
    {
        OGS_FATAL("D.o.f. table or mesh have not been set.");
    }

    for (unsigned global_component = 0; global_component < _abstols.size();
         ++global_component)
    {
        // TODO short cut if tol <= 0.0
        auto error_dx = norm(minus_delta_x, global_component, _norm_type,
                             *_dof_table, *_mesh);
        auto norm_x =
            norm(x, global_component, _norm_type, *_dof_table, *_mesh);

        INFO(
            "Convergence criterion, component %u: |dx|=%.4e, |x|=%.4e, "
            "|dx|/|x|=%.4e",
            error_dx, global_component, norm_x,
            (norm_x == 0. ? std::numeric_limits<double>::quiet_NaN()
                          : (error_dx / norm_x)));
    }
}

static std::unique_ptr<ProcessData> makeProcessData(
    std::unique_ptr<NumLib::TimeStepAlgorithm>&& timestepper,
    NumLib::NonlinearSolverBase& nonlinear_solver,
    Process& process,
    std::unique_ptr<NumLib::TimeDiscretization>&& time_disc,
    std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit,
    ProcessOutput&& process_output)
{
    using Tag = NumLib::NonlinearSolverTag;

    if (auto* nonlinear_solver_picard =
            dynamic_cast<NumLib::NonlinearSolver<Tag::Picard>*>(
                &nonlinear_solver))
    {
        return std::make_unique<ProcessData>(
            std::move(timestepper), *nonlinear_solver_picard,
            std::move(conv_crit), std::move(time_disc), process,
            std::move(process_output));
    }
    if (auto* nonlinear_solver_newton =
            dynamic_cast<NumLib::NonlinearSolver<Tag::Newton>*>(
                &nonlinear_solver))
    {
        return std::make_unique<ProcessData>(
            std::move(timestepper), *nonlinear_solver_newton,
            std::move(conv_crit), std::move(time_disc), process,
            std::move(process_output));
    }

    OGS_FATAL("Encountered unknown nonlinear solver type. Aborting");
}

std::pair<MatVec*, bool>
SimpleMatrixVectorProvider::
get_(std::size_t& id,
     std::map<std::size_t, MatVec*>& unused_map,
     std::map<MatVec*, std::size_t>& used_map,
     Args&&... args)
{
    if (id >= _next_id) {
        OGS_FATAL("An obviously uninitialized id argument has been passed."
            " This might not be a serious error for the current implementation,"
            " but it might become one in the future."
            " Hence, I will abort now.");
    }

    if (do_search)
    {
        auto it = unused_map.find(id);
        if (it != unused_map.end()) // unused matrix/vector found
            return { ::detail::transfer(unused_map, used_map, it), false };
    }

    // not searched or not found, so create a new one
    id = _next_id++;
    auto res = used_map.emplace(
        MathLib::MatrixVectorTraits<MatVec>::newInstance(std::forward<Args>(args)...).release(),
        id);
    assert(res.second && "Emplacement failed.");
    return { res.first->first, true };
}

std::unique_ptr<TESFEMReactionAdaptor> TESFEMReactionAdaptor::newInstance(
    TESLocalAssemblerData const& data)
{
    auto const* ads = data.ap.react_sys.get();
    if (dynamic_cast<Adsorption::AdsorptionReaction const*>(ads) != nullptr)
    {
        return std::unique_ptr<TESFEMReactionAdaptor>(
            new TESFEMReactionAdaptorAdsorption(data));
    }
    else if (dynamic_cast<Adsorption::ReactionInert const*>(ads) != nullptr)
    {
        return std::unique_ptr<TESFEMReactionAdaptor>(
            new TESFEMReactionAdaptorInert(data));
    }
    else if (dynamic_cast<Adsorption::ReactionSinusoidal const*>(ads) !=
             nullptr)
    {
        return std::unique_ptr<TESFEMReactionAdaptor>(
            new TESFEMReactionAdaptorSinusoidal(data));
    }
    else if (dynamic_cast<Adsorption::ReactionCaOH2 const*>(ads) != nullptr)
    {
        return std::unique_ptr<TESFEMReactionAdaptor>(
            new TESFEMReactionAdaptorCaOH2(data));
    }

    OGS_FATAL("No suitable TESFEMReactionAdaptor found. Aborting.");
    return std::unique_ptr<TESFEMReactionAdaptor>(nullptr);
}

std::unique_ptr<TimeDiscretization> createTimeDiscretization(
    BaseLib::ConfigTree const& config)
{
    //! \ogs_file_param{prj__time_loop__processes__process__time_discretization__type}
    auto const type = config.getConfigParameter<std::string>("type");

    //! \ogs_file_param_special{prj__time_loop__processes__process__time_discretization__BackwardEuler}
    if (type == "BackwardEuler")
    {
        return std::make_unique<BackwardEuler>();
    }
    //! \ogs_file_param_special{prj__time_loop__processes__process__time_discretization__ForwardEuler}
    if (type == "ForwardEuler")
    {
        return std::make_unique<ForwardEuler>();
    }
    //! \ogs_file_param_special{prj__time_loop__processes__process__time_discretization__CrankNicolson}
    if (type == "CrankNicolson")
    {
        //! \ogs_file_param{prj__time_loop__processes__process__time_discretization__CrankNicolson__theta}
        auto const theta = config.getConfigParameter<double>("theta");
        return std::make_unique<CrankNicolson>(theta);
    }
    //! \ogs_file_param_special{prj__time_loop__processes__process__time_discretization__BackwardDifferentiationFormula}
    if (type == "BackwardDifferentiationFormula")
    {
        //! \ogs_file_param{prj__time_loop__processes__process__time_discretization__BackwardDifferentiationFormula__order}
        auto const order = config.getConfigParameter<unsigned>("order");
        return std::make_unique<BackwardDifferentiationFormula>(order);
    }

    OGS_FATAL("Unrecognized time discretization type `%s'", type.c_str());
}

LocalToGlobalIndexMap::LocalToGlobalIndexMap(
    std::vector<MeshLib::MeshSubset>&& mesh_subsets,
    std::vector<int> const& global_component_ids,
    std::vector<int> const& variable_component_offsets,
    std::vector<MeshLib::Element*> const& elements,
    NumLib::MeshComponentMap&& mesh_component_map,
    LocalToGlobalIndexMap::ConstructorTag)
    : _mesh_subsets(std::move(mesh_subsets)),
      _mesh_component_map(std::move(mesh_component_map)),
      _variable_component_offsets(variable_component_offsets)
{
    // Each subset in the mesh_subsets represents a single component.
    if (_mesh_subsets.size() != global_component_ids.size())
        OGS_FATAL(
            "Number of mesh subsets is not equal to number of components. "
            "There are %d mesh subsets and %d components.",
            _mesh_subsets.size(), global_component_ids.size());

    for (int i = 0; i < static_cast<int>(global_component_ids.size()); ++i)
    {
        auto const& ms = _mesh_subsets[i];

        // For all MeshSubset in mesh_subsets and each element of that
        // MeshSubset save a line of global indices.
        std::size_t const mesh_id = ms.getMeshID();

        findGlobalIndices(elements.cbegin(), elements.cend(), ms.getNodes(),
                          mesh_id, global_component_ids[i], i);
    }
}