From 2d891b068ba0864ba8826818085de8c2cab808b5 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Thu, 7 Jul 2016 17:22:15 +0200
Subject: [PATCH] [PL] moved index computation to local assembler interface

---
 .../GroundwaterFlow/GroundwaterFlowFEM.h      | 18 +++----
 ProcessLib/LocalAssemblerInterface.cpp        | 51 +++++++++++++++++++
 ProcessLib/LocalAssemblerInterface.h          | 34 ++++++++-----
 ProcessLib/TES/TESLocalAssembler-impl.h       | 20 +++-----
 ProcessLib/TES/TESLocalAssembler.h            |  9 ++--
 5 files changed, 90 insertions(+), 42 deletions(-)
 create mode 100644 ProcessLib/LocalAssemblerInterface.cpp

diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
index 1696df68044..6f309631377 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -12,7 +12,6 @@
 
 #include <vector>
 
-#include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
@@ -75,16 +74,13 @@ public:
         assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
     }
 
-    void assemble(std::size_t const mesh_item_id,
-                  NumLib::LocalToGlobalIndexMap const& dof_table,
-                  double const /*t*/, GlobalVector const& x,
-                  GlobalMatrix& /*M*/, GlobalMatrix& K,
-                  GlobalVector& b) override
+    void assembleConcrete(
+        double const /*t*/, std::vector<double> const& local_x,
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+        GlobalMatrix& /*M*/, GlobalMatrix& K, GlobalVector& b) override
     {
         _localA.setZero();
         _localRhs.setZero();
-        auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
-        auto const local_x = x.get(indices);
 
         IntegrationMethod integration_method(_integration_order);
         unsigned const n_integration_points = integration_method.getNumberOfPoints();
@@ -108,10 +104,8 @@ public:
             }
         }
 
-        auto const r_c_indices =
-            NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
-        K.add(r_c_indices, _localA);
-        b.add(indices, _localRhs);
+        K.add(indices, _localA);
+        b.add(indices.rows, _localRhs);
     }
 
     Eigen::Map<const Eigen::RowVectorXd>
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
new file mode 100644
index 00000000000..6b78e5ed89e
--- /dev/null
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -0,0 +1,51 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "LocalAssemblerInterface.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+
+namespace ProcessLib
+{
+void LocalAssemblerInterface::assemble(
+    const std::size_t mesh_item_id,
+    const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
+    const GlobalVector& x, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
+{
+    auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
+    auto const local_x = x.get(indices);
+    auto const r_c_indices =
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
+
+    assembleConcrete(t, local_x, r_c_indices, M, K, b);
+}
+
+void LocalAssemblerInterface::assembleJacobian(
+    const std::size_t mesh_item_id,
+    const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
+    const GlobalVector& x, GlobalMatrix& Jac)
+{
+    auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
+    auto const local_x = x.get(indices);
+    auto const r_c_indices =
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
+
+    assembleJacobianConcrete(t, local_x, r_c_indices, Jac);
+}
+
+void LocalAssemblerInterface::assembleJacobianConcrete(
+        double const /*t*/, std::vector<double> const& /*local_x*/,
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& /*indices*/,
+        GlobalMatrix& /*Jac*/)
+{
+    OGS_FATAL(
+        "The assembleJacobian() function is not implemented in the local "
+        "assembler.");
+}
+
+}  // namespace ProcessLib
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 3b4b0e7cac5..37b25e48beb 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -10,6 +10,9 @@
 #ifndef PROCESSLIB_LOCALASSEMBLERINTERFACE_H
 #define PROCESSLIB_LOCALASSEMBLERINTERFACE_H
 
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/NumericsConfig.h"
+
 namespace ProcessLib
 {
 
@@ -23,27 +26,32 @@ class LocalAssemblerInterface
 public:
     virtual ~LocalAssemblerInterface() = default;
 
-    virtual void assemble(std::size_t const mesh_item_id,
+    void assemble(std::size_t const mesh_item_id,
+                  NumLib::LocalToGlobalIndexMap const& dof_table,
+                  double const t, GlobalVector const& x, GlobalMatrix& M,
+                  GlobalMatrix& K, GlobalVector& b);
+
+    void assembleJacobian(std::size_t const mesh_item_id,
                           NumLib::LocalToGlobalIndexMap const& dof_table,
                           double const t, GlobalVector const& x,
-                          GlobalMatrix& M, GlobalMatrix& K,
-                          GlobalVector& b) = 0;
-
-    virtual void assembleJacobian(
-        std::size_t const /*id*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/, double const /*t*/,
-        GlobalVector const& /*x*/, GlobalMatrix& /*Jac*/)
-    {
-        OGS_FATAL(
-            "The assembleJacobian() function is not implemented in the local "
-            "assembler.");
-    }
+                          GlobalMatrix& Jac);
 
     virtual void preTimestep(std::vector<double> const& /*local_x*/,
                              double const /*t*/, double const /*delta_t*/)
     {
     }
     virtual void postTimestep(std::vector<double> const& /*local_x*/) {}
+
+protected:
+    virtual void assembleConcrete(
+            double const t, std::vector<double> const& local_x,
+            NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+            GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) = 0;
+
+    virtual void assembleJacobianConcrete(
+        double const t, std::vector<double> const& local_x,
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+        GlobalMatrix& Jac);
 };
 
 } // namespace ProcessLib
diff --git a/ProcessLib/TES/TESLocalAssembler-impl.h b/ProcessLib/TES/TESLocalAssembler-impl.h
index 7e96f6adb00..331f67b6352 100644
--- a/ProcessLib/TES/TESLocalAssembler-impl.h
+++ b/ProcessLib/TES/TESLocalAssembler-impl.h
@@ -10,7 +10,6 @@
 #define PROCESS_LIB_TES_FEM_IMPL_H_
 
 #include "MaterialsLib/Adsorption/Adsorption.h"
-#include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "NumLib/Function/Interpolation.h"
@@ -133,16 +132,15 @@ TESLocalAssembler<
 
 template <typename ShapeFunction_, typename IntegrationMethod_,
           unsigned GlobalDim>
-void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::assemble(
-    std::size_t const mesh_item_id,
-    NumLib::LocalToGlobalIndexMap const& dof_table, double const /*t*/,
-    GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
+void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
+    assembleConcrete(
+        double const /*t*/, std::vector<double> const& local_x,
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+        GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     _local_M.setZero();
     _local_K.setZero();
     _local_b.setZero();
-    auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
-    auto const local_x = x.get(indices);
 
     IntegrationMethod_ integration_method(_integration_order);
     unsigned const n_integration_points = integration_method.getNumberOfPoints();
@@ -187,11 +185,9 @@ void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::assemble(
         std::printf("\n");
     }
 
-    auto const r_c_indices =
-        NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
-    M.add(r_c_indices, _local_M);
-    K.add(r_c_indices, _local_K);
-    b.add(indices, _local_b);
+    M.add(indices, _local_M);
+    K.add(indices, _local_K);
+    b.add(indices.rows, _local_b);
 }
 
 template <typename ShapeFunction_, typename IntegrationMethod_,
diff --git a/ProcessLib/TES/TESLocalAssembler.h b/ProcessLib/TES/TESLocalAssembler.h
index f21ae2b37d7..37ce25c3045 100644
--- a/ProcessLib/TES/TESLocalAssembler.h
+++ b/ProcessLib/TES/TESLocalAssembler.h
@@ -49,11 +49,10 @@ public:
                       unsigned const integration_order,
                       AssemblyParams const& asm_params);
 
-    void assemble(std::size_t const mesh_item_id,
-                  NumLib::LocalToGlobalIndexMap const& dof_table,
-                  double const t, GlobalVector const& x,
-                  GlobalMatrix& M, GlobalMatrix& K,
-                  GlobalVector& b) override;
+    void assembleConcrete(
+        double const t, std::vector<double> const& local_x,
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+        GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override;
 
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
-- 
GitLab