From 92f9d06eab08f1de4dd3cec8bf0db0942e244c80 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Sun, 21 Jul 2024 12:33:01 +0200
Subject: [PATCH] [PL] Separate matrix output for M,K,b and J,b

Same as 3e9a70cd397af307a97011096814bcc3e4deae2a but for global matrices.

A little code duplication in favor not to deal with pointers.
---
 ProcessLib/Assembly/MatrixOutput.cpp          | 47 ++++++++++++++-----
 ProcessLib/Assembly/MatrixOutput.h            |  8 ++--
 .../ParallelVectorMatrixAssembler.cpp         |  2 +-
 .../LargeDeformationProcess.cpp               |  4 +-
 .../SmallDeformationProcess.cpp               |  4 +-
 .../ThermalTwoPhaseFlowWithPPProcess.cpp      |  4 +-
 6 files changed, 47 insertions(+), 22 deletions(-)

diff --git a/ProcessLib/Assembly/MatrixOutput.cpp b/ProcessLib/Assembly/MatrixOutput.cpp
index 56ed505c69f..5aabd69dc9f 100644
--- a/ProcessLib/Assembly/MatrixOutput.cpp
+++ b/ProcessLib/Assembly/MatrixOutput.cpp
@@ -175,10 +175,9 @@ GlobalMatrixOutput::GlobalMatrixOutput()
 }
 
 void GlobalMatrixOutput::operator()(double const t, int const process_id,
-                                    GlobalMatrix const* M,
-                                    GlobalMatrix const* K,
-                                    GlobalVector const& b,
-                                    GlobalMatrix const* const Jac)
+                                    GlobalMatrix const& M,
+                                    GlobalMatrix const& K,
+                                    GlobalVector const& b)
 {
     if (!do_output_)
     {
@@ -188,22 +187,20 @@ void GlobalMatrixOutput::operator()(double const t, int const process_id,
 #ifndef USE_PETSC
     ++counter_;
 
-    if (M)
     {
         auto fh = openGlobalMatrixOutputFile(filenamePrefix_, counter_, t,
                                              process_id, "M", "mat");
 
         fh << "M ";
-        outputGlobalMatrix(*M, fh);
+        outputGlobalMatrix(M, fh);
     }
 
-    if (K)
     {
         auto fh = openGlobalMatrixOutputFile(filenamePrefix_, counter_, t,
                                              process_id, "K", "mat");
 
         fh << "K ";
-        outputGlobalMatrix(*K, fh);
+        outputGlobalMatrix(K, fh);
     }
 
     {
@@ -213,21 +210,47 @@ void GlobalMatrixOutput::operator()(double const t, int const process_id,
         fh << "b ";
         outputGlobalVector(b, fh);
     }
+#else
+    // do nothing, warning message already printed in the constructor
+    (void)t;
+    (void)process_id;
+    (void)M;
+    (void)K;
+    (void)b;
+#endif
+}
+
+void GlobalMatrixOutput::operator()(double const t, int const process_id,
+                                    GlobalVector const& b,
+                                    GlobalMatrix const& Jac)
+{
+    if (!do_output_)
+    {
+        return;
+    }
+
+#ifndef USE_PETSC
+    ++counter_;
+
+    {
+        auto fh = openGlobalMatrixOutputFile(filenamePrefix_, counter_, t,
+                                             process_id, "b", "vec");
+
+        fh << "b ";
+        outputGlobalVector(b, fh);
+    }
 
-    if (Jac)
     {
         auto fh = openGlobalMatrixOutputFile(filenamePrefix_, counter_, t,
                                              process_id, "Jac", "mat");
 
         fh << "Jac ";
-        outputGlobalMatrix(*Jac, fh);
+        outputGlobalMatrix(Jac, fh);
     }
 #else
     // do nothing, warning message already printed in the constructor
     (void)t;
     (void)process_id;
-    (void)M;
-    (void)K;
     (void)b;
     (void)Jac;
 #endif
diff --git a/ProcessLib/Assembly/MatrixOutput.h b/ProcessLib/Assembly/MatrixOutput.h
index 2b3d285f7b1..855dbf33475 100644
--- a/ProcessLib/Assembly/MatrixOutput.h
+++ b/ProcessLib/Assembly/MatrixOutput.h
@@ -24,9 +24,11 @@ struct GlobalMatrixOutput
 {
     GlobalMatrixOutput();
 
-    void operator()(double const t, int const process_id, GlobalMatrix const* M,
-                    GlobalMatrix const* K, GlobalVector const& b,
-                    GlobalMatrix const* const Jac = nullptr);
+    void operator()(double const t, int const process_id, GlobalMatrix const& M,
+                    GlobalMatrix const& K, GlobalVector const& b);
+
+    void operator()(double const t, int const process_id, GlobalVector const& b,
+                    GlobalMatrix const& Jac);
 
 private:
     std::string filenamePrefix_;
diff --git a/ProcessLib/Assembly/ParallelVectorMatrixAssembler.cpp b/ProcessLib/Assembly/ParallelVectorMatrixAssembler.cpp
index e799f5153ce..5b070cb55db 100644
--- a/ProcessLib/Assembly/ParallelVectorMatrixAssembler.cpp
+++ b/ProcessLib/Assembly/ParallelVectorMatrixAssembler.cpp
@@ -249,7 +249,7 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
 
     stats->print();
 
-    global_matrix_output_(t, process_id, nullptr, nullptr, b, &Jac);
+    global_matrix_output_(t, process_id, b, Jac);
     exception.rethrow();
 }
 }  // namespace ProcessLib::Assembly
diff --git a/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp b/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp
index bf5a13101ff..f6c58710518 100644
--- a/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp
+++ b/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp
@@ -145,7 +145,7 @@ void LargeDeformationProcess<DisplacementDim>::assembleConcreteProcess(
         getActiveElementIDs(), dof_table, t, dt, x, x_prev, process_id, &M, &K,
         &b);
 
-    _global_output(t, process_id, &M, &K, b);
+    _global_output(t, process_id, M, K, b);
 }
 
 template <int DisplacementDim>
@@ -168,7 +168,7 @@ void LargeDeformationProcess<DisplacementDim>::
     transformVariableFromGlobalVector(b, 0, *_local_to_global_index_map,
                                       *_nodal_forces, std::negate<double>());
 
-    _global_output(t, process_id, nullptr, nullptr, b, &Jac);
+    _global_output(t, process_id, b, Jac);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
index b588aef7a42..8e9ea109e51 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
@@ -147,7 +147,7 @@ void SmallDeformationProcess<DisplacementDim>::assembleConcreteProcess(
         getActiveElementIDs(), dof_table, t, dt, x, x_prev, process_id, &M, &K,
         &b);
 
-    _global_output(t, process_id, &M, &K, b);
+    _global_output(t, process_id, M, K, b);
 }
 
 template <int DisplacementDim>
@@ -170,7 +170,7 @@ void SmallDeformationProcess<DisplacementDim>::
     transformVariableFromGlobalVector(b, 0, *_local_to_global_index_map,
                                       *_nodal_forces, std::negate<double>());
 
-    _global_output(t, process_id, nullptr, nullptr, b, &Jac);
+    _global_output(t, process_id, b, Jac);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
index 20165029fdd..a7df17686ff 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
@@ -96,7 +96,7 @@ void ThermalTwoPhaseFlowWithPPProcess::assembleConcreteProcess(
         getActiveElementIDs(), dof_table, t, dt, x, x_prev, process_id, &M, &K,
         &b);
 
-    _global_output(t, process_id, &M, &K, b);
+    _global_output(t, process_id, M, K, b);
 }
 
 void ThermalTwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
@@ -115,7 +115,7 @@ void ThermalTwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
         _local_assemblers, getActiveElementIDs(), dof_table, t, dt, x, x_prev,
         process_id, &b, &Jac);
 
-    _global_output(t, process_id, nullptr, nullptr, b, &Jac);
+    _global_output(t, process_id, b, Jac);
 }
 
 void ThermalTwoPhaseFlowWithPPProcess::preTimestepConcreteProcess(
-- 
GitLab