From 351b7a760ef08666d1ed34eee7d09030974919ea Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 08:26:15 +0200 Subject: [PATCH 01/14] [PL/RM] Moved sigma_sw and eps_m to stateful data --- .../ConstitutiveRelations/ConstitutiveData.h | 6 +- .../RichardsMechanics/IntegrationPointData.h | 25 +- .../RichardsMechanicsFEM-impl.h | 244 +++++++++++++----- 3 files changed, 199 insertions(+), 76 deletions(-) diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h index 7783b02fc7d..d32cba9ac66 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h @@ -30,7 +30,11 @@ template using StatefulData = std::tuple< StrainData, ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature:: - EffectiveStressData>; + EffectiveStressData, + ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature:: + SwellingDataStateful, + ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature:: + MechanicalStrainData>; template using StatefulDataPrev = ProcessLib::ConstitutiveRelations::PrevStateOf< diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h index 6363455b5f7..189fdcf4b47 100644 --- a/ProcessLib/RichardsMechanics/IntegrationPointData.h +++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h @@ -29,19 +29,8 @@ struct IntegrationPointData final material_state_variables( solid_material.createMaterialStateVariables()) { - // Initialize current time step values - static const int kelvin_vector_size = - MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim); - sigma_sw.setZero(kelvin_vector_size); - eps_m.setZero(kelvin_vector_size); - - // Previous time step values are not initialized and are set later. - eps_m_prev.resize(kelvin_vector_size); } - typename BMatricesType::KelvinVectorType sigma_sw, sigma_sw_prev; - typename BMatricesType::KelvinVectorType eps_m, eps_m_prev; - typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u; typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u; @@ -74,8 +63,6 @@ struct IntegrationPointData final void pushBackState() { - eps_m_prev = eps_m; - sigma_sw_prev = sigma_sw; saturation_prev = saturation; saturation_m_prev = saturation_m; porosity_prev = porosity; @@ -132,13 +119,21 @@ struct IntegrationPointData final DisplacementDim>& sigma_eff, PrevState> const& sigma_eff_prev) + DisplacementDim>> const& sigma_eff_prev, + ProcessLib::ThermoRichardsMechanics:: + ConstitutiveStress_StrainTemperature::MechanicalStrainData< + DisplacementDim> const& + /*eps_m*/, + PrevState< + ProcessLib::ThermoRichardsMechanics:: + ConstitutiveStress_StrainTemperature::MechanicalStrainData< + DisplacementDim>> const& eps_m_prev) { MaterialPropertyLib::VariableArray variable_array_prev; variable_array_prev.stress = sigma_eff_prev->sigma_eff; variable_array_prev.mechanical_strain .emplace>( - eps_m_prev); + eps_m_prev->eps_m); variable_array_prev.temperature = temperature; auto&& solution = solid_material.integrateStress( diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index 13728943eda..e952e68bd78 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -41,19 +41,22 @@ void updateSwellingStressAndVolumetricStrain( double const alpha, double const phi, double const p_cap_ip, MPL::VariableArray& variables, MPL::VariableArray& variables_prev, ParameterLib::SpatialPosition const& x_position, double const t, - double const dt) + double const dt, + ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature:: + SwellingDataStateful& sigma_sw, + PrevState> const& sigma_sw_prev) { auto const& identity2 = MathLib::KelvinVector::Invariants< MathLib::KelvinVector::kelvin_vector_dimensions( DisplacementDim)>::identity2; - auto& sigma_sw = ip_data.sigma_sw; - auto const& sigma_sw_prev = ip_data.sigma_sw_prev; if (!medium.hasProperty(MPL::PropertyType::saturation_micro)) { // If there is swelling, compute it. Update volumetric strain rate, // s.t. it corresponds to the mechanical part only. - sigma_sw = sigma_sw_prev; + sigma_sw = *sigma_sw_prev; if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)) { auto const sigma_sw_dot = @@ -62,14 +65,15 @@ void updateSwellingStressAndVolumetricStrain( solid_phase[MPL::PropertyType::swelling_stress_rate] .value(variables, variables_prev, x_position, t, dt))); - sigma_sw += sigma_sw_dot * dt; + sigma_sw.sigma_sw += sigma_sw_dot * dt; // !!! Misusing volumetric strain for mechanical volumetric // strain just to update the transport porosity !!! variables.volumetric_strain += - identity2.transpose() * C_el.inverse() * sigma_sw; - variables_prev.volumetric_strain += - identity2.transpose() * C_el.inverse() * sigma_sw_prev; + identity2.transpose() * C_el.inverse() * sigma_sw.sigma_sw; + variables_prev.volumetric_strain += identity2.transpose() * + C_el.inverse() * + sigma_sw_prev->sigma_sw; } } @@ -109,7 +113,7 @@ void updateSwellingStressAndVolumetricStrain( medium.property(MPL::PropertyType::saturation_micro) .template value(variables, x_position, t, dt); } - sigma_sw = sigma_sw_prev + delta_sigma_sw; + sigma_sw.sigma_sw = sigma_sw_prev->sigma_sw + delta_sigma_sw; } } @@ -248,7 +252,15 @@ std::size_t RichardsMechanicsLocalAssembler< if (name == "swelling_stress") { return ProcessLib::setIntegrationPointKelvinVectorData( - values, _ip_data, &IpData::sigma_sw); + values, this->current_states_, + [](auto& tuple) -> auto& + { + return std::get>( + tuple) + .sigma_sw; + }); } if (name == "epsilon") { @@ -350,9 +362,20 @@ void RichardsMechanicsLocalAssembler>(this->current_states_[ip]) .eps; - auto& sigma_sw = _ip_data[ip].sigma_sw; + auto const& sigma_sw = + std::get>( + this->current_states_[ip]) + .sigma_sw; + auto& eps_m_prev = + std::get>>( + this->prev_states_[ip]) + ->eps_m; - _ip_data[ip].eps_m_prev.noalias() = + eps_m_prev.noalias() = solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate) ? eps + C_el.inverse() * sigma_sw : eps; @@ -450,7 +473,6 @@ void RichardsMechanicsLocalAssembler< auto& eps = std::get>(this->current_states_[ip]); - auto& eps_m = _ip_data[ip].eps_m; eps.eps.noalias() = B * u; auto& S_L = _ip_data[ip].saturation; @@ -546,8 +568,17 @@ void RichardsMechanicsLocalAssembler< // Swelling and possibly volumetric strain rate update. { - auto& sigma_sw = _ip_data[ip].sigma_sw; - auto const& sigma_sw_prev = _ip_data[ip].sigma_sw_prev; + auto& sigma_sw = + std::get>( + this->current_states_[ip]) + .sigma_sw; + auto const& sigma_sw_prev = std::get>>(this->prev_states_[ip]) + ->sigma_sw; // If there is swelling, compute it. Update volumetric strain rate, // s.t. it corresponds to the mechanical part only. @@ -595,7 +626,12 @@ void RichardsMechanicsLocalAssembler< liquid_phase.property(MPL::PropertyType::viscosity) .template value(variables, x_position, t, dt); - auto const& sigma_sw = _ip_data[ip].sigma_sw; + auto const& sigma_sw = + std::get>( + this->current_states_[ip]) + .sigma_sw; auto const& sigma_eff = std::get>( + { + auto& eps_m = + std::get>( + this->current_states_[ip]) + .eps_m; + eps_m.noalias() = + solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate) + ? eps.eps + C_el.inverse() * sigma_sw + : eps.eps; + variables.mechanical_strain.emplace< + MathLib::KelvinVector::KelvinVectorType>( eps_m); + } { auto& SD = this->current_states_[ip]; @@ -648,10 +692,19 @@ void RichardsMechanicsLocalAssembler< ConstitutiveStress_StrainTemperature:: EffectiveStressData>>( SD_prev); + auto const& eps_m = + std::get>(SD); + auto& eps_m_prev = std::get< + PrevState>>( + SD_prev); - _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, - dt, temperature, sigma_eff, - sigma_eff_prev); + _ip_data[ip].updateConstitutiveRelation( + variables, t, x_position, dt, temperature, sigma_eff, + sigma_eff_prev, eps_m, eps_m_prev); } // p_SR @@ -747,7 +800,6 @@ void RichardsMechanicsLocalAssembler>(SD); - auto& eps_m = ip_data.eps_m; auto& S_L = ip_data.saturation; auto const S_L_prev = ip_data.saturation_prev; auto const alpha = @@ -849,11 +901,24 @@ void RichardsMechanicsLocalAssembler(CD) = mu; - // Swelling and possibly volumetric strain rate update. - updateSwellingStressAndVolumetricStrain( - ip_data, *medium, solid_phase, C_el, rho_LR, mu, - this->process_data_.micro_porosity_parameters, alpha, phi, p_cap_ip, - variables, variables_prev, x_position, t, dt); + { + // Swelling and possibly volumetric strain rate update. + auto& sigma_sw = + std::get>(SD); + auto const& sigma_sw_prev = + std::get>>( + SD_prev); + + updateSwellingStressAndVolumetricStrain( + ip_data, *medium, solid_phase, C_el, rho_LR, mu, + this->process_data_.micro_porosity_parameters, alpha, phi, p_cap_ip, + variables, variables_prev, x_position, t, dt, sigma_sw, + sigma_sw_prev); + } if (medium->hasProperty(MPL::PropertyType::transport_porosity)) { @@ -912,15 +977,27 @@ void RichardsMechanicsLocalAssembler>( - eps_m); + { + auto& sigma_sw = + std::get>(SD) + .sigma_sw; + + auto& eps_m = + std::get>(SD) + .eps_m; + eps_m.noalias() = + solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate) + ? eps.eps + C_el.inverse() * sigma_sw + : eps.eps; + variables.mechanical_strain + .emplace>( + eps_m); + } { auto& sigma_eff = @@ -932,10 +1009,19 @@ void RichardsMechanicsLocalAssembler>>( SD_prev); + auto const& eps_m = + std::get>(SD); + auto& eps_m_prev = + std::get>>( + SD_prev); - auto C = ip_data.updateConstitutiveRelation(variables, t, x_position, - dt, temperature, sigma_eff, - sigma_eff_prev); + auto C = ip_data.updateConstitutiveRelation( + variables, t, x_position, dt, temperature, sigma_eff, + sigma_eff_prev, eps_m, eps_m_prev); *std::get>(CD) = std::move(C); } @@ -1430,7 +1516,12 @@ std::vector const& RichardsMechanicsLocalAssembler< for (unsigned ip = 0; ip < n_integration_points; ++ip) { - auto const& sigma_sw = _ip_data[ip].sigma_sw; + auto const& sigma_sw = + std::get>( + this->current_states_[ip]) + .sigma_sw; cache_mat.col(ip) = MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_sw); } @@ -1795,7 +1886,6 @@ void RichardsMechanicsLocalAssembler>(this->current_states_[ip]) .eps; eps.noalias() = B * u; - auto& eps_m = _ip_data[ip].eps_m; auto& S_L = _ip_data[ip].saturation; auto const S_L_prev = _ip_data[ip].saturation_prev; S_L = medium->property(MPL::PropertyType::saturation) @@ -1849,11 +1939,25 @@ void RichardsMechanicsLocalAssembler(variables, x_position, t, dt); - // Swelling and possibly volumetric strain rate update. - updateSwellingStressAndVolumetricStrain( - _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu, - this->process_data_.micro_porosity_parameters, alpha, phi, p_cap_ip, - variables, variables_prev, x_position, t, dt); + { + // Swelling and possibly volumetric strain rate update. + auto& sigma_sw = + std::get>( + this->current_states_[ip]); + auto const& sigma_sw_prev = std::get< + PrevState>>( + this->prev_states_[ip]); + + updateSwellingStressAndVolumetricStrain( + _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu, + this->process_data_.micro_porosity_parameters, alpha, phi, + p_cap_ip, variables, variables_prev, x_position, t, dt, + sigma_sw, sigma_sw_prev); + } if (medium->hasProperty(MPL::PropertyType::transport_porosity)) { @@ -1914,15 +2018,26 @@ void RichardsMechanicsLocalAssembler(variables, x_position, t, dt); _ip_data[ip].dry_density_solid = (1 - phi) * rho_SR; - auto& sigma_sw = _ip_data[ip].sigma_sw; - - eps_m.noalias() = - solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate) - ? eps + C_el.inverse() * sigma_sw - : eps; - variables.mechanical_strain - .emplace>( + { + auto& SD = this->current_states_[ip]; + auto const& sigma_sw = + std::get>(SD) + .sigma_sw; + auto& eps_m = + std::get>(SD) + .eps_m; + eps_m.noalias() = + solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate) + ? eps + C_el.inverse() * sigma_sw + : eps; + variables.mechanical_strain.emplace< + MathLib::KelvinVector::KelvinVectorType>( eps_m); + } { auto& SD = this->current_states_[ip]; @@ -1936,10 +2051,19 @@ void RichardsMechanicsLocalAssembler>>( SD_prev); + auto const& eps_m = + std::get>(SD); + auto const& eps_m_prev = std::get< + PrevState>>( + SD_prev); - _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, - dt, temperature, sigma_eff, - sigma_eff_prev); + _ip_data[ip].updateConstitutiveRelation( + variables, t, x_position, dt, temperature, sigma_eff, + sigma_eff_prev, eps_m, eps_m_prev); } auto const& b = this->process_data_.specific_body_force; From e9702a2e1b8166d83e2df32ee2eef1190dd77161 Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 08:51:13 +0200 Subject: [PATCH 02/14] [PL/RM] Made method static ... to exclude access to data members via the this pointer. --- .../RichardsMechanicsFEM-impl.h | 21 ++++++++++++------- .../RichardsMechanics/RichardsMechanicsFEM.h | 6 ++++-- 2 files changed, 18 insertions(+), 9 deletions(-) diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index e952e68bd78..ddf0947c7ab 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -786,7 +786,8 @@ void RichardsMechanicsLocalAssembler const& p_cap_data, ConstitutiveData& CD, StatefulData& SD, - StatefulDataPrev const& SD_prev) + StatefulDataPrev const& SD_prev, + std::optional const& micro_porosity_parameters) { auto const& liquid_phase = medium->phase("AqueousLiquid"); auto const& solid_phase = medium->phase("Solid"); @@ -889,11 +890,18 @@ void RichardsMechanicsLocalAssembler(*x_position.getElementID()) + : static_cast(-1); + auto const ip = + x_position.getIntegrationPoint() + ? static_cast(*x_position.getIntegrationPoint()) + : static_cast(-1); OGS_FATAL( "RichardsMechanics: Biot-coefficient {} is smaller than " "porosity {} in element/integration point {}/{}.", - alpha, phi, this->element_.getID(), - -1 /*ip*/ /* TODO (CL) re-enable */); + alpha, phi, eid, ip); } auto const mu = liquid_phase.property(MPL::PropertyType::viscosity) @@ -915,9 +923,8 @@ void RichardsMechanicsLocalAssembler( ip_data, *medium, solid_phase, C_el, rho_LR, mu, - this->process_data_.micro_porosity_parameters, alpha, phi, p_cap_ip, - variables, variables_prev, x_position, t, dt, sigma_sw, - sigma_sw_prev); + micro_porosity_parameters, alpha, phi, p_cap_ip, variables, + variables_prev, x_position, t, dt, sigma_sw, sigma_sw_prev); } if (medium->hasProperty(MPL::PropertyType::transport_porosity)) @@ -1183,7 +1190,7 @@ void RichardsMechanicsLocalAssembler{ p_cap_ip, p_cap_prev_ip, Eigen::Vector::Zero()}, - CD, SD, SD_prev); + CD, SD, SD_prev, this->process_data_.micro_porosity_parameters); { auto const& C = *std::get>(CD); diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h index 16b82107462..d18799cba7d 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h @@ -342,7 +342,7 @@ class RichardsMechanicsLocalAssembler getMaterialStateVariablesAt(unsigned integration_point) const override; private: - void assembleWithJacobianEvalConstitutiveSetting( + static void assembleWithJacobianEvalConstitutiveSetting( double const t, double const dt, ParameterLib::SpatialPosition const& x_position, IpData& ip_data, MPL::VariableArray& variables, MPL::VariableArray& variables_prev, @@ -350,7 +350,9 @@ class RichardsMechanicsLocalAssembler CapillaryPressureData const& p_cap_data, ConstitutiveData& CD, StatefulData& SD, - StatefulDataPrev const& SD_prev); + StatefulDataPrev const& SD_prev, + std::optional const& + micro_porosity_parameters); std::vector> _ip_data; From 11d9ce6180ccbfcd4ea715a53a1d23ebc602d0bf Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 09:12:11 +0200 Subject: [PATCH 03/14] [PL/RM] Added missing includes --- ProcessLib/RichardsMechanics/IntegrationPointData.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h index 189fdcf4b47..78838eee745 100644 --- a/ProcessLib/RichardsMechanics/IntegrationPointData.h +++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h @@ -12,7 +12,10 @@ #include +#include "ConstitutiveRelations/Base.h" +#include "MaterialLib/SolidModels/MechanicsBase.h" #include "MathLib/KelvinVector.h" +#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/SolidMechanics.h" namespace ProcessLib { From 0b059fd0ac245b2fd7f287282c8c01139d8e3e81 Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 10:04:20 +0200 Subject: [PATCH 04/14] [PL/RM] Moved porosity, transport_porosity and saturation to stateful data --- .../ConstitutiveRelations/ConstitutiveData.h | 6 +- .../RichardsMechanics/IntegrationPointData.h | 10 +- .../RichardsMechanicsFEM-impl.h | 266 ++++++++++++++---- 3 files changed, 212 insertions(+), 70 deletions(-) diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h index d32cba9ac66..68fe7ac8e93 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h @@ -19,6 +19,7 @@ #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Porosity.h" #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Saturation.h" #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/SolidCompressibilityData.h" +#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TransportPorosity.h" #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/SolidMechanics.h" #include "SaturationSecantDerivative.h" #include "StiffnessTensor.h" @@ -34,7 +35,10 @@ using StatefulData = std::tuple< ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature:: SwellingDataStateful, ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature:: - MechanicalStrainData>; + MechanicalStrainData, + ProcessLib::ThermoRichardsMechanics::SaturationData, + ProcessLib::ThermoRichardsMechanics::PorosityData, + ProcessLib::ThermoRichardsMechanics::TransportPorosityData>; template using StatefulDataPrev = ProcessLib::ConstitutiveRelations::PrevStateOf< diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h index 78838eee745..9bf248410e7 100644 --- a/ProcessLib/RichardsMechanics/IntegrationPointData.h +++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h @@ -44,14 +44,9 @@ struct IntegrationPointData final double liquid_pressure_m = std::numeric_limits::quiet_NaN(); double liquid_pressure_m_prev = std::numeric_limits::quiet_NaN(); - double saturation = std::numeric_limits::quiet_NaN(); - double saturation_prev = std::numeric_limits::quiet_NaN(); double saturation_m = std::numeric_limits::quiet_NaN(); double saturation_m_prev = std::numeric_limits::quiet_NaN(); - double porosity = std::numeric_limits::quiet_NaN(); - double porosity_prev = std::numeric_limits::quiet_NaN(); - double transport_porosity = std::numeric_limits::quiet_NaN(); - double transport_porosity_prev = std::numeric_limits::quiet_NaN(); + double dry_density_solid = std::numeric_limits::quiet_NaN(); double dry_density_pellet_saturated = std::numeric_limits::quiet_NaN(); @@ -66,10 +61,7 @@ struct IntegrationPointData final void pushBackState() { - saturation_prev = saturation; saturation_m_prev = saturation_m; - porosity_prev = porosity; - transport_porosity_prev = transport_porosity; liquid_pressure_m_prev = liquid_pressure_m; material_state_variables->pushBackState(); } diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index ddf0947c7ab..46b40360f9b 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -46,7 +46,11 @@ void updateSwellingStressAndVolumetricStrain( SwellingDataStateful& sigma_sw, PrevState> const& sigma_sw_prev) + DisplacementDim>> const& sigma_sw_prev, + PrevState + phi_M_prev, + PrevState phi_prev, + ProcessLib::ThermoRichardsMechanics::TransportPorosityData& phi_M) { auto const& identity2 = MathLib::KelvinVector::Invariants< MathLib::KelvinVector::kelvin_vector_dimensions( @@ -81,9 +85,7 @@ void updateSwellingStressAndVolumetricStrain( // the micro_porosity_parameters. if (medium.hasProperty(MPL::PropertyType::saturation_micro)) { - double const phi_M_prev = ip_data.transport_porosity_prev; - double const phi_prev = ip_data.porosity_prev; - double const phi_m_prev = phi_prev - phi_M_prev; + double const phi_m_prev = phi_prev->phi - phi_M_prev->phi; double const p_L_m_prev = ip_data.liquid_pressure_m_prev; auto const S_L_m_prev = ip_data.saturation_m_prev; @@ -96,10 +98,9 @@ void updateSwellingStressAndVolumetricStrain( medium.property(MPL::PropertyType::saturation_micro), solid_phase.property(MPL::PropertyType::swelling_stress_rate)); - auto& phi_M = ip_data.transport_porosity; - phi_M = phi - (phi_m_prev + delta_phi_m); - variables_prev.transport_porosity = phi_M_prev; - variables.transport_porosity = phi_M; + phi_M.phi = phi - (phi_m_prev + delta_phi_m); + variables_prev.transport_porosity = phi_M_prev->phi; + variables.transport_porosity = phi_M.phi; auto& p_L_m = ip_data.liquid_pressure_m; p_L_m = p_L_m_prev + delta_p_L_m; @@ -174,17 +175,25 @@ RichardsMechanicsLocalAssemblerproperty(MPL::porosity) - .template initialValue( - x_position, - std::numeric_limits< - double>::quiet_NaN() /* t independent */); - - ip_data.transport_porosity = ip_data.porosity; + auto& porosity = + std::get( + this->current_states_[ip]) + .phi; + porosity = medium->property(MPL::porosity) + .template initialValue( + x_position, + std::numeric_limits< + double>::quiet_NaN() /* t independent */); + + auto& transport_porosity = + std::get< + ProcessLib::ThermoRichardsMechanics::TransportPorosityData>( + this->current_states_[ip]) + .phi; + transport_porosity = porosity; if (medium->hasProperty(MPL::PropertyType::transport_porosity)) { - ip_data.transport_porosity = + transport_porosity = medium->property(MPL::transport_porosity) .template initialValue( x_position, @@ -236,18 +245,38 @@ std::size_t RichardsMechanicsLocalAssembler< if (name == "saturation") { - return ProcessLib::setIntegrationPointScalarData(values, _ip_data, - &IpData::saturation); + return ProcessLib::setIntegrationPointScalarData( + values, this->current_states_, + [](auto& tuple) -> auto& + { + return std::get< + ProcessLib::ThermoRichardsMechanics::SaturationData>( + tuple) + .S_L; + }); } if (name == "porosity") { - return ProcessLib::setIntegrationPointScalarData(values, _ip_data, - &IpData::porosity); + return ProcessLib::setIntegrationPointScalarData( + values, this->current_states_, + [](auto& tuple) -> auto& + { + return std::get< + ProcessLib::ThermoRichardsMechanics::PorosityData>( + tuple) + .phi; + }); } if (name == "transport_porosity") { return ProcessLib::setIntegrationPointScalarData( - values, _ip_data, &IpData::transport_porosity); + values, this->current_states_, + [](auto& tuple) -> auto& + { + return std::get(tuple) + .phi; + }); } if (name == "swelling_stress") { @@ -341,9 +370,13 @@ void RichardsMechanicsLocalAssembler(variables, x_position, t, dt); variables.temperature = temperature; - _ip_data[ip].saturation_prev = - medium->property(MPL::PropertyType::saturation) - .template value(variables, x_position, t, dt); + auto& S_L_prev = + std::get< + PrevState>( + this->prev_states_[ip]) + ->S_L; + S_L_prev = medium->property(MPL::PropertyType::saturation) + .template value(variables, x_position, t, dt); if (medium->hasProperty(MPL::PropertyType::saturation_micro)) { @@ -475,8 +508,15 @@ void RichardsMechanicsLocalAssembler< std::get>(this->current_states_[ip]); eps.eps.noalias() = B * u; - auto& S_L = _ip_data[ip].saturation; - auto const S_L_prev = _ip_data[ip].saturation_prev; + auto& S_L = + std::get( + this->current_states_[ip]) + .S_L; + auto const S_L_prev = + std::get< + PrevState>( + this->prev_states_[ip]) + ->S_L; double p_cap_ip; NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip); @@ -549,9 +589,15 @@ void RichardsMechanicsLocalAssembler< variables.volumetric_strain = Invariants::trace(eps.eps); variables_prev.volumetric_strain = Invariants::trace(B * u_prev); - auto& phi = _ip_data[ip].porosity; + auto& phi = std::get( + this->current_states_[ip]) + .phi; { // Porosity update - variables_prev.porosity = _ip_data[ip].porosity_prev; + auto const phi_prev = std::get>( + this->prev_states_[ip]) + ->phi; + variables_prev.porosity = phi_prev; phi = medium->property(MPL::PropertyType::porosity) .template value(variables, variables_prev, x_position, t, dt); @@ -604,14 +650,23 @@ void RichardsMechanicsLocalAssembler< if (medium->hasProperty(MPL::PropertyType::transport_porosity)) { - variables_prev.transport_porosity = - _ip_data[ip].transport_porosity_prev; - - _ip_data[ip].transport_porosity = + auto& transport_porosity = + std::get( + this->current_states_[ip]) + .phi; + auto const transport_porosity_prev = + std::get>( + this->prev_states_[ip]) + ->phi; + variables_prev.transport_porosity = transport_porosity_prev; + + transport_porosity = medium->property(MPL::PropertyType::transport_porosity) .template value(variables, variables_prev, x_position, t, dt); - variables.transport_porosity = _ip_data[ip].transport_porosity; + variables.transport_porosity = transport_porosity; } else { @@ -801,8 +856,13 @@ void RichardsMechanicsLocalAssembler>(SD); - auto& S_L = ip_data.saturation; - auto const S_L_prev = ip_data.saturation_prev; + auto& S_L = + std::get(SD).S_L; + auto const S_L_prev = + std::get< + PrevState>( + SD_prev) + ->S_L; auto const alpha = medium->property(MPL::PropertyType::biot_coefficient) .template value(variables, x_position, t, dt); @@ -877,10 +937,15 @@ void RichardsMechanicsLocalAssembler>>(SD_prev)->eps); - auto& phi = ip_data.porosity; + auto& phi = + std::get(SD).phi; { // Porosity update - - variables_prev.porosity = ip_data.porosity_prev; + auto const phi_prev = + std::get< + PrevState>( + SD_prev) + ->phi; + variables_prev.porosity = phi_prev; phi = medium->property(MPL::PropertyType::porosity) .template value(variables, variables_prev, x_position, t, dt); @@ -920,24 +985,42 @@ void RichardsMechanicsLocalAssembler>>( SD_prev); + auto const transport_porosity_prev = std::get>( + SD_prev); + auto const phi_prev = std::get< + PrevState>( + SD_prev); + auto& transport_porosity = std::get< + ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(SD); updateSwellingStressAndVolumetricStrain( ip_data, *medium, solid_phase, C_el, rho_LR, mu, micro_porosity_parameters, alpha, phi, p_cap_ip, variables, - variables_prev, x_position, t, dt, sigma_sw, sigma_sw_prev); + variables_prev, x_position, t, dt, sigma_sw, sigma_sw_prev, + transport_porosity_prev, phi_prev, transport_porosity); } if (medium->hasProperty(MPL::PropertyType::transport_porosity)) { if (!medium->hasProperty(MPL::PropertyType::saturation_micro)) { - variables_prev.transport_porosity = ip_data.transport_porosity_prev; - - ip_data.transport_porosity = + auto& transport_porosity = + std::get< + ProcessLib::ThermoRichardsMechanics::TransportPorosityData>( + SD) + .phi; + auto const transport_porosity_prev = std::get>( + SD_prev) + ->phi; + variables_prev.transport_porosity = transport_porosity_prev; + + transport_porosity = medium->property(MPL::PropertyType::transport_porosity) .template value(variables, variables_prev, x_position, t, dt); - variables.transport_porosity = ip_data.transport_porosity; + variables.transport_porosity = transport_porosity; } } else @@ -1284,7 +1367,10 @@ void RichardsMechanicsLocalAssembler( + this->current_states_[ip]) + .S_L; if (this->process_data_.explicit_hm_coupling_in_unsaturated_zone) { double const chi_S_L_prev = std::get>( + this->prev_states_[ip]) + ->S_L; storage_p_a_S_Jpp.noalias() -= N_p.transpose() * rho_LR * ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap + @@ -1642,7 +1732,15 @@ std::vector const& RichardsMechanicsLocalAssembler< std::vector& cache) const { return ProcessLib::getIntegrationPointScalarData( - _ip_data, &IpData::saturation, cache); + this->current_states_, + [](auto& tuple) -> auto& + { + return std::get< + ProcessLib::ThermoRichardsMechanics::SaturationData>( + tuple) + .S_L; + }, + cache); } template const& RichardsMechanicsLocalAssembler< std::vector const& /*dof_table*/, std::vector& cache) const { - return ProcessLib::getIntegrationPointScalarData(_ip_data, - &IpData::porosity, cache); + return ProcessLib::getIntegrationPointScalarData( + this->current_states_, + [](auto& tuple) -> auto& + { + return std::get( + tuple) + .phi; + }, + cache); } template const& RichardsMechanicsLocalAssembler< std::vector& cache) const { return ProcessLib::getIntegrationPointScalarData( - _ip_data, &IpData::transport_porosity, cache); + this->current_states_, + [](auto& tuple) -> auto& + { + return std:: + get( + tuple) + .phi; + }, + cache); } template >(this->current_states_[ip]) .eps; eps.noalias() = B * u; - auto& S_L = _ip_data[ip].saturation; - auto const S_L_prev = _ip_data[ip].saturation_prev; + auto& S_L = + std::get( + this->current_states_[ip]) + .S_L; + auto const S_L_prev = + std::get< + PrevState>( + this->prev_states_[ip]) + ->S_L; S_L = medium->property(MPL::PropertyType::saturation) .template value(variables, x_position, t, dt); variables.liquid_saturation = S_L; @@ -1929,9 +2049,15 @@ void RichardsMechanicsLocalAssembler( + this->current_states_[ip]) + .phi; { // Porosity update - variables_prev.porosity = _ip_data[ip].porosity_prev; + auto const phi_prev = std::get>( + this->prev_states_[ip]) + ->phi; + variables_prev.porosity = phi_prev; phi = medium->property(MPL::PropertyType::porosity) .template value(variables, variables_prev, x_position, t, dt); @@ -1958,26 +2084,46 @@ void RichardsMechanicsLocalAssembler>>( this->prev_states_[ip]); + auto const transport_porosity_prev = std::get>( + this->prev_states_[ip]); + auto const phi_prev = std::get< + PrevState>( + this->prev_states_[ip]); + auto& transport_porosity = std::get< + ProcessLib::ThermoRichardsMechanics::TransportPorosityData>( + this->current_states_[ip]); updateSwellingStressAndVolumetricStrain( _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu, this->process_data_.micro_porosity_parameters, alpha, phi, p_cap_ip, variables, variables_prev, x_position, t, dt, - sigma_sw, sigma_sw_prev); + sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev, + transport_porosity); } if (medium->hasProperty(MPL::PropertyType::transport_porosity)) { if (!medium->hasProperty(MPL::PropertyType::saturation_micro)) { - variables_prev.transport_porosity = - _ip_data[ip].transport_porosity_prev; - - _ip_data[ip].transport_porosity = + auto& transport_porosity = + std::get( + this->current_states_[ip]) + .phi; + auto const transport_porosity_prev = + std::get>( + this->prev_states_[ip]) + ->phi; + + variables_prev.transport_porosity = transport_porosity_prev; + + transport_porosity = medium->property(MPL::PropertyType::transport_porosity) .template value(variables, variables_prev, x_position, t, dt); - variables.transport_porosity = _ip_data[ip].transport_porosity; + variables.transport_porosity = transport_porosity; } } else From 958adffd55728f097a8c5df9840b4fe1d946ae36 Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 10:42:17 +0200 Subject: [PATCH 05/14] [PL/RM] Moved micro saturation and pressure to stateful data --- .../ConstitutiveRelations/ConstitutiveData.h | 5 +- .../ConstitutiveRelations/MicroPressure.h | 18 +++++ .../ConstitutiveRelations/MicroSaturation.h | 18 +++++ .../RichardsMechanics/IntegrationPointData.h | 12 +-- .../RichardsMechanicsFEM-impl.h | 79 ++++++++++++------- 5 files changed, 93 insertions(+), 39 deletions(-) create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h index 68fe7ac8e93..242707fb1ef 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h @@ -12,6 +12,8 @@ #include "Density.h" #include "LiquidDensity.h" +#include "MicroPressure.h" +#include "MicroSaturation.h" #include "ProcessLib/ConstitutiveRelations/Base.h" #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Biot.h" #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/LiquidViscosity.h" @@ -38,7 +40,8 @@ using StatefulData = std::tuple< MechanicalStrainData, ProcessLib::ThermoRichardsMechanics::SaturationData, ProcessLib::ThermoRichardsMechanics::PorosityData, - ProcessLib::ThermoRichardsMechanics::TransportPorosityData>; + ProcessLib::ThermoRichardsMechanics::TransportPorosityData, MicroPressure, + MicroSaturation>; template using StatefulDataPrev = ProcessLib::ConstitutiveRelations::PrevStateOf< diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h new file mode 100644 index 00000000000..541be0780ca --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h @@ -0,0 +1,18 @@ +/** + * \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 "BaseLib/StrongType.h" + +namespace ProcessLib::RichardsMechanics +{ +using MicroPressure = BaseLib::StrongType; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h new file mode 100644 index 00000000000..2456caa8ed3 --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h @@ -0,0 +1,18 @@ +/** + * \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 "BaseLib/StrongType.h" + +namespace ProcessLib::RichardsMechanics +{ +using MicroSaturation = BaseLib::StrongType; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h index 9bf248410e7..5c3908bdd70 100644 --- a/ProcessLib/RichardsMechanics/IntegrationPointData.h +++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h @@ -42,11 +42,6 @@ struct IntegrationPointData final typename ShapeMatricesTypePressure::GlobalDimVectorType v_darcy; - double liquid_pressure_m = std::numeric_limits::quiet_NaN(); - double liquid_pressure_m_prev = std::numeric_limits::quiet_NaN(); - double saturation_m = std::numeric_limits::quiet_NaN(); - double saturation_m_prev = std::numeric_limits::quiet_NaN(); - double dry_density_solid = std::numeric_limits::quiet_NaN(); double dry_density_pellet_saturated = std::numeric_limits::quiet_NaN(); @@ -59,12 +54,7 @@ struct IntegrationPointData final material_state_variables; double integration_weight = std::numeric_limits::quiet_NaN(); - void pushBackState() - { - saturation_m_prev = saturation_m; - liquid_pressure_m_prev = liquid_pressure_m; - material_state_variables->pushBackState(); - } + void pushBackState() { material_state_variables->pushBackState(); } MathLib::KelvinVector::KelvinMatrixType computeElasticTangentStiffness( diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index 46b40360f9b..d1ff1775d73 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -47,10 +47,13 @@ void updateSwellingStressAndVolumetricStrain( PrevState> const& sigma_sw_prev, - PrevState + PrevState const phi_M_prev, - PrevState phi_prev, - ProcessLib::ThermoRichardsMechanics::TransportPorosityData& phi_M) + PrevState const phi_prev, + ProcessLib::ThermoRichardsMechanics::TransportPorosityData& phi_M, + PrevState const p_L_m_prev, + PrevState const S_L_m_prev, MicroPressure& p_L_m, + MicroSaturation& S_L_m) { auto const& identity2 = MathLib::KelvinVector::Invariants< MathLib::KelvinVector::kelvin_vector_dimensions( @@ -86,15 +89,12 @@ void updateSwellingStressAndVolumetricStrain( if (medium.hasProperty(MPL::PropertyType::saturation_micro)) { double const phi_m_prev = phi_prev->phi - phi_M_prev->phi; - double const p_L_m_prev = ip_data.liquid_pressure_m_prev; - - auto const S_L_m_prev = ip_data.saturation_m_prev; auto const [delta_phi_m, delta_e_sw, delta_p_L_m, delta_sigma_sw] = computeMicroPorosity( identity2.transpose() * C_el.inverse(), rho_LR, mu, - *micro_porosity_parameters, alpha, phi, -p_cap_ip, p_L_m_prev, - variables_prev, S_L_m_prev, phi_m_prev, x_position, t, dt, + *micro_porosity_parameters, alpha, phi, -p_cap_ip, **p_L_m_prev, + variables_prev, **S_L_m_prev, phi_m_prev, x_position, t, dt, medium.property(MPL::PropertyType::saturation_micro), solid_phase.property(MPL::PropertyType::swelling_stress_rate)); @@ -102,17 +102,15 @@ void updateSwellingStressAndVolumetricStrain( variables_prev.transport_porosity = phi_M_prev->phi; variables.transport_porosity = phi_M.phi; - auto& p_L_m = ip_data.liquid_pressure_m; - p_L_m = p_L_m_prev + delta_p_L_m; + *p_L_m = **p_L_m_prev + delta_p_L_m; { // Update micro saturation. MPL::VariableArray variables_prev; - variables_prev.capillary_pressure = -p_L_m_prev; + variables_prev.capillary_pressure = -**p_L_m_prev; MPL::VariableArray variables; - variables.capillary_pressure = -p_L_m; + variables.capillary_pressure = -*p_L_m; - ip_data.saturation_m = - medium.property(MPL::PropertyType::saturation_micro) - .template value(variables, x_position, t, dt); + *S_L_m = medium.property(MPL::PropertyType::saturation_micro) + .template value(variables, x_position, t, dt); } sigma_sw.sigma_sw = sigma_sw_prev->sigma_sw + delta_sigma_sw; } @@ -362,8 +360,14 @@ void RichardsMechanicsLocalAssembler(this->current_states_[ip]); + auto& p_L_m_prev = + std::get>(this->prev_states_[ip]); + **p_L_m_prev = -p_cap_ip; + *p_L_m = -p_cap_ip; + } auto const temperature = medium->property(MPL::PropertyType::reference_temperature) @@ -382,10 +386,14 @@ void RichardsMechanicsLocalAssemblerproperty(MPL::PropertyType::saturation_micro) - .template value(vars, x_position, t, dt); - _ip_data[ip].saturation_m_prev = S_L_m; + + auto& S_L_m = std::get(this->current_states_[ip]); + auto& S_L_m_prev = + std::get>(this->prev_states_[ip]); + + *S_L_m = medium->property(MPL::PropertyType::saturation_micro) + .template value(vars, x_position, t, dt); + *S_L_m_prev = S_L_m; } // Set eps_m_prev from potentially non-zero eps and sigma_sw from @@ -993,12 +1001,17 @@ void RichardsMechanicsLocalAssembler(SD); + auto& p_L_m = std::get(SD); + auto const p_L_m_prev = std::get>(SD_prev); + auto& S_L_m = std::get(SD); + auto const S_L_m_prev = std::get>(SD_prev); updateSwellingStressAndVolumetricStrain( ip_data, *medium, solid_phase, C_el, rho_LR, mu, micro_porosity_parameters, alpha, phi, p_cap_ip, variables, variables_prev, x_position, t, dt, sigma_sw, sigma_sw_prev, - transport_porosity_prev, phi_prev, transport_porosity); + transport_porosity_prev, phi_prev, transport_porosity, p_L_m_prev, + S_L_m_prev, p_L_m, S_L_m); } if (medium->hasProperty(MPL::PropertyType::transport_porosity)) @@ -1490,7 +1503,8 @@ void RichardsMechanicsLocalAssemblerprocess_data_.micro_porosity_parameters ->mass_exchange_coefficient; - auto const p_L_m = _ip_data[ip].liquid_pressure_m; + auto const p_L_m = + *std::get(this->current_states_[ip]); local_rhs.template segment(pressure_index) .noalias() -= N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w; @@ -1501,7 +1515,8 @@ void RichardsMechanicsLocalAssembler>( + this->prev_states_[ip]); local_Jac .template block( pressure_index, pressure_index) @@ -1765,7 +1780,9 @@ std::vector const& RichardsMechanicsLocalAssembler< std::vector& cache) const { return ProcessLib::getIntegrationPointScalarData( - _ip_data, &IpData::saturation_m, cache); + this->current_states_, + [](auto& tuple) -> auto& { return *std::get(tuple); }, + cache); } template const& RichardsMechanicsLocalAssembler< std::vector& cache) const { return ProcessLib::getIntegrationPointScalarData( - _ip_data, &IpData::liquid_pressure_m, cache); + this->current_states_, + [](auto& tuple) -> auto& { return *std::get(tuple); }, + cache); } template ( this->current_states_[ip]); + auto& p_L_m = std::get(this->current_states_[ip]); + auto const p_L_m_prev = + std::get>(this->prev_states_[ip]); + auto& S_L_m = std::get(this->current_states_[ip]); + auto const S_L_m_prev = + std::get>(this->prev_states_[ip]); updateSwellingStressAndVolumetricStrain( _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu, this->process_data_.micro_porosity_parameters, alpha, phi, p_cap_ip, variables, variables_prev, x_position, t, dt, sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev, - transport_porosity); + transport_porosity, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m); } if (medium->hasProperty(MPL::PropertyType::transport_porosity)) From 9e938fe2de40f04f0f4215be2d3ec923fd2c5ff3 Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 10:47:39 +0200 Subject: [PATCH 06/14] [PL/RM] Removed obsolete function argument --- .../RichardsMechanics/RichardsMechanicsFEM-impl.h | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index d1ff1775d73..c76007b6d9b 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -31,9 +31,9 @@ namespace ProcessLib { namespace RichardsMechanics { -template +template void updateSwellingStressAndVolumetricStrain( - IPData& ip_data, MaterialPropertyLib::Medium const& medium, + MaterialPropertyLib::Medium const& medium, MaterialPropertyLib::Phase const& solid_phase, MathLib::KelvinVector::KelvinMatrixType const& C_el, double const rho_LR, double const mu, @@ -1007,11 +1007,10 @@ void RichardsMechanicsLocalAssembler>(SD_prev); updateSwellingStressAndVolumetricStrain( - ip_data, *medium, solid_phase, C_el, rho_LR, mu, - micro_porosity_parameters, alpha, phi, p_cap_ip, variables, - variables_prev, x_position, t, dt, sigma_sw, sigma_sw_prev, - transport_porosity_prev, phi_prev, transport_porosity, p_L_m_prev, - S_L_m_prev, p_L_m, S_L_m); + *medium, solid_phase, C_el, rho_LR, mu, micro_porosity_parameters, + alpha, phi, p_cap_ip, variables, variables_prev, x_position, t, dt, + sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev, + transport_porosity, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m); } if (medium->hasProperty(MPL::PropertyType::transport_porosity)) @@ -2120,7 +2119,7 @@ void RichardsMechanicsLocalAssembler>(this->prev_states_[ip]); updateSwellingStressAndVolumetricStrain( - _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu, + *medium, solid_phase, C_el, rho_LR, mu, this->process_data_.micro_porosity_parameters, alpha, phi, p_cap_ip, variables, variables_prev, x_position, t, dt, sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev, From 75b938cb237a6bd89c2fbc203cd466e01eea126d Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 11:03:27 +0200 Subject: [PATCH 07/14] [PL/RM] Moved Darcy velocity to output data --- .../ConstitutiveRelations/ConstitutiveData.h | 4 +++- ProcessLib/RichardsMechanics/IntegrationPointData.h | 2 -- .../RichardsMechanics/RichardsMechanicsFEM-impl.h | 10 +++++++--- 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h index 242707fb1ef..417faf5bc92 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h @@ -16,6 +16,7 @@ #include "MicroSaturation.h" #include "ProcessLib/ConstitutiveRelations/Base.h" #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Biot.h" +#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/DarcyLaw.h" #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/LiquidViscosity.h" #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityData.h" #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Porosity.h" @@ -49,7 +50,8 @@ using StatefulDataPrev = ProcessLib::ConstitutiveRelations::PrevStateOf< /// Data that is needed for output purposes, but not directly for the assembly. template -using OutputData = std::tuple<>; +using OutputData = std::tuple< + ProcessLib::ThermoRichardsMechanics::DarcyLawData>; /// Data that is needed for the equation system assembly. template diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h index 5c3908bdd70..ee813f5b7ba 100644 --- a/ProcessLib/RichardsMechanics/IntegrationPointData.h +++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h @@ -40,8 +40,6 @@ struct IntegrationPointData final typename ShapeMatricesTypePressure::NodalRowVectorType N_p; typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p; - typename ShapeMatricesTypePressure::GlobalDimVectorType v_darcy; - double dry_density_solid = std::numeric_limits::quiet_NaN(); double dry_density_pellet_saturated = std::numeric_limits::quiet_NaN(); diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index c76007b6d9b..99c9e940d94 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -1718,7 +1718,9 @@ std::vector const& RichardsMechanicsLocalAssembler< for (unsigned ip = 0; ip < n_integration_points; ip++) { - cache_matrix.col(ip).noalias() = _ip_data[ip].v_darcy; + cache_matrix.col(ip).noalias() = *std::get< + ProcessLib::ThermoRichardsMechanics::DarcyLawData>( + this->output_data_[ip]); } return cache; @@ -2247,8 +2249,10 @@ void RichardsMechanicsLocalAssembler>( + this->output_data_[ip]) + ->noalias() = -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b; saturation_avg += S_L; porosity_avg += phi; From 9cae9efbc734e0c1e3ae2f836c78ec397315a69f Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 11:09:26 +0200 Subject: [PATCH 08/14] [PL/RM] Removed unused data members --- ProcessLib/RichardsMechanics/IntegrationPointData.h | 4 ---- 1 file changed, 4 deletions(-) diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h index ee813f5b7ba..23076399f5f 100644 --- a/ProcessLib/RichardsMechanics/IntegrationPointData.h +++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h @@ -41,10 +41,6 @@ struct IntegrationPointData final typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p; double dry_density_solid = std::numeric_limits::quiet_NaN(); - double dry_density_pellet_saturated = - std::numeric_limits::quiet_NaN(); - double dry_density_pellet_unsaturated = - std::numeric_limits::quiet_NaN(); MaterialLib::Solids::MechanicsBase const& solid_material; std::unique_ptr Date: Tue, 11 Jun 2024 11:20:27 +0200 Subject: [PATCH 09/14] [PL/RM] Moved dry solid density to output data --- .../ConstitutiveRelations/ConstitutiveData.h | 4 +++- .../ConstitutiveRelations/DrySolidDensity.h | 19 +++++++++++++++++++ .../RichardsMechanics/IntegrationPointData.h | 2 -- .../RichardsMechanicsFEM-impl.h | 6 ++++-- 4 files changed, 26 insertions(+), 5 deletions(-) create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h index 417faf5bc92..2641de6abc8 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h @@ -11,6 +11,7 @@ #pragma once #include "Density.h" +#include "DrySolidDensity.h" #include "LiquidDensity.h" #include "MicroPressure.h" #include "MicroSaturation.h" @@ -51,7 +52,8 @@ using StatefulDataPrev = ProcessLib::ConstitutiveRelations::PrevStateOf< /// Data that is needed for output purposes, but not directly for the assembly. template using OutputData = std::tuple< - ProcessLib::ThermoRichardsMechanics::DarcyLawData>; + ProcessLib::ThermoRichardsMechanics::DarcyLawData, + DrySolidDensity>; /// Data that is needed for the equation system assembly. template diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h new file mode 100644 index 00000000000..a2af87eae9f --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h @@ -0,0 +1,19 @@ +/** + * \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 "BaseLib/StrongType.h" + +namespace ProcessLib::RichardsMechanics +{ +// Apparent dry solid density +using DrySolidDensity = BaseLib::StrongType; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h index 23076399f5f..a1276108f9f 100644 --- a/ProcessLib/RichardsMechanics/IntegrationPointData.h +++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h @@ -40,8 +40,6 @@ struct IntegrationPointData final typename ShapeMatricesTypePressure::NodalRowVectorType N_p; typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p; - double dry_density_solid = std::numeric_limits::quiet_NaN(); - MaterialLib::Solids::MechanicsBase const& solid_material; std::unique_ptr::MaterialStateVariables> diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index 99c9e940d94..4192ce317b5 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -1889,7 +1889,9 @@ std::vector const& RichardsMechanicsLocalAssembler< std::vector& cache) const { return ProcessLib::getIntegrationPointScalarData( - _ip_data, &IpData::dry_density_solid, cache); + this->output_data_, + [](auto& tuple) -> auto& { return *std::get(tuple); }, + cache); } template (variables, x_position, t, dt); - _ip_data[ip].dry_density_solid = (1 - phi) * rho_SR; + *std::get(this->output_data_[ip]) = (1 - phi) * rho_SR; { auto& SD = this->current_states_[ip]; From 40a8bb45dd27c2ecf89c8d00713b7c788e4c195e Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 11:42:02 +0200 Subject: [PATCH 10/14] [PL/RM] Use reflection for output --- .../ConstitutiveRelations/DrySolidDensity.h | 7 + .../ConstitutiveRelations/MicroPressure.h | 7 + .../ConstitutiveRelations/MicroSaturation.h | 7 + .../LocalAssemblerInterface.h | 84 +---- .../RichardsMechanicsFEM-impl.h | 311 ------------------ .../RichardsMechanics/RichardsMechanicsFEM.h | 69 ---- .../RichardsMechanicsProcess.cpp | 82 +---- 7 files changed, 39 insertions(+), 528 deletions(-) diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h index a2af87eae9f..c021158782c 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h @@ -10,10 +10,17 @@ #pragma once +#include + #include "BaseLib/StrongType.h" namespace ProcessLib::RichardsMechanics { // Apparent dry solid density using DrySolidDensity = BaseLib::StrongType; + +constexpr std::string_view ioName(struct DrySolidDensityTag*) +{ + return "dry_density_solid"; +} } // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h index 541be0780ca..96ace800113 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h @@ -10,9 +10,16 @@ #pragma once +#include + #include "BaseLib/StrongType.h" namespace ProcessLib::RichardsMechanics { using MicroPressure = BaseLib::StrongType; + +constexpr std::string_view ioName(struct MicroPressureTag*) +{ + return "micro_pressure"; +} } // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h index 2456caa8ed3..4048fdd58fe 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h @@ -10,9 +10,16 @@ #pragma once +#include + #include "BaseLib/StrongType.h" namespace ProcessLib::RichardsMechanics { using MicroSaturation = BaseLib::StrongType; + +constexpr std::string_view ioName(struct MicroSaturationTag*) +{ + return "micro_saturation"; +} } // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h index 529002d1f65..a9a57f99fdf 100644 --- a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h +++ b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h @@ -66,82 +66,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, std::string_view const name, double const* values, int const integration_order) = 0; - virtual std::vector getSigma() const = 0; - - virtual std::vector const& getIntPtSigma( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector getSwellingStress() const = 0; - - virtual std::vector const& getIntPtSwellingStress( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector getEpsilon() const = 0; - - virtual std::vector const& getIntPtEpsilon( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector const& getIntPtDarcyVelocity( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector getSaturation() const = 0; - - virtual std::vector const& getIntPtSaturation( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector getMicroSaturation() const = 0; - - virtual std::vector const& getIntPtMicroSaturation( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector getMicroPressure() const = 0; - - virtual std::vector const& getIntPtMicroPressure( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector getPorosity() const = 0; - - virtual std::vector const& getIntPtPorosity( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector getTransportPorosity() const = 0; - - virtual std::vector const& getIntPtTransportPorosity( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector const& getIntPtDryDensitySolid( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - virtual std::vector getMaterialStateVariableInternalState( std::function( typename MaterialLib::Solids::MechanicsBase:: @@ -157,6 +81,14 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, DisplacementDim>::MaterialStateVariables const& getMaterialStateVariablesAt(unsigned /*integration_point*/) const = 0; + static auto getReflectionDataForOutput() + { + using Self = LocalAssemblerInterface; + + return ProcessLib::Reflection::reflectWithoutName( + &Self::current_states_, &Self::output_data_); + } + protected: RichardsMechanicsProcessData& process_data_; NumLib::GenericIntegrationMethod const& integration_method_; diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index 4192ce317b5..f1bd168a25d 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -1557,121 +1557,6 @@ void RichardsMechanicsLocalAssembler -std::vector RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, - DisplacementDim>::getSigma() const -{ - constexpr int kelvin_vector_size = - MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim); - - return transposeInPlace( - [this](std::vector& values) - { return getIntPtSigma(0, {}, {}, values); }); -} - -template -std::vector const& RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtSigma( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - return ProcessLib::getIntegrationPointKelvinVectorData( - this->current_states_, - [](auto& tuple) -> auto& { - return std::get>(tuple) - .sigma_eff; - }, - cache); -} - -template -std::vector RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, - DisplacementDim>::getSwellingStress() const -{ - constexpr int kelvin_vector_size = - MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim); - - return transposeInPlace( - [this](std::vector& values) - { return getIntPtSwellingStress(0, {}, {}, values); }); -} - -template -std::vector const& RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtSwellingStress( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - constexpr int kelvin_vector_size = - MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim); - auto const n_integration_points = _ip_data.size(); - - cache.clear(); - auto cache_mat = MathLib::createZeroedMatrix>( - cache, kelvin_vector_size, n_integration_points); - - for (unsigned ip = 0; ip < n_integration_points; ++ip) - { - auto const& sigma_sw = - std::get>( - this->current_states_[ip]) - .sigma_sw; - cache_mat.col(ip) = - MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_sw); - } - - return cache; -} - -template -std::vector RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, - DisplacementDim>::getEpsilon() const -{ - constexpr int kelvin_vector_size = - MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim); - - return transposeInPlace( - [this](std::vector& values) - { return getIntPtEpsilon(0, {}, {}, values); }); -} - -template -std::vector const& RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtEpsilon( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - return ProcessLib::getIntegrationPointKelvinVectorData( - this->current_states_, - [](auto& tuple) -> auto& { - return std::get>(tuple).eps; - }, - cache); -} - template int RichardsMechanicsLocalAssembler RichardsMechanicsLocalAssembler< n_components); } -template -std::vector const& RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtDarcyVelocity( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - unsigned const n_integration_points = - this->integration_method_.getNumberOfPoints(); - - cache.clear(); - auto cache_matrix = MathLib::createZeroedMatrix>( - cache, DisplacementDim, n_integration_points); - - for (unsigned ip = 0; ip < n_integration_points; ip++) - { - cache_matrix.col(ip).noalias() = *std::get< - ProcessLib::ThermoRichardsMechanics::DarcyLawData>( - this->output_data_[ip]); - } - - return cache; -} - -template -std::vector RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, - DisplacementDim>::getSaturation() const -{ - std::vector result; - getIntPtSaturation(0, {}, {}, result); - return result; -} - -template -std::vector const& RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtSaturation( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - return ProcessLib::getIntegrationPointScalarData( - this->current_states_, - [](auto& tuple) -> auto& - { - return std::get< - ProcessLib::ThermoRichardsMechanics::SaturationData>( - tuple) - .S_L; - }, - cache); -} - -template -std::vector RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, - DisplacementDim>::getMicroSaturation() const -{ - std::vector result; - getIntPtMicroSaturation(0, {}, {}, result); - return result; -} - -template -std::vector const& RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtMicroSaturation( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - return ProcessLib::getIntegrationPointScalarData( - this->current_states_, - [](auto& tuple) -> auto& { return *std::get(tuple); }, - cache); -} - -template -std::vector RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, - DisplacementDim>::getMicroPressure() const -{ - std::vector result; - getIntPtMicroPressure(0, {}, {}, result); - return result; -} - -template -std::vector const& RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtMicroPressure( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - return ProcessLib::getIntegrationPointScalarData( - this->current_states_, - [](auto& tuple) -> auto& { return *std::get(tuple); }, - cache); -} - -template -std::vector RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, - DisplacementDim>::getPorosity() const -{ - std::vector result; - getIntPtPorosity(0, {}, {}, result); - return result; -} - -template -std::vector const& RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtPorosity( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - return ProcessLib::getIntegrationPointScalarData( - this->current_states_, - [](auto& tuple) -> auto& - { - return std::get( - tuple) - .phi; - }, - cache); -} - -template -std::vector RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, - DisplacementDim>::getTransportPorosity() const -{ - std::vector result; - getIntPtTransportPorosity(0, {}, {}, result); - return result; -} - -template -std::vector const& RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtTransportPorosity( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - return ProcessLib::getIntegrationPointScalarData( - this->current_states_, - [](auto& tuple) -> auto& - { - return std:: - get( - tuple) - .phi; - }, - cache); -} - -template -std::vector const& RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtDryDensitySolid( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - return ProcessLib::getIntegrationPointScalarData( - this->output_data_, - [](auto& tuple) -> auto& { return *std::get(tuple); }, - cache); -} - template void RichardsMechanicsLocalAssembler(N_u.data(), N_u.size()); } - std::vector getSigma() const override; - - std::vector const& getIntPtDarcyVelocity( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - - std::vector getSaturation() const override; - std::vector const& getIntPtSaturation( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - - std::vector getMicroSaturation() const override; - std::vector const& getIntPtMicroSaturation( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - - std::vector getMicroPressure() const override; - std::vector const& getIntPtMicroPressure( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - - std::vector getPorosity() const override; - std::vector const& getIntPtPorosity( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - - std::vector getTransportPorosity() const override; - std::vector const& getIntPtTransportPorosity( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - - std::vector const& getIntPtSigma( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - - std::vector getSwellingStress() const override; - std::vector const& getIntPtSwellingStress( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - - std::vector getEpsilon() const override; - std::vector const& getIntPtEpsilon( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - int getMaterialID() const override; std::vector getMaterialStateVariableInternalState( @@ -267,12 +204,6 @@ class RichardsMechanicsLocalAssembler MaterialStateVariables&)> const& get_values_span, int const& n_components) const override; - std::vector const& getIntPtDryDensitySolid( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - private: /** * Assemble local matrices and vectors arise from the linearized discretized diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp index ff36291ec52..4ab85ce74c2 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp @@ -16,6 +16,8 @@ #include "MeshLib/Utils/getOrCreateMeshProperty.h" #include "NumLib/DOF/ComputeSparsityPattern.h" #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h" +#include "ProcessLib/Reflection/ReflectionForExtrapolation.h" +#include "ProcessLib/Reflection/ReflectionForIPWriters.h" #include "ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h" #include "ProcessLib/Utils/SetIPDataInitialConditions.h" #include "RichardsMechanicsFEM.h" @@ -48,43 +50,10 @@ RichardsMechanicsProcess::RichardsMechanicsProcess( _hydraulic_flow = MeshLib::getOrCreateMeshProperty( mesh, "MassFlowRate", MeshLib::MeshItemType::Node, 1); - // TODO (naumov) remove ip suffix. Probably needs modification of the mesh - // properties, s.t. there is no "overlapping" with cell/point data. - // See getOrCreateMeshProperty. - _integration_point_writer.emplace_back( - std::make_unique( - "sigma_ip", - static_cast(mesh.getDimension() == 2 ? 4 : 6) /*n components*/, - integration_order, _local_assemblers, &LocalAssemblerIF::getSigma)); - - _integration_point_writer.emplace_back( - std::make_unique( - "saturation_ip", 1 /*n components*/, integration_order, - _local_assemblers, &LocalAssemblerIF::getSaturation)); - - _integration_point_writer.emplace_back( - std::make_unique( - "porosity_ip", 1 /*n components*/, integration_order, - _local_assemblers, &LocalAssemblerIF::getPorosity)); - - _integration_point_writer.emplace_back( - std::make_unique( - "transport_porosity_ip", 1 /*n components*/, integration_order, - _local_assemblers, &LocalAssemblerIF::getTransportPorosity)); - - _integration_point_writer.emplace_back( - std::make_unique( - "swelling_stress_ip", - static_cast(mesh.getDimension() == 2 ? 4 : 6) /*n components*/, - integration_order, _local_assemblers, - &LocalAssemblerIF::getSwellingStress)); - - _integration_point_writer.emplace_back( - std::make_unique( - "epsilon_ip", - static_cast(mesh.getDimension() == 2 ? 4 : 6) /*n components*/, - integration_order, _local_assemblers, - &LocalAssemblerIF::getEpsilon)); + ProcessLib::Reflection::addReflectedIntegrationPointWriters< + DisplacementDim>(LocalAssemblerIF::getReflectionDataForOutput(), + _integration_point_writer, integration_order, + _local_assemblers); } template @@ -202,6 +171,10 @@ void RichardsMechanicsProcess::initializeConcreteProcess( NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(), _process_data); + ProcessLib::Reflection::addReflectedSecondaryVariables( + LocalAssemblerIF::getReflectionDataForOutput(), _secondary_variables, + getExtrapolator(), _local_assemblers); + auto add_secondary_variable = [&](std::string const& name, int const num_components, auto get_ip_values_function) @@ -213,41 +186,6 @@ void RichardsMechanicsProcess::initializeConcreteProcess( std::move(get_ip_values_function))); }; - add_secondary_variable("sigma", - MathLib::KelvinVector::KelvinVectorType< - DisplacementDim>::RowsAtCompileTime, - &LocalAssemblerIF::getIntPtSigma); - - add_secondary_variable("swelling_stress", - MathLib::KelvinVector::KelvinVectorType< - DisplacementDim>::RowsAtCompileTime, - &LocalAssemblerIF::getIntPtSwellingStress); - - add_secondary_variable("epsilon", - MathLib::KelvinVector::KelvinVectorType< - DisplacementDim>::RowsAtCompileTime, - &LocalAssemblerIF::getIntPtEpsilon); - - add_secondary_variable("velocity", DisplacementDim, - &LocalAssemblerIF::getIntPtDarcyVelocity); - - add_secondary_variable("saturation", 1, - &LocalAssemblerIF::getIntPtSaturation); - - add_secondary_variable("micro_saturation", 1, - &LocalAssemblerIF::getIntPtMicroSaturation); - - add_secondary_variable("micro_pressure", 1, - &LocalAssemblerIF::getIntPtMicroPressure); - - add_secondary_variable("porosity", 1, &LocalAssemblerIF::getIntPtPorosity); - - add_secondary_variable("transport_porosity", 1, - &LocalAssemblerIF::getIntPtTransportPorosity); - - add_secondary_variable("dry_density_solid", 1, - &LocalAssemblerIF::getIntPtDryDensitySolid); - // // enable output of internal variables defined by material models // From a361f5bd09fd91983feca8806a9f9a70000ddb1e Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 12:35:38 +0200 Subject: [PATCH 11/14] [PL/RM] Moved solid material (state) to local assembler interface --- .../RichardsMechanics/IntegrationPointData.h | 29 +++----- .../LocalAssemblerInterface.h | 14 ++-- .../RichardsMechanicsFEM-impl.h | 73 +++++++++++-------- .../RichardsMechanics/RichardsMechanicsFEM.h | 18 +++-- 4 files changed, 70 insertions(+), 64 deletions(-) diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h index a1276108f9f..bef36fc89ee 100644 --- a/ProcessLib/RichardsMechanics/IntegrationPointData.h +++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h @@ -25,35 +25,24 @@ template struct IntegrationPointData final { - explicit IntegrationPointData( - MaterialLib::Solids::MechanicsBase const& - solid_material) - : solid_material(solid_material), - material_state_variables( - solid_material.createMaterialStateVariables()) - { - } - typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u; typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u; typename ShapeMatricesTypePressure::NodalRowVectorType N_p; typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p; - MaterialLib::Solids::MechanicsBase const& solid_material; - std::unique_ptr::MaterialStateVariables> - material_state_variables; double integration_weight = std::numeric_limits::quiet_NaN(); - void pushBackState() { material_state_variables->pushBackState(); } - MathLib::KelvinVector::KelvinMatrixType computeElasticTangentStiffness( double const t, ParameterLib::SpatialPosition const& x_position, double const dt, - double const temperature) + double const temperature, + MaterialLib::Solids::MechanicsBase const& + solid_material, + typename MaterialLib::Solids::MechanicsBase:: + MaterialStateVariables const& material_state_variables) { namespace MPL = MaterialPropertyLib; @@ -72,7 +61,7 @@ struct IntegrationPointData final auto&& solution = solid_material.integrateStress( variable_array_prev, variable_array, t, x_position, dt, - *material_state_variables); + material_state_variables); if (!solution) { @@ -104,7 +93,11 @@ struct IntegrationPointData final PrevState< ProcessLib::ThermoRichardsMechanics:: ConstitutiveStress_StrainTemperature::MechanicalStrainData< - DisplacementDim>> const& eps_m_prev) + DisplacementDim>> const& eps_m_prev, + MaterialLib::Solids::MechanicsBase const& + solid_material, + std::unique_ptr::MaterialStateVariables>& material_state_variables) { MaterialPropertyLib::VariableArray variable_array_prev; variable_array_prev.stress = sigma_eff_prev->sigma_eff; diff --git a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h index a9a57f99fdf..d9463fcc5fc 100644 --- a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h +++ b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h @@ -35,7 +35,10 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, : process_data_(process_data), integration_method_(integration_method), element_(e), - is_axially_symmetric_(is_axially_symmetric) + is_axially_symmetric_(is_axially_symmetric), + solid_material_(MaterialLib::Solids::selectSolidConstitutiveRelation( + process_data_.solid_materials, process_data_.material_ids, + e.getID())) { unsigned const n_integration_points = integration_method_.getNumberOfPoints(); @@ -46,15 +49,10 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, material_states_.reserve(n_integration_points); - auto const& solid_material = - MaterialLib::Solids::selectSolidConstitutiveRelation( - process_data_.solid_materials, process_data_.material_ids, - e.getID()); - for (unsigned ip = 0; ip < n_integration_points; ++ip) { material_states_.emplace_back( - solid_material.createMaterialStateVariables()); + solid_material_.createMaterialStateVariables()); // Set initial strain field to zero. std::get>(current_states_[ip]).eps = @@ -95,6 +93,8 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, MeshLib::Element const& element_; bool const is_axially_symmetric_; + MaterialLib::Solids::MechanicsBase const& solid_material_; + std::vector> current_states_; std::vector> prev_states_; std::vector< diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index f1bd168a25d..080271c44a2 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -132,7 +132,7 @@ RichardsMechanicsLocalAssemblerintegration_method_.getNumberOfPoints(); - _ip_data.reserve(n_integration_points); + _ip_data.resize(n_integration_points); _secondary_data.N_u.resize(n_integration_points); auto const shape_matrices_u = @@ -146,11 +146,6 @@ RichardsMechanicsLocalAssembler( e, is_axially_symmetric, this->integration_method_); - auto const& solid_material = - MaterialLib::Solids::selectSolidConstitutiveRelation( - this->process_data_.solid_materials, - this->process_data_.material_ids, e.getID()); - auto const& medium = this->process_data_.media_map.getMedium(this->element_.getID()); @@ -159,7 +154,6 @@ RichardsMechanicsLocalAssemblersolid_material_, + *this->material_states_[ip].material_state_variables); auto& eps = std::get>(this->current_states_[ip]) .eps; @@ -547,11 +542,11 @@ void RichardsMechanicsLocalAssembler< medium->property(MPL::PropertyType::biot_coefficient) .template value(variables, x_position, t, dt); auto const C_el = _ip_data[ip].computeElasticTangentStiffness( - t, x_position, dt, temperature); + t, x_position, dt, temperature, this->solid_material_, + *this->material_states_[ip].material_state_variables); - auto const beta_SR = - (1 - alpha) / - _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el); + auto const beta_SR = (1 - alpha) / this->solid_material_.getBulkModulus( + t, x_position, &C_el); variables.grain_compressibility = beta_SR; auto const rho_LR = @@ -715,7 +710,8 @@ void RichardsMechanicsLocalAssembler< } variables.equivalent_plastic_strain = - _ip_data[ip].material_state_variables->getEquivalentPlasticStrain(); + this->material_states_[ip] + .material_state_variables->getEquivalentPlasticStrain(); auto const K_intrinsic = MPL::formEigenTensor( medium->property(MPL::PropertyType::permeability) @@ -767,7 +763,8 @@ void RichardsMechanicsLocalAssembler< _ip_data[ip].updateConstitutiveRelation( variables, t, x_position, dt, temperature, sigma_eff, - sigma_eff_prev, eps_m, eps_m_prev); + sigma_eff_prev, eps_m, eps_m_prev, this->solid_material_, + this->material_states_[ip].material_state_variables); } // p_SR @@ -850,7 +847,11 @@ void RichardsMechanicsLocalAssembler& CD, StatefulData& SD, StatefulDataPrev const& SD_prev, - std::optional const& micro_porosity_parameters) + std::optional const& micro_porosity_parameters, + MaterialLib::Solids::MechanicsBase const& + solid_material, + ProcessLib::ThermoRichardsMechanics::MaterialStateData& + material_state_data) { auto const& liquid_phase = medium->phase("AqueousLiquid"); auto const& solid_phase = medium->phase("Solid"); @@ -876,11 +877,12 @@ void RichardsMechanicsLocalAssembler(variables, x_position, t, dt); *std::get(CD) = alpha; - auto const C_el = - ip_data.computeElasticTangentStiffness(t, x_position, dt, temperature); + auto const C_el = ip_data.computeElasticTangentStiffness( + t, x_position, dt, temperature, solid_material, + *material_state_data.material_state_variables); - auto const beta_SR = (1 - alpha) / ip_data.solid_material.getBulkModulus( - t, x_position, &C_el); + auto const beta_SR = + (1 - alpha) / solid_material.getBulkModulus(t, x_position, &C_el); variables.grain_compressibility = beta_SR; std::get(CD) .beta_SR = beta_SR; @@ -1057,7 +1059,8 @@ void RichardsMechanicsLocalAssemblergetEquivalentPlasticStrain(); + material_state_data.material_state_variables + ->getEquivalentPlasticStrain(); double const k_rel = medium->property(MPL::PropertyType::relative_permeability) @@ -1123,7 +1126,8 @@ void RichardsMechanicsLocalAssembler>(CD) = std::move(C); } @@ -1239,7 +1243,7 @@ void RichardsMechanicsLocalAssemblercurrent_states_[ip]; auto const& SD_prev = this->prev_states_[ip]; [[maybe_unused]] auto models = createConstitutiveModels( - this->process_data_, _ip_data[ip].solid_material); + this->process_data_, this->solid_material_); x_position.setIntegrationPoint(ip); auto const& w = _ip_data[ip].integration_weight; @@ -1285,7 +1289,8 @@ void RichardsMechanicsLocalAssembler{ p_cap_ip, p_cap_prev_ip, Eigen::Vector::Zero()}, - CD, SD, SD_prev, this->process_data_.micro_porosity_parameters); + CD, SD, SD_prev, this->process_data_.micro_porosity_parameters, + this->solid_material_, this->material_states_[ip]); { auto const& C = *std::get>(CD); @@ -1579,8 +1584,10 @@ std::vector RichardsMechanicsLocalAssembler< int const& n_components) const { return ProcessLib::getIntegrationPointDataMaterialStateVariables( - _ip_data, &IpData::material_state_variables, get_values_span, - n_components); + this->material_states_, + &ProcessLib::ThermoRichardsMechanics::MaterialStateData< + DisplacementDim>::material_state_variables, + get_values_span, n_components); } template (variables, x_position, t, dt); auto const C_el = _ip_data[ip].computeElasticTangentStiffness( - t, x_position, dt, temperature); + t, x_position, dt, temperature, this->solid_material_, + *this->material_states_[ip].material_state_variables); - auto const beta_SR = - (1 - alpha) / - _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el); + auto const beta_SR = (1 - alpha) / this->solid_material_.getBulkModulus( + t, x_position, &C_el); variables.grain_compressibility = beta_SR; variables.effective_pore_pressure = -chi_S_L * p_cap_ip; @@ -1867,7 +1874,8 @@ void RichardsMechanicsLocalAssemblergetEquivalentPlasticStrain(); + this->material_states_[ip] + .material_state_variables->getEquivalentPlasticStrain(); auto const K_intrinsic = MPL::formEigenTensor( medium->property(MPL::PropertyType::permeability) @@ -1933,7 +1941,8 @@ void RichardsMechanicsLocalAssemblersolid_material_, + this->material_states_[ip].material_state_variables); } auto const& b = this->process_data_.specific_body_force; @@ -1986,7 +1995,7 @@ RichardsMechanicsLocalAssembler:: getMaterialStateVariablesAt(unsigned integration_point) const { - return *_ip_data[integration_point].material_state_variables; + return *this->material_states_[integration_point].material_state_variables; } } // namespace RichardsMechanics } // namespace ProcessLib diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h index e00e0c864c9..f07fb69b937 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h @@ -154,10 +154,11 @@ class RichardsMechanicsLocalAssembler } double const t = 0; // TODO (naumov) pass t from top - ip_data.solid_material.initializeInternalStateVariables( - t, x_position, *ip_data.material_state_variables); + this->solid_material_.initializeInternalStateVariables( + t, x_position, + *this->material_states_[ip].material_state_variables); - ip_data.pushBackState(); + this->material_states_[ip].pushBackState(); this->prev_states_[ip] = SD; } @@ -171,9 +172,9 @@ class RichardsMechanicsLocalAssembler unsigned const n_integration_points = this->integration_method_.getNumberOfPoints(); - for (unsigned ip = 0; ip < n_integration_points; ip++) + for (auto& s : this->material_states_) { - _ip_data[ip].pushBackState(); + s.pushBackState(); } // TODO move to the local assembler interface @@ -282,8 +283,11 @@ class RichardsMechanicsLocalAssembler ConstitutiveData& CD, StatefulData& SD, StatefulDataPrev const& SD_prev, - std::optional const& - micro_porosity_parameters); + std::optional const& micro_porosity_parameters, + MaterialLib::Solids::MechanicsBase const& + solid_material, + ProcessLib::ThermoRichardsMechanics::MaterialStateData& + material_state_data); std::vector> _ip_data; From 3c1b7944ae21479f0490565e491d28625d805d3f Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 12:38:06 +0200 Subject: [PATCH 12/14] [PL/RM] Use reflection for IP data input --- .../LocalAssemblerInterface.h | 56 ++++++++- .../RichardsMechanicsFEM-impl.h | 119 ------------------ .../RichardsMechanics/RichardsMechanicsFEM.h | 6 - 3 files changed, 53 insertions(+), 128 deletions(-) diff --git a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h index d9463fcc5fc..a9f01c33a63 100644 --- a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h +++ b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h @@ -16,7 +16,9 @@ #include "NumLib/Extrapolation/ExtrapolatableElement.h" #include "NumLib/Fem/Integration/GenericIntegrationMethod.h" #include "ProcessLib/LocalAssemblerInterface.h" +#include "ProcessLib/Reflection/ReflectionSetIPData.h" #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/MaterialState.h" +#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h" #include "RichardsMechanicsProcessData.h" namespace ProcessLib @@ -60,9 +62,57 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, } } - virtual std::size_t setIPDataInitialConditions( - std::string_view const name, double const* values, - int const integration_order) = 0; + std::size_t setIPDataInitialConditions(std::string_view name, + double const* values, + int const integration_order) + { + if (integration_order != + static_cast(integration_method_.getIntegrationOrder())) + { + OGS_FATAL( + "Setting integration point initial conditions; The integration " + "order of the local assembler for element {:d} is different " + "from the integration order in the initial condition.", + element_.getID()); + } + + if (name == "sigma" && process_data_.initial_stress != nullptr) + { + OGS_FATAL( + "Setting initial conditions for stress from integration " + "point data and from a parameter '{:s}' is not possible " + "simultaneously.", + process_data_.initial_stress->name); + } + + if (name.starts_with("material_state_variable_")) + { + name.remove_prefix(24); + + auto const& internal_variables = + solid_material_.getInternalVariables(); + if (auto const iv = std::find_if( + begin(internal_variables), end(internal_variables), + [&name](auto const& iv) { return iv.name == name; }); + iv != end(internal_variables)) + { + DBUG("Setting material state variable '{:s}'", name); + return ProcessLib:: + setIntegrationPointDataMaterialStateVariables( + values, material_states_, + &ProcessLib::ThermoRichardsMechanics::MaterialStateData< + DisplacementDim>::material_state_variables, + iv->reference); + } + return 0; + } + + // TODO this logic could be pulled out of the local assembler into the + // process. That might lead to a slightly better performance due to less + // string comparisons. + return ProcessLib::Reflection::reflectSetIPData( + name, values, current_states_); + } virtual std::vector getMaterialStateVariableInternalState( std::function( diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index 080271c44a2..fbb5818088c 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -197,125 +197,6 @@ RichardsMechanicsLocalAssembler -std::size_t RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, - DisplacementDim>::setIPDataInitialConditions(std::string_view name, - double const* values, - int const integration_order) -{ - if (integration_order != - static_cast(this->integration_method_.getIntegrationOrder())) - { - OGS_FATAL( - "Setting integration point initial conditions; The integration " - "order of the local assembler for element {:d} is different " - "from the integration order in the initial condition.", - this->element_.getID()); - } - - if (name == "sigma") - { - if (this->process_data_.initial_stress != nullptr) - { - OGS_FATAL( - "Setting initial conditions for stress from integration " - "point data and from a parameter '{:s}' is not possible " - "simultaneously.", - this->process_data_.initial_stress->name); - } - return ProcessLib::setIntegrationPointKelvinVectorData( - values, this->current_states_, [](auto& tuple) -> auto& { - return std::get>( - tuple) - .sigma_eff; - }); - } - - if (name == "saturation") - { - return ProcessLib::setIntegrationPointScalarData( - values, this->current_states_, - [](auto& tuple) -> auto& - { - return std::get< - ProcessLib::ThermoRichardsMechanics::SaturationData>( - tuple) - .S_L; - }); - } - if (name == "porosity") - { - return ProcessLib::setIntegrationPointScalarData( - values, this->current_states_, - [](auto& tuple) -> auto& - { - return std::get< - ProcessLib::ThermoRichardsMechanics::PorosityData>( - tuple) - .phi; - }); - } - if (name == "transport_porosity") - { - return ProcessLib::setIntegrationPointScalarData( - values, this->current_states_, - [](auto& tuple) -> auto& - { - return std::get(tuple) - .phi; - }); - } - if (name == "swelling_stress") - { - return ProcessLib::setIntegrationPointKelvinVectorData( - values, this->current_states_, - [](auto& tuple) -> auto& - { - return std::get>( - tuple) - .sigma_sw; - }); - } - if (name == "epsilon") - { - return ProcessLib::setIntegrationPointKelvinVectorData( - values, this->current_states_, [](auto& tuple) -> auto& { - return std::get>(tuple).eps; - }); - } - if (name.starts_with("material_state_variable_")) - { - name.remove_prefix(24); - - // Using first ip data for solid material. TODO (naumov) move solid - // material into element, store only material state in IPs. - auto const& internal_variables = - _ip_data[0].solid_material.getInternalVariables(); - if (auto const iv = std::find_if( - begin(internal_variables), end(internal_variables), - [&name](auto const& iv) { return iv.name == name; }); - iv != end(internal_variables)) - { - DBUG("Setting material state variable '{:s}'", name); - return ProcessLib::setIntegrationPointDataMaterialStateVariables( - values, _ip_data, &IpData::material_state_variables, - iv->reference); - } - - ERR("Could not find variable {:s} in solid material model's internal " - "variables.", - name); - } - return 0; -} - template void RichardsMechanicsLocalAssembler& process_data); - /// \return the number of read integration points. - std::size_t setIPDataInitialConditions( - std::string_view const name, - double const* values, - int const integration_order) override; - void setInitialConditionsConcrete(Eigen::VectorXd const local_x, double const t, int const process_id) override; From 30608a97d350e7b485ba847af9787382367bed15 Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 12:47:53 +0200 Subject: [PATCH 13/14] [PL/SD] Moved methods to local assembler interface --- .../LocalAssemblerInterface.h | 49 ++++++++++++++++--- .../RichardsMechanicsFEM-impl.h | 48 ------------------ .../RichardsMechanics/RichardsMechanicsFEM.h | 34 ------------- 3 files changed, 43 insertions(+), 88 deletions(-) diff --git a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h index a9f01c33a63..fd8d870dafe 100644 --- a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h +++ b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h @@ -114,20 +114,38 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, name, values, current_states_); } - virtual std::vector getMaterialStateVariableInternalState( + std::vector getMaterialStateVariableInternalState( std::function( typename MaterialLib::Solids::MechanicsBase:: MaterialStateVariables&)> const& get_values_span, - int const& n_components) const = 0; + int const& n_components) const + { + return ProcessLib::getIntegrationPointDataMaterialStateVariables( + material_states_, + &ProcessLib::ThermoRichardsMechanics::MaterialStateData< + DisplacementDim>::material_state_variables, + get_values_span, n_components); + } // TODO move to NumLib::ExtrapolatableElement - virtual unsigned getNumberOfIntegrationPoints() const = 0; + unsigned getNumberOfIntegrationPoints() const + { + return integration_method_.getNumberOfPoints(); + } - virtual int getMaterialID() const = 0; + int getMaterialID() const + { + return process_data_.material_ids == nullptr + ? 0 + : (*process_data_.material_ids)[element_.getID()]; + } - virtual typename MaterialLib::Solids::MechanicsBase< + typename MaterialLib::Solids::MechanicsBase< DisplacementDim>::MaterialStateVariables const& - getMaterialStateVariablesAt(unsigned /*integration_point*/) const = 0; + getMaterialStateVariablesAt(unsigned integration_point) const + { + return *material_states_[integration_point].material_state_variables; + } static auto getReflectionDataForOutput() { @@ -137,6 +155,25 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, &Self::current_states_, &Self::output_data_); } + void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/, + Eigen::VectorXd const& /*local_x_prev*/, + double const /*t*/, double const /*dt*/, + int const /*process_id*/) override final + { + unsigned const n_integration_points = + integration_method_.getNumberOfPoints(); + + for (auto& s : material_states_) + { + s.pushBackState(); + } + + for (unsigned ip = 0; ip < n_integration_points; ip++) + { + prev_states_[ip] = current_states_[ip]; + } + } + protected: RichardsMechanicsProcessData& process_data_; NumLib::GenericIntegrationMethod const& integration_method_; diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index fbb5818088c..b225b114f76 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -1443,34 +1443,6 @@ void RichardsMechanicsLocalAssembler -int RichardsMechanicsLocalAssembler::getMaterialID() const -{ - return this->process_data_.material_ids == nullptr - ? 0 - : (*this->process_data_.material_ids)[this->element_.getID()]; -} - -template -std::vector RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getMaterialStateVariableInternalState( - std::function( - typename MaterialLib::Solids::MechanicsBase:: - MaterialStateVariables&)> const& get_values_span, - int const& n_components) const -{ - return ProcessLib::getIntegrationPointDataMaterialStateVariables( - this->material_states_, - &ProcessLib::ThermoRichardsMechanics::MaterialStateData< - DisplacementDim>::material_state_variables, - get_values_span, n_components); -} - template void RichardsMechanicsLocalAssembler(this->element_, this->is_axially_symmetric_, p_L, *this->process_data_.pressure_interpolated); } - -template -unsigned RichardsMechanicsLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, - DisplacementDim>::getNumberOfIntegrationPoints() const -{ - return this->integration_method_.getNumberOfPoints(); -} - -template -typename MaterialLib::Solids::MechanicsBase< - DisplacementDim>::MaterialStateVariables const& -RichardsMechanicsLocalAssembler:: - getMaterialStateVariablesAt(unsigned integration_point) const -{ - return *this->material_states_[integration_point].material_state_variables; -} } // namespace RichardsMechanics } // namespace ProcessLib diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h index 553bfb996ee..f11da1d504d 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h @@ -158,26 +158,6 @@ class RichardsMechanicsLocalAssembler } } - void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/, - Eigen::VectorXd const& /*local_x_prev*/, - double const /*t*/, double const /*dt*/, - int const /*process_id*/) override - { - unsigned const n_integration_points = - this->integration_method_.getNumberOfPoints(); - - for (auto& s : this->material_states_) - { - s.pushBackState(); - } - - // TODO move to the local assembler interface - for (unsigned ip = 0; ip < n_integration_points; ip++) - { - this->prev_states_[ip] = this->current_states_[ip]; - } - } - void computeSecondaryVariableConcrete( double const t, double const dt, Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev) override; @@ -191,14 +171,6 @@ class RichardsMechanicsLocalAssembler return Eigen::Map(N_u.data(), N_u.size()); } - int getMaterialID() const override; - - std::vector getMaterialStateVariableInternalState( - std::function( - typename MaterialLib::Solids::MechanicsBase:: - MaterialStateVariables&)> const& get_values_span, - int const& n_components) const override; - private: /** * Assemble local matrices and vectors arise from the linearized discretized @@ -261,12 +233,6 @@ class RichardsMechanicsLocalAssembler std::vector& local_K_data, std::vector& local_b_data, std::vector& local_Jac_data); - unsigned getNumberOfIntegrationPoints() const override; - - typename MaterialLib::Solids::MechanicsBase< - DisplacementDim>::MaterialStateVariables const& - getMaterialStateVariablesAt(unsigned integration_point) const override; - private: static void assembleWithJacobianEvalConstitutiveSetting( double const t, double const dt, From 58b0bff7e49a4f2a69f9ba31652830d3448af000 Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 11 Jun 2024 12:51:41 +0200 Subject: [PATCH 14/14] [PL/RM] Made two methods static --- .../RichardsMechanics/IntegrationPointData.h | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h index bef36fc89ee..d2dd04447eb 100644 --- a/ProcessLib/RichardsMechanics/IntegrationPointData.h +++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h @@ -33,16 +33,16 @@ struct IntegrationPointData final double integration_weight = std::numeric_limits::quiet_NaN(); - MathLib::KelvinVector::KelvinMatrixType - computeElasticTangentStiffness( - double const t, - ParameterLib::SpatialPosition const& x_position, - double const dt, - double const temperature, - MaterialLib::Solids::MechanicsBase const& - solid_material, - typename MaterialLib::Solids::MechanicsBase:: - MaterialStateVariables const& material_state_variables) + MathLib::KelvinVector:: + KelvinMatrixType static computeElasticTangentStiffness( + double const t, + ParameterLib::SpatialPosition const& x_position, + double const dt, + double const temperature, + MaterialLib::Solids::MechanicsBase const& + solid_material, + typename MaterialLib::Solids::MechanicsBase:: + MaterialStateVariables const& material_state_variables) { namespace MPL = MaterialPropertyLib; @@ -74,7 +74,7 @@ struct IntegrationPointData final return C; } - typename BMatricesType::KelvinMatrixType updateConstitutiveRelation( + static typename BMatricesType::KelvinMatrixType updateConstitutiveRelation( MaterialPropertyLib::VariableArray const& variable_array, double const t, ParameterLib::SpatialPosition const& x_position,