Add vtu output to reduced lung

Author: bgoderbauer

src/reduced_lung/4C_reduced_lung_helpers.cpp

@@ -176,6 +176,61 @@ namespace ReducedLung
 
     return column_map;
   }
+
+  void collect_runtime_output_data(
+      Core::IO::DiscretizationVisualizationWriterMesh& visualization_writer,
+      const std::vector<Airway>& airways, const std::vector<TerminalUnit>& terminal_units,
+      const Core::LinAlg::Vector<double>& dofs, const Core::LinAlg::Map* element_row_map)
+  {
+    Core::LinAlg::Vector<double> pressure_in(*element_row_map, true);
+    Core::LinAlg::Vector<double> pressure_out(*element_row_map, true);
+    Core::LinAlg::Vector<double> flow_in(*element_row_map, true);
+    for (const auto& airway : airways)
+    {
+      [[maybe_unused]] int err = pressure_in.replace_local_value(
+          airway.local_element_id, 0, dofs[airway.local_dof_ids[p_in]]);
+      FOUR_C_ASSERT(err == 0,
+          "Internal error: replace_local_value for runtime output (p_in from airways) did not "
+          "work.");
+      err = pressure_out.replace_local_value(
+          airway.local_element_id, 0, dofs[airway.local_dof_ids[p_out]]);
+      FOUR_C_ASSERT(err == 0,
+          "Internal error: replace_local_value for runtime output (p_out from airways) did not "
+          "work.");
+      err =
+          flow_in.replace_local_value(airway.local_element_id, 0, dofs[airway.local_dof_ids[q_in]]);
+      FOUR_C_ASSERT(err == 0,
+          "Internal error: replace_local_value for runtime output (q_in from airways) did not "
+          "work.");
+    }
+    for (const auto& terminal_unit : terminal_units)
+    {
+      [[maybe_unused]] int err = pressure_in.replace_local_value(
+          terminal_unit.local_element_id, 0, dofs[terminal_unit.local_dof_ids[p_in]]);
+      FOUR_C_ASSERT(err == 0,
+          "Internal error: replace_local_value for runtime output (p_in from terminal units) did "
+          "not "
+          "work.");
+      err = pressure_out.replace_local_value(
+          terminal_unit.local_element_id, 0, dofs[terminal_unit.local_dof_ids[p_out]]);
+      FOUR_C_ASSERT(err == 0,
+          "Internal error: replace_local_value for runtime output (p_out from terminal units) did "
+          "not "
+          "work.");
+      err = flow_in.replace_local_value(
+          terminal_unit.local_element_id, 0, dofs[terminal_unit.local_dof_ids[q_in]]);
+      FOUR_C_ASSERT(err == 0,
+          "Internal error: replace_local_value for runtime output (q_in from terminal units) did "
+          "not "
+          "work.");
+    }
+    visualization_writer.append_result_data_vector_with_context(
+        pressure_in, Core::IO::OutputEntity::element, {"p_1"});
+    visualization_writer.append_result_data_vector_with_context(
+        pressure_out, Core::IO::OutputEntity::element, {"p_2"});
+    visualization_writer.append_result_data_vector_with_context(
+        flow_in, Core::IO::OutputEntity::element, {"q_in"});
+  }
 }  // namespace ReducedLung
 
 FOUR_C_NAMESPACE_CLOSE

src/reduced_lung/4C_reduced_lung_helpers.hpp

@@ -12,6 +12,7 @@
 
 #include "4C_comm_mpi_utils.hpp"
 #include "4C_global_data.hpp"
+#include "4C_io_discretization_visualization_writer_mesh.hpp"
 #include "4C_linalg_map.hpp"
 
 FOUR_C_NAMESPACE_OPEN
@@ -98,8 +99,8 @@ namespace ReducedLung
 
   struct Airway
   {
-    int global_equation_id;
-    int local_equation_id;
+    int global_element_id;
+    int local_element_id;
     int local_airway_id;
     AirwayType airway_type;
     // dofs: {p1, p2, q} for resistive airways; {p1, p2, q1, q2} for compliant airways
@@ -117,8 +118,8 @@ namespace ReducedLung
   // Save rheological model and (non-)linear elasticity separately
   struct TerminalUnit
   {
-    int global_equation_id;
-    int local_equation_id;
+    int global_element_id;
+    int local_element_id;
     int local_terminal_unit_id;
     TerminalUnitType tu_type;
     double E;
@@ -203,6 +204,11 @@ namespace ReducedLung
       const std::map<int, int>& first_global_dof_of_ele, const std::vector<Connection>& connections,
       const std::vector<Bifurcation>& bifurcations,
       const std::vector<BoundaryCondition>& boundary_conditions);
+
+  void collect_runtime_output_data(
+      Core::IO::DiscretizationVisualizationWriterMesh& visualization_writer,
+      const std::vector<Airway>& airways, const std::vector<TerminalUnit>& terminal_units,
+      const Core::LinAlg::Vector<double>& dofs, const Core::LinAlg::Map* element_row_map);
 }  // namespace ReducedLung
 
 FOUR_C_NAMESPACE_CLOSE

src/reduced_lung/4C_reduced_lung_main.cpp

@@ -13,6 +13,7 @@
 #include "4C_fem_general_node.hpp"
 #include "4C_global_data.hpp"
 #include "4C_io.hpp"
+#include "4C_io_discretization_visualization_writer_mesh.hpp"
 #include "4C_linalg_utils_sparse_algebra_manipulation.hpp"
 #include "4C_linear_solver_method_linalg.hpp"
 #include "4C_mat_maxwell_0d_acinus_NeoHookean.hpp"
@@ -47,6 +48,12 @@ namespace ReducedLung
         Teuchos::getIntegralValue<Core::IO::Verbositylevel>(
             Global::Problem::instance()->io_params(), "VERBOSITY"));
     actdis->compute_null_space_if_necessary(solver->params());
+    // Create runtime output writer
+    const auto visualization_writer =
+        std::make_unique<Core::IO::DiscretizationVisualizationWriterMesh>(
+            actdis, Core::IO::visualization_parameters_factory(
+                        Global::Problem::instance()->io_params().sublist("RUNTIME VTK OUTPUT"),
+                        *Global::Problem::instance()->output_control_file(), 0));
     // The existing mpi communicator is recycled for the new data layout.
     const auto& comm = actdis->get_comm();
 
@@ -122,17 +129,17 @@ namespace ReducedLung
     // Assign every local element its associated global dof ids.
     for (auto& airway : airways)
     {
-      int first_dof_gid = first_global_dof_of_ele[airway.global_equation_id];
-      int n_dof = dof_per_ele[airway.global_equation_id];
+      int first_dof_gid = first_global_dof_of_ele[airway.global_element_id];
+      int n_dof = dof_per_ele[airway.global_element_id];
       for (int i = 0; i < n_dof; i++)
       {
         airway.global_dof_ids.insert(airway.global_dof_ids.end(), first_dof_gid + i);
       }
     }
     for (auto& terminal_unit : terminal_units)
     {
-      int first_dof_gid = first_global_dof_of_ele[terminal_unit.global_equation_id];
-      int n_dof = dof_per_ele[terminal_unit.global_equation_id];
+      int first_dof_gid = first_global_dof_of_ele[terminal_unit.global_element_id];
+      int n_dof = dof_per_ele[terminal_unit.global_element_id];
       for (int i = 0; i < n_dof; i++)
       {
         terminal_unit.global_dof_ids.insert(terminal_unit.global_dof_ids.end(), first_dof_gid + i);
@@ -434,7 +441,7 @@ namespace ReducedLung
     {
       const int airway_id = airway.local_airway_id;
       const auto* airway_ele =
-          static_cast<Discret::Elements::RedAirway*>(actdis->g_element(airway.global_equation_id));
+          static_cast<Discret::Elements::RedAirway*>(actdis->g_element(airway.global_element_id));
       const Discret::ReducedLung::AirwayParams airway_params = airway_ele->get_airway_params();
       const std::vector<double>& coords_1 = airway_ele->nodes()[0]->x();
       const std::vector<double>& coords_2 = airway_ele->nodes()[1]->x();
@@ -457,7 +464,7 @@ namespace ReducedLung
     {
       const int terminal_unit_id = terminal_unit.local_terminal_unit_id;
       const auto* terminal_unit_ele = static_cast<Discret::Elements::RedAcinus*>(
-          actdis->g_element(terminal_unit.global_equation_id));
+          actdis->g_element(terminal_unit.global_element_id));
       const std::vector<double>& coords_1 = terminal_unit_ele->nodes()[0]->x();
       const std::vector<double>& coords_2 = terminal_unit_ele->nodes()[1]->x();
       const double radius = std::sqrt((coords_1[0] - coords_2[0]) * (coords_1[0] - coords_2[0]) +
@@ -484,6 +491,7 @@ namespace ReducedLung
     auto sysmat = Epetra_CrsMatrix(
         Copy, row_map.get_epetra_map(), locally_relevant_dof_map.get_epetra_map(), 3);
 
+    const int results_every = rawdyn.get<int>("RESULTSEVERY");
     // Time integration parameters.
     const double dt = rawdyn.get<double>("TIMESTEP");
     const int n_timesteps = rawdyn.get<int>("NUMSTEP");
@@ -526,15 +534,15 @@ namespace ReducedLung
         if (!sysmat.Filled())
         {
           err = sysmat.InsertMyValues(
-              airway.local_equation_id, vals.size(), vals.data(), airway.local_dof_ids.data());
+              airway.local_element_id, vals.size(), vals.data(), airway.local_dof_ids.data());
         }
         else
         {
           err = sysmat.ReplaceMyValues(
-              airway.local_equation_id, vals.size(), vals.data(), airway.local_dof_ids.data());
+              airway.local_element_id, vals.size(), vals.data(), airway.local_dof_ids.data());
         }
         FOUR_C_ASSERT(err == 0, "Internal error: Airway equation assembly did not work.");
-        err = rhs.replace_local_value(airway.local_equation_id, 0, res);
+        err = rhs.replace_local_value(airway.local_element_id, 0, res);
         FOUR_C_ASSERT(err == 0, "Internal error: Airway equation calculation did not work.");
       }
 
@@ -554,16 +562,16 @@ namespace ReducedLung
             E * (V_tu - V0_tu) / V0_tu;
         if (!sysmat.Filled())
         {
-          err = sysmat.InsertMyValues(terminal_unit.local_equation_id, vals.size(), vals.data(),
+          err = sysmat.InsertMyValues(terminal_unit.local_element_id, vals.size(), vals.data(),
               terminal_unit.local_dof_ids.data());
         }
         else
         {
-          err = sysmat.ReplaceMyValues(terminal_unit.local_equation_id, vals.size(), vals.data(),
+          err = sysmat.ReplaceMyValues(terminal_unit.local_element_id, vals.size(), vals.data(),
               terminal_unit.local_dof_ids.data());
         }
         FOUR_C_ASSERT(err == 0, "Internal error: Terminal Unit equation assembly did not work.");
-        err = rhs.replace_local_value(terminal_unit.local_equation_id, 0, res);
+        err = rhs.replace_local_value(terminal_unit.local_element_id, 0, res);
         FOUR_C_ASSERT(err == 0, "Internal error: Terminal Unit equation calculation did not work.");
       }
 
@@ -735,6 +743,15 @@ namespace ReducedLung
         // Backwards Euler: V_n+1 = V_n + q_n+1 * dt.
         volume[tu_id] += locally_relevant_dofs[q_id] * dt;
       }
+
+      // Runtime output
+      if (n % results_every == 0)
+      {
+        visualization_writer->reset();
+        collect_runtime_output_data(
+            *visualization_writer, airways, terminal_units, dofs, actdis->element_row_map());
+        visualization_writer->write_to_disk(dt * n, n);
+      }
     }
     // Print time monitor
     Teuchos::TimeMonitor::summarize(std::cout, false, true, false);