Merge pull request ufz#2145 from TomFischer/ConstraintBCs

First version of constraint boundary conditions

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.

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.

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.

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.

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.

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.

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

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,

ProcessLib/BoundaryCondition/BoundaryCondition.h

@@ -27,6 +27,7 @@ namespace ProcessLib
 {
 struct BoundaryConditionConfig;
 struct ParameterBase;
+class Process;
 
 class BoundaryCondition
 {
@@ -65,11 +66,11 @@ class BoundaryConditionBuilder
     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 @@ class BoundaryConditionBuilder
         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,

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));

ProcessLib/BoundaryCondition/BoundaryConditionCollection.h

@@ -42,7 +42,7 @@ class BoundaryConditionCollection final
         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)
     {