From a814fff0dc64f87c86497b563ce6a432b6e7eabd Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Fri, 4 May 2018 20:57:34 +0200
Subject: [PATCH] [PL] LIE/HM: Extract copyVariableFromGlobalVector.

---
 .../HydroMechanics/HydroMechanicsProcess.cpp  | 44 ++++-----------
 ProcessLib/Utils/GlobalVectorUtils.h          | 53 +++++++++++++++++++
 2 files changed, 63 insertions(+), 34 deletions(-)
 create mode 100644 ProcessLib/Utils/GlobalVectorUtils.h

diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
index 27b6b32663e..87fed2424bc 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
@@ -20,6 +20,7 @@
 
 #include "ProcessLib/LIE/BoundaryCondition/BoundaryConditionBuilder.h"
 #include "ProcessLib/LIE/Common/MeshUtils.h"
+#include "ProcessLib/Utils/GlobalVectorUtils.h"
 
 #include "LocalAssembler/CreateLocalAssemblers.h"
 #include "LocalAssembler/HydroMechanicsLocalAssemblerFracture.h"
@@ -547,45 +548,20 @@ void HydroMechanicsProcess<GlobalDim>::assembleWithJacobianConcreteProcess(
         _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
         dx_dx, M, K, b, Jac, _coupled_solutions);
 
-    auto copy_variable_from_global_vector = [this, &b](int const variable_id,
-                                                       int const n_components,
-                                                       auto& output_vector) {
-        std::fill(output_vector.begin(), output_vector.end(),
-                  std::numeric_limits<double>::quiet_NaN());
-        for (int component = 0; component < n_components; ++component)
-        {
-            auto const& mesh_subsets =
-                _local_to_global_index_map->getMeshSubsets(variable_id,
-                                                           component);
-            for (auto const& ms : mesh_subsets)
-            {
-                auto const mesh_id = _mesh.getID();
-                assert(ms->getMeshID() ==
-                       mesh_id);  // Multiple meshes not supported.
-                for (auto const& node : ms->getNodes())
-                {
-                    auto const node_id = node->getID();
-                    MeshLib::Location const l(
-                        mesh_id, MeshLib::MeshItemType::Node, node_id);
-                    output_vector.getComponent(node_id, component) =
-                        -b[_local_to_global_index_map->getGlobalIndex(
-                            l, variable_id, component)];
-                }
-            }
-        }
+    auto copyRhs = [&](int const variable_id, auto& output_vector) {
+        transformVariableFromGlobalVector(b, variable_id,
+                                          *_local_to_global_index_map,
+                                          output_vector, std::negate<double>());
     };
-    copy_variable_from_global_vector(0, 1,
-                                     *_process_data.mesh_prop_hydraulic_flow);
-    copy_variable_from_global_vector(1, GlobalDim,
-                                     *_process_data.mesh_prop_nodal_forces);
-    copy_variable_from_global_vector(
-        2, GlobalDim, *_process_data.mesh_prop_nodal_forces_jump);
+    copyRhs(0, *_process_data.mesh_prop_hydraulic_flow);
+    copyRhs(1, *_process_data.mesh_prop_nodal_forces);
+    copyRhs(2, *_process_data.mesh_prop_nodal_forces_jump);
 }
 
 template <int GlobalDim>
 void HydroMechanicsProcess<GlobalDim>::preTimestepConcreteProcess(
-    GlobalVector const& x, double const t,
-    double const dt, const int /*process_id*/)
+    GlobalVector const& x, double const t, double const dt,
+    const int /*process_id*/)
 {
     DBUG("PreTimestep HydroMechanicsProcess.");
 
diff --git a/ProcessLib/Utils/GlobalVectorUtils.h b/ProcessLib/Utils/GlobalVectorUtils.h
new file mode 100644
index 00000000000..906feffe169
--- /dev/null
+++ b/ProcessLib/Utils/GlobalVectorUtils.h
@@ -0,0 +1,53 @@
+/**
+ * \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
+ *
+ */
+
+#include "MeshLib/Properties.h"
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+
+namespace ProcessLib
+{
+/// Copies part of a global vector for the given variable into output vector
+/// while applying a function to each value.
+///
+/// \attention The output_vector is accessed through node id and component,
+/// therefore multiple meshes are not supported.
+template <typename Functor>
+void transformVariableFromGlobalVector(
+    GlobalVector const& input_vector, int const variable_id,
+    NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
+    MeshLib::PropertyVector<double>& output_vector, Functor mapFunction)
+{
+    std::fill(begin(output_vector), end(output_vector),
+              std::numeric_limits<double>::quiet_NaN());
+
+    int const n_components =
+        local_to_global_index_map.getNumberOfVariableComponents(variable_id);
+    for (int component = 0; component < n_components; ++component)
+    {
+        auto const& mesh_subsets =
+            local_to_global_index_map.getMeshSubsets(variable_id, component);
+        assert(mesh_subsets.size() ==
+               1);  // Multiple meshes are not supported by the output_vector.
+        for (auto const& ms : mesh_subsets)
+        {
+            auto const mesh_id = ms->getMeshID();
+            for (auto const& node : ms->getNodes())
+            {
+                auto const node_id = node->getID();
+                MeshLib::Location const l(mesh_id, MeshLib::MeshItemType::Node,
+                                          node_id);
+                output_vector.getComponent(node_id, component) = mapFunction(
+                    input_vector[local_to_global_index_map.getGlobalIndex(
+                        l, variable_id, component)]);
+            }
+        }
+    }
+}
+}  // namespace ProcessLib
-- 
GitLab