10 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr)
12 commonDataPtr(common_data_ptr) {
21 const size_t nb_gauss_pts =
commonDataPtr->mStressPtr->size2();
23 getFTensor2SymmetricFromMat<SPACE_DIM>(*(
commonDataPtr->mStressPtr));
28 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->plasticFlow);
37 t_flow(
i,
j) = t_flow_tmp(
i,
j);
47 boost::shared_ptr<CommonData> common_data_ptr,
48 boost::shared_ptr<MatrixDouble> m_D_ptr)
50 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
62 auto t_tau_dot = getFTensor0FromVec(
commonDataPtr->plasticTauDot);
65 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->plasticFlow);
66 auto t_plastic_strain_dot =
67 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->plasticStrainDot);
69 getFTensor2SymmetricFromMat<SPACE_DIM>(*(
commonDataPtr->mStressPtr));
72 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
commonDataPtr->mDPtr);
73 auto t_D_Op = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
mDPtr);
80 t_flow_dir_dstress(
i,
j,
k,
l) =
81 1.5 * (t_diff_deviator(
M,
N,
i,
j) * t_diff_deviator(
M,
N,
k,
l));
82 t_flow_dir_dstrain(
i,
j,
k,
l) =
83 t_flow_dir_dstress(
i,
j,
m,
n) * t_D_Op(
m,
n,
k,
l);
106 auto t_res_c_dtau = getFTensor0FromVec(
commonDataPtr->resCdTau);
107 auto t_res_c_dstrain =
108 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resCdStrain);
109 auto t_res_c_dstrain_dot =
110 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resCdStrainDot);
112 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resFlow);
113 auto t_res_flow_dtau =
114 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resFlowDtau);
115 auto t_res_flow_dstrain = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(
117 auto t_res_flow_dstrain_dot = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(
125 ++t_plastic_strain_dot;
130 ++t_res_c_dstrain_dot;
133 ++t_res_flow_dstrain;
134 ++t_res_flow_dstrain_dot;
138 auto get_avtive_pts = [&]() {
139 int nb_points_avtive_on_elem = 0;
140 int nb_points_on_elem = 0;
143 auto t_tau_dot = getFTensor0FromVec(
commonDataPtr->plasticTauDot);
144 auto t_f = getFTensor0FromVec(
commonDataPtr->plasticSurface);
145 auto t_plastic_strain_dot =
146 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->plasticStrainDot);
150 const auto ww =
w(eqiv, t_tau_dot, t_f,
hardening(t_tau));
155 ++nb_points_avtive_on_elem;
161 ++t_plastic_strain_dot;
171 nb_points += nb_points_on_elem;
172 if (nb_points_avtive_on_elem > 0) {
174 active_points += nb_points_avtive_on_elem;
175 if (nb_points_avtive_on_elem == nb_points_on_elem) {
180 if (nb_points_avtive_on_elem != nb_points_on_elem)
186 if (
getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
194 t_diff_plastic_strain);
199 auto ww =
w(eqiv, t_tau_dot, t_f, sigma_y);
203 auto c =
constraint(eqiv, t_tau_dot, t_f, sigma_y, abs_ww);
211 t_flow_dir(
k,
l) = 1.5 * (t_dev_stress(
I,
J) * t_diff_deviator(
I,
J,
k,
l));
213 t_flow_dstrain(
i,
j) = t_flow(
k,
l) * t_D_Op(
k,
l,
i,
j);
215 auto get_res_c = [&]() {
return c; };
217 auto get_res_c_dstrain = [&](
auto &t_diff_res) {
218 t_diff_res(
i,
j) = c_f * t_flow_dstrain(
i,
j);
221 auto get_res_c_dstrain_dot = [&](
auto &t_diff_res) {
222 t_diff_res(
i,
j) = (
getTSa() * c_equiv) * t_diff_eqiv(
i,
j);
225 auto get_res_c_dtau = [&]() {
226 return getTSa() * c_dot_tau + c_sigma_y * d_sigma_y;
229 auto get_res_flow = [&](
auto &t_res_flow) {
230 const auto a = sigma_y;
231 const auto b = t_tau_dot;
232 t_res_flow(
k,
l) =
a * t_plastic_strain_dot(
k,
l) - b * t_flow_dir(
k,
l);
235 auto get_res_flow_dtau = [&](
auto &t_res_flow_dtau) {
236 const auto da = d_sigma_y;
238 t_res_flow_dtau(
k,
l) =
239 da * t_plastic_strain_dot(
k,
l) - db * t_flow_dir(
k,
l);
242 auto get_res_flow_dstrain = [&](
auto &t_res_flow_dstrain) {
243 const auto b = t_tau_dot;
244 t_res_flow_dstrain(
m,
n,
k,
l) = -t_flow_dir_dstrain(
m,
n,
k,
l) * b;
247 auto get_res_flow_dstrain_dot = [&](
auto &t_res_flow_dstrain_dot) {
248 const auto a = sigma_y;
249 t_res_flow_dstrain_dot(
m,
n,
k,
l) =
250 (
a *
getTSa()) * t_diff_plastic_strain(
m,
n,
k,
l);
253 t_res_c = get_res_c();
254 get_res_flow(t_res_flow);
256 if (
getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
257 t_res_c_dtau = get_res_c_dtau();
258 get_res_c_dstrain(t_res_c_dstrain);
259 get_res_c_dstrain_dot(t_res_c_dstrain_dot);
260 get_res_flow_dtau(t_res_flow_dtau);
261 get_res_flow_dstrain(t_res_flow_dstrain);
262 get_res_flow_dstrain_dot(t_res_flow_dstrain_dot);
272 boost::shared_ptr<CommonData> common_data_ptr,
273 boost::shared_ptr<MatrixDouble> m_D_ptr,
276 commonDataPtr(common_data_ptr), scaleStress(
scale), mDPtr(m_D_ptr) {
285 const size_t nb_gauss_pts =
commonDataPtr->mStrainPtr->size2();
288 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
mDPtr);
290 getFTensor2SymmetricFromMat<SPACE_DIM>(*(
commonDataPtr->mStrainPtr));
291 auto t_plastic_strain =
292 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->plasticStrain);
294 getFTensor2SymmetricFromMat<SPACE_DIM>(*(
commonDataPtr->mStressPtr));
296 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
298 t_D(
i,
j,
k,
l) * (t_strain(
k,
l) - t_plastic_strain(
k,
l));
310 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
311 boost::shared_ptr<MatrixDouble> m_D_ptr)
313 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
319 const auto nb_integration_pts = getGaussPts().size2();
320 const auto nb_base_functions = data.getN().size2();
323 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resFlow);
327 auto next = [&]() { ++t_res_flow; };
329 auto t_w = getFTensor0IntegrationWeight();
330 auto t_base = data.getFTensor0N();
332 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
333 double alpha = getMeasure() * t_w;
337 t_rhs(
L) = alpha * (t_res_flow(
i,
j) * t_L(
i,
j,
L));
340 auto t_nf = getFTensor1FromArray<size_symm, size_symm>(nf);
343 t_nf(
L) += t_base * t_rhs(
L);
347 for (; bb < nb_base_functions; ++bb)
355 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
356 boost::shared_ptr<MatrixDouble> m_D_ptr)
358 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
364 const size_t nb_integration_pts = getGaussPts().size2();
365 const size_t nb_base_functions = data.getN().size2();
369 auto next = [&]() { ++t_res_c; };
371 auto t_w = getFTensor0IntegrationWeight();
373 auto t_base = data.getFTensor0N();
374 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
375 const double alpha = getMeasure() * t_w;
377 const auto res = alpha * t_res_c;
382 nf[bb] += t_base * res;
385 for (; bb < nb_base_functions; ++bb)
393 const std::string row_field_name,
const std::string col_field_name,
394 boost::shared_ptr<CommonData> common_data_ptr,
395 boost::shared_ptr<MatrixDouble> m_D_ptr)
398 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
405 &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 0, 2),
406 &mat(3 * rr + 1, 0), &mat(3 * rr + 1, 1), &mat(3 * rr + 1, 2),
407 &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1), &mat(3 * rr + 2, 2)};
413 &mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
414 &mat(6 * rr + 0, 3), &mat(6 * rr + 0, 4), &mat(6 * rr + 0, 5),
415 &mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
416 &mat(6 * rr + 1, 3), &mat(6 * rr + 1, 4), &mat(6 * rr + 1, 5),
417 &mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
418 &mat(6 * rr + 2, 3), &mat(6 * rr + 2, 4), &mat(6 * rr + 2, 5),
419 &mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
420 &mat(6 * rr + 3, 3), &mat(6 * rr + 3, 4), &mat(6 * rr + 3, 5),
421 &mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
422 &mat(6 * rr + 4, 3), &mat(6 * rr + 4, 4), &mat(6 * rr + 4, 5),
423 &mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2),
424 &mat(6 * rr + 5, 3), &mat(6 * rr + 5, 4), &mat(6 * rr + 5, 5)};
429 EntitiesFieldData::EntData &col_data) {
434 const auto nb_integration_pts = getGaussPts().size2();
435 const auto nb_row_base_functions = row_data.getN().size2();
437 auto t_res_flow_dstrain = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(
439 auto t_res_flow_dstrain_dot = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(
444 ++t_res_flow_dstrain;
445 ++t_res_flow_dstrain_dot;
448 auto t_w = getFTensor0IntegrationWeight();
449 auto t_row_base = row_data.getFTensor0N();
450 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
451 double alpha = getMeasure() * t_w;
456 alpha * (t_L(
i,
j,
O) * ((t_res_flow_dstrain_dot(
i,
j,
k,
l) -
457 t_res_flow_dstrain(
i,
j,
k,
l)) *
465 auto t_col_base = col_data.getFTensor0N(gg, 0);
467 t_mat(
O,
L) += ((t_row_base * t_col_base) * t_res_mat(
O,
L));
475 for (; rr < nb_row_base_functions; ++rr)
483 const std::string row_field_name,
const std::string col_field_name,
484 boost::shared_ptr<CommonData> common_data_ptr,
485 boost::shared_ptr<MatrixDouble> m_D_ptr)
488 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
495 &mat(3 * rr + 0, 0), &mat(3 * rr + 1, 0), &mat(3 * rr + 2, 0)};
501 &mat(6 * rr + 0, 0), &mat(6 * rr + 1, 0), &mat(6 * rr + 2, 0),
502 &mat(6 * rr + 3, 0), &mat(6 * rr + 4, 0), &mat(6 * rr + 5, 0)};
506 EntitiesFieldData::EntData &row_data,
507 EntitiesFieldData::EntData &col_data) {
510 const auto nb_integration_pts = getGaussPts().size2();
511 const size_t nb_row_base_functions = row_data.getN().size2();
514 const auto type = getFEType();
515 const auto nb_nodes = moab::CN::VerticesPerEntity(type);
517 auto t_res_flow_dtau =
518 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resFlowDtau);
522 auto next = [&]() { ++t_res_flow_dtau; };
524 auto t_w = getFTensor0IntegrationWeight();
525 auto t_row_base = row_data.getFTensor0N();
526 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
527 double alpha = getMeasure() * t_w;
530 t_res_vec(
L) = alpha * (t_res_flow_dtau(
i,
j) * t_L(
i,
j,
L));
537 auto t_col_base = col_data.getFTensor0N(gg, 0);
539 t_mat(
L) += t_row_base * t_col_base * t_res_vec(
L);
545 for (; rr != nb_row_base_functions; ++rr)
553 const std::string row_field_name,
const std::string col_field_name,
554 boost::shared_ptr<CommonData> common_data_ptr,
555 boost::shared_ptr<MatrixDouble> m_D_ptr)
558 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
564 &mat(0, 0), &mat(0, 1), &mat(0, 2)};
569 &mat(0, 0), &mat(0, 1), &mat(0, 2), &mat(0, 3), &mat(0, 4), &mat(0, 5)};
574 EntitiesFieldData::EntData &col_data) {
577 const auto nb_integration_pts = getGaussPts().size2();
578 const auto nb_row_base_functions = row_data.getN().size2();
581 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resCdStrain);
582 auto t_c_dstrain_dot =
583 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->resCdStrainDot);
592 auto t_w = getFTensor0IntegrationWeight();
593 auto t_row_base = row_data.getFTensor0N();
594 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
595 const double alpha = getMeasure() * t_w;
599 t_res_vec(
L) = t_L(
i,
j,
L) * (t_c_dstrain_dot(
i,
j) - t_c_dstrain(
i,
j));
606 const auto row_base = alpha * t_row_base;
607 auto t_col_base = col_data.getFTensor0N(gg, 0);
609 t_mat(
L) += (row_base * t_col_base) * t_res_vec(
L);
615 for (; rr != nb_row_base_functions; ++rr)
623 const std::string row_field_name,
const std::string col_field_name,
624 boost::shared_ptr<CommonData> common_data_ptr)
627 commonDataPtr(common_data_ptr) {
632 EntitiesFieldData::EntData &row_data,
633 EntitiesFieldData::EntData &col_data) {
636 const auto nb_integration_pts = getGaussPts().size2();
637 const auto nb_row_base_functions = row_data.getN().size2();
639 auto t_res_c_dtau = getFTensor0FromVec(
commonDataPtr->resCdTau);
640 auto next = [&]() { ++t_res_c_dtau; };
642 auto t_w = getFTensor0IntegrationWeight();
643 auto t_row_base = row_data.getFTensor0N();
644 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
645 const double alpha = getMeasure() * t_w;
648 const auto res = alpha * (t_res_c_dtau);
651 auto mat_ptr =
locMat.data().begin();
654 auto t_col_base = col_data.getFTensor0N(gg, 0);
656 *mat_ptr += t_row_base * t_col_base * res;
662 for (; rr < nb_row_base_functions; ++rr)
#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
double hardening(double tau)
double hardening_dtau(double tau)
constexpr auto field_name
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
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
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