Merge pull request #2272 from wenqing/time_dependent_D_BC

Dirichlet boundary condition within a time interval

BaseLib/TimeInterval.cpp

@@ -0,0 +1,36 @@
+/**
+ * \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:   TimeInterval.cpp
+ *
+ * Created on November 27, 2018, 5:06 PM
+ *
+ */
+
+#include "TimeInterval.h"
+
+#include "BaseLib/ConfigTree.h"
+
+namespace BaseLib
+{
+std::unique_ptr<TimeInterval> createTimeInterval(
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{prj__time_loop__processes__process__time_interval}
+    auto const& time_interval_config = config.getConfigSubtree("time_interval");
+
+    const double start_time =
+        //! \ogs_file_param{prj__time_loop__processes__process__time_interval__start}
+        time_interval_config.getConfigParameter<double>("start");
+
+    const double end_time =
+        //! \ogs_file_param{prj__time_loop__processes__process__time_interval__end}
+        time_interval_config.getConfigParameter<double>("end");
+
+    return std::make_unique<BaseLib::TimeInterval>(start_time, end_time);
+}
+}  // end of namespace

BaseLib/TimeInterval.h

@@ -0,0 +1,46 @@
+/**
+ *
+ * \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:   TimeInterval.h
+ *
+ * Created on November 26, 2018, 4:44 PM
+ */
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+
+/*!
+ * Class for a time interval, which has a member to check whether the given time
+ * is in this time interval.
+ */
+class TimeInterval final
+{
+public:
+    TimeInterval(const double start_time, const double end_time)
+        : _start_time(start_time), _end_time(end_time)
+    {
+    }
+
+    bool contains(const double current_time) const
+    {
+        return (current_time >= _start_time && current_time <= _end_time);
+    }
+
+private:
+    const double _start_time;
+    const double _end_time;
+};
+
+std::unique_ptr<TimeInterval> createTimeInterval(
+    BaseLib::ConfigTree const& config);
+
+}  // end of namespace

Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/DirichletWithinTimeInterval/c_DirichletWithinTimeInterval.md

@@ -0,0 +1,4 @@
+It defines a Dirichlet boundary condition that exists only in a time interval.
+
+Benchmark:
+Tests/Data/Parabolic/LiquidFlow/TimeIntervalDirichletBC/TimeIntervalDirichletBC.prj

Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/DirichletWithinTimeInterval/t_parameter.md

@@ -0,0 +1,2 @@
+The value of the Dirichlet boundary condition within the
+ specified time interval.

Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/DirichletWithinTimeInterval/t_time_interval.md

@@ -0,0 +1 @@
+\copydoc ogs_file_param__prj__time_loop__processes__process__time_interval

Documentation/ProjectFile/prj/time_loop/processes/process/time_interval/i_time_interval.md

@@ -0,0 +1 @@
+It defines a time interval.

Documentation/ProjectFile/prj/time_loop/processes/process/time_interval/t_end.md

@@ -0,0 +1 @@
+The end time of the time interval.

Documentation/ProjectFile/prj/time_loop/processes/process/time_interval/t_start.md

@@ -0,0 +1 @@
+The start time of the time interval.

ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp

@@ -13,13 +13,17 @@
 #include "BoundaryConditionConfig.h"
 #include "ConstraintDirichletBoundaryCondition.h"
 #include "DirichletBoundaryCondition.h"
+#include "DirichletBoundaryConditionWithinTimeInterval.h"
 #include "NeumannBoundaryCondition.h"
 #include "NonuniformDirichletBoundaryCondition.h"
 #include "NonuniformNeumannBoundaryCondition.h"
 #include "NonuniformVariableDependentNeumannBoundaryCondition.h"
 #include "NormalTractionBoundaryCondition.h"
 #include "PhaseFieldIrreversibleDamageOracleBoundaryCondition.h"
 #include "RobinBoundaryCondition.h"
+
+#include "BaseLib/TimeInterval.h"
+
 #ifdef OGS_USE_PYTHON
 #include "Python/PythonBoundaryCondition.h"
 #endif
@@ -54,6 +58,12 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
             config.config, config.boundary_mesh, dof_table, variable_id,
             *config.component_id, parameters);
     }
+    if (type == "DirichletWithinTimeInterval")
+    {
+        return ProcessLib::createDirichletBoundaryConditionWithinTimeInterval(
+            config.config, config.boundary_mesh, dof_table, variable_id,
+            *config.component_id, parameters);
+    }
     if (type == "Neumann")
     {
         return ProcessLib::createNeumannBoundaryCondition(

ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp

@@ -12,46 +12,44 @@
 #include <algorithm>
 #include <logog/include/logog.hpp>
 #include <vector>
+
+#include "DirichletBoundaryConditionAuxiliaryFunctions.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "NumLib/IndexValueVector.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/Parameter/Parameter.h"
 #include "ProcessLib/Utils/ProcessUtils.h"
 
 namespace ProcessLib
 {
-void DirichletBoundaryCondition::getEssentialBCValues(
-    const double t, GlobalVector const& /*x*/,
-    NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
+DirichletBoundaryCondition::DirichletBoundaryCondition(
+    Parameter<double> const& parameter, MeshLib::Mesh const& bc_mesh,
+    NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
+    int const component_id)
+    : _parameter(parameter),
+      _bc_mesh(bc_mesh),
+      _variable_id(variable_id),
+      _component_id(component_id)
 {
-    SpatialPosition pos;
+    checkParametersOfDirichletBoundaryCondition(_bc_mesh, dof_table_bulk,
+                                                _variable_id, _component_id);
 
-    bc_values.ids.clear();
-    bc_values.values.clear();
+    std::vector<MeshLib::Node*> const& bc_nodes = bc_mesh.getNodes();
+    MeshLib::MeshSubset bc_mesh_subset(_bc_mesh, bc_nodes);
 
-    // 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());
-        }
-    }
+    // 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 DirichletBoundaryCondition::getEssentialBCValues(
+    const double t, GlobalVector const& x,
+    NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
+{
+    getEssentialBCValuesLocal(_parameter, _bc_mesh, *_dof_table_boundary,
+                              _variable_id, _component_id, t, x, bc_values);
 }
 
 std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(
@@ -70,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.

ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h

@@ -9,13 +9,18 @@
 
 #pragma once
 
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "NumLib/IndexValueVector.h"
-#include "ProcessLib/Parameter/Parameter.h"
 #include "BoundaryCondition.h"
 
+namespace BaseLib
+{
+class ConfigTree;
+}
+
 namespace ProcessLib
 {
+template <typename T>
+struct Parameter;
+
 // TODO docu
 /// The DirichletBoundaryCondition class describes a constant in space
 /// and time Dirichlet boundary condition.
@@ -27,47 +32,7 @@ class DirichletBoundaryCondition final : public BoundaryCondition
     DirichletBoundaryCondition(
         Parameter<double> const& parameter, MeshLib::Mesh const& bc_mesh,
         NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
-        int const variable_id, int const component_id)
-        : _parameter(parameter),
-          _bc_mesh(bc_mesh),
-          _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);
-
-        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)));
-    }
+        int const variable_id, int const component_id);
 
     void getEssentialBCValues(
         const double t, GlobalVector const& x,

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