From 557f0e66d84f41bb139810fc0f1f14444406681f Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Wed, 28 Nov 2018 12:09:25 +0100
Subject: [PATCH] [BC] Added two auxiliary functions to avoid the way of
 inheritance for reducing source code duplication.

---
 .../CreateBoundaryCondition.cpp               |  3 +
 .../DirichletBoundaryCondition.cpp            | 77 ++-------------
 .../DirichletBoundaryCondition.h              | 12 +--
 ...letBoundaryConditionAuxiliaryFunctions.cpp | 95 +++++++++++++++++++
 ...chletBoundaryConditionAuxiliaryFunctions.h | 51 ++++++++++
 ...letBoundaryConditionWithinTimeInterval.cpp | 28 ++++--
 ...chletBoundaryConditionWithinTimeInterval.h | 21 +++-
 7 files changed, 202 insertions(+), 85 deletions(-)
 create mode 100644 ProcessLib/BoundaryCondition/DirichletBoundaryConditionAuxiliaryFunctions.cpp
 create mode 100644 ProcessLib/BoundaryCondition/DirichletBoundaryConditionAuxiliaryFunctions.h

diff --git a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
index 7115b756d4d..02407f68562 100644
--- a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
@@ -21,6 +21,9 @@
 #include "NormalTractionBoundaryCondition.h"
 #include "PhaseFieldIrreversibleDamageOracleBoundaryCondition.h"
 #include "RobinBoundaryCondition.h"
+
+#include "NumLib/TimeStepping/TimeInterval.h"
+
 #ifdef OGS_USE_PYTHON
 #include "Python/PythonBoundaryCondition.h"
 #endif
diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
index 8e3274407de..473ecf504ef 100644
--- a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
@@ -13,6 +13,8 @@
 #include <logog/include/logog.hpp>
 #include <vector>
 
+#include "DirichletBoundaryConditionAuxiliaryFunctions.h"
+
 #include "BaseLib/ConfigTree.h"
 #include "NumLib/IndexValueVector.h"
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
@@ -30,34 +32,10 @@ DirichletBoundaryCondition::DirichletBoundaryCondition(
       _variable_id(variable_id),
       _component_id(component_id)
 {
-    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));
-    }
-
-    if (!_bc_mesh.getProperties().existsPropertyVector<std::size_t>(
-            "bulk_node_ids"))
-    {
-        OGS_FATAL(
-            "The required bulk node ids map does not exist in the boundary "
-            "mesh '%s'.",
-            _bc_mesh.getName().c_str());
-    }
-
-    std::vector<MeshLib::Node*> const& bc_nodes = _bc_mesh.getNodes();
-    DBUG(
-        "Found %d nodes for Dirichlet BCs for the variable %d and "
-        "component "
-        "%d",
-        bc_nodes.size(), variable_id, component_id);
+    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
@@ -70,45 +48,8 @@ void DirichletBoundaryCondition::getEssentialBCValues(
     const double t, GlobalVector const& x,
     NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
 {
-    getEssentialBCValuesLocal(t, x, bc_values);
-}
-
-void DirichletBoundaryCondition::getEssentialBCValuesLocal(
-    const double t, GlobalVector const& /*x*/,
-    NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
-{
-    SpatialPosition pos;
-
-    bc_values.ids.clear();
-    bc_values.values.clear();
-
-    // convert mesh node ids to global index for the given component
-    bc_values.ids.reserve(bc_values.ids.size() + _bc_mesh.getNumberOfNodes());
-    bc_values.values.reserve(bc_values.values.size() +
-                             _bc_mesh.getNumberOfNodes());
-    for (auto const* const node : _bc_mesh.getNodes())
-    {
-        auto const id = node->getID();
-        pos.setNodeID(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)
-        {
-            bc_values.ids.emplace_back(global_index);
-            bc_values.values.emplace_back(_parameter(t, pos).front());
-        }
-    }
+    getEssentialBCValuesLocal(_parameter, _bc_mesh, *_dof_table_boundary,
+                              _variable_id, _component_id, t, x, bc_values);
 }
 
 std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(
@@ -127,8 +68,8 @@ std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(
 
     auto& param = findParameter<double>(param_name, parameters, 1);
 
-    // In case of partitioned mesh the boundary could be empty, i.e. there is no
-    // boundary condition.
+// 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.
diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h
index b913ca62694..6680f59dc6e 100644
--- a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h
@@ -18,14 +18,15 @@ class ConfigTree;
 
 namespace ProcessLib
 {
-template <typename T> struct Parameter;
+template <typename T>
+struct Parameter;
 
 // TODO docu
 /// The DirichletBoundaryCondition class describes a constant in space
 /// and time Dirichlet boundary condition.
 /// The expected parameter in the passed configuration is "value" which, when
 /// not present defaults to zero.
-class DirichletBoundaryCondition  : public BoundaryCondition
+class DirichletBoundaryCondition final : public BoundaryCondition
 {
 public:
     DirichletBoundaryCondition(
@@ -37,18 +38,13 @@ public:
         const double t, GlobalVector const& x,
         NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
 
-protected:
+private:
     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;
-
-    void getEssentialBCValuesLocal(
-        const double t, GlobalVector const& x,
-        NumLib::IndexValueVector<GlobalIndexType>& bc_values) const;
-
 };
 
 std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(
diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryConditionAuxiliaryFunctions.cpp b/ProcessLib/BoundaryCondition/DirichletBoundaryConditionAuxiliaryFunctions.cpp
new file mode 100644
index 00000000000..19ada570ff3
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/DirichletBoundaryConditionAuxiliaryFunctions.cpp
@@ -0,0 +1,95 @@
+/**
+ *
+ * \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
+ *
+ * File:   DirichletBoundaryConditionAuxiliaryFunctions.cpp
+ *
+ * Created on November 28, 2018, 11:26 AM
+ */
+#include "DirichletBoundaryConditionAuxiliaryFunctions.h"
+
+#include "NumLib/IndexValueVector.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/Parameter/Parameter.h"
+
+namespace ProcessLib
+{
+void checkParametersOfDirichletBoundaryCondition(
+    MeshLib::Mesh const& bc_mesh,
+    NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
+    int const variable_id,
+    int const component_id)
+{
+    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));
+    }
+
+    if (!bc_mesh.getProperties().existsPropertyVector<std::size_t>(
+            "bulk_node_ids"))
+    {
+        OGS_FATAL(
+            "The required bulk node ids map does not exist in the boundary "
+            "mesh '%s'.",
+            bc_mesh.getName().c_str());
+    }
+
+    DBUG(
+        "Found %d nodes for Dirichlet BCs for the variable %d and "
+        "component "
+        "%d",
+        bc_mesh.getNodes().size(), variable_id, component_id);
+}
+
+void getEssentialBCValuesLocal(
+    Parameter<double> const& parameter, MeshLib::Mesh const& bc_mesh,
+    NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
+    int const variable_id, int const component_id, const double t,
+    GlobalVector const& /*x*/,
+    NumLib::IndexValueVector<GlobalIndexType>& bc_values)
+{
+    SpatialPosition pos;
+
+    bc_values.ids.clear();
+    bc_values.values.clear();
+
+    // convert mesh node ids to global index for the given component
+    bc_values.ids.reserve(bc_values.ids.size() + bc_mesh.getNumberOfNodes());
+    bc_values.values.reserve(bc_values.values.size() +
+                             bc_mesh.getNumberOfNodes());
+    for (auto const* const node : bc_mesh.getNodes())
+    {
+        auto const id = node->getID();
+        pos.setNodeID(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)
+        {
+            bc_values.ids.emplace_back(global_index);
+            bc_values.values.emplace_back(parameter(t, pos).front());
+        }
+    }
+}
+} // end of name space
diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryConditionAuxiliaryFunctions.h b/ProcessLib/BoundaryCondition/DirichletBoundaryConditionAuxiliaryFunctions.h
new file mode 100644
index 00000000000..a3392491c30
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/DirichletBoundaryConditionAuxiliaryFunctions.h
@@ -0,0 +1,51 @@
+/**
+ *
+ * \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
+ *
+ * \brief
+ *  Defines functions that are shared by DirichletBoundaryCondition
+ *  and DirichletBoundaryConditionWithinTimeInterval, which avoid the way of
+ *  inheritance for reducing source code duplication.
+ *
+ * File:   DirichletBoundaryConditionAuxiliaryFunctions.h
+ *
+ * Created on November 28, 2018, 11:26 AM
+ */
+#pragma once
+
+#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+template <typename>
+struct IndexValueVector;
+}
+
+namespace ProcessLib
+{
+template <typename T>
+struct Parameter;
+
+void checkParametersOfDirichletBoundaryCondition(
+    MeshLib::Mesh const& bc_mesh,
+    NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
+    int const variable_id,
+    int const component_id);
+
+void getEssentialBCValuesLocal(
+    Parameter<double> const& parameter, MeshLib::Mesh const& bc_mesh,
+    NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
+    int const variable_id, int const component_id, const double t,
+    GlobalVector const& x,
+    NumLib::IndexValueVector<GlobalIndexType>& bc_values);
+}  // end of name space
diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryConditionWithinTimeInterval.cpp b/ProcessLib/BoundaryCondition/DirichletBoundaryConditionWithinTimeInterval.cpp
index 4e51b27e171..209c787e6d4 100644
--- a/ProcessLib/BoundaryCondition/DirichletBoundaryConditionWithinTimeInterval.cpp
+++ b/ProcessLib/BoundaryCondition/DirichletBoundaryConditionWithinTimeInterval.cpp
@@ -12,6 +12,9 @@
  */
 #include "DirichletBoundaryConditionWithinTimeInterval.h"
 
+#include "DirichletBoundaryCondition.h"
+#include "DirichletBoundaryConditionAuxiliaryFunctions.h"
+
 #include "BaseLib/ConfigTree.h"
 
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
@@ -29,10 +32,22 @@ DirichletBoundaryConditionWithinTimeInterval::
         Parameter<double> const& parameter, MeshLib::Mesh const& bc_mesh,
         NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
         int const variable_id, int const component_id)
-    : DirichletBoundaryCondition(parameter, bc_mesh, dof_table_bulk,
-                                 variable_id, component_id),
+    : _parameter(parameter),
+      _bc_mesh(bc_mesh),
+      _variable_id(variable_id),
+      _component_id(component_id),
       _time_interval(std::move(time_interval))
 {
+    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 DirichletBoundaryConditionWithinTimeInterval::getEssentialBCValues(
@@ -41,7 +56,8 @@ void DirichletBoundaryConditionWithinTimeInterval::getEssentialBCValues(
 {
     if (_time_interval->contains(t))
     {
-        getEssentialBCValuesLocal(t, x, bc_values);
+        getEssentialBCValuesLocal(_parameter, _bc_mesh, *_dof_table_boundary,
+                                  _variable_id, _component_id, t, x, bc_values);
         return;
     }
 
@@ -49,7 +65,7 @@ void DirichletBoundaryConditionWithinTimeInterval::getEssentialBCValues(
     bc_values.values.clear();
 }
 
-std::unique_ptr<DirichletBoundaryCondition>
+std::unique_ptr<DirichletBoundaryConditionWithinTimeInterval>
 createDirichletBoundaryConditionWithinTimeInterval(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
     NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
@@ -87,8 +103,8 @@ createDirichletBoundaryConditionWithinTimeInterval(
     config.peekConfigParameter<std::string>("time_interval");
 
     return std::make_unique<DirichletBoundaryConditionWithinTimeInterval>(
-        NumLib::createTimeInterval(config), param, bc_mesh,
-        dof_table_bulk, variable_id, component_id);
+        NumLib::createTimeInterval(config), param, bc_mesh, dof_table_bulk,
+        variable_id, component_id);
 }
 
 }  // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryConditionWithinTimeInterval.h b/ProcessLib/BoundaryCondition/DirichletBoundaryConditionWithinTimeInterval.h
index b96932fb5bd..118a4e8e124 100644
--- a/ProcessLib/BoundaryCondition/DirichletBoundaryConditionWithinTimeInterval.h
+++ b/ProcessLib/BoundaryCondition/DirichletBoundaryConditionWithinTimeInterval.h
@@ -14,7 +14,12 @@
 
 #include <memory>
 
-#include "DirichletBoundaryCondition.h"
+#include "BoundaryCondition.h"
+
+namespace BaseLib
+{
+class ConfigTree;
+}
 
 namespace NumLib
 {
@@ -23,8 +28,11 @@ class TimeInterval;
 
 namespace ProcessLib
 {
+template <typename T>
+struct Parameter;
+
 class DirichletBoundaryConditionWithinTimeInterval final
-    : public DirichletBoundaryCondition
+    : public BoundaryCondition
 {
 public:
     DirichletBoundaryConditionWithinTimeInterval(
@@ -38,10 +46,17 @@ public:
         NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
 
 private:
+    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;
+
     std::unique_ptr<NumLib::TimeInterval> _time_interval;
 };
 
-std::unique_ptr<DirichletBoundaryCondition>
+std::unique_ptr<DirichletBoundaryConditionWithinTimeInterval>
 createDirichletBoundaryConditionWithinTimeInterval(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
     NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
-- 
GitLab