fixup HeatConduction

ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp

@@ -20,6 +20,7 @@ namespace HeatConduction
 {
 std::unique_ptr<Process> createHeatConductionProcess(
     MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
     std::vector<ProcessVariable> const& variables,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     BaseLib::ConfigTree const& config)
@@ -72,7 +73,8 @@ std::unique_ptr<Process> createHeatConductionProcess(
                                         named_function_caller);
 
     return std::unique_ptr<Process>{new HeatConductionProcess{
-        mesh, parameters, std::move(process_variables), std::move(process_data),
+        mesh, std::move(jacobian_assembler), parameters,
+        std::move(process_variables), std::move(process_data),
         std::move(secondary_variables), std::move(named_function_caller)}};
 }
 

ProcessLib/HeatConduction/CreateHeatConductionProcess.h

@@ -19,6 +19,7 @@ namespace HeatConduction
 {
 std::unique_ptr<Process> createHeatConductionProcess(
     MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
     std::vector<ProcessVariable> const& variables,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     BaseLib::ConfigTree const& config);

ProcessLib/HeatConduction/HeatConductionFEM.h

@@ -12,6 +12,8 @@
 
 #include <vector>
 
+#include "HeatConductionProcessData.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
@@ -19,7 +21,6 @@
 #include "ProcessLib/LocalAssemblerTraits.h"
 #include "ProcessLib/Parameter/Parameter.h"
 #include "ProcessLib/Utils/InitShapeMatrices.h"
-#include "HeatConductionProcessData.h"
 
 namespace ProcessLib
 {
@@ -65,10 +66,6 @@ public:
                        HeatConductionProcessData const& process_data)
         : _element(element),
           _process_data(process_data),
-          _localK(local_matrix_size,
-                  local_matrix_size),  // TODO narrowing conversion
-          _localM(local_matrix_size, local_matrix_size),
-          _localRhs(local_matrix_size),
           _integration_method(integration_order),
           _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
                                             IntegrationMethod, GlobalDim>(
@@ -79,14 +76,19 @@ public:
         assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
     }
 
-    void assembleConcrete(
-        double const t, std::vector<double> const& local_x,
-        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
-        GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override
+    void assemble(double const t, std::vector<double> const& local_x,
+                  std::vector<double>& local_M_data,
+                  std::vector<double>& local_K_data,
+                  std::vector<double>& /*local_b_data*/) override
     {
-        _localK.setZero();
-        _localM.setZero();
-        _localRhs.setZero();
+        auto const local_matrix_size = local_x.size();
+        // This assertion is valid only if all nodal d.o.f. use the same shape matrices.
+        assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+
+        auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
+            local_M_data, local_matrix_size, local_matrix_size);
+        auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
+            local_K_data, local_matrix_size, local_matrix_size);
 
         unsigned const n_integration_points =
             _integration_method.getNumberOfPoints();
@@ -103,9 +105,9 @@ public:
             auto const heat_capacity = _process_data.heat_capacity(t, pos)[0];
             auto const density = _process_data.density(t, pos)[0];
 
-            _localK.noalias() +=
+            local_K.noalias() +=
                 sm.dNdx.transpose() * k * sm.dNdx * sm.detJ * wp.getWeight();
-            _localM.noalias() += sm.N.transpose() * density * heat_capacity *
+            local_M.noalias() += sm.N.transpose() * density * heat_capacity *
                                  sm.N * sm.detJ * wp.getWeight();
             // heat flux only computed for output.
             GlobalDimVectorType const heat_flux =
@@ -117,10 +119,6 @@ public:
                 _heat_fluxes[d][ip] = heat_flux[d];
             }
         }
-
-        K.add(indices, _localK);
-        M.add(indices, _localM);
-        b.add(indices.rows, _localRhs);
     }
 
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
@@ -157,10 +155,6 @@ private:
     MeshLib::Element const& _element;
     HeatConductionProcessData const& _process_data;
 
-    NodalMatrixType _localK;
-    NodalMatrixType _localM;
-    NodalVectorType _localRhs;
-
     IntegrationMethod const _integration_method;
     std::vector<ShapeMatrices> _shape_matrices;
 

ProcessLib/HeatConduction/HeatConductionProcess.cpp

@@ -19,13 +19,15 @@ namespace HeatConduction
 {
 HeatConductionProcess::HeatConductionProcess(
     MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
     HeatConductionProcessData&& process_data,
     SecondaryVariableCollection&& secondary_variables,
     NumLib::NamedFunctionCaller&& named_function_caller)
-    : Process(mesh, parameters, std::move(process_variables),
-              std::move(secondary_variables), std::move(named_function_caller)),
+    : Process(mesh, std::move(jacobian_assembler), parameters,
+              std::move(process_variables), std::move(secondary_variables),
+              std::move(named_function_caller)),
       _process_data(std::move(process_data))
 {
 }
@@ -72,10 +74,24 @@ void HeatConductionProcess::assembleConcreteProcess(const double t,
     DBUG("Assemble HeatConductionProcess.");
 
     // Call global assembler for each local assembly item.
-    GlobalExecutor::executeMemberOnDereferenced(
-        &HeatConductionLocalAssemblerInterface::assemble, _local_assemblers,
+    GlobalExecutor::executeMemberDereferenced(
+        _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         *_local_to_global_index_map, t, x, M, K, b);
 }
 
+void HeatConductionProcess::assembleWithJacobianConcreteProcess(
+    const double t, GlobalVector const& x, GlobalVector const& xdot,
+    const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
+    GlobalVector& b, GlobalMatrix& Jac)
+{
+    DBUG("AssembleWithJacobian HeatConductionProcess.");
+
+    // Call global assembler for each local assembly item.
+    GlobalExecutor::executeMemberDereferenced(
+        _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac);
+}
+
 }  // namespace HeatConduction
 }  // namespace ProcessLib

ProcessLib/HeatConduction/HeatConductionProcess.h

@@ -25,6 +25,8 @@ class HeatConductionProcess final : public Process
 public:
     HeatConductionProcess(
         MeshLib::Mesh& mesh,
+        std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+            jacobian_assembler,
         std::vector<std::unique_ptr<ParameterBase>> const& parameters,
         std::vector<std::reference_wrapper<ProcessVariable>>&&
             process_variables,
@@ -48,6 +50,11 @@ private:
                                  GlobalMatrix& M, GlobalMatrix& K,
                                  GlobalVector& b) override;
 
+    void assembleWithJacobianConcreteProcess(
+        const double t, GlobalVector const& x, GlobalVector const& xdot,
+        const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+
     HeatConductionProcessData _process_data;
 
     std::vector<std::unique_ptr<HeatConductionLocalAssemblerInterface>>