https://gitlab.opengeosys.org/ogs/ogs.git
Tip revision: 67fde4f3bd01e44239b00ac37051b9ce91a6780d authored by Tom Fischer on 08 November 2021, 11:17:29 UTC
Fix typo.
Fix typo.
Tip revision: 67fde4f
SmallDeformationFEM.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 <limits>
#include <memory>
#include <vector>
#include "LocalAssemblerInterface.h"
#include "MaterialLib/PhysicalConstant.h"
#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "ParameterLib/Parameter.h"
#include "ProcessLib/Deformation/BMatrixPolicy.h"
#include "ProcessLib/Deformation/GMatrixPolicy.h"
#include "ProcessLib/Deformation/LinearBMatrix.h"
#include "ProcessLib/LocalAssemblerInterface.h"
#include "ProcessLib/LocalAssemblerTraits.h"
#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
#include "SmallDeformationProcessData.h"
namespace ProcessLib
{
namespace SmallDeformation
{
namespace MPL = MaterialPropertyLib;
template <typename BMatricesType, typename ShapeMatricesType,
int DisplacementDim>
struct IntegrationPointData final
{
explicit IntegrationPointData(
MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
solid_material)
: solid_material(solid_material),
material_state_variables(
solid_material.createMaterialStateVariables())
{
}
typename BMatricesType::KelvinVectorType sigma, sigma_prev;
typename BMatricesType::KelvinVectorType eps, eps_prev;
double free_energy_density = 0;
MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
DisplacementDim>::MaterialStateVariables>
material_state_variables;
double integration_weight;
typename ShapeMatricesType::NodalRowVectorType N;
typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
void pushBackState()
{
eps_prev = eps;
sigma_prev = sigma;
material_state_variables->pushBackState();
}
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;
};
template <typename ShapeFunction, typename IntegrationMethod,
int DisplacementDim>
class SmallDeformationLocalAssembler
: public SmallDeformationLocalAssemblerInterface<DisplacementDim>
{
public:
using ShapeMatricesType =
ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>;
using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
using BMatricesType = BMatrixPolicyType<ShapeFunction, DisplacementDim>;
using BMatrixType = typename BMatricesType::BMatrixType;
using StiffnessMatrixType = typename BMatricesType::StiffnessMatrixType;
using NodalForceVectorType = typename BMatricesType::NodalForceVectorType;
using NodalDisplacementVectorType =
typename BMatricesType::NodalForceVectorType;
using GMatricesType = GMatrixPolicyType<ShapeFunction, DisplacementDim>;
using GradientVectorType = typename GMatricesType::GradientVectorType;
using GradientMatrixType = typename GMatricesType::GradientMatrixType;
using IpData =
IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>;
SmallDeformationLocalAssembler(SmallDeformationLocalAssembler const&) =
delete;
SmallDeformationLocalAssembler(SmallDeformationLocalAssembler&&) = delete;
SmallDeformationLocalAssembler(
MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
bool const is_axially_symmetric,
unsigned const integration_order,
SmallDeformationProcessData<DisplacementDim>& process_data)
: _process_data(process_data),
_integration_method(integration_order),
_element(e),
_is_axially_symmetric(is_axially_symmetric)
{
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
_ip_data.reserve(n_integration_points);
_secondary_data.N.resize(n_integration_points);
auto const shape_matrices =
NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
DisplacementDim>(e, is_axially_symmetric,
_integration_method);
auto& solid_material =
MaterialLib::Solids::selectSolidConstitutiveRelation(
_process_data.solid_materials,
_process_data.material_ids,
e.getID());
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
_ip_data.emplace_back(solid_material);
auto& ip_data = _ip_data[ip];
auto const& sm = shape_matrices[ip];
_ip_data[ip].integration_weight =
_integration_method.getWeightedPoint(ip).getWeight() *
sm.integralMeasure * sm.detJ;
ip_data.N = sm.N;
ip_data.dNdx = sm.dNdx;
static const int kelvin_vector_size =
MathLib::KelvinVector::kelvin_vector_dimensions(
DisplacementDim);
// Initialize current time step values
ip_data.sigma.setZero(kelvin_vector_size);
ip_data.eps.setZero(kelvin_vector_size);
// Previous time step values are not initialized and are set later.
ip_data.sigma_prev.resize(kelvin_vector_size);
ip_data.eps_prev.resize(kelvin_vector_size);
_secondary_data.N[ip] = shape_matrices[ip].N;
}
}
/// Returns number of read integration points.
std::size_t setIPDataInitialConditions(std::string const& name,
double const* values,
int const integration_order) override
{
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 setSigma(values);
}
return 0;
}
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<ShapeFunction,
ShapeMatricesType>(
_element, ip_data.N))};
ip_data.sigma =
MathLib::KelvinVector::symmetricTensorToKelvinVector<
DisplacementDim>((*_process_data.initial_stress)(
std::numeric_limits<
double>::quiet_NaN() /* time independent */,
x_position));
}
ip_data.pushBackState();
}
}
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_b_data*/) override
{
OGS_FATAL(
"SmallDeformationLocalAssembler: 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_b_data,
std::vector<double>& local_Jac_data) override
{
auto const local_matrix_size = local_x.size();
auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
local_Jac_data, local_matrix_size, local_matrix_size);
auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
local_b_data, local_matrix_size);
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
MPL::VariableArray variables_prev;
MPL::VariableArray variables;
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
auto const& b = _process_data.specific_body_force;
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
auto const& w = _ip_data[ip].integration_weight;
auto const& N = _ip_data[ip].N;
auto const& dNdx = _ip_data[ip].dNdx;
typename ShapeMatricesType::template MatrixType<DisplacementDim,
displacement_size>
N_u_op = ShapeMatricesType::template MatrixType<
DisplacementDim,
displacement_size>::Zero(DisplacementDim,
displacement_size);
for (int i = 0; i < DisplacementDim; ++i)
{
N_u_op
.template block<1, displacement_size / DisplacementDim>(
i, i * displacement_size / DisplacementDim)
.noalias() = N;
}
auto const x_coord =
NumLib::interpolateXCoordinate<ShapeFunction,
ShapeMatricesType>(_element, N);
auto const B = LinearBMatrix::computeBMatrix<
DisplacementDim, ShapeFunction::NPOINTS,
typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
_is_axially_symmetric);
auto const& eps_prev = _ip_data[ip].eps_prev;
auto const& sigma_prev = _ip_data[ip].sigma_prev;
auto& eps = _ip_data[ip].eps;
auto& sigma = _ip_data[ip].sigma;
auto& state = _ip_data[ip].material_state_variables;
eps.noalias() =
B *
Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
variables_prev[static_cast<int>(MPL::Variable::stress)]
.emplace<
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
sigma_prev);
variables_prev[static_cast<int>(MPL::Variable::mechanical_strain)]
.emplace<
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
eps_prev);
double const T_ref =
_process_data.reference_temperature
? (*_process_data.reference_temperature)(t, x_position)[0]
: std::numeric_limits<double>::quiet_NaN();
variables_prev[static_cast<int>(MPL::Variable::temperature)]
.emplace<double>(T_ref);
variables[static_cast<int>(MPL::Variable::mechanical_strain)]
.emplace<
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
eps);
variables[static_cast<int>(MPL::Variable::temperature)]
.emplace<double>(T_ref);
auto&& solution = _ip_data[ip].solid_material.integrateStress(
variables_prev, variables, t, x_position, dt, *state);
if (!solution)
{
OGS_FATAL("Computation of local constitutive relation failed.");
}
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C;
std::tie(sigma, state, C) = std::move(*solution);
auto const rho = _process_data.solid_density(t, x_position)[0];
local_b.noalias() -=
(B.transpose() * sigma - N_u_op.transpose() * rho * b) * w;
local_Jac.noalias() += B.transpose() * C * B * w;
}
}
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();
}
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
auto& ip_data = _ip_data[ip];
// Update free energy density needed for material forces.
ip_data.free_energy_density =
ip_data.solid_material.computeFreeEnergyDensity(
t, x_position, dt, ip_data.eps, ip_data.sigma,
*ip_data.material_state_variables);
}
}
std::vector<double> const& getMaterialForces(
std::vector<double> const& local_x,
std::vector<double>& nodal_values) override
{
return ProcessLib::SmallDeformation::getMaterialForces<
DisplacementDim, ShapeFunction, ShapeMatricesType,
typename BMatricesType::NodalForceVectorType,
NodalDisplacementVectorType, GradientVectorType,
GradientMatrixType>(local_x, nodal_values, _integration_method,
_ip_data, _element, _is_axially_symmetric);
}
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
auto const& N = _secondary_data.N[integration_point];
// assumes N is stored contiguously in memory
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
}
std::vector<double> const& getIntPtFreeEnergyDensity(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& cache) const override
{
cache.clear();
cache.reserve(_ip_data.size());
transform(
cbegin(_ip_data), cend(_ip_data), back_inserter(cache),
[](auto const& ip_data) { return ip_data.free_energy_density; });
return cache;
}
std::size_t setSigma(double const* values)
{
return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
values, _ip_data, &IpData::sigma);
}
// 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);
}
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, 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);
}
unsigned getNumberOfIntegrationPoints() const override
{
return _integration_method.getNumberOfPoints();
}
typename MaterialLib::Solids::MechanicsBase<
DisplacementDim>::MaterialStateVariables const&
getMaterialStateVariablesAt(unsigned integration_point) const override
{
return *_ip_data[integration_point].material_state_variables;
}
void computeSecondaryVariableConcrete(
double const /*t*/, double const /*dt*/, Eigen::VectorXd const& /*x*/,
Eigen::VectorXd const& /*x_dot*/) override
{
int const elem_id = _element.getID();
ParameterLib::SpatialPosition x_position;
x_position.setElementID(elem_id);
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
auto sigma_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
Eigen::Matrix<double, 3, 3>::Zero());
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
auto const& sigma = _ip_data[ip].sigma;
sigma_sum += sigma;
}
Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
MathLib::KelvinVector::kelvinVectorToTensor(sigma_sum) /
n_integration_points;
Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(
sigma_avg);
Eigen::Map<Eigen::Vector3d>(
&(*_process_data.principal_stress_values)[elem_id * 3], 3) =
e_s.eigenvalues();
auto eigen_vectors = e_s.eigenvectors();
for (auto i = 0; i < 3; i++)
{
Eigen::Map<Eigen::Vector3d>(
&(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
eigen_vectors.col(i);
}
}
private:
SmallDeformationProcessData<DisplacementDim>& _process_data;
std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
IntegrationMethod _integration_method;
MeshLib::Element const& _element;
SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
bool const _is_axially_symmetric;
static const int displacement_size =
ShapeFunction::NPOINTS * DisplacementDim;
};
} // namespace SmallDeformation
} // namespace ProcessLib