diff --git a/ProcessLib/HydroMechanics/Tests.cmake b/ProcessLib/HydroMechanics/Tests.cmake
index 74e681640abe97b2f01efb89d1c8a5a7e83f7cd6..14293466cd679ddc7b4dfc0a675100b5965cf546 100644
--- a/ProcessLib/HydroMechanics/Tests.cmake
+++ b/ProcessLib/HydroMechanics/Tests.cmake
@@ -277,7 +277,7 @@ AddTest(
REQUIREMENTS NOT OGS_USE_MPI
DIFF_DATA
expected_square_1e2_UC_late_ts_10_t_1000.000000.vtu square_1e2_UC_late_ts_10_t_1000.000000.vtu displacement displacement 1e-13 1e-16
- expected_square_1e2_UC_late_ts_10_t_1000.000000.vtu square_1e2_UC_late_ts_10_t_1000.000000.vtu pressure pressure 1e-13 1e-16
+ expected_square_1e2_UC_late_ts_10_t_1000.000000.vtu square_1e2_UC_late_ts_10_t_1000.000000.vtu pressure pressure 2e-13 1e-16
)
AddTest(
diff --git a/ProcessLib/Utils/SetOrGetIntegrationPointData.h b/ProcessLib/Utils/SetOrGetIntegrationPointData.h
index cc355456bb4aa8c8d989a94f8691b3e3a690f6f2..451d04ba06bd67286308f44c9eb79a7de78e421f 100644
--- a/ProcessLib/Utils/SetOrGetIntegrationPointData.h
+++ b/ProcessLib/Utils/SetOrGetIntegrationPointData.h
@@ -16,7 +16,7 @@
#include "BaseLib/DynamicSpan.h"
#include "MathLib/KelvinVector.h"
-#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "TransposeInPlace.h"
namespace ProcessLib
{
@@ -65,17 +65,25 @@ std::vector<double> const& getIntegrationPointKelvinVectorData(
return cache;
}
+//! Overload without \c cache argument.
+//!
+//! \note This function returns the data in transposed storage order compared to
+//! the overloads that have a \c cache argument.
template <int DisplacementDim, typename IntegrationPointDataVector,
typename MemberType>
std::vector<double> getIntegrationPointKelvinVectorData(
IntegrationPointDataVector const& ip_data_vector, MemberType member)
{
- std::vector<double> ip_kelvin_vector_values;
-
- getIntegrationPointKelvinVectorData<DisplacementDim>(
- ip_data_vector, member, ip_kelvin_vector_values);
+ constexpr int kelvin_vector_size =
+ MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
- return ip_kelvin_vector_values;
+ return transposeInPlace<kelvin_vector_size>(
+ [&](std::vector<double>& values)
+ {
+ return getIntegrationPointKelvinVectorData<DisplacementDim>(
+ ip_data_vector, member, values);
+ ;
+ });
}
template <int DisplacementDim, typename IntegrationPointDataVector,
@@ -207,9 +215,9 @@ std::vector<double> getIntegrationPointDataMaterialStateVariables(
std::vector<double> result;
result.reserve(ip_data_vector.size() * n_components);
- for (auto& ip_data_vector : ip_data_vector)
+ for (auto& ip_data : ip_data_vector)
{
- auto const values_span = get_values_span(*(ip_data_vector.*member));
+ auto const values_span = get_values_span(*(ip_data.*member));
assert(values_span.size() == static_cast<std::size_t>(n_components));
result.insert(end(result), values_span.begin(), values_span.end());
@@ -230,7 +238,7 @@ std::size_t setIntegrationPointDataMaterialStateVariables(
auto const n_integration_points = ip_data_vector.size();
std::size_t position = 0;
- for (auto& ip_data : ip_data_vector)
+ for (auto const& ip_data : ip_data_vector)
{
auto const values_span = get_values_span(*(ip_data.*member));
std::copy_n(values + position, values_span.size(), values_span.begin());
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloan/slope.prj b/Tests/Data/Mechanics/MohrCoulombAbboSloan/slope.prj
index d39414f39bff882533b05a1be937d83e685217b7..b25f8278a0449c1c5b1284ba46438e019215c4a5 100644
--- a/Tests/Data/Mechanics/MohrCoulombAbboSloan/slope.prj
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloan/slope.prj
@@ -41,7 +41,7 @@
<convergence_criterion>
<type>DeltaX</type>
<norm_type>NORM2</norm_type>
- <reltol>1e-14</reltol>
+ <reltol>1e-13</reltol>
</convergence_criterion>
<time_discretization>
<type>BackwardEuler</type>
diff --git a/Tests/ProcessLib/TestIPDataAccess.cpp b/Tests/ProcessLib/TestIPDataAccess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5f8a3d162b58e468068668088b3a2f457d6bd3fb
--- /dev/null
+++ b/Tests/ProcessLib/TestIPDataAccess.cpp
@@ -0,0 +1,229 @@
+/**
+ * \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 "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
+
+template <int DisplacementDim>
+struct IPData
+{
+ using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
+
+ KV kelvin;
+ double scalar;
+};
+
+template <class Dim>
+struct ProcessLib_IPDataAccess : ::testing::Test
+{
+ static constexpr int dim = Dim::value;
+ static constexpr int kv_size =
+ MathLib::KelvinVector::kelvin_vector_dimensions(dim);
+
+ static std::vector<IPData<Dim::value>> getIPData()
+ {
+ using KV = typename IPData<dim>::KV;
+
+ constexpr int off_diag_size = dim == 2 ? 1 : 3;
+
+ constexpr std::size_t num_int_pts = 10;
+
+ std::vector<IPData<dim>> ip_data(num_int_pts);
+
+ for (std::size_t i = 0; i < num_int_pts; ++i)
+ {
+ ip_data[i].kelvin =
+ KV::Constant(10. * i) + KV::LinSpaced(0., kv_size - 1.);
+
+ // compensate Kelvin vector <-> symmetric tensor conversion
+ ip_data[i].kelvin.template tail<off_diag_size>() *= std::sqrt(2.0);
+
+ ip_data[i].scalar = 10. * num_int_pts + i;
+ }
+
+ return ip_data;
+ }
+
+ static std::vector<IPData<Dim::value>> getIPDataNaNs()
+ {
+ using KV = typename IPData<dim>::KV;
+
+ constexpr std::size_t num_int_pts = 10;
+ constexpr double nan = std::numeric_limits<double>::quiet_NaN();
+
+ std::vector<IPData<dim>> ip_data(num_int_pts);
+
+ for (std::size_t i = 0; i < num_int_pts; ++i)
+ {
+ ip_data[i].kelvin = KV::Constant(nan);
+ ip_data[i].scalar = nan;
+ }
+
+ return ip_data;
+ }
+
+ static std::vector<double> getScalarData()
+ {
+ return {100, 101, 102, 103, 104, 105, 106, 107, 108, 109};
+ }
+
+ static std::vector<double> getKVDataDefaultOrder()
+ {
+ if constexpr (dim == 2)
+ {
+ return {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, // 1st comp
+ 1, 11, 21, 31, 41, 51, 61, 71, 81, 91, // 2nd comp
+ 2, 12, 22, 32, 42, 52, 62, 72, 82, 92, // 3rd comp
+ 3, 13, 23, 33, 43, 53, 63, 73, 83, 93}; // 4th comp
+ }
+ else if constexpr (dim == 3)
+ {
+ return {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, //
+ 1, 11, 21, 31, 41, 51, 61, 71, 81, 91, //
+ 2, 12, 22, 32, 42, 52, 62, 72, 82, 92, //
+ 3, 13, 23, 33, 43, 53, 63, 73, 83, 93, //
+ 4, 14, 24, 34, 44, 54, 64, 74, 84, 94, //
+ 5, 15, 25, 35, 45, 55, 65, 75, 85, 95};
+ };
+ }
+
+ static std::vector<double> getKVDataTransposedOrder()
+ {
+ if constexpr (dim == 2)
+ {
+ return {0, 1, 2, 3, // 1st IP
+ 10, 11, 12, 13, // 2nd IP
+ 20, 21, 22, 23, // 3rd IP
+ 30, 31, 32, 33, // ...
+ 40, 41, 42, 43, //
+ 50, 51, 52, 53, //
+ 60, 61, 62, 63, //
+ 70, 71, 72, 73, //
+ 80, 81, 82, 83, //
+ 90, 91, 92, 93};
+ }
+ else if constexpr (dim == 3)
+ {
+ return {
+ 0, 1, 2, 3, 4, 5, //
+ 10, 11, 12, 13, 14, 15, //
+ 20, 21, 22, 23, 24, 25, //
+ 30, 31, 32, 33, 34, 35, //
+ 40, 41, 42, 43, 44, 45, //
+ 50, 51, 52, 53, 54, 55, //
+ 60, 61, 62, 63, 64, 65, //
+ 70, 71, 72, 73, 74, 75, //
+ 80, 81, 82, 83, 84, 85, //
+ 90, 91, 92, 93, 94, 95,
+ };
+ };
+ }
+};
+
+using ProcessLib_IPDataAccess_TestCases =
+ ::testing::Types<std::integral_constant<int, 2>,
+ std::integral_constant<int, 3>>;
+
+TYPED_TEST_SUITE(ProcessLib_IPDataAccess, ProcessLib_IPDataAccess_TestCases);
+
+TYPED_TEST(ProcessLib_IPDataAccess, GetScalarData)
+{
+ constexpr int dim = TypeParam::value;
+
+ auto const ip_data = this->getIPData();
+
+ std::vector<double> cache;
+
+ ProcessLib::getIntegrationPointScalarData(
+ ip_data, &IPData<dim>::scalar, cache);
+
+ ASSERT_THAT(cache,
+ testing::Pointwise(testing::DoubleEq(), this->getScalarData()));
+}
+
+TYPED_TEST(ProcessLib_IPDataAccess, GetKelvinVectorDataDefaultOrder)
+{
+ constexpr int dim = TypeParam::value;
+
+ auto const ip_data = this->getIPData();
+
+ std::vector<double> cache;
+
+ ProcessLib::getIntegrationPointKelvinVectorData<dim>(
+ ip_data, &IPData<dim>::kelvin, cache);
+
+ ASSERT_THAT(
+ cache,
+ testing::Pointwise(testing::DoubleEq(), this->getKVDataDefaultOrder()));
+}
+
+TYPED_TEST(ProcessLib_IPDataAccess, GetKelvinVectorDataTransposedOrder)
+{
+ constexpr int dim = TypeParam::value;
+
+ auto const ip_data = this->getIPData();
+
+ const std::vector<double> cache =
+ ProcessLib::getIntegrationPointKelvinVectorData<dim>(
+ ip_data,
+ &IPData<dim>::kelvin); // pretty subtle: if no cache argument is
+ // passed, data is returned transposed
+
+ ASSERT_THAT(cache,
+ testing::Pointwise(testing::DoubleEq(),
+ this->getKVDataTransposedOrder()));
+}
+
+TYPED_TEST(ProcessLib_IPDataAccess, SetScalarData)
+{
+ constexpr int dim = TypeParam::value;
+
+ auto ip_data = this->getIPDataNaNs();
+
+ auto const cache = this->getScalarData();
+
+ auto const num_read = ProcessLib::setIntegrationPointScalarData(
+ &cache[0], ip_data, &IPData<dim>::scalar);
+
+ ASSERT_EQ(ip_data.size(), num_read);
+
+ auto const ip_data_expected = this->getIPData();
+
+ for (std::size_t i = 0; i < ip_data_expected.size(); ++i)
+ {
+ EXPECT_DOUBLE_EQ(ip_data_expected[i].scalar, ip_data[i].scalar)
+ << "Values at integration point " << i << " differ.";
+ }
+}
+
+TYPED_TEST(ProcessLib_IPDataAccess, SetKelvinVectorData)
+{
+ constexpr int dim = TypeParam::value;
+
+ auto ip_data = this->getIPDataNaNs();
+
+ auto const cache = this->getKVDataTransposedOrder();
+
+ auto const num_read = ProcessLib::setIntegrationPointKelvinVectorData<dim>(
+ &cache[0], ip_data, &IPData<dim>::kelvin);
+
+ ASSERT_EQ(ip_data.size(), num_read);
+
+ auto const ip_data_expected = this->getIPData();
+
+ for (std::size_t i = 0; i < ip_data_expected.size(); ++i)
+ {
+ EXPECT_THAT(
+ ip_data[i].kelvin,
+ testing::Pointwise(testing::DoubleEq(), ip_data_expected[i].kelvin))
+ << "Values at integration point " << i << " differ.";
+ }
+}