From 86b35a15d6bda1404836a4db233e4268e424b477 Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Mon, 19 Apr 2021 12:35:36 +0200
Subject: [PATCH] [PL/Output] Distinguish between bulk and non-bulk meshes.

---
 ProcessLib/Output/Output.cpp        | 26 ++++++++--------
 ProcessLib/Output/ProcessOutput.cpp | 47 +++++++++++++++++++++++++++++
 ProcessLib/Output/ProcessOutput.h   |  3 +-
 3 files changed, 62 insertions(+), 14 deletions(-)

diff --git a/ProcessLib/Output/Output.cpp b/ProcessLib/Output/Output.cpp
index 8249b314ba9..4da2b755d66 100644
--- a/ProcessLib/Output/Output.cpp
+++ b/ProcessLib/Output/Output.cpp
@@ -284,7 +284,7 @@ void Output::doOutputAlways(Process const& process,
     bool output_secondary_variable = true;
     // Need to add variables of process to vtu even no output takes place.
     addProcessDataToMesh(t, x, process_id, process.getMesh(), dof_tables,
-                         process.getProcessVariables(process_id),
+                         dof_tables, process.getProcessVariables(process_id),
                          process.getSecondaryVariables(),
                          output_secondary_variable,
                          process.getIntegrationPointWriter(process.getMesh()),
@@ -339,16 +339,16 @@ void Output::doOutputAlways(Process const& process,
             continue;
         }
         // mesh related output
-        auto& mesh = *BaseLib::findElementOrError(
+        auto& non_bulk_mesh = *BaseLib::findElementOrError(
             begin(_meshes), end(_meshes),
             [&mesh_output_name](auto const& m) {
                 return m->getName() == mesh_output_name;
             },
             "Need mesh '" + mesh_output_name + "' for the output.");
 
-        std::vector<MeshLib::Node*> const& nodes = mesh.getNodes();
+        std::vector<MeshLib::Node*> const& nodes = non_bulk_mesh.getNodes();
         DBUG("Found {:d} nodes for output at mesh '{:s}'.", nodes.size(),
-             mesh.getName());
+             non_bulk_mesh.getName());
 
         std::vector<std::unique_ptr<NumLib::LocalToGlobalIndexMap>>
             mesh_dof_tables;
@@ -357,7 +357,7 @@ void Output::doOutputAlways(Process const& process,
         {
             mesh_dof_tables.push_back(
                 process.getDOFTable(i).deriveBoundaryConstrainedMap(
-                    MeshLib::MeshSubset{mesh, nodes}));
+                    MeshLib::MeshSubset{non_bulk_mesh, nodes}));
         }
         std::vector<NumLib::LocalToGlobalIndexMap const*>
             mesh_dof_table_pointers;
@@ -369,18 +369,18 @@ void Output::doOutputAlways(Process const& process,
                   });
 
         output_secondary_variable = false;
-        addProcessDataToMesh(t, x, process_id, mesh, mesh_dof_table_pointers,
-                             process.getProcessVariables(process_id),
-                             process.getSecondaryVariables(),
-                             output_secondary_variable,
-                             process.getIntegrationPointWriter(mesh),
-                             _output_data_specification);
+        addProcessDataToMesh(
+            t, x, process_id, non_bulk_mesh, dof_tables,
+            mesh_dof_table_pointers, process.getProcessVariables(process_id),
+            process.getSecondaryVariables(), output_secondary_variable,
+            process.getIntegrationPointWriter(non_bulk_mesh),
+            _output_data_specification);
 
         // TODO (TomFischer): add pvd support here. This can be done if the
         // output is mesh related instead of process related. This would also
         // allow for merging bulk mesh output and arbitrary mesh output.
 
-        output_bulk_mesh(mesh);
+        output_bulk_mesh(non_bulk_mesh);
     }
     INFO("[time] Output of timestep {:d} took {:g} s.", timestep,
          time_output.elapsed());
@@ -444,7 +444,7 @@ void Output::doOutputNonlinearIteration(Process const& process,
 
     bool const output_secondary_variable = true;
     addProcessDataToMesh(t, x, process_id, process.getMesh(), dof_tables,
-                         process.getProcessVariables(process_id),
+                         dof_tables, process.getProcessVariables(process_id),
                          process.getSecondaryVariables(),
                          output_secondary_variable,
                          process.getIntegrationPointWriter(process.getMesh()),
diff --git a/ProcessLib/Output/ProcessOutput.cpp b/ProcessLib/Output/ProcessOutput.cpp
index 1a82f4e16ba..ae800fe77d8 100644
--- a/ProcessLib/Output/ProcessOutput.cpp
+++ b/ProcessLib/Output/ProcessOutput.cpp
@@ -19,6 +19,9 @@
 #include "InfoLib/GitInfo.h"
 #include "IntegrationPointWriter.h"
 #include "MathLib/LinAlg/LinAlg.h"
+#ifdef USE_PETSC
+#include "MeshLib/NodePartitionedMesh.h"
+#endif
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
 #include "MeshLib/IO/XDMF/XdmfHdfWriter.h"
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
@@ -136,6 +139,8 @@ namespace ProcessLib
 void addProcessDataToMesh(
     const double t, std::vector<GlobalVector*> const& x, int const process_id,
     MeshLib::Mesh& mesh,
+    [[maybe_unused]] std::vector<NumLib::LocalToGlobalIndexMap const*> const&
+        bulk_dof_table,
     std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
     std::vector<std::reference_wrapper<ProcessVariable>> const&
         process_variables,
@@ -207,6 +212,48 @@ void addProcessDataToMesh(
             auto const mesh_id = mesh_subset.getMeshID();
             for (auto const* node : mesh_subset.getNodes())
             {
+#ifdef USE_PETSC
+                if (bulk_dof_table[process_id] != dof_table[process_id])
+                {
+                    if (!mesh.getProperties().existsPropertyVector<std::size_t>(
+                            "bulk_node_ids"))
+                    {
+                        OGS_FATAL(
+                            "The required bulk node ids map does not exist in "
+                            "the boundary mesh '{:s}' or has the wrong data "
+                            "type (should be equivalent to C++ data type "
+                            "std::size_t which is an unsigned integer of size "
+                            "{:d} or UInt64 in vtk terminology).",
+                            mesh.getName(), sizeof(std::size_t));
+                    }
+                    auto const bulk_node_id_map =
+                        *mesh.getProperties().getPropertyVector<std::size_t>(
+                            "bulk_node_ids");
+
+                    if (static_cast<MeshLib::NodePartitionedMesh const&>(mesh)
+                            .isGhostNode(node->getID()))
+                    {
+                        auto const bulk_node_id =
+                            bulk_node_id_map[node->getID()];
+                        // use bulk_dof_table to find information
+                        std::size_t const bulk_mesh_id = 0;
+                        MeshLib::Location const l(bulk_mesh_id,
+                                                  MeshLib::MeshItemType::Node,
+                                                  bulk_node_id);
+                        auto const global_component_id =
+                            global_component_offset + component_id;
+                        auto const index =
+                            bulk_dof_table[process_id]->getLocalIndex(
+                                l, global_component_id,
+                                x[process_id]->getRangeBegin(),
+                                x[process_id]->getRangeEnd());
+
+                        output_data[node->getID() * n_components +
+                                    component_id] = x_copy[index];
+                        continue;
+                    }
+                }
+#endif
                 MeshLib::Location const l(mesh_id, MeshLib::MeshItemType::Node,
                                           node->getID());
 
diff --git a/ProcessLib/Output/ProcessOutput.h b/ProcessLib/Output/ProcessOutput.h
index 89b4cc631d1..028451aff6e 100644
--- a/ProcessLib/Output/ProcessOutput.h
+++ b/ProcessLib/Output/ProcessOutput.h
@@ -33,7 +33,8 @@ struct OutputDataSpecification final
 void addProcessDataToMesh(
     const double t, std::vector<GlobalVector*> const& x, int const process_id,
     MeshLib::Mesh& mesh,
-    std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+    std::vector<NumLib::LocalToGlobalIndexMap const*> const& bulk_dof_tables,
+    std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
     std::vector<std::reference_wrapper<ProcessVariable>> const&
         process_variables,
     SecondaryVariableCollection const& secondary_variables,
-- 
GitLab