[PL] Impl. of ConstraintDirichletBoundaryCondition.

ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "ConstraintDirichletBoundaryCondition.h"
+
+#include <algorithm>
+#include <vector>
+#include <logog/include/logog.hpp>
+
+#include "MeshLib/Node.h"
+#include "MeshLib/MeshSearch/NodeSearch.h"  // for getUniqueNodes
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace ProcessLib
+{
+ConstraintDirichletBoundaryCondition::ConstraintDirichletBoundaryCondition(
+    Parameter<double> const& parameter,
+    NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
+    int const component_id, MeshLib::Mesh const& bc_mesh,
+    unsigned const integration_order, MeshLib::Mesh const& bulk_mesh,
+    double const constraint_threshold, bool const lower,
+    std::function<Eigen::Vector3d(std::size_t const, MathLib::Point3d const&,
+                                  double const, GlobalVector const&)>
+        getFlux)
+    : _parameter(parameter),
+      _variable_id(variable_id),
+      _component_id(component_id),
+      _bc_mesh(bc_mesh),
+      _integration_order(integration_order),
+      _constraint_threshold(constraint_threshold),
+      _lower(lower),
+      _bulk_mesh(bulk_mesh),
+      _getFlux(getFlux)
+{
+    if (variable_id >=
+            static_cast<int>(dof_table_bulk.getNumberOfVariables()) ||
+        component_id >=
+            dof_table_bulk.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_bulk.getNumberOfVariables(),
+            dof_table_bulk.getNumberOfVariableComponents(variable_id));
+    }
+
+    std::vector<MeshLib::Node*> const& bc_nodes = _bc_mesh.getNodes();
+    DBUG(
+        "Found %d nodes for constraint Dirichlet BCs for the variable %d and "
+        "component %d",
+        bc_nodes.size(), variable_id, component_id);
+
+    MeshLib::MeshSubset bc_mesh_subset{_bc_mesh, bc_nodes};
+
+    // Create local DOF table from intersected mesh subsets for the given
+    // variable and component ids.
+    _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
+        variable_id, {component_id}, std::move(bc_mesh_subset)));
+
+    auto const& bc_elements(_bc_mesh.getElements());
+    _local_assemblers.resize(bc_elements.size());
+    _flux_values.resize(bc_elements.size());
+    // create _bulk_ids vector
+    auto const* bulk_element_ids =
+        _bc_mesh.getProperties().getPropertyVector<std::size_t>(
+            "bulk_element_ids");
+    if (!bulk_element_ids)
+    {
+        OGS_FATAL(
+            "The boundary mesh '%s' doesn't contain the needed property "
+            "'bulk_element_ids'.",
+            _bc_mesh.getName().c_str());
+    }
+    auto const* bulk_node_ids =
+        _bc_mesh.getProperties().getPropertyVector<std::size_t>(
+            "bulk_node_ids");
+    if (!bulk_node_ids)
+    {
+        OGS_FATAL(
+            "The boundary mesh '%s' doesn't contain the needed property "
+            "'bulk_node_ids'.",
+            _bc_mesh.getName().c_str());
+    }
+    auto const& bulk_nodes = bulk_mesh.getNodes();
+
+    auto get_bulk_element_face_id =
+        [&](auto const bulk_element_id, MeshLib::Element const* bc_elem) {
+            auto const* bulk_elem = _bulk_mesh.getElement(bulk_element_id);
+            std::array<MeshLib::Node*, 3> nodes{
+                {bulk_nodes[(*bulk_node_ids)[bc_elem->getNode(0)->getID()]],
+                 bulk_nodes[(*bulk_node_ids)[bc_elem->getNode(1)->getID()]],
+                 bulk_nodes[(*bulk_node_ids)[bc_elem->getNode(2)->getID()]]}};
+            return bulk_elem->identifyFace(nodes.data());
+        };
+
+    _bulk_ids.reserve(bc_elements.size());
+    std::transform(begin(bc_elements), end(bc_elements),
+                   std::back_inserter(_bulk_ids), [&](auto const* bc_element) {
+                       auto const bulk_element_id =
+                           (*bulk_element_ids)[bc_element->getID()];
+                       return std::make_pair(bulk_element_id,
+                                             get_bulk_element_face_id(
+                                                 bulk_element_id, bc_element));
+                   });
+
+    const int shape_function_order = 1;
+
+    ProcessLib::createLocalAssemblers<
+        ConstraintDirichletBoundaryConditionLocalAssembler>(
+        _bulk_mesh.getDimension(), _bc_mesh.getElements(), *_dof_table_boundary,
+        shape_function_order, _local_assemblers, _bc_mesh.isAxiallySymmetric(),
+        _integration_order, _bulk_ids);
+}
+
+void ConstraintDirichletBoundaryCondition::preTimestep(double t,
+                                                       GlobalVector const& x)
+{
+    DBUG(
+        "ConstraintDirichletBoundaryCondition::preTimestep: computing flux "
+        "constraints");
+    for (auto const* boundary_element : _bc_mesh.getElements())
+    {
+        _flux_values[boundary_element->getID()] =
+            _local_assemblers[boundary_element->getID()]->integrate(
+                x, t, _bulk_mesh,
+                [this](std::size_t const element_id,
+                       MathLib::Point3d const& pnt, double const t,
+                       GlobalVector const& x) {
+                    return _getFlux(element_id, pnt, t, x);
+                });
+    }
+}
+
+void ConstraintDirichletBoundaryCondition::getEssentialBCValues(
+    const double t, NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
+{
+    SpatialPosition pos;
+
+    bc_values.ids.clear();
+    bc_values.values.clear();
+
+    std::vector<std::pair<GlobalIndexType, double>> tmp_bc_values;
+
+    auto isFlux = [&](const std::size_t element_id) {
+        return _lower ? _flux_values[element_id] < _constraint_threshold
+                      : _flux_values[element_id] > _constraint_threshold;
+    };
+
+    for (auto const* boundary_element : _bc_mesh.getElements())
+    {
+        // check if the boundary element is active
+        if (isFlux(boundary_element->getID()))
+        {
+            continue;
+        }
+
+        // loop over the boundary element nodes and set the dirichlet value
+        unsigned const number_nodes = boundary_element->getNumberOfNodes();
+        for (unsigned i = 0; i < number_nodes; ++i)
+        {
+            auto const id = boundary_element->getNode(i)->getID();
+            pos.setNodeID(id);
+
+            MeshLib::Location l(_bulk_mesh.getID(), MeshLib::MeshItemType::Node,
+                                id);
+            // TODO: that might be slow, but only done once
+            const auto g_idx = _dof_table_boundary->getGlobalIndex(
+                l, _variable_id, _component_id);
+            if (g_idx == NumLib::MeshComponentMap::nop)
+                continue;
+            // For the DDC approach (e.g. with PETSc option), the negative
+            // index of g_idx means that the entry by that index is a ghost one,
+            // which should be dropped. Especially for PETSc routines
+            // MatZeroRows and MatZeroRowsColumns, which are called to apply the
+            // Dirichlet BC, the negative index is not accepted like other
+            // matrix or vector PETSc routines. Therefore, the following
+            // if-condition is applied.
+            if (g_idx >= 0)
+            {
+                tmp_bc_values.emplace_back(g_idx, _parameter(t,pos).front());
+            }
+        }
+    }
+
+    if (tmp_bc_values.empty())
+    {
+        DBUG("The domain on which the boundary is defined is empty");
+        return;
+    }
+
+    // store unique pairs of node id and value with respect to the node id in
+    // the bc_values vector. The values related to the same node id are
+    // averaged.
+    // first: sort the (node id, value) pairs according to the node id
+    std::sort(tmp_bc_values.begin(), tmp_bc_values.end(),
+              [](std::pair<GlobalIndexType, double> const& a,
+                 std::pair<GlobalIndexType, double> const& b) {
+                  return a.first < b.first;
+              });
+    // second: average the values over equal node id ranges
+    unsigned cnt = 1;
+    GlobalIndexType current_id = tmp_bc_values.begin()->first;
+    double sum = tmp_bc_values.begin()->second;
+    for (auto const& tmp_bc_value : tmp_bc_values)
+    {
+        if (tmp_bc_value.first == current_id)
+        {
+            cnt++;
+            sum += tmp_bc_value.second;
+        }
+        else
+        {
+            bc_values.ids.emplace_back(current_id);
+            bc_values.values.emplace_back(sum/cnt);
+            cnt = 1;
+            sum = tmp_bc_value.second;
+            current_id = tmp_bc_value.first;
+        }
+    }
+    bc_values.ids.emplace_back(current_id);
+    bc_values.values.emplace_back(sum / cnt);
+
+    DBUG("Found %d dirichlet values.", bc_values.ids.size());
+    for (unsigned i = 0; i < bc_values.ids.size(); ++i)
+    {
+        DBUG("\tid: %d, value: %e", bc_values.ids[i], bc_values.values[i]);
+    }
+}
+
+std::unique_ptr<ConstraintDirichletBoundaryCondition>
+createConstraintDirichletBoundaryCondition(
+    BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
+    NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
+    unsigned const integration_order, int const component_id,
+    std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+    Process const& constraining_process)
+{
+    DBUG("Constructing ConstraintDirichletBoundaryCondition from config.");
+    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
+    config.checkConfigParameter("type", "ConstraintDirichlet");
+
+    auto const constraint_type =
+        config.getConfigParameter<std::string>("constraint_type");
+    if (constraint_type != "Flux")
+    {
+        OGS_FATAL("The constraint type is '%s', but has to be 'Flux'.",
+                  constraint_type.c_str());
+    }
+
+    auto const constraining_process_variable =
+        config.getConfigParameter<std::string>("constraining_process_variable");
+
+    if (!constraining_process.isMonolithicSchemeUsed())
+    {
+        OGS_FATAL(
+            "The constraint dirichlet boundary condition is implemented only "
+            "for monolithic implemented processes.");
+    }
+    const int process_id = 0;
+    auto process_variables =
+        constraining_process.getProcessVariables(process_id);
+    auto constraining_pv =
+        std::find_if(process_variables.cbegin(), process_variables.cend(),
+                     [&constraining_process_variable](ProcessVariable const& pv) {
+                         return pv.getName() == constraining_process_variable;
+                     });
+    if (constraining_pv == std::end(process_variables))
+    {
+        auto const& constraining_process_variable_name =
+            process_variables[variable_id].get().getName();
+        OGS_FATAL(
+            "<constraining_process_variable> in process variable name '%s' at "
+            "geometry 'TODO' : The constraining process variable is set as "
+            "'%s', but this is not specified in the project file.",
+            constraining_process_variable_name.c_str(),
+            constraining_process_variable.c_str());
+    }
+
+    auto const constraint_threshold =
+        config.getConfigParameter<double>("constraint_threshold");
+
+    auto const constraint_direction_string =
+        config.getConfigParameter<std::string>("constraint_direction");
+    if (constraint_direction_string != "greater" &&
+        constraint_direction_string != "lower")
+    {
+        OGS_FATAL(
+            "The constraint direction is '%s', but has to be either 'greater' "
+            "or 'lower'.",
+            constraint_direction_string.c_str());
+    }
+    bool const lower = constraint_direction_string == "lower";
+
+
+    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__parameter}
+    auto const param_name = config.getConfigParameter<std::string>("parameter");
+    DBUG("Using parameter %s", param_name.c_str());
+
+    auto& param = findParameter<double>(param_name, parameters, 1);
+
+    return std::make_unique<ConstraintDirichletBoundaryCondition>(
+        param, dof_table_bulk, variable_id, component_id, bc_mesh,
+        integration_order, constraining_process.getMesh(), constraint_threshold,
+        lower,
+        [&constraining_process](std::size_t const element_id,
+                                MathLib::Point3d const& pnt, double const t,
+                                GlobalVector const& x) {
+            return constraining_process.getFlux(element_id, pnt, t, x);
+        });
+}
+
+}  // namespace ProcessLib

ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/IndexValueVector.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "BoundaryCondition.h"
+
+#include "ConstraintDirichletBoundaryConditionLocalAssembler.h"
+
+namespace ProcessLib
+{
+/// The ConstraintDirichletBoundaryCondition class describes a Dirichlet-type
+/// boundary condition that is constant in space and time where the domain can
+/// shrink and grow within the simulation. The expected parameter in the passed
+/// configuration is "value" which, when not present defaults to zero.
+class ConstraintDirichletBoundaryCondition final : public BoundaryCondition
+{
+public:
+    /// @param parameter Used for setting the values for the boundary condition.
+    /// @param dof_table_bulk The bulk local to global index map is used to
+    /// derive the local to global index map for the boundary.
+    /// @param variable_id The variable id is needed to determine the global
+    /// index.
+    /// @param component_id The component id is needed to determine the global
+    /// index.
+    /// @param bc_mesh Lower dimensional mesh the boundary condition is defined
+    /// on. The bc_mesh must have the two PropertyVector objects
+    /// 'bulk_element_ids' and 'bulk_node_ids' containing the corresponding
+    /// information.
+    /// @param integration_order Order the order of integration used to compute
+    /// the constraint.
+    /// @param bulk_mesh The FE mesh for the simulation.
+    /// @param constraint_threshold The threshold value used for the switch
+    /// off/on decision.
+    /// @param lower Boolean value used for the calculation of the constraint
+    /// criterion, i.e., if lower is set to true the criterion 'calculated_value
+    /// < constraint_threshold' is evaluated to switch on/off the boundary
+    /// condition, else 'calculated_value > constraint_threshold' is evaluated.
+    /// @param getFlux The function used for the flux calculation.
+    /// @note The function has to be stored by value, else the process value is
+    /// not captured properly.
+    ConstraintDirichletBoundaryCondition(
+        Parameter<double> const& parameter,
+        NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
+        int const variable_id, int const component_id,
+        MeshLib::Mesh const& bc_mesh, unsigned const integration_order,
+        MeshLib::Mesh const& bulk_mesh, double const constraint_threshold,
+        bool const lower,
+        std::function<Eigen::Vector3d(std::size_t const,
+                                      MathLib::Point3d const&, double const,
+                                      GlobalVector const&)>
+            getFlux);
+
+    void preTimestep(double t, GlobalVector const& x) override;
+
+    void getEssentialBCValues(
+        const double t,
+        NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
+
+private:
+    Parameter<double> const& _parameter;
+
+    /// Local dof table, a subset of the global one restricted to the
+    /// participating number of elements of the boundary condition.
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table_boundary;
+
+    int const _variable_id;
+    int const _component_id;
+
+    /// Vector of (lower-dimensional) boundary elements on which the boundary
+    /// condition is defined.
+    MeshLib::Mesh const& _bc_mesh;
+
+    /// Integration order for integration over the lower-dimensional elements
+    unsigned const _integration_order;
+
+    /// The first item of the pair is the element id in the bulk mesh, the
+    /// second item is the face id of the bulk element that is part of the
+    /// boundary
+    std::vector<std::pair<std::size_t, unsigned>> _bulk_ids;
+
+    /// Stores the results of the flux computations per boundary element.
+    std::vector<double> _flux_values;
+
+    /// Local assemblers for each boundary element.
+    std::vector<std::unique_ptr<
+        ConstraintDirichletBoundaryConditionLocalAssemblerInterface>>
+        _local_assemblers;
+
+    /// The threshold value used to the switch off/on the Dirichlet-type
+    /// boundary condition.
+    double const _constraint_threshold;
+
+    /// The boolean value lower is used for the calculation of the constraint
+    /// criterion, i.e., if lower is set to true the criterion 'calculated_value
+    /// < constraint_threshold' is evaluated to switch on/off the boundary
+    /// condition, else 'calculated_value > constraint_threshold' is evaluated.
+    bool const _lower;
+
+    /// The mesh _bulk_mesh is the discretized domain the process(es) are
+    /// defined on. It is needed to get values for the constraint calculation.
+    MeshLib::Mesh const& _bulk_mesh;
+
+    /// The function _getFlux calculates the flux through the boundary element.
+    std::function<Eigen::Vector3d(
+            std::size_t const, MathLib::Point3d const&, double const,
+            GlobalVector const&)> _getFlux;
+};
+
+/// The function parses the config tree and creates a
+/// ConstraintDirichletBoundaryCondition.
+std::unique_ptr<ConstraintDirichletBoundaryCondition>
+createConstraintDirichletBoundaryCondition(
+    BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
+    NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
+    unsigned const integration_order, int const component_id,
+    std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+    Process const& constraining_process);
+}  // namespace ProcessLib

ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+
+#include "NumLib/DOF/DOFTableUtil.h"
+
+#include "ProcessLib/Process.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "MeshLib/Elements/MapBulkElementPoint.h"
+#include "MeshLib/Elements/FaceRule.h"
+#include "MeshLib/Elements/Elements.h"
+
+namespace ProcessLib
+{
+class ConstraintDirichletBoundaryConditionLocalAssemblerInterface
+{
+public:
+    virtual ~ConstraintDirichletBoundaryConditionLocalAssemblerInterface() =
+        default;
+
+    virtual double integrate(
+        GlobalVector const& x, double const t, MeshLib::Mesh const& bulk_mesh,
+        std::function<Eigen::Vector3d(std::size_t const,
+                                      MathLib::Point3d const&, double const,
+                                      GlobalVector const&)> const& getFlux) = 0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+class ConstraintDirichletBoundaryConditionLocalAssembler final
+    : public ConstraintDirichletBoundaryConditionLocalAssemblerInterface
+{
+protected:
+    using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+    using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+    using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+
+public:
+    /// Precomputes the shape matrices for a given surface element.
+    /// @param surface_element The surface element used for precomputing the
+    /// @param local_matrix_size Unused parameter, only necessary for the
+    /// creation scheme of local assemblers
+    /// shape matrices used later on for the integration.
+    /// @param is_axially_symmetric Corrects integration measure for cylinder
+    /// coordinates.
+    /// @param integration_order The order of the integration.
+    /// @param bulk_ids Pairs of bulk element ids and bulk element face ids.
+    ConstraintDirichletBoundaryConditionLocalAssembler(
+        MeshLib::Element const& surface_element,
+        std::size_t const local_matrix_size,
+        bool const is_axially_symmetric, unsigned const integration_order,
+        std::vector<std::pair<std::size_t, unsigned>> bulk_ids)
+        : _surface_element(surface_element),
+          _integration_method(integration_order),
+          _bulk_element_id(bulk_ids[_surface_element.getID()].first),
+          _bulk_face_id(bulk_ids[_surface_element.getID()].second)
+    {
+        (void)local_matrix_size; // unused, but needed for the interface
+        using FemType =
+            NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
+
+        FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
+            &_surface_element));
+
+        std::size_t const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        std::vector<
+            typename ShapeMatricesType::ShapeMatrices,
+            Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
+            shape_matrices;
+        shape_matrices.reserve(n_integration_points);
+        _detJ_times_integralMeasure_times_weight.reserve(n_integration_points);
+        for (std::size_t ip = 0; ip < n_integration_points; ++ip)
+        {
+            shape_matrices.emplace_back(ShapeFunction::DIM, GlobalDim,
+                                        ShapeFunction::NPOINTS);
+            fe.template computeShapeFunctions<NumLib::ShapeMatrixType::N_J>(
+                _integration_method.getWeightedPoint(ip).getCoords(),
+                shape_matrices[ip], GlobalDim, is_axially_symmetric);
+            auto const& wp = _integration_method.getWeightedPoint(ip);
+            _detJ_times_integralMeasure_times_weight.push_back(
+                shape_matrices[ip].detJ * shape_matrices[ip].integralMeasure *
+                wp.getWeight());
+        }
+    }
+
+    /// Integration for the element with the id \c element_id.
+    /// @param x The global vector containing the values for numerical
+    /// integration.
+    /// @param t The point in time the the integration will be performed.
+    /// @param bulk_mesh The bulk mesh the process is defined on.
+    /// @param getFlux The function of the constraining process used to
+    /// calculate the flux.
+    double integrate(
+        GlobalVector const& x, double const t, MeshLib::Mesh const& bulk_mesh,
+        std::function<Eigen::Vector3d(
+            std::size_t const, MathLib::Point3d const&, double const,
+            GlobalVector const&)> const& getFlux) override
+    {
+        MathLib::Vector3 surface_element_normal;
+        if (_surface_element.getDimension() < 2)
+        {
+            auto const bulk_element_normal =
+                MeshLib::FaceRule::getSurfaceNormal(
+                    bulk_mesh.getElement(_bulk_element_id));
+            MathLib::Vector3 const edge_vector(*_surface_element.getNode(0),
+                                               *_surface_element.getNode(1));
+            surface_element_normal =
+                MathLib::crossProduct(bulk_element_normal, edge_vector);
+        } else {
+            surface_element_normal =
+                MeshLib::FaceRule::getSurfaceNormal(&_surface_element);
+        }
+
+        surface_element_normal.normalize();
+        // At the moment (2018-04-26) the surface normal is not oriented
+        // according to the right hand rule
+        // for correct results it is necessary to multiply the normal with -1
+        surface_element_normal *= -1;
+
+        auto const n_integration_points =
+            _integration_method.getNumberOfPoints();
+        // integrated_value +=
+        //   int_{\Gamma_e} \alpha \cdot flux \cdot normal \d \Gamma
+        double integrated_value = 0;
+        for (auto ip(0); ip < n_integration_points; ip++)
+        {
+            auto const& wp = _integration_method.getWeightedPoint(ip);
+            auto const bulk_element_point = MeshLib::getBulkElementPoint(
+                bulk_mesh, _bulk_element_id, _bulk_face_id, wp);
+            auto const bulk_flux =
+                getFlux(_bulk_element_id, bulk_element_point, t, x);
+
+            // TODO find solution for 2d case
+            double const bulk_grad_times_normal(
+                Eigen::Map<Eigen::RowVectorXd const>(bulk_flux.data(),
+                                                     bulk_flux.size())
+                    .dot(Eigen::Map<Eigen::RowVectorXd const>(
+                        surface_element_normal.getCoords(), 3)));
+
+            integrated_value += bulk_grad_times_normal *
+                                _detJ_times_integralMeasure_times_weight[ip];
+        }
+        return integrated_value;
+    }
+
+private:
+    MeshLib::Element const& _surface_element;
+
+    std::vector<double> _detJ_times_integralMeasure_times_weight;
+
+    IntegrationMethod const _integration_method;
+    std::size_t const _bulk_element_id;
+    unsigned const _bulk_face_id;
+};
+
+}  // ProcessLib