diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h index 7783b02fc7d..2641de6abc8 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h @@ -11,14 +11,19 @@ #pragma once #include "Density.h" +#include "DrySolidDensity.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/DarcyLaw.h" #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/LiquidViscosity.h" #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityData.h" #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" @@ -30,7 +35,15 @@ template using StatefulData = std::tuple< StrainData, ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature:: - EffectiveStressData>; + EffectiveStressData, + ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature:: + SwellingDataStateful, + ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature:: + MechanicalStrainData, + ProcessLib::ThermoRichardsMechanics::SaturationData, + ProcessLib::ThermoRichardsMechanics::PorosityData, + ProcessLib::ThermoRichardsMechanics::TransportPorosityData, MicroPressure, + MicroSaturation>; template using StatefulDataPrev = ProcessLib::ConstitutiveRelations::PrevStateOf< @@ -38,7 +51,9 @@ 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, + 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..c021158782c --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h @@ -0,0 +1,26 @@ +/** + * \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 + +#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 new file mode 100644 index 00000000000..96ace800113 --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h @@ -0,0 +1,25 @@ +/** + * \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 + +#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 new file mode 100644 index 00000000000..4048fdd58fe --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h @@ -0,0 +1,25 @@ +/** + * \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 + +#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/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h index 6363455b5f7..d2dd04447eb 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 { @@ -22,74 +25,24 @@ template struct IntegrationPointData final { - explicit IntegrationPointData( - MaterialLib::Solids::MechanicsBase const& - solid_material) - : solid_material(solid_material), - 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; typename ShapeMatricesTypePressure::NodalRowVectorType N_p; typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p; - 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 = 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(); - double dry_density_pellet_unsaturated = - std::numeric_limits::quiet_NaN(); - - MaterialLib::Solids::MechanicsBase const& solid_material; - std::unique_ptr::MaterialStateVariables> - material_state_variables; double integration_weight = std::numeric_limits::quiet_NaN(); - void pushBackState() - { - eps_m_prev = eps_m; - sigma_sw_prev = sigma_sw; - 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(); - } - - MathLib::KelvinVector::KelvinMatrixType - computeElasticTangentStiffness( - double const t, - ParameterLib::SpatialPosition const& x_position, - double const dt, - double const temperature) + 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; @@ -108,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) { @@ -121,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, @@ -132,13 +85,25 @@ 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, + 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; 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/LocalAssemblerInterface.h b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h index 529002d1f65..fd8d870dafe 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 @@ -35,7 +37,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 +51,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 = @@ -62,100 +62,117 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, } } - virtual std::size_t setIPDataInitialConditions( - 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; + 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()); + } - virtual std::vector getTransportPorosity() const = 0; + 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); + } - virtual std::vector const& getIntPtTransportPorosity( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; + 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; + } - virtual std::vector const& getIntPtDryDensitySolid( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 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::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() + { + using Self = LocalAssemblerInterface; + + return ProcessLib::Reflection::reflectWithoutName( + &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_; @@ -163,6 +180,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 13728943eda..b225b114f76 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, @@ -41,19 +41,29 @@ 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, + PrevState const + phi_M_prev, + 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( 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 +72,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; } } @@ -77,39 +88,31 @@ 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 p_L_m_prev = ip_data.liquid_pressure_m_prev; - - auto const S_L_m_prev = ip_data.saturation_m_prev; + double const phi_m_prev = phi_prev->phi - phi_M_prev->phi; 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)); - 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; + *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_prev + delta_sigma_sw; + sigma_sw.sigma_sw = sigma_sw_prev->sigma_sw + delta_sigma_sw; } } @@ -129,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 = @@ -143,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()); @@ -156,7 +154,6 @@ 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, @@ -192,97 +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, _ip_data, - &IpData::saturation); - } - if (name == "porosity") - { - return ProcessLib::setIntegrationPointScalarData(values, _ip_data, - &IpData::porosity); - } - if (name == "transport_porosity") - { - return ProcessLib::setIntegrationPointScalarData( - values, _ip_data, &IpData::transport_porosity); - } - if (name == "swelling_stress") - { - return ProcessLib::setIntegrationPointKelvinVectorData( - values, _ip_data, &IpData::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(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) .template value(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)) { MPL::VariableArray vars; vars.capillary_pressure = p_cap_ip; - double const S_L_m = - medium->property(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 // restart. 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& eps = std::get>(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,11 +390,17 @@ 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; - 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); @@ -477,11 +423,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 = @@ -527,9 +473,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); @@ -546,8 +498,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. @@ -573,14 +534,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 { @@ -595,7 +565,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::getgetEquivalentPlasticStrain(); + this->material_states_[ip] + .material_state_variables->getEquivalentPlasticStrain(); auto const K_intrinsic = MPL::formEigenTensor( medium->property(MPL::PropertyType::permeability) @@ -628,13 +604,21 @@ void RichardsMechanicsLocalAssembler< // // displacement equation, displacement part // - eps_m.noalias() = - solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate) - ? eps.eps + C_el.inverse() * sigma_sw - : eps.eps; - variables.mechanical_strain - .emplace>( + { + 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 +632,20 @@ 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, this->solid_material_, + this->material_states_[ip].material_state_variables); } // p_SR @@ -733,7 +727,12 @@ void RichardsMechanicsLocalAssembler const& p_cap_data, ConstitutiveData& CD, StatefulData& SD, - StatefulDataPrev const& SD_prev) + StatefulDataPrev const& SD_prev, + 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"); @@ -747,19 +746,24 @@ 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& 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); *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; @@ -824,10 +828,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); @@ -837,11 +846,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) @@ -849,23 +865,57 @@ 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); + 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); + 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( + *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)) { 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 @@ -890,7 +940,8 @@ void RichardsMechanicsLocalAssemblergetEquivalentPlasticStrain(); + material_state_data.material_state_variables + ->getEquivalentPlasticStrain(); double const k_rel = medium->property(MPL::PropertyType::relative_permeability) @@ -912,15 +963,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 +995,20 @@ 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, solid_material, + material_state_data.material_state_variables); *std::get>(CD) = std::move(C); } @@ -1051,7 +1124,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; @@ -1097,7 +1170,8 @@ 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, + this->solid_material_, this->material_states_[ip]); { auto const& C = *std::get>(CD); @@ -1191,7 +1265,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 + @@ -1307,7 +1388,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; @@ -1318,7 +1400,8 @@ void RichardsMechanicsLocalAssembler>( + this->prev_states_[ip]); local_Jac .template block( pressure_index, pressure_index) @@ -1360,307 +1443,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 = _ip_data[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::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( - _ip_data, &IpData::material_state_variables, get_values_span, - 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() = _ip_data[ip].v_darcy; - } - - 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( - _ip_data, &IpData::saturation, 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( - _ip_data, &IpData::saturation_m, 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( - _ip_data, &IpData::liquid_pressure_m, 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(_ip_data, - &IpData::porosity, 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( - _ip_data, &IpData::transport_porosity, 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( - _ip_data, &IpData::dry_density_solid, cache); -} - template 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; + 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; @@ -1818,11 +1606,11 @@ void RichardsMechanicsLocalAssembler(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; @@ -1832,9 +1620,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); @@ -1849,24 +1643,64 @@ 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]); + 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]); + 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( + *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, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m); + } 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 @@ -1893,7 +1727,8 @@ void RichardsMechanicsLocalAssemblergetEquivalentPlasticStrain(); + this->material_states_[ip] + .material_state_variables->getEquivalentPlasticStrain(); auto const K_intrinsic = MPL::formEigenTensor( medium->property(MPL::PropertyType::permeability) @@ -1912,17 +1747,28 @@ 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; + *std::get(this->output_data_[ip]) = (1 - phi) * rho_SR; - 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,18 +1782,30 @@ 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, this->solid_material_, + this->material_states_[ip].material_state_variables); } auto const& b = this->process_data_.specific_body_force; // Compute the velocity auto const& dNdx_p = _ip_data[ip].dNdx_p; - _ip_data[ip].v_darcy.noalias() = - -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b; + std::get< + ProcessLib::ThermoRichardsMechanics::DarcyLawData>( + 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; @@ -1972,25 +1830,5 @@ 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 *_ip_data[integration_point].material_state_variables; -} } // namespace RichardsMechanics } // namespace ProcessLib diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h index 16b82107462..f11da1d504d 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h @@ -90,12 +90,6 @@ class RichardsMechanicsLocalAssembler bool const is_axially_symmetric, RichardsMechanicsProcessData& 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; @@ -154,35 +148,16 @@ 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; } } - 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 (unsigned ip = 0; ip < n_integration_points; ip++) - { - _ip_data[ip].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; @@ -196,83 +171,6 @@ class RichardsMechanicsLocalAssembler return Eigen::Map(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( - std::function( - typename MaterialLib::Solids::MechanicsBase:: - 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 @@ -335,14 +233,8 @@ 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: - 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 +242,12 @@ 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, + MaterialLib::Solids::MechanicsBase const& + solid_material, + ProcessLib::ThermoRichardsMechanics::MaterialStateData& + material_state_data); std::vector> _ip_data; 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 //