From cf6e008bb785a2f6bb08a1020c65d55dbbf7a795 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Tue, 4 Dec 2018 14:40:01 +0100
Subject: [PATCH] [deactivation] Added _deactivated_subdomains to
 ProcessVariable as a member

---
 ProcessLib/Process.cpp         | 17 ++++++
 ProcessLib/ProcessVariable.cpp | 96 ++++++++++++++++++++++++++++------
 ProcessLib/ProcessVariable.h   | 27 ++++++++--
 3 files changed, 119 insertions(+), 21 deletions(-)

diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index cc96693cbd7..d9107b04979 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -59,6 +59,23 @@ Process::Process(
           return pcs_sts;
       }(_process_variables.size()))
 {
+    // If there are deactivated subdomains, check whether MaterialIDs exist in
+    // mesh data.
+    for (auto const& per_process_process_variables : _process_variables)
+    {
+        for (auto const& process_variable : per_process_process_variables)
+        {
+            if ((!process_variable.get()
+                      .getDeactivatedSubdomains()
+                      .empty()) &&
+                (!materialIDs(mesh)))
+            {
+                OGS_FATAL(
+                    "The mesh does not contain matertialIDs for the "
+                    "deactivation of subdomains. The program terminates now.");
+            }
+        }
+    }
 }
 
 void Process::initializeProcessBoundaryConditionsAndSourceTerms(
diff --git a/ProcessLib/ProcessVariable.cpp b/ProcessLib/ProcessVariable.cpp
index dda4286dbcc..0b7fb131da6 100644
--- a/ProcessLib/ProcessVariable.cpp
+++ b/ProcessLib/ProcessVariable.cpp
@@ -9,14 +9,22 @@
 
 #include "ProcessVariable.h"
 
+#include <algorithm>
 #include <utility>
+
 #include <logog/include/logog.hpp>
 
 #include "BaseLib/Algorithm.h"
+#include "BaseLib/TimeInterval.h"
+
 #include "MeshGeoToolsLib/ConstructMeshesFromGeometries.h"
 #include "MeshLib/Mesh.h"
+#include "MeshLib/Node.h"
+
 #include "ProcessLib/BoundaryCondition/BoundaryCondition.h"
 #include "ProcessLib/BoundaryCondition/CreateBoundaryCondition.h"
+#include "ProcessLib/BoundaryCondition/DirichletBoundaryConditionWithinTimeInterval.h"
+#include "ProcessLib/DeactivatedSubdomain.h"
 #include "ProcessLib/SourceTerms/CreateSourceTerm.h"
 #include "ProcessLib/SourceTerms/SourceTerm.h"
 #include "ProcessLib/Utils/ProcessUtils.h"
@@ -97,6 +105,7 @@ ProcessVariable::ProcessVariable(
       _n_components(config.getConfigParameter<int>("components")),
       //! \ogs_file_param{prj__process_variables__process_variable__order}
       _shapefunction_order(config.getConfigParameter<unsigned>("order")),
+      _deactivated_subdomains(createDeactivatedSubdomains(config, mesh)),
       _initial_condition(findParameter<double>(
           //! \ogs_file_param{prj__process_variables__process_variable__initial_condition}
           config.getConfigParameter<std::string>("initial_condition"),
@@ -105,15 +114,18 @@ ProcessVariable::ProcessVariable(
     DBUG("Constructing process variable %s", _name.c_str());
 
     if (_shapefunction_order < 1 || 2 < _shapefunction_order)
-        OGS_FATAL("The given shape function order %d is not supported", _shapefunction_order);
+        OGS_FATAL("The given shape function order %d is not supported",
+                  _shapefunction_order);
 
     // Boundary conditions
     //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions}
-    if (auto bcs_config = config.getConfigSubtreeOptional("boundary_conditions"))
+    if (auto bcs_config =
+            config.getConfigSubtreeOptional("boundary_conditions"))
     {
-        for (auto bc_config :
-             //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition}
-             bcs_config->getConfigSubtreeList("boundary_condition"))
+        for (
+            auto bc_config :
+            //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition}
+            bcs_config->getConfigSubtreeList("boundary_condition"))
         {
             auto const& mesh = findMeshInConfig(bc_config, meshes);
             auto component_id =
@@ -126,7 +138,9 @@ ProcessVariable::ProcessVariable(
 
             _bc_configs.emplace_back(std::move(bc_config), mesh, component_id);
         }
-    } else {
+    }
+    else
+    {
         INFO("No boundary conditions found.");
     }
 
@@ -134,9 +148,10 @@ ProcessVariable::ProcessVariable(
     //! \ogs_file_param{prj__process_variables__process_variable__source_terms}
     if (auto sts_config = config.getConfigSubtreeOptional("source_terms"))
     {
-        for (auto st_config :
-             //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term}
-             sts_config->getConfigSubtreeList("source_term"))
+        for (
+            auto st_config :
+            //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term}
+            sts_config->getConfigSubtreeList("source_term"))
         {
             MeshLib::Mesh const& mesh = findMeshInConfig(st_config, meshes);
             auto component_id =
@@ -150,7 +165,9 @@ ProcessVariable::ProcessVariable(
             _source_term_configs.emplace_back(std::move(st_config), mesh,
                                               component_id);
         }
-    } else {
+    }
+    else
+    {
         INFO("No source terms found.");
     }
 }
@@ -160,6 +177,7 @@ ProcessVariable::ProcessVariable(ProcessVariable&& other)
       _mesh(other._mesh),
       _n_components(other._n_components),
       _shapefunction_order(other._shapefunction_order),
+      _deactivated_subdomains(std::move(other._deactivated_subdomains)),
       _initial_condition(std::move(other._initial_condition)),
       _bc_configs(std::move(other._bc_configs)),
       _source_term_configs(std::move(other._source_term_configs))
@@ -189,10 +207,9 @@ ProcessVariable::createBoundaryConditions(
 
     for (auto& config : _bc_configs)
     {
-        auto bc = createBoundaryCondition(config, dof_table, _mesh, variable_id,
-                                          integration_order,
-                                          _shapefunction_order, parameters,
-                                          process);
+        auto bc = createBoundaryCondition(
+            config, dof_table, _mesh, variable_id, integration_order,
+            _shapefunction_order, parameters, process);
 #ifdef USE_PETSC
         if (bc == nullptr)
         {
@@ -202,11 +219,58 @@ ProcessVariable::createBoundaryConditions(
         bcs.push_back(std::move(bc));
     }
 
+    if (_deactivated_subdomains.empty())
+        return bcs;
+
+    createBoundaryConditionsForDeactivatedSubDomains(dof_table, variable_id,
+                                                     parameters, bcs);
     return bcs;
 }
 
-std::vector<std::unique_ptr<SourceTerm>>
-ProcessVariable::createSourceTerms(
+void ProcessVariable::createBoundaryConditionsForDeactivatedSubDomains(
+    const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    std::vector<std::unique_ptr<BoundaryCondition>>& bcs)
+{
+    auto& parameter = findParameter<double>(
+        DeactivatedSubdomain::name_of_paramater_of_zero, parameters, 1);
+
+    for (auto const& deactivated_subdomain : _deactivated_subdomains)
+    {
+        // Copy the time interval.
+        std::unique_ptr<BaseLib::TimeInterval> time_interval =
+            std::make_unique<BaseLib::TimeInterval>(
+                (*(*deactivated_subdomain).time_interval));
+
+        auto const& deactivated_sudomain_meshes =
+            (*deactivated_subdomain).deactivated_sudomain_meshes;
+        for (auto const& deactivated_sudomain_mesh :
+             deactivated_sudomain_meshes)
+        {
+            for (int component_id = 0;
+                 component_id < dof_table.getNumberOfComponents();
+                 component_id++)
+            {
+                auto bc = std::make_unique<
+                    DirichletBoundaryConditionWithinTimeInterval>(
+                    std::move(time_interval), parameter,
+                    (*(*deactivated_sudomain_mesh).mesh),
+                    (*deactivated_sudomain_mesh).inactive_nodes, dof_table,
+                    variable_id, component_id);
+
+#ifdef USE_PETSC
+                if (bc == nullptr)
+                {
+                    continue;
+                }
+#endif  // USE_PETSC
+                bcs.push_back(std::move(bc));
+            }
+        }
+    }
+}
+
+std::vector<std::unique_ptr<SourceTerm>> ProcessVariable::createSourceTerms(
     const NumLib::LocalToGlobalIndexMap& dof_table,
     const int variable_id,
     unsigned const integration_order,
diff --git a/ProcessLib/ProcessVariable.h b/ProcessLib/ProcessVariable.h
index 51cea3dd78e..2d38f678cad 100644
--- a/ProcessLib/ProcessVariable.h
+++ b/ProcessLib/ProcessVariable.h
@@ -9,6 +9,8 @@
 
 #pragma once
 
+#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
+
 #include "ProcessLib/BoundaryCondition/BoundaryConditionConfig.h"
 #include "ProcessLib/Parameter/Parameter.h"
 #include "ProcessLib/SourceTerms/SourceTermConfig.h"
@@ -16,7 +18,8 @@
 namespace MeshLib
 {
 class Mesh;
-template <typename T> class PropertyVector;
+template <typename T>
+class PropertyVector;
 }
 namespace NumLib
 {
@@ -32,14 +35,15 @@ class Process;
 
 namespace ProcessLib
 {
+struct DeactivatedSubdomain;
+
 /// A named process variable. Its properties includes the mesh, and the initial
 /// and boundary conditions as well as the source terms.
 class ProcessVariable
 {
 public:
     ProcessVariable(
-        BaseLib::ConfigTree const& config,
-        MeshLib::Mesh& mesh,
+        BaseLib::ConfigTree const& config, MeshLib::Mesh& mesh,
         std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
         std::vector<std::unique_ptr<ParameterBase>> const& parameters);
 
@@ -50,9 +54,14 @@ public:
     /// Returns a mesh on which the process variable is defined.
     MeshLib::Mesh const& getMesh() const;
 
+    std::vector<std::unique_ptr<DeactivatedSubdomain const>> const&
+    getDeactivatedSubdomains() const
+    {
+        return _deactivated_subdomains;
+    }
+
     /// Returns the number of components of the process variable.
     int getNumberOfComponents() const { return _n_components; }
-
     std::vector<std::unique_ptr<BoundaryCondition>> createBoundaryConditions(
         const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id,
         unsigned const integration_order,
@@ -70,7 +79,6 @@ public:
     }
 
     unsigned getShapeFunctionOrder() const { return _shapefunction_order; }
-
 private:
     std::string const _name;
     MeshLib::Mesh& _mesh;
@@ -89,6 +97,15 @@ private:
     ///
     /// \sa MeshLib::CellRule MeshLib::FaceRule MeshLib::EdgeRule.
     unsigned _shapefunction_order;
+
+    std::vector<std::unique_ptr<DeactivatedSubdomain const>>
+        _deactivated_subdomains;
+
+    void createBoundaryConditionsForDeactivatedSubDomains(
+        const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id,
+        std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+        std::vector<std::unique_ptr<BoundaryCondition>>& bcs);
+
     Parameter<double> const& _initial_condition;
 
     std::vector<BoundaryConditionConfig> _bc_configs;
-- 
GitLab