[PL] multi-component secondary variable output

ProcessLib/ProcessOutput.cpp

@@ -41,7 +41,6 @@ ProcessOutput::ProcessOutput(BaseLib::ConfigTree const& output_config)
 }
 
 void doProcessOutput(std::string const& file_name,
-                     bool const compress_output,
                      GlobalVector const& x,
                      MeshLib::Mesh& mesh,
                      NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -119,6 +118,45 @@ void doProcessOutput(std::string const& file_name,
 #ifndef USE_PETSC
     // the following section is for the output of secondary variables
 
+    auto count_mesh_items = [](MeshLib::Mesh const& mesh,
+                               MeshLib::MeshItemType type) -> std::size_t {
+        switch (type)
+        {
+            case MeshLib::MeshItemType::Cell:
+                return mesh.getNumberOfElements();
+            case MeshLib::MeshItemType::Node:
+                return mesh.getNumberOfNodes();
+            default:
+                break;  // avoid compiler warning
+        }
+        return 0;
+    };
+
+    auto get_or_create_mesh_property = [&mesh, &count_mesh_items](
+        std::string const& property_name, MeshLib::MeshItemType type,
+        unsigned const num_components) {
+        // Get or create a property vector for results.
+        MeshLib::PropertyVector<double>* result = nullptr;
+
+        auto const N = count_mesh_items(mesh, type) * num_components;
+
+        if (mesh.getProperties().hasPropertyVector(property_name))
+        {
+            result = mesh.getProperties().template getPropertyVector<double>(
+                property_name);
+        }
+        else
+        {
+            result =
+                mesh.getProperties().template createNewPropertyVector<double>(
+                    property_name, type, num_components);
+            result->resize(N);
+        }
+        assert(result && result->size() == N);
+
+        return result;
+    };
+
     auto add_secondary_var = [&](SecondaryVariable const& var,
                                  std::string const& output_name) {
         assert(var.fcts.num_components == 1);  // TODO implement other cases
@@ -126,18 +164,22 @@ void doProcessOutput(std::string const& file_name,
         {
             DBUG("  secondary variable %s", output_name.c_str());
 
-            auto result = MeshLib::getOrCreateMeshProperty<double>(
-                mesh, output_name, MeshLib::MeshItemType::Node, 1);
+            auto& result = *get_or_create_mesh_property(
+                output_name, MeshLib::MeshItemType::Node,
+                var.fcts.num_components);
+            assert(result.size() ==
+                   mesh.getNumberOfNodes() * var.fcts.num_components);
 
             std::unique_ptr<GlobalVector> result_cache;
             auto const& nodal_values =
                     var.fcts.eval_field(x, dof_table, result_cache);
+            assert(nodal_values.size() == result.size());
 
             // Copy result
             for (GlobalIndexType i = 0; i < nodal_values.size(); ++i)
             {
                 assert(!std::isnan(nodal_values[i]));
-                (*result)[i] = nodal_values[i];
+                result[i] = nodal_values[i];
             }
         }
 
@@ -146,18 +188,23 @@ void doProcessOutput(std::string const& file_name,
             DBUG("  secondary variable %s residual", output_name.c_str());
             auto const& property_name_res = output_name + "_residual";
 
-            auto result = MeshLib::getOrCreateMeshProperty<double>(
-                mesh, property_name_res, MeshLib::MeshItemType::Cell, 1);
+            auto& result = *get_or_create_mesh_property(
+                property_name_res, MeshLib::MeshItemType::Cell,
+                var.fcts.num_components);
+            assert(result.size() ==
+                   mesh.getNumberOfElements() * var.fcts.num_components);
 
             std::unique_ptr<GlobalVector> result_cache;
             auto const& residuals =
                     var.fcts.eval_residuals(x, dof_table, result_cache);
 
+            assert(residuals.size() == result.size());
+
             // Copy result
             for (GlobalIndexType i = 0; i < residuals.size(); ++i)
             {
                 assert(!std::isnan(residuals[i]));
-                (*result)[i] = residuals[i];
+                result[i] = residuals[i];
             }
         }
     };