diff --git a/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h
index c22a3cbc758311ea47d09ddc2b89c22a16d976a5..426265183b3371eb7bd609389429c6247a520f71 100644
--- a/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h
+++ b/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h
@@ -11,8 +11,14 @@
#pragma once
#include "MaterialLib/SolidModels/MechanicsBase.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/Integration/GenericIntegrationMethod.h"
#include "ProcessLib/LocalAssemblerInterface.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveSetting.h"
+#include "ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h"
+#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
+#include "ProcessLib/Utils/TransposeInPlace.h"
namespace ProcessLib::ThermoRichardsMechanics
{
@@ -20,88 +26,389 @@ template <int DisplacementDim>
struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
public NumLib::ExtrapolatableElement
{
- virtual std::size_t setIPDataInitialConditions(
- std::string const& name, double const* values,
- int const integration_order) = 0;
-
- virtual std::vector<double> getSigma() const = 0;
-
- virtual std::vector<double> const& getIntPtSigma(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> getSwellingStress() const = 0;
-
- virtual std::vector<double> const& getIntPtSwellingStress(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> getEpsilon() const = 0;
-
- virtual std::vector<double> const& getIntPtEpsilon(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> const& getIntPtDarcyVelocity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> getSaturation() const = 0;
-
- virtual std::vector<double> const& getIntPtSaturation(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> getPorosity() const = 0;
-
- virtual std::vector<double> const& getIntPtPorosity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> getTransportPorosity() const = 0;
-
- virtual std::vector<double> const& getIntPtTransportPorosity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> const& getIntPtDryDensitySolid(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> const& getIntPtLiquidDensity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> const& getIntPtViscosity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
+ LocalAssemblerInterface(
+ MeshLib::Element const& e,
+ NumLib::GenericIntegrationMethod const& integration_method,
+ bool const is_axially_symmetric,
+ ThermoRichardsMechanicsProcessData<DisplacementDim>& 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,
+ e.getID()))
+ {
+ unsigned const n_integration_points =
+ integration_method_.getNumberOfPoints();
+
+ current_states_.resize(n_integration_points);
+ prev_states_.resize(n_integration_points);
+ output_data_.resize(n_integration_points);
+
+ material_states_.reserve(n_integration_points);
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ material_states_.emplace_back(solid_material_);
+ }
+
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(e.getID());
+
+ auto const& medium = process_data_.media_map->getMedium(e.getID());
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ namespace MPL = MaterialPropertyLib;
+
+ x_position.setIntegrationPoint(ip);
+
+ auto& current_state = current_states_[ip];
+
+ // Initial porosity. Could be read from integration point data or
+ // mesh.
+ current_state.poro_data.phi =
+ medium->property(MPL::porosity)
+ .template initialValue<double>(
+ x_position,
+ std::numeric_limits<
+ double>::quiet_NaN() /* t independent */);
+
+ if (medium->hasProperty(MPL::PropertyType::transport_porosity))
+ {
+ current_state.transport_poro_data.phi =
+ medium->property(MPL::transport_porosity)
+ .template initialValue<double>(
+ x_position,
+ std::numeric_limits<
+ double>::quiet_NaN() /* t independent */);
+ }
+ else
+ {
+ current_state.transport_poro_data.phi =
+ current_state.poro_data.phi;
+ }
+ }
+ }
+
+ std::size_t setIPDataInitialConditions(std::string const& name,
+ double const* values,
+ int const integration_order)
+ {
+ if (integration_order !=
+ static_cast<int>(integration_method_.getIntegrationOrder()))
+ {
+ OGS_FATAL(
+ "Setting integration point initial conditions; The integration "
+ "order of the local assembler for element {:d} is different "
+ "from the integration order in the initial condition.",
+ element_.getID());
+ }
+
+ if (name == "sigma_ip")
+ {
+ if (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);
+ }
+ return ProcessLib::setIntegrationPointKelvinVectorData<
+ DisplacementDim>(
+ values, current_states_,
+ [](auto& cs) -> auto& { return cs.s_mech_data.sigma_eff; });
+ }
+
+ if (name == "saturation_ip")
+ {
+ return ProcessLib::setIntegrationPointScalarData(
+ values, current_states_,
+ [](auto& state) -> auto& { return state.S_L_data.S_L; });
+ }
+ if (name == "porosity_ip")
+ {
+ return ProcessLib::setIntegrationPointScalarData(
+ values, current_states_,
+ [](auto& state) -> auto& { return state.poro_data.phi; });
+ }
+ if (name == "transport_porosity_ip")
+ {
+ return ProcessLib::setIntegrationPointScalarData(
+ values, current_states_, [](auto& state) -> auto& {
+ return state.transport_poro_data.phi;
+ });
+ }
+ if (name == "swelling_stress_ip")
+ {
+ return ProcessLib::setIntegrationPointKelvinVectorData<
+ DisplacementDim>(
+ values, current_states_,
+ [](auto& cs) -> auto& { return cs.swelling_data.sigma_sw; });
+ }
+ if (name == "epsilon_ip")
+ {
+ return ProcessLib::setIntegrationPointKelvinVectorData<
+ DisplacementDim>(
+ values, current_states_,
+ [](auto& cs) -> auto& { return cs.eps_data.eps; });
+ }
+ return 0;
+ }
+
+ std::vector<double> getSigma() const
+ {
+ constexpr int kelvin_vector_size =
+ MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+
+ return transposeInPlace<kelvin_vector_size>(
+ [this](std::vector<double>& values)
+ { return getIntPtSigma(0, {}, {}, values); });
+ }
+
+ std::vector<double> const& getIntPtSigma(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+ {
+ return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
+ current_states_, [](auto const& cs) -> auto const& {
+ return cs.s_mech_data.sigma_eff;
+ },
+ cache);
+ }
+
+ std::vector<double> getSwellingStress() const
+ {
+ constexpr int kelvin_vector_size =
+ MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+
+ return transposeInPlace<kelvin_vector_size>(
+ [this](std::vector<double>& values)
+ { return getIntPtSwellingStress(0, {}, {}, values); });
+ }
+
+ std::vector<double> const& getIntPtSwellingStress(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+ {
+ constexpr int kelvin_vector_size =
+ MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+ auto const n_integration_points = current_states_.size();
+
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, kelvin_vector_size, n_integration_points);
+
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ auto const& sigma_sw = current_states_[ip].swelling_data.sigma_sw;
+ cache_mat.col(ip) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_sw);
+ }
+
+ return cache;
+ }
+
+ std::vector<double> getEpsilon() const
+ {
+ constexpr int kelvin_vector_size =
+ MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+
+ return transposeInPlace<kelvin_vector_size>(
+ [this](std::vector<double>& values)
+ { return getIntPtEpsilon(0, {}, {}, values); });
+ }
+
+ std::vector<double> const& getIntPtEpsilon(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+ {
+ return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
+ current_states_,
+ [](auto const& cs) -> auto const& { return cs.eps_data.eps; },
+ cache);
+ }
+
+ std::vector<double> const& getIntPtDarcyVelocity(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+ {
+ unsigned const n_integration_points =
+ integration_method_.getNumberOfPoints();
+
+ cache.clear();
+ auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, DisplacementDim, n_integration_points);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ cache_matrix.col(ip).noalias() =
+ output_data_[ip].darcy_data.v_darcy;
+ }
+
+ return cache;
+ }
+
+ std::vector<double> getSaturation() const
+ {
+ std::vector<double> result;
+ getIntPtSaturation(0, {}, {}, result);
+ return result;
+ }
+
+ std::vector<double> const& getIntPtSaturation(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+ {
+ return ProcessLib::getIntegrationPointScalarData(
+ current_states_,
+ [](auto const& state) -> auto const& { return state.S_L_data.S_L; },
+ cache);
+ }
+
+ std::vector<double> getPorosity() const
+ {
+ std::vector<double> result;
+ getIntPtPorosity(0, {}, {}, result);
+ return result;
+ }
+
+ std::vector<double> const& getIntPtPorosity(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+ {
+ return ProcessLib::getIntegrationPointScalarData(
+ current_states_, [](auto const& state) -> auto const& {
+ return state.poro_data.phi;
+ },
+ cache);
+ }
+
+ std::vector<double> getTransportPorosity() const
+ {
+ std::vector<double> result;
+ getIntPtTransportPorosity(0, {}, {}, result);
+ return result;
+ }
+
+ std::vector<double> const& getIntPtTransportPorosity(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+ {
+ return ProcessLib::getIntegrationPointScalarData(
+ current_states_,
+ [](auto const& state) -> auto const& {
+ return state.transport_poro_data.phi;
+ },
+ cache);
+ }
+
+ std::vector<double> const& getIntPtDryDensitySolid(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+ {
+ return ProcessLib::getIntegrationPointScalarData(
+ output_data_,
+ [](auto const& out) -> auto const& {
+ return out.rho_S_data.dry_density_solid;
+ },
+ cache);
+ }
+
+ std::vector<double> const& getIntPtLiquidDensity(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+ {
+ return ProcessLib::getIntegrationPointScalarData(
+ output_data_,
+ [](auto const& out) -> auto const& {
+ return out.rho_L_data.rho_LR;
+ },
+ cache);
+ }
+
+ std::vector<double> const& getIntPtViscosity(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+ {
+ return ProcessLib::getIntegrationPointScalarData(
+ output_data_, [](auto const& out) -> auto const& {
+ return out.mu_L_data.viscosity;
+ },
+ cache);
+ }
// TODO move to NumLib::ExtrapolatableElement
- virtual unsigned getNumberOfIntegrationPoints() const = 0;
+ unsigned getNumberOfIntegrationPoints() const
+ {
+ return integration_method_.getNumberOfPoints();
+ }
- 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;
+ }
+
+ void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
+ double const /*t*/,
+ double const /*dt*/) override
+ {
+ unsigned const n_integration_points =
+ integration_method_.getNumberOfPoints();
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ // TODO re-evaluate part of the assembly in order to be consistent?
+ material_states_[ip].pushBackState();
+ }
+
+ prev_states_ = current_states_;
+ }
+
+protected:
+ ThermoRichardsMechanicsProcessData<DisplacementDim>& process_data_;
+
+ std::vector<StatefulData<DisplacementDim>>
+ current_states_; // TODO maybe do not store but rather re-evaluate for
+ // state update
+ std::vector<StatefulData<DisplacementDim>> prev_states_;
+
+ // Material state is special, because it contains both the current and the
+ // old state.
+ std::vector<MaterialStateData<DisplacementDim>> material_states_;
+
+ NumLib::GenericIntegrationMethod const& integration_method_;
+ MeshLib::Element const& element_;
+ bool const is_axially_symmetric_;
+
+ MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material_;
+
+ std::vector<OutputData<DisplacementDim>> output_data_;
};
} // namespace ProcessLib::ThermoRichardsMechanics
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
index 90783f6de16f34db8eca8683b9cd6c69bf4c05a7..fa06f3b21868dac80f6f39d5c381a9c0e83e5f23 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
@@ -21,9 +21,9 @@
#include "MaterialLib/MPL/Utils/FormKelvinVectorFromThermalExpansivity.h"
#include "MaterialLib/MPL/Utils/GetLiquidThermalExpansivity.h"
#include "MaterialLib/PhysicalConstant.h"
-#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
#include "MathLib/KelvinVector.h"
#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
#include "ProcessLib/Utils/TransposeInPlace.h"
#include "ThermoRichardsMechanicsFEM.h"
@@ -42,50 +42,34 @@ ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
NumLib::GenericIntegrationMethod const& integration_method,
bool const is_axially_symmetric,
ThermoRichardsMechanicsProcessData<DisplacementDim>& 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, e.getID()))
+ : LocalAssemblerInterface<DisplacementDim>(
+ e, integration_method, is_axially_symmetric, process_data)
{
unsigned const n_integration_points =
- integration_method_.getNumberOfPoints();
+ this->integration_method_.getNumberOfPoints();
- current_states_.resize(n_integration_points);
- prev_states_.resize(n_integration_points);
ip_data_.resize(n_integration_points);
- secondary_data_.N_u.resize(n_integration_points);
- output_data_.resize(n_integration_points);
-
- material_states_.reserve(n_integration_points);
- for (unsigned ip = 0; ip < n_integration_points; ++ip)
- {
- material_states_.emplace_back(solid_material_);
- }
auto const shape_matrices_u =
NumLib::initShapeMatrices<ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement,
DisplacementDim>(e, is_axially_symmetric,
- integration_method_);
+ integration_method);
auto const shape_matrices =
NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
DisplacementDim>(e, is_axially_symmetric,
- integration_method_);
-
- auto const& medium = process_data_.media_map->getMedium(element_.getID());
+ integration_method);
ParameterLib::SpatialPosition x_position;
- x_position.setElementID(element_.getID());
+ x_position.setElementID(e.getID());
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
auto& ip_data = ip_data_[ip];
auto const& sm_u = shape_matrices_u[ip];
ip_data_[ip].integration_weight =
- integration_method_.getWeightedPoint(ip).getWeight() *
+ integration_method.getWeightedPoint(ip).getWeight() *
sm_u.integralMeasure * sm_u.detJ;
ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
@@ -105,101 +89,9 @@ ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
// ip_data.N_p and ip_data.dNdx_p are used for both p and T variables
ip_data.N_p = shape_matrices[ip].N;
ip_data.dNdx_p = shape_matrices[ip].dNdx;
-
- auto& current_state = current_states_[ip];
-
- // Initial porosity. Could be read from integration point data or mesh.
- current_state.poro_data.phi =
- medium->property(MPL::porosity)
- .template initialValue<double>(
- x_position,
- std::numeric_limits<
- double>::quiet_NaN() /* t independent */);
-
- if (medium->hasProperty(MPL::PropertyType::transport_porosity))
- {
- current_state.transport_poro_data.phi =
- medium->property(MPL::transport_porosity)
- .template initialValue<double>(
- x_position,
- std::numeric_limits<
- double>::quiet_NaN() /* t independent */);
- }
- else
- {
- current_state.transport_poro_data.phi = current_state.poro_data.phi;
- }
-
- secondary_data_.N_u[ip] = shape_matrices_u[ip].N;
}
}
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::size_t ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction,
- DisplacementDim>::setIPDataInitialConditions(std::string const& name,
- double const* values,
- int const integration_order)
-{
- if (integration_order !=
- static_cast<int>(integration_method_.getIntegrationOrder()))
- {
- OGS_FATAL(
- "Setting integration point initial conditions; The integration "
- "order of the local assembler for element {:d} is different "
- "from the integration order in the initial condition.",
- element_.getID());
- }
-
- if (name == "sigma_ip")
- {
- if (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);
- }
- return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
- values, current_states_,
- [](auto& cs) -> auto& { return cs.s_mech_data.sigma_eff; });
- }
-
- if (name == "saturation_ip")
- {
- return ProcessLib::setIntegrationPointScalarData(
- values, current_states_,
- [](auto& state) -> auto& { return state.S_L_data.S_L; });
- }
- if (name == "porosity_ip")
- {
- return ProcessLib::setIntegrationPointScalarData(
- values, current_states_,
- [](auto& state) -> auto& { return state.poro_data.phi; });
- }
- if (name == "transport_porosity_ip")
- {
- return ProcessLib::setIntegrationPointScalarData(
- values, current_states_,
- [](auto& state) -> auto& { return state.transport_poro_data.phi; });
- }
- if (name == "swelling_stress_ip")
- {
- return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
- values, current_states_,
- [](auto& cs) -> auto& { return cs.swelling_data.sigma_sw; });
- }
- if (name == "epsilon_ip")
- {
- return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
- values, current_states_,
- [](auto& cs) -> auto& { return cs.eps_data.eps; });
- }
- return 0;
-}
-
template <typename ShapeFunctionDisplacement, typename ShapeFunction,
int DisplacementDim>
void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
@@ -221,16 +113,17 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
temperature_size);
constexpr double dt = std::numeric_limits<double>::quiet_NaN();
- auto const& medium = process_data_.media_map->getMedium(element_.getID());
+ auto const& medium =
+ this->process_data_.media_map->getMedium(this->element_.getID());
MPL::VariableArray variables;
ParameterLib::SpatialPosition x_position;
- x_position.setElementID(element_.getID());
+ x_position.setElementID(this->element_.getID());
auto const& solid_phase = medium->phase("Solid");
unsigned const n_integration_points =
- integration_method_.getNumberOfPoints();
+ this->integration_method_.getNumberOfPoints();
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
@@ -249,7 +142,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
NumLib::shapeFunctionInterpolate(T, N, T_ip);
variables[static_cast<int>(MPL::Variable::temperature)] = T_ip;
- prev_states_[ip].S_L_data.S_L =
+ this->prev_states_[ip].S_L_data.S_L =
medium->property(MPL::PropertyType::saturation)
.template value<double>(variables, x_position, t, dt);
@@ -257,12 +150,12 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
// restart.
SpaceTimeData const x_t{x_position, t, dt};
ElasticTangentStiffnessData<DisplacementDim> C_el_data;
- ElasticTangentStiffnessModel<DisplacementDim>{solid_material_}.eval(
- x_t, {T_ip, 0, {}}, C_el_data);
+ ElasticTangentStiffnessModel<DisplacementDim>{this->solid_material_}
+ .eval(x_t, {T_ip, 0, {}}, C_el_data);
- auto const& eps = current_states_[ip].eps_data.eps;
- auto const& sigma_sw = current_states_[ip].swelling_data.sigma_sw;
- prev_states_[ip].s_mech_data.eps_m.noalias() =
+ auto const& eps = this->current_states_[ip].eps_data.eps;
+ auto const& sigma_sw = this->current_states_[ip].swelling_data.sigma_sw;
+ this->prev_states_[ip].s_mech_data.eps_m.noalias() =
solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
? eps + C_el_data.C_el.inverse() * sigma_sw
: eps;
@@ -281,30 +174,33 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
std::vector<double>& local_rhs_data,
std::vector<double>& local_Jac_data)
{
- auto& medium = *process_data_.media_map->getMedium(element_.getID());
+ auto& medium =
+ *this->process_data_.media_map->getMedium(this->element_.getID());
ParameterLib::SpatialPosition x_position;
- x_position.setElementID(element_.getID());
+ x_position.setElementID(this->element_.getID());
LocalMatrices loc_mat;
loc_mat.setZero();
LocalMatrices loc_mat_current_ip;
loc_mat_current_ip.setZero(); // only to set the right matrix sizes
- ConstitutiveSetting<DisplacementDim> constitutive_setting(solid_material_,
- process_data_);
+ ConstitutiveSetting<DisplacementDim> constitutive_setting(
+ this->solid_material_, this->process_data_);
- for (unsigned ip = 0; ip < integration_method_.getNumberOfPoints(); ++ip)
+ for (unsigned ip = 0; ip < this->integration_method_.getNumberOfPoints();
+ ++ip)
{
x_position.setIntegrationPoint(ip);
- assembleWithJacobianSingleIP(t, dt, x_position, //
- local_x, local_xdot, //
- ip_data_[ip], constitutive_setting,
- medium, //
- loc_mat_current_ip, //
- current_states_[ip], prev_states_[ip],
- material_states_[ip], output_data_[ip]);
+ assembleWithJacobianSingleIP(
+ t, dt, x_position, //
+ local_x, local_xdot, //
+ ip_data_[ip], constitutive_setting,
+ medium, //
+ loc_mat_current_ip, //
+ this->current_states_[ip], this->prev_states_[ip],
+ this->material_states_[ip], this->output_data_[ip]);
loc_mat += loc_mat_current_ip;
}
@@ -322,7 +218,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
ShapeFunctionDisplacement, ShapeFunction,
DisplacementDim>::LocalMatrices& loc_mat) const
{
- if (process_data_.apply_mass_lumping)
+ if (this->process_data_.apply_mass_lumping)
{
loc_mat.storage_p_a_p =
loc_mat.storage_p_a_p.colwise().sum().eval().asDiagonal();
@@ -421,20 +317,21 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const x_coord =
NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
- ShapeMatricesTypeDisplacement>(element_,
- N_u);
+ ShapeMatricesTypeDisplacement>(
+ this->element_, N_u);
auto const B =
LinearBMatrix::computeBMatrix<DisplacementDim,
ShapeFunctionDisplacement::NPOINTS,
typename BMatricesType::BMatrixType>(
- dNdx_u, N_u, x_coord, is_axially_symmetric_);
+ dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
auto const [T, p_L, u] = localDOF(local_x);
auto const [T_dot, p_L_dot, u_dot] = localDOF(local_xdot);
GlobalDimVectorType const grad_T_ip = dNdx * T;
- ConstitutiveModels<DisplacementDim> models(process_data_, solid_material_);
+ ConstitutiveModels<DisplacementDim> models(this->process_data_,
+ this->solid_material_);
ConstitutiveTempData<DisplacementDim> tmp;
ConstitutiveData<DisplacementDim> CD;
@@ -546,265 +443,6 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
out *= ip_data.integration_weight;
}
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction, DisplacementDim>::getSigma() const
-{
- constexpr int kelvin_vector_size =
- MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-
- return transposeInPlace<kelvin_vector_size>(
- [this](std::vector<double>& values)
- { return getIntPtSigma(0, {}, {}, values); });
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction, DisplacementDim>::
- getIntPtSigma(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
- current_states_,
- [](auto const& cs) -> auto const& { return cs.s_mech_data.sigma_eff; },
- cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction,
- DisplacementDim>::getSwellingStress() const
-{
- constexpr int kelvin_vector_size =
- MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-
- return transposeInPlace<kelvin_vector_size>(
- [this](std::vector<double>& values)
- { return getIntPtSwellingStress(0, {}, {}, values); });
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction, DisplacementDim>::
- getIntPtSwellingStress(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- constexpr int kelvin_vector_size =
- MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
- auto const n_integration_points = ip_data_.size();
-
- cache.clear();
- auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
- double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, kelvin_vector_size, n_integration_points);
-
- for (unsigned ip = 0; ip < n_integration_points; ++ip)
- {
- auto const& sigma_sw = current_states_[ip].swelling_data.sigma_sw;
- cache_mat.col(ip) =
- MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_sw);
- }
-
- return cache;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double>
-ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
- DisplacementDim>::getEpsilon() const
-{
- constexpr int kelvin_vector_size =
- MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-
- return transposeInPlace<kelvin_vector_size>(
- [this](std::vector<double>& values)
- { return getIntPtEpsilon(0, {}, {}, values); });
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction, DisplacementDim>::
- getIntPtEpsilon(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
- current_states_,
- [](auto const& cs) -> auto const& { return cs.eps_data.eps; }, cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction, DisplacementDim>::
- getIntPtDarcyVelocity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- unsigned const n_integration_points =
- integration_method_.getNumberOfPoints();
-
- cache.clear();
- auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
- double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, DisplacementDim, n_integration_points);
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- cache_matrix.col(ip).noalias() = output_data_[ip].darcy_data.v_darcy;
- }
-
- return cache;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double>
-ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
- DisplacementDim>::getSaturation() const
-{
- std::vector<double> result;
- getIntPtSaturation(0, {}, {}, result);
- return result;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction, DisplacementDim>::
- getIntPtSaturation(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- return ProcessLib::getIntegrationPointScalarData(
- current_states_,
- [](auto const& state) -> auto const& { return state.S_L_data.S_L; },
- cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double>
-ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
- DisplacementDim>::getPorosity() const
-{
- std::vector<double> result;
- getIntPtPorosity(0, {}, {}, result);
- return result;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction, DisplacementDim>::
- getIntPtPorosity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- return ProcessLib::getIntegrationPointScalarData(
- current_states_,
- [](auto const& state) -> auto const& { return state.poro_data.phi; },
- cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction,
- DisplacementDim>::getTransportPorosity() const
-{
- std::vector<double> result;
- getIntPtTransportPorosity(0, {}, {}, result);
- return result;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction, DisplacementDim>::
- getIntPtTransportPorosity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- return ProcessLib::getIntegrationPointScalarData(
- current_states_, [](auto const& state) -> auto const& {
- return state.transport_poro_data.phi;
- },
- cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction, DisplacementDim>::
- getIntPtDryDensitySolid(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- return ProcessLib::getIntegrationPointScalarData(
- output_data_,
- [](auto const& out) -> auto const& {
- return out.rho_S_data.dry_density_solid;
- },
- cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction, DisplacementDim>::
- getIntPtLiquidDensity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- return ProcessLib::getIntegrationPointScalarData(
- output_data_,
- [](auto const& out) -> auto const& { return out.rho_L_data.rho_LR; },
- cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction, DisplacementDim>::
- getIntPtViscosity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- return ProcessLib::getIntegrationPointScalarData(
- output_data_,
- [](auto const& out) -> auto const& { return out.mu_L_data.viscosity; },
- cache);
-}
-
template <typename ShapeFunctionDisplacement, typename ShapeFunction,
int DisplacementDim>
void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
@@ -821,14 +459,15 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const p_L_dot = block_p(local_x_dot);
auto const u_dot = block_u(local_x_dot);
- auto const e_id = element_.getID();
- auto& medium = *process_data_.media_map->getMedium(e_id);
+ auto const e_id = this->element_.getID();
+ auto const& process_data = this->process_data_;
+ auto& medium = *process_data.media_map->getMedium(e_id);
ParameterLib::SpatialPosition x_position;
x_position.setElementID(e_id);
unsigned const n_integration_points =
- integration_method_.getNumberOfPoints();
+ this->integration_method_.getNumberOfPoints();
double saturation_avg = 0;
double porosity_avg = 0;
@@ -838,15 +477,19 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
KV sigma_avg = KV::Zero();
- ConstitutiveSetting<DisplacementDim> constitutive_setting(solid_material_,
- process_data_);
+ ConstitutiveSetting<DisplacementDim> constitutive_setting(
+ this->solid_material_, process_data);
- ConstitutiveModels<DisplacementDim> models(process_data_, solid_material_);
+ ConstitutiveModels<DisplacementDim> models(process_data,
+ this->solid_material_);
ConstitutiveTempData<DisplacementDim> tmp;
ConstitutiveData<DisplacementDim> CD;
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
+ auto& current_state = this->current_states_[ip];
+ auto& output_data = this->output_data_[ip];
+
x_position.setIntegrationPoint(ip);
auto const& ip_data = ip_data_[ip];
@@ -860,12 +503,12 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const x_coord =
NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement>(
- element_, N_u);
+ this->element_, N_u);
auto const B =
LinearBMatrix::computeBMatrix<DisplacementDim,
ShapeFunctionDisplacement::NPOINTS,
typename BMatricesType::BMatrixType>(
- dNdx_u, N_u, x_coord, is_axially_symmetric_);
+ dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
double const T_ip = N * T;
double const T_dot_ip = N * T_dot;
@@ -880,21 +523,21 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
// rate computation
KelvinVectorType eps_prev = eps - B * (u_dot * dt);
- constitutive_setting.eval(models, //
- t, dt, x_position, //
- medium, //
- {T_ip, T_dot_ip, grad_T_ip}, //
- {p_cap_ip, p_cap_dot_ip, grad_p_cap_ip}, //
- eps, eps_prev, current_states_[ip],
- prev_states_[ip], material_states_[ip], tmp,
- output_data_[ip], CD);
-
- saturation_avg += current_states_[ip].S_L_data.S_L;
- porosity_avg += current_states_[ip].poro_data.phi;
-
- liquid_density_avg += output_data_[ip].rho_L_data.rho_LR;
- viscosity_avg += output_data_[ip].mu_L_data.viscosity;
- sigma_avg += current_states_[ip].s_mech_data.sigma_eff;
+ constitutive_setting.eval(
+ models, //
+ t, dt, x_position, //
+ medium, //
+ {T_ip, T_dot_ip, grad_T_ip}, //
+ {p_cap_ip, p_cap_dot_ip, grad_p_cap_ip}, //
+ eps, eps_prev, current_state, this->prev_states_[ip],
+ this->material_states_[ip], tmp, output_data, CD);
+
+ saturation_avg += current_state.S_L_data.S_L;
+ porosity_avg += current_state.poro_data.phi;
+
+ liquid_density_avg += output_data.rho_L_data.rho_LR;
+ viscosity_avg += output_data.mu_L_data.viscosity;
+ sigma_avg += current_state.s_mech_data.sigma_eff;
}
saturation_avg /= n_integration_points;
porosity_avg /= n_integration_points;
@@ -902,43 +545,23 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
liquid_density_avg /= n_integration_points;
sigma_avg /= n_integration_points;
- (*process_data_.element_saturation)[e_id] = saturation_avg;
- (*process_data_.element_porosity)[e_id] = porosity_avg;
- (*process_data_.element_liquid_density)[e_id] = liquid_density_avg;
- (*process_data_.element_viscosity)[e_id] = viscosity_avg;
+ (*process_data.element_saturation)[e_id] = saturation_avg;
+ (*process_data.element_porosity)[e_id] = porosity_avg;
+ (*process_data.element_liquid_density)[e_id] = liquid_density_avg;
+ (*process_data.element_viscosity)[e_id] = viscosity_avg;
Eigen::Map<KV>(
- &(*process_data_.element_stresses)[e_id * KV::RowsAtCompileTime]) =
+ &(*process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) =
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_avg);
NumLib::interpolateToHigherOrderNodes<
ShapeFunction, typename ShapeFunctionDisplacement::MeshElement,
- DisplacementDim>(element_, is_axially_symmetric_, p_L,
- *process_data_.pressure_interpolated);
+ DisplacementDim>(this->element_, this->is_axially_symmetric_, p_L,
+ *process_data.pressure_interpolated);
NumLib::interpolateToHigherOrderNodes<
ShapeFunction, typename ShapeFunctionDisplacement::MeshElement,
- DisplacementDim>(element_, is_axially_symmetric_, T,
- *process_data_.temperature_interpolated);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-unsigned ThermoRichardsMechanicsLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunction,
- DisplacementDim>::getNumberOfIntegrationPoints() const
-{
- return integration_method_.getNumberOfPoints();
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunction,
- int DisplacementDim>
-typename MaterialLib::Solids::MechanicsBase<
- DisplacementDim>::MaterialStateVariables const&
-ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
- DisplacementDim>::
- getMaterialStateVariablesAt(unsigned integration_point) const
-{
- return *material_states_[integration_point].material_state_variables;
+ DisplacementDim>(this->element_, this->is_axially_symmetric_, T,
+ *process_data.temperature_interpolated);
}
} // namespace ThermoRichardsMechanics
} // namespace ProcessLib
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h
index 88393a3ea427abcc0bb5a64913e59c6d1940d80b..e90d17447464798b97df90ccb7953216f124155b 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h
@@ -16,17 +16,11 @@
#include "ConstitutiveSetting.h"
#include "IntegrationPointData.h"
#include "LocalAssemblerInterface.h"
-#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
#include "MathLib/KelvinVector.h"
-#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
-#include "NumLib/DOF/DOFTableUtil.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"
-#include "ProcessLib/Deformation/LinearBMatrix.h"
-#include "ProcessLib/LocalAssemblerTraits.h"
#include "ThermoRichardsMechanicsProcessData.h"
namespace ProcessLib
@@ -35,14 +29,6 @@ namespace ThermoRichardsMechanics
{
namespace MPL = MaterialPropertyLib;
-/// Used for the extrapolation of the integration point values. It is ordered
-/// (and stored) by integration points.
-template <typename ShapeMatrixType>
-struct SecondaryData
-{
- std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N_u;
-};
-
template <typename ShapeFunctionDisplacement, typename ShapeFunction,
int DisplacementDim>
class ThermoRichardsMechanicsLocalAssembler
@@ -93,12 +79,6 @@ public:
bool const is_axially_symmetric,
ThermoRichardsMechanicsProcessData<DisplacementDim>& process_data);
- /// \return the number of read integration points.
- std::size_t setIPDataInitialConditions(
- std::string const& name,
- double const* values,
- int const integration_order) override;
-
void setInitialConditionsConcrete(std::vector<double> const& local_x,
double const t,
bool const use_monolithic_scheme,
@@ -274,48 +254,32 @@ public:
void initializeConcrete() override
{
unsigned const n_integration_points =
- integration_method_.getNumberOfPoints();
+ this->integration_method_.getNumberOfPoints();
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
// Set initial stress from parameter.
- if (process_data_.initial_stress != nullptr)
+ if (this->process_data_.initial_stress != nullptr)
{
ParameterLib::SpatialPosition const x_position{
- std::nullopt, element_.getID(), ip,
+ std::nullopt, this->element_.getID(), ip,
MathLib::Point3d(NumLib::interpolateCoordinates<
ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement>(
- element_, ip_data_[ip].N_u))};
+ this->element_, ip_data_[ip].N_u))};
- current_states_[ip].s_mech_data.sigma_eff =
+ this->current_states_[ip].s_mech_data.sigma_eff =
MathLib::KelvinVector::symmetricTensorToKelvinVector<
- DisplacementDim>((*process_data_.initial_stress)(
+ DisplacementDim>((*this->process_data_.initial_stress)(
std::numeric_limits<
double>::quiet_NaN() /* time independent */,
x_position));
}
- material_states_[ip].pushBackState();
+ this->material_states_[ip].pushBackState();
}
- prev_states_ = current_states_;
- }
-
- void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
- double const /*t*/,
- double const /*dt*/) override
- {
- unsigned const n_integration_points =
- integration_method_.getNumberOfPoints();
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- // TODO re-evaluate part of the assembly in order to be consistent?
- material_states_[ip].pushBackState();
- }
-
- prev_states_ = current_states_;
+ this->prev_states_ = this->current_states_;
}
void computeSecondaryVariableConcrete(
@@ -325,111 +289,15 @@ public:
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
- auto const& N_u = secondary_data_.N_u[integration_point];
+ auto const& N_u = ip_data_[integration_point].N_u;
// assumes N is stored contiguously in memory
return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
}
- std::vector<double> getSigma() const override;
-
- std::vector<double> const& getIntPtDarcyVelocity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
-
- std::vector<double> getSaturation() const override;
- std::vector<double> const& getIntPtSaturation(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
-
- std::vector<double> getPorosity() const override;
- std::vector<double> const& getIntPtPorosity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
-
- std::vector<double> getTransportPorosity() const override;
- std::vector<double> const& getIntPtTransportPorosity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
-
- std::vector<double> const& getIntPtSigma(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
-
- std::vector<double> getSwellingStress() const override;
- std::vector<double> const& getIntPtSwellingStress(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
-
- std::vector<double> getEpsilon() const override;
- std::vector<double> const& getIntPtEpsilon(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
-
- std::vector<double> const& getIntPtDryDensitySolid(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
-
- std::vector<double> const& getIntPtLiquidDensity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
-
- std::vector<double> const& getIntPtViscosity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
-
-private:
- unsigned getNumberOfIntegrationPoints() const override;
-
- typename MaterialLib::Solids::MechanicsBase<
- DisplacementDim>::MaterialStateVariables const&
- getMaterialStateVariablesAt(unsigned integration_point) const override;
-
private:
- ThermoRichardsMechanicsProcessData<DisplacementDim>& process_data_;
-
- std::vector<StatefulData<DisplacementDim>>
- current_states_; // TODO maybe do not store but rather re-evaluate for
- // state update
- std::vector<StatefulData<DisplacementDim>> prev_states_;
-
- // Material state is special, because it contains both the current and the
- // old state.
- std::vector<MaterialStateData<DisplacementDim>> material_states_;
-
std::vector<IpData> ip_data_;
- NumLib::GenericIntegrationMethod const& integration_method_;
- MeshLib::Element const& element_;
- bool const is_axially_symmetric_;
- SecondaryData<
- typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
- secondary_data_;
-
- MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material_;
-
- std::vector<OutputData<DisplacementDim>> output_data_;
-
static auto block_uu(auto& mat)
{
return mat.template block<displacement_size, displacement_size>(
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
index 0c495508502fe2032e82b5351af9eb619993e8fb..9d48203acffd20e50d9bd0caca32b34fdcbafe38 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
@@ -14,7 +14,7 @@
#include "BaseLib/Error.h"
#include "MeshLib/Elements/Utils.h"
-#include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "NumLib/DOF/DOFTableUtil.h"
#include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
#include "ProcessLib/Process.h"
#include "ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h"