diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/c_ConstraintDirichletBoundaryCondition.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/c_ConstraintDirichletBoundaryCondition.md
new file mode 100644
index 0000000000000000000000000000000000000000..02834b84b65e939b0cc6ef27e2e8fdf2174e9f80
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/c_ConstraintDirichletBoundaryCondition.md
@@ -0,0 +1,2 @@
+For the constraint Dirichlet-type boundary condition the type has to be
+ConstraintDirichlet.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_constraining_process_variable.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_constraining_process_variable.md
new file mode 100644
index 0000000000000000000000000000000000000000..e333bfb00a055bae38529eaa85f8f2f33d62e903
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_constraining_process_variable.md
@@ -0,0 +1,2 @@
+This tag specifies the process variable whose flux values constrain the current
+process variable.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_constraint_direction.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_constraint_direction.md
new file mode 100644
index 0000000000000000000000000000000000000000..f1b350eb2b77253eed7e129d131281e75faba924
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_constraint_direction.md
@@ -0,0 +1,6 @@
+The value of the tag determines the evaluation of the constraint condition.
+Possible values for the tag 'constraint_directions' are 'greater' or 'lower'.
+
+If the value 'greater' is given the condition 'calculated_flux_value >
+constraint_threshold' is evaluated. If the value 'less' is given the condition
+'calculated_flux_value < constraint_threshold' is evaluated.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_constraint_threshold.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_constraint_threshold.md
new file mode 100644
index 0000000000000000000000000000000000000000..3f6784b236bd9a2be11be553ad51e2d1353e539c
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_constraint_threshold.md
@@ -0,0 +1 @@
+Threshold value used in the constraint condition.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_constraint_type.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_constraint_type.md
new file mode 100644
index 0000000000000000000000000000000000000000..11c0e76288a4a1eca33a93c2c070b5a704b6f2b8
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_constraint_type.md
@@ -0,0 +1,3 @@
+The constraint type has to be 'Flux', i.e., the constraint is based on the
+secondary variable. It is planned to add constraints based on the value of the
+primary variables.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_parameter.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_parameter.md
new file mode 100644
index 0000000000000000000000000000000000000000..f36803a9dcae4c1602b4fac0fb36f961aaf361c6
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/ConstraintDirichletBoundaryCondition/t_parameter.md
@@ -0,0 +1 @@
+The name of the parameter that defines the Dirichlet-type condition values.
diff --git a/MeshLib/Elements/Utils.h b/MeshLib/Elements/Utils.h
index 70e1e20eded12046e0e66d0d2fb6f61a054ba202..df3fa27524dd29cab7fa4895307ae6dd0a773420 100644
--- a/MeshLib/Elements/Utils.h
+++ b/MeshLib/Elements/Utils.h
@@ -16,6 +16,7 @@
#include "MeshLib/Node.h"
#include "Element.h"
+#include "FaceRule.h"
namespace MeshLib
{
@@ -39,4 +40,34 @@ inline std::vector<Node*> getBaseNodes(std::vector<Element*> const& elements)
return base_nodes;
}
+inline MathLib::Vector3 calculateNormalizedSurfaceNormal(
+ MeshLib::Element const& surface_element,
+ MeshLib::Element const& bulk_element)
+{
+ MathLib::Vector3 surface_element_normal;
+ if (surface_element.getDimension() < 2)
+ {
+ auto const bulk_element_normal = MeshLib::FaceRule::getSurfaceNormal(
+ &bulk_element);
+ 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;
+
+ return surface_element_normal;
+}
+
} // namespace MeshLib
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index ba1385d9ef2619cb1b5500c2ca58852536de2edf..389f371d8c499b9d8e98a8f575e9fd160f66d158 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -10,11 +10,13 @@
#include "BoundaryCondition.h"
#include "BoundaryConditionConfig.h"
#include "DirichletBoundaryCondition.h"
+#include "ConstraintDirichletBoundaryCondition.h"
#include "NeumannBoundaryCondition.h"
#include "NonuniformDirichletBoundaryCondition.h"
#include "NonuniformNeumannBoundaryCondition.h"
#include "NormalTractionBoundaryCondition.h"
#include "PhaseFieldIrreversibleDamageOracleBoundaryCondition.h"
+#include "ProcessLib/Process.h"
#include "RobinBoundaryCondition.h"
namespace ProcessLib
@@ -25,11 +27,18 @@ BoundaryConditionBuilder::createBoundaryCondition(
const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
const int variable_id, const unsigned integration_order,
const unsigned shapefunction_order,
- const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters)
+ const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters,
+ Process const& process)
{
//! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
auto const type = config.config.peekConfigParameter<std::string>("type");
+ if (type == "ConstraintDirichlet")
+ {
+ return createConstraintDirichletBoundaryCondition(
+ config, dof_table, mesh, variable_id, integration_order,
+ parameters, process);
+ }
if (type == "Dirichlet")
{
return createDirichletBoundaryCondition(
@@ -90,6 +99,19 @@ BoundaryConditionBuilder::createDirichletBoundaryCondition(
*config.component_id, parameters);
}
+std::unique_ptr<BoundaryCondition>
+BoundaryConditionBuilder::createConstraintDirichletBoundaryCondition(
+ const BoundaryConditionConfig& config,
+ const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
+ const int variable_id, const unsigned integration_order,
+ const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters,
+ Process const& process)
+{
+ return ProcessLib::createConstraintDirichletBoundaryCondition(
+ config.config, config.mesh, dof_table, variable_id, integration_order,
+ *config.component_id, parameters, process);
+}
+
std::unique_ptr<BoundaryCondition>
BoundaryConditionBuilder::createNeumannBoundaryCondition(
const BoundaryConditionConfig& config,
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.h b/ProcessLib/BoundaryCondition/BoundaryCondition.h
index ef2638956fccab40080a66fb3c084af439628735..43fc1db7f3816bec3e1ac2ab0f863e7cc91d7aa3 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.h
@@ -27,6 +27,7 @@ namespace ProcessLib
{
struct BoundaryConditionConfig;
struct ParameterBase;
+class Process;
class BoundaryCondition
{
@@ -65,11 +66,11 @@ public:
virtual std::unique_ptr<BoundaryCondition> createBoundaryCondition(
const BoundaryConditionConfig& config,
const NumLib::LocalToGlobalIndexMap& dof_table,
- const MeshLib::Mesh& mesh,
- const int variable_id, const unsigned integration_order,
- const unsigned shapefunction_order,
+ const MeshLib::Mesh& mesh, const int variable_id,
+ const unsigned integration_order, const unsigned shapefunction_order,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
- parameters);
+ parameters,
+ Process const& process);
protected:
virtual std::unique_ptr<BoundaryCondition> createDirichletBoundaryCondition(
@@ -81,6 +82,16 @@ protected:
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
parameters);
+ virtual std::unique_ptr<BoundaryCondition>
+ createConstraintDirichletBoundaryCondition(
+ const BoundaryConditionConfig& config,
+ const NumLib::LocalToGlobalIndexMap& dof_table,
+ const MeshLib::Mesh& mesh, const int variable_id,
+ const unsigned integration_order,
+ const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
+ parameters,
+ Process const& process);
+
virtual std::unique_ptr<BoundaryCondition> createNeumannBoundaryCondition(
const BoundaryConditionConfig& config,
const NumLib::LocalToGlobalIndexMap& dof_table,
diff --git a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp
index d2c6d1783f8fa0b4b763412f06da72eeb0f7b3f4..63293a99ab7efe5da2acb545b45f7fd6c6549afc 100644
--- a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp
@@ -24,15 +24,15 @@ void BoundaryConditionCollection::addBCsForProcessVariables(
std::vector<std::reference_wrapper<ProcessVariable>> const&
process_variables,
NumLib::LocalToGlobalIndexMap const& dof_table,
- unsigned const integration_order)
+ unsigned const integration_order, Process const& process)
{
for (int variable_id = 0;
variable_id < static_cast<int>(process_variables.size());
++variable_id)
{
ProcessVariable& pv = process_variables[variable_id];
- auto bcs = pv.createBoundaryConditions(dof_table, variable_id,
- integration_order, _parameters);
+ auto bcs = pv.createBoundaryConditions(
+ dof_table, variable_id, integration_order, _parameters, process);
std::move(bcs.begin(), bcs.end(),
std::back_inserter(_boundary_conditions));
diff --git a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
index 679f76579600208e9c06ea9d56425fb3bfa7ce2d..b511506028c9900542f20d89cc72ae225fb67468 100644
--- a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
+++ b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
@@ -42,7 +42,7 @@ public:
std::vector<std::reference_wrapper<ProcessVariable>> const&
process_variables,
NumLib::LocalToGlobalIndexMap const& dof_table,
- unsigned const integration_order);
+ unsigned const integration_order, Process const& process);
void preTimestep(const double t, GlobalVector const& x)
{
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..acce3726620c45b75644b929f8c681063e7d7717
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
@@ -0,0 +1,325 @@
+/**
+ * \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_mesh, _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,
+ [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 =
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__ConstraintDirichletBoundaryCondition__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());
+ }
+
+ // Todo (TF) Open question: How to specify which getFlux function should be
+ // used for the constraint calculation?
+ auto const constraining_process_variable =
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__ConstraintDirichletBoundaryCondition__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 =
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__ConstraintDirichletBoundaryCondition__constraint_threshold}
+ config.getConfigParameter<double>("constraint_threshold");
+
+ 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 != "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__ConstraintDirichletBoundaryCondition__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..f87b56b25de86e6334b20e4ba72d81b235c2d7e4
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
@@ -0,0 +1,173 @@
+/**
+ * \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/Elements.h"
+#include "MeshLib/Elements/Utils.h"
+
+namespace ProcessLib
+{
+
+struct IntegrationPointData final
+{
+ IntegrationPointData(double const& detJ,
+ double const& integral_measure,
+ double const& integration_weight,
+ MathLib::Point3d&& bulk_element_point_)
+ : detJ_times_integralMeasure_times_weight(detJ * integral_measure *
+ integration_weight),
+ bulk_element_point(bulk_element_point_)
+ {
+ }
+
+ double const detJ_times_integralMeasure_times_weight;
+ MathLib::Point3d bulk_element_point;
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+class ConstraintDirichletBoundaryConditionLocalAssemblerInterface
+{
+public:
+ virtual ~ConstraintDirichletBoundaryConditionLocalAssemblerInterface() =
+ default;
+
+ virtual double integrate(
+ GlobalVector const& x, double const t,
+ 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_mesh The bulk mesh the process is defined on.
+ /// @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, MeshLib::Mesh const& bulk_mesh,
+ 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),
+ _surface_element_normal(MeshLib::calculateNormalizedSurfaceNormal(
+ _surface_element, *(bulk_mesh.getElements()[_bulk_element_id])))
+ {
+ (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();
+
+ auto const bulk_face_id = bulk_ids[_surface_element.getID()].second;
+ std::vector<
+ typename ShapeMatricesType::ShapeMatrices,
+ Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
+ shape_matrices;
+ shape_matrices.reserve(n_integration_points);
+ _ip_data.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);
+ auto bulk_element_point = MeshLib::getBulkElementPoint(
+ bulk_mesh, _bulk_element_id, bulk_face_id, wp);
+ _ip_data.emplace_back(shape_matrices[ip].detJ,
+ shape_matrices[ip].integralMeasure,
+ wp.getWeight(),
+ std::move(bulk_element_point));
+ }
+ }
+
+ /// 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 getFlux The function of the constraining process used to
+ /// calculate the flux.
+ double integrate(
+ GlobalVector const& x, double const t,
+ std::function<Eigen::Vector3d(
+ std::size_t const, MathLib::Point3d const&, double const,
+ GlobalVector const&)> const& getFlux) override
+ {
+ 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 bulk_flux = getFlux(
+ _bulk_element_id, _ip_data[ip].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 *
+ _ip_data[ip].detJ_times_integralMeasure_times_weight;
+ }
+ return integrated_value;
+ }
+
+private:
+ MeshLib::Element const& _surface_element;
+
+ std::vector<IntegrationPointData> _ip_data;
+
+ IntegrationMethod const _integration_method;
+ std::size_t const _bulk_element_id;
+ MathLib::Vector3 const _surface_element_normal;
+};
+
+} // ProcessLib
diff --git a/ProcessLib/HT/Tests.cmake b/ProcessLib/HT/Tests.cmake
index 872fa49e43098cacf4b8a552de2f3c4db72925cd..4036ed06b3321decafba92b098eb40121d8cfe5a 100644
--- a/ProcessLib/HT/Tests.cmake
+++ b/ProcessLib/HT/Tests.cmake
@@ -349,6 +349,18 @@ AddTest(
CoupledPressureParabolicTemperatureParabolic_ts_10_expected.vtu CoupledPressureParabolicTemperatureParabolicStaggered_pcs_1_ts_10_t_1.000000.vtu darcy_velocity darcy_velocity 1e-10 1e-10
)
+AddTest(
+ NAME HT_SimpleSynthetics_constraint_dirichlet_bc
+ PATH Parabolic/HT/SimpleSynthetics
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS constraint_bc_1e3.prj
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ DIFF_DATA
+ GLOB ConstViscosityThermalConvection_pcs_0_ts_*.vtu T T 1e-14 1e-13
+ GLOB ConstViscosityThermalConvection_pcs_0_ts_*.vtu p p 1e-15 1e-14
+)
+
# MPI/PETSc tests
AddTest(
NAME Parallel_LARGE_2D_ThermalConvection_constviscosity
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 1dfeac7f53628f8c6390a76e9e89551ea1d3a4b8..7debb42fcfbdc5fea92bdcaad9f7da432aa399b3 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -38,11 +38,11 @@ Process::Process(
_coupled_solutions(nullptr),
_integration_order(integration_order),
_process_variables(std::move(process_variables)),
- _boundary_conditions([&](const std::size_t number_of_processes)
+ _boundary_conditions([&](const std::size_t number_of_process_variables)
-> std::vector<BoundaryConditionCollection> {
std::vector<BoundaryConditionCollection> pcs_BCs;
- pcs_BCs.reserve(number_of_processes);
- for (std::size_t i = 0; i < number_of_processes; i++)
+ pcs_BCs.reserve(number_of_process_variables);
+ for (std::size_t i = 0; i < number_of_process_variables; i++)
{
pcs_BCs.emplace_back(BoundaryConditionCollection(parameters));
}
@@ -68,7 +68,7 @@ void Process::initializeProcessBoundaryConditionsAndSourceTerms(
auto& per_process_BCs = _boundary_conditions[process_id];
per_process_BCs.addBCsForProcessVariables(per_process_variables, dof_table,
- _integration_order);
+ _integration_order, *this);
auto& per_process_sts = _source_term_collections[process_id];
per_process_sts.addSourceTermsForProcessVariables(
diff --git a/ProcessLib/ProcessVariable.cpp b/ProcessLib/ProcessVariable.cpp
index fcb7c71bb62b6071ac858e8acd24b0110bdaf6e5..87a84d7d38b89b430b700bdb17637cd9c80f7e71 100644
--- a/ProcessLib/ProcessVariable.cpp
+++ b/ProcessLib/ProcessVariable.cpp
@@ -173,7 +173,8 @@ ProcessVariable::createBoundaryConditions(
const NumLib::LocalToGlobalIndexMap& dof_table,
const int variable_id,
unsigned const integration_order,
- std::vector<std::unique_ptr<ParameterBase>> const& parameters)
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ Process const& process)
{
std::vector<std::unique_ptr<BoundaryCondition>> bcs;
bcs.reserve(_bc_configs.size());
@@ -182,7 +183,7 @@ ProcessVariable::createBoundaryConditions(
{
auto bc = _bc_builder->createBoundaryCondition(
config, dof_table, _mesh, variable_id, integration_order,
- _shapefunction_order, parameters);
+ _shapefunction_order, parameters, process);
bcs.push_back(std::move(bc));
}
diff --git a/ProcessLib/ProcessVariable.h b/ProcessLib/ProcessVariable.h
index 11343c425d81274966bab713b81b6c3c7adfaa75..9cd33a18caa6ab1fe2a219235012f80f54518064 100644
--- a/ProcessLib/ProcessVariable.h
+++ b/ProcessLib/ProcessVariable.h
@@ -26,7 +26,7 @@ namespace ProcessLib
class BoundaryConditionBuilder;
/// A named process variable. Its properties includes the mesh, and the initial
-/// and boundary conditions.
+/// and boundary conditions as well as the source terms.
class ProcessVariable
{
public:
@@ -53,7 +53,8 @@ public:
std::vector<std::unique_ptr<BoundaryCondition>> createBoundaryConditions(
const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id,
unsigned const integration_order,
- std::vector<std::unique_ptr<ParameterBase>> const& parameters);
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ Process const& process);
std::vector<std::unique_ptr<NodalSourceTerm>> createSourceTerms(
const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id,
diff --git a/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_0_t_0.000000.vtu b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_0_t_0.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..8459471b96979f888c475ed60f4a27db6f49261c
--- /dev/null
+++ b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_0_t_0.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a71f2792610cb5301f403cb9db81ff7ffb292a07fa9838a4bfa9280cb974944f
+size 38479
diff --git a/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_1_t_0.000000.vtu b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_1_t_0.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..3f91edb6d91a1cbedad285a5e0d68afc0e6765cf
--- /dev/null
+++ b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_1_t_0.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f877b776ec381e0bb2565f21da42aa78c3df860a9ec2ef67d78de7c06c4c7a83
+size 53644
diff --git a/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_2_t_0.000010.vtu b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_2_t_0.000010.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..10918fcdaa2ec3cc8b536adc0ae0af766d985b1b
--- /dev/null
+++ b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_2_t_0.000010.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:15236a721f365aef10fae0fc2d8002cee0aa66d80f0713ae780b16bf17551d81
+size 54049
diff --git a/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_3_t_0.001010.vtu b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_3_t_0.001010.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..0c4003561a97bb497a91e2adf250ed6baa0fafb7
--- /dev/null
+++ b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_3_t_0.001010.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6227f1e0703c4468020492ab37d3d49444264d15eb688837240ba3f495ebbbe1
+size 53972
diff --git a/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_4_t_0.101010.vtu b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_4_t_0.101010.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..4532a88fa7f8e4c0597d41fd6a29c8d95fb17d47
--- /dev/null
+++ b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_4_t_0.101010.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:888c4a960d53d26817702c30e4f241be67c81b44c57fad54234ec0d362317590
+size 54273
diff --git a/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_5_t_1.101010.vtu b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_5_t_1.101010.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..c4e6f9a697e239fb20ab70148afdf12210bc9890
--- /dev/null
+++ b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_5_t_1.101010.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2b1ed670e29384ded16d9fdb55d6fe44fac2a710f2e5a5bd28a25b84d01e8039
+size 53867
diff --git a/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_6_t_10.000000.vtu b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_6_t_10.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..513e31a8397d38dff9f5d854cecaf569647930e9
--- /dev/null
+++ b/Tests/Data/Parabolic/HT/SimpleSynthetics/ConstViscosityThermalConvection_pcs_0_ts_6_t_10.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:975dec6ad7bcaa8bcd2a8c8669fcc35b5c44172638738b4d6e0945a74f0edbd0
+size 53011
diff --git a/Tests/Data/Parabolic/HT/SimpleSynthetics/constraint_bc_1e3.prj b/Tests/Data/Parabolic/HT/SimpleSynthetics/constraint_bc_1e3.prj
new file mode 100644
index 0000000000000000000000000000000000000000..e76ad2b33bbaddaf405b03e57a37414e3c5a2895
--- /dev/null
+++ b/Tests/Data/Parabolic/HT/SimpleSynthetics/constraint_bc_1e3.prj
@@ -0,0 +1,405 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>cube_1x1x1_hex_1e3.vtu</mesh>
+ <geometry>cube_1x1x1.gml</geometry>
+ <processes>
+ <process>
+ <name>SimpleThermalConvection</name>
+ <type>HT</type>
+ <integration_order>2</integration_order>
+ <process_variables>
+ <temperature>T</temperature>
+ <pressure>p</pressure>
+ </process_variables>
+ <fluid>
+ <density>
+ <type>TemperatureDependent</type>
+ <rho0>1000</rho0>
+ <temperature0>20</temperature0>
+ <beta>4.3e-4</beta>
+ </density>
+ <viscosity>
+ <type>Constant</type>
+ <value>1.0e-3</value>
+ </viscosity>
+ <specific_heat_capacity>
+ <type>Constant</type>
+ <value>4200</value>
+ </specific_heat_capacity>
+ </fluid>
+ <porous_medium>
+ <porous_medium id="0">
+ <permeability>
+ <permeability_tensor_entries>kappa1</permeability_tensor_entries>
+ <type>Constant</type>
+ </permeability>
+ <porosity>
+ <type>Constant</type>
+ <porosity_parameter>constant_porosity_parameter</porosity_parameter>
+ </porosity>
+ <storage>
+ <type>Constant</type>
+ <value>0.0</value>
+ </storage>
+ </porous_medium>
+ </porous_medium>
+ <density_solid>rho_solid</density_solid>
+ <specific_heat_capacity_solid>c_solid</specific_heat_capacity_solid>
+ <thermal_conductivity_solid>lambda_solid</thermal_conductivity_solid>
+ <thermal_conductivity_fluid>lambda_fluid</thermal_conductivity_fluid>
+ <thermal_dispersivity>
+ <longitudinal>alpha_l</longitudinal>
+ <transversal>alpha_t</transversal>
+ </thermal_dispersivity>
+ <specific_body_force>0 0 0</specific_body_force>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="darcy_velocity" output_name="darcy_velocity"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="SimpleThermalConvection">
+ <nonlinear_solver>basic_picard</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1.e-3</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <output>
+ <variables>
+ <variable>T</variable>
+ <variable>p</variable>
+ <variable>darcy_velocity</variable>
+ </variables>
+ </output>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial> 0.0 </t_initial>
+ <t_end> 10 </t_end>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <delta_t>1e-7</delta_t>
+ </pair>
+ <pair>
+ <repeat>1</repeat>
+ <delta_t>1e-5</delta_t>
+ </pair>
+ <pair>
+ <repeat>1</repeat>
+ <delta_t>1e-3</delta_t>
+ </pair>
+ <pair>
+ <repeat>1</repeat>
+ <delta_t>1e-1</delta_t>
+ </pair>
+ <pair>
+ <repeat>1</repeat>
+ <delta_t>1</delta_t>
+ </pair>
+ <pair>
+ <repeat>1</repeat>
+ <delta_t>10</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>ConstViscosityThermalConvection</prefix>
+ <timesteps>
+ <pair>
+ <repeat> 10 </repeat>
+ <each_steps> 1 </each_steps>
+ </pair>
+ </timesteps>
+ </output>
+ </time_loop>
+ <parameters>
+ <parameter>
+ <name>rho_solid</name>
+ <type>Constant</type>
+ <value>0.0</value>
+ </parameter>
+ <parameter>
+ <name>rho_fluid</name>
+ <type>Constant</type>
+ <value>1000</value>
+ </parameter>
+ <parameter>
+ <name>lambda_solid</name>
+ <type>Constant</type>
+ <value>3.0</value>
+ </parameter>
+ <parameter>
+ <name>lambda_fluid</name>
+ <type>Constant</type>
+ <value>0.65</value>
+ </parameter>
+ <parameter>
+ <name>alpha_l</name>
+ <type>Constant</type>
+ <value>0.0</value>
+ </parameter>
+ <parameter>
+ <name>alpha_t</name>
+ <type>Constant</type>
+ <value>0.0</value>
+ </parameter>
+ <parameter>
+ <name>c_solid</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>T0</name>
+ <type>Constant</type>
+ <value>20</value>
+ </parameter>
+ <parameter>
+ <name>P0</name>
+ <type>Constant</type>
+ <value>2.99e7</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_left0</name>
+ <type>Constant</type>
+ <value>3e7</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_left1</name>
+ <type>Constant</type>
+ <value>2.999e7</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_left2</name>
+ <type>Constant</type>
+ <value>2.998e7</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_left3</name>
+ <type>Constant</type>
+ <value>2.997e7</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_left4</name>
+ <type>Constant</type>
+ <value>2.996e7</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_left5</name>
+ <type>Constant</type>
+ <value>2.995e7</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_left6</name>
+ <type>Constant</type>
+ <value>2.994e7</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_left7</name>
+ <type>Constant</type>
+ <value>2.993e7</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_left8</name>
+ <type>Constant</type>
+ <value>2.992e7</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_left9</name>
+ <type>Constant</type>
+ <value>2.991e7</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_right</name>
+ <type>Constant</type>
+ <value>2.99e7</value>
+ </parameter>
+ <parameter>
+ <name>t_Dirichlet_left</name>
+ <type>Constant</type>
+ <value>20</value>
+ </parameter>
+ <parameter>
+ <name>t_Dirichlet_right</name>
+ <type>Constant</type>
+ <value>10</value>
+ </parameter>
+ <parameter>
+ <name>constant_porosity_parameter</name>
+ <type>Constant</type>
+ <value>0.001</value>
+ </parameter>
+ <parameter>
+ <name>kappa1</name>
+ <type>Constant</type>
+ <values>1.e-14 0 0 0 1.e-14 0 0 0 1.e-14</values>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>T</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>T0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Dirichlet</type>
+ <parameter>t_Dirichlet_left</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Dirichlet</type>
+ <parameter>t_Dirichlet_right</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ <process_variable>
+ <name>p</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>P0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left0</geometry>
+ <type>ConstraintDirichlet</type>
+ <constraint_type>Flux</constraint_type>
+ <constraining_process_variable>p</constraining_process_variable>
+ <constraint_threshold>0.0</constraint_threshold>
+ <constraint_direction>greater</constraint_direction>
+ <parameter>p_Dirichlet_left0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left1</geometry>
+ <type>ConstraintDirichlet</type>
+ <constraint_type>Flux</constraint_type>
+ <constraining_process_variable>p</constraining_process_variable>
+ <constraint_threshold>0.0</constraint_threshold>
+ <constraint_direction>greater</constraint_direction>
+ <parameter>p_Dirichlet_left1</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left2</geometry>
+ <type>ConstraintDirichlet</type>
+ <parameter>p_Dirichlet_left2</parameter>
+ <constraint_type>Flux</constraint_type>
+ <constraining_process_variable>p</constraining_process_variable>
+ <constraint_threshold>0.0</constraint_threshold>
+ <constraint_direction>greater</constraint_direction>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left3</geometry>
+ <type>ConstraintDirichlet</type>
+ <constraint_type>Flux</constraint_type>
+ <constraining_process_variable>p</constraining_process_variable>
+ <constraint_threshold>0.0</constraint_threshold>
+ <constraint_direction>greater</constraint_direction>
+ <parameter>p_Dirichlet_left3</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left4</geometry>
+ <type>ConstraintDirichlet</type>
+ <constraint_type>Flux</constraint_type>
+ <constraining_process_variable>p</constraining_process_variable>
+ <constraint_threshold>0.0</constraint_threshold>
+ <constraint_direction>greater</constraint_direction>
+ <parameter>p_Dirichlet_left4</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left5</geometry>
+ <type>ConstraintDirichlet</type>
+ <constraint_type>Flux</constraint_type>
+ <constraining_process_variable>p</constraining_process_variable>
+ <constraint_threshold>0.0</constraint_threshold>
+ <constraint_direction>greater</constraint_direction>
+ <parameter>p_Dirichlet_left5</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left6</geometry>
+ <type>ConstraintDirichlet</type>
+ <constraint_type>Flux</constraint_type>
+ <constraining_process_variable>p</constraining_process_variable>
+ <constraint_threshold>0.0</constraint_threshold>
+ <constraint_direction>greater</constraint_direction>
+ <parameter>p_Dirichlet_left6</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left7</geometry>
+ <type>ConstraintDirichlet</type>
+ <constraint_type>Flux</constraint_type>
+ <constraining_process_variable>p</constraining_process_variable>
+ <constraint_threshold>0.0</constraint_threshold>
+ <constraint_direction>greater</constraint_direction>
+ <parameter>p_Dirichlet_left7</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left8</geometry>
+ <type>ConstraintDirichlet</type>
+ <constraint_type>Flux</constraint_type>
+ <constraining_process_variable>p</constraining_process_variable>
+ <constraint_threshold>0.0</constraint_threshold>
+ <constraint_direction>greater</constraint_direction>
+ <parameter>p_Dirichlet_left8</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left9</geometry>
+ <type>ConstraintDirichlet</type>
+ <constraint_type>Flux</constraint_type>
+ <constraining_process_variable>p</constraining_process_variable>
+ <constraint_threshold>0.0</constraint_threshold>
+ <constraint_direction>greater</constraint_direction>
+ <parameter>p_Dirichlet_left9</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Dirichlet</type>
+ <parameter>p_Dirichlet_right</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_picard</name>
+ <type>Picard</type>
+ <max_iter>1000</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i bicgstab -p jacobi -tol 1e-12 -maxiter 10000</lis>
+ <petsc>
+ <prefix>T</prefix>
+ <parameters> -T_ksp_type bcgs -T_pc_type bjacobi -T_ksp_rtol 1e-16 -T_ksp_max_it 4000</parameters>
+ </petsc>
+ <eigen>
+ <solver_type>BiCGSTAB</solver_type>
+ <precon_type>ILUT</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-12</error_tolerance>
+ </eigen>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Parabolic/HT/SimpleSynthetics/cube_1x1x1.gml b/Tests/Data/Parabolic/HT/SimpleSynthetics/cube_1x1x1.gml
index 46b181cb3906ac5b1fda28eae10ea426cb27ba1b..fc62efe078c5ea66fba547223d2706b5553861fb 100644
--- a/Tests/Data/Parabolic/HT/SimpleSynthetics/cube_1x1x1.gml
+++ b/Tests/Data/Parabolic/HT/SimpleSynthetics/cube_1x1x1.gml
@@ -1,3 +1,3 @@
version https://git-lfs.github.com/spec/v1
-oid sha256:b2a6b8f99867442bfc6f06c8b81cd5aa29b208d07842327cb0076ebc923e087a
-size 1690
+oid sha256:32af7a16c6e10edd667a7634dc4b355328bb8c5448f21b71cf6d55a8001f2970
+size 3998