www.mooseframework.org
MultiPlasticityRawComponentAssembler.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
8 
9 template <>
10 InputParameters
12 {
13  InputParameters params = emptyInputParameters();
14  MooseEnum specialIC("none rock joint", "none");
15  params.addParam<MooseEnum>("specialIC",
16  specialIC,
17  "For certain combinations of plastic models, the set of active "
18  "constraints can be initialized optimally. 'none': no special "
19  "initialization is performed. For all other choices, the "
20  "plastic_models must be chosen to have the following types. 'rock': "
21  "'TensileMulti MohrCoulombMulti'. 'joint': 'WeakPlaneTensile "
22  "WeakPlaneShear'.");
23  params.addParam<std::vector<UserObjectName>>(
24  "plastic_models",
25  "List of names of user objects that define the plastic models that could "
26  "be active for this material. If no plastic_models are provided, only "
27  "elasticity will be used.");
28  params.addClassDescription("RawComponentAssembler class to calculate yield functions, etc, used "
29  "in multi-surface finite-strain plasticity");
30  return params;
31 }
32 
34  const MooseObject * moose_object)
35  : UserObjectInterface(moose_object),
36  _params(moose_object->parameters()),
37  _num_models(_params.get<std::vector<UserObjectName>>("plastic_models").size()),
38  _num_surfaces(0),
39  _specialIC(_params.get<MooseEnum>("specialIC"))
40 {
41  _f.resize(_num_models);
42  for (unsigned model = 0; model < _num_models; ++model)
43  _f[model] = &getUserObjectByName<TensorMechanicsPlasticModel>(
44  _params.get<std::vector<UserObjectName>>("plastic_models")[model]);
45 
46  for (unsigned model = 0; model < _num_models; ++model)
47  _num_surfaces += _f[model]->numberSurfaces();
48 
51  unsigned int surface = 0;
52  for (unsigned model = 0; model < _num_models; ++model)
53  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
54  {
55  _model_given_surface[surface] = model;
56  _model_surface_given_surface[surface] = model_surface;
57  surface++;
58  }
59 
60  _surfaces_given_model.resize(_num_models);
61  surface = 0;
62  for (unsigned model = 0; model < _num_models; ++model)
63  {
64  _surfaces_given_model[model].resize(_f[model]->numberSurfaces());
65  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
66  _surfaces_given_model[model][model_surface] = surface++;
67  }
68 
69  // check the plastic_models for specialIC
70  if (_specialIC == "rock")
71  {
72  if (_num_models != 2)
73  mooseError("Choosing specialIC=rock, you must have plasticity models of type 'TensileMulti "
74  "MohrCoulombMulti'\n");
75  if (!(_f[0]->modelName().compare("TensileMulti") == 0 &&
76  _f[1]->modelName().compare("MohrCoulombMulti") == 0))
77  mooseError("Choosing specialIC=rock, you must have plasticity models of type 'TensileMulti "
78  "MohrCoulombMulti'\n");
79  }
80  if (_specialIC == "joint")
81  {
82  if (_num_models != 2)
83  mooseError("Choosing specialIC=joint, you must have plasticity models of type "
84  "'WeakPlaneTensile WeakPlaneShear'\n");
85  if (!(_f[0]->modelName().compare("WeakPlaneTensile") == 0 &&
86  _f[1]->modelName().compare("WeakPlaneShear") == 0))
87  mooseError("Choosing specialIC=joint, you must have plasticity models of type "
88  "'WeakPlaneTensile WeakPlaneShear'\n");
89  }
90 }
91 
92 void
94  const std::vector<Real> & intnl,
95  const std::vector<bool> & active,
96  std::vector<Real> & f)
97 {
98  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
99  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
100 
101  f.resize(0);
102  std::vector<unsigned int> active_surfaces_of_model;
103  std::vector<unsigned int>::iterator active_surface;
104  std::vector<Real> model_f;
105  for (unsigned model = 0; model < _num_models; ++model)
106  {
107  activeModelSurfaces(model, active, active_surfaces_of_model);
108  if (active_surfaces_of_model.size() > 0)
109  {
110  _f[model]->yieldFunctionV(stress, intnl[model], model_f);
111  for (active_surface = active_surfaces_of_model.begin();
112  active_surface != active_surfaces_of_model.end();
113  ++active_surface)
114  f.push_back(model_f[*active_surface]);
115  }
116  }
117 }
118 
119 void
121  const RankTwoTensor & stress,
122  const std::vector<Real> & intnl,
123  const std::vector<bool> & active,
124  std::vector<RankTwoTensor> & df_dstress)
125 {
126  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
127  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
128 
129  df_dstress.resize(0);
130  std::vector<unsigned int> active_surfaces_of_model;
131  std::vector<unsigned int>::iterator active_surface;
132  std::vector<RankTwoTensor> model_df_dstress;
133  for (unsigned model = 0; model < _num_models; ++model)
134  {
135  activeModelSurfaces(model, active, active_surfaces_of_model);
136  if (active_surfaces_of_model.size() > 0)
137  {
138  _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
139  for (active_surface = active_surfaces_of_model.begin();
140  active_surface != active_surfaces_of_model.end();
141  ++active_surface)
142  df_dstress.push_back(model_df_dstress[*active_surface]);
143  }
144  }
145 }
146 
147 void
149  const std::vector<Real> & intnl,
150  const std::vector<bool> & active,
151  std::vector<Real> & df_dintnl)
152 {
153  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
154  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
155 
156  df_dintnl.resize(0);
157  std::vector<unsigned int> active_surfaces_of_model;
158  std::vector<unsigned int>::iterator active_surface;
159  std::vector<Real> model_df_dintnl;
160  for (unsigned model = 0; model < _num_models; ++model)
161  {
162  activeModelSurfaces(model, active, active_surfaces_of_model);
163  if (active_surfaces_of_model.size() > 0)
164  {
165  _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
166  for (active_surface = active_surfaces_of_model.begin();
167  active_surface != active_surfaces_of_model.end();
168  ++active_surface)
169  df_dintnl.push_back(model_df_dintnl[*active_surface]);
170  }
171  }
172 }
173 
174 void
176  const std::vector<Real> & intnl,
177  const std::vector<bool> & active,
178  std::vector<RankTwoTensor> & r)
179 {
180  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
181  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
182 
183  r.resize(0);
184  std::vector<unsigned int> active_surfaces_of_model;
185  std::vector<unsigned int>::iterator active_surface;
186  std::vector<RankTwoTensor> model_r;
187  for (unsigned model = 0; model < _num_models; ++model)
188  {
189  activeModelSurfaces(model, active, active_surfaces_of_model);
190  if (active_surfaces_of_model.size() > 0)
191  {
192  _f[model]->flowPotentialV(stress, intnl[model], model_r);
193  for (active_surface = active_surfaces_of_model.begin();
194  active_surface != active_surfaces_of_model.end();
195  ++active_surface)
196  r.push_back(model_r[*active_surface]);
197  }
198  }
199 }
200 
201 void
203  const RankTwoTensor & stress,
204  const std::vector<Real> & intnl,
205  const std::vector<bool> & active,
206  std::vector<RankFourTensor> & dr_dstress)
207 {
208  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
209  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
210 
211  dr_dstress.resize(0);
212  std::vector<unsigned int> active_surfaces_of_model;
213  std::vector<unsigned int>::iterator active_surface;
214  std::vector<RankFourTensor> model_dr_dstress;
215  for (unsigned model = 0; model < _num_models; ++model)
216  {
217  activeModelSurfaces(model, active, active_surfaces_of_model);
218  if (active_surfaces_of_model.size() > 0)
219  {
220  _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
221  for (active_surface = active_surfaces_of_model.begin();
222  active_surface != active_surfaces_of_model.end();
223  ++active_surface)
224  dr_dstress.push_back(model_dr_dstress[*active_surface]);
225  }
226  }
227 }
228 
229 void
231  const std::vector<Real> & intnl,
232  const std::vector<bool> & active,
233  std::vector<RankTwoTensor> & dr_dintnl)
234 {
235  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
236  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
237 
238  dr_dintnl.resize(0);
239  std::vector<unsigned int> active_surfaces_of_model;
240  std::vector<unsigned int>::iterator active_surface;
241  std::vector<RankTwoTensor> model_dr_dintnl;
242  for (unsigned model = 0; model < _num_models; ++model)
243  {
244  activeModelSurfaces(model, active, active_surfaces_of_model);
245  if (active_surfaces_of_model.size() > 0)
246  {
247  _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
248  for (active_surface = active_surfaces_of_model.begin();
249  active_surface != active_surfaces_of_model.end();
250  ++active_surface)
251  dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
252  }
253  }
254 }
255 
256 void
258  const std::vector<Real> & intnl,
259  const std::vector<bool> & active,
260  std::vector<Real> & h)
261 {
262  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
263  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
264 
265  h.resize(0);
266  std::vector<unsigned int> active_surfaces_of_model;
267  std::vector<unsigned int>::iterator active_surface;
268  std::vector<Real> model_h;
269  for (unsigned model = 0; model < _num_models; ++model)
270  {
271  activeModelSurfaces(model, active, active_surfaces_of_model);
272  if (active_surfaces_of_model.size() > 0)
273  {
274  _f[model]->hardPotentialV(stress, intnl[model], model_h);
275  for (active_surface = active_surfaces_of_model.begin();
276  active_surface != active_surfaces_of_model.end();
277  ++active_surface)
278  h.push_back(model_h[*active_surface]);
279  }
280  }
281 }
282 
283 void
285  const RankTwoTensor & stress,
286  const std::vector<Real> & intnl,
287  const std::vector<bool> & active,
288  std::vector<RankTwoTensor> & dh_dstress)
289 {
290  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
291  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
292 
293  dh_dstress.resize(0);
294  std::vector<unsigned int> active_surfaces_of_model;
295  std::vector<unsigned int>::iterator active_surface;
296  std::vector<RankTwoTensor> model_dh_dstress;
297  for (unsigned model = 0; model < _num_models; ++model)
298  {
299  activeModelSurfaces(model, active, active_surfaces_of_model);
300  if (active_surfaces_of_model.size() > 0)
301  {
302  _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
303  for (active_surface = active_surfaces_of_model.begin();
304  active_surface != active_surfaces_of_model.end();
305  ++active_surface)
306  dh_dstress.push_back(model_dh_dstress[*active_surface]);
307  }
308  }
309 }
310 
311 void
313  const std::vector<Real> & intnl,
314  const std::vector<bool> & active,
315  std::vector<Real> & dh_dintnl)
316 {
317  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
318  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
319 
320  dh_dintnl.resize(0);
321  std::vector<unsigned int> active_surfaces_of_model;
322  std::vector<unsigned int>::iterator active_surface;
323  std::vector<Real> model_dh_dintnl;
324  for (unsigned model = 0; model < _num_models; ++model)
325  {
326  activeModelSurfaces(model, active, active_surfaces_of_model);
327  if (active_surfaces_of_model.size() > 0)
328  {
329  _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
330  for (active_surface = active_surfaces_of_model.begin();
331  active_surface != active_surfaces_of_model.end();
332  ++active_surface)
333  dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
334  }
335  }
336 }
337 
338 void
340  const RankTwoTensor & stress,
341  const std::vector<Real> & intnl,
342  const RankFourTensor & Eijkl,
343  std::vector<bool> & act)
344 {
345  mooseAssert(f.size() == _num_surfaces,
346  "buildActiveConstraints called with f.size = " << f.size() << " while there are "
347  << _num_surfaces
348  << " surfaces");
349  mooseAssert(intnl.size() == _num_models,
350  "buildActiveConstraints called with intnl.size = " << intnl.size()
351  << " while there are "
352  << _num_models
353  << " 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
597 MultiPlasticityRawComponentAssembler::returnMapAll(const RankTwoTensor & trial_stress,
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
789 MultiPlasticityRawComponentAssembler::anyActiveSurfaces(int model, const std::vector<bool> & active)
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
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.
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)
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.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
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.
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.
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...
InputParameters validParams< MultiPlasticityRawComponentAssembler >()
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...
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=...