diff --git a/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp b/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp index f19707a2dc2..724f100c1ce 100644 --- a/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp +++ b/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp @@ -18,6 +18,7 @@ #include "MeshLib/Utils/getOrCreateMeshProperty.h" #include "NumLib/DOF/DOFTableUtil.h" #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h" +#include "ProcessLib/Output/CellAverageAlgorithm.h" #include "ProcessLib/Process.h" #include "ProcessLib/Reflection/ReflectionForExtrapolation.h" #include "ProcessLib/Reflection/ReflectionForIPWriters.h" @@ -207,6 +208,8 @@ void LargeDeformationProcess::computeSecondaryVariableConcrete( GlobalExecutor::executeSelectedMemberOnDereferenced( &LocalAssemblerInterface::computeSecondaryVariable, _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, x_prev, process_id); + + computeCellAverages(cell_average_data_, _local_assemblers); } template class LargeDeformationProcess<2>; template class LargeDeformationProcess<3>; diff --git a/ProcessLib/Output/CellAverageAlgorithm.h b/ProcessLib/Output/CellAverageAlgorithm.h new file mode 100644 index 00000000000..1bb82cc45f4 --- /dev/null +++ b/ProcessLib/Output/CellAverageAlgorithm.h @@ -0,0 +1,63 @@ +/** + * \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 "CellAverageData.h" +#include "ProcessLib/Reflection/ReflectionIPData.h" + +namespace ProcessLib +{ +namespace detail +{ +void computeCellAverages(CellAverageData& cell_average_data, + std::string const& name, unsigned const num_comp, + auto&& flattened_ip_data_accessor, + auto const& local_assemblers) +{ + auto& prop_vec = + cell_average_data.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(ip_data.size() / num_comp); + Eigen::Map ip_data_mapped{ip_data.data(), + num_comp, num_ips}; + + Eigen::Map{&prop_vec[i * num_comp], num_comp} = + ip_data_mapped.rowwise().mean(); + } +} +} // namespace detail + +template +void computeCellAverages( + CellAverageData& cell_average_data, + std::vector> const& local_assemblers) +{ + auto const callback = [&cell_average_data, &local_assemblers]( + std::string const& name, + unsigned const num_comp, + auto&& flattened_ip_data_accessor) + { + detail::computeCellAverages(cell_average_data, name, num_comp, + flattened_ip_data_accessor, + local_assemblers); + }; + + ProcessLib::Reflection::forEachReflectedFlattenedIPDataAccessor( + LAIntf::getReflectionDataForOutput(), callback); +} +} // namespace ProcessLib 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& 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(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( + const_cast(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..19dda7699d1 --- /dev/null +++ b/ProcessLib/Output/CellAverageData.h @@ -0,0 +1,29 @@ +/** + * \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" + +namespace ProcessLib +{ +struct CellAverageData +{ + explicit CellAverageData(MeshLib::Mesh& mesh) : mesh_{mesh} {} + + MeshLib::PropertyVector& getOrCreatePropertyVector( + std::string const& name, unsigned const num_comp); + +private: + MeshLib::Mesh const& mesh_; + std::map*> cell_averages_; +}; +} // namespace ProcessLib diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp index 98ed6566750..9e9277e0169 100644 --- a/ProcessLib/Process.cpp +++ b/ProcessLib/Process.cpp @@ -54,6 +54,7 @@ Process::Process( : name(std::move(name_)), _mesh(mesh), _secondary_variables(std::move(secondary_variables)), + cell_average_data_(mesh), _jacobian_assembler(std::move(jacobian_assembler)), _global_assembler(*_jacobian_assembler), _use_monolithic_scheme(use_monolithic_scheme), diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h index 1816b836668..be39d2bcb13 100644 --- a/ProcessLib/Process.h +++ b/ProcessLib/Process.h @@ -21,6 +21,7 @@ #include "ParameterLib/Parameter.h" #include "ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.h" #include "ProcessLib/BoundaryConditionAndSourceTerm/SourceTermCollection.h" +#include "ProcessLib/Output/CellAverageData.h" #include "ProcessLib/Output/ExtrapolatorData.h" #include "ProcessLib/Output/SecondaryVariable.h" #include "ProcessVariable.h" @@ -361,6 +362,8 @@ class Process SecondaryVariableCollection _secondary_variables; + CellAverageData cell_average_data_; + // TODO (CL) switch to the parallel vector matrix assembler here, once // feature complete, or use the assembly mixin and remove these two members. std::unique_ptr _jacobian_assembler; diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp index 2b6223c5524..d4147358251 100644 --- a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp +++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp @@ -15,6 +15,7 @@ #include "MeshLib/Utils/IntegrationPointWriter.h" #include "MeshLib/Utils/getOrCreateMeshProperty.h" #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h" +#include "ProcessLib/Output/CellAverageAlgorithm.h" #include "ProcessLib/Process.h" #include "ProcessLib/Reflection/ReflectionForExtrapolation.h" #include "ProcessLib/Reflection/ReflectionForIPWriters.h" @@ -216,6 +217,8 @@ void SmallDeformationProcess::computeSecondaryVariableConcrete( GlobalExecutor::executeSelectedMemberOnDereferenced( &LocalAssemblerInterface::computeSecondaryVariable, _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, x_prev, process_id); + + computeCellAverages(cell_average_data_, _local_assemblers); } template class SmallDeformationProcess<2>; template class SmallDeformationProcess<3>; diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h index 154248a716b..42d127276f9 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -1910,23 +1910,10 @@ void TH2MLocalAssemblerprocess_data_.temperature_interpolated); - unsigned const n_integration_points = - this->integration_method_.getNumberOfPoints(); - - double saturation_avg = 0; - ConstitutiveRelations::ConstitutiveModels const models{ this->solid_material_, *this->process_data_.phase_transition_model_}; updateConstitutiveVariables(local_x, local_x_prev, t, dt, models); - - for (unsigned ip = 0; ip < n_integration_points; ip++) - { - saturation_avg += this->current_states_[ip].S_L_data.S_L; - } - saturation_avg /= n_integration_points; - (*this->process_data_.element_saturation)[this->element_.getID()] = - saturation_avg; } } // namespace TH2M diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp index b3289ca9d14..376369b107e 100644 --- a/ProcessLib/TH2M/TH2MProcess.cpp +++ b/ProcessLib/TH2M/TH2MProcess.cpp @@ -19,6 +19,7 @@ #include "MeshLib/Utils/getOrCreateMeshProperty.h" #include "NumLib/DOF/DOFTableUtil.h" #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h" +#include "ProcessLib/Output/CellAverageAlgorithm.h" #include "ProcessLib/Process.h" #include "ProcessLib/Reflection/ReflectionForExtrapolation.h" #include "ProcessLib/Reflection/ReflectionForIPWriters.h" @@ -175,10 +176,6 @@ void TH2MProcess::initializeConcreteProcess( _process_data.solid_materials, local_assemblers_, _integration_point_writer, integration_order); - _process_data.element_saturation = MeshLib::getOrCreateMeshProperty( - const_cast(mesh), "saturation_avg", - MeshLib::MeshItemType::Cell, 1); - _process_data.gas_pressure_interpolated = MeshLib::getOrCreateMeshProperty( const_cast(mesh), "gas_pressure_interpolated", @@ -320,6 +317,8 @@ void TH2MProcess::computeSecondaryVariableConcrete( &LocalAssemblerInterface::computeSecondaryVariable, local_assemblers_, pv.getActiveElementIDs(), getDOFTables(x.size()), t, dt, x, x_prev, process_id); + + computeCellAverages(cell_average_data_, local_assemblers_); } template diff --git a/ProcessLib/TH2M/TH2MProcessData.h b/ProcessLib/TH2M/TH2MProcessData.h index 8aa138d9758..e23f42d8cff 100644 --- a/ProcessLib/TH2M/TH2MProcessData.h +++ b/ProcessLib/TH2M/TH2MProcessData.h @@ -52,7 +52,6 @@ struct TH2MProcessData const bool use_TaylorHood_elements; - MeshLib::PropertyVector* element_saturation = nullptr; MeshLib::PropertyVector* gas_pressure_interpolated = nullptr; MeshLib::PropertyVector* capillary_pressure_interpolated = nullptr; MeshLib::PropertyVector* temperature_interpolated = nullptr; diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h index f76508b2bc1..3a9c9776320 100644 --- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h +++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h @@ -541,14 +541,6 @@ void ThermoRichardsMechanicsLocalAssemblerintegration_method_.getNumberOfPoints(); - double saturation_avg = 0; - double porosity_avg = 0; - double liquid_density_avg = 0; - double viscosity_avg = 0; - - using KV = MathLib::KelvinVector::KelvinVectorType; - KV sigma_avg = KV::Zero(); - typename ConstitutiveTraits::ConstitutiveSetting constitutive_setting; auto models = ConstitutiveTraits::createConstitutiveModels( @@ -603,29 +595,7 @@ void ThermoRichardsMechanicsLocalAssemblerprev_states_[ip], this->material_states_[ip], tmp, output_data, CD); - - saturation_avg += std::get(current_state).S_L; - porosity_avg += std::get(current_state).phi; - - liquid_density_avg += std::get(output_data).rho_LR; - viscosity_avg += std::get(output_data).viscosity; - sigma_avg += ConstitutiveTraits::ConstitutiveSetting::statefulStress( - current_state); } - saturation_avg /= n_integration_points; - porosity_avg /= n_integration_points; - viscosity_avg /= n_integration_points; - liquid_density_avg /= n_integration_points; - sigma_avg /= n_integration_points; - - (*process_data.element_saturation)[e_id] = saturation_avg; - (*process_data.element_porosity)[e_id] = porosity_avg; - (*process_data.element_liquid_density)[e_id] = liquid_density_avg; - (*process_data.element_viscosity)[e_id] = viscosity_avg; - - Eigen::Map( - &(*process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) = - MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_avg); NumLib::interpolateToHigherOrderNodes< ShapeFunction, typename ShapeFunctionDisplacement::MeshElement, diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp index b2961f64854..f00e5791de9 100644 --- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp +++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp @@ -18,6 +18,7 @@ #include "MeshLib/Utils/getOrCreateMeshProperty.h" #include "NumLib/DOF/DOFTableUtil.h" #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h" +#include "ProcessLib/Output/CellAverageAlgorithm.h" #include "ProcessLib/Process.h" #include "ProcessLib/Reflection/ReflectionForExtrapolation.h" #include "ProcessLib/Reflection/ReflectionForIPWriters.h" @@ -156,29 +157,6 @@ void ThermoRichardsMechanicsProcess:: process_data_.solid_materials, local_assemblers_, _integration_point_writer, integration_order); - process_data_.element_saturation = MeshLib::getOrCreateMeshProperty( - const_cast(mesh), "saturation_avg", - MeshLib::MeshItemType::Cell, 1); - - process_data_.element_porosity = MeshLib::getOrCreateMeshProperty( - const_cast(mesh), "porosity_avg", - MeshLib::MeshItemType::Cell, 1); - - process_data_.element_liquid_density = - MeshLib::getOrCreateMeshProperty( - const_cast(mesh), "liquid_density_avg", - MeshLib::MeshItemType::Cell, 1); - - process_data_.element_viscosity = MeshLib::getOrCreateMeshProperty( - const_cast(mesh), "viscosity_avg", - MeshLib::MeshItemType::Cell, 1); - - process_data_.element_stresses = MeshLib::getOrCreateMeshProperty( - const_cast(mesh), "stress_avg", - MeshLib::MeshItemType::Cell, - MathLib::KelvinVector::KelvinVectorType< - DisplacementDim>::RowsAtCompileTime); - process_data_.pressure_interpolated = MeshLib::getOrCreateMeshProperty( const_cast(mesh), "pressure_interpolated", @@ -309,6 +287,8 @@ void ThermoRichardsMechanicsProcess:: &LocalAssemblerIF::computeSecondaryVariable, local_assemblers_, pv.getActiveElementIDs(), getDOFTables(x.size()), t, dt, x, x_prev, process_id); + + computeCellAverages(cell_average_data_, local_assemblers_); } template diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h index b7933599c8f..b01f9d61280 100644 --- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h +++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h @@ -55,11 +55,6 @@ struct ThermoRichardsMechanicsProcessData InitializePorosityFromMediumProperty const initialize_porosity_from_medium_property; - MeshLib::PropertyVector* element_saturation = nullptr; - MeshLib::PropertyVector* element_porosity = nullptr; - MeshLib::PropertyVector* element_liquid_density = nullptr; - MeshLib::PropertyVector* element_viscosity = nullptr; - MeshLib::PropertyVector* element_stresses = nullptr; MeshLib::PropertyVector* temperature_interpolated = nullptr; MeshLib::PropertyVector* pressure_interpolated = nullptr; diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small.prj b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small.prj index 02b646e0280..ae7609e2109 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small.prj +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small.prj @@ -205,6 +205,7 @@ NodalForces pressure_interpolated saturation_avg + sigma_avg _ts_{:timestep}_t_{:time} @@ -354,6 +355,12 @@ 1e-8 0 + + RichardsFlow_2d_small_ts_.*.vtu + sigma_avg + 2.8e-12 + 0 + RichardsFlow_2d_small_ts_.*.vtu epsilon diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_0_t_0.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_0_t_0.000000.vtu index ce7d33c10c9..b9dd3ced866 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_0_t_0.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_0_t_0.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_104_t_2000.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_104_t_2000.000000.vtu index 91e4affb979..784c259c3e4 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_104_t_2000.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_104_t_2000.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_19_t_318.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_19_t_318.000000.vtu index a9f53bc77b8..339fa819ca5 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_19_t_318.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_19_t_318.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_1_t_1.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_1_t_1.000000.vtu index 703fec56f30..83ab412ef92 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_1_t_1.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_1_t_1.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_29_t_518.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_29_t_518.000000.vtu index b6202102d14..8cd523c95c1 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_29_t_518.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_29_t_518.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_2_t_3.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_2_t_3.000000.vtu index 264abfd4282..8d6a382bdfa 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_2_t_3.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_2_t_3.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_39_t_718.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_39_t_718.000000.vtu index cce79dd592d..39c90e840f4 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_39_t_718.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_39_t_718.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_3_t_8.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_3_t_8.000000.vtu index 47475355ff2..03367285fba 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_3_t_8.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_3_t_8.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_49_t_918.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_49_t_918.000000.vtu index 6f67c0b5486..19d052c3012 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_49_t_918.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_49_t_918.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_4_t_18.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_4_t_18.000000.vtu index a993960c709..ea244a19a65 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_4_t_18.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_4_t_18.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_59_t_1118.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_59_t_1118.000000.vtu index 6a6c5a26731..22d5a98f6bd 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_59_t_1118.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_59_t_1118.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_5_t_38.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_5_t_38.000000.vtu index 89bef788e68..4ee3963156d 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_5_t_38.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_5_t_38.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_69_t_1318.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_69_t_1318.000000.vtu index 38a1c9d9cd6..879a21c0feb 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_69_t_1318.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_69_t_1318.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_6_t_58.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_6_t_58.000000.vtu index 99d0546d156..0848e3a65dd 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_6_t_58.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_6_t_58.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_79_t_1518.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_79_t_1518.000000.vtu index f27991cbb38..f4f39ef92a8 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_79_t_1518.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_79_t_1518.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_7_t_78.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_7_t_78.000000.vtu index 287a168431e..017f00e2fab 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_7_t_78.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_7_t_78.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_89_t_1718.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_89_t_1718.000000.vtu index e71fddbc39c..d012b085e7a 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_89_t_1718.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_89_t_1718.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_8_t_98.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_8_t_98.000000.vtu index 0b00240c916..bd511e3f761 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_8_t_98.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_8_t_98.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_99_t_1918.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_99_t_1918.000000.vtu index 69eacf9c7f6..df0b16698e0 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_99_t_1918.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_99_t_1918.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_9_t_118.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_9_t_118.000000.vtu index 8add72d94db..47c3c4b5601 100644 --- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_9_t_118.000000.vtu +++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_9_t_118.000000.vtu @@ -27,7 +27,7 @@ - + diff --git a/Tests/ProcessLib/TestReflectIPData.cpp b/Tests/ProcessLib/TestReflectIPData.cpp index 5ace325219c..ad6fbe060b6 100644 --- a/Tests/ProcessLib/TestReflectIPData.cpp +++ b/Tests/ProcessLib/TestReflectIPData.cpp @@ -14,6 +14,9 @@ #include #include +#include "MeshToolsLib/MeshGenerators/MeshGenerator.h" +#include "ProcessLib/Output/CellAverageAlgorithm.h" +#include "ProcessLib/Reflection/ReflectionForIPWriters.h" #include "ProcessLib/Reflection/ReflectionIPData.h" #include "ProcessLib/Reflection/ReflectionSetIPData.h" @@ -130,6 +133,8 @@ struct LocAsmIF makeReflectionData(&LocAsmIF::ip_data_level1), makeReflectionData(&LocAsmIF::ip_data_level1b)}; } + + static auto getReflectionDataForOutput() { return reflect(); } }; template @@ -251,6 +256,7 @@ std::vector initKelvin(LocAsmIF& loc_asm, for (std::size_t ip = 0; ip < num_int_pts; ++ip) { ip_data_accessor(loc_asm, ip) = + // the local assembler stores Kelvin vector data... MathLib::KelvinVector::symmetricTensorToKelvinVector( Eigen::Vector::LinSpaced( kv_size, ip * kv_size + start_value, @@ -268,6 +274,7 @@ std::vector initKelvin(LocAsmIF& loc_asm, } // prepare reference data + // ... the reference data used in the tests is not Kelvin-mapped! std::vector vector_expected(num_int_pts * kv_size); iota(begin(vector_expected), end(vector_expected), start_value); return vector_expected; @@ -329,6 +336,10 @@ std::vector initMatrix(LocAsmIF& loc_asm, template struct ReferenceData { +private: + ReferenceData() = default; + +public: std::vector scalar; std::vector vector; std::vector kelvin; @@ -348,6 +359,8 @@ struct ReferenceData std::vector scalar2b; std::vector scalar3b; + // Computes reference data and initializes the internal (integration point) + // data of the passed \c loc_asm. static ReferenceData create(LocAsmIF& loc_asm, bool const for_read_test) { @@ -484,6 +497,45 @@ struct ReferenceData } }; +template +void checkAll(ReferenceData const& ref, auto&& checker, + bool const check_level_0_data = true) +{ + auto constexpr kv_size = + MathLib::KelvinVector::kelvin_vector_dimensions(dim); + + // The storage of "level 0 data" (i.e., non-reflected data directly in the + // local assembler, e.g. a std::vector of integration point data of + // a scalar value) is not compatible with the implemented IP data input + // logic. Therefore we don't test level 0 data in that case. + // This incompatibility is not a problem in practice, as such data is not + // used ATM in OGS's local assemblers. + if (check_level_0_data) + { + checker("scalar", 1, ref.scalar); + checker("vector", dim, ref.vector); + checker("kelvin", kv_size, ref.kelvin); + } + + // level 1 + checker("scalar1", 1, ref.scalar1); + checker("vector1", dim, ref.vector1); + checker("kelvin1", kv_size, ref.kelvin1); + + // level 3 + checker("scalar3", 1, ref.scalar3); + checker("vector3", dim, ref.vector3); + checker("kelvin3", kv_size, ref.kelvin3); + checker("matrix3", dim * 4, ref.matrix3); + checker("matrix3_1", dim * 4, ref.matrix3_1); + checker("matrix3_2", 4, ref.matrix3_2); + + // b levels + checker("scalar1b", 1, ref.scalar1b); + checker("scalar2b", 1, ref.scalar2b); + checker("scalar3b", 1, ref.scalar3b); +} + template struct ProcessLib_ReflectIPData : ::testing::Test { @@ -499,8 +551,6 @@ 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; @@ -548,28 +598,7 @@ TYPED_TEST(ProcessLib_ReflectIPData, ReadTest) << "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); - check("matrix3", dim * 4, ref.matrix3); - check("matrix3_1", dim * 4, ref.matrix3_1); - check("matrix3_2", 4, ref.matrix3_2); - - // b levels - check("scalar1b", 1, ref.scalar1b); - check("scalar2b", 1, ref.scalar2b); - check("scalar3b", 1, ref.scalar3b); + checkAll(ref, check); } TYPED_TEST(ProcessLib_ReflectIPData, WriteTest) @@ -583,7 +612,8 @@ TYPED_TEST(ProcessLib_ReflectIPData, WriteTest) auto const ref = ReferenceData::create(loc_asm, false); - // data getters - used for checks ////////////////////////////////////////// + // set up data getters - used for checks /////////////////////////////////// + // that critically relies on the read test above being successful! std::map> map_name_to_num_comp_and_function; @@ -603,7 +633,7 @@ TYPED_TEST(ProcessLib_ReflectIPData, WriteTest) // checks ////////////////////////////////////////////////////////////////// auto check = [&map_name_to_num_comp_and_function, &loc_asm]( - std::string const& name, std::size_t const size_expected, + std::string const& name, unsigned const num_comp, std::vector const& values_plain) { auto const it = map_name_to_num_comp_and_function.find(name); @@ -611,19 +641,34 @@ TYPED_TEST(ProcessLib_ReflectIPData, WriteTest) 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; + auto const& [num_comp_2, fct] = it->second; + + // consistency checks + ASSERT_EQ(num_comp, num_comp_2); + ASSERT_EQ(0, values_plain.size() % num_comp); + auto const num_int_pts_expected = values_plain.size() / num_comp; EXPECT_THAT(fct(loc_asm), testing::Each(testing::IsNan())) - << "All values must be initialize to NaN in this unit test. Check " + << "All values must be initialized to NaN in this unit test. Check " "failed for ip data with name '" << name << "'"; - auto const size = ProcessLib::Reflection::reflectSetIPData( - name, values_plain.data(), loc_asm.ip_data_level1); + // function under test ///////////////////////////////////////////////// - EXPECT_EQ(size_expected, size) - << "Unexpected size obtained for ip data with name '" << name - << "'"; + auto const num_int_pts_actual = + ProcessLib::Reflection::reflectSetIPData( + name, values_plain.data(), loc_asm.ip_data_level1); + auto const num_int_pts_actual_1 = + ProcessLib::Reflection::reflectSetIPData( + name, values_plain.data(), loc_asm.ip_data_level1b); + + // end function under test ///////////////////////////////////////////// + + ASSERT_EQ(num_int_pts_expected, num_int_pts_actual) + << "Unexpected number of integration points obtained for ip data " + "with name '" + << name << "'"; + ASSERT_EQ(num_int_pts_expected, num_int_pts_actual_1); // check set values via round-trip with getter tested in previous unit // test @@ -633,19 +678,7 @@ TYPED_TEST(ProcessLib_ReflectIPData, WriteTest) << "'"; }; - check("scalar1", num_int_pts, ref.scalar1); - check("vector1", num_int_pts, ref.vector1); - check("kelvin1", num_int_pts, ref.kelvin1); - - check("scalar3", num_int_pts, ref.scalar3); - check("vector3", num_int_pts, ref.vector3); - check("kelvin3", num_int_pts, ref.kelvin3); - check("matrix3", num_int_pts, ref.matrix3); - check("matrix3_1", num_int_pts, ref.matrix3_1); - check("matrix3_2", num_int_pts, ref.matrix3_2); - - check("scalar2b", num_int_pts, ref.scalar2b); - check("scalar3b", num_int_pts, ref.scalar3b); + checkAll(ref, check, false); } TYPED_TEST(ProcessLib_ReflectIPData, RawDataTypes) @@ -687,3 +720,147 @@ TYPED_TEST(ProcessLib_ReflectIPData, RawDataTypes) static_assert( !PRD::is_raw_data>::value); } + +TYPED_TEST(ProcessLib_ReflectIPData, IPWriterTest) +{ + constexpr int dim = TypeParam::value; + + using LocAsm = LocAsmIF; + + std::size_t const num_int_pts = 8; + std::vector> loc_asms; + // emplace...new is a workaround for an MSVC compiler warning: + // C:\Program Files (x86)\Microsoft Visual + // Studio\2019\Community\VC\Tools\MSVC\14.29.30133\include\memory(3382): + // warning C4244: 'argument': conversion from 'const uintptr_t' to 'const + // unsigned int', possible loss of data + // C:/Users/gitlab/gitlab/_b/bg4d5s_d/0/ogs/ogs/Tests/ProcessLib/ + // TestReflectIPData.cpp(792): + // note: see reference to function template instantiation + // 'std::unique_ptr> + // std::make_unique(const size_t &)' being compiled + loc_asms.emplace_back(new LocAsm{num_int_pts}); + + std::unique_ptr mesh{ + MeshToolsLib::MeshGenerator::generateLineMesh(1.0, 1)}; + + auto const ref = ReferenceData::create(*loc_asms.front(), true); + + // function under test ///////////////////////////////////////////////////// + + std::vector> + integration_point_writers; + constexpr int integration_order = 0xc0ffee; // dummy value + ProcessLib::Reflection::addReflectedIntegrationPointWriters( + LocAsm::getReflectionDataForOutput(), integration_point_writers, + integration_order, loc_asms); + + // this finally invokes the IP writers + MeshLib::addIntegrationPointDataToMesh(*mesh, integration_point_writers); + + // checks ////////////////////////////////////////////////////////////////// + + auto check = [&mesh](std::string const& name, + unsigned const num_comp, + std::vector const& values_expected) + { + auto const& props = mesh->getProperties(); + + ASSERT_TRUE(props.existsPropertyVector(name + "_ip")) + << "Property vector '" << name + "_ip" + << "' does not exist in the mesh."; + + ASSERT_EQ(0, values_expected.size() % num_comp); + auto const num_int_pts = values_expected.size() / num_comp; + + auto const& ip_data_actual = + *props.getPropertyVector(name + "_ip"); + + ASSERT_EQ(MeshLib::MeshItemType::IntegrationPoint, + ip_data_actual.getMeshItemType()); + ASSERT_EQ(mesh->getNumberOfElements() * num_int_pts, + ip_data_actual.getNumberOfTuples()); + ASSERT_EQ(num_comp, ip_data_actual.getNumberOfGlobalComponents()); + + EXPECT_THAT(ip_data_actual, + testing::Pointwise(testing::DoubleEq(), values_expected)) + << "Values differ for ip data with name '" << name << "'"; + }; + + checkAll(ref, check); +} + +TYPED_TEST(ProcessLib_ReflectIPData, CellAverageTest) +{ + constexpr int dim = TypeParam::value; + + using LocAsm = LocAsmIF; + + std::size_t const num_int_pts = 8; + std::vector> loc_asms; + // emplace...new is a workaround for an MSVC compiler warning: + // C:\Program Files (x86)\Microsoft Visual + // Studio\2019\Community\VC\Tools\MSVC\14.29.30133\include\memory(3382): + // warning C4244: 'argument': conversion from 'const uintptr_t' to 'const + // unsigned int', possible loss of data + // C:/Users/gitlab/gitlab/_b/bg4d5s_d/0/ogs/ogs/Tests/ProcessLib/ + // TestReflectIPData.cpp(792): + // note: see reference to function template instantiation + // 'std::unique_ptr> + // std::make_unique(const size_t &)' being compiled + loc_asms.emplace_back(new LocAsm{num_int_pts}); + + std::unique_ptr mesh{ + MeshToolsLib::MeshGenerator::generateLineMesh(1.0, 1)}; + + auto const ref = ReferenceData::create(*loc_asms.front(), true); + + // function under test ///////////////////////////////////////////////////// + + ProcessLib::CellAverageData cell_average_data{*mesh}; + computeCellAverages(cell_average_data, loc_asms); + + // checks ////////////////////////////////////////////////////////////////// + + auto check = [&mesh](std::string const& name, + unsigned const num_comp_expected, + std::vector const& ip_data) + { + ASSERT_EQ(num_int_pts * num_comp_expected, ip_data.size()); + + // TODO the implementation here in the unit test computing cell average + // reference data might be too similar to the production code. In fact, + // it's almost the same code. + + // each integration point corresponds to a column in the mapped + // matrix, vector components/matrix entries are stored contiguously + // in memory. + Eigen::Map> + ip_data_mat{ip_data.data(), num_comp_expected, num_int_pts}; + + Eigen::VectorXd const cell_avg_expected = ip_data_mat.rowwise().mean(); + + auto const& props = mesh->getProperties(); + + ASSERT_TRUE(props.existsPropertyVector(name + "_avg")) + << "Property vector '" << name + "_avg" + << "' does not exist in the mesh."; + + auto const& cell_avg_actual = + *props.getPropertyVector(name + "_avg"); + + ASSERT_EQ(MeshLib::MeshItemType::Cell, + cell_avg_actual.getMeshItemType()); + ASSERT_EQ(mesh->getNumberOfElements(), + cell_avg_actual.getNumberOfTuples()); + ASSERT_EQ(num_comp_expected, + cell_avg_actual.getNumberOfGlobalComponents()); + + EXPECT_THAT(cell_avg_actual, + testing::Pointwise(testing::DoubleEq(), cell_avg_expected)) + << "Values differ for cell average data with name '" << name << "'"; + }; + + checkAll(ref, check); +}