diff --git a/MeshLib/Utils/IntegrationPointWriter.h b/MeshLib/Utils/IntegrationPointWriter.h
index 10df6f5fd86ad2f617c0f4b85c663fae9f10263a..e4eae288c6240074c0bfe4bcac56acc41ad05fa0 100644
--- a/MeshLib/Utils/IntegrationPointWriter.h
+++ b/MeshLib/Utils/IntegrationPointWriter.h
@@ -22,6 +22,10 @@ class Properties;
struct IntegrationPointWriter final
{
+ /// Constructor taking a member function of the \c LocalAssemblerInterface
+ /// as a callback.
+ ///
+ /// The callback \c getIpData must return a \c std::vector<double>.
template <typename LocalAssemblerInterface, typename... Args>
IntegrationPointWriter(
std::string const& name,
@@ -53,6 +57,40 @@ struct IntegrationPointWriter final
return result;
};
}
+
+ /// Constructor taking an \c accessor function as a callback.
+ ///
+ /// The accessor function must take a const reference to a
+ /// \c LocalAssemblerInterface object as its only argument.
+ /// The callback \c accessor must return a \c std::vector<double>.
+ template <typename LocalAssemblerInterface, typename Accessor>
+ IntegrationPointWriter(
+ std::string const& name,
+ int const n_components,
+ int const integration_order,
+ std::vector<std::unique_ptr<LocalAssemblerInterface>> const&
+ local_assemblers,
+ Accessor accessor)
+ : _name(name),
+ _n_components(n_components),
+ _integration_order(integration_order)
+ {
+ _callback = [&local_assemblers, accessor]
+ {
+ // Result containing integration point data for each local
+ // assembler.
+ std::vector<std::vector<double>> result;
+ result.reserve(local_assemblers.size());
+
+ std::transform(begin(local_assemblers), end(local_assemblers),
+ std::back_inserter(result),
+ [&accessor](auto const& la)
+ { return accessor(*la); });
+
+ return result;
+ };
+ }
+
int numberOfComponents() const { return _n_components; }
int integrationOrder() const { return _integration_order; }
std::string name() const { return _name; }
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index 97f952ea02e1f591b11cc20e321cd446d19461d1..100d9684a49782916e4cb5133389d14b6fed0e59 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -13,6 +13,7 @@ append_source_files(SOURCES BoundaryConditionAndSourceTerm/Utils)
append_source_files(SOURCES SurfaceFlux)
append_source_files(SOURCES Output)
append_source_files(SOURCES Utils)
+append_source_files(SOURCES Reflection)
ogs_add_library(ProcessLib GENERATE_EXPORT_HEADER ${SOURCES})
@@ -29,6 +30,7 @@ target_link_libraries(
MeshGeoToolsLib
MeshLib
NumLib
+ Boost::mp11
range-v3
$<$<TARGET_EXISTS:ProcessLibBoundaryConditionAndSourceTermPython>:ProcessLibBoundaryConditionAndSourceTermPython>
$<$<TARGET_EXISTS:petsc>:petsc>
diff --git a/ProcessLib/Output/SecondaryVariable.h b/ProcessLib/Output/SecondaryVariable.h
index 3448165720fa4c914210f2b0529643a2894e9772..25b21e199381aa69849628bbb7100b8fa8783c0f 100644
--- a/ProcessLib/Output/SecondaryVariable.h
+++ b/ProcessLib/Output/SecondaryVariable.h
@@ -13,6 +13,7 @@
#include "BaseLib/Algorithm.h"
#include "NumLib/Extrapolation/ExtrapolatableElementCollection.h"
#include "NumLib/Extrapolation/Extrapolator.h"
+#include "ProcessLib/Utils/TransposeInPlace.h"
namespace NumLib
{
@@ -174,8 +175,9 @@ SecondaryVariableFunctions makeExtrapolator(
const double t,
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::unique_ptr<GlobalVector> & /*result_cache*/
- ) -> GlobalVector const& {
+ std::unique_ptr<GlobalVector>& /*result_cache*/
+ ) -> GlobalVector const&
+ {
auto const extrapolatables = NumLib::makeExtrapolatable(
local_assemblers, integration_point_values_method);
extrapolator.extrapolate(num_components, extrapolatables, t, x,
@@ -189,8 +191,9 @@ SecondaryVariableFunctions makeExtrapolator(
const double t,
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::unique_ptr<GlobalVector> & /*result_cache*/
- ) -> GlobalVector const& {
+ std::unique_ptr<GlobalVector>& /*result_cache*/
+ ) -> GlobalVector const&
+ {
auto const extrapolatables = NumLib::makeExtrapolatable(
local_assemblers, integration_point_values_method);
extrapolator.calculateResiduals(num_components, extrapolatables, t, x,
@@ -200,4 +203,57 @@ SecondaryVariableFunctions makeExtrapolator(
return {num_components, eval_field, eval_residuals};
}
+//! A variant that takes an \c accessor function as a callback.
+//!
+//! The \c accessor must take a local assembler by const reference as its only
+//! argument and must return a <tt>std::vector\<double\></tt> of integration
+//! point data. The data returned by the accessor must be in "integration point
+//! writer order", which is transposed compared to "extrapolator order".
+template <typename LocalAssemblerCollection, typename IPDataAccessor>
+SecondaryVariableFunctions makeExtrapolator2(
+ const unsigned num_components,
+ NumLib::Extrapolator& extrapolator,
+ LocalAssemblerCollection const& local_assemblers,
+ IPDataAccessor&& accessor)
+{
+ using LocalAssemblerInterface = std::remove_cvref_t<
+ decltype(*std::declval<LocalAssemblerCollection>()[0])>;
+ static_assert(std::is_invocable_r_v<std::vector<double>, IPDataAccessor,
+ LocalAssemblerInterface const&>);
+
+ if (num_components == 1)
+ {
+ auto method_wrapped =
+ [accessor](
+ LocalAssemblerInterface const& loc_asm, const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<
+ NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) -> std::vector<double> const&
+ {
+ cache = accessor(loc_asm);
+ return cache;
+ };
+
+ return makeExtrapolator(num_components, extrapolator, local_assemblers,
+ method_wrapped);
+ }
+
+ auto method_wrapped =
+ [accessor, num_components](
+ LocalAssemblerInterface const& loc_asm, const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<
+ NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) -> std::vector<double> const&
+ {
+ cache = accessor(loc_asm);
+ transposeInPlace(cache, cache.size() / num_components);
+ return cache;
+ };
+
+ return makeExtrapolator(num_components, extrapolator, local_assemblers,
+ method_wrapped);
+}
+
} // namespace ProcessLib
diff --git a/ProcessLib/Reflection/ReflectionData.h b/ProcessLib/Reflection/ReflectionData.h
new file mode 100644
index 0000000000000000000000000000000000000000..eecdc2b39d533e62cf22a557dc64d13115b2c2e8
--- /dev/null
+++ b/ProcessLib/Reflection/ReflectionData.h
@@ -0,0 +1,44 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, 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 <string>
+#include <tuple>
+
+namespace ProcessLib::Reflection
+{
+template <typename Class, typename Member>
+struct ReflectionData
+{
+ ReflectionData(std::string name, Member Class::*field)
+ : name(std::move(name)), field(field)
+ {
+ }
+
+ explicit ReflectionData(Member Class::*field) : field(field) {}
+
+ std::string name;
+ Member Class::*field;
+};
+
+template <typename Class, typename Member>
+std::tuple<ReflectionData<Class, Member>> reflectWithName(
+ std::string const& name, Member Class::*field)
+{
+ return {{name, field}};
+}
+
+template <typename Class, typename... Members>
+auto reflectWithoutName(Members Class::*... members)
+{
+ return std::tuple{ReflectionData<Class, Members>{"", members}...};
+}
+} // namespace ProcessLib::Reflection
diff --git a/ProcessLib/Reflection/ReflectionForExtrapolation.h b/ProcessLib/Reflection/ReflectionForExtrapolation.h
new file mode 100644
index 0000000000000000000000000000000000000000..a8dac205cd1db2329791c999a86b78d2eb062ad3
--- /dev/null
+++ b/ProcessLib/Reflection/ReflectionForExtrapolation.h
@@ -0,0 +1,40 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, 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 "ProcessLib/Output/SecondaryVariable.h"
+#include "ReflectionIPData.h"
+
+namespace ProcessLib::Reflection
+{
+/// Adds secondary variables for all IP data obtained recursively from the given
+/// \c reflection_data to the given secondary variable collection.
+template <int Dim, typename LocAsmIF, typename ReflData>
+void addReflectedSecondaryVariables(
+ ReflData const& reflection_data,
+ SecondaryVariableCollection& secondary_variables,
+ NumLib::Extrapolator& extrapolator,
+ std::vector<std::unique_ptr<LocAsmIF>> const& local_assemblers)
+{
+ forEachReflectedFlattenedIPDataAccessor<Dim, LocAsmIF>(
+ reflection_data,
+ [&secondary_variables, &local_assemblers, &extrapolator](
+ std::string const& name,
+ unsigned const num_comp,
+ auto&& flattened_ip_data_accessor)
+ {
+ secondary_variables.addSecondaryVariable(
+ name,
+ makeExtrapolator2(num_comp, extrapolator, local_assemblers,
+ flattened_ip_data_accessor));
+ });
+}
+} // namespace ProcessLib::Reflection
diff --git a/ProcessLib/Reflection/ReflectionForIPWriters.h b/ProcessLib/Reflection/ReflectionForIPWriters.h
new file mode 100644
index 0000000000000000000000000000000000000000..dc7a202b2a09ec849630040cc09bd39d1ed24bee
--- /dev/null
+++ b/ProcessLib/Reflection/ReflectionForIPWriters.h
@@ -0,0 +1,42 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, 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/Utils/IntegrationPointWriter.h"
+#include "ReflectionIPData.h"
+
+namespace ProcessLib::Reflection
+{
+/// Adds IP data writers for all IP data obtained recursively from the given
+/// \c reflection_data to the given IP writer vector.
+template <int Dim, typename LocAsmIF, typename ReflData>
+void addReflectedIntegrationPointWriters(
+ ReflData const& reflection_data,
+ std::vector<std::unique_ptr<MeshLib::IntegrationPointWriter>>&
+ integration_point_writers,
+ unsigned const integration_order,
+ std::vector<std::unique_ptr<LocAsmIF>> const& local_assemblers)
+{
+ forEachReflectedFlattenedIPDataAccessor<Dim, LocAsmIF>(
+ reflection_data,
+ [&integration_point_writers, integration_order, &local_assemblers](
+ std::string const& name,
+ unsigned const num_comp,
+ auto&& flattened_ip_data_accessor)
+ {
+ // TODO check if writer with such a name already exists.
+ integration_point_writers.emplace_back(
+ std::make_unique<MeshLib::IntegrationPointWriter>(
+ name + "_ip", num_comp, integration_order, local_assemblers,
+ flattened_ip_data_accessor));
+ });
+}
+} // namespace ProcessLib::Reflection
diff --git a/ProcessLib/Reflection/ReflectionIPData.h b/ProcessLib/Reflection/ReflectionIPData.h
new file mode 100644
index 0000000000000000000000000000000000000000..a71ce46fd9ed542691b22d1b29cfd79d59d9a6bb
--- /dev/null
+++ b/ProcessLib/Reflection/ReflectionIPData.h
@@ -0,0 +1,310 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, 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 <boost/mp11.hpp>
+
+#include "MathLib/KelvinVector.h"
+#include "ReflectionData.h"
+
+namespace ProcessLib::Reflection
+{
+namespace detail
+{
+// Used in metaprogramming to check if the type T has a member "reflect".
+template <typename T>
+constexpr bool has_reflect = requires
+{
+ T::reflect();
+};
+
+template <typename T>
+struct NumberOfComponents;
+
+template <>
+struct NumberOfComponents<double> : std::integral_constant<unsigned, 1>
+{
+};
+
+template <int N>
+struct NumberOfComponents<Eigen::Matrix<double, N, 1, Eigen::ColMajor, N, 1>>
+ : std::integral_constant<unsigned, N>
+{
+};
+
+/** A function object taking a local assembler as its argument and returning a
+ * <tt>std::vector\<double\></tt> of some specific "flattened" integration point
+ * (IP) data.
+ *
+ * \tparam Dim the space dimension
+ * \tparam Accessor_IPDataVecInLocAsm see below
+ * \tparam Accessor_CurrentLevelFromIPDataVecElement see below
+ *
+ * In OGS IP data is usually stored in the local assembler in the following way:
+ *
+ * \code
+ * struct LocAsm
+ * {
+ * std::vector<IPData1> ip_data1;
+ * std::vector<IPData2> ip_data2;
+ * };
+ * \endcode
+ *
+ * The types \c IPData1 and \c IPData2 might directly contain the IP data or
+ * might have \c struct members who contain the IP data, e.g.:
+ *
+ * \code
+ * struct IPData1
+ * {
+ * double scalar1; // IPData1 directly holds IP data
+ * Eigen::Vector<double, 3> vector1;
+ * };
+ *
+ * struct IPData_Level2
+ * {
+ * double scalar2;
+ * Eigen::Vector<double, 6> kelvin2;
+ * };
+ *
+ * struct IPData2
+ * {
+ * IPData_Level2 level2; // IPData2 does not directly hold IP data
+ * };
+ * \endcode
+ *
+ * \c Accessor_IPDataVecInLocAsm is a function object with signature
+ * <tt>LocAsm const\& -\> std::vector\<IPData\> const\&</tt>.
+ *
+ * \c Accessor_CurrentLevelFromIPDataVecElement is a function object with
+ * signature <tt>IPData const\& -\> double (or Eigen::Vector)</tt>,
+ * where \c IPData is the "top level" struct contained in the
+ * <tt>std::vector\<IPData\></tt>.
+ *
+ * I.e. the first accessor takes us from the local assembler to the IP data
+ * vector and the second accessor takes us from an IP data vector element to the
+ * final IP data of type \c double or <tt>Eigen::Vector</tt>.
+ *
+ * \note This function object transforms Kelvin vector typed IP data to a
+ * ParaView compatible symmetric tensor representation.
+ */
+template <int Dim, typename Accessor_IPDataVecInLocAsm,
+ typename Accessor_CurrentLevelFromIPDataVecElement>
+struct GetFlattenedIPDataFromLocAsm
+{
+ static_assert(!std::is_reference_v<Accessor_IPDataVecInLocAsm>);
+ static_assert(
+ !std::is_reference_v<Accessor_CurrentLevelFromIPDataVecElement>);
+
+ Accessor_IPDataVecInLocAsm accessor_ip_data_vec_in_loc_asm;
+ Accessor_CurrentLevelFromIPDataVecElement
+ accessor_current_level_from_ip_data_vec_element;
+
+ template <typename LocAsm>
+ std::vector<double> operator()(LocAsm const& loc_asm) const
+ {
+ using IPDataVector = std::remove_cvref_t<
+ std::invoke_result_t<Accessor_IPDataVecInLocAsm, LocAsm const&>>;
+ using IPDataVectorElement = typename IPDataVector::value_type;
+
+ // the concrete IP data, e.g. double or Eigen::Vector
+ using ConcreteIPData = std::remove_cvref_t<
+ std::invoke_result_t<Accessor_CurrentLevelFromIPDataVecElement,
+ IPDataVectorElement const&>>;
+
+ constexpr unsigned num_comp = NumberOfComponents<ConcreteIPData>::value;
+ auto const& ip_data_vector = accessor_ip_data_vec_in_loc_asm(loc_asm);
+ auto const num_ips = ip_data_vector.size();
+
+ std::vector<double> result(num_comp * num_ips);
+
+ for (std::size_t ip = 0; ip < num_ips; ++ip)
+ {
+ auto const& ip_data_vector_element = ip_data_vector[ip];
+ auto const& ip_data =
+ accessor_current_level_from_ip_data_vec_element(
+ ip_data_vector_element);
+
+ if constexpr (num_comp == 1)
+ {
+ // TODO optimization if nothing needs to be copied?
+ result[ip] = ip_data;
+ }
+ else if constexpr (num_comp ==
+ MathLib::KelvinVector::kelvin_vector_dimensions(
+ Dim))
+ {
+ auto const converted =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
+ ip_data);
+
+ for (unsigned comp = 0; comp < num_comp; ++comp)
+ {
+ result[ip * num_comp + comp] = converted[comp];
+ }
+ }
+ else
+ {
+ for (unsigned comp = 0; comp < num_comp; ++comp)
+ {
+ result[ip * num_comp + comp] = ip_data[comp];
+ }
+ }
+ }
+ return result;
+ }
+};
+
+// Convenience function for template argument deduction with
+// GetFlattenedIPDataFromLocAsm
+template <int Dim, typename Accessor_IPDataVecInLocAsm,
+ typename Accessor_CurrentLevelFromIPDataVecElement>
+GetFlattenedIPDataFromLocAsm<
+ Dim, std::remove_cvref_t<Accessor_IPDataVecInLocAsm>,
+ std::remove_cvref_t<Accessor_CurrentLevelFromIPDataVecElement>>
+getFlattenedIPDataFromLocAsm(
+ Accessor_IPDataVecInLocAsm accessor_ip_data_vec_in_loc_asm,
+ Accessor_CurrentLevelFromIPDataVecElement
+ accessor_current_level_from_ip_data_vec_element)
+{
+ return {std::forward<Accessor_IPDataVecInLocAsm>(
+ accessor_ip_data_vec_in_loc_asm),
+ std::forward<Accessor_CurrentLevelFromIPDataVecElement>(
+ accessor_current_level_from_ip_data_vec_element)};
+}
+
+// Convenience function for template argument deduction with
+// GetFlattenedIPDataFromLocAsm. Overload of the function above with less
+// arguments.
+template <int Dim, typename Accessor_IPDataVecInLocAsm>
+auto getFlattenedIPDataFromLocAsm(
+ Accessor_IPDataVecInLocAsm&& accessor_ip_data_vec_in_loc_asm)
+{
+ return getFlattenedIPDataFromLocAsm<Dim>(
+ std::forward<Accessor_IPDataVecInLocAsm>(
+ accessor_ip_data_vec_in_loc_asm),
+ std::identity{});
+}
+
+/// Calls the given \c callback for each flattened IP data accessor obtained
+/// recursively from the given \c reflection data.
+///
+/// The \c callback function must take (i) the name of the IP data (\c
+/// std::string), (ii) their number of components (\c unsigned) and a flattened
+/// IP data accessor of type
+/// ProcessLib::Reflection::detail::GetFlattenedIPDataFromLocAsm as arguments.
+///
+/// \see ProcessLib::Reflection::detail::GetFlattenedIPDataFromLocAsm for
+/// details, also on the template parameters.
+template <int Dim, typename Callback, typename ReflectionDataTuple,
+ typename Accessor_IPDataVecInLocAsm,
+ typename Accessor_CurrentLevelFromIPDataVecElement>
+void forEachReflectedFlattenedIPDataAccessor(
+ Callback const& callback, ReflectionDataTuple const& reflection_data,
+ Accessor_IPDataVecInLocAsm const& accessor_ip_data_vec_in_loc_asm,
+ Accessor_CurrentLevelFromIPDataVecElement const&
+ accessor_current_level_from_ip_data_vec_element)
+{
+ boost::mp11::tuple_for_each(
+ reflection_data,
+ [&accessor_ip_data_vec_in_loc_asm,
+ &accessor_current_level_from_ip_data_vec_element,
+ &callback]<typename Class, typename Member>(
+ ReflectionData<Class, Member> const& refl_data)
+ {
+ auto accessor_member_from_ip_data_vec_element =
+ [level = refl_data.field,
+ accessor_current_level_from_ip_data_vec_element](
+ auto const& ip_data_vec_element) -> Member const&
+ {
+ return accessor_current_level_from_ip_data_vec_element(
+ ip_data_vec_element).*
+ level;
+ };
+
+ if constexpr (has_reflect<Member>)
+ {
+ forEachReflectedFlattenedIPDataAccessor<Dim>(
+ callback, Member::reflect(),
+ accessor_ip_data_vec_in_loc_asm,
+ accessor_member_from_ip_data_vec_element);
+ }
+ else
+ {
+ constexpr unsigned num_comp = NumberOfComponents<Member>::value;
+
+ callback(refl_data.name, num_comp,
+ getFlattenedIPDataFromLocAsm<Dim>(
+ accessor_ip_data_vec_in_loc_asm,
+ accessor_member_from_ip_data_vec_element));
+ }
+ });
+}
+
+// Overload of the function above with less arguments
+template <int Dim, typename Callback, typename ReflectionDataTuple,
+ typename Accessor_IPDataVecInLocAsm>
+void forEachReflectedFlattenedIPDataAccessor(
+ Callback const& callback,
+ ReflectionDataTuple const& reflection_data,
+ Accessor_IPDataVecInLocAsm const& accessor_ip_data_vec_in_loc_asm)
+{
+ forEachReflectedFlattenedIPDataAccessor<Dim>(
+ callback, reflection_data, accessor_ip_data_vec_in_loc_asm,
+ std::identity{});
+}
+
+} // namespace detail
+
+/// Calls the passed \c callback for each IP data accessor for the given
+/// \c LocAsmIF class.
+///
+/// The IP data accessors are obtained via reflection (i.e., via the static
+/// \c reflect() method.
+///
+/// The IP data accessors provide IP data as a flat \c std::vector<double>.
+///
+/// The \c callback must accept name, number of components and a function object
+/// with signature <tt>LocAsmIF const\& -\> std::vector\<double\></tt> as its
+/// arguments.
+template <int Dim, typename LocAsmIF, typename Callback, typename ReflData>
+void forEachReflectedFlattenedIPDataAccessor(ReflData const& reflection_data,
+ Callback const& callback)
+{
+ boost::mp11::tuple_for_each(
+ reflection_data,
+ [&callback]<typename Class, typename Member>(
+ ReflectionData<Class, std::vector<Member>> const& refl_data)
+ {
+ auto accessor_ip_data_vec_in_loc_asm = [ip_data_vector =
+ refl_data.field](
+ LocAsmIF const& loc_asm) -> auto const&
+ {
+ return loc_asm.*ip_data_vector;
+ };
+
+ if constexpr (detail::has_reflect<Member>)
+ {
+ detail::forEachReflectedFlattenedIPDataAccessor<Dim>(
+ callback, Member::reflect(),
+ accessor_ip_data_vec_in_loc_asm);
+ }
+ else
+ {
+ constexpr unsigned num_comp =
+ detail::NumberOfComponents<Member>::value;
+
+ callback(refl_data.name, num_comp,
+ detail::getFlattenedIPDataFromLocAsm<Dim>(
+ accessor_ip_data_vec_in_loc_asm));
+ }
+ });
+}
+} // namespace ProcessLib::Reflection
diff --git a/ProcessLib/ThermoRichardsMechanics/Constitutive/Base.h b/ProcessLib/ThermoRichardsMechanics/Constitutive/Base.h
index 60d229558fe9fca170c7be472a9d23b88bcc394a..10a8ffc179aca6860aab1e6ffa47e28be079ffda 100644
--- a/ProcessLib/ThermoRichardsMechanics/Constitutive/Base.h
+++ b/ProcessLib/ThermoRichardsMechanics/Constitutive/Base.h
@@ -13,6 +13,7 @@
#include "MaterialLib/MPL/Medium.h"
#include "MathLib/KelvinVector.h"
#include "ParameterLib/SpatialPosition.h"
+#include "ProcessLib/Reflection/ReflectionData.h"
namespace ProcessLib::ThermoRichardsMechanics
{
@@ -116,6 +117,13 @@ struct StrainData
{
// TODO Move initialization to the local assembler.
KelvinVector<DisplacementDim> eps = KVzero<DisplacementDim>();
+
+ static auto reflect()
+ {
+ using Self = StrainData<DisplacementDim>;
+
+ return ProcessLib::Reflection::reflectWithName("epsilon", &Self::eps);
+ }
};
} // namespace ProcessLib::ThermoRichardsMechanics
diff --git a/ProcessLib/ThermoRichardsMechanics/Constitutive/DarcyLaw.h b/ProcessLib/ThermoRichardsMechanics/Constitutive/DarcyLaw.h
index 74508febd7f33b641add3116b440d6ee21589d7f..1fc24f14182eafb0ef1034f51f1d7d1a955b062f 100644
--- a/ProcessLib/ThermoRichardsMechanics/Constitutive/DarcyLaw.h
+++ b/ProcessLib/ThermoRichardsMechanics/Constitutive/DarcyLaw.h
@@ -11,7 +11,6 @@
#pragma once
#include "LiquidDensity.h"
-#include "LiquidViscosity.h"
#include "Permeability.h"
#include "ThermoOsmosis.h"
@@ -21,6 +20,14 @@ template <int DisplacementDim>
struct DarcyLawData
{
Eigen::Vector<double, DisplacementDim> v_darcy;
+
+ static auto reflect()
+ {
+ using Self = DarcyLawData<DisplacementDim>;
+
+ return ProcessLib::Reflection::reflectWithName("velocity",
+ &Self::v_darcy);
+ }
};
template <int DisplacementDim>
diff --git a/ProcessLib/ThermoRichardsMechanics/Constitutive/LiquidDensity.h b/ProcessLib/ThermoRichardsMechanics/Constitutive/LiquidDensity.h
index c4de0275b4bdee8d0f073aa57cf2d411074f93e2..9ca326252dba2d2ac459802ca3346367890e419e 100644
--- a/ProcessLib/ThermoRichardsMechanics/Constitutive/LiquidDensity.h
+++ b/ProcessLib/ThermoRichardsMechanics/Constitutive/LiquidDensity.h
@@ -19,6 +19,12 @@ struct LiquidDensityData
double rho_LR;
double drho_LR_dp;
double drho_LR_dT;
+
+ static auto reflect()
+ {
+ return ProcessLib::Reflection::reflectWithName(
+ "liquid_density", &LiquidDensityData::rho_LR);
+ }
};
template <int DisplacementDim>
diff --git a/ProcessLib/ThermoRichardsMechanics/Constitutive/LiquidViscosity.cpp b/ProcessLib/ThermoRichardsMechanics/Constitutive/LiquidViscosity.cpp
index 7c6e0e7ff5ac09b1da5d51a19474657a7a9cf52a..d8997365742f74a168f7c42bc7fb5bcf3d643679 100644
--- a/ProcessLib/ThermoRichardsMechanics/Constitutive/LiquidViscosity.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/Constitutive/LiquidViscosity.cpp
@@ -10,8 +10,6 @@
#include "LiquidViscosity.h"
-#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
-
namespace ProcessLib::ThermoRichardsMechanics
{
template <int DisplacementDim>
diff --git a/ProcessLib/ThermoRichardsMechanics/Constitutive/LiquidViscosity.h b/ProcessLib/ThermoRichardsMechanics/Constitutive/LiquidViscosity.h
index 2f60ad5aaf7a1a6241dc0762dc5388560942c494..6ccdcd370ac95548118a175fa6e2ea0737153291 100644
--- a/ProcessLib/ThermoRichardsMechanics/Constitutive/LiquidViscosity.h
+++ b/ProcessLib/ThermoRichardsMechanics/Constitutive/LiquidViscosity.h
@@ -10,13 +10,19 @@
#pragma once
-#include "Porosity.h"
+#include "ProcessLib/ThermoRichardsMechanics/Constitutive/Base.h"
namespace ProcessLib::ThermoRichardsMechanics
{
struct LiquidViscosityData
{
double viscosity;
+
+ static auto reflect()
+ {
+ return ProcessLib::Reflection::reflectWithName(
+ "viscosity", &LiquidViscosityData::viscosity);
+ }
};
template <int DisplacementDim>
diff --git a/ProcessLib/ThermoRichardsMechanics/Constitutive/Permeability.cpp b/ProcessLib/ThermoRichardsMechanics/Constitutive/Permeability.cpp
index cb7494151c2f9fd83c4dafc70cf529415ab68ab5..03f4229158df7866cfa3cdb7c6f56efa9bdf7f99 100644
--- a/ProcessLib/ThermoRichardsMechanics/Constitutive/Permeability.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/Constitutive/Permeability.cpp
@@ -23,8 +23,8 @@ void PermeabilityModel<DisplacementDim>::eval(
CapillaryPressureData<DisplacementDim> const& p_cap_data,
TemperatureData<DisplacementDim> const& T_data,
PorosityData const& poro_data, LiquidViscosityData const& mu_L_data,
- PorosityData& transport_poro_data,
- PorosityData const& transport_poro_data_prev,
+ TransportPorosityData& transport_poro_data,
+ TransportPorosityData const& transport_poro_data_prev,
SolidMechanicsDataStateless<DisplacementDim> const& s_mech_data,
StrainData<DisplacementDim> const& eps_data,
StrainData<DisplacementDim> const& eps_prev_data,
diff --git a/ProcessLib/ThermoRichardsMechanics/Constitutive/Permeability.h b/ProcessLib/ThermoRichardsMechanics/Constitutive/Permeability.h
index e3a252b373d9d19a731bb92623ed1928454253fd..ae3eb3c5b52e33337595fdd38498e81b5781b94b 100644
--- a/ProcessLib/ThermoRichardsMechanics/Constitutive/Permeability.h
+++ b/ProcessLib/ThermoRichardsMechanics/Constitutive/Permeability.h
@@ -36,8 +36,8 @@ struct PermeabilityModel
PorosityData const& poro_data,
LiquidViscosityData const& mu_L_data,
// TODO evaluate transport porosity evolution separately
- PorosityData& transport_poro_data,
- PorosityData const& transport_poro_data_prev,
+ TransportPorosityData& transport_poro_data,
+ TransportPorosityData const& transport_poro_data_prev,
SolidMechanicsDataStateless<DisplacementDim> const& s_mech_data,
StrainData<DisplacementDim> const& eps_data,
StrainData<DisplacementDim> const& eps_prev_data,
diff --git a/ProcessLib/ThermoRichardsMechanics/Constitutive/Porosity.h b/ProcessLib/ThermoRichardsMechanics/Constitutive/Porosity.h
index 42dfa5c813fab638336cb292d44b539295082605..7bafb13b290f185e1137eb8b106046d33bee9b55 100644
--- a/ProcessLib/ThermoRichardsMechanics/Constitutive/Porosity.h
+++ b/ProcessLib/ThermoRichardsMechanics/Constitutive/Porosity.h
@@ -11,7 +11,6 @@
#pragma once
#include "Bishops.h"
-#include "MathLib/KelvinVector.h"
#include "Saturation.h"
#include "SolidCompressibility.h"
@@ -20,6 +19,23 @@ namespace ProcessLib::ThermoRichardsMechanics
struct PorosityData
{
double phi;
+
+ static auto reflect()
+ {
+ return ProcessLib::Reflection::reflectWithName("porosity",
+ &PorosityData::phi);
+ }
+};
+
+struct TransportPorosityData
+{
+ double phi;
+
+ static auto reflect()
+ {
+ return ProcessLib::Reflection::reflectWithName(
+ "transport_porosity", &TransportPorosityData::phi);
+ }
};
template <int DisplacementDim>
diff --git a/ProcessLib/ThermoRichardsMechanics/Constitutive/Saturation.h b/ProcessLib/ThermoRichardsMechanics/Constitutive/Saturation.h
index eb1c1229a9b72649cae70226555d2b525f31e42b..589c87ecb4cc417c28b8669d8a6b688ea9c7fb45 100644
--- a/ProcessLib/ThermoRichardsMechanics/Constitutive/Saturation.h
+++ b/ProcessLib/ThermoRichardsMechanics/Constitutive/Saturation.h
@@ -17,6 +17,12 @@ namespace ProcessLib::ThermoRichardsMechanics
struct SaturationData
{
double S_L;
+
+ static auto reflect()
+ {
+ return ProcessLib::Reflection::reflectWithName("saturation",
+ &SaturationData::S_L);
+ }
};
struct SaturationDataDeriv
diff --git a/ProcessLib/ThermoRichardsMechanics/Constitutive/SolidDensity.h b/ProcessLib/ThermoRichardsMechanics/Constitutive/SolidDensity.h
index 4ea8bedb59f72f7014af5c5ba04c2f934ccdd181..fb13c344fa11b36143bcba33171881c159705d4f 100644
--- a/ProcessLib/ThermoRichardsMechanics/Constitutive/SolidDensity.h
+++ b/ProcessLib/ThermoRichardsMechanics/Constitutive/SolidDensity.h
@@ -19,6 +19,12 @@ struct SolidDensityData
{
double rho_SR;
double dry_density_solid;
+
+ static auto reflect()
+ {
+ return ProcessLib::Reflection::reflectWithName(
+ "dry_density_solid", &SolidDensityData::dry_density_solid);
+ }
};
template <int DisplacementDim>
diff --git a/ProcessLib/ThermoRichardsMechanics/Constitutive/SolidMechanics.h b/ProcessLib/ThermoRichardsMechanics/Constitutive/SolidMechanics.h
index 3def597ab16aba9fba969076280c2ddf463eb608..ec591359a8b72c10606a7dd359f4f8ad090a5d20 100644
--- a/ProcessLib/ThermoRichardsMechanics/Constitutive/SolidMechanics.h
+++ b/ProcessLib/ThermoRichardsMechanics/Constitutive/SolidMechanics.h
@@ -25,6 +25,15 @@ struct SolidMechanicsDataStateful
// TODO it seems fragile that some data have to be initialized that way.
KelvinVector<DisplacementDim> sigma_eff = KVzero<DisplacementDim>();
KelvinVector<DisplacementDim> eps_m = KVzero<DisplacementDim>();
+
+ static auto reflect()
+ {
+ using Self = SolidMechanicsDataStateful<DisplacementDim>;
+
+ // TODO add eps_m?
+ return ProcessLib::Reflection::reflectWithName("sigma",
+ &Self::sigma_eff);
+ }
};
template <int DisplacementDim>
diff --git a/ProcessLib/ThermoRichardsMechanics/Constitutive/Swelling.h b/ProcessLib/ThermoRichardsMechanics/Constitutive/Swelling.h
index b3a14ad19160fccac81194e450d06f75c98b2a20..54091d39898e7329a46e1322bbf4c384924035b3 100644
--- a/ProcessLib/ThermoRichardsMechanics/Constitutive/Swelling.h
+++ b/ProcessLib/ThermoRichardsMechanics/Constitutive/Swelling.h
@@ -20,6 +20,14 @@ template <int DisplacementDim>
struct SwellingDataStateful
{
KelvinVector<DisplacementDim> sigma_sw = KVzero<DisplacementDim>();
+
+ static auto reflect()
+ {
+ using Self = SwellingDataStateful<DisplacementDim>;
+
+ return ProcessLib::Reflection::reflectWithName("swelling_stress",
+ &Self::sigma_sw);
+ }
};
template <int DisplacementDim>
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveSetting.h b/ProcessLib/ThermoRichardsMechanics/ConstitutiveSetting.h
index 77307c26c2863431794cd491d85799d5e80d463e..2df2974fdbaea1d4fa505a0833f13f7c2a1ef3c0 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveSetting.h
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveSetting.h
@@ -30,16 +30,30 @@
namespace ProcessLib::ThermoRichardsMechanics
{
+
/// Data whose state must be tracked by the TRM process.
template <int DisplacementDim>
struct StatefulData
{
SaturationData S_L_data;
PorosityData poro_data;
- PorosityData transport_poro_data;
+ TransportPorosityData transport_poro_data;
StrainData<DisplacementDim> eps_data;
SwellingDataStateful<DisplacementDim> swelling_data;
SolidMechanicsDataStateful<DisplacementDim> s_mech_data;
+
+ static auto reflect()
+ {
+ using Self = StatefulData<DisplacementDim>;
+
+ return ProcessLib::Reflection::reflectWithoutName(
+ &Self::S_L_data,
+ &Self::poro_data,
+ &Self::transport_poro_data,
+ &Self::eps_data,
+ &Self::swelling_data,
+ &Self::s_mech_data);
+ }
};
/// Data that is needed for output purposes, but not directly for the assembly.
@@ -50,6 +64,16 @@ struct OutputData
LiquidDensityData rho_L_data;
LiquidViscosityData mu_L_data;
SolidDensityData rho_S_data;
+
+ static auto reflect()
+ {
+ using Self = OutputData<DisplacementDim>;
+
+ return ProcessLib::Reflection::reflectWithoutName(&Self::darcy_data,
+ &Self::rho_L_data,
+ &Self::mu_L_data,
+ &Self::rho_S_data);
+ }
};
/// Data that is needed for the equation system assembly.
diff --git a/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h
index 7e87c3ab59923d33fc66e3f187a273c6e42a8cd3..965c03a15ea422fb0c30a040228e4ec8607c78da 100644
--- a/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h
+++ b/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h
@@ -18,7 +18,6 @@
#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveSetting.h"
#include "ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h"
#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
-#include "ProcessLib/Utils/TransposeInPlace.h"
namespace ProcessLib::ThermoRichardsMechanics
{
@@ -119,212 +118,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
return 0;
}
-private:
- std::vector<double> getSigma() const
- {
- constexpr int kelvin_vector_size =
- MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-
- return transposeInPlace<kelvin_vector_size>(
- [this](std::vector<double>& values)
- { return getIntPtSigma(0, {}, {}, values); });
- }
-
- std::vector<double> const& getIntPtSigma(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
- {
- return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
- current_states_, [](auto const& cs) -> auto const& {
- return cs.s_mech_data.sigma_eff;
- },
- cache);
- }
-
- std::vector<double> getSwellingStress() const
- {
- constexpr int kelvin_vector_size =
- MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-
- return transposeInPlace<kelvin_vector_size>(
- [this](std::vector<double>& values)
- { return getIntPtSwellingStress(0, {}, {}, values); });
- }
-
- std::vector<double> const& getIntPtSwellingStress(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
- {
- constexpr int kelvin_vector_size =
- MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
- auto const n_integration_points = current_states_.size();
-
- cache.clear();
- auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
- double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, kelvin_vector_size, n_integration_points);
-
- for (unsigned ip = 0; ip < n_integration_points; ++ip)
- {
- auto const& sigma_sw = current_states_[ip].swelling_data.sigma_sw;
- cache_mat.col(ip) =
- MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_sw);
- }
-
- return cache;
- }
-
- std::vector<double> getEpsilon() const
- {
- constexpr int kelvin_vector_size =
- MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-
- return transposeInPlace<kelvin_vector_size>(
- [this](std::vector<double>& values)
- { return getIntPtEpsilon(0, {}, {}, values); });
- }
-
- std::vector<double> const& getIntPtEpsilon(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
- {
- return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
- current_states_,
- [](auto const& cs) -> auto const& { return cs.eps_data.eps; },
- cache);
- }
-
- std::vector<double> const& getIntPtDarcyVelocity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
- {
- unsigned const n_integration_points =
- integration_method_.getNumberOfPoints();
-
- cache.clear();
- auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
- double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, DisplacementDim, n_integration_points);
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- cache_matrix.col(ip).noalias() =
- output_data_[ip].darcy_data.v_darcy;
- }
-
- return cache;
- }
-
- std::vector<double> getSaturation() const
- {
- std::vector<double> result;
- getIntPtSaturation(0, {}, {}, result);
- return result;
- }
-
- std::vector<double> const& getIntPtSaturation(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
- {
- return ProcessLib::getIntegrationPointScalarData(
- current_states_,
- [](auto const& state) -> auto const& { return state.S_L_data.S_L; },
- cache);
- }
-
- std::vector<double> getPorosity() const
- {
- std::vector<double> result;
- getIntPtPorosity(0, {}, {}, result);
- return result;
- }
-
- std::vector<double> const& getIntPtPorosity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
- {
- return ProcessLib::getIntegrationPointScalarData(
- current_states_, [](auto const& state) -> auto const& {
- return state.poro_data.phi;
- },
- cache);
- }
-
- std::vector<double> getTransportPorosity() const
- {
- std::vector<double> result;
- getIntPtTransportPorosity(0, {}, {}, result);
- return result;
- }
-
- std::vector<double> const& getIntPtTransportPorosity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
- {
- return ProcessLib::getIntegrationPointScalarData(
- current_states_,
- [](auto const& state) -> auto const& {
- return state.transport_poro_data.phi;
- },
- cache);
- }
-
- std::vector<double> const& getIntPtDryDensitySolid(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
- {
- return ProcessLib::getIntegrationPointScalarData(
- output_data_,
- [](auto const& out) -> auto const& {
- return out.rho_S_data.dry_density_solid;
- },
- cache);
- }
-
- std::vector<double> const& getIntPtLiquidDensity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
- {
- return ProcessLib::getIntegrationPointScalarData(
- output_data_,
- [](auto const& out) -> auto const& {
- return out.rho_L_data.rho_LR;
- },
- cache);
- }
-
- std::vector<double> const& getIntPtViscosity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
- {
- return ProcessLib::getIntegrationPointScalarData(
- output_data_, [](auto const& out) -> auto const& {
- return out.mu_L_data.viscosity;
- },
- cache);
- }
-
-public:
// TODO move to NumLib::ExtrapolatableElement
unsigned getNumberOfIntegrationPoints() const
{
@@ -354,65 +147,12 @@ public:
prev_states_ = current_states_;
}
- struct IPDataAccessorForExtrapolation
- {
- std::string name;
- unsigned num_comp;
- std::function<std::vector<double> const&(
- LocalAssemblerInterface<DisplacementDim> const&,
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<
- NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& /*cache*/)>
- ip_data_accessor;
- };
-
- static std::vector<IPDataAccessorForExtrapolation>
- getIPDataAccessorsForExtrapolation()
- {
- using Self = LocalAssemblerInterface<DisplacementDim>;
- constexpr auto kv_size =
- MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-
- return {{"sigma", kv_size, &Self::getIntPtSigma},
- {"swelling_stress", kv_size, &Self::getIntPtSwellingStress},
- {"epsilon", kv_size, &Self::getIntPtEpsilon},
- {"velocity", DisplacementDim, &Self::getIntPtDarcyVelocity},
- {"saturation", 1, &Self::getIntPtSaturation},
- {"porosity", 1, &Self::getIntPtPorosity},
- {"transport_porosity", 1, &Self::getIntPtTransportPorosity},
- {"dry_density_solid", 1, &Self::getIntPtDryDensitySolid},
- {"liquid_density", 1, &Self::getIntPtLiquidDensity},
- {"viscosity", 1, &Self::getIntPtViscosity}};
- }
-
- struct IPDataAccessorForIPWriter
- {
- std::string name;
- unsigned num_comp;
- std::vector<double> (LocalAssemblerInterface<DisplacementDim>::*
- ip_data_accessor)() const;
- };
-
- static std::vector<IPDataAccessorForIPWriter>
- getIPDataAccessorsForIPWriter()
+ static auto getReflectionDataForOutput()
{
using Self = LocalAssemblerInterface<DisplacementDim>;
- constexpr auto kv_size =
- MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
- // TODO (naumov) remove ip suffix. Probably needs modification of the
- // mesh properties, s.t. there is no "overlapping" with cell/point data.
- // See getOrCreateMeshProperty.
- return {
- {"sigma_ip", kv_size, &Self::getSigma},
- {"saturation_ip", 1, &Self::getSaturation},
- {"porosity_ip", 1, &Self::getPorosity},
- {"transport_porosity_ip", 1, &Self::getTransportPorosity},
- {"swelling_stress_ip", kv_size, &Self::getSwellingStress},
- {"epsilon_ip", kv_size, &Self::getEpsilon},
- };
+ return ProcessLib::Reflection::reflectWithoutName(
+ &Self::current_states_, &Self::output_data_);
}
protected:
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
index fc1d1b0dbd1c8567249042589f7f8cef084cdaca..439a664a56c5cc73a74e6cb3fdc9cab4ef9bc986 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
@@ -17,6 +17,8 @@
#include "NumLib/DOF/DOFTableUtil.h"
#include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
#include "ProcessLib/Process.h"
+#include "ProcessLib/Reflection/ReflectionForExtrapolation.h"
+#include "ProcessLib/Reflection/ReflectionForIPWriters.h"
#include "ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h"
#include "ProcessLib/Utils/SetIPDataInitialConditions.h"
#include "ThermoRichardsMechanicsFEM.h"
@@ -51,14 +53,10 @@ ThermoRichardsMechanicsProcess<DisplacementDim>::ThermoRichardsMechanicsProcess(
heat_flux_ = MeshLib::getOrCreateMeshProperty<double>(
mesh, "HeatFlowRate", MeshLib::MeshItemType::Node, 1);
- for (auto& [name, num_comp, ip_data_accessor] :
- LocalAssemblerIF::getIPDataAccessorsForIPWriter())
- {
- _integration_point_writer.emplace_back(
- std::make_unique<MeshLib::IntegrationPointWriter>(
- name, num_comp, integration_order, local_assemblers_,
- ip_data_accessor));
- }
+ ProcessLib::Reflection::addReflectedIntegrationPointWriters<
+ DisplacementDim>(LocalAssemblerIF::getReflectionDataForOutput(),
+ _integration_point_writer, integration_order,
+ local_assemblers_);
}
template <int DisplacementDim>
@@ -143,11 +141,9 @@ void ThermoRichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
std::move(get_ip_values_function)));
};
- for (auto& [name, num_comp, fct] : LocalAssemblerInterface<
- DisplacementDim>::getIPDataAccessorsForExtrapolation())
- {
- add_secondary_variable(name, num_comp, fct);
- };
+ ProcessLib::Reflection::addReflectedSecondaryVariables<DisplacementDim>(
+ LocalAssemblerIF::getReflectionDataForOutput(), _secondary_variables,
+ getExtrapolator(), local_assemblers_);
//
// enable output of internal variables defined by material models
diff --git a/ProcessLib/Utils/TransposeInPlace.h b/ProcessLib/Utils/TransposeInPlace.h
index 163ebf9713f2800ae344744777a130f1d04926b3..6b8b338ef5a58f1d940a9ca8885465aee3397336 100644
--- a/ProcessLib/Utils/TransposeInPlace.h
+++ b/ProcessLib/Utils/TransposeInPlace.h
@@ -37,4 +37,32 @@ std::vector<double> transposeInPlace(
return result;
}
+template <int Components>
+void transposeInPlace(std::vector<double>& values)
+{
+ MathLib::toMatrix<
+ Eigen::Matrix<double, Eigen::Dynamic, Components, Eigen::RowMajor>>(
+ values, values.size() / Components, Components) =
+ MathLib::toMatrix<
+ Eigen::Matrix<double, Components, Eigen::Dynamic, Eigen::RowMajor>>(
+ values, Components, values.size() / Components)
+ .transpose()
+ .eval();
+}
+
+inline void transposeInPlace(std::vector<double>& values,
+ unsigned const num_components)
+{
+ MathLib::toMatrix<
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
+ values, values.size() / num_components, num_components) =
+ MathLib::toMatrix<Eigen::Matrix<double,
+ Eigen::Dynamic,
+ Eigen::Dynamic,
+ Eigen::RowMajor>>(
+ values, num_components, values.size() / num_components)
+ .transpose()
+ .eval();
+}
+
} // namespace ProcessLib
diff --git a/Tests/ProcessLib/TestReflectIPData.cpp b/Tests/ProcessLib/TestReflectIPData.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3dd1f8b47c7bcc5a40e32dacf39c3891451267f8
--- /dev/null
+++ b/Tests/ProcessLib/TestReflectIPData.cpp
@@ -0,0 +1,487 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, 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 <gmock/gmock-matchers.h>
+#include <gtest/gtest.h>
+
+#include <numeric>
+#include <random>
+
+#include "ProcessLib/Reflection/ReflectionIPData.h"
+
+template <int Dim>
+struct Level3
+{
+ MathLib::KelvinVector::KelvinVectorType<Dim> kelvin3;
+ Eigen::Vector<double, Dim> vector3;
+ double scalar3;
+
+ static auto reflect()
+ {
+ using namespace ProcessLib::Reflection;
+ return std::tuple{ReflectionData{"kelvin3", &Level3::kelvin3},
+ ReflectionData{"vector3", &Level3::vector3},
+ ReflectionData{"scalar3", &Level3::scalar3}};
+ }
+};
+
+template <int Dim>
+struct Level3b
+{
+ double scalar3b;
+
+ static auto reflect()
+ {
+ using namespace ProcessLib::Reflection;
+ return std::tuple{ReflectionData{"scalar3b", &Level3b::scalar3b}};
+ }
+};
+
+template <int Dim>
+struct Level2
+{
+ Level3<Dim> level3;
+ Level3b<Dim> level3b;
+
+ static auto reflect()
+ {
+ using namespace ProcessLib::Reflection;
+ return std::tuple{ReflectionData{&Level2::level3},
+ ReflectionData{&Level2::level3b}};
+ }
+};
+
+template <int Dim>
+struct Level2b
+{
+ double scalar2b;
+
+ static auto reflect()
+ {
+ using namespace ProcessLib::Reflection;
+ return std::tuple{ReflectionData{"scalar2b", &Level2b::scalar2b}};
+ }
+};
+
+template <int Dim>
+struct Level1
+{
+ MathLib::KelvinVector::KelvinVectorType<Dim> kelvin1;
+ Eigen::Vector<double, Dim> vector1;
+ double scalar1;
+ Level2<Dim> level2;
+ Level2b<Dim> level2b;
+
+ static auto reflect()
+ {
+ using namespace ProcessLib::Reflection;
+ return std::tuple{ReflectionData{"kelvin1", &Level1::kelvin1},
+ ReflectionData{"vector1", &Level1::vector1},
+ ReflectionData{"scalar1", &Level1::scalar1},
+ ReflectionData{&Level1::level2},
+ ReflectionData{&Level1::level2b}};
+ }
+};
+
+template <int Dim>
+struct Level1b
+{
+ double scalar1b;
+
+ static auto reflect()
+ {
+ using namespace ProcessLib::Reflection;
+ return std::tuple{ReflectionData{"scalar1b", &Level1b::scalar1b}};
+ }
+};
+
+template <int Dim>
+struct LocAsmIF
+{
+ explicit LocAsmIF(unsigned const num_ips)
+ : ip_data_scalar(num_ips),
+ ip_data_vector(num_ips),
+ ip_data_kelvin(num_ips),
+ ip_data_level1(num_ips),
+ ip_data_level1b(num_ips)
+ {
+ }
+
+ std::size_t numIPs() const { return ip_data_scalar.size(); }
+
+ std::vector<double> ip_data_scalar;
+ std::vector<Eigen::Vector<double, Dim>> ip_data_vector;
+ std::vector<MathLib::KelvinVector::KelvinVectorType<Dim>> ip_data_kelvin;
+ std::vector<Level1<Dim>> ip_data_level1;
+ std::vector<Level1b<Dim>> ip_data_level1b;
+
+ static auto reflect()
+ {
+ using namespace ProcessLib::Reflection;
+ return std::tuple{ReflectionData{"scalar", &LocAsmIF::ip_data_scalar},
+ ReflectionData{"vector", &LocAsmIF::ip_data_vector},
+ ReflectionData{"kelvin", &LocAsmIF::ip_data_kelvin},
+ ReflectionData{&LocAsmIF::ip_data_level1},
+ ReflectionData{&LocAsmIF::ip_data_level1b}};
+ }
+};
+
+template <int dim>
+struct NumCompAndFunction
+{
+ unsigned num_comp;
+ std::function<std::vector<double>(LocAsmIF<dim> const&)> function;
+};
+
+// Prepares scalar IP data for the passed local assembler.
+//
+// The IP data are a sequence of double values starting at the passed start
+// value and incremented by one for each integration point.
+//
+// The location of the prepared data is specified by the IP data accessor
+// callback function.
+//
+// Returns the expected data for use in unit test checks.
+template <int dim>
+std::vector<double> initScalar(LocAsmIF<dim>& loc_asm,
+ double const start_value,
+ auto const ip_data_accessor,
+ bool const for_read_test)
+{
+ auto const num_int_pts = loc_asm.numIPs();
+
+ // init ip data in the local assembler
+ if (for_read_test)
+ {
+ for (std::size_t ip = 0; ip < num_int_pts; ++ip)
+ {
+ ip_data_accessor(loc_asm, ip) = start_value + ip;
+ }
+ }
+ else
+ {
+ for (std::size_t ip = 0; ip < num_int_pts; ++ip)
+ {
+ ip_data_accessor(loc_asm, ip) =
+ std::numeric_limits<double>::quiet_NaN();
+ }
+ }
+
+ // prepare reference data
+ std::vector<double> scalar_expected(num_int_pts);
+ iota(begin(scalar_expected), end(scalar_expected), start_value);
+ return scalar_expected;
+}
+
+// Prepares vector valued IP data for the passed local assembler.
+//
+// The IP data are a sequence of double values starting at the passed start
+// value and incremented by one for each integration point and vector
+// component.
+//
+// The location of the prepared data is specified by the IP data accessor
+// callback function.
+//
+// Returns the expected data for use in unit test checks.
+template <int dim>
+std::vector<double> initVector(LocAsmIF<dim>& loc_asm,
+ double const start_value,
+ auto const ip_data_accessor,
+ bool const for_read_test)
+{
+ auto const num_int_pts = loc_asm.numIPs();
+
+ // init ip data in the local assembler
+ if (for_read_test)
+ {
+ for (std::size_t ip = 0; ip < num_int_pts; ++ip)
+ {
+ ip_data_accessor(loc_asm, ip) =
+ Eigen::Vector<double, dim>::LinSpaced(
+ dim, ip * dim + start_value,
+ ip * dim + start_value - 1 + dim);
+ }
+ }
+ else
+ {
+ for (std::size_t ip = 0; ip < num_int_pts; ++ip)
+ {
+ ip_data_accessor(loc_asm, ip) =
+ Eigen::Vector<double, dim>::Constant(
+ std::numeric_limits<double>::quiet_NaN());
+ }
+ }
+
+ // prepare reference data
+ std::vector<double> vector_expected(num_int_pts * dim);
+ iota(begin(vector_expected), end(vector_expected), start_value);
+ return vector_expected;
+}
+
+// Prepares Kelvin vector valued IP data for the passed local assembler.
+//
+// The IP data are a sequence of double values starting at the passed start
+// value and incremented by one for each integration point and Kelvin vector
+// component.
+//
+// The location of the prepared data is specified by the IP data accessor
+// callback function.
+//
+// Returns the expected data for use in unit test checks.
+template <int dim>
+std::vector<double> initKelvin(LocAsmIF<dim>& loc_asm,
+ double const start_value,
+ auto const ip_data_accessor,
+ bool const for_read_test)
+{
+ auto constexpr kv_size =
+ MathLib::KelvinVector::kelvin_vector_dimensions(dim);
+
+ auto const num_int_pts = loc_asm.numIPs();
+
+ // init ip data in the local assembler
+ if (for_read_test)
+ {
+ for (std::size_t ip = 0; ip < num_int_pts; ++ip)
+ {
+ ip_data_accessor(loc_asm, ip) =
+ MathLib::KelvinVector::symmetricTensorToKelvinVector(
+ Eigen::Vector<double, kv_size>::LinSpaced(
+ kv_size, ip * kv_size + start_value,
+ ip * kv_size + start_value - 1 + kv_size));
+ }
+ }
+ else
+ {
+ for (std::size_t ip = 0; ip < num_int_pts; ++ip)
+ {
+ ip_data_accessor(loc_asm, ip) =
+ Eigen::Vector<double, kv_size>::Constant(
+ std::numeric_limits<double>::quiet_NaN());
+ }
+ }
+
+ // prepare reference data
+ std::vector<double> vector_expected(num_int_pts * kv_size);
+ iota(begin(vector_expected), end(vector_expected), start_value);
+ return vector_expected;
+}
+
+template <int dim>
+struct ReferenceData
+{
+ std::vector<double> scalar;
+ std::vector<double> vector;
+ std::vector<double> kelvin;
+
+ std::vector<double> scalar1;
+ std::vector<double> vector1;
+ std::vector<double> kelvin1;
+
+ std::vector<double> scalar3;
+ std::vector<double> vector3;
+ std::vector<double> kelvin3;
+
+ std::vector<double> scalar1b;
+ std::vector<double> scalar2b;
+ std::vector<double> scalar3b;
+
+ static ReferenceData<dim> create(LocAsmIF<dim>& loc_asm,
+ bool const for_read_test)
+ {
+ std::random_device ran_dev;
+ std::mt19937 ran_gen(ran_dev());
+ std::uniform_real_distribution<> ran_dist(1.0, 2.0);
+ auto start_value = [&]() { return ran_dist(ran_gen); };
+
+ ReferenceData<dim> ref;
+
+ // level 0 - data preparation //////////////////////////////////////////
+
+ ref.scalar = initScalar(
+ loc_asm, start_value(),
+ [](auto& loc_asm, unsigned const ip) -> auto& {
+ return loc_asm.ip_data_scalar[ip];
+ },
+ for_read_test);
+
+ ref.vector = initVector(
+ loc_asm, start_value(),
+ [](auto& loc_asm, unsigned const ip) -> auto& {
+ return loc_asm.ip_data_vector[ip];
+ },
+ for_read_test);
+
+ ref.kelvin = initKelvin(
+ loc_asm, start_value(),
+ [](auto& loc_asm, unsigned const ip) -> auto& {
+ return loc_asm.ip_data_kelvin[ip];
+ },
+ for_read_test);
+
+ // level 1 - data preparation //////////////////////////////////////////
+
+ ref.scalar1 = initScalar(
+ loc_asm, start_value(),
+ [](auto& loc_asm, unsigned const ip) -> auto& {
+ return loc_asm.ip_data_level1[ip].scalar1;
+ },
+ for_read_test);
+
+ ref.vector1 = initVector(
+ loc_asm, start_value(),
+ [](auto& loc_asm, unsigned const ip) -> auto& {
+ return loc_asm.ip_data_level1[ip].vector1;
+ },
+ for_read_test);
+
+ ref.kelvin1 = initKelvin(
+ loc_asm, start_value(),
+ [](auto& loc_asm, unsigned const ip) -> auto& {
+ return loc_asm.ip_data_level1[ip].kelvin1;
+ },
+ for_read_test);
+
+ // level 3 - data preparation //////////////////////////////////////////
+
+ ref.scalar3 = initScalar(
+ loc_asm, start_value(),
+ [](auto& loc_asm, unsigned const ip) -> auto& {
+ return loc_asm.ip_data_level1[ip].level2.level3.scalar3;
+ },
+ for_read_test);
+
+ ref.vector3 = initVector(
+ loc_asm, start_value(),
+ [](auto& loc_asm, unsigned const ip) -> auto& {
+ return loc_asm.ip_data_level1[ip].level2.level3.vector3;
+ },
+ for_read_test);
+
+ ref.kelvin3 = initKelvin(
+ loc_asm, start_value(),
+ [](auto& loc_asm, unsigned const ip) -> auto& {
+ return loc_asm.ip_data_level1[ip].level2.level3.kelvin3;
+ },
+ for_read_test);
+
+ // b levels - data preparation /////////////////////////////////////////
+ // b levels test that the reflection implementation recurses on multiple
+ // members, not only on one.
+
+ ref.scalar1b = initScalar(
+ loc_asm, start_value(),
+ [](auto& loc_asm, unsigned const ip) -> auto& {
+ return loc_asm.ip_data_level1b[ip].scalar1b;
+ },
+ for_read_test);
+
+ ref.scalar2b = initScalar(
+ loc_asm, start_value(),
+ [](auto& loc_asm, unsigned const ip) -> auto& {
+ return loc_asm.ip_data_level1[ip].level2b.scalar2b;
+ },
+ for_read_test);
+
+ ref.scalar3b = initScalar(
+ loc_asm, start_value(),
+ [](auto& loc_asm, unsigned const ip) -> auto& {
+ return loc_asm.ip_data_level1[ip].level2.level3b.scalar3b;
+ },
+ for_read_test);
+
+ return ref;
+ }
+};
+
+template <class Dim>
+struct ProcessLib_ReflectIPData : ::testing::Test
+{
+ static constexpr auto dim = Dim::value;
+};
+
+using ProcessLib_ReflectIPData_TestCases =
+ ::testing::Types<std::integral_constant<int, 2>,
+ std::integral_constant<int, 3>>;
+
+TYPED_TEST_SUITE(ProcessLib_ReflectIPData, ProcessLib_ReflectIPData_TestCases);
+
+TYPED_TEST(ProcessLib_ReflectIPData, ReadTest)
+{
+ constexpr int dim = TypeParam::value;
+ auto constexpr kv_size =
+ MathLib::KelvinVector::kelvin_vector_dimensions(dim);
+
+ using LocAsm = LocAsmIF<dim>;
+
+ std::size_t const num_int_pts = 8;
+ LocAsm loc_asm(num_int_pts);
+
+ auto const ref = ReferenceData<dim>::create(loc_asm, true);
+
+ // function under test /////////////////////////////////////////////////////
+
+ std::map<std::string, NumCompAndFunction<dim>>
+ map_name_to_num_comp_and_function;
+
+ ProcessLib::Reflection::forEachReflectedFlattenedIPDataAccessor<dim,
+ LocAsm>(
+ LocAsm::reflect(),
+ [&map_name_to_num_comp_and_function](std::string const& name,
+ unsigned const num_comp,
+ auto&& double_vec_from_loc_asm)
+ {
+ EXPECT_FALSE(map_name_to_num_comp_and_function.contains(name));
+ map_name_to_num_comp_and_function[name] = {
+ num_comp, std::move(double_vec_from_loc_asm)};
+ });
+
+ // checks //////////////////////////////////////////////////////////////////
+
+ auto check = [&map_name_to_num_comp_and_function, &loc_asm](
+ std::string const& name,
+ unsigned const num_comp_expected,
+ std::vector<double> const& values_expected)
+ {
+ auto const it = map_name_to_num_comp_and_function.find(name);
+
+ ASSERT_NE(map_name_to_num_comp_and_function.end(), it)
+ << "No accessor found for ip data with name '" << name << "'";
+
+ auto const& [num_comp, fct] = it->second;
+
+ EXPECT_EQ(num_comp_expected, num_comp)
+ << "Number of components differs for ip data with name '" << name
+ << "'";
+ EXPECT_THAT(fct(loc_asm),
+ testing::Pointwise(testing::DoubleEq(), values_expected))
+ << "Values differ for ip data with name '" << name << "'";
+ };
+
+ // level 0
+ check("scalar", 1, ref.scalar);
+ check("vector", dim, ref.vector);
+ check("kelvin", kv_size, ref.kelvin);
+
+ // level 1
+ check("scalar1", 1, ref.scalar1);
+ check("vector1", dim, ref.vector1);
+ check("kelvin1", kv_size, ref.kelvin1);
+
+ // level 3
+ check("scalar3", 1, ref.scalar3);
+ check("vector3", dim, ref.vector3);
+ check("kelvin3", kv_size, ref.kelvin3);
+
+ // b levels
+ check("scalar1b", 1, ref.scalar1b);
+ check("scalar2b", 1, ref.scalar2b);
+ check("scalar3b", 1, ref.scalar3b);
+}
diff --git a/Tests/ProcessLib/TestTransposeInplace.cpp b/Tests/ProcessLib/TestTransposeInplace.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d35bc7abe6e9ccf84fbc8e8f52c1ff399d772acc
--- /dev/null
+++ b/Tests/ProcessLib/TestTransposeInplace.cpp
@@ -0,0 +1,59 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2022, 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 <gmock/gmock.h>
+#include <gtest/gtest.h>
+
+#include "ProcessLib/Utils/TransposeInPlace.h"
+
+TEST(ProcessLib, TransposeInplaceFixed)
+{
+ std::vector<double> const values_in{
+ 1, 2, 3, 4, //
+ 5, 6, 7, 8, //
+ 9, 10, 11, 12 //
+ };
+ constexpr unsigned num_rows = 3;
+
+ std::vector<double> const values_out_expected{
+ 1, 5, 9, //
+ 2, 6, 10, //
+ 3, 7, 11, //
+ 4, 8, 12 //
+ };
+
+ std::vector<double> values_out_actual(values_in);
+
+ ProcessLib::transposeInPlace<num_rows>(values_out_actual);
+
+ ASSERT_THAT(values_out_actual, testing::ContainerEq(values_out_expected));
+}
+
+TEST(ProcessLib, TransposeInplaceDynamic)
+{
+ std::vector<double> const values_in{
+ 1, 2, 3, 4, //
+ 5, 6, 7, 8, //
+ 9, 10, 11, 12 //
+ };
+ unsigned const num_rows = 3;
+
+ std::vector<double> const values_out_expected{
+ 1, 5, 9, //
+ 2, 6, 10, //
+ 3, 7, 11, //
+ 4, 8, 12 //
+ };
+
+ std::vector<double> values_out_actual(values_in);
+
+ ProcessLib::transposeInPlace(values_out_actual, num_rows);
+
+ ASSERT_THAT(values_out_actual, testing::ContainerEq(values_out_expected));
+}