diff --git a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
index 650c53e4a09f2a89c53557a8615620cd567cf9f5..0766abc63371ebf994dce48e57230b3fb96f9d18 100644
--- a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
@@ -19,6 +19,7 @@
#include "NeumannBoundaryCondition.h"
#include "NormalTractionBoundaryCondition.h"
#include "PhaseFieldIrreversibleDamageOracleBoundaryCondition.h"
+#include "PrimaryVariableConstraintDirichletBoundaryCondition.h"
#include "RobinBoundaryCondition.h"
#include "VariableDependentNeumannBoundaryCondition.h"
@@ -107,6 +108,12 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
config.config, config.boundary_mesh, dof_table, variable_id,
integration_order, *config.component_id, parameters, process);
}
+ if (type == "PrimaryVariableConstraintDirichlet")
+ {
+ return createPrimaryVariableConstraintDirichletBoundaryCondition(
+ config.config, config.boundary_mesh, dof_table, variable_id,
+ *config.component_id, parameters);
+ }
if (type == "HCNonAdvectiveFreeComponentFlowBoundary")
{
return createHCNonAdvectiveFreeComponentFlowBoundaryCondition(
diff --git a/ProcessLib/BoundaryCondition/PrimaryVariableConstraintDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/PrimaryVariableConstraintDirichletBoundaryCondition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8a08182ebcec30439762853e13f3bc0e98e70de7
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/PrimaryVariableConstraintDirichletBoundaryCondition.cpp
@@ -0,0 +1,170 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "PrimaryVariableConstraintDirichletBoundaryCondition.h"
+
+#include <algorithm>
+#include <vector>
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/Logging.h"
+#include "DirichletBoundaryConditionAuxiliaryFunctions.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/IndexValueVector.h"
+#include "ParameterLib/Parameter.h"
+#include "ParameterLib/Utils.h"
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+PrimaryVariableConstraintDirichletBoundaryCondition::
+ PrimaryVariableConstraintDirichletBoundaryCondition(
+ ParameterLib::Parameter<double> const& parameter,
+ MeshLib::Mesh const& bc_mesh,
+ NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
+ int const variable_id, int const component_id,
+ ParameterLib::Parameter<double> const& constraint_threshold_parameter,
+ bool const less)
+ : _parameter(parameter),
+ _bc_mesh(bc_mesh),
+ _variable_id(variable_id),
+ _component_id(component_id),
+ _constraint_threshold_parameter(constraint_threshold_parameter),
+ _less(less)
+{
+ checkParametersOfDirichletBoundaryCondition(_bc_mesh, dof_table_bulk,
+ _variable_id, _component_id);
+
+ std::vector<MeshLib::Node*> const& bc_nodes = bc_mesh.getNodes();
+ MeshLib::MeshSubset bc_mesh_subset(_bc_mesh, bc_nodes);
+
+ // Create local DOF table from the BC mesh subset for the given variable
+ // and component id.
+ _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
+ variable_id, {component_id}, std::move(bc_mesh_subset)));
+}
+
+void PrimaryVariableConstraintDirichletBoundaryCondition::getEssentialBCValues(
+ const double t, GlobalVector const& x,
+ NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
+{
+ ParameterLib::SpatialPosition pos;
+
+ bc_values.ids.clear();
+ bc_values.values.clear();
+
+ auto const& nodes_in_bc_mesh = _bc_mesh.getNodes();
+ // convert mesh node ids to global index for the given component
+ bc_values.ids.reserve(bc_values.ids.size() + nodes_in_bc_mesh.size());
+ bc_values.values.reserve(bc_values.values.size() + nodes_in_bc_mesh.size());
+ for (auto const* const node : nodes_in_bc_mesh)
+ {
+ auto const id = node->getID();
+ // TODO: that might be slow, but only done once
+ auto const global_index = _dof_table_boundary->getGlobalIndex(
+ {_bc_mesh.getID(), MeshLib::MeshItemType::Node, id}, _variable_id,
+ _component_id);
+ if (global_index == NumLib::MeshComponentMap::nop)
+ {
+ continue;
+ }
+ // For the DDC approach (e.g. with PETSc option), the negative
+ // index of global_index 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 (global_index >= 0)
+ {
+ // fetch the value of the primary variable
+ auto const local_x = x.get(std::vector{global_index});
+
+ bc_values.ids.emplace_back(global_index);
+ pos.setNodeID(id);
+ bc_values.values.emplace_back(_parameter(t, pos).front());
+ }
+ }
+}
+
+std::unique_ptr<PrimaryVariableConstraintDirichletBoundaryCondition>
+createPrimaryVariableConstraintDirichletBoundaryCondition(
+ BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
+ NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
+ int const component_id,
+ const std::vector<std::unique_ptr<ParameterLib::ParameterBase>>& parameters)
+{
+ DBUG(
+ "Constructing PrimaryVariableConstraintDirichletBoundaryCondition from "
+ "config.");
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
+ config.checkConfigParameter("type", "PrimaryVariableConstraintDirichlet");
+
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__PrimaryVariableConstraintDirichlet__parameter}
+ auto const param_name = config.getConfigParameter<std::string>("parameter");
+ DBUG("Using parameter {:s}", param_name);
+
+ auto& parameter = ParameterLib::findParameter<double>(
+ param_name, parameters, 1, &bc_mesh);
+
+ auto const constraint_type =
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__ConstraintDirichletBoundaryCondition__constraint_type}
+ config.getConfigParameter<std::string>("constraint_type");
+ if (constraint_type != "PrimaryVariable")
+ {
+ OGS_FATAL(
+ "The constraint type is '{:s}', but has to be 'PrimaryVariable'.",
+ constraint_type);
+ }
+
+ auto const constraint_parameter_name = config.getConfigParameter<
+ std::string>(
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__PrimaryVariableConstraintDirichletBoundaryCondition__constraint_threshold_parameter}
+ "constraint_threshold_parameter");
+ DBUG("Using parameter {:s} as constraint_threshold_parameter",
+ constraint_parameter_name);
+
+ auto& constraint_threshold_parameter = ParameterLib::findParameter<double>(
+ constraint_parameter_name, parameters, 1, &bc_mesh);
+
+ auto const constraint_direction_string =
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__ConstraintDirichletBoundaryCondition__constraint_direction}
+ config.getConfigParameter<std::string>("constraint_direction");
+ if (constraint_direction_string != "greater" &&
+ constraint_direction_string != "less")
+ {
+ OGS_FATAL(
+ "The constraint direction is '{:s}', but has to be either "
+ "'greater' or 'less'.",
+ constraint_direction_string);
+ }
+ bool const less = constraint_direction_string == "less";
+
+// 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 (bc_mesh.getDimension() == 0 && bc_mesh.getNumberOfNodes() == 0 &&
+ bc_mesh.getNumberOfElements() == 0)
+ {
+ return nullptr;
+ }
+#endif // USE_PETSC
+
+ return std::make_unique<
+ PrimaryVariableConstraintDirichletBoundaryCondition>(
+ parameter, bc_mesh, dof_table_bulk, variable_id, component_id,
+ constraint_threshold_parameter, less);
+}
+
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/PrimaryVariableConstraintDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/PrimaryVariableConstraintDirichletBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..eb2b44aad607e6ea42fa54648c13f19d0b5b5581
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/PrimaryVariableConstraintDirichletBoundaryCondition.h
@@ -0,0 +1,75 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "BoundaryCondition.h"
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+namespace ParameterLib
+{
+template <typename T>
+struct Parameter;
+}
+
+namespace ProcessLib
+{
+// TODO docu
+/// The PrimaryVariableConstraintDirichletBoundaryCondition class describes a
+/// constant in space and time PrimaryVariableConstraintDirichlet boundary
+/// condition. The expected parameter in the passed configuration is "value"
+/// which, when not present defaults to zero.
+class PrimaryVariableConstraintDirichletBoundaryCondition final
+ : public BoundaryCondition
+{
+public:
+ PrimaryVariableConstraintDirichletBoundaryCondition(
+ ParameterLib::Parameter<double> const& parameter,
+ MeshLib::Mesh const& bc_mesh,
+ NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
+ int const variable_id, int const component_id,
+ ParameterLib::Parameter<double> const& constraint_threshold_parameter,
+ bool const less);
+
+ void getEssentialBCValues(
+ const double t, GlobalVector const& x,
+ NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
+
+private:
+ ParameterLib::Parameter<double> const& _parameter;
+
+ MeshLib::Mesh const& _bc_mesh;
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap const> _dof_table_boundary;
+ int const _variable_id;
+ int const _component_id;
+
+ /// The threshold value used to the switch off/on the Dirichlet-type
+ /// boundary condition.
+ ParameterLib::Parameter<double> const& _constraint_threshold_parameter;
+
+ /// 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 _less;
+};
+
+std::unique_ptr<PrimaryVariableConstraintDirichletBoundaryCondition>
+createPrimaryVariableConstraintDirichletBoundaryCondition(
+ BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
+ NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
+ int const component_id,
+ const std::vector<std::unique_ptr<ParameterLib::ParameterBase>>&
+ parameters);
+
+} // namespace ProcessLib