From 0e5c9184df34048c94e2d329ceb2aba27c7ac155 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Tue, 23 Apr 2024 12:20:08 +0200
Subject: [PATCH] [PL] Added generic cell average output functionalty
---
ProcessLib/Output/CellAverageData.cpp | 59 ++++++++++++++++++
ProcessLib/Output/CellAverageData.h | 87 +++++++++++++++++++++++++++
2 files changed, 146 insertions(+)
create mode 100644 ProcessLib/Output/CellAverageData.cpp
create mode 100644 ProcessLib/Output/CellAverageData.h
diff --git a/ProcessLib/Output/CellAverageData.cpp b/ProcessLib/Output/CellAverageData.cpp
new file mode 100644
index 00000000000..1e46fa1a4cf
--- /dev/null
+++ b/ProcessLib/Output/CellAverageData.cpp
@@ -0,0 +1,59 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CellAverageData.h"
+
+#include "MeshLib/Utils/getOrCreateMeshProperty.h"
+
+namespace ProcessLib
+{
+MeshLib::PropertyVector<double>& CellAverageData::getOrCreatePropertyVector(
+ std::string const& name, unsigned const num_comp)
+{
+ if (auto const it = cell_averages_.find(name); it != cell_averages_.end())
+ {
+ auto& prop_vec = *it->second;
+ auto const num_comp_mesh = prop_vec.getNumberOfGlobalComponents();
+ if (num_comp_mesh == static_cast<int>(num_comp))
+ {
+ return prop_vec;
+ }
+
+ OGS_FATAL(
+ "The requested property '{}' has {} components, but the one "
+ "present in the mesh has {} components.",
+ name, num_comp, num_comp_mesh);
+ }
+
+ auto const name_in_mesh = name + "_avg";
+ auto [it, emplaced] = cell_averages_.emplace(
+ name, MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh_), name_in_mesh,
+ MeshLib::MeshItemType::Cell, num_comp));
+
+ if (!it->second)
+ {
+ OGS_FATAL("The cell property '{}' could not be added to the mesh.",
+ name_in_mesh);
+ }
+
+ if (!emplaced)
+ {
+ OGS_FATAL(
+ "Internal logic error. Something very bad happened. The cell "
+ "property '{}' was not added to the list of cell averages to "
+ "compute. There is some very strange inconsistency in the "
+ "code. Trouble ahead!",
+ name_in_mesh);
+ }
+
+ return *it->second;
+}
+} // namespace ProcessLib
diff --git a/ProcessLib/Output/CellAverageData.h b/ProcessLib/Output/CellAverageData.h
new file mode 100644
index 00000000000..3f8be132048
--- /dev/null
+++ b/ProcessLib/Output/CellAverageData.h
@@ -0,0 +1,87 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "MeshLib/Mesh.h"
+#include "MeshLib/PropertyVector.h"
+#include "ProcessLib/Reflection/ReflectionIPData.h"
+
+namespace ProcessLib
+{
+struct CellAverageData
+{
+ explicit CellAverageData(MeshLib::Mesh& mesh) : mesh_{mesh} {}
+
+ template <typename LAIntf>
+ void computeSecondaryVariable(
+ int const dim,
+ std::vector<std::unique_ptr<LAIntf>> const& local_assemblers)
+ {
+ auto const callback =
+ [this, &local_assemblers](std::string const& name,
+ unsigned const num_comp,
+ auto&& flattened_ip_data_accessor)
+ {
+ computeCellAverages(name, num_comp, flattened_ip_data_accessor,
+ local_assemblers);
+ };
+
+ if (dim == 2)
+ {
+ ProcessLib::Reflection::
+ forEachReflectedFlattenedIPDataAccessor<2, LAIntf>(
+ LAIntf::getReflectionDataForOutput(), callback);
+ }
+ else if (dim == 3)
+ {
+ ProcessLib::Reflection::
+ forEachReflectedFlattenedIPDataAccessor<3, LAIntf>(
+ LAIntf::getReflectionDataForOutput(), callback);
+ }
+ else
+ {
+ OGS_FATAL(
+ "Generic averaged output is only implemented for dimensions 2 "
+ "and 3 ATM. You requested dim = {}.",
+ dim);
+ }
+ }
+
+private:
+ MeshLib::PropertyVector<double>& getOrCreatePropertyVector(
+ std::string const& name, unsigned const num_comp);
+
+ void computeCellAverages(std::string const& name,
+ unsigned const num_comp,
+ auto&& flattened_ip_data_accessor,
+ auto const& local_assemblers)
+ {
+ auto& prop_vec = getOrCreatePropertyVector(name, num_comp);
+
+ for (std::size_t i = 0; i < local_assemblers.size(); ++i)
+ {
+ auto const& loc_asm = *local_assemblers[i];
+ auto const& ip_data = flattened_ip_data_accessor(loc_asm);
+ assert(ip_data.size() % num_comp == 0);
+ auto const num_ips =
+ static_cast<Eigen::Index>(ip_data.size() / num_comp);
+ Eigen::Map<const Eigen::MatrixXd> ip_data_mapped{ip_data.data(),
+ num_comp, num_ips};
+
+ Eigen::Map<Eigen::VectorXd>{&prop_vec[i * num_comp], num_comp} =
+ ip_data_mapped.rowwise().mean();
+ }
+ }
+
+ MeshLib::Mesh const& mesh_;
+ std::map<std::string, MeshLib::PropertyVector<double>*> cell_averages_;
+};
+} // namespace ProcessLib
--
GitLab