13#ifndef __ADOLCPLASTICITY_HPP_
14#define __ADOLCPLASTICITY_HPP_
17#error "MoFEM need to be compiled with ADOL-C"
32 return FTensor::Dg<double, 3, 6>{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
33 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0,
34 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
35 0.0, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
42 return FTensor::Dg<double, 3, 6>{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
43 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
44 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
45 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
52 return FTensor::Dg<double, 3, 6>{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
53 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
54 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
55 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
61struct CommonData : boost::enable_shared_from_this<CommonData> {
106 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
gradAtGaussPts);
117 PetscBool b_bar =
bBar ? PETSC_TRUE : PETSC_FALSE;
118 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, PETSC_NULLPTR,
"-b_bar", &b_bar,
125 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
stressMatrix);
128 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
168 return getVectorAdaptor(&(
gradientW[6]), 6);
188 ublas::symmetric_matrix<double, ublas::lower>
C,
D;
208 MoFEMErrorCode
playW();
212 MoFEMErrorCode
playY();
217 MoFEMErrorCode
playH();
268 virtual MoFEMErrorCode
270 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
271 std::string block_name, Sev sev) {
283 bool &recalculate_elastic_tangent) {
329 SmartPetscObj<Mat>
A;
330 SmartPetscObj<Vec>
R;
333 ublas::matrix<double, ublas::column_major>
dataA;
359template <
int DIM, StrainType STRAIN>
362 boost::shared_ptr<ClosestPointProjection> cp_ptr,
bool calc_lhs);
367 boost::shared_ptr<ClosestPointProjection> cp_ptr,
bool calc_lhs);
372 boost::shared_ptr<ClosestPointProjection> cp_ptr,
bool calc_lhs);
377 boost::shared_ptr<ClosestPointProjection> cp_ptr,
bool calc_lhs);
382 boost::shared_ptr<ClosestPointProjection> cp_ptr,
bool calc_lhs);
392template <
int DIM, IntegrationType I,
typename AssemblyDomainEleOp>
403template <
int DIM, IntegrationType I,
typename AssemblyDomainEleOp>
413 template <
int DIM, IntegrationType I>
416 template <
int DIM, IntegrationType I>
421using Pip = boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator>;
439template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
442 std::string block_name, boost::shared_ptr<CommonData> common_data_ptr,
443 boost::shared_ptr<ClosestPointProjection> cp_ptr, Sev sev = Sev::inform) {
446 CHKERR cp_ptr->addMatBlockOps(m_field, pip, block_name, sev);
449 pip.push_back(
new typename P::template Assembly<A>::template
OpRhs<DIM, I>(
469template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
472 std::string block_name,
473 boost::shared_ptr<CommonData> common_data_ptr,
474 boost::shared_ptr<ClosestPointProjection> cp_ptr) {
477 CHKERR cp_ptr->addMatBlockOps(m_field, pip, block_name, Sev::noisy);
480 pip.push_back(
new typename P::template Assembly<A>::template
OpLhs<DIM, I>(
500 std::string block_name,
501 boost::shared_ptr<CommonData> common_data_ptr,
502 boost::shared_ptr<ClosestPointProjection> cp_ptr);
510 std::string block_name,
511 boost::shared_ptr<CommonData> common_data_ptr,
512 boost::shared_ptr<ClosestPointProjection> cp_ptr);
520 std::string block_name,
521 boost::shared_ptr<CommonData> common_data_ptr,
522 boost::shared_ptr<ClosestPointProjection> cp_ptr);
534 boost::shared_ptr<FEMethod> fe_ptr);
543 MatrixDouble &storage,
const bool b_bar,
544 const int nb_integration_pts,
double *w_ptr,
549 MatrixDouble &storage,
const bool b_bar,
550 const int nb_integration_pts,
double *w_ptr,
555template <
int DIM,
typename AssemblyDomainEleOp>
558 boost::shared_ptr<CommonData> common_data_ptr);
559 MoFEMErrorCode
iNtegrate(EntitiesFieldData::EntData &data);
566template <
int DIM,
typename AssemblyDomainEleOp>
569 boost::shared_ptr<CommonData> common_data_ptr);
570 MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data,
571 DataForcesAndSourcesCore::EntData &col_data);
579template <
int DIM,
typename AssemblyDomainEleOp>
581 string field_name, boost::shared_ptr<CommonData> common_data_ptr)
583 commonDataPtr(common_data_ptr) {}
586template <
int DIM,
typename AssemblyDomainEleOp>
588 EntitiesFieldData::EntData &data) {
593 double *w_ptr = &(OP::getGaussPts()(DIM, 0));
596 data, baseStorage, commonDataPtr->bBar, OP::getGaussPts().size2(), w_ptr,
599 auto t_stress = getFTensor2SymmetricFromMat<3>(commonDataPtr->stressMatrix);
601 const auto vol = OP::getMeasure();
602 auto t_w = OP::getFTensor0IntegrationWeight();
603 for (
int gg = 0; gg != OP::nbIntegrationPts; gg++) {
604 double alpha = vol * t_w;
605 for (
int bb = 0; bb != OP::nbRows; ++bb) {
606 OP::locF[bb] += alpha * (t_stress(
i,
j) * t_diff(
i,
j));
616template <
int DIM,
typename AssemblyDomainEleOp>
618 string field_name, boost::shared_ptr<CommonData> common_data_ptr)
621 commonDataPtr(common_data_ptr) {
625template <
int DIM,
typename AssemblyDomainEleOp>
627 DataForcesAndSourcesCore::EntData &row_data,
628 DataForcesAndSourcesCore::EntData &col_data) {
632 double *w_ptr = &(OP::getGaussPts()(DIM, 0));
634 row_data, baseRowStorage, commonDataPtr->bBar, OP::getGaussPts().size2(),
637 col_data, baseColStorage, commonDataPtr->bBar, OP::getGaussPts().size2(),
640 auto get_ftensor2_symmetric = [&](
auto &storage,
const auto gg,
643 &storage(gg, 6 * rr + 0), &storage(gg, 6 * rr + 1),
644 &storage(gg, 6 * rr + 2), &storage(gg, 6 * rr + 3),
645 &storage(gg, 6 * rr + 4), &storage(gg, 6 * rr + 5)};
654 getFTensor4DdgFromMat<3, 3, 1>(commonDataPtr->materialTangentOperator);
655 const auto vol = OP::getMeasure();
656 auto t_w = OP::getFTensor0IntegrationWeight();
657 for (
int gg = 0; gg != OP::nbIntegrationPts; ++gg) {
658 const double alpha = vol * t_w;
661 for (
auto rr = 0; rr != OP::nbRows; ++rr) {
663 t_rowCep(
k,
l) = t_Cep(
i,
j,
k,
l) * t_diff_row(
i,
j);
664 auto t_diff_col = get_ftensor2_symmetric(baseColStorage, gg, 0);
665 for (
auto cc = 0; cc != OP::nbCols; ++cc) {
666 OP::locMat(rr, cc) += alpha * (t_rowCep(
k,
l) * t_diff_col(
k,
l));
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, string field_name, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr)
Assemble the left hand side, i.e. tangent matrix.
MoFEMErrorCode opFactoryDomainUpdate(MoFEM::Interface &m_field, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr)
Push operators to update history variables.
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, string field_name, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, Sev sev=Sev::inform)
Assemble the left hand side, i.e. tangent matrix.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress< 3, SMALL_STRAIN >(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress< 2, SMALL_STRAIN >(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
FTensor::Dg< double, 3, 6 > strain_to_voight_op()
Op convert strain tensor to Vight strain vector.
FTensor::Dg< double, 3, 6 > voight_to_stress_op()
Op convert Vight stress vector to stress tensor.
MoFEMErrorCode opFactoryDomainUpdate< 2 >(MoFEM::Interface &m_field, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr)
Push operators to update history variables.
MoFEMErrorCode ADOLCPlasticityJac(SNES, Vec, Mat, Mat, void *ctx)
Internal SNES function used at integration points to calulate tangent matrix.
FTensor::Dg< double, 3, 6 > voight_to_strain_op()
Op convert Vight strain vector to strain tensor.
MoFEMErrorCode ADOLCPlasticityRes(SNES snes, Vec chi, Vec r, void *ctx)
Internal SNES function used at integration points to calulate stress.
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
Get opreator to calulate stress.
boost::shared_ptr< TSUpdate > createTSUpdate(std::string fe_name, boost::shared_ptr< FEMethod > fe_ptr)
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress< 2, LARGE_STRAIN >(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress< 3, LARGE_STRAIN >(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > Pip
MoFEMErrorCode opFactoryDomainUpdate< 3 >(MoFEM::Interface &m_field, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr)
Push operators to update history variables.
constexpr auto field_name
Closest point projection algorithm.
VectorAdaptor activeVariablesW
MatrixDouble partial2HSigmaFlux
VectorAdaptor getInternalFluxes()
MoFEMErrorCode playY_NoGradient()
ublas::vector< adouble > a_internalFluxes
MoFEMErrorCode recordTapes()
Record tapes.
ublas::vector< adouble > a_plasticStrain
std::array< int, LAST_TAPE > tapesTags
MoFEMErrorCode solveClosestProjection()
Solve nonlinear system of equations to find stress, internal fluxes, and Lagrange plastic multiplier.
VectorAdaptor partialYFlux
MatrixDouble partialWStrainPlasticStrain
MoFEMErrorCode recordH()
Record flow potential.
MoFEMErrorCode evaluatePotentials()
friend MoFEMErrorCode ADOLCPlasticityRes(SNES, Vec, Vec, void *ctx)
Function executed by nested SENES at evaluationg residual.
MoFEMErrorCode recordY()
Record yield function.
boost::function< int(int, int, int)> integrationRule
virtual MoFEMErrorCode flowPotential()=0
Set flow potential.
VectorAdaptor getActiveVariablesYH()
ublas::symmetric_matrix< double, ublas::lower > partial2HFlux
MoFEMErrorCode recordW()
Record strain energy.
virtual MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, Sev sev)
Get model parameters from blocks.
VectorAdaptor partialHFlux
ublas::symmetric_matrix< double, ublas::lower > D
MoFEMErrorCode calculateA()
VectorDouble internalVariables0
VectorAdaptor partialHSigma
MoFEMErrorCode playH_NoHessian()
ublas::vector< adouble > a_internalVariables
VectorAdaptor getTotalStrain()
ublas::symmetric_matrix< double, ublas::lower > C
virtual MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent)
Set parameters for ADOL-C of constitutive relations.
virtual MoFEMErrorCode codedHessianW(vector< double * >)
VectorAdaptor partialYSigma
MoFEMErrorCode playPotentials_NoHessian()
SmartPetscObj< Vec > dChi
SmartPetscObj< SNES > sNes
virtual MoFEMErrorCode yieldFunction()=0
Set yield function.
VectorAdaptor getInternalVariables()
MoFEMErrorCode calculateR(Vec R)
MoFEMErrorCode createMatAVecR()
ublas::vector< adouble > a_sTress
MoFEMErrorCode consistentTangent()
Calculate consistent tangent matrix.
MoFEMErrorCode playW_NoHessian()
ublas::symmetric_matrix< double, ublas::lower > partial2HSigma
VectorAdaptor getPlasticStrain()
ublas::vector< adouble > a_sTrain
friend MoFEMErrorCode ADOLCPlasticityJac(SNES, Vec, Mat, Mat, void *ctx)
Function executed by nested SENES at evaluationg tangent matrix.
ublas::matrix< double, ublas::column_major > dataA
VectorDouble plasticStrain0
virtual MoFEMErrorCode freeHemholtzEnergy()=0
Set Hemholtz energy.
MoFEMErrorCode playPotentials()
VectorAdaptor getStress()
MoFEMErrorCode snesCreate()
Create nested snes.
common data used by volume elements
MatrixDouble materialTangentOperator
MoFEMErrorCode getDefaultMaterialParameters()
double * plasticStrainPtr
auto getPlasticStrain(int gg=0)
MatrixDouble activeVariablesW
MatrixDouble plasticStrainMatrix
MatrixDouble stressMatrix
MatrixDouble gradAtGaussPts
auto getInternalVariables(int gg=0)
auto getFTensor1StressAtGaussPts(int gg=0)
auto getFTensor1PlasticStrainAtGaussPts(int gg=0)
boost::shared_ptr< MatrixDouble > getStressMatrixPtr()
double * internalVariablesPtr
auto getGradAtGaussPtsPtr()
int internalVariablesSize
vector< double > deltaGamma
boost::shared_ptr< MatrixDouble > getPlasticStrainMatrixPtr()
Calculate tensorial base functions. Apply bBar method when needed.
static FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 6 >, 3 > getFTensor2SymmetricDiffBase(DataForcesAndSourcesCore::EntData &data, MatrixDouble &storage, const bool b_bar, const int nb_integration_pts, double *w_ptr, FTensor::Number< 2 >)
MatrixDouble baseRowStorage
Store tensorial base functions.
OpLhsImpl(std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
MatrixDouble baseColStorage
Store tensorial base functions.
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
MatrixDouble baseStorage
Store tensorial base functions.
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Assemble right hand side.
Update internal fluxes (update history variables)
virtual MoFEMErrorCode postProcess(TS ts)=0
Deprecated interface functions.
FormsIntegrators< DomainEleOp >::Assembly< A >::OpBase AssemblyDomainEleOp