const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
auto t_stress =
getFTensor2SymmetricFromMat<SPACE_DIM>(*(
commonDataPtr->mStressPtr));
auto t_flow =
getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->plasticFlow);
t_flow(
i,
j) = t_flow_tmp(
i,
j);
++t_flow;
++t_stress;
}
}
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
auto t_tau_dot = getFTensor0FromVec(
commonDataPtr->plasticTauDot);
auto t_flow =
getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->plasticFlow);
auto t_plastic_strain_dot =
getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->plasticStrainDot);
auto t_stress =
getFTensor2SymmetricFromMat<SPACE_DIM>(*(
commonDataPtr->mStressPtr));
auto t_D =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
commonDataPtr->mDPtr);
auto t_D_Op = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
mDPtr);
t_flow_dir_dstress(
i,
j,
k,
l) =
1.5 * (t_diff_deviator(
M,
N,
i,
j) * t_diff_deviator(
M,
N,
k,
l));
t_flow_dir_dstrain(
i,
j,
k,
l) =
t_flow_dir_dstress(
i,
j,
m,
n) * t_D_Op(
m,
n,
k,
l);
false);
false);
auto t_res_c_dstrain =
getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resCdStrain);
auto t_res_c_dstrain_dot =
getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resCdStrainDot);
auto t_res_flow =
auto t_res_flow_dtau =
getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resFlowDtau);
auto t_res_flow_dstrain = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(
auto t_res_flow_dstrain_dot = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(
auto next = [&]() {
++t_tau;
++t_tau_dot;
++t_f;
++t_flow;
++t_plastic_strain_dot;
++t_stress;
++t_res_c;
++t_res_c_dtau;
++t_res_c_dstrain;
++t_res_c_dstrain_dot;
++t_res_flow;
++t_res_flow_dtau;
++t_res_flow_dstrain;
++t_res_flow_dstrain_dot;
++t_w;
};
auto get_avtive_pts = [&]() {
int nb_points_avtive_on_elem = 0;
int nb_points_on_elem = 0;
auto t_tau_dot = getFTensor0FromVec(
commonDataPtr->plasticTauDot);
auto t_plastic_strain_dot =
getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->plasticStrainDot);
const auto ww =
w(eqiv, t_tau_dot, t_f,
hardening(t_tau));
++nb_points_on_elem;
if (sign_ww > 0) {
++nb_points_avtive_on_elem;
}
++t_tau;
++t_tau_dot;
++t_f;
++t_plastic_strain_dot;
}
++nb_elements;
nb_points += nb_points_on_elem;
if (nb_points_avtive_on_elem > 0) {
++avtive_elems;
active_points += nb_points_avtive_on_elem;
if (nb_points_avtive_on_elem == nb_points_on_elem) {
++avtive_full_elems;
}
}
if (nb_points_avtive_on_elem != nb_points_on_elem)
return 1;
else
return 0;
};
if (
getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
get_avtive_pts();
}
t_diff_plastic_strain);
auto ww =
w(eqiv, t_tau_dot, t_f, sigma_y);
auto c =
constraint(eqiv, t_tau_dot, t_f, sigma_y, abs_ww);
t_flow_dir(
k,
l) = 1.5 * (t_dev_stress(
I,
J) * t_diff_deviator(
I,
J,
k,
l));
t_flow_dstrain(
i,
j) = t_flow(
k,
l) * t_D_Op(
k,
l,
i,
j);
auto get_res_c = [&]() {
return c; };
auto get_res_c_dstrain = [&](auto &t_diff_res) {
t_diff_res(
i,
j) = c_f * t_flow_dstrain(
i,
j);
};
auto get_res_c_dstrain_dot = [&](auto &t_diff_res) {
t_diff_res(
i,
j) = (
getTSa() * c_equiv) * t_diff_eqiv(
i,
j);
};
auto get_res_c_dtau = [&]() {
return getTSa() * c_dot_tau + c_sigma_y * d_sigma_y;
};
auto get_res_flow = [&](auto &t_res_flow) {
const auto b = t_tau_dot;
t_res_flow(
k,
l) =
a * t_plastic_strain_dot(
k,
l) - b * t_flow_dir(
k,
l);
};
auto get_res_flow_dtau = [&](auto &t_res_flow_dtau) {
const auto da = d_sigma_y;
da * t_plastic_strain_dot(
k,
l) - db * t_flow_dir(
k,
l);
};
auto get_res_flow_dstrain = [&](auto &t_res_flow_dstrain) {
const auto b = t_tau_dot;
t_res_flow_dstrain(
m,
n,
k,
l) = -t_flow_dir_dstrain(
m,
n,
k,
l) * b;
};
auto get_res_flow_dstrain_dot = [&](auto &t_res_flow_dstrain_dot) {
t_res_flow_dstrain_dot(
m,
n,
k,
l) =
};
t_res_c = get_res_c();
get_res_flow(t_res_flow);
if (
getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
t_res_c_dtau = get_res_c_dtau();
get_res_c_dstrain(t_res_c_dstrain);
get_res_c_dstrain_dot(t_res_c_dstrain_dot);
get_res_flow_dtau(t_res_flow_dtau);
get_res_flow_dstrain(t_res_flow_dstrain);
get_res_flow_dstrain_dot(t_res_flow_dstrain_dot);
}
next();
}
}
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr,
commonDataPtr(common_data_ptr), scaleStress(
scale), mDPtr(m_D_ptr) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
nb_gauss_pts);
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
mDPtr);
auto t_strain =
getFTensor2SymmetricFromMat<SPACE_DIM>(*(
commonDataPtr->mStrainPtr));
auto t_plastic_strain =
getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->plasticStrain);
auto t_stress =
getFTensor2SymmetricFromMat<SPACE_DIM>(*(
commonDataPtr->mStressPtr));
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
t_D(
i,
j,
k,
l) * (t_strain(
k,
l) - t_plastic_strain(
k,
l));
++t_strain;
++t_plastic_strain;
++t_stress;
}
}
const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
MoFEMErrorCode
const auto nb_integration_pts = getGaussPts().size2();
const auto nb_base_functions = data.getN().size2();
auto t_res_flow =
auto next = [&]() { ++t_res_flow; };
auto t_w = getFTensor0IntegrationWeight();
auto t_base = data.getFTensor0N();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
double alpha = getMeasure() * t_w;
++t_w;
t_rhs(
L) = alpha * (t_res_flow(
i,
j) * t_L(
i,
j,
L));
next();
auto t_nf = getFTensor1FromArray<size_symm, size_symm>(nf);
size_t bb = 0;
t_nf(
L) += t_base * t_rhs(
L);
++t_base;
++t_nf;
}
for (; bb < nb_base_functions; ++bb)
++t_base;
}
}
const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
MoFEMErrorCode
const size_t nb_integration_pts = getGaussPts().size2();
const size_t nb_base_functions = data.getN().size2();
auto next = [&]() { ++t_res_c; };
auto t_w = getFTensor0IntegrationWeight();
auto t_base = data.getFTensor0N();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = getMeasure() * t_w;
++t_w;
const auto res = alpha * t_res_c;
next();
size_t bb = 0;
nf[bb] += t_base * res;
++t_base;
}
for (; bb < nb_base_functions; ++bb)
++t_base;
}
}
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
sYmm = false;
}
&mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 0, 2),
&mat(3 * rr + 1, 0), &mat(3 * rr + 1, 1), &mat(3 * rr + 1, 2),
&mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1), &mat(3 * rr + 2, 2)};
}
&mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
&mat(6 * rr + 0, 3), &mat(6 * rr + 0, 4), &mat(6 * rr + 0, 5),
&mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
&mat(6 * rr + 1, 3), &mat(6 * rr + 1, 4), &mat(6 * rr + 1, 5),
&mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
&mat(6 * rr + 2, 3), &mat(6 * rr + 2, 4), &mat(6 * rr + 2, 5),
&mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
&mat(6 * rr + 3, 3), &mat(6 * rr + 3, 4), &mat(6 * rr + 3, 5),
&mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
&mat(6 * rr + 4, 3), &mat(6 * rr + 4, 4), &mat(6 * rr + 4, 5),
&mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2),
&mat(6 * rr + 5, 3), &mat(6 * rr + 5, 4), &mat(6 * rr + 5, 5)};
}
MoFEMErrorCode
EntitiesFieldData::EntData &col_data) {
const auto nb_integration_pts = getGaussPts().size2();
const auto nb_row_base_functions = row_data.getN().size2();
auto t_res_flow_dstrain = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(
auto t_res_flow_dstrain_dot = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(
auto next = [&]() {
++t_res_flow_dstrain;
++t_res_flow_dstrain_dot;
};
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
double alpha = getMeasure() * t_w;
++t_w;
alpha * (t_L(
i,
j,
O) * ((t_res_flow_dstrain_dot(
i,
j,
k,
l) -
t_res_flow_dstrain(
i,
j,
k,
l)) *
next();
size_t rr = 0;
auto t_col_base = col_data.getFTensor0N(gg, 0);
t_mat(
O,
L) += ((t_row_base * t_col_base) * t_res_mat(
O,
L));
++t_mat;
++t_col_base;
}
++t_row_base;
}
for (; rr < nb_row_base_functions; ++rr)
++t_row_base;
}
}
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
sYmm = false;
}
&mat(3 * rr + 0, 0), &mat(3 * rr + 1, 0), &mat(3 * rr + 2, 0)};
}
&mat(6 * rr + 0, 0), &mat(6 * rr + 1, 0), &mat(6 * rr + 2, 0),
&mat(6 * rr + 3, 0), &mat(6 * rr + 4, 0), &mat(6 * rr + 5, 0)};
}
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const auto nb_integration_pts = getGaussPts().size2();
const size_t nb_row_base_functions = row_data.getN().size2();
const auto type = getFEType();
const auto nb_nodes = moab::CN::VerticesPerEntity(type);
auto t_res_flow_dtau =
getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resFlowDtau);
auto next = [&]() { ++t_res_flow_dtau; };
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
double alpha = getMeasure() * t_w;
++t_w;
t_res_vec(
L) = alpha * (t_res_flow_dtau(
i,
j) * t_L(
i,
j,
L));
next();
size_t rr = 0;
auto t_mat =
auto t_col_base = col_data.getFTensor0N(gg, 0);
t_mat(
L) += t_row_base * t_col_base * t_res_vec(
L);
++t_mat;
++t_col_base;
}
++t_row_base;
}
for (; rr != nb_row_base_functions; ++rr)
++t_row_base;
}
}
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr)
commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
sYmm = false;
}
&mat(0, 0), &mat(0, 1), &mat(0, 2)};
}
&mat(0, 0), &mat(0, 1), &mat(0, 2), &mat(0, 3), &mat(0, 4), &mat(0, 5)};
}
MoFEMErrorCode
EntitiesFieldData::EntData &col_data) {
const auto nb_integration_pts = getGaussPts().size2();
const auto nb_row_base_functions = row_data.getN().size2();
auto t_c_dstrain =
getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resCdStrain);
auto t_c_dstrain_dot =
getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resCdStrainDot);
auto next = [&]() {
++t_c_dstrain;
++t_c_dstrain_dot;
};
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = getMeasure() * t_w;
++t_w;
t_res_vec(
L) = t_L(
i,
j,
L) * (t_c_dstrain_dot(
i,
j) - t_c_dstrain(
i,
j));
next();
auto t_mat =
size_t rr = 0;
const auto row_base = alpha * t_row_base;
auto t_col_base = col_data.getFTensor0N(gg, 0);
t_mat(
L) += (row_base * t_col_base) * t_res_vec(
L);
++t_mat;
++t_col_base;
}
++t_row_base;
}
for (; rr != nb_row_base_functions; ++rr)
++t_row_base;
}
}
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {
sYmm = false;
}
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const auto nb_integration_pts = getGaussPts().size2();
const auto nb_row_base_functions = row_data.getN().size2();
auto next = [&]() { ++t_res_c_dtau; };
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = getMeasure() * t_w;
++t_w;
const auto res = alpha * (t_res_c_dtau);
next();
auto mat_ptr =
locMat.data().begin();
size_t rr = 0;
auto t_col_base = col_data.getFTensor0N(gg, 0);
*mat_ptr += t_row_base * t_col_base * res;
++t_col_base;
++mat_ptr;
}
++t_row_base;
}
for (; rr < nb_row_base_functions; ++rr)
++t_row_base;
}
}
};
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
#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()
const double c
speed of light (cm/ns)
static auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
auto deviator(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_stress, double trace)
static auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
double constrian_sign(double x)
double platsic_surface(FTensor::Tensor2_symmetric< double, 3 > &&t_stress_deviator)
auto get_mat_scalar_dtensor_sym(MatrixDouble &mat, FTensor::Number< 2 >)
FTensor::Index< 'j', SPACE_DIM > j
FTensor::Index< 'M', 3 > M
double diff_constrain_eqiv(double sign, double eqiv, double dot_tau)
auto plastic_flow(long double f, FTensor::Tensor2_symmetric< double, 3 > &&t_dev_stress, FTensor::Ddg< double, 3, SPACE_DIM > &&t_diff_deviator)
FTensor::Index< 'O', size_symm > O
FTensor::Index< 'L', size_symm > L
FTensor::Index< 'l', SPACE_DIM > l
FTensor::Index< 'k', SPACE_DIM > k
auto diff_tensor()
[Operators definitions]
auto diff_constrain_df(double sign)
double w(double eqiv, double dot_tau, double f, double sigma_y)
FTensor::Index< 'N', 3 > N
auto diff_constrain_dsigma_y(double sign)
FTensor::Index< 'I', 3 > I
auto diff_deviator(FTensor::Ddg< double, SPACE_DIM, SPACE_DIM > &&t_diff_stress)
FTensor::Index< 'i', SPACE_DIM > i
[Common data]
FTensor::Index< 'm', SPACE_DIM > m
double constrain_abs(double x)
FTensor::Index< 'J', 3 > J
auto diff_equivalent_strain_dot(const T1 eqiv, T2 &t_plastic_strain_dot, T3 &t_diff_plastic_strain)
double trace(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_stress)
double constraint(double eqiv, double dot_tau, double f, double sigma_y, double abs_w)
double diff_constrain_ddot_tau(double sign, double eqiv)
auto equivalent_strain_dot(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_plastic_strain_dot)
FTensor::Index< 'n', SPACE_DIM > n
FormsIntegrators< DomainEleOp >::Assembly< A >::OpBase AssemblyDomainEleOp
double hardening(double tau)
double hardening_dtau(double tau)
constexpr auto field_name
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
auto getFTensor0IntegrationWeight()
Get integration weights.
const TSMethod::TSContext getTSCtx() const
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
OpCalculateConstraintsLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > mat_D_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
OpCalculateConstraintsLhs_dTAU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
OpCalculateConstraintsRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
boost::shared_ptr< CommonData > commonDataPtr
OpCalculatePlasticFlowLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
OpCalculatePlasticFlowLhs_dTAU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
OpCalculatePlasticFlowRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
[Calculate stress]
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculatePlasticSurface(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
OpCalculatePlasticity(const std::string field_name, MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > mDPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
OpPlasticStress(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > mDPtr, const double scale=1)
boost::shared_ptr< MatrixDouble > mDPtr