Merge branch 'ExtrapolatePostReactionSolutionToNodes' into 'master'

[CL] Chemical calculation shifted from nodes to int_pts (Part I): Extrapolate post reaction solution to nodes

See merge request ogs/ogs!3017

ProcessLib/ComponentTransport/ChemicalProcessData.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <vector>
+
+#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
+
+namespace ProcessLib
+{
+namespace ComponentTransport
+{
+struct ChemicalProcessData
+{
+    ChemicalProcessData(
+        std::vector<std::vector<GlobalIndexType>>& chemical_system_index_map_)
+        : chemical_system_index_map(chemical_system_index_map_)
+    {
+    }
+
+    std::vector<std::vector<GlobalIndexType>>& chemical_system_index_map;
+};
+}  // namespace ComponentTransport
+}  // namespace ProcessLib

ProcessLib/ComponentTransport/ComponentTransportFEM.h

@@ -70,6 +70,12 @@ public:
         std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
         std::vector<double>& cache) const = 0;
 
+    virtual std::vector<double> const& getInterpolatedLocalSolution(
+        const double /*t*/,
+        std::vector<GlobalVector*> const& int_pt_x,
+        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
 protected:
     CoupledSolutionsForStaggeredScheme* _coupled_solutions{nullptr};
 };
@@ -891,6 +897,24 @@ public:
         return flux;
     }
 
+    std::vector<double> const& getInterpolatedLocalSolution(
+        const double /*t*/,
+        std::vector<GlobalVector*> const& int_pt_x,
+        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+        std::vector<double>& cache) const override
+    {
+        assert(_process_data.chemical_process_data);
+        assert(int_pt_x.size() == 1);
+
+        cache.clear();
+        auto const ele_id = _element.getID();
+        auto const& indices = _process_data.chemical_process_data
+                                  ->chemical_system_index_map[ele_id];
+        cache = int_pt_x[0]->get(indices);
+
+        return cache;
+    }
+
 private:
     MeshLib::Element const& _element;
     ComponentTransportProcessData const& _process_data;

ProcessLib/ComponentTransport/ComponentTransportProcess.cpp

@@ -174,6 +174,35 @@ void ComponentTransportProcess::
         _local_assemblers, pv.getActiveElementIDs(), _coupled_solutions);
 }
 
+void ComponentTransportProcess::extrapolateIntegrationPointValuesToNodes(
+    const double t,
+    std::vector<GlobalVector*> const& integration_point_values_vectors,
+    std::vector<GlobalVector*>& nodal_values_vectors)
+{
+    auto& extrapolator = getExtrapolator();
+    auto const extrapolatables =
+        NumLib::makeExtrapolatable(_local_assemblers,
+                                   &ComponentTransportLocalAssemblerInterface::
+                                       getInterpolatedLocalSolution);
+
+    for (unsigned transport_process_id = 0;
+         transport_process_id < integration_point_values_vectors.size();
+         ++transport_process_id)
+    {
+        auto const& pv = _process_variables[transport_process_id + 1][0].get();
+        auto const& int_pt_C =
+            integration_point_values_vectors[transport_process_id];
+
+        extrapolator.extrapolate(pv.getNumberOfComponents(), extrapolatables, t,
+                                 {int_pt_C},
+                                 {_local_to_global_index_map.get()});
+
+        auto const& nodal_values = extrapolator.getNodalValues();
+        MathLib::LinAlg::copy(nodal_values,
+                              *nodal_values_vectors[transport_process_id + 1]);
+    }
+}
+
 void ComponentTransportProcess::preTimestepConcreteProcess(
     std::vector<GlobalVector*> const& x, const double /*t*/,
     const double /*delta_t*/, int const process_id)

ProcessLib/ComponentTransport/ComponentTransportProcess.h

@@ -117,6 +117,11 @@ public:
     void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
         int const process_id) override;
 
+    void extrapolateIntegrationPointValuesToNodes(
+        const double t,
+        std::vector<GlobalVector*> const& integration_point_values_vectors,
+        std::vector<GlobalVector*>& nodal_values_vectors) override;
+
     void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
                                     const double /*t*/,
                                     const double /*delta_t*/,

ProcessLib/ComponentTransport/ComponentTransportProcessData.h

@@ -12,6 +12,7 @@
 
 #include <memory>
 
+#include "ChemicalProcessData.h"
 #include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 
@@ -43,6 +44,7 @@ struct ComponentTransportProcessData
     Eigen::VectorXd const specific_body_force;
     bool const has_gravity;
     bool const non_advective_form;
+    std::unique_ptr<ChemicalProcessData> chemical_process_data;
 };
 
 }  // namespace ComponentTransport

ProcessLib/Process.h

@@ -100,6 +100,14 @@ public:
         int const /*process_id*/)
     {
     }
+
+    virtual void extrapolateIntegrationPointValuesToNodes(
+        const double /*t*/,
+        std::vector<GlobalVector*> const& /*integration_point_values_vectors*/,
+        std::vector<GlobalVector*>& /*nodal_values_vectors*/)
+    {
+    }
+
     void preAssemble(const double t, double const dt,
                      GlobalVector const& x) final;
     void assemble(const double t, double const dt,