https://gitlab.opengeosys.org/ogs/ogs.git
Tip revision: 691eea29eab979ef0bdfbaa2814b2ac32122bc53 authored by wenqing on 09 July 2021, 19:21:52 UTC
Merge branch 'fixing' into 'master'
Merge branch 'fixing' into 'master'
Tip revision: 691eea2
ThermoHydroMechanicsFEM.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 "IntegrationPointData.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"
#include "ThermoHydroMechanicsProcessData.h"
namespace ProcessLib
{
namespace ThermoHydroMechanics
{
namespace MPL = MaterialPropertyLib;
/// Used by for 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 ThermoHydroMechanicsLocalAssembler : public LocalAssemblerInterface
{
public:
using ShapeMatricesTypeDisplacement =
ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
// Types for pressure.
using ShapeMatricesTypePressure =
ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>;
using GlobalDimMatrixType =
typename ShapeMatricesTypePressure::GlobalDimMatrixType;
using GlobalDimVectorType =
typename ShapeMatricesTypePressure::GlobalDimVectorType;
static int const KelvinVectorSize =
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
ThermoHydroMechanicsLocalAssembler(
ThermoHydroMechanicsLocalAssembler const&) = delete;
ThermoHydroMechanicsLocalAssembler(ThermoHydroMechanicsLocalAssembler&&) =
delete;
ThermoHydroMechanicsLocalAssembler(
MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
bool const is_axially_symmetric,
unsigned const integration_order,
ThermoHydroMechanicsProcessData<DisplacementDim>& process_data);
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(
"ThermoHydroMechanicsLocalAssembler: 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 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_x,
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;
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::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;
private:
std::size_t setSigma(double const* values)
{
return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
values, _ip_data, &IpData::sigma_eff);
}
// 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
{
return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
_ip_data, &IpData::sigma_eff);
}
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);
}
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 override
{
return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
_ip_data, &IpData::eps, cache);
}
private:
ThermoHydroMechanicsProcessData<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;
// The shape function of pressure has the same form with the shape function
// of temperature
static const int temperature_index = 0;
static const int temperature_size = ShapeFunctionPressure::NPOINTS;
static const int pressure_index = ShapeFunctionPressure::NPOINTS;
static const int pressure_size = ShapeFunctionPressure::NPOINTS;
static const int displacement_index = ShapeFunctionPressure::NPOINTS * 2;
static const int displacement_size =
ShapeFunctionDisplacement::NPOINTS * DisplacementDim;
};
} // namespace ThermoHydroMechanics
} // namespace ProcessLib
#include "ThermoHydroMechanicsFEM-impl.h"