From 7b2af516c3498a7c0bdfa7ff9dc315c0f9136e89 Mon Sep 17 00:00:00 2001
From: renchao_lu <renchao.lu@gmail.com>
Date: Wed, 23 Dec 2020 12:14:21 +0100
Subject: [PATCH] [PL/CT] Set chemical_system_id.

---
 .../ComponentTransportFEM.h                   | 34 +++++++++++--------
 .../ComponentTransportProcess.cpp             |  4 +++
 2 files changed, 23 insertions(+), 15 deletions(-)

diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 42524cd0a60..957e5c8bbc8 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -48,6 +48,8 @@ struct IntegrationPointData final
     GlobalDimNodalMatrixType const dNdx;
     double const integration_weight;
 
+    GlobalIndexType chemical_system_id = 0;
+
     double porosity = std::numeric_limits<double>::quiet_NaN();
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
 };
@@ -66,6 +68,8 @@ public:
         _coupled_solutions = coupling_term;
     }
 
+    virtual void setChemicalSystemID(std::size_t const /*mesh_item_id*/) = 0;
+
     virtual std::vector<double> const& getIntPtDarcyVelocity(
         const double t,
         std::vector<GlobalVector*> const& x,
@@ -159,23 +163,23 @@ public:
                 medium->property(MaterialPropertyLib::PropertyType::porosity)
                     .template initialValue<double>(pos, 0.);
         }
+    }
 
-        if (_process_data.chemical_process_data)
-        {
-            // chemical system index map
-            auto& chemical_system_index_map =
-                _process_data.chemical_process_data->chemical_system_index_map;
-
-            GlobalIndexType const start_value =
-                chemical_system_index_map.empty()
-                    ? 0
-                    : chemical_system_index_map.back().back() + 1;
-
-            std::vector<GlobalIndexType> indices(
-                _integration_method.getNumberOfPoints());
-            std::iota(indices.begin(), indices.end(), start_value);
+    void setChemicalSystemID(std::size_t const /*mesh_item_id*/) override
+    {
+        assert(_process_data.chemical_solver_interface);
+        // chemical system index map
+        auto& chemical_system_index_map =
+            _process_data.chemical_solver_interface->chemical_system_index_map;
 
-            chemical_system_index_map.push_back(std::move(indices));
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            _ip_data[ip].chemical_system_id = chemical_system_index_map.empty()
+                                     ? 0
+                                     : chemical_system_index_map.back() + 1;
+            chemical_system_index_map.push_back(_ip_data[ip].chemical_system_id);
         }
     }
 
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index b0f217f60f5..f6c880a3537 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -86,6 +86,10 @@ void ComponentTransportProcess::initializeConcreteProcess(
 
     if (_chemical_solver_interface)
     {
+        GlobalExecutor::executeSelectedMemberOnDereferenced(
+            &ComponentTransportLocalAssemblerInterface::setChemicalSystemID,
+            _local_assemblers, pv.getActiveElementIDs());
+
         _chemical_solver_interface->initialize();
     }
 
-- 
GitLab