Commit 094a91d1 authored by Tom Fischer's avatar Tom Fischer
Browse files

Merge branch 'trm-move-code-to-superclass' into 'master'

TRM: Moved some members and methods to the local assembler interface.

See merge request ogs/ogs!4230
parents 04ac8593 aeb873dd
......@@ -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
......@@ -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();