swh:1:snp:f521c49ab17ef7db6ec70b2430e1ed203f50383f
Tip revision: f8eb3fa67df5b9b5c7b36af69a10c90258b3eb2b authored by Dmitri Naumov on 09 March 2021, 11:17:46 UTC
[T/RM] Relax tolerances for failing CI job.
[T/RM] Relax tolerances for failing CI job.
Tip revision: f8eb3fa
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;
void pushBackState()
{
eps_prev = eps;
sigma_eff_prev = sigma_eff;
material_state_variables->pushBackState();
}
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::KelvinVectorDimensions<DisplacementDim>::value;
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++)
{
_ip_data[ip].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());
}
std::size_t setSigma(double const* values);
// 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::size_t setEpsilon(double const* values);
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"