17 #include "libmesh/utility.h" 29 "max_NR_iterations>0",
30 "Maximum number of Newton-Raphson iterations allowed");
32 "The Newton-Raphson process is only deemed " 33 "converged if the plastic strain increment " 34 "constraints have L2 norm less than this.");
38 "min_stepsize>0 & min_stepsize<=1",
39 "If ordinary Newton-Raphson + line-search fails, then the applied strain increment is " 40 "subdivided, and the return-map is tried again. This parameter is the minimum fraction of " 41 "applied strain increment that may be applied before the algorithm gives up entirely");
44 "max_stepsize_for_dumb>0 & max_stepsize_for_dumb<=1",
45 "If your deactivation_scheme is 'something_to_dumb', then " 46 "'dumb' will only be used if the stepsize falls below this " 47 "value. This parameter is useful because the 'dumb' scheme is " 48 "computationally expensive");
49 MooseEnum deactivation_scheme(
"optimized safe dumb optimized_to_safe safe_to_dumb " 50 "optimized_to_safe_to_dumb optimized_to_dumb",
53 "deactivation_scheme",
55 "Scheme by which constraints are deactivated. (NOTE: This is irrelevant if there is only " 56 "one yield surface.) safe: return to the yield surface and then deactivate constraints with " 57 "negative plasticity multipliers. optimized: deactivate a constraint as soon as its " 58 "plasticity multiplier becomes negative. dumb: iteratively try all combinations of active " 59 "constraints until the solution is found. You may specify fall-back options. Eg " 60 "optimized_to_safe: first use 'optimized', and if that fails, try the return with 'safe'.");
62 "transverse_direction",
63 "If this parameter is provided, before the return-map algorithm is " 64 "called a rotation is performed so that the 'z' axis in the new " 65 "frame lies along the transverse_direction in the original frame. " 66 "After returning, the inverse rotation is performed. The " 67 "transverse_direction will itself rotate with large strains. This " 68 "is so that transversely-isotropic plasticity models may be easily " 69 "defined in the frame where the isotropy holds in the x-y plane.");
70 params.
addParam<
bool>(
"ignore_failures",
72 "The return-map algorithm will return with the best admissible " 73 "stresses and internal parameters that it can, even if they don't " 74 "fully correspond to the applied strain increment. To speed " 75 "computations, this flag can be set to true, the max_NR_iterations " 76 "set small, and the min_stepsize large.");
77 MooseEnum tangent_operator(
"elastic linear nonlinear",
"nonlinear");
80 "Type of tangent operator to return. 'elastic': return the " 81 "elasticity tensor. 'linear': return the consistent tangent operator " 82 "that is correct for plasticity with yield functions linear in " 83 "stress. 'nonlinear': return the full, general consistent tangent " 84 "operator. The calculations assume the hardening potentials are " 85 "independent of stress and hardening parameters.");
86 params.
addParam<
bool>(
"perform_finite_strain_rotations",
88 "Tensors are correctly rotated in " 89 "finite-strain simulations. For " 90 "optimal performance you can set " 91 "this to 'false' if you are only " 92 "ever using small strains");
100 _max_iter(getParam<unsigned
int>(
"max_NR_iterations")),
101 _min_stepsize(getParam<
Real>(
"min_stepsize")),
102 _max_stepsize_for_dumb(getParam<
Real>(
"max_stepsize_for_dumb")),
103 _ignore_failures(getParam<bool>(
"ignore_failures")),
107 _epp_tol(getParam<
Real>(
"ep_plastic_tolerance")),
115 _n_supplied(parameters.isParamValid(
"transverse_direction")),
119 _perform_finite_strain_rotations(getParam<bool>(
"perform_finite_strain_rotations")),
121 _elasticity_tensor_name(_base_name +
"elasticity_tensor"),
122 _elasticity_tensor(getMaterialPropertyByName<
RankFourTensor>(_elasticity_tensor_name)),
123 _plastic_strain(declareProperty<
RankTwoTensor>(
"plastic_strain")),
124 _plastic_strain_old(getMaterialPropertyOld<
RankTwoTensor>(
"plastic_strain")),
125 _intnl(declareProperty<
std::vector<
Real>>(
"plastic_internal_parameter")),
126 _intnl_old(getMaterialPropertyOld<
std::vector<
Real>>(
"plastic_internal_parameter")),
127 _yf(declareProperty<
std::vector<
Real>>(
"plastic_yield_function")),
128 _iter(declareProperty<
Real>(
"plastic_NR_iterations")),
130 _linesearch_needed(declareProperty<
Real>(
"plastic_linesearch_needed")),
134 _ld_encountered(declareProperty<
Real>(
135 "plastic_linear_dependence_encountered")),
137 _constraints_added(declareProperty<
Real>(
"plastic_constraints_added")),
142 _n_old(getMaterialPropertyOld<
RealVectorValue>(
"plastic_transverse_direction")),
144 _strain_increment(getMaterialPropertyByName<
RankTwoTensor>(_base_name +
"strain_increment")),
145 _total_strain_old(getMaterialPropertyOldByName<
RankTwoTensor>(_base_name +
"total_strain")),
147 getMaterialPropertyByName<
RankTwoTensor>(_base_name +
"rotation_increment")),
149 _stress_old(getMaterialPropertyOld<
RankTwoTensor>(_base_name +
"stress")),
150 _elastic_strain_old(getMaterialPropertyOld<
RankTwoTensor>(_base_name +
"elastic_strain")),
155 hasMaterialProperty<
RankFourTensor>(
"elastic_flexural_rigidity_tensor")),
156 _curvature(_cosserat ? &getMaterialPropertyByName<
RankTwoTensor>(
"curvature") : nullptr),
157 _elastic_flexural_rigidity_tensor(
158 _cosserat ? &getMaterialPropertyByName<
RankFourTensor>(
"elastic_flexural_rigidity_tensor")
160 _couple_stress(_cosserat ? &declareProperty<
RankTwoTensor>(
"couple_stress") : nullptr),
161 _couple_stress_old(_cosserat ? &getMaterialPropertyOld<
RankTwoTensor>(
"couple_stress")
163 _Jacobian_mult_couple(_cosserat ? &declareProperty<
RankFourTensor>(
"couple_Jacobian_mult")
172 mooseError(
"ComputeMultiPlasticityStress: ep_plastic_tolerance must be positive");
179 "ComputeMultiPlasticityStress: transverse_direction vector must not have zero length");
214 mooseError(
"Finite-differencing completed. Exiting with no error");
235 mooseError(
"Finite-differencing completed. Exiting with no error");
240 unsigned int number_iterations;
241 bool linesearch_needed =
false;
242 bool ld_encountered =
false;
243 bool constraints_added =
false;
289 _iter[
_qp] = 1.0 * number_iterations;
353 (*_Jacobian_mult_couple)[
_qp].rotate(
_rot);
355 (*_couple_stress)[
_qp].rotate(
_rot);
371 const std::vector<Real> & intnl_old,
372 std::vector<Real> & intnl,
373 std::vector<Real> & pm,
374 std::vector<Real> & cumulative_pm,
379 std::vector<Real> & yf,
380 unsigned int & iterations,
387 unsigned num_plastic_returns;
392 unsigned custom_model = 0;
393 bool successful_return =
returnMapAll(stress_old + E_ijkl * strain_increment,
408 if (num_plastic_returns == 0)
415 plastic_strain = plastic_strain_old;
416 if (successful_return && final_step)
419 consistent_tangent_operator = E_ijkl;
422 consistent_tangent_operator =
425 return successful_return;
427 else if (num_plastic_returns == 1 && successful_return)
432 plastic_strain = plastic_strain_old + delta_dp;
438 consistent_tangent_operator = E_ijkl;
441 std::vector<Real> custom_model_pm;
442 for (
unsigned surface = 0; surface <
_f[custom_model]->numberSurfaces(); ++surface)
444 consistent_tangent_operator =
445 _f[custom_model]->consistentTangentOperator(stress_old + E_ijkl * strain_increment,
446 intnl_old[custom_model],
456 consistent_tangent_operator =
462 mooseError(
"ComputeMultiPlasticityStress::quickStep should not get here!");
468 const std::vector<Real> & intnl_old,
469 std::vector<Real> & intnl,
474 std::vector<Real> & yf,
475 unsigned int & iterations,
476 bool & linesearch_needed,
477 bool & ld_encountered,
478 bool & constraints_added,
496 bool return_successful =
false;
499 unsigned int iter = 0;
501 Real step_size = 1.0;
502 Real time_simulated = 0.0;
521 unsigned int num_consecutive_successes = 0;
525 return_successful =
returnMap(stress_good,
539 time_simulated + step_size >= 1,
540 consistent_tangent_operator,
544 if (return_successful)
546 num_consecutive_successes += 1;
547 time_simulated += step_size;
549 if (time_simulated < 1.0)
552 stress_good = stress;
553 plastic_strain_good = plastic_strain;
556 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
557 yf_good[surface] = yf[surface];
558 if (num_consecutive_successes >= 2)
561 step_size = std::min(step_size, 1.0 - time_simulated);
566 num_consecutive_successes = 0;
567 stress = stress_good;
568 plastic_strain = plastic_strain_good;
572 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
573 yf[surface] = yf_good[surface];
574 dep = step_size * this_strain_increment;
578 if (!return_successful)
582 stress = stress_good;
583 plastic_strain = plastic_strain_good;
586 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
587 yf[surface] = yf_good[surface];
591 Moose::out <<
"After reducing the stepsize to " << step_size
592 <<
" with original strain increment with L2norm " << this_strain_increment.
L2norm()
593 <<
" the returnMap algorithm failed" << std::endl;
609 return return_successful;
615 const std::vector<Real> & intnl_old,
616 std::vector<Real> & intnl,
621 std::vector<Real> &
f,
623 bool can_revert_to_dumb,
624 bool & linesearch_needed,
625 bool & ld_encountered,
626 bool & constraints_added,
629 std::vector<Real> & cumulative_pm)
635 std::vector<Real> pm;
638 bool successful_return =
quickStep(stress_old,
650 consistent_tangent_operator,
654 if (successful_return)
655 return successful_return;
685 std::vector<Real> ic;
696 std::vector<bool> act;
710 bool single_step_success =
true;
716 std::vector<bool> initial_act;
724 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
725 initial_act[surface] = act[surface];
732 int dumb_iteration = 0;
733 std::vector<unsigned int> dumb_order;
750 std::set<unsigned int> actives_tried;
757 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
759 nr_res2 += 0.5 * Utility::pow<2>(
f[surface] /
_f[
modelNumber(surface)]->_f_tol);
761 successful_return =
false;
763 bool still_finding_solution =
true;
764 while (still_finding_solution)
766 single_step_success =
true;
767 unsigned int local_iter = 0;
770 while (nr_res2 > 0.5 && local_iter++ <
_max_iter && single_step_success)
786 bool nr_good = (nr_res2 <= 0.5 && local_iter <=
_max_iter && single_step_success);
800 bool change_scheme =
false;
801 bool increment_dumb =
false;
803 if (!change_scheme && deact_scheme ==
dumb)
806 still_finding_solution = (change_scheme || increment_dumb);
821 if (!still_finding_solution)
824 successful_return =
false;
829 bool kt_good =
false;
836 if (deact_scheme !=
dumb)
841 still_finding_solution =
843 if (!still_finding_solution)
849 still_finding_solution =
true;
863 bool increment_dumb =
false;
865 still_finding_solution = increment_dumb;
870 if (!still_finding_solution)
874 successful_return =
false;
880 bool admissible =
false;
881 if (nr_good && kt_good)
884 std::vector<Real> all_f;
900 constraints_added =
true;
903 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
905 (!act[surface] && (all_f[surface] >
_f[
modelNumber(surface)]->_f_tol)))
906 act_plus[surface] =
true;
910 constraints_added =
true;
911 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
912 act[surface] = act_plus[surface];
915 add_constraints =
false;
918 bool change_scheme =
false;
919 bool increment_dumb =
false;
921 if (!add_constraints)
924 if (!add_constraints && !change_scheme && deact_scheme ==
dumb)
927 still_finding_solution = (add_constraints || change_scheme || increment_dumb);
942 if (!still_finding_solution)
945 successful_return =
false;
951 successful_return = (nr_good && admissible && kt_good);
952 if (successful_return)
955 if (still_finding_solution)
957 stress = initial_stress;
966 successful_return =
false;
978 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
982 nr_res2 += 0.5 * Utility::pow<2>(
f[ind] /
_f[
modelNumber(surface)]->_f_tol);
989 if (successful_return)
991 plastic_strain += delta_dp;
993 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
994 cumulative_pm[surface] += pm[surface];
997 consistent_tangent_operator =
1009 return successful_return;
1014 const std::vector<Real> & all_f)
1016 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1017 if (!act[surface] && (all_f[surface] >
_f[
modelNumber(surface)]->_f_tol))
1024 bool can_revert_to_dumb)
1049 bool can_revert_to_dumb,
1051 const std::vector<Real> & intnl_old,
1053 std::vector<bool> & act,
1054 int & dumb_iteration,
1055 std::vector<unsigned int> & dumb_order)
1057 if (current_deactivation_scheme ==
optimized &&
1061 current_deactivation_scheme =
safe;
1062 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1063 act[surface] = initial_act[surface];
1065 else if ((current_deactivation_scheme ==
safe &&
1068 can_revert_to_dumb) ||
1070 can_revert_to_dumb))
1072 current_deactivation_scheme =
dumb;
1082 const std::vector<Real> & intnl_old,
1083 std::vector<Real> & intnl,
1084 std::vector<Real> & pm,
1087 std::vector<Real> &
f,
1089 std::vector<Real> & ic,
1090 std::vector<bool> & active,
1092 bool & linesearch_needed,
1093 bool & ld_encountered)
1095 bool successful_step;
1097 Real nr_res2_before_step = nr_res2;
1099 std::vector<Real> intnl_before_step;
1100 std::vector<Real> pm_before_step;
1106 stress_before_step = stress;
1111 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1112 pm_before_step[surface] = pm[surface];
1113 delta_dp_before_step = delta_dp;
1120 std::vector<Real> dpm;
1127 std::vector<bool> deact_ld;
1135 bool constraints_changing =
true;
1136 bool reinstated_actives;
1137 while (constraints_changing)
1140 nrStep(stress, intnl_old, intnl, pm, E_inv, delta_dp, dstress, dpm, dintnl, active, deact_ld);
1142 for (
unsigned surface = 0; surface < deact_ld.size(); ++surface)
1143 if (deact_ld[surface])
1145 ld_encountered =
true;
1168 if (!successful_step)
1170 return successful_step;
1173 constraints_changing =
false;
1176 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1177 if (active[surface] && pm[surface] < 0.0)
1178 constraints_changing =
true;
1181 if (constraints_changing)
1183 stress = stress_before_step;
1184 delta_dp = delta_dp_before_step;
1185 nr_res2 = nr_res2_before_step;
1186 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1188 if (active[surface] && pm[surface] < 0.0)
1191 active[surface] =
false;
1192 pm_before_step[surface] = 0.0;
1197 pm[surface] = pm_before_step[surface];
1202 successful_step =
false;
1203 return successful_step;
1213 if (reinstated_actives)
1215 bool completely_converged =
true;
1216 if (successful_step && nr_res2 < 0.5)
1230 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1231 if (active[surface])
1233 completely_converged =
false;
1236 completely_converged =
false;
1238 if (!completely_converged)
1239 nr_res2 =
residual2(pm,
f, epp, ic, active, deact_ld);
1242 return successful_step;
1247 std::vector<bool> & deactivated_due_to_ld)
1249 bool reinstated_actives =
false;
1250 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1251 if (deactivated_due_to_ld[surface])
1252 reinstated_actives =
true;
1255 return reinstated_actives;
1261 unsigned num_active = 0;
1262 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1263 if (active[surface])
1270 const std::vector<Real> & intnl,
1271 std::vector<Real> & all_f)
1273 std::vector<bool> act;
1278 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1287 const std::vector<Real> & pm,
1288 const std::vector<bool> & active)
1291 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1293 if (active[surface])
1296 if (pm[surface] != 0)
1299 else if (pm[surface] != 0)
1300 mooseError(
"Crash due to plastic multiplier not being zero. This occurred because of poor " 1303 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1304 if (pm[surface] < 0)
1312 const std::vector<Real> & pm,
1313 std::vector<bool> & active)
1315 bool turned_off =
false;
1319 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1321 if (active[surface])
1324 if (pm[surface] != 0)
1327 active[surface] =
false;
1330 else if (pm[surface] != 0)
1331 mooseError(
"Crash due to plastic multiplier not being zero. This occurred because of poor " 1338 int surface_to_turn_off = -1;
1340 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1341 if (pm[surface] < min_pm)
1343 min_pm = pm[surface];
1344 surface_to_turn_off = surface;
1346 if (surface_to_turn_off >= 0)
1347 active[surface_to_turn_off] =
false;
1353 const std::vector<Real> &
f,
1355 const std::vector<Real> & ic,
1356 const std::vector<bool> & active,
1357 const std::vector<bool> & deactivated_due_to_ld)
1362 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1363 if (active[surface])
1365 if (!deactivated_due_to_ld[surface])
1367 if (!(pm[surface] == 0 &&
f[ind] <= 0))
1368 nr_res2 += 0.5 * Utility::pow<2>(
f[ind] /
_f[
modelNumber(surface)]->_f_tol);
1370 else if (deactivated_due_to_ld[surface] &&
f[ind] > 0)
1371 nr_res2 += 0.5 * Utility::pow<2>(
f[ind] /
_f[
modelNumber(surface)]->_f_tol);
1378 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1379 active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
1383 nr_res2 += 0.5 * Utility::pow<2>(ic[ind++] /
_f[
model]->_ic_tol);
1391 const std::vector<Real> & intnl_old,
1392 std::vector<Real> & intnl,
1393 std::vector<Real> & pm,
1397 const std::vector<Real> & dpm,
1398 const std::vector<Real> & dintnl,
1399 std::vector<Real> &
f,
1401 std::vector<Real> & ic,
1402 const std::vector<bool> & active,
1403 const std::vector<bool> & deactivated_due_to_ld,
1404 bool & linesearch_needed)
1417 Real slope = -2 * nr_res2;
1424 std::vector<Real> ls_pm;
1425 ls_pm.resize(pm.size());
1431 std::vector<Real> ls_intnl;
1432 ls_intnl.resize(intnl.size());
1438 std::vector<RankTwoTensor> r;
1445 ls_delta_dp = delta_dp - E_inv * dstress * lam;
1446 for (
unsigned a = 0;
a < intnl.size(); ++
a)
1447 ls_intnl[
a] = intnl[
a] + dintnl[
a] * lam;
1448 ls_stress = stress + dstress * lam;
1454 nr_res2 =
residual2(ls_pm,
f, epp, ic, active, deactivated_due_to_ld);
1456 if (nr_res2 < f0 + 1E-4 * lam * slope)
1458 else if (lam < lam_min)
1464 ls_delta_dp = delta_dp;
1465 for (
unsigned a = 0;
a < intnl.size(); ++
a)
1466 ls_intnl[
a] = intnl[
a];
1469 ls_stress, intnl_old, ls_intnl, ls_pm, ls_delta_dp,
f, r, epp, ic, active);
1470 nr_res2 =
residual2(ls_pm,
f, epp, ic, active, deactivated_due_to_ld);
1473 else if (lam == 1.0)
1476 tmp_lam = -slope / 2.0 / (nr_res2 - f0 - slope);
1481 Real rhs1 = nr_res2 - f0 - lam * slope;
1482 Real rhs2 = f2 - f0 - lam2 * slope;
1483 Real a = (rhs1 / Utility::pow<2>(lam) - rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
1485 (-lam2 * rhs1 / Utility::pow<2>(lam) + lam * rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
1487 tmp_lam = -slope / (2.0 *
b);
1490 Real disc = Utility::pow<2>(
b) - 3 *
a * slope;
1492 tmp_lam = 0.5 * lam;
1494 tmp_lam = (-
b + std::sqrt(disc)) / (3.0 *
a);
1496 tmp_lam = -slope / (
b + std::sqrt(disc));
1498 if (tmp_lam > 0.5 * lam)
1499 tmp_lam = 0.5 * lam;
1503 lam = std::max(tmp_lam, 0.1 * lam);
1507 linesearch_needed =
true;
1513 delta_dp = ls_delta_dp;
1514 for (
unsigned a = 0;
a < intnl.size(); ++
a)
1515 intnl[
a] = ls_intnl[
a];
1523 const std::vector<Real> & intnl,
1524 std::vector<unsigned int> & dumb_order)
1526 if (dumb_order.size() != 0)
1529 std::vector<bool> act;
1532 std::vector<Real>
f;
1534 std::vector<RankTwoTensor> df_dstress;
1537 typedef std::pair<Real, unsigned> pair_for_sorting;
1539 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1541 dist[surface].first =
f[surface] / df_dstress[surface].L2norm();
1542 dist[surface].second = surface;
1544 std::sort(dist.begin(), dist.end());
1554 const std::vector<unsigned int> & dumb_order,
1555 std::vector<bool> & act)
1557 dumb_iteration += 1;
1558 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1559 act[dumb_order[surface]] =
1575 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1577 num += (1 << surface);
1584 const std::vector<Real> & intnl,
1586 const std::vector<Real> & pm_this_step,
1587 const std::vector<Real> & cumulative_pm)
1597 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1598 act_at_some_step[surface] = (cumulative_pm[surface] > 0);
1605 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1606 act_vary[surface] = (pm_this_step[surface] > 0);
1608 std::vector<RankTwoTensor> df_dstress;
1610 std::vector<Real> df_dintnl;
1612 std::vector<RankTwoTensor> r;
1614 std::vector<RankFourTensor> dr_dstress_at_some_step;
1616 std::vector<RankTwoTensor> dr_dintnl_at_some_step;
1618 std::vector<Real> h;
1627 std::vector<RankTwoTensor> r_minus_stuff;
1629 for (
unsigned surface1 = 0; surface1 <
_num_surfaces; ++surface1)
1630 if (act_vary[surface1])
1632 r_minus_stuff.push_back(r[ind1]);
1634 for (
unsigned surface2 = 0; surface2 <
_num_surfaces; ++surface2)
1635 if (act_at_some_step[surface2])
1639 r_minus_stuff.back() -=
1640 cumulative_pm[surface2] * dr_dintnl_at_some_step[ind2] * h[ind1];
1647 unsigned int num_currently_active = 0;
1648 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
1649 if (act_vary[surface])
1650 num_currently_active += 1;
1658 std::vector<PetscScalar> zzz;
1659 zzz.assign(num_currently_active * num_currently_active, 0.0);
1663 for (
unsigned surface1 = 0; surface1 <
_num_surfaces; ++surface1)
1664 if (act_vary[surface1])
1667 for (
unsigned surface2 = 0; surface2 <
_num_surfaces; ++surface2)
1668 if (act_vary[surface2])
1670 r2 = df_dstress[ind1] * (E_ijkl * r_minus_stuff[ind2]);
1671 zzz[ind1 * num_currently_active + ind2] += r2(0, 0) + r2(1, 1) + r2(2, 2);
1673 zzz[ind1 * num_currently_active + ind2] += df_dintnl[ind1] * h[ind2];
1679 if (num_currently_active > 0)
1695 for (
unsigned surface1 = 0; surface1 <
_num_surfaces; ++surface1)
1696 if (act_vary[surface1])
1700 for (
unsigned surface2 = 0; surface2 <
_num_surfaces; ++surface2)
1701 if (act_vary[surface2])
1704 for (
unsigned i = 0; i < 3; i++)
1705 for (
unsigned j = 0;
j < 3;
j++)
1706 for (
unsigned k = 0;
k < 3;
k++)
1707 for (
unsigned l = 0; l < 3; l++)
1708 strain_coeff(i,
j,
k, l) -=
1709 part1(i,
j) * part2(
k, l) * zzz[ind1 * num_currently_active + ind2];
1716 return strain_coeff;
1722 for (
unsigned surface1 = 0; surface1 <
_num_surfaces; ++surface1)
1723 if (act_at_some_step[surface1])
1725 part3 += cumulative_pm[surface1] * E_ijkl * dr_dstress_at_some_step[ind1];
1729 stress_coeff += part3;
1735 for (
unsigned surface1 = 0; surface1 <
_num_surfaces; ++surface1)
1736 if (act_vary[surface1])
1740 for (
unsigned surface2 = 0; surface2 <
_num_surfaces; ++surface2)
1741 if (act_vary[surface2])
1744 for (
unsigned i = 0; i < 3; i++)
1745 for (
unsigned j = 0;
j < 3;
j++)
1746 for (
unsigned k = 0;
k < 3;
k++)
1747 for (
unsigned l = 0; l < 3; l++)
1748 stress_coeff(i,
j,
k, l) -=
1749 part1(i,
j) * part2(
k, l) * zzz[ind1 * num_currently_active + ind2];
1766 s_inv = stress_coeff.
invSymm();
1770 return strain_coeff;
1774 return s_inv * strain_coeff;
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
bool _perform_finite_strain_rotations
whether to perform the rotations necessary in finite-strain simulations
virtual bool reinstateLinearDependentConstraints(std::vector< bool > &deactivated_due_to_ld)
makes all deactivated_due_to_ld false, and if >0 of them were initially true, returns true ...
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
ComputeStressBase is the base class for stress tensors computed from MOOSE's strain calculators...
auto norm() const -> decltype(std::norm(Real()))
bool canAddConstraints(const std::vector< bool > &act, const std::vector< Real > &all_f)
initIdentitySymmetricFour
virtual void dyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
The derivative of the active yield function(s) with respect to stress.
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
quickStep_called_from_t
The functions from which quickStep can be called.
MooseEnum _fspb_debug
none - don't do any debugging crash - currently inactive jacobian - check the jacobian entries jacobi...
virtual bool singleStep(Real &nr_res2, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, RankTwoTensor &delta_dp, const RankFourTensor &E_inv, std::vector< Real > &f, RankTwoTensor &epp, std::vector< Real > &ic, std::vector< bool > &active, DeactivationSchemeEnum deactivation_scheme, bool &linesearch_needed, bool &ld_encountered)
Performs a single Newton-Raphson + linesearch step Constraints are deactivated and the step is re-don...
static InputParameters validParams()
virtual void buildActiveConstraints(const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
Constructs a set of active constraints, given the yield functions, f.
const MaterialProperty< RankTwoTensor > & _rotation_increment
Rotation increment (coming from ComputeIncrementalStrain, for example)
virtual void preReturnMap()
const MaterialProperty< RealVectorValue > & _n_old
old value of transverse direction
MaterialProperty< Real > & _constraints_added
Whether constraints were added in during the latest Newton-Raphson process (1 if true, 0 otherwise)
MaterialProperty< std::vector< Real > > & _intnl
internal parameters
bool returnMapAll(const RankTwoTensor &trial_stress, const std::vector< Real > &intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &stress, std::vector< Real > &intnl, std::vector< Real > &pm, std::vector< Real > &cumulative_pm, RankTwoTensor &delta_dp, std::vector< Real > &yf, unsigned &num_successful_plastic_returns, unsigned &custom_model)
Performs a returnMap for each plastic model using their inbuilt returnMap functions.
bool canChangeScheme(DeactivationSchemeEnum current_deactivation_scheme, bool can_revert_to_dumb)
RankFourTensor _my_flexural_rigidity_tensor
Flexual rigidity tensor that can be rotated by this class (ie, its not const)
static constexpr std::size_t dim
ComputeMultiPlasticityStress performs the return-map algorithm and associated stress updates for plas...
const MaterialProperty< RankTwoTensor > & _strain_increment
strain increment (coming from ComputeIncrementalStrain, for example)
registerMooseObject("SolidMechanicsApp", ComputeMultiPlasticityStress)
TangentOperatorEnum
The type of tangent operator to return. tangent operator = d(stress_rate)/d(strain_rate).
static InputParameters validParams()
virtual void initQpStatefulProperties()
const MaterialProperty< RankTwoTensor > & _stress_old
Old value of stress.
unsigned int _max_iter
Maximum number of Newton-Raphson iterations allowed.
RankTwoTensor _my_strain_increment
Strain increment that can be rotated by this class, and split into multiple increments (ie...
unsigned int activeCombinationNumber(const std::vector< bool > &act)
virtual void applyKuhnTucker(const std::vector< Real > &f, const std::vector< Real > &pm, std::vector< bool > &active)
Checks Kuhn-Tucker conditions, and alters "active" if appropriate.
virtual void yieldFunction(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
The active yield function(s)
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
virtual void postReturnMap()
void rotate(const RankTwoTensorTempl< Real > &R)
bool _ignore_failures
Even if the returnMap fails, return the best values found for stress and internal parameters...
RankFourTensor _my_elasticity_tensor
Elasticity tensor that can be rotated by this class (ie, its not const)
MaterialProperty< Real > & _iter
Number of Newton-Raphson iterations used in the return-map.
const MaterialProperty< RankTwoTensor > & _elastic_strain_old
Old value of elastic strain.
TensorValue< Real > RealTensorValue
unsigned int _num_surfaces
Number of surfaces within the plastic models.
virtual void initQpStatefulProperties() override
Real f(Real x)
Test function for Brents method.
bool _n_supplied
User supplied the transverse direction vector.
MaterialProperty< Real > & _ld_encountered
Whether linear-dependence was encountered in the latest Newton-Raphson process (1 if true...
void changeScheme(const std::vector< bool > &initial_act, bool can_revert_to_dumb, const RankTwoTensor &initial_stress, const std::vector< Real > &intnl_old, DeactivationSchemeEnum ¤t_deactivation_scheme, std::vector< bool > &act, int &dumb_iteration, std::vector< unsigned int > &dumb_order)
virtual bool checkKuhnTucker(const std::vector< Real > &f, const std::vector< Real > &pm, const std::vector< bool > &active)
Checks Kuhn-Tucker conditions, and alters "active" if appropriate.
virtual bool lineSearch(Real &nr_res2, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, const RankFourTensor &E_inv, RankTwoTensor &delta_dp, const RankTwoTensor &dstress, const std::vector< Real > &dpm, const std::vector< Real > &dintnl, std::vector< Real > &f, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, bool &linesearch_needed)
Performs a line search.
RankTwoTensor rot(const RankTwoTensor &tens)
MaterialProperty< std::vector< Real > > & _yf
yield functions
void assign(const TypeTensor< T2 > &)
RankTwoTensorTempl< Real > rotated(const RankTwoTensorTempl< Real > &R) const
virtual void dflowPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dr_dintnl)
The derivative of the active flow potentials with respect to the active internal parameters The UserO...
void checkSolution(const RankFourTensor &E_inv)
Checks that Ax does equal b in the NR procedure.
void checkDerivatives()
Checks the derivatives, eg dyieldFunction_dstress by using finite difference approximations.
MaterialProperty< RealVectorValue > & _n
current value of transverse direction
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to 'active' ...
std::vector< Real > _dummy_pm
dummy "consistency parameters" (plastic multipliers) used in quickStep when called from computeQpStre...
RankFourTensorTempl< T > transposeMajor() const
virtual void incrementDumb(int &dumb_iteration, const std::vector< unsigned int > &dumb_order, std::vector< bool > &act)
Increments "dumb_iteration" by 1, and sets "act" appropriately (act[alpha] = true iff alpha_th bit of...
MaterialProperty< Real > & _linesearch_needed
Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise) ...
virtual void flowPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
The active flow potential(s) - one for each yield function.
void rotate(const TypeTensor< T > &R)
bool canIncrementDumb(int dumb_iteration)
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
RealVectorValue _n_input
the supplied transverse direction vector
Real _max_stepsize_for_dumb
"dumb" deactivation will only be used if the stepsize falls below this quantity
const MaterialProperty< std::vector< Real > > & _intnl_old
old values of internal parameters
enum ComputeMultiPlasticityStress::TangentOperatorEnum _tangent_operator_type
virtual bool checkAdmissible(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< Real > &all_f)
Checks whether the yield functions are in the admissible region.
GenericRealTensorValue< is_ad > rotVecToZ(GenericRealVectorValue< is_ad > vec)
virtual bool returnMap(const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &f, unsigned int &iter, bool can_revert_to_dumb, bool &linesearch_needed, bool &ld_encountered, bool &constraints_added, bool final_step, RankFourTensor &consistent_tangent_operator, std::vector< Real > &cumulative_pm)
Implements the return map.
virtual Real residual2(const std::vector< Real > &pm, const std::vector< Real > &f, const RankTwoTensor &epp, const std::vector< Real > &ic, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld)
The residual-squared.
ComputeMultiPlasticityStress(const InputParameters ¶meters)
MaterialProperty< RankTwoTensor > & _plastic_strain
plastic strain
unsigned int _num_models
Number of plastic models for this material.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void computeQpStress()
Compute the stress and store it in the _stress material property for the current quadrature point...
static const std::string alpha
virtual void dflowPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankFourTensor > &dr_dstress)
The derivative of the active flow potential(s) with respect to stress.
Real _epp_tol
Tolerance on the plastic strain increment ("direction") constraint.
RankFourTensor consistentTangentOperator(const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &E_ijkl, const std::vector< Real > &pm_this_step, const std::vector< Real > &cumulative_pm)
Computes the consistent tangent operator (another name for the jacobian = d(stress_rate)/d(strain_rat...
MaterialProperty< RankTwoTensor > & _elastic_strain
Elastic strain material property.
IntRange< T > make_range(T beg, T end)
Real _min_stepsize
Minimum fraction of applied strain that may be applied during adaptive stepsizing.
void mooseError(Args &&... args) const
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.
MultiPlasticityDebugger computes various finite-difference things to help developers remove bugs in t...
virtual void hardPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &h)
The active hardening potentials (one for each internal parameter and for each yield function) by assu...
RankTwoTensor _my_curvature
Curvature that can be rotated by this class, and split into multiple increments (ie, its not const)
std::vector< Real > _cumulative_pm
the sum of the plastic multipliers over all the sub-steps.
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
MaterialProperty< RankTwoTensor > & _stress
Stress material property.
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
Old value of plastic strain.
virtual unsigned int numberActive(const std::vector< bool > &active)
counts the number of active constraints
RealTensorValue _rot
rotation matrix that takes _n to (0, 0, 1)
virtual bool quickStep(const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, std::vector< Real > &cumulative_pm, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &yf, unsigned int &iterations, RankFourTensor &consistent_tangent_operator, const quickStep_called_from_t called_from, bool final_step)
Attempts to find an admissible (stress, intnl) by using the customized return-map algorithms defined ...
bool _cosserat
whether Cosserat mechanics should be used
void buildDumbOrder(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< unsigned int > &dumb_order)
Builds the order which "dumb" activation will take.
virtual void dyieldFunction_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &df_dintnl)
The derivative of active yield function(s) with respect to their internal parameters (the user object...
RankFourTensorTempl< T > invSymm() const
enum ComputeMultiPlasticityStress::DeactivationSchemeEnum _deactivation_scheme
virtual bool plasticStep(const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &yf, unsigned int &iterations, bool &linesearch_needed, bool &ld_encountered, bool &constraints_added, RankFourTensor &consistent_tangent_operator)
performs a plastic step
const MaterialProperty< RankFourTensor > & _elasticity_tensor
Elasticity tensor material property.
virtual void nrStep(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const RankTwoTensor &delta_dp, RankTwoTensor &dstress, std::vector< Real > &dpm, std::vector< Real > &dintnl, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs one Newton-Raphson step.
static const std::string k
virtual void calculateConstraints(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &f, std::vector< RankTwoTensor > &r, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active)
The constraints.
void ErrorVector unsigned int
void checkJacobian(const RankFourTensor &E_inv, const std::vector< Real > &intnl_old)
Checks the full Jacobian, which is just certain linear combinations of the dyieldFunction_dstress, etc, by using finite difference approximations.
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
static InputParameters validParams()