diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ec717c6c710c26651a4d1f3f08d610c551c7c50b
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
@@ -0,0 +1,319 @@
+/**
+ * \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
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..253abfceb7fddd23c46f41ad64de4d8e63e7b1f8
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h
@@ -0,0 +1,128 @@
+/**
+ * \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
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..8fa1e191c5928f43279b3019133c1be73f6dccb4
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
@@ -0,0 +1,169 @@
+/**
+ * \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