www.mooseframework.org
MultiPlasticityRawComponentAssembler.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
11 #include "RankFourTensor.h"
12 
15 {
17  MooseEnum specialIC("none rock joint", "none");
18  params.addParam<MooseEnum>("specialIC",
19  specialIC,
20  "For certain combinations of plastic models, the set of active "
21  "constraints can be initialized optimally. 'none': no special "
22  "initialization is performed. For all other choices, the "
23  "plastic_models must be chosen to have the following types. 'rock': "
24  "'TensileMulti MohrCoulombMulti'. 'joint': 'WeakPlaneTensile "
25  "WeakPlaneShear'.");
26  params.addParam<std::vector<UserObjectName>>(
27  "plastic_models",
28  "List of names of user objects that define the plastic models that could "
29  "be active for this material. If no plastic_models are provided, only "
30  "elasticity will be used.");
31  params.addClassDescription("RawComponentAssembler class to calculate yield functions, etc, used "
32  "in multi-surface finite-strain plasticity");
33  return params;
34 }
35 
37  const MooseObject * moose_object)
38  : UserObjectInterface(moose_object),
39  _params(moose_object->parameters()),
40  _num_models(_params.get<std::vector<UserObjectName>>("plastic_models").size()),
41  _num_surfaces(0),
42  _specialIC(_params.get<MooseEnum>("specialIC"))
43 {
44  _f.resize(_num_models);
45  for (unsigned model = 0; model < _num_models; ++model)
46  _f[model] = &getUserObjectByName<SolidMechanicsPlasticModel>(
47  _params.get<std::vector<UserObjectName>>("plastic_models")[model]);
48 
49  for (unsigned model = 0; model < _num_models; ++model)
50  _num_surfaces += _f[model]->numberSurfaces();
51 
54  unsigned int surface = 0;
55  for (unsigned model = 0; model < _num_models; ++model)
56  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
57  {
58  _model_given_surface[surface] = model;
59  _model_surface_given_surface[surface] = model_surface;
60  surface++;
61  }
62 
64  surface = 0;
65  for (unsigned model = 0; model < _num_models; ++model)
66  {
67  _surfaces_given_model[model].resize(_f[model]->numberSurfaces());
68  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
69  _surfaces_given_model[model][model_surface] = surface++;
70  }
71 
72  // check the plastic_models for specialIC
73  if (_specialIC == "rock")
74  {
75  if (_num_models != 2)
76  mooseError("Choosing specialIC=rock, you must have plasticity models of type 'TensileMulti "
77  "MohrCoulombMulti'\n");
78  if (!(_f[0]->modelName().compare("TensileMulti") == 0 &&
79  _f[1]->modelName().compare("MohrCoulombMulti") == 0))
80  mooseError("Choosing specialIC=rock, you must have plasticity models of type 'TensileMulti "
81  "MohrCoulombMulti'\n");
82  }
83  if (_specialIC == "joint")
84  {
85  if (_num_models != 2)
86  mooseError("Choosing specialIC=joint, you must have plasticity models of type "
87  "'WeakPlaneTensile WeakPlaneShear'\n");
88  if (!(_f[0]->modelName().compare("WeakPlaneTensile") == 0 &&
89  _f[1]->modelName().compare("WeakPlaneShear") == 0))
90  mooseError("Choosing specialIC=joint, you must have plasticity models of type "
91  "'WeakPlaneTensile WeakPlaneShear'\n");
92  }
93 }
94 
95 void
97  const std::vector<Real> & intnl,
98  const std::vector<bool> & active,
99  std::vector<Real> & f)
100 {
101  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
102  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
103 
104  f.resize(0);
105  std::vector<unsigned int> active_surfaces_of_model;
106  std::vector<unsigned int>::iterator active_surface;
107  std::vector<Real> model_f;
108  for (unsigned model = 0; model < _num_models; ++model)
109  {
110  activeModelSurfaces(model, active, active_surfaces_of_model);
111  if (active_surfaces_of_model.size() > 0)
112  {
113  _f[model]->yieldFunctionV(stress, intnl[model], model_f);
114  for (active_surface = active_surfaces_of_model.begin();
115  active_surface != active_surfaces_of_model.end();
116  ++active_surface)
117  f.push_back(model_f[*active_surface]);
118  }
119  }
120 }
121 
122 void
124  const RankTwoTensor & stress,
125  const std::vector<Real> & intnl,
126  const std::vector<bool> & active,
127  std::vector<RankTwoTensor> & df_dstress)
128 {
129  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
130  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
131 
132  df_dstress.resize(0);
133  std::vector<unsigned int> active_surfaces_of_model;
134  std::vector<unsigned int>::iterator active_surface;
135  std::vector<RankTwoTensor> model_df_dstress;
136  for (unsigned model = 0; model < _num_models; ++model)
137  {
138  activeModelSurfaces(model, active, active_surfaces_of_model);
139  if (active_surfaces_of_model.size() > 0)
140  {
141  _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
142  for (active_surface = active_surfaces_of_model.begin();
143  active_surface != active_surfaces_of_model.end();
144  ++active_surface)
145  df_dstress.push_back(model_df_dstress[*active_surface]);
146  }
147  }
148 }
149 
150 void
152  const std::vector<Real> & intnl,
153  const std::vector<bool> & active,
154  std::vector<Real> & df_dintnl)
155 {
156  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
157  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
158 
159  df_dintnl.resize(0);
160  std::vector<unsigned int> active_surfaces_of_model;
161  std::vector<unsigned int>::iterator active_surface;
162  std::vector<Real> model_df_dintnl;
163  for (unsigned model = 0; model < _num_models; ++model)
164  {
165  activeModelSurfaces(model, active, active_surfaces_of_model);
166  if (active_surfaces_of_model.size() > 0)
167  {
168  _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
169  for (active_surface = active_surfaces_of_model.begin();
170  active_surface != active_surfaces_of_model.end();
171  ++active_surface)
172  df_dintnl.push_back(model_df_dintnl[*active_surface]);
173  }
174  }
175 }
176 
177 void
179  const std::vector<Real> & intnl,
180  const std::vector<bool> & active,
181  std::vector<RankTwoTensor> & r)
182 {
183  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
184  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
185 
186  r.resize(0);
187  std::vector<unsigned int> active_surfaces_of_model;
188  std::vector<unsigned int>::iterator active_surface;
189  std::vector<RankTwoTensor> model_r;
190  for (unsigned model = 0; model < _num_models; ++model)
191  {
192  activeModelSurfaces(model, active, active_surfaces_of_model);
193  if (active_surfaces_of_model.size() > 0)
194  {
195  _f[model]->flowPotentialV(stress, intnl[model], model_r);
196  for (active_surface = active_surfaces_of_model.begin();
197  active_surface != active_surfaces_of_model.end();
198  ++active_surface)
199  r.push_back(model_r[*active_surface]);
200  }
201  }
202 }
203 
204 void
206  const RankTwoTensor & stress,
207  const std::vector<Real> & intnl,
208  const std::vector<bool> & active,
209  std::vector<RankFourTensor> & dr_dstress)
210 {
211  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
212  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
213 
214  dr_dstress.resize(0);
215  std::vector<unsigned int> active_surfaces_of_model;
216  std::vector<unsigned int>::iterator active_surface;
217  std::vector<RankFourTensor> model_dr_dstress;
218  for (unsigned model = 0; model < _num_models; ++model)
219  {
220  activeModelSurfaces(model, active, active_surfaces_of_model);
221  if (active_surfaces_of_model.size() > 0)
222  {
223  _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
224  for (active_surface = active_surfaces_of_model.begin();
225  active_surface != active_surfaces_of_model.end();
226  ++active_surface)
227  dr_dstress.push_back(model_dr_dstress[*active_surface]);
228  }
229  }
230 }
231 
232 void
234  const std::vector<Real> & intnl,
235  const std::vector<bool> & active,
236  std::vector<RankTwoTensor> & dr_dintnl)
237 {
238  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
239  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
240 
241  dr_dintnl.resize(0);
242  std::vector<unsigned int> active_surfaces_of_model;
243  std::vector<unsigned int>::iterator active_surface;
244  std::vector<RankTwoTensor> model_dr_dintnl;
245  for (unsigned model = 0; model < _num_models; ++model)
246  {
247  activeModelSurfaces(model, active, active_surfaces_of_model);
248  if (active_surfaces_of_model.size() > 0)
249  {
250  _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
251  for (active_surface = active_surfaces_of_model.begin();
252  active_surface != active_surfaces_of_model.end();
253  ++active_surface)
254  dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
255  }
256  }
257 }
258 
259 void
261  const std::vector<Real> & intnl,
262  const std::vector<bool> & active,
263  std::vector<Real> & h)
264 {
265  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
266  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
267 
268  h.resize(0);
269  std::vector<unsigned int> active_surfaces_of_model;
270  std::vector<unsigned int>::iterator active_surface;
271  std::vector<Real> model_h;
272  for (unsigned model = 0; model < _num_models; ++model)
273  {
274  activeModelSurfaces(model, active, active_surfaces_of_model);
275  if (active_surfaces_of_model.size() > 0)
276  {
277  _f[model]->hardPotentialV(stress, intnl[model], model_h);
278  for (active_surface = active_surfaces_of_model.begin();
279  active_surface != active_surfaces_of_model.end();
280  ++active_surface)
281  h.push_back(model_h[*active_surface]);
282  }
283  }
284 }
285 
286 void
288  const RankTwoTensor & stress,
289  const std::vector<Real> & intnl,
290  const std::vector<bool> & active,
291  std::vector<RankTwoTensor> & dh_dstress)
292 {
293  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
294  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
295 
296  dh_dstress.resize(0);
297  std::vector<unsigned int> active_surfaces_of_model;
298  std::vector<unsigned int>::iterator active_surface;
299  std::vector<RankTwoTensor> model_dh_dstress;
300  for (unsigned model = 0; model < _num_models; ++model)
301  {
302  activeModelSurfaces(model, active, active_surfaces_of_model);
303  if (active_surfaces_of_model.size() > 0)
304  {
305  _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
306  for (active_surface = active_surfaces_of_model.begin();
307  active_surface != active_surfaces_of_model.end();
308  ++active_surface)
309  dh_dstress.push_back(model_dh_dstress[*active_surface]);
310  }
311  }
312 }
313 
314 void
316  const std::vector<Real> & intnl,
317  const std::vector<bool> & active,
318  std::vector<Real> & dh_dintnl)
319 {
320  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
321  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
322 
323  dh_dintnl.resize(0);
324  std::vector<unsigned int> active_surfaces_of_model;
325  std::vector<unsigned int>::iterator active_surface;
326  std::vector<Real> model_dh_dintnl;
327  for (unsigned model = 0; model < _num_models; ++model)
328  {
329  activeModelSurfaces(model, active, active_surfaces_of_model);
330  if (active_surfaces_of_model.size() > 0)
331  {
332  _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
333  for (active_surface = active_surfaces_of_model.begin();
334  active_surface != active_surfaces_of_model.end();
335  ++active_surface)
336  dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
337  }
338  }
339 }
340 
341 void
343  const RankTwoTensor & stress,
344  const std::vector<Real> & intnl,
345  const RankFourTensor & Eijkl,
346  std::vector<bool> & act)
347 {
348  mooseAssert(f.size() == _num_surfaces,
349  "buildActiveConstraints called with f.size = " << f.size() << " while there are "
350  << _num_surfaces << " surfaces");
351  mooseAssert(intnl.size() == _num_models,
352  "buildActiveConstraints called with intnl.size = "
353  << intnl.size() << " while there are " << _num_models << " models");
354 
355  if (_specialIC == "rock")
356  buildActiveConstraintsRock(f, stress, intnl, Eijkl, act);
357  else if (_specialIC == "joint")
358  buildActiveConstraintsJoint(f, stress, intnl, Eijkl, act);
359  else // no specialIC
360  {
361  act.resize(0);
362  unsigned ind = 0;
363  for (unsigned model = 0; model < _num_models; ++model)
364  {
365  std::vector<Real> model_f(0);
366  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
367  model_f.push_back(f[ind++]);
368  std::vector<bool> model_act;
369  RankTwoTensor returned_stress;
370  _f[model]->activeConstraints(
371  model_f, stress, intnl[model], Eijkl, model_act, returned_stress);
372  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
373  act.push_back(model_act[model_surface]);
374  }
375  }
376 }
377 
378 void
380  const RankTwoTensor & stress,
381  const std::vector<Real> & intnl,
382  const RankFourTensor & Eijkl,
383  std::vector<bool> & act)
384 {
385  act.assign(2, false);
386 
387  RankTwoTensor returned_stress;
388  std::vector<bool> active_tensile;
389  std::vector<bool> active_shear;
390  std::vector<Real> f_single;
391 
392  // first try tensile alone
393  f_single.assign(1, 0);
394  f_single[0] = f[0];
395  _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
396  _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
397  if (f_single[0] <= _f[1]->_f_tol)
398  {
399  act[0] = active_tensile[0];
400  return;
401  }
402 
403  // next try shear alone
404  f_single.assign(1, 0);
405  f_single[0] = f[1];
406  _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_shear, returned_stress);
407  _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
408  if (f_single[0] <= _f[0]->_f_tol)
409  {
410  act[1] = active_shear[0];
411  return;
412  }
413 
414  // must be mixed
415  act[0] = act[1] = true;
416  return;
417 }
418 
419 void
421  const RankTwoTensor & stress,
422  const std::vector<Real> & intnl,
423  const RankFourTensor & Eijkl,
424  std::vector<bool> & act)
425 {
426  act.assign(9, false);
427 
428  RankTwoTensor returned_stress;
429  std::vector<bool> active_tensile;
430  std::vector<bool> active_MC;
431  std::vector<Real> f_single;
432 
433  // first try tensile alone
434  f_single.assign(3, 0);
435  f_single[0] = f[0];
436  f_single[1] = f[1];
437  f_single[2] = f[2];
438  _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
439  _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
440  if (f_single[0] <= _f[1]->_f_tol && f_single[1] <= _f[1]->_f_tol &&
441  f_single[2] <= _f[1]->_f_tol && f_single[3] <= _f[1]->_f_tol &&
442  f_single[4] <= _f[1]->_f_tol && f_single[5] <= _f[1]->_f_tol)
443  {
444  act[0] = active_tensile[0];
445  act[1] = active_tensile[1];
446  act[2] = active_tensile[2];
447  return;
448  }
449 
450  // next try MC alone
451  f_single.assign(6, 0);
452  f_single[0] = f[3];
453  f_single[1] = f[4];
454  f_single[2] = f[5];
455  f_single[3] = f[6];
456  f_single[4] = f[7];
457  f_single[5] = f[8];
458  _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_MC, returned_stress);
459  _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
460  if (f_single[0] <= _f[0]->_f_tol && f_single[1] <= _f[0]->_f_tol && f_single[2] <= _f[0]->_f_tol)
461  {
462  act[3] = active_MC[0];
463  act[4] = active_MC[1];
464  act[5] = active_MC[2];
465  act[6] = active_MC[3];
466  act[7] = active_MC[4];
467  act[8] = active_MC[5];
468  return;
469  }
470 
471  // must be a mix.
472  // The possibilities are enumerated below.
473 
474  // tensile=edge, MC=tip (two possibilities)
475  if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
476  active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
477  active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
478  {
479  act[1] = act[2] = act[6] = true;
480  act[4] = true;
481  return;
482  }
483  if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
484  active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
485  active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
486  {
487  act[1] = act[2] = act[6] = true; // i don't think act[4] is necessary, is it?!
488  return;
489  }
490 
491  // tensile = edge, MC=edge (two possibilities)
492  if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
493  active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
494  active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
495  {
496  act[1] = act[2] = act[4] = act[6] = true;
497  return;
498  }
499  if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
500  active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
501  active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
502  {
503  act[1] = act[2] = act[4] = act[6] = true;
504  return;
505  }
506 
507  // tensile = edge, MC=face
508  if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
509  active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
510  active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
511  {
512  act[1] = act[2] = act[6] = true;
513  return;
514  }
515 
516  // tensile = face, MC=tip (two possibilities)
517  if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
518  active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
519  active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
520  {
521  act[2] = act[6] = true;
522  act[4] = true;
523  act[8] = true;
524  return;
525  }
526  if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
527  active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
528  active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
529  {
530  act[2] = act[6] = true;
531  act[8] = true;
532  return;
533  }
534 
535  // tensile = face, MC=face
536  if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
537  active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
538  active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
539  {
540  act[1] = act[2] = act[6] = true;
541  return;
542  }
543 
544  // tensile = face, MC=edge (two possibilites).
545  act[2] = true; // tensile face
546  act[3] = active_MC[0];
547  act[4] = active_MC[1];
548  act[5] = active_MC[2];
549  act[6] = active_MC[3];
550  act[7] = active_MC[4];
551  act[8] = active_MC[5];
552  return;
553 }
554 
596 bool
598  const std::vector<Real> & intnl_old,
599  const RankFourTensor & E_ijkl,
600  Real ep_plastic_tolerance,
601  RankTwoTensor & stress,
602  std::vector<Real> & intnl,
603  std::vector<Real> & pm,
604  std::vector<Real> & cumulative_pm,
605  RankTwoTensor & delta_dp,
606  std::vector<Real> & yf,
607  unsigned & num_successful_plastic_returns,
608  unsigned & custom_model)
609 {
610  mooseAssert(intnl_old.size() == _num_models,
611  "returnMapAll: Incorrect size of internal parameters");
612  mooseAssert(intnl.size() == _num_models, "returnMapAll: Incorrect size of internal parameters");
613  mooseAssert(pm.size() == _num_surfaces, "returnMapAll: Incorrect size of pm");
614 
615  num_successful_plastic_returns = 0;
616  yf.resize(0);
617  pm.assign(_num_surfaces, 0.0);
618 
619  RankTwoTensor returned_stress; // each model will give a returned_stress. if only one model is
620  // plastically active, i set stress=returned_stress, so as to
621  // record this returned value
622  std::vector<Real> model_f;
623  RankTwoTensor model_delta_dp;
624  std::vector<Real> model_pm;
625  bool trial_stress_inadmissible;
626  bool successful_return = true;
627  unsigned the_single_plastic_model = 0;
628  bool using_custom_return_map = true;
629 
630  // run through all the plastic models, performing their
631  // returnMap algorithms.
632  // If one finds (trial_stress, intnl) inadmissible and
633  // successfully returns, break from the loop to evaluate
634  // all the others at that returned stress
635  for (unsigned model = 0; model < _num_models; ++model)
636  {
637  if (using_custom_return_map)
638  {
639  model_pm.assign(_f[model]->numberSurfaces(), 0.0);
640  bool model_returned = _f[model]->returnMap(trial_stress,
641  intnl_old[model],
642  E_ijkl,
643  ep_plastic_tolerance,
644  returned_stress,
645  intnl[model],
646  model_pm,
647  model_delta_dp,
648  model_f,
649  trial_stress_inadmissible);
650  if (!trial_stress_inadmissible)
651  {
652  // in the elastic zone: record the yield-function values (returned_stress, intnl, model_pm
653  // and model_delta_dp are undefined)
654  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
655  ++model_surface)
656  yf.push_back(model_f[model_surface]);
657  }
658  else if (trial_stress_inadmissible && !model_returned)
659  {
660  // in the plastic zone, and the customized returnMap failed
661  // for some reason (or wasn't implemented). The coder
662  // should have correctly returned model_f(trial_stress, intnl_old)
663  // so record them
664  // (returned_stress, intnl, model_pm and model_delta_dp are undefined)
665  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
666  ++model_surface)
667  yf.push_back(model_f[model_surface]);
668  // now there's almost zero point in using the custom
669  // returnMap functions
670  using_custom_return_map = false;
671  successful_return = false;
672  }
673  else
674  {
675  // in the plastic zone, and the customized returnMap
676  // succeeded.
677  // record the first returned_stress and delta_dp if everything is going OK
678  // as they could be the actual answer
679  if (trial_stress_inadmissible)
680  num_successful_plastic_returns++;
681  the_single_plastic_model = model;
682  stress = returned_stress;
683  // note that i break here, and don't push_back
684  // model_f to yf. So now yf contains only the values of
685  // model_f from previous models to the_single_plastic_model
686  // also i don't set delta_dp = model_delta_dp yet, because
687  // i might find problems later on
688  // also, don't increment cumulative_pm for the same reason
689 
690  break;
691  }
692  }
693  else
694  {
695  // not using custom returnMap functions because one
696  // has already failed and that one said trial_stress
697  // was inadmissible. So now calculate the yield functions
698  // at the trial stress
699  _f[model]->yieldFunctionV(trial_stress, intnl_old[model], model_f);
700  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
701  yf.push_back(model_f[model_surface]);
702  }
703  }
704 
705  if (num_successful_plastic_returns == 0)
706  {
707  // here either all the models were elastic (successful_return=true),
708  // or some were plastic and either the customized returnMap failed
709  // or wasn't implemented (successful_return=false).
710  // In either case, have to set the following:
711  stress = trial_stress;
712  for (unsigned model = 0; model < _num_models; ++model)
713  intnl[model] = intnl_old[model];
714  return successful_return;
715  }
716 
717  // Now we know that num_successful_plastic_returns == 1 and all the other
718  // models (with model number < the_single_plastic_model) must have been
719  // admissible at (trial_stress, intnl). However, all models might
720  // not be admissible at (trial_stress, intnl), so must check that
721  std::vector<Real> yf_at_returned_stress(0);
722  bool all_admissible = true;
723  for (unsigned model = 0; model < _num_models; ++model)
724  {
725  if (model == the_single_plastic_model)
726  {
727  // no need to spend time calculating the yield function: we know its admissible
728  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
729  yf_at_returned_stress.push_back(model_f[model_surface]);
730  continue;
731  }
732  _f[model]->yieldFunctionV(stress, intnl_old[model], model_f);
733  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
734  {
735  if (model_f[model_surface] > _f[model]->_f_tol)
736  // bummer, this model is not admissible at the returned_stress
737  all_admissible = false;
738  yf_at_returned_stress.push_back(model_f[model_surface]);
739  }
740  if (!all_admissible)
741  // no point in continuing computing yield functions
742  break;
743  }
744 
745  if (!all_admissible)
746  {
747  // we tried using the returned value of stress predicted by
748  // the_single_plastic_model, but it wasn't admissible according
749  // to other plastic models. We need to set:
750  stress = trial_stress;
751  for (unsigned model = 0; model < _num_models; ++model)
752  intnl[model] = intnl_old[model];
753  // and calculate the remainder of the yield functions at trial_stress
754  for (unsigned model = the_single_plastic_model; model < _num_models; ++model)
755  {
756  _f[model]->yieldFunctionV(trial_stress, intnl[model], model_f);
757  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
758  yf.push_back(model_f[model_surface]);
759  }
760  num_successful_plastic_returns = 0;
761  return false;
762  }
763 
764  // So the customized returnMap algorithm can provide a returned
765  // (stress, intnl) configuration, and that is admissible according
766  // to all plastic models
767  yf.resize(0);
768  for (unsigned surface = 0; surface < yf_at_returned_stress.size(); ++surface)
769  yf.push_back(yf_at_returned_stress[surface]);
770  delta_dp = model_delta_dp;
771  for (unsigned model_surface = 0; model_surface < _f[the_single_plastic_model]->numberSurfaces();
772  ++model_surface)
773  {
774  cumulative_pm[_surfaces_given_model[the_single_plastic_model][model_surface]] +=
775  model_pm[model_surface];
776  pm[_surfaces_given_model[the_single_plastic_model][model_surface]] = model_pm[model_surface];
777  }
778  custom_model = the_single_plastic_model;
779  return true;
780 }
781 
782 unsigned int
784 {
785  return _model_given_surface[surface];
786 }
787 
788 bool
790 {
791  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
792  if (active[_surfaces_given_model[model][model_surface]])
793  return true;
794  return false;
795 }
796 
797 void
799  const std::vector<bool> & active,
800  std::vector<unsigned int> & active_surfaces)
801 {
802  active_surfaces.resize(0);
803  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
804  if (active[_surfaces_given_model[model][model_surface]])
805  active_surfaces.push_back(_surfaces_given_model[model][model_surface]);
806 }
807 
808 void
810  int model,
811  const std::vector<bool> & active,
812  std::vector<unsigned int> & active_surfaces_of_model)
813 {
814  active_surfaces_of_model.resize(0);
815  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
816  if (active[_surfaces_given_model[model][model_surface]])
817  active_surfaces_of_model.push_back(model_surface);
818 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
std::vector< unsigned int > _model_surface_given_surface
given a surface number, this returns the corresponding-model&#39;s internal surface number ...
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.
MooseEnum _specialIC
Allows initial set of active constraints to be chosen optimally.
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
void mooseError(Args &&... args)
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.
std::vector< std::pair< R1, R2 > > get(const std::string &param1, const std::string &param2) const
virtual void dhardPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &dh_dintnl)
The derivative of the active hardening potentials with respect to the active 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.
std::vector< unsigned int > _model_given_surface
given a surface number, this returns the model number
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)
InputParameters emptyInputParameters()
virtual void dhardPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dh_dstress)
The derivative of the active hardening potentials with respect to stress By assumption in the Userobj...
void buildActiveConstraintsRock(const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
"Rock" version Constructs a set of active constraints, given the yield functions, f...
void buildActiveConstraintsJoint(const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
"Joint" version Constructs a set of active constraints, given the yield functions, f.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
Real f(Real x)
Test function for Brents method.
void assign(const TypeTensor< T2 > &)
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...
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
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.
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-")
unsigned int _num_models
Number of plastic models for this material.
MultiPlasticityRawComponentAssembler(const MooseObject *moose_object)
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.
void addClassDescription(const std::string &doc_string)
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...
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...
const Elem & get(const ElemType type_in)
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
void activeSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces)
Returns the external surface number(s) of the active surfaces of the given model This may be of size=...