From fdf4056488f163d2e9161333d8178cbbd09de61c Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 22 Sep 2023 17:35:46 +0200 Subject: [PATCH 01/10] [PL/SD] Move elem.-type independent members to LA Especially the solid_material is now stored once per element and not in every integration point. --- .../LocalAssemblerInterface.h | 24 ++++ .../SmallDeformation/SmallDeformationFEM.h | 111 ++++++++---------- 2 files changed, 71 insertions(+), 64 deletions(-) diff --git a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h index bb525d62d31..4459c25070a 100644 --- a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h +++ b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h @@ -13,9 +13,12 @@ #include #include "MaterialLib/SolidModels/MechanicsBase.h" +#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h" #include "NumLib/Extrapolation/ExtrapolatableElement.h" +#include "NumLib/Fem/Integration/GenericIntegrationMethod.h" #include "ProcessLib/Deformation/MaterialForces.h" #include "ProcessLib/LocalAssemblerInterface.h" +#include "SmallDeformationProcessData.h" namespace ProcessLib { @@ -27,6 +30,20 @@ struct SmallDeformationLocalAssemblerInterface public MaterialForcesInterface, public NumLib::ExtrapolatableElement { + SmallDeformationLocalAssemblerInterface( + MeshLib::Element const& e, + NumLib::GenericIntegrationMethod const& integration_method, + bool const is_axially_symmetric, + SmallDeformationProcessData& process_data) + : process_data_(process_data), + integration_method_(integration_method), + element_(e), + is_axially_symmetric_(is_axially_symmetric), + solid_material_(MaterialLib::Solids::selectSolidConstitutiveRelation( + process_data_.solid_materials, process_data_.material_ids, + element_.getID())) + { + } virtual std::size_t setIPDataInitialConditions( std::string const& name, double const* values, int const integration_order) = 0; @@ -64,6 +81,13 @@ struct SmallDeformationLocalAssemblerInterface virtual typename MaterialLib::Solids::MechanicsBase< DisplacementDim>::MaterialStateVariables const& getMaterialStateVariablesAt(unsigned /*integration_point*/) const = 0; + +protected: + SmallDeformationProcessData& process_data_; + NumLib::GenericIntegrationMethod const& integration_method_; + MeshLib::Element const& element_; + bool const is_axially_symmetric_; + MaterialLib::Solids::MechanicsBase const& solid_material_; }; } // namespace SmallDeformation diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h index d69d1ac738d..79e50dec43d 100644 --- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h +++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h @@ -17,14 +17,12 @@ #include "LocalAssemblerInterface.h" #include "MaterialLib/PhysicalConstant.h" -#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h" #include "MathLib/EigenBlockMatrixView.h" #include "MathLib/LinAlg/Eigen/EigenMapTools.h" #include "NumLib/DOF/LocalDOF.h" #include "NumLib/Extrapolation/ExtrapolatableElement.h" #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h" #include "NumLib/Fem/InitShapeMatrices.h" -#include "NumLib/Fem/Integration/GenericIntegrationMethod.h" #include "NumLib/Fem/ShapeMatrixPolicy.h" #include "ParameterLib/Parameter.h" #include "ProcessLib/Deformation/BMatrixPolicy.h" @@ -33,7 +31,6 @@ #include "ProcessLib/LocalAssemblerInterface.h" #include "ProcessLib/LocalAssemblerTraits.h" #include "ProcessLib/Utils/SetOrGetIntegrationPointData.h" -#include "SmallDeformationProcessData.h" namespace ProcessLib { @@ -48,8 +45,7 @@ struct IntegrationPointData final explicit IntegrationPointData( MaterialLib::Solids::MechanicsBase const& solid_material) - : solid_material(solid_material), - material_state_variables( + : material_state_variables( solid_material.createMaterialStateVariables()) { } @@ -58,7 +54,6 @@ struct IntegrationPointData final typename BMatricesType::KelvinVectorType eps; double free_energy_density = 0; - MaterialLib::Solids::MechanicsBase const& solid_material; std::unique_ptr::MaterialStateVariables> material_state_variables; @@ -121,35 +116,27 @@ class SmallDeformationLocalAssembler NumLib::GenericIntegrationMethod const& integration_method, bool const is_axially_symmetric, SmallDeformationProcessData& process_data) - : _process_data(process_data), - _integration_method(integration_method), - _element(e), - _is_axially_symmetric(is_axially_symmetric) + : SmallDeformationLocalAssemblerInterface( + e, integration_method, is_axially_symmetric, process_data) { unsigned const n_integration_points = - _integration_method.getNumberOfPoints(); + this->integration_method_.getNumberOfPoints(); _ip_data.reserve(n_integration_points); _secondary_data.N.resize(n_integration_points); auto const shape_matrices = NumLib::initShapeMatrices(e, is_axially_symmetric, - _integration_method); - - auto& solid_material = - MaterialLib::Solids::selectSolidConstitutiveRelation( - _process_data.solid_materials, - _process_data.material_ids, - e.getID()); + DisplacementDim>( + e, is_axially_symmetric, this->integration_method_); for (unsigned ip = 0; ip < n_integration_points; ip++) { - _ip_data.emplace_back(solid_material); + _ip_data.emplace_back(this->solid_material_); auto& ip_data = _ip_data[ip]; auto const& sm = shape_matrices[ip]; _ip_data[ip].integration_weight = - _integration_method.getWeightedPoint(ip).getWeight() * + this->integration_method_.getWeightedPoint(ip).getWeight() * sm.integralMeasure * sm.detJ; ip_data.N = sm.N; @@ -175,24 +162,24 @@ class SmallDeformationLocalAssembler int const integration_order) override { if (integration_order != - static_cast(_integration_method.getIntegrationOrder())) + 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.", - _element.getID()); + this->element_.getID()); } if (name == "sigma_ip") { - if (_process_data.initial_stress != nullptr) + 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.", - _process_data.initial_stress->name); + this->process_data_.initial_stress->name); } return setSigma(values); } @@ -203,10 +190,8 @@ class SmallDeformationLocalAssembler name.substr(24, name.size() - 24 - 3); DBUG("Setting material state variable '{:s}'", variable_name); - // 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(); + this->solid_material_.getInternalVariables(); if (auto const iv = std::find_if( begin(internal_variables), end(internal_variables), [&variable_name](auto const& iv) @@ -231,31 +216,31 @@ class SmallDeformationLocalAssembler void initializeConcrete() override { unsigned const n_integration_points = - _integration_method.getNumberOfPoints(); + this->integration_method_.getNumberOfPoints(); for (unsigned ip = 0; ip < n_integration_points; ip++) { auto& ip_data = _ip_data[ip]; ParameterLib::SpatialPosition const x_position{ - std::nullopt, _element.getID(), ip, + std::nullopt, this->element_.getID(), ip, MathLib::Point3d( NumLib::interpolateCoordinates( - _element, ip_data.N))}; + this->element_, ip_data.N))}; /// Set initial stress from parameter. - if (_process_data.initial_stress != nullptr) + if (this->process_data_.initial_stress != nullptr) { ip_data.sigma = MathLib::KelvinVector::symmetricTensorToKelvinVector< - DisplacementDim>((*_process_data.initial_stress)( + DisplacementDim>((*this->process_data_.initial_stress)( std::numeric_limits< double>::quiet_NaN() /* time independent */, x_position)); } double const t = 0; // TODO (naumov) pass t from top - ip_data.solid_material.initializeInternalStateVariables( + this->solid_material_.initializeInternalStateVariables( t, x_position, *ip_data.material_state_variables); ip_data.pushBackState(); @@ -274,7 +259,7 @@ class SmallDeformationLocalAssembler ip_data) const { auto const& solid_phase = - this->_process_data.media_map.getMedium(this->_element.getID()) + this->process_data_.media_map.getMedium(this->element_.getID()) ->phase("Solid"); MPL::VariableArray variables_prev; @@ -290,12 +275,12 @@ class SmallDeformationLocalAssembler auto const& dNdx = ip_data.dNdx; auto const x_coord = NumLib::interpolateXCoordinate( - _element, N); + this->element_, N); auto const B = LinearBMatrix::computeBMatrix( - dNdx, N, x_coord, _is_axially_symmetric); + dNdx, N, x_coord, this->is_axially_symmetric_); eps.noalias() = B * u; @@ -307,8 +292,8 @@ class SmallDeformationLocalAssembler B * u_prev); double const T_ref = - _process_data.reference_temperature - ? (*_process_data.reference_temperature)(t, x_position)[0] + this->process_data_.reference_temperature + ? (*this->process_data_.reference_temperature)(t, x_position)[0] : std::numeric_limits::quiet_NaN(); variables_prev.temperature = T_ref; @@ -321,7 +306,7 @@ class SmallDeformationLocalAssembler solid_phase[MPL::PropertyType::density].template value( variables, x_position, t, dt); - auto&& solution = ip_data.solid_material.integrateStress( + auto&& solution = this->solid_material_.integrateStress( variables_prev, variables, t, x_position, dt, *state); if (!solution) @@ -367,12 +352,12 @@ class SmallDeformationLocalAssembler auto [u_prev] = localDOF(local_x_prev); unsigned const n_integration_points = - _integration_method.getNumberOfPoints(); + this->integration_method_.getNumberOfPoints(); ParameterLib::SpatialPosition x_position; - x_position.setElementID(_element.getID()); + x_position.setElementID(this->element_.getID()); - auto const& b = _process_data.specific_body_force; + auto const& b = this->process_data_.specific_body_force; for (unsigned ip = 0; ip < n_integration_points; ip++) { @@ -383,11 +368,12 @@ class SmallDeformationLocalAssembler auto const x_coord = NumLib::interpolateXCoordinate(_element, N); + ShapeMatricesType>( + this->element_, N); auto const B = LinearBMatrix::computeBMatrix< DisplacementDim, ShapeFunction::NPOINTS, - typename BMatricesType::BMatrixType>(dNdx, N, x_coord, - _is_axially_symmetric); + typename BMatricesType::BMatrixType>( + dNdx, N, x_coord, this->is_axially_symmetric_); auto const& sigma = _ip_data[ip].sigma; @@ -407,10 +393,10 @@ class SmallDeformationLocalAssembler int const /*process_id*/) override { unsigned const n_integration_points = - _integration_method.getNumberOfPoints(); + this->integration_method_.getNumberOfPoints(); ParameterLib::SpatialPosition x_position; - x_position.setElementID(_element.getID()); + x_position.setElementID(this->element_.getID()); for (unsigned ip = 0; ip < n_integration_points; ip++) { @@ -425,7 +411,7 @@ class SmallDeformationLocalAssembler // Update free energy density needed for material forces. _ip_data[ip].free_energy_density = - _ip_data[ip].solid_material.computeFreeEnergyDensity( + this->solid_material_.computeFreeEnergyDensity( t, x_position, dt, eps, sigma, *state); _ip_data[ip].pushBackState(); @@ -440,8 +426,9 @@ class SmallDeformationLocalAssembler DisplacementDim, ShapeFunction, ShapeMatricesType, typename BMatricesType::NodalForceVectorType, NodalDisplacementVectorType, GradientVectorType, - GradientMatrixType>(local_x, nodal_values, _integration_method, - _ip_data, _element, _is_axially_symmetric); + GradientMatrixType>(local_x, nodal_values, + this->integration_method_, _ip_data, + this->element_, this->is_axially_symmetric_); } Eigen::Map getShapeMatrix( @@ -506,14 +493,15 @@ class SmallDeformationLocalAssembler unsigned getNumberOfIntegrationPoints() const override { - return _integration_method.getNumberOfPoints(); + return this->integration_method_.getNumberOfPoints(); } int getMaterialID() const override { - return _process_data.material_ids == nullptr + return this->process_data_.material_ids == nullptr ? 0 - : (*_process_data.material_ids)[_element.getID()]; + : (*this->process_data_ + .material_ids)[this->element_.getID()]; } std::vector getMaterialStateVariableInternalState( @@ -538,11 +526,11 @@ class SmallDeformationLocalAssembler double const /*t*/, double const /*dt*/, Eigen::VectorXd const& /*x*/, Eigen::VectorXd const& /*x_prev*/) override { - int const elem_id = _element.getID(); + int const elem_id = this->element_.getID(); ParameterLib::SpatialPosition x_position; x_position.setElementID(elem_id); unsigned const n_integration_points = - _integration_method.getNumberOfPoints(); + this->integration_method_.getNumberOfPoints(); auto sigma_sum = MathLib::KelvinVector::tensorToKelvin( Eigen::Matrix::Zero()); @@ -562,7 +550,7 @@ class SmallDeformationLocalAssembler sigma_avg); Eigen::Map( - &(*_process_data.principal_stress_values)[elem_id * 3], 3) = + &(*this->process_data_.principal_stress_values)[elem_id * 3], 3) = e_s.eigenvalues(); auto eigen_vectors = e_s.eigenvectors(); @@ -570,8 +558,8 @@ class SmallDeformationLocalAssembler for (auto i = 0; i < 3; i++) { Eigen::Map( - &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) = - eigen_vectors.col(i); + &(*this->process_data_.principal_stress_vector[i])[elem_id * 3], + 3) = eigen_vectors.col(i); } } @@ -583,14 +571,9 @@ class SmallDeformationLocalAssembler } private: - SmallDeformationProcessData& _process_data; - std::vector> _ip_data; - NumLib::GenericIntegrationMethod const& _integration_method; - MeshLib::Element const& _element; SecondaryData _secondary_data; - bool const _is_axially_symmetric; }; } // namespace SmallDeformation From bdb94b63e3ca91aa40dd1888627af34bb3b0b5e1 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 22 Sep 2023 18:41:12 +0200 Subject: [PATCH 02/10] [PL/SD] Move solid material states in LA Separating the solid material states from integration point data. Also move all but one functions into the local assembler, since they are no longer accessing the integration point data. --- .../ConstitutiveRelations/MaterialState.h | 32 ++++++++ .../LocalAssemblerInterface.h | 51 ++++++++++-- .../SmallDeformation/SmallDeformationFEM.h | 82 +++++-------------- 3 files changed, 95 insertions(+), 70 deletions(-) create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/MaterialState.h diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/MaterialState.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/MaterialState.h new file mode 100644 index 00000000000..d0ddb14fb9b --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/MaterialState.h @@ -0,0 +1,32 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 "MaterialLib/SolidModels/MechanicsBase.h" + +namespace ProcessLib::SmallDeformation +{ +template +class MaterialStateData +{ + using MSV = typename MaterialLib::Solids::MechanicsBase< + DisplacementDim>::MaterialStateVariables; + +public: + explicit MaterialStateData(std::unique_ptr&& material_state_variables) + : material_state_variables(std::move(material_state_variables)) + { + } + + void pushBackState() { material_state_variables->pushBackState(); } + + std::unique_ptr material_state_variables; +}; +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h index 4459c25070a..9954f7afc01 100644 --- a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h +++ b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h @@ -12,12 +12,14 @@ #include +#include "ConstitutiveRelations/MaterialState.h" #include "MaterialLib/SolidModels/MechanicsBase.h" #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h" #include "NumLib/Extrapolation/ExtrapolatableElement.h" #include "NumLib/Fem/Integration/GenericIntegrationMethod.h" #include "ProcessLib/Deformation/MaterialForces.h" #include "ProcessLib/LocalAssemblerInterface.h" +#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h" #include "SmallDeformationProcessData.h" namespace ProcessLib @@ -43,6 +45,15 @@ struct SmallDeformationLocalAssemblerInterface process_data_.solid_materials, process_data_.material_ids, element_.getID())) { + unsigned const n_integration_points = + integration_method_.getNumberOfPoints(); + + material_states_.reserve(n_integration_points); + for (unsigned ip = 0; ip < n_integration_points; ++ip) + { + material_states_.emplace_back( + solid_material_.createMaterialStateVariables()); + } } virtual std::size_t setIPDataInitialConditions( std::string const& name, double const* values, @@ -67,23 +78,45 @@ struct SmallDeformationLocalAssemblerInterface std::vector const& dof_table, std::vector& cache) const = 0; - virtual std::vector getMaterialStateVariableInternalState( + // TODO move to NumLib::ExtrapolatableElement + unsigned getNumberOfIntegrationPoints() const + { + return integration_method_.getNumberOfPoints(); + } + + int getMaterialID() const + { + return process_data_.material_ids == nullptr + ? 0 + : (*process_data_.material_ids)[element_.getID()]; + } + + std::vector getMaterialStateVariableInternalState( std::function( typename MaterialLib::Solids::MechanicsBase:: MaterialStateVariables&)> const& get_values_span, - int const& n_components) const = 0; - - // TODO move to NumLib::ExtrapolatableElement - virtual unsigned getNumberOfIntegrationPoints() const = 0; - - virtual int getMaterialID() const = 0; + int const& n_components) const + { + return ProcessLib::getIntegrationPointDataMaterialStateVariables( + material_states_, + &MaterialStateData::material_state_variables, + get_values_span, n_components); + } - 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; + } protected: SmallDeformationProcessData& process_data_; + + // Material state is special, because it contains both the current and the + // old state. + std::vector> material_states_; + NumLib::GenericIntegrationMethod const& integration_method_; MeshLib::Element const& element_; bool const is_axially_symmetric_; diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h index 79e50dec43d..c01136a7b84 100644 --- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h +++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h @@ -30,7 +30,6 @@ #include "ProcessLib/Deformation/LinearBMatrix.h" #include "ProcessLib/LocalAssemblerInterface.h" #include "ProcessLib/LocalAssemblerTraits.h" -#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h" namespace ProcessLib { @@ -42,31 +41,15 @@ template struct IntegrationPointData final { - explicit IntegrationPointData( - MaterialLib::Solids::MechanicsBase const& - solid_material) - : material_state_variables( - solid_material.createMaterialStateVariables()) - { - } - typename BMatricesType::KelvinVectorType sigma, sigma_prev; typename BMatricesType::KelvinVectorType eps; double free_energy_density = 0; - std::unique_ptr::MaterialStateVariables> - material_state_variables; - double integration_weight; typename ShapeMatricesType::NodalRowVectorType N; typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx; - void pushBackState() - { - sigma_prev = sigma; - material_state_variables->pushBackState(); - } + void pushBackState() { sigma_prev = sigma; } EIGEN_MAKE_ALIGNED_OPERATOR_NEW; }; @@ -122,7 +105,7 @@ class SmallDeformationLocalAssembler unsigned const n_integration_points = this->integration_method_.getNumberOfPoints(); - _ip_data.reserve(n_integration_points); + _ip_data.resize(n_integration_points); _secondary_data.N.resize(n_integration_points); auto const shape_matrices = @@ -132,7 +115,6 @@ class SmallDeformationLocalAssembler for (unsigned ip = 0; ip < n_integration_points; ip++) { - _ip_data.emplace_back(this->solid_material_); auto& ip_data = _ip_data[ip]; auto const& sm = shape_matrices[ip]; _ip_data[ip].integration_weight = @@ -200,7 +182,9 @@ class SmallDeformationLocalAssembler { return ProcessLib:: setIntegrationPointDataMaterialStateVariables( - values, _ip_data, &IpData::material_state_variables, + values, this->material_states_, + &MaterialStateData< + DisplacementDim>::material_state_variables, iv->reference); } @@ -240,10 +224,12 @@ class SmallDeformationLocalAssembler } double const t = 0; // TODO (naumov) pass t from top + auto& material_state = this->material_states_[ip]; this->solid_material_.initializeInternalStateVariables( - t, x_position, *ip_data.material_state_variables); + t, x_position, *material_state.material_state_variables); ip_data.pushBackState(); + material_state.pushBackState(); } } @@ -256,7 +242,8 @@ class SmallDeformationLocalAssembler ParameterLib::SpatialPosition const& x_position, double const t, double const dt, IntegrationPointData& - ip_data) const + ip_data, + MaterialStateData& material_state) const { auto const& solid_phase = this->process_data_.media_map.getMedium(this->element_.getID()) @@ -269,7 +256,6 @@ class SmallDeformationLocalAssembler auto& eps = ip_data.eps; auto& sigma = ip_data.sigma; - auto& state = ip_data.material_state_variables; auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; @@ -307,7 +293,8 @@ class SmallDeformationLocalAssembler variables, x_position, t, dt); auto&& solution = this->solid_material_.integrateStress( - variables_prev, variables, t, x_position, dt, *state); + variables_prev, variables, t, x_position, dt, + *material_state.material_state_variables); if (!solution) { @@ -315,7 +302,8 @@ class SmallDeformationLocalAssembler } MathLib::KelvinVector::KelvinMatrixType C; - std::tie(sigma, state, C) = std::move(*solution); + std::tie(sigma, material_state.material_state_variables, C) = + std::move(*solution); return {C, rho}; } @@ -378,7 +366,8 @@ class SmallDeformationLocalAssembler auto const& sigma = _ip_data[ip].sigma; auto const [C, rho] = updateConstitutiveRelations( - u, u_prev, x_position, t, dt, _ip_data[ip]); + u, u_prev, x_position, t, dt, _ip_data[ip], + this->material_states_[ip]); local_b.noalias() -= (B.transpose() * sigma - N_u_op(N).transpose() * rho * b) * w; @@ -403,18 +392,20 @@ class SmallDeformationLocalAssembler x_position.setIntegrationPoint(ip); updateConstitutiveRelations(local_x, local_x_prev, x_position, t, - dt, _ip_data[ip]); + dt, _ip_data[ip], + this->material_states_[ip]); auto& eps = _ip_data[ip].eps; auto& sigma = _ip_data[ip].sigma; - auto& state = _ip_data[ip].material_state_variables; + auto& state = *this->material_states_[ip].material_state_variables; // Update free energy density needed for material forces. _ip_data[ip].free_energy_density = this->solid_material_.computeFreeEnergyDensity( - t, x_position, dt, eps, sigma, *state); + t, x_position, dt, eps, sigma, state); _ip_data[ip].pushBackState(); + this->material_states_[ip].pushBackState(); } } @@ -491,37 +482,6 @@ class SmallDeformationLocalAssembler _ip_data, &IpData::eps, cache); } - unsigned getNumberOfIntegrationPoints() const override - { - return this->integration_method_.getNumberOfPoints(); - } - - int getMaterialID() const override - { - return this->process_data_.material_ids == nullptr - ? 0 - : (*this->process_data_ - .material_ids)[this->element_.getID()]; - } - - std::vector getMaterialStateVariableInternalState( - std::function( - typename MaterialLib::Solids::MechanicsBase:: - MaterialStateVariables&)> const& get_values_span, - int const& n_components) const override - { - return ProcessLib::getIntegrationPointDataMaterialStateVariables( - _ip_data, &IpData::material_state_variables, get_values_span, - n_components); - } - - typename MaterialLib::Solids::MechanicsBase< - DisplacementDim>::MaterialStateVariables const& - getMaterialStateVariablesAt(unsigned integration_point) const override - { - return *_ip_data[integration_point].material_state_variables; - } - void computeSecondaryVariableConcrete( double const /*t*/, double const /*dt*/, Eigen::VectorXd const& /*x*/, Eigen::VectorXd const& /*x_prev*/) override From 4a4af87e717ceb067c8bb45a6955360bced8e464 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 22 Sep 2023 15:38:15 +0200 Subject: [PATCH 03/10] [PL/SD] Add ConstitutiveRelations; not yet used --- ProcessLib/SmallDeformation/CMakeLists.txt | 1 + .../ConstitutiveRelations/Base.h | 120 ++++++++++++++++++ .../ConstitutiveRelations/ConstitutiveData.h | 81 ++++++++++++ .../ConstitutiveModels.h | 38 ++++++ .../ConstitutiveSetting.cpp | 62 +++++++++ .../ConstitutiveSetting.h | 39 ++++++ .../CreateConstitutiveSetting.cpp | 33 +++++ .../CreateConstitutiveSetting.h | 46 +++++++ .../ConstitutiveRelations/Gravity.cpp | 26 ++++ .../ConstitutiveRelations/Gravity.h | 41 ++++++ .../ConstitutiveRelations/Invoke.h | 54 ++++++++ .../ConstitutiveRelations/SolidDensity.cpp | 31 +++++ .../ConstitutiveRelations/SolidDensity.h | 31 +++++ .../ConstitutiveRelations/SolidMechanics.cpp | 67 ++++++++++ .../ConstitutiveRelations/SolidMechanics.h | 61 +++++++++ .../ConstitutiveRelations/StressData.h | 32 +++++ 16 files changed, 763 insertions(+) create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/Base.h create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveModels.h create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.h create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/CreateConstitutiveSetting.cpp create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/CreateConstitutiveSetting.h create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.cpp create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.h create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/Invoke.h create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.cpp create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.h create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/StressData.h diff --git a/ProcessLib/SmallDeformation/CMakeLists.txt b/ProcessLib/SmallDeformation/CMakeLists.txt index 37d9043240d..72f111ef3e8 100644 --- a/ProcessLib/SmallDeformation/CMakeLists.txt +++ b/ProcessLib/SmallDeformation/CMakeLists.txt @@ -1,4 +1,5 @@ get_source_files(SOURCES) +append_source_files(SOURCES ConstitutiveRelations) ogs_add_library(SmallDeformation ${SOURCES}) target_link_libraries(SmallDeformation PUBLIC ProcessLib PRIVATE ParameterLib) diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/Base.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/Base.h new file mode 100644 index 00000000000..b92cb296b4b --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/Base.h @@ -0,0 +1,120 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 "MaterialLib/MPL/Medium.h" +#include "MathLib/KelvinVector.h" +#include "ParameterLib/SpatialPosition.h" +#include "ProcessLib/Reflection/ReflectionData.h" + +namespace ProcessLib::SmallDeformation +{ +template +using KelvinVector = MathLib::KelvinVector::KelvinVectorType; + +template +using KelvinMatrix = MathLib::KelvinVector::KelvinMatrixType; + +template +using GlobalDimVector = Eigen::Vector; + +/// Convenience alias for not a number. +static constexpr double nan = std::numeric_limits::quiet_NaN(); + +/// Used to set a Kelvin vector to all not-a-number. +template +constexpr KelvinVector KVnan() +{ + return KelvinVector::Constant(nan); +} + +/// Used to set a Kelvin matrix to all not-a-number. +template +constexpr KelvinMatrix KMnan() +{ + return KelvinMatrix::Constant(nan); +} + +/// Used to set a Kelvin vector to all zero. +template +constexpr KelvinVector KVzero() +{ + return KelvinVector::Zero(); +} + +/// Represents a previous state of type T. +template +struct PrevState +{ + PrevState() = default; + explicit PrevState(T const& t) : t{t} {} + explicit PrevState(T&& t) : t{std::move(t)} {} + + PrevState& operator=(T const& u) + { + t = u; + return *this; + } + + PrevState& operator=(T&& u) + { + t = std::move(u); + return *this; + } + + T& operator*() { return t; } + T const& operator*() const { return t; } + + T* operator->() { return &t; } + T const* operator->() const { return &t; } + +private: + T t; +}; + +struct SpaceTimeData +{ + ParameterLib::SpatialPosition x; + double t; + double dt; +}; + +struct MediaData +{ + explicit MediaData(MaterialPropertyLib::Medium const& medium) + : medium{medium}, solid{medium.phase("Solid")} + { + } + + MaterialPropertyLib::Medium const& medium; + MaterialPropertyLib::Phase const& solid; +}; + +template +struct TemperatureData +{ + double T; +}; + +template +struct StrainData +{ + // TODO Move initialization to the local assembler. + KelvinVector eps = KVnan(); + + static auto reflect() + { + using Self = StrainData; + + return ProcessLib::Reflection::reflectWithName("epsilon", &Self::eps); + } +}; + +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h new file mode 100644 index 00000000000..4c7c1745777 --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h @@ -0,0 +1,81 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 "Gravity.h" +#include "ProcessLib/Reflection/ReflectionData.h" +#include "SolidDensity.h" +#include "SolidMechanics.h" + +namespace ProcessLib::SmallDeformation +{ +namespace ConstitutiveRelations +{ +/// Data whose state must be tracked by the process. +template +struct StatefulData +{ + StressData stress_data; + + static auto reflect() + { + using Self = StatefulData; + + return Reflection::reflectWithoutName(&Self::stress_data); + } +}; + +/// Data whose state must be tracked by the process. +template +struct StatefulDataPrev +{ + PrevState> stress_data; + + StatefulDataPrev& operator=( + StatefulData const& state) + { + stress_data = state.stress_data; + + return *this; + } +}; + +/// Data that is needed for output purposes, but not directly for the assembly. +template +struct OutputData +{ + StrainData eps_data; + + static auto reflect() + { + using Self = OutputData; + + return Reflection::reflectWithoutName(&Self::eps_data); + } +}; + +/// Data that is needed for the equation system assembly. +template +struct ConstitutiveData +{ + SolidMechanicsDataStateless s_mech_data; + GravityData grav_data; +}; + +/// Data that stores intermediate values, which are not needed outside the +/// constitutive setting. +template +struct ConstitutiveTempData +{ + PrevState> eps_data_prev; + SolidDensityData rho_S_data; +}; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveModels.h new file mode 100644 index 00000000000..7daed425a5d --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveModels.h @@ -0,0 +1,38 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 "Gravity.h" +#include "SolidDensity.h" +#include "SolidMechanics.h" + +namespace ProcessLib::SmallDeformation +{ +namespace ConstitutiveRelations +{ +/// Constitutive models used for assembly. +template +struct ConstitutiveModels +{ + template + explicit ConstitutiveModels( + TRMProcessData const& process_data, + SolidConstitutiveRelation const& solid_material) + : s_mech_model(solid_material), + grav_model(process_data.specific_body_force) + { + } + + SolidMechanicsModel s_mech_model; + SolidDensityModel rho_S_model; + GravityModel grav_model; +}; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp new file mode 100644 index 00000000000..f691c5efbee --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp @@ -0,0 +1,62 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + */ + +#include "ConstitutiveSetting.h" + +#include "Invoke.h" + +namespace ProcessLib::SmallDeformation +{ +namespace ConstitutiveRelations +{ +template +void ConstitutiveSetting::eval( + ConstitutiveModels& models, double const t, + double const dt, ParameterLib::SpatialPosition const& x_position, + MaterialPropertyLib::Medium const& medium, double const T_ref, + KelvinVector const& eps, + KelvinVector const& eps_prev, + StatefulData& state, + StatefulDataPrev const& prev_state, + MaterialStateData& mat_state, + ConstitutiveTempData& tmp, + OutputData& out, + ConstitutiveData& cd) const +{ + namespace MPL = MaterialPropertyLib; + + auto& eps_data = out.eps_data; + eps_data.eps = eps; + auto& eps_data_prev = tmp.eps_data_prev; + eps_data_prev->eps = eps_prev; + auto& rho_S_data = tmp.rho_S_data; + + auto& s_mech_data = cd.s_mech_data; + auto& grav_data = cd.grav_data; + + TemperatureData const T_data{T_ref}; + SpaceTimeData const x_t{x_position, t, dt}; + MediaData const media_data{medium}; + + assertEvalArgsUnique(models.s_mech_model); + models.s_mech_model.eval(x_t, T_data, eps_data, eps_data_prev, mat_state, + prev_state.stress_data, state.stress_data, + s_mech_data); + + assertEvalArgsUnique(models.rho_S_model); + models.rho_S_model.eval(x_t, media_data, T_data, rho_S_data); + + assertEvalArgsUnique(models.grav_model); + models.grav_model.eval(rho_S_data, grav_data); +} + +template struct ConstitutiveSetting<2>; +template struct ConstitutiveSetting<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.h new file mode 100644 index 00000000000..d9b275c0b46 --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.h @@ -0,0 +1,39 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 "ConstitutiveData.h" +#include "ConstitutiveModels.h" + +namespace ProcessLib::SmallDeformation +{ +namespace ConstitutiveRelations +{ +template +struct ConstitutiveSetting +{ + /// Evaluate the constitutive setting. + void eval(ConstitutiveModels& models, double const t, + double const dt, ParameterLib::SpatialPosition const& x_position, + MaterialPropertyLib::Medium const& medium, double const T_ref, + KelvinVector const& eps, + KelvinVector const& eps_prev, + StatefulData& state, + StatefulDataPrev const& prev_state, + MaterialStateData& mat_state, + ConstitutiveTempData& tmp, + OutputData& out, + ConstitutiveData& cd) const; +}; + +extern template struct ConstitutiveSetting<2>; +extern template struct ConstitutiveSetting<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/CreateConstitutiveSetting.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/CreateConstitutiveSetting.cpp new file mode 100644 index 00000000000..85162b2c2d3 --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/CreateConstitutiveSetting.cpp @@ -0,0 +1,33 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + */ + +#include "CreateConstitutiveSetting.h" + +#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h" + +namespace ProcessLib::SmallDeformation +{ +namespace ConstitutiveRelations +{ +template +std::map>> +CreateConstitutiveSetting::createSolidConstitutiveRelations( + std::vector> const& parameters, + std::optional const& + local_coordinate_system, + BaseLib::ConfigTree const& config) +{ + return MaterialLib::Solids::createConstitutiveRelations( + parameters, local_coordinate_system, config); +} + +template struct CreateConstitutiveSetting<2>; +template struct CreateConstitutiveSetting<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/CreateConstitutiveSetting.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/CreateConstitutiveSetting.h new file mode 100644 index 00000000000..76bedc5d4bf --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/CreateConstitutiveSetting.h @@ -0,0 +1,46 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 +#include +#include + +#include "SolidMechanics.h" + +namespace BaseLib +{ +class ConfigTree; +} +namespace ParameterLib +{ +struct ParameterBase; +struct CoordinateSystem; +} // namespace ParameterLib + +namespace ProcessLib::SmallDeformation +{ +namespace ConstitutiveRelations +{ +template +struct CreateConstitutiveSetting +{ + static std::map>> + createSolidConstitutiveRelations( + std::vector> const& + parameters, + std::optional const& + local_coordinate_system, + BaseLib::ConfigTree const& config); +}; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.cpp new file mode 100644 index 00000000000..2b25e1b3ac7 --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.cpp @@ -0,0 +1,26 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + */ + +#include "Gravity.h" + +namespace ProcessLib::SmallDeformation +{ +template +void GravityModel::eval( + SolidDensityData const& rho_S_data, GravityData& out) const +{ + auto const rho_SR = rho_S_data.rho_SR; + auto const b = specific_body_force_; + + out.volumetric_body_force = rho_SR * b; +} + +template struct GravityModel<2>; +template struct GravityModel<3>; +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.h new file mode 100644 index 00000000000..6fa2e365c1a --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.h @@ -0,0 +1,41 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 "SolidDensity.h" + +namespace ProcessLib::SmallDeformation +{ +template +struct GravityData +{ + GlobalDimVector volumetric_body_force; +}; + +template +struct GravityModel +{ + explicit GravityModel( + Eigen::Vector const& specific_body_force) + : specific_body_force_(specific_body_force) + { + } + + void eval(SolidDensityData const& rho_S_data, + GravityData& out) const; + +private: + // TODO (naumov) Do we need to store this for each integration point? + Eigen::Vector const specific_body_force_; +}; + +extern template struct GravityModel<2>; +extern template struct GravityModel<3>; +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/Invoke.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/Invoke.h new file mode 100644 index 00000000000..94309e15e5d --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/Invoke.h @@ -0,0 +1,54 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 + +namespace ProcessLib::SmallDeformation +{ +namespace detail +{ +template +constexpr bool areEvalArgumentTypesUnique(Result (Class::*)(Args...)) +{ + using namespace boost::mp11; + return mp_is_set< + mp_transform>>::value; +} + +template +constexpr bool areEvalArgumentTypesUnique(Result (Class::*)(Args...) const) +{ + using namespace boost::mp11; + return mp_is_set< + mp_transform>>::value; +} +} // namespace detail + +/// Checks whether the argument types of the eval() method of the given type T +/// are unique. +/// +/// Argument types differing only in constness, reference or volatility are +/// considered equal. +template +constexpr bool areEvalArgumentTypesUnique() +{ + return detail::areEvalArgumentTypesUnique(&T::eval); +} + +/// Statically asserts that the argument types of the passed Model's eval() +/// method are unique. +template +constexpr void assertEvalArgsUnique(Model const&) +{ + static_assert(areEvalArgumentTypesUnique>()); +} + +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.cpp new file mode 100644 index 00000000000..fd77ffb2b9c --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.cpp @@ -0,0 +1,31 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + */ + +#include "SolidDensity.h" + +namespace ProcessLib::SmallDeformation +{ +template +void SolidDensityModel::eval( + SpaceTimeData const& x_t, + MediaData const& media_data, + TemperatureData const& T_data, + SolidDensityData& out) const +{ + namespace MPL = MaterialPropertyLib; + MPL::VariableArray variables; + variables.temperature = T_data.T; + + out.rho_SR = media_data.solid.property(MPL::PropertyType::density) + .template value(variables, x_t.x, x_t.t, x_t.dt); +} + +template struct SolidDensityModel<2>; +template struct SolidDensityModel<3>; +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.h new file mode 100644 index 00000000000..95bea57e630 --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.h @@ -0,0 +1,31 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 "Base.h" + +namespace ProcessLib::SmallDeformation +{ +struct SolidDensityData +{ + double rho_SR; +}; + +template +struct SolidDensityModel +{ + void eval(SpaceTimeData const& x_t, MediaData const& media_data, + TemperatureData const& T_data, + SolidDensityData& out) const; +}; + +extern template struct SolidDensityModel<2>; +extern template struct SolidDensityModel<3>; +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp new file mode 100644 index 00000000000..a7b889eaf4b --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp @@ -0,0 +1,67 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + */ + +#include "SolidMechanics.h" + +namespace ProcessLib::SmallDeformation +{ +namespace ConstitutiveRelations +{ +template +void SolidMechanicsModel::eval( + SpaceTimeData const& x_t, + TemperatureData const& T_data, + StrainData const& eps_data, + PrevState> const& eps_data_prev, + MaterialStateData& mat_state, + PrevState> const& stress_data_prev, + StressData& stress_data, + SolidMechanicsDataStateless& current_stateless) const +{ + namespace MPL = MaterialPropertyLib; + + // current state + MPL::VariableArray variables; + { + // thermodynamic forces + variables.stress = stress_data.sigma; + variables.mechanical_strain = eps_data.eps; + + // external state variables + variables.temperature = T_data.T; + } + + // previous state + MPL::VariableArray variables_prev; + { + // thermodynamic forces + variables_prev.stress = stress_data_prev->sigma; + variables_prev.mechanical_strain = eps_data_prev->eps; + + // external state variables + variables_prev.temperature = T_data.T; + } + + auto solution = solid_material_.integrateStress( + variables_prev, variables, x_t.t, x_t.x, x_t.dt, + *mat_state.material_state_variables); + + if (!solution) + { + OGS_FATAL("Computation of local constitutive relation failed."); + } + + std::tie(stress_data.sigma, mat_state.material_state_variables, + current_stateless.stiffness_tensor) = std::move(*solution); +} + +template struct SolidMechanicsModel<2>; +template struct SolidMechanicsModel<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h new file mode 100644 index 00000000000..1fc9fc4d2f5 --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h @@ -0,0 +1,61 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 "MaterialLib/SolidModels/MechanicsBase.h" +#include "MaterialState.h" +#include "StressData.h" + +namespace ProcessLib::SmallDeformation +{ +namespace ConstitutiveRelations +{ +template +struct SolidMechanicsDataStateless +{ + KelvinMatrix stiffness_tensor = KMnan(); +}; + +template +using SolidConstitutiveRelation = + MaterialLib::Solids::MechanicsBase; + +template +struct SolidMechanicsModel +{ + explicit SolidMechanicsModel( + SolidConstitutiveRelation const& solid_material) + : solid_material_(solid_material) + { + } + + void eval( + SpaceTimeData const& x_t, + TemperatureData const& T_data, + StrainData const& eps_data, + PrevState> const& eps_data_prev, + MaterialStateData& mat_state, + PrevState> const& stress_data_prev, + StressData& stress_data, + SolidMechanicsDataStateless& current_stateless) const; + + auto getInternalVariables() const + { + return solid_material_.getInternalVariables(); + } + +private: + SolidConstitutiveRelation const& solid_material_; +}; + +extern template struct SolidMechanicsModel<2>; +extern template struct SolidMechanicsModel<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/StressData.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/StressData.h new file mode 100644 index 00000000000..da7039d601e --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/StressData.h @@ -0,0 +1,32 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 "Base.h" + +namespace ProcessLib::SmallDeformation +{ + +template +struct StressData +{ + // Stress is stateful for some constitutive settings and therefore must be + // initialized to something valid, e.g., zero. + // TODO find a better solution for that. + KelvinVector sigma = KVzero(); + + static auto reflect() + { + using Self = StressData; + + return ProcessLib::Reflection::reflectWithName("sigma", &Self::sigma); + } +}; +} // namespace ProcessLib::SmallDeformation From cd49e2e8eaebea1cbf585357a4854ea8a60459f0 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Tue, 10 Oct 2023 14:46:41 +0200 Subject: [PATCH 04/10] [PL/SD] Use constitutive relations eval in FEM - Removing sigma and eps from IPData, now stored in current state and prev state. - Remove associated ip getter and setter, now handled through reflections - Stress, stiffness tensor, and volumetric body force are now fetched from newly updated constitutive data (CD). --- ProcessLib/Deformation/MaterialForces.h | 7 +- .../LocalAssemblerInterface.h | 37 ++-- .../SmallDeformation/SmallDeformationFEM.h | 205 +++++++----------- .../SmallDeformationProcess.cpp | 32 ++- 4 files changed, 114 insertions(+), 167 deletions(-) diff --git a/ProcessLib/Deformation/MaterialForces.h b/ProcessLib/Deformation/MaterialForces.h index 0ae7b363836..3852be91f46 100644 --- a/ProcessLib/Deformation/MaterialForces.h +++ b/ProcessLib/Deformation/MaterialForces.h @@ -36,12 +36,13 @@ struct MaterialForcesInterface template std::vector const& getMaterialForces( std::vector const& local_x, std::vector& nodal_values, IntegrationMethod const& _integration_method, IPData const& _ip_data, - MeshLib::Element const& element, bool const is_axially_symmetric) + StressData const& stress_data, MeshLib::Element const& element, + bool const is_axially_symmetric) { unsigned const n_integration_points = _integration_method.getNumberOfPoints(); @@ -52,7 +53,7 @@ std::vector const& getMaterialForces( for (unsigned ip = 0; ip < n_integration_points; ip++) { - auto const& sigma = _ip_data[ip].sigma; + auto const& sigma = stress_data[ip].stress_data.sigma; auto const& N = _ip_data[ip].N; auto const& dNdx = _ip_data[ip].dNdx; diff --git a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h index 9954f7afc01..237dcc590cd 100644 --- a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h +++ b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h @@ -12,6 +12,8 @@ #include +#include "ConstitutiveRelations/ConstitutiveData.h" +#include "ConstitutiveRelations/ConstitutiveSetting.h" #include "ConstitutiveRelations/MaterialState.h" #include "MaterialLib/SolidModels/MechanicsBase.h" #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h" @@ -54,6 +56,10 @@ struct SmallDeformationLocalAssemblerInterface material_states_.emplace_back( solid_material_.createMaterialStateVariables()); } + + current_states_.resize(n_integration_points); + prev_states_.resize(n_integration_points); + output_data_.resize(n_integration_points); } virtual std::size_t setIPDataInitialConditions( std::string const& name, double const* values, @@ -65,19 +71,6 @@ struct SmallDeformationLocalAssemblerInterface std::vector const& dof_table, std::vector& cache) const = 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 const& getIntPtEpsilon( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - // TODO move to NumLib::ExtrapolatableElement unsigned getNumberOfIntegrationPoints() const { @@ -110,17 +103,35 @@ struct SmallDeformationLocalAssemblerInterface return *material_states_[integration_point].material_state_variables; } + static auto getReflectionDataForOutput() + { + using Self = SmallDeformationLocalAssemblerInterface; + + return ProcessLib::Reflection::reflectWithoutName( + &Self::current_states_, &Self::output_data_); + } + protected: SmallDeformationProcessData& process_data_; // Material state is special, because it contains both the current and the // old state. std::vector> material_states_; + std::vector> + current_states_; // TODO maybe do not store but rather re-evaluate for + // state update + std::vector< + typename ConstitutiveRelations::StatefulDataPrev> + prev_states_; + std::vector> + output_data_; NumLib::GenericIntegrationMethod const& integration_method_; MeshLib::Element const& element_; bool const is_axially_symmetric_; MaterialLib::Solids::MechanicsBase const& solid_material_; + typename ConstitutiveRelations::ConstitutiveSetting + constitutive_setting; }; } // namespace SmallDeformation diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h index c01136a7b84..9c99e5ba552 100644 --- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h +++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h @@ -30,6 +30,7 @@ #include "ProcessLib/Deformation/LinearBMatrix.h" #include "ProcessLib/LocalAssemblerInterface.h" #include "ProcessLib/LocalAssemblerTraits.h" +#include "ProcessLib/Reflection/ReflectionSetIPData.h" namespace ProcessLib { @@ -41,16 +42,12 @@ template struct IntegrationPointData final { - typename BMatricesType::KelvinVectorType sigma, sigma_prev; - typename BMatricesType::KelvinVectorType eps; double free_energy_density = 0; double integration_weight; typename ShapeMatricesType::NodalRowVectorType N; typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx; - void pushBackState() { sigma_prev = sigma; } - EIGEN_MAKE_ALIGNED_OPERATOR_NEW; }; @@ -124,16 +121,6 @@ class SmallDeformationLocalAssembler ip_data.N = sm.N; ip_data.dNdx = sm.dNdx; - static const int kelvin_vector_size = - MathLib::KelvinVector::kelvin_vector_dimensions( - DisplacementDim); - // Initialize current time step values - ip_data.sigma.setZero(kelvin_vector_size); - ip_data.eps.setZero(kelvin_vector_size); - - // Previous time step values are not initialized and are set later. - ip_data.sigma_prev.resize(kelvin_vector_size); - _secondary_data.N[ip] = shape_matrices[ip].N; } } @@ -153,23 +140,9 @@ class SmallDeformationLocalAssembler this->element_.getID()); } - if (name == "sigma_ip") + if (name.starts_with("material_state_variable_")) { - 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 setSigma(values); - } - if (name.starts_with("material_state_variable_") && - name.ends_with("_ip")) - { - std::string const variable_name = - name.substr(24, name.size() - 24 - 3); + std::string const variable_name = name.substr(24, name.size() - 24); DBUG("Setting material state variable '{:s}'", variable_name); auto const& internal_variables = @@ -192,9 +165,14 @@ class SmallDeformationLocalAssembler "Could not find variable {:s} in solid material model's " "internal variables.", variable_name); + return 0; } - 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, this->current_states_); } void initializeConcrete() override @@ -203,7 +181,7 @@ class SmallDeformationLocalAssembler this->integration_method_.getNumberOfPoints(); for (unsigned ip = 0; ip < n_integration_points; ip++) { - auto& ip_data = _ip_data[ip]; + auto const& ip_data = _ip_data[ip]; ParameterLib::SpatialPosition const x_position{ std::nullopt, this->element_.getID(), ip, @@ -215,7 +193,7 @@ class SmallDeformationLocalAssembler /// Set initial stress from parameter. if (this->process_data_.initial_stress != nullptr) { - ip_data.sigma = + this->current_states_[ip].stress_data.sigma.noalias() = MathLib::KelvinVector::symmetricTensorToKelvinVector< DisplacementDim>((*this->process_data_.initial_stress)( std::numeric_limits< @@ -228,14 +206,16 @@ class SmallDeformationLocalAssembler this->solid_material_.initializeInternalStateVariables( t, x_position, *material_state.material_state_variables); - ip_data.pushBackState(); material_state.pushBackState(); } + + for (unsigned ip = 0; ip < n_integration_points; ip++) + { + this->prev_states_[ip] = this->current_states_[ip]; + } } - /// Updates sigma, eps, and state through passed integration point data. - /// \returns tangent stiffness and density. - std::tuple, double> + typename ConstitutiveRelations::ConstitutiveData updateConstitutiveRelations( Eigen::Ref const& u, Eigen::Ref const& u_prev, @@ -243,20 +223,17 @@ class SmallDeformationLocalAssembler double const dt, IntegrationPointData& ip_data, - MaterialStateData& material_state) const + typename ConstitutiveRelations::ConstitutiveSetting& + CS, + MaterialPropertyLib::Medium const& medium, + typename ConstitutiveRelations::StatefulData& + current_state, + typename ConstitutiveRelations::StatefulDataPrev const& + prev_state, + MaterialStateData& material_state, + typename ConstitutiveRelations::OutputData& + output_data) const { - auto const& solid_phase = - this->process_data_.media_map.getMedium(this->element_.getID()) - ->phase("Solid"); - - MPL::VariableArray variables_prev; - MPL::VariableArray variables; - - auto const& sigma_prev = ip_data.sigma_prev; - - auto& eps = ip_data.eps; - auto& sigma = ip_data.sigma; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; auto const x_coord = @@ -268,44 +245,24 @@ class SmallDeformationLocalAssembler typename BMatricesType::BMatrixType>( dNdx, N, x_coord, this->is_axially_symmetric_); - eps.noalias() = B * u; - - variables_prev.stress - .emplace>( - sigma_prev); - variables_prev.mechanical_strain - .emplace>( - B * u_prev); - double const T_ref = this->process_data_.reference_temperature ? (*this->process_data_.reference_temperature)(t, x_position)[0] : std::numeric_limits::quiet_NaN(); - variables_prev.temperature = T_ref; - variables.mechanical_strain - .emplace>( - eps); - variables.temperature = T_ref; - - auto const rho = - solid_phase[MPL::PropertyType::density].template value( - variables, x_position, t, dt); + typename ConstitutiveRelations::ConstitutiveModels + models{this->process_data_, this->solid_material_}; + typename ConstitutiveRelations::ConstitutiveTempData + tmp; + typename ConstitutiveRelations::ConstitutiveData CD; - auto&& solution = this->solid_material_.integrateStress( - variables_prev, variables, t, x_position, dt, - *material_state.material_state_variables); + CS.eval(models, t, dt, x_position, // + medium, // + T_ref, B * u, B * u_prev, // + current_state, prev_state, material_state, tmp, output_data, + CD); - if (!solution) - { - OGS_FATAL("Computation of local constitutive relation failed."); - } - - MathLib::KelvinVector::KelvinMatrixType C; - std::tie(sigma, material_state.material_state_variables, C) = - std::move(*solution); - - return {C, rho}; + return CD; } void assemble(double const /*t*/, double const /*dt*/, @@ -345,7 +302,10 @@ class SmallDeformationLocalAssembler ParameterLib::SpatialPosition x_position; x_position.setElementID(this->element_.getID()); - auto const& b = this->process_data_.specific_body_force; + typename ConstitutiveRelations::ConstitutiveSetting + constitutive_setting; + auto const& medium = + *this->process_data_.media_map.getMedium(this->element_.getID()); for (unsigned ip = 0; ip < n_integration_points; ip++) { @@ -363,14 +323,18 @@ class SmallDeformationLocalAssembler typename BMatricesType::BMatrixType>( dNdx, N, x_coord, this->is_axially_symmetric_); - auto const& sigma = _ip_data[ip].sigma; - - auto const [C, rho] = updateConstitutiveRelations( + auto const CD = updateConstitutiveRelations( u, u_prev, x_position, t, dt, _ip_data[ip], - this->material_states_[ip]); + constitutive_setting, medium, this->current_states_[ip], + this->prev_states_[ip], this->material_states_[ip], + this->output_data_[ip]); + + auto const& sigma = this->current_states_[ip].stress_data.sigma; + auto const& b = CD.grav_data.volumetric_body_force; + auto const& C = CD.s_mech_data.stiffness_tensor; local_b.noalias() -= - (B.transpose() * sigma - N_u_op(N).transpose() * rho * b) * w; + (B.transpose() * sigma - N_u_op(N).transpose() * b) * w; local_Jac.noalias() += B.transpose() * C * B * w; } } @@ -387,16 +351,24 @@ class SmallDeformationLocalAssembler ParameterLib::SpatialPosition x_position; x_position.setElementID(this->element_.getID()); + typename ConstitutiveRelations::ConstitutiveSetting + constitutive_setting; + + auto const& medium = + *this->process_data_.media_map.getMedium(this->element_.getID()); + for (unsigned ip = 0; ip < n_integration_points; ip++) { x_position.setIntegrationPoint(ip); - updateConstitutiveRelations(local_x, local_x_prev, x_position, t, - dt, _ip_data[ip], - this->material_states_[ip]); + updateConstitutiveRelations( + local_x, local_x_prev, x_position, t, dt, _ip_data[ip], + constitutive_setting, medium, this->current_states_[ip], + this->prev_states_[ip], this->material_states_[ip], + this->output_data_[ip]); - auto& eps = _ip_data[ip].eps; - auto& sigma = _ip_data[ip].sigma; + auto const& eps = this->output_data_[ip].eps_data.eps; + auto const& sigma = this->current_states_[ip].stress_data.sigma; auto& state = *this->material_states_[ip].material_state_variables; // Update free energy density needed for material forces. @@ -404,9 +376,13 @@ class SmallDeformationLocalAssembler this->solid_material_.computeFreeEnergyDensity( t, x_position, dt, eps, sigma, state); - _ip_data[ip].pushBackState(); this->material_states_[ip].pushBackState(); } + + for (unsigned ip = 0; ip < n_integration_points; ip++) + { + this->prev_states_[ip] = this->current_states_[ip]; + } } std::vector const& getMaterialForces( @@ -417,9 +393,9 @@ class SmallDeformationLocalAssembler DisplacementDim, ShapeFunction, ShapeMatricesType, typename BMatricesType::NodalForceVectorType, NodalDisplacementVectorType, GradientVectorType, - GradientMatrixType>(local_x, nodal_values, - this->integration_method_, _ip_data, - this->element_, this->is_axially_symmetric_); + GradientMatrixType>( + local_x, nodal_values, this->integration_method_, _ip_data, + this->current_states_, this->element_, this->is_axially_symmetric_); } Eigen::Map getShapeMatrix( @@ -447,41 +423,6 @@ class SmallDeformationLocalAssembler return cache; } - std::size_t setSigma(double const* values) - { - return ProcessLib::setIntegrationPointKelvinVectorData( - values, _ip_data, &IpData::sigma); - } - - // TODO (naumov) This method is same as getIntPtSigma but for arguments and - // the ordering of the cache_mat. - // There should be only one. - std::vector getSigma() const override - { - return ProcessLib::getIntegrationPointKelvinVectorData( - _ip_data, &IpData::sigma); - } - - std::vector const& getIntPtSigma( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const override - { - return ProcessLib::getIntegrationPointKelvinVectorData( - _ip_data, &IpData::sigma, cache); - } - - std::vector const& getIntPtEpsilon( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const override - { - return ProcessLib::getIntegrationPointKelvinVectorData( - _ip_data, &IpData::eps, cache); - } - void computeSecondaryVariableConcrete( double const /*t*/, double const /*dt*/, Eigen::VectorXd const& /*x*/, Eigen::VectorXd const& /*x_prev*/) override @@ -498,7 +439,7 @@ class SmallDeformationLocalAssembler for (unsigned ip = 0; ip < n_integration_points; ip++) { x_position.setIntegrationPoint(ip); - auto const& sigma = _ip_data[ip].sigma; + auto const& sigma = this->current_states_[ip].stress_data.sigma; sigma_sum += sigma; } diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp index fef518075af..e9ae759b7ff 100644 --- a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp +++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp @@ -16,6 +16,8 @@ #include "MeshLib/Utils/getOrCreateMeshProperty.h" #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h" #include "ProcessLib/Process.h" +#include "ProcessLib/Reflection/ReflectionForExtrapolation.h" +#include "ProcessLib/Reflection/ReflectionForIPWriters.h" #include "ProcessLib/SmallDeformation/CreateLocalAssemblers.h" #include "ProcessLib/Utils/SetIPDataInitialConditions.h" #include "SmallDeformationFEM.h" @@ -66,15 +68,11 @@ SmallDeformationProcess::SmallDeformationProcess( const_cast(mesh), "principal_stress_values", MeshLib::MeshItemType::Cell, 3); - // 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, - &LocalAssemblerInterface::getSigma)); + ProcessLib::Reflection::addReflectedIntegrationPointWriters< + DisplacementDim>(SmallDeformationLocalAssemblerInterface< + DisplacementDim>::getReflectionDataForOutput(), + _integration_point_writer, integration_order, + _local_assemblers); } template @@ -110,15 +108,10 @@ void SmallDeformationProcess::initializeConcreteProcess( 1, &LocalAssemblerInterface::getIntPtFreeEnergyDensity); - add_secondary_variable("sigma", - MathLib::KelvinVector::KelvinVectorType< - DisplacementDim>::RowsAtCompileTime, - &LocalAssemblerInterface::getIntPtSigma); - - add_secondary_variable("epsilon", - MathLib::KelvinVector::KelvinVectorType< - DisplacementDim>::RowsAtCompileTime, - &LocalAssemblerInterface::getIntPtEpsilon); + ProcessLib::Reflection::addReflectedSecondaryVariables( + SmallDeformationLocalAssemblerInterface< + DisplacementDim>::getReflectionDataForOutput(), + _secondary_variables, getExtrapolator(), _local_assemblers); // // enable output of internal variables defined by material models @@ -132,8 +125,9 @@ void SmallDeformationProcess::initializeConcreteProcess( _process_data.solid_materials, _local_assemblers, _integration_point_writer, integration_order); + bool const remove_name_suffix = true; setIPDataInitialConditions(_integration_point_writer, mesh.getProperties(), - _local_assemblers); + _local_assemblers, remove_name_suffix); // Initialize local assemblers after all variables have been set. GlobalExecutor::executeMemberOnDereferenced( From 08dce7a21eea771de7737a09deb41ad04e0eac48 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Tue, 10 Oct 2023 14:58:04 +0200 Subject: [PATCH 05/10] [PL/SD] Move shape function independent fct in I/F --- .../LocalAssemblerInterface.h | 95 ++++++++++++++++++- .../SmallDeformation/SmallDeformationFEM.h | 93 ------------------ 2 files changed, 92 insertions(+), 96 deletions(-) diff --git a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h index 237dcc590cd..5ac074a06ee 100644 --- a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h +++ b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h @@ -10,6 +10,7 @@ #pragma once +#include #include #include "ConstitutiveRelations/ConstitutiveData.h" @@ -21,6 +22,7 @@ #include "NumLib/Fem/Integration/GenericIntegrationMethod.h" #include "ProcessLib/Deformation/MaterialForces.h" #include "ProcessLib/LocalAssemblerInterface.h" +#include "ProcessLib/Reflection/ReflectionSetIPData.h" #include "ProcessLib/Utils/SetOrGetIntegrationPointData.h" #include "SmallDeformationProcessData.h" @@ -61,9 +63,96 @@ struct SmallDeformationLocalAssemblerInterface prev_states_.resize(n_integration_points); output_data_.resize(n_integration_points); } - virtual std::size_t setIPDataInitialConditions( - std::string const& name, double const* values, - int const integration_order) = 0; + /// Returns number of read integration points. + std::size_t setIPDataInitialConditions(std::string const& 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.starts_with("material_state_variable_")) + { + std::string const variable_name = name.substr(24, name.size() - 24); + DBUG("Setting material state variable '{:s}'", variable_name); + + auto const& internal_variables = + solid_material_.getInternalVariables(); + if (auto const iv = std::find_if( + begin(internal_variables), end(internal_variables), + [&variable_name](auto const& iv) + { return iv.name == variable_name; }); + iv != end(internal_variables)) + { + return ProcessLib:: + setIntegrationPointDataMaterialStateVariables( + values, material_states_, + &MaterialStateData< + DisplacementDim>::material_state_variables, + iv->reference); + } + + WARN( + "Could not find variable {:s} in solid material model's " + "internal variables.", + variable_name); + 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_); + } + + void computeSecondaryVariableConcrete( + double const /*t*/, double const /*dt*/, Eigen::VectorXd const& /*x*/, + Eigen::VectorXd const& /*x_prev*/) override + { + int const elem_id = element_.getID(); + ParameterLib::SpatialPosition x_position; + x_position.setElementID(elem_id); + unsigned const n_integration_points = + integration_method_.getNumberOfPoints(); + + auto sigma_sum = MathLib::KelvinVector::tensorToKelvin( + Eigen::Matrix::Zero()); + + for (unsigned ip = 0; ip < n_integration_points; ip++) + { + x_position.setIntegrationPoint(ip); + auto const& sigma = current_states_[ip].stress_data.sigma; + sigma_sum += sigma; + } + + Eigen::Matrix const sigma_avg = + MathLib::KelvinVector::kelvinVectorToTensor(sigma_sum) / + n_integration_points; + + Eigen::SelfAdjointEigenSolver> e_s( + sigma_avg); + + Eigen::Map( + &(*process_data_.principal_stress_values)[elem_id * 3], 3) = + e_s.eigenvalues(); + + auto eigen_vectors = e_s.eigenvectors(); + + for (auto i = 0; i < 3; i++) + { + Eigen::Map( + &(*process_data_.principal_stress_vector[i])[elem_id * 3], 3) = + eigen_vectors.col(i); + } + } virtual std::vector const& getIntPtFreeEnergyDensity( const double t, diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h index 9c99e5ba552..7e72690a8bf 100644 --- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h +++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h @@ -10,7 +10,6 @@ #pragma once -#include #include #include #include @@ -30,7 +29,6 @@ #include "ProcessLib/Deformation/LinearBMatrix.h" #include "ProcessLib/LocalAssemblerInterface.h" #include "ProcessLib/LocalAssemblerTraits.h" -#include "ProcessLib/Reflection/ReflectionSetIPData.h" namespace ProcessLib { @@ -125,56 +123,6 @@ class SmallDeformationLocalAssembler } } - /// Returns number of read integration points. - std::size_t setIPDataInitialConditions(std::string const& name, - double const* values, - int const integration_order) override - { - 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.starts_with("material_state_variable_")) - { - std::string const variable_name = name.substr(24, name.size() - 24); - DBUG("Setting material state variable '{:s}'", variable_name); - - auto const& internal_variables = - this->solid_material_.getInternalVariables(); - if (auto const iv = std::find_if( - begin(internal_variables), end(internal_variables), - [&variable_name](auto const& iv) - { return iv.name == variable_name; }); - iv != end(internal_variables)) - { - return ProcessLib:: - setIntegrationPointDataMaterialStateVariables( - values, this->material_states_, - &MaterialStateData< - DisplacementDim>::material_state_variables, - iv->reference); - } - - WARN( - "Could not find variable {:s} in solid material model's " - "internal variables.", - variable_name); - 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, this->current_states_); - } - void initializeConcrete() override { unsigned const n_integration_points = @@ -423,47 +371,6 @@ class SmallDeformationLocalAssembler return cache; } - void computeSecondaryVariableConcrete( - double const /*t*/, double const /*dt*/, Eigen::VectorXd const& /*x*/, - Eigen::VectorXd const& /*x_prev*/) override - { - int const elem_id = this->element_.getID(); - ParameterLib::SpatialPosition x_position; - x_position.setElementID(elem_id); - unsigned const n_integration_points = - this->integration_method_.getNumberOfPoints(); - - auto sigma_sum = MathLib::KelvinVector::tensorToKelvin( - Eigen::Matrix::Zero()); - - for (unsigned ip = 0; ip < n_integration_points; ip++) - { - x_position.setIntegrationPoint(ip); - auto const& sigma = this->current_states_[ip].stress_data.sigma; - sigma_sum += sigma; - } - - Eigen::Matrix const sigma_avg = - MathLib::KelvinVector::kelvinVectorToTensor(sigma_sum) / - n_integration_points; - - Eigen::SelfAdjointEigenSolver> e_s( - sigma_avg); - - Eigen::Map( - &(*this->process_data_.principal_stress_values)[elem_id * 3], 3) = - e_s.eigenvalues(); - - auto eigen_vectors = e_s.eigenvectors(); - - for (auto i = 0; i < 3; i++) - { - Eigen::Map( - &(*this->process_data_.principal_stress_vector[i])[elem_id * 3], - 3) = eigen_vectors.col(i); - } - } - private: static constexpr auto localDOF(std::vector const& x) { From 42ed73c8e66ba22b026f27a92a3b17f5c38c4ce4 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Wed, 11 Oct 2023 15:37:02 +0200 Subject: [PATCH 06/10] [PL/SD] Use StrongType for solid density --- .../ConstitutiveRelations/ConstitutiveData.h | 2 +- .../ConstitutiveRelations/ConstitutiveModels.h | 2 +- .../ConstitutiveSetting.cpp | 6 +++--- .../ConstitutiveRelations/Gravity.cpp | 5 ++--- .../ConstitutiveRelations/Gravity.h | 2 +- .../ConstitutiveRelations/SolidDensity.cpp | 17 ++++++----------- .../ConstitutiveRelations/SolidDensity.h | 12 +++--------- 7 files changed, 17 insertions(+), 29 deletions(-) diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h index 4c7c1745777..4823c300e0c 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h @@ -75,7 +75,7 @@ template struct ConstitutiveTempData { PrevState> eps_data_prev; - SolidDensityData rho_S_data; + SolidDensity rho_SR; }; } // namespace ConstitutiveRelations } // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveModels.h index 7daed425a5d..1ae91962f1f 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveModels.h +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveModels.h @@ -31,7 +31,7 @@ struct ConstitutiveModels } SolidMechanicsModel s_mech_model; - SolidDensityModel rho_S_model; + SolidDensityModel rho_S_model; GravityModel grav_model; }; } // namespace ConstitutiveRelations diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp index f691c5efbee..b3067c70af9 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp @@ -35,7 +35,7 @@ void ConstitutiveSetting::eval( eps_data.eps = eps; auto& eps_data_prev = tmp.eps_data_prev; eps_data_prev->eps = eps_prev; - auto& rho_S_data = tmp.rho_S_data; + auto& rho_SR = tmp.rho_SR; auto& s_mech_data = cd.s_mech_data; auto& grav_data = cd.grav_data; @@ -50,10 +50,10 @@ void ConstitutiveSetting::eval( s_mech_data); assertEvalArgsUnique(models.rho_S_model); - models.rho_S_model.eval(x_t, media_data, T_data, rho_S_data); + models.rho_S_model.eval(x_t, media_data, T_data, rho_SR); assertEvalArgsUnique(models.grav_model); - models.grav_model.eval(rho_S_data, grav_data); + models.grav_model.eval(rho_SR, grav_data); } template struct ConstitutiveSetting<2>; diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.cpp index 2b25e1b3ac7..d5d6fbe6d38 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.cpp +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.cpp @@ -13,12 +13,11 @@ namespace ProcessLib::SmallDeformation { template void GravityModel::eval( - SolidDensityData const& rho_S_data, GravityData& out) const + SolidDensity const& rho_SR, GravityData& out) const { - auto const rho_SR = rho_S_data.rho_SR; auto const b = specific_body_force_; - out.volumetric_body_force = rho_SR * b; + out.volumetric_body_force = *rho_SR * b; } template struct GravityModel<2>; diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.h index 6fa2e365c1a..8ab8883fd3d 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.h +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.h @@ -28,7 +28,7 @@ struct GravityModel { } - void eval(SolidDensityData const& rho_S_data, + void eval(SolidDensity const& rho_SR, GravityData& out) const; private: diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.cpp index fd77ffb2b9c..aa39cda70d8 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.cpp +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.cpp @@ -11,21 +11,16 @@ namespace ProcessLib::SmallDeformation { -template -void SolidDensityModel::eval( - SpaceTimeData const& x_t, - MediaData const& media_data, - TemperatureData const& T_data, - SolidDensityData& out) const +void SolidDensityModel::eval(SpaceTimeData const& x_t, + MediaData const& media_data, + TemperatureData const& T_data, + SolidDensity& out) const { namespace MPL = MaterialPropertyLib; MPL::VariableArray variables; variables.temperature = T_data.T; - out.rho_SR = media_data.solid.property(MPL::PropertyType::density) - .template value(variables, x_t.x, x_t.t, x_t.dt); + *out = media_data.solid.property(MPL::PropertyType::density) + .template value(variables, x_t.x, x_t.t, x_t.dt); } - -template struct SolidDensityModel<2>; -template struct SolidDensityModel<3>; } // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.h index 95bea57e630..6bfcd8c48c3 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.h +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.h @@ -10,22 +10,16 @@ #pragma once #include "Base.h" +#include "BaseLib/StrongType.h" namespace ProcessLib::SmallDeformation { -struct SolidDensityData -{ - double rho_SR; -}; +using SolidDensity = BaseLib::StrongType; -template struct SolidDensityModel { void eval(SpaceTimeData const& x_t, MediaData const& media_data, TemperatureData const& T_data, - SolidDensityData& out) const; + SolidDensity& out) const; }; - -extern template struct SolidDensityModel<2>; -extern template struct SolidDensityModel<3>; } // namespace ProcessLib::SmallDeformation From 5c468812e1f41886b4fdd1487937807e2571a782 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Wed, 11 Oct 2023 15:47:48 +0200 Subject: [PATCH 07/10] [PL/SD] Use StrongType for gravity --- .../ConstitutiveRelations/ConstitutiveData.h | 2 +- .../ConstitutiveRelations/ConstitutiveModels.h | 4 ++-- .../ConstitutiveRelations/ConstitutiveSetting.cpp | 6 +++--- .../SmallDeformation/ConstitutiveRelations/Gravity.cpp | 6 ++---- .../SmallDeformation/ConstitutiveRelations/Gravity.h | 9 ++++----- ProcessLib/SmallDeformation/SmallDeformationFEM.h | 2 +- 6 files changed, 13 insertions(+), 16 deletions(-) diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h index 4823c300e0c..b5e3e124151 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h @@ -66,7 +66,7 @@ template struct ConstitutiveData { SolidMechanicsDataStateless s_mech_data; - GravityData grav_data; + VolumetricBodyForce volumetric_body_force; }; /// Data that stores intermediate values, which are not needed outside the diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveModels.h index 1ae91962f1f..859250115d2 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveModels.h +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveModels.h @@ -26,13 +26,13 @@ struct ConstitutiveModels TRMProcessData const& process_data, SolidConstitutiveRelation const& solid_material) : s_mech_model(solid_material), - grav_model(process_data.specific_body_force) + gravity_model(process_data.specific_body_force) { } SolidMechanicsModel s_mech_model; SolidDensityModel rho_S_model; - GravityModel grav_model; + GravityModel gravity_model; }; } // namespace ConstitutiveRelations } // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp index b3067c70af9..21c2486f7ad 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp @@ -38,7 +38,7 @@ void ConstitutiveSetting::eval( auto& rho_SR = tmp.rho_SR; auto& s_mech_data = cd.s_mech_data; - auto& grav_data = cd.grav_data; + auto& volumetric_body_force = cd.volumetric_body_force; TemperatureData const T_data{T_ref}; SpaceTimeData const x_t{x_position, t, dt}; @@ -52,8 +52,8 @@ void ConstitutiveSetting::eval( assertEvalArgsUnique(models.rho_S_model); models.rho_S_model.eval(x_t, media_data, T_data, rho_SR); - assertEvalArgsUnique(models.grav_model); - models.grav_model.eval(rho_SR, grav_data); + assertEvalArgsUnique(models.gravity_model); + models.gravity_model.eval(rho_SR, volumetric_body_force); } template struct ConstitutiveSetting<2>; diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.cpp index d5d6fbe6d38..c35e7bcd06c 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.cpp +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.cpp @@ -13,11 +13,9 @@ namespace ProcessLib::SmallDeformation { template void GravityModel::eval( - SolidDensity const& rho_SR, GravityData& out) const + SolidDensity const& rho_SR, VolumetricBodyForce& out) const { - auto const b = specific_body_force_; - - out.volumetric_body_force = *rho_SR * b; + *out = *rho_SR * specific_body_force_; } template struct GravityModel<2>; diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.h index 8ab8883fd3d..e0f93caa19e 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.h +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/Gravity.h @@ -9,15 +9,14 @@ #pragma once +#include "BaseLib/StrongType.h" #include "SolidDensity.h" namespace ProcessLib::SmallDeformation { template -struct GravityData -{ - GlobalDimVector volumetric_body_force; -}; +using VolumetricBodyForce = + BaseLib::StrongType, struct GravityTag>; template struct GravityModel @@ -29,7 +28,7 @@ struct GravityModel } void eval(SolidDensity const& rho_SR, - GravityData& out) const; + VolumetricBodyForce& out) const; private: // TODO (naumov) Do we need to store this for each integration point? diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h index 7e72690a8bf..c9c5899d3d3 100644 --- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h +++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h @@ -278,7 +278,7 @@ class SmallDeformationLocalAssembler this->output_data_[ip]); auto const& sigma = this->current_states_[ip].stress_data.sigma; - auto const& b = CD.grav_data.volumetric_body_force; + auto const& b = *CD.volumetric_body_force; auto const& C = CD.s_mech_data.stiffness_tensor; local_b.noalias() -= From bc984dc6850af84124bad5643e8933f3a6839d67 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Wed, 11 Oct 2023 15:57:28 +0200 Subject: [PATCH 08/10] [PL/SD] Use StrongType for temperature --- ProcessLib/SmallDeformation/ConstitutiveRelations/Base.h | 7 ++----- .../ConstitutiveRelations/ConstitutiveSetting.cpp | 6 +++--- .../ConstitutiveRelations/SolidDensity.cpp | 4 ++-- .../SmallDeformation/ConstitutiveRelations/SolidDensity.h | 3 +-- .../ConstitutiveRelations/SolidMechanics.cpp | 6 +++--- .../ConstitutiveRelations/SolidMechanics.h | 2 +- 6 files changed, 12 insertions(+), 16 deletions(-) diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/Base.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/Base.h index b92cb296b4b..1739e12ebe0 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/Base.h +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/Base.h @@ -9,6 +9,7 @@ #pragma once +#include "BaseLib/StrongType.h" #include "MaterialLib/MPL/Medium.h" #include "MathLib/KelvinVector.h" #include "ParameterLib/SpatialPosition.h" @@ -97,11 +98,7 @@ struct MediaData MaterialPropertyLib::Phase const& solid; }; -template -struct TemperatureData -{ - double T; -}; +using Temperature = BaseLib::StrongType; template struct StrainData diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp index 21c2486f7ad..d0eafec3f1e 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp @@ -40,17 +40,17 @@ void ConstitutiveSetting::eval( auto& s_mech_data = cd.s_mech_data; auto& volumetric_body_force = cd.volumetric_body_force; - TemperatureData const T_data{T_ref}; + Temperature const T{T_ref}; SpaceTimeData const x_t{x_position, t, dt}; MediaData const media_data{medium}; assertEvalArgsUnique(models.s_mech_model); - models.s_mech_model.eval(x_t, T_data, eps_data, eps_data_prev, mat_state, + models.s_mech_model.eval(x_t, T, eps_data, eps_data_prev, mat_state, prev_state.stress_data, state.stress_data, s_mech_data); assertEvalArgsUnique(models.rho_S_model); - models.rho_S_model.eval(x_t, media_data, T_data, rho_SR); + models.rho_S_model.eval(x_t, media_data, T, rho_SR); assertEvalArgsUnique(models.gravity_model); models.gravity_model.eval(rho_SR, volumetric_body_force); diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.cpp index aa39cda70d8..eb743e2fcd3 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.cpp +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.cpp @@ -13,12 +13,12 @@ namespace ProcessLib::SmallDeformation { void SolidDensityModel::eval(SpaceTimeData const& x_t, MediaData const& media_data, - TemperatureData const& T_data, + Temperature const& temperature, SolidDensity& out) const { namespace MPL = MaterialPropertyLib; MPL::VariableArray variables; - variables.temperature = T_data.T; + variables.temperature = *temperature; *out = media_data.solid.property(MPL::PropertyType::density) .template value(variables, x_t.x, x_t.t, x_t.dt); diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.h index 6bfcd8c48c3..6550d6db8b5 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.h +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidDensity.h @@ -19,7 +19,6 @@ using SolidDensity = BaseLib::StrongType; struct SolidDensityModel { void eval(SpaceTimeData const& x_t, MediaData const& media_data, - TemperatureData const& T_data, - SolidDensity& out) const; + Temperature const& temperature, SolidDensity& out) const; }; } // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp index a7b889eaf4b..00eb985ae1f 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp @@ -16,7 +16,7 @@ namespace ConstitutiveRelations template void SolidMechanicsModel::eval( SpaceTimeData const& x_t, - TemperatureData const& T_data, + Temperature const& temperature, StrainData const& eps_data, PrevState> const& eps_data_prev, MaterialStateData& mat_state, @@ -34,7 +34,7 @@ void SolidMechanicsModel::eval( variables.mechanical_strain = eps_data.eps; // external state variables - variables.temperature = T_data.T; + variables.temperature = *temperature; } // previous state @@ -45,7 +45,7 @@ void SolidMechanicsModel::eval( variables_prev.mechanical_strain = eps_data_prev->eps; // external state variables - variables_prev.temperature = T_data.T; + variables_prev.temperature = *temperature; } auto solution = solid_material_.integrateStress( diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h index 1fc9fc4d2f5..572ac2391de 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h @@ -38,7 +38,7 @@ struct SolidMechanicsModel void eval( SpaceTimeData const& x_t, - TemperatureData const& T_data, + Temperature const& temperature, StrainData const& eps_data, PrevState> const& eps_data_prev, MaterialStateData& mat_state, From 19867ce55078dfb024871c7140c2f0d8d4e1ef79 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 13 Oct 2023 13:33:08 +0200 Subject: [PATCH 09/10] [PL/SD] Move free energy density into solid model --- ProcessLib/Deformation/MaterialForces.h | 9 ++--- .../ConstitutiveRelations/ConstitutiveData.h | 5 ++- .../ConstitutiveSetting.cpp | 4 ++- .../ConstitutiveRelations/FreeEnergyDensity.h | 26 ++++++++++++++ .../ConstitutiveRelations/SolidMechanics.cpp | 8 ++++- .../ConstitutiveRelations/SolidMechanics.h | 19 ++++++----- .../LocalAssemblerInterface.h | 6 ---- .../SmallDeformation/SmallDeformationFEM.h | 34 +++---------------- .../SmallDeformationProcess.cpp | 4 --- 9 files changed, 59 insertions(+), 56 deletions(-) create mode 100644 ProcessLib/SmallDeformation/ConstitutiveRelations/FreeEnergyDensity.h diff --git a/ProcessLib/Deformation/MaterialForces.h b/ProcessLib/Deformation/MaterialForces.h index 3852be91f46..89b5c89bbdb 100644 --- a/ProcessLib/Deformation/MaterialForces.h +++ b/ProcessLib/Deformation/MaterialForces.h @@ -37,12 +37,12 @@ template + typename OutputData, typename IntegrationMethod> std::vector const& getMaterialForces( std::vector const& local_x, std::vector& nodal_values, IntegrationMethod const& _integration_method, IPData const& _ip_data, - StressData const& stress_data, MeshLib::Element const& element, - bool const is_axially_symmetric) + StressData const& stress_data, OutputData const& output_data, + MeshLib::Element const& element, bool const is_axially_symmetric) { unsigned const n_integration_points = _integration_method.getNumberOfPoints(); @@ -57,7 +57,8 @@ std::vector const& getMaterialForces( auto const& N = _ip_data[ip].N; auto const& dNdx = _ip_data[ip].dNdx; - auto const& psi = _ip_data[ip].free_energy_density; + auto const& psi = + output_data[ip].free_energy_density_data.free_energy_density; auto const x_coord = NumLib::interpolateXCoordinate( diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h index b5e3e124151..f6525c7d18f 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h @@ -9,6 +9,7 @@ #pragma once +#include "FreeEnergyDensity.h" #include "Gravity.h" #include "ProcessLib/Reflection/ReflectionData.h" #include "SolidDensity.h" @@ -52,12 +53,14 @@ template struct OutputData { StrainData eps_data; + FreeEnergyDensityData free_energy_density_data; static auto reflect() { using Self = OutputData; - return Reflection::reflectWithoutName(&Self::eps_data); + return Reflection::reflectWithoutName(&Self::eps_data, + &Self::free_energy_density_data); } }; diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp index d0eafec3f1e..f5071e1c1a3 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp @@ -40,6 +40,8 @@ void ConstitutiveSetting::eval( auto& s_mech_data = cd.s_mech_data; auto& volumetric_body_force = cd.volumetric_body_force; + auto& free_energy_density_data = out.free_energy_density_data; + Temperature const T{T_ref}; SpaceTimeData const x_t{x_position, t, dt}; MediaData const media_data{medium}; @@ -47,7 +49,7 @@ void ConstitutiveSetting::eval( assertEvalArgsUnique(models.s_mech_model); models.s_mech_model.eval(x_t, T, eps_data, eps_data_prev, mat_state, prev_state.stress_data, state.stress_data, - s_mech_data); + s_mech_data, free_energy_density_data); assertEvalArgsUnique(models.rho_S_model); models.rho_S_model.eval(x_t, media_data, T, rho_SR); diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/FreeEnergyDensity.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/FreeEnergyDensity.h new file mode 100644 index 00000000000..c52ee5ae847 --- /dev/null +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/FreeEnergyDensity.h @@ -0,0 +1,26 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 "Base.h" + +namespace ProcessLib::SmallDeformation +{ +struct FreeEnergyDensityData +{ + double free_energy_density; + + static auto reflect() + { + return ProcessLib::Reflection::reflectWithName( + "free_energy_density", &FreeEnergyDensityData::free_energy_density); + } +}; +} // namespace ProcessLib::SmallDeformation diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp index 00eb985ae1f..283e1957681 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp @@ -22,7 +22,8 @@ void SolidMechanicsModel::eval( MaterialStateData& mat_state, PrevState> const& stress_data_prev, StressData& stress_data, - SolidMechanicsDataStateless& current_stateless) const + SolidMechanicsDataStateless& current_stateless, + FreeEnergyDensityData& free_energy_density_data) const { namespace MPL = MaterialPropertyLib; @@ -59,6 +60,11 @@ void SolidMechanicsModel::eval( std::tie(stress_data.sigma, mat_state.material_state_variables, current_stateless.stiffness_tensor) = std::move(*solution); + + free_energy_density_data.free_energy_density = + solid_material_.computeFreeEnergyDensity( + x_t.t, x_t.x, x_t.dt, eps_data.eps, stress_data.sigma, + *mat_state.material_state_variables); } template struct SolidMechanicsModel<2>; diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h index 572ac2391de..95eafc247ab 100644 --- a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h +++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h @@ -9,6 +9,7 @@ #pragma once +#include "FreeEnergyDensity.h" #include "MaterialLib/SolidModels/MechanicsBase.h" #include "MaterialState.h" #include "StressData.h" @@ -36,15 +37,15 @@ struct SolidMechanicsModel { } - void eval( - SpaceTimeData const& x_t, - Temperature const& temperature, - StrainData const& eps_data, - PrevState> const& eps_data_prev, - MaterialStateData& mat_state, - PrevState> const& stress_data_prev, - StressData& stress_data, - SolidMechanicsDataStateless& current_stateless) const; + void eval(SpaceTimeData const& x_t, + Temperature const& temperature, + StrainData const& eps_data, + PrevState> const& eps_data_prev, + MaterialStateData& mat_state, + PrevState> const& stress_data_prev, + StressData& stress_data, + SolidMechanicsDataStateless& current_stateless, + FreeEnergyDensityData& free_energy_density_data) const; auto getInternalVariables() const { diff --git a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h index 5ac074a06ee..50626d70315 100644 --- a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h +++ b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h @@ -154,12 +154,6 @@ struct SmallDeformationLocalAssemblerInterface } } - virtual std::vector const& getIntPtFreeEnergyDensity( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - // TODO move to NumLib::ExtrapolatableElement unsigned getNumberOfIntegrationPoints() const { diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h index c9c5899d3d3..888044b7fa3 100644 --- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h +++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h @@ -40,8 +40,6 @@ template struct IntegrationPointData final { - double free_energy_density = 0; - double integration_weight; typename ShapeMatricesType::NodalRowVectorType N; typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx; @@ -315,15 +313,6 @@ class SmallDeformationLocalAssembler this->prev_states_[ip], this->material_states_[ip], this->output_data_[ip]); - auto const& eps = this->output_data_[ip].eps_data.eps; - auto const& sigma = this->current_states_[ip].stress_data.sigma; - auto& state = *this->material_states_[ip].material_state_variables; - - // Update free energy density needed for material forces. - _ip_data[ip].free_energy_density = - this->solid_material_.computeFreeEnergyDensity( - t, x_position, dt, eps, sigma, state); - this->material_states_[ip].pushBackState(); } @@ -341,9 +330,10 @@ class SmallDeformationLocalAssembler DisplacementDim, ShapeFunction, ShapeMatricesType, typename BMatricesType::NodalForceVectorType, NodalDisplacementVectorType, GradientVectorType, - GradientMatrixType>( - local_x, nodal_values, this->integration_method_, _ip_data, - this->current_states_, this->element_, this->is_axially_symmetric_); + GradientMatrixType>(local_x, nodal_values, + this->integration_method_, _ip_data, + this->current_states_, this->output_data_, + this->element_, this->is_axially_symmetric_); } Eigen::Map getShapeMatrix( @@ -355,22 +345,6 @@ class SmallDeformationLocalAssembler return Eigen::Map(N.data(), N.size()); } - std::vector const& getIntPtFreeEnergyDensity( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const override - { - cache.clear(); - cache.reserve(_ip_data.size()); - - transform( - cbegin(_ip_data), cend(_ip_data), back_inserter(cache), - [](auto const& ip_data) { return ip_data.free_energy_density; }); - - return cache; - } - private: static constexpr auto localDOF(std::vector const& x) { diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp index e9ae759b7ff..bed27c80494 100644 --- a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp +++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp @@ -104,10 +104,6 @@ void SmallDeformationProcess::initializeConcreteProcess( std::move(get_ip_values_function))); }; - add_secondary_variable("free_energy_density", - 1, - &LocalAssemblerInterface::getIntPtFreeEnergyDensity); - ProcessLib::Reflection::addReflectedSecondaryVariables( SmallDeformationLocalAssemblerInterface< DisplacementDim>::getReflectionDataForOutput(), From 69bbd6f1845c529c763d34750ecf679ce7d6116c Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 13 Oct 2023 14:13:21 +0200 Subject: [PATCH 10/10] [PL/SD] Move SD-specific MaterialForces.h --- ProcessLib/SmallDeformation/LocalAssemblerInterface.h | 2 +- ProcessLib/{Deformation => SmallDeformation}/MaterialForces.h | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) rename ProcessLib/{Deformation => SmallDeformation}/MaterialForces.h (99%) diff --git a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h index 50626d70315..d3e0ee3e1e9 100644 --- a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h +++ b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h @@ -16,11 +16,11 @@ #include "ConstitutiveRelations/ConstitutiveData.h" #include "ConstitutiveRelations/ConstitutiveSetting.h" #include "ConstitutiveRelations/MaterialState.h" +#include "MaterialForces.h" #include "MaterialLib/SolidModels/MechanicsBase.h" #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h" #include "NumLib/Extrapolation/ExtrapolatableElement.h" #include "NumLib/Fem/Integration/GenericIntegrationMethod.h" -#include "ProcessLib/Deformation/MaterialForces.h" #include "ProcessLib/LocalAssemblerInterface.h" #include "ProcessLib/Reflection/ReflectionSetIPData.h" #include "ProcessLib/Utils/SetOrGetIntegrationPointData.h" diff --git a/ProcessLib/Deformation/MaterialForces.h b/ProcessLib/SmallDeformation/MaterialForces.h similarity index 99% rename from ProcessLib/Deformation/MaterialForces.h rename to ProcessLib/SmallDeformation/MaterialForces.h index 89b5c89bbdb..c273669c4a8 100644 --- a/ProcessLib/Deformation/MaterialForces.h +++ b/ProcessLib/SmallDeformation/MaterialForces.h @@ -1,6 +1,5 @@ /** * \file - * * \copyright * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org) * Distributed under a Modified BSD License. @@ -182,7 +181,8 @@ void writeMaterialForces( LocalAssemblerInterface& local_assembler, const NumLib::LocalToGlobalIndexMap& dof_table, GlobalVector const& x, - GlobalVector& node_values) { + GlobalVector& node_values) + { auto const indices = NumLib::getIndices(mesh_item_id, dof_table); std::vector local_data; auto const local_x = x.get(indices);