Revision c8a94924e1587005284bada9a0eac0af0d37746e authored by Dmitry Yu. Naumov on 13 July 2021, 14:14:43 UTC, committed by Dmitry Yu. Naumov on 13 July 2021, 14:14:43 UTC
Enabled using MFront in HydroMechanics Closes #3162 See merge request ogs/ogs!3704
HydroMechanicsFEM.h
/**
* \file
* \copyright
* Copyright (c) 2012-2021, 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 <memory>
#include <vector>
#include "HydroMechanicsProcessData.h"
#include "LocalAssemblerInterface.h"
#include "MaterialLib/PhysicalConstant.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/ShapeMatrixPolicy.h"
#include "ParameterLib/Parameter.h"
#include "ProcessLib/Deformation/BMatrixPolicy.h"
#include "ProcessLib/Deformation/LinearBMatrix.h"
#include "ProcessLib/LocalAssemblerTraits.h"
#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
namespace ProcessLib
{
namespace HydroMechanics
{
namespace MPL = MaterialPropertyLib;
template <typename BMatricesType, typename ShapeMatrixTypeDisplacement,
typename ShapeMatricesTypePressure, int DisplacementDim, int NPoints>
struct IntegrationPointData final
{
explicit IntegrationPointData(
MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
solid_material)
: solid_material(solid_material),
material_state_variables(
solid_material.createMaterialStateVariables())
{
}
typename ShapeMatrixTypeDisplacement::template MatrixType<
DisplacementDim, NPoints * DisplacementDim>
N_u_op; /**< for interpolation of the displacement vector, whereas N_u interpolates one component (scalar) */
typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
typename BMatricesType::KelvinVectorType eps, eps_prev;
typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
DisplacementDim>::MaterialStateVariables>
material_state_variables;
double integration_weight;
// TODO disable in monolithic scheme to save memory
double coupling_pressure = std::numeric_limits<double>::quiet_NaN(); /**<
needed for staggered scheme to store value from last coupling iteration */
void pushBackState()
{
eps_prev = eps;
sigma_eff_prev = sigma_eff;
material_state_variables->pushBackState();
}
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>
computeElasticTangentStiffness(
double const t,
ParameterLib::SpatialPosition const& x_position,
double const dt,
double const temperature)
{
namespace MPL = MaterialPropertyLib;
MPL::VariableArray variable_array;
MPL::VariableArray variable_array_prev;
auto const null_state = solid_material.createMaterialStateVariables();
using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
variable_array[static_cast<int>(MPL::Variable::stress)].emplace<KV>(
KV::Zero());
variable_array[static_cast<int>(MPL::Variable::mechanical_strain)]
.emplace<KV>(KV::Zero());
variable_array[static_cast<int>(MPL::Variable::temperature)]
.emplace<double>(temperature);
variable_array_prev[static_cast<int>(MPL::Variable::stress)]
.emplace<KV>(KV::Zero());
variable_array_prev[static_cast<int>(MPL::Variable::mechanical_strain)]
.emplace<KV>(KV::Zero());
variable_array_prev[static_cast<int>(MPL::Variable::temperature)]
.emplace<double>(temperature);
auto&& solution =
solid_material.integrateStress(variable_array_prev, variable_array,
t, x_position, dt, *null_state);
if (!solution)
{
OGS_FATAL("Computation of elastic tangent stiffness failed.");
}
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C =
std::move(std::get<2>(*solution));
return C;
}
template <typename DisplacementVectorType>
typename BMatricesType::KelvinMatrixType updateConstitutiveRelation(
MPL::VariableArray const& variable_array,
double const t,
ParameterLib::SpatialPosition const& x_position,
double const dt,
DisplacementVectorType const& /*u*/,
double const T)
{
MaterialPropertyLib::VariableArray variable_array_prev;
variable_array_prev[static_cast<int>(
MaterialPropertyLib::Variable::stress)]
.emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
sigma_eff_prev);
variable_array_prev[static_cast<int>(MaterialPropertyLib::Variable::
mechanical_strain)]
.emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
eps_prev);
variable_array_prev[static_cast<int>(
MaterialPropertyLib::Variable::temperature)]
.emplace<double>(T);
auto&& solution = solid_material.integrateStress(
variable_array_prev, variable_array, t, x_position, dt,
*material_state_variables);
if (!solution)
{
OGS_FATAL("Computation of local constitutive relation failed.");
}
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C;
std::tie(sigma_eff, material_state_variables, C) = std::move(*solution);
return C;
}
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
/// 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 ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
class HydroMechanicsLocalAssembler
: public LocalAssemblerInterface<DisplacementDim>
{
public:
using ShapeMatricesTypeDisplacement =
ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
// Types for pressure.
using ShapeMatricesTypePressure =
ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>;
static int const KelvinVectorSize =
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler const&) = delete;
HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler&&) = delete;
HydroMechanicsLocalAssembler(
MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
bool const is_axially_symmetric,
unsigned const integration_order,
HydroMechanicsProcessData<DisplacementDim>& process_data);
/// Returns number of read integration points.
std::size_t setIPDataInitialConditions(
std::string const& name,
double const* values,
int const integration_order) override;
void assemble(double const /*t*/, double const /*dt*/,
std::vector<double> const& /*local_x*/,
std::vector<double> const& /*local_xdot*/,
std::vector<double>& /*local_M_data*/,
std::vector<double>& /*local_K_data*/,
std::vector<double>& /*local_rhs_data*/) override
{
OGS_FATAL(
"HydroMechanicsLocalAssembler: assembly without jacobian is not "
"implemented.");
}
void assembleWithJacobian(double const t, double const dt,
std::vector<double> const& local_x,
std::vector<double> const& local_xdot,
const double /*dxdot_dx*/, const double /*dx_dx*/,
std::vector<double>& /*local_M_data*/,
std::vector<double>& /*local_K_data*/,
std::vector<double>& local_rhs_data,
std::vector<double>& local_Jac_data) override;
void assembleWithJacobianForStaggeredScheme(
const double t, double const dt, Eigen::VectorXd const& local_x,
Eigen::VectorXd const& local_xdot, const double dxdot_dx,
const double dx_dx, int const process_id,
std::vector<double>& local_M_data, std::vector<double>& local_K_data,
std::vector<double>& local_b_data,
std::vector<double>& local_Jac_data) override;
void initializeConcrete() override
{
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
auto& ip_data = _ip_data[ip];
/// Set initial stress from parameter.
if (_process_data.initial_stress != nullptr)
{
ParameterLib::SpatialPosition const x_position{
std::nullopt, _element.getID(), ip,
MathLib::Point3d(NumLib::interpolateCoordinates<
ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement>(
_element, ip_data.N_u))};
ip_data.sigma_eff =
MathLib::KelvinVector::symmetricTensorToKelvinVector<
DisplacementDim>((*_process_data.initial_stress)(
std::numeric_limits<
double>::quiet_NaN() /* time independent */,
x_position));
}
ip_data.pushBackState();
}
}
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++)
{
_ip_data[ip].pushBackState();
}
}
void computeSecondaryVariableConcrete(
double const t, double const dt, Eigen::VectorXd const& local_xs,
Eigen::VectorXd const& local_x_dot) override;
void postNonLinearSolverConcrete(std::vector<double> const& local_x,
std::vector<double> const& local_xdot,
double const t, double const dt,
bool const use_monolithic_scheme,
int const process_id) override;
void setInitialConditionsConcrete(std::vector<double> const& local_x,
double const t,
bool const use_monolithic_scheme,
int const process_id) override;
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
auto const& N_u = _secondary_data.N_u[integration_point];
// assumes N is stored contiguously in memory
return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
}
// 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<double> getSigma() const override;
std::vector<double> getEpsilon() 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> const& getIntPtSigma(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& cache) const override
{
return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
_ip_data, &IpData::sigma_eff, cache);
}
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
{
return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
_ip_data, &IpData::eps, cache);
}
private:
/**
* Assemble local matrices and vectors arise from the linearized discretized
* weak form of the residual of the momentum balance equation,
* \f[
* \nabla (\sigma - \alpha_b p \mathrm{I}) = f
* \f]
* where \f$ \sigma\f$ is the effective stress tensor, \f$p\f$ is the pore
* pressure, \f$\alpha_b\f$ is the Biot constant, \f$\mathrm{I}\f$ is the
* identity tensor, and \f$f\f$ is the body force.
*
* @param t Time
* @param dt Time increment
* @param local_x Nodal values of \f$x\f$ of an element of all
* coupled processes.
* @param local_b_data Right hand side vector of an element.
* @param local_Jac_data Element Jacobian matrix for the Newton-Raphson
* method.
*/
void assembleWithJacobianForDeformationEquations(
const double t, double const dt, Eigen::VectorXd const& local_x,
std::vector<double>& local_b_data, std::vector<double>& local_Jac_data);
/**
* Assemble local matrices and vectors arise from the linearized discretized
* weak form of the residual of the mass balance equation of single phase
* flow,
* \f[
* \alpha \cdot{p} - \nabla (K (\nabla p + \rho g \nabla z) +
* \alpha_b \nabla \cdot \dot{u} = Q
* \f]
* where \f$ alpha\f$ is a coefficient may depend on storage or the fluid
* density change, \f$ \rho\f$ is the fluid density, \f$g\f$ is the
* gravitational acceleration, \f$z\f$ is the vertical coordinate, \f$u\f$
* is the displacement, and \f$Q\f$ is the source/sink term.
*
* @param t Time
* @param dt Time increment
* @param local_x Nodal values of \f$x\f$ of an element of all
* coupled processes.
* @param local_xdot Nodal values of \f$\dot{x}\f$ of an element of
* all coupled processes.
* @param local_b_data Right hand side vector of an element.
* @param local_Jac_data Element Jacobian matrix for the Newton-Raphson
* method.
*/
void assembleWithJacobianForPressureEquations(
const double t, double const dt, Eigen::VectorXd const& local_x,
Eigen::VectorXd const& local_xdot, std::vector<double>& local_b_data,
std::vector<double>& local_Jac_data);
unsigned getNumberOfIntegrationPoints() const override;
typename MaterialLib::Solids::MechanicsBase<
DisplacementDim>::MaterialStateVariables const&
getMaterialStateVariablesAt(unsigned integration_point) const override;
private:
HydroMechanicsProcessData<DisplacementDim>& _process_data;
using BMatricesType =
BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
using IpData =
IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement,
ShapeMatricesTypePressure, DisplacementDim,
ShapeFunctionDisplacement::NPOINTS>;
std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
IntegrationMethod _integration_method;
MeshLib::Element const& _element;
bool const _is_axially_symmetric;
SecondaryData<
typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
_secondary_data;
static const int pressure_index = 0;
static const int pressure_size = ShapeFunctionPressure::NPOINTS;
static const int displacement_index = ShapeFunctionPressure::NPOINTS;
static const int displacement_size =
ShapeFunctionDisplacement::NPOINTS * DisplacementDim;
};
} // namespace HydroMechanics
} // namespace ProcessLib
#include "HydroMechanicsFEM-impl.h"
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...