Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 :
10 : #include "MultiPlasticityRawComponentAssembler.h"
11 : #include "RankFourTensor.h"
12 :
13 : InputParameters
14 4174 : MultiPlasticityRawComponentAssembler::validParams()
15 : {
16 4174 : InputParameters params = emptyInputParameters();
17 8348 : MooseEnum specialIC("none rock joint", "none");
18 8348 : 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 8348 : 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 4174 : params.addClassDescription("RawComponentAssembler class to calculate yield functions, etc, used "
32 : "in multi-surface finite-strain plasticity");
33 4174 : return params;
34 4174 : }
35 :
36 3122 : MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(
37 3122 : const MooseObject * moose_object)
38 : : UserObjectInterface(moose_object),
39 3122 : _params(moose_object->parameters()),
40 3122 : _num_models(_params.get<std::vector<UserObjectName>>("plastic_models").size()),
41 3122 : _num_surfaces(0),
42 3122 : _specialIC(_params.get<MooseEnum>("specialIC"))
43 : {
44 3122 : _f.resize(_num_models);
45 7942 : for (unsigned model = 0; model < _num_models; ++model)
46 4820 : _f[model] = &getUserObjectByName<SolidMechanicsPlasticModel>(
47 4820 : _params.get<std::vector<UserObjectName>>("plastic_models")[model]);
48 :
49 7942 : for (unsigned model = 0; model < _num_models; ++model)
50 4820 : _num_surfaces += _f[model]->numberSurfaces();
51 :
52 3122 : _model_given_surface.resize(_num_surfaces);
53 3122 : _model_surface_given_surface.resize(_num_surfaces);
54 : unsigned int surface = 0;
55 7942 : for (unsigned model = 0; model < _num_models; ++model)
56 11200 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
57 : {
58 6380 : _model_given_surface[surface] = model;
59 6380 : _model_surface_given_surface[surface] = model_surface;
60 6380 : surface++;
61 : }
62 :
63 3122 : _surfaces_given_model.resize(_num_models);
64 : surface = 0;
65 7942 : for (unsigned model = 0; model < _num_models; ++model)
66 : {
67 4820 : _surfaces_given_model[model].resize(_f[model]->numberSurfaces());
68 11200 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
69 6380 : _surfaces_given_model[model][model_surface] = surface++;
70 : }
71 :
72 : // check the plastic_models for specialIC
73 3122 : if (_specialIC == "rock")
74 : {
75 48 : if (_num_models != 2)
76 0 : mooseError("Choosing specialIC=rock, you must have plasticity models of type 'TensileMulti "
77 : "MohrCoulombMulti'\n");
78 144 : if (!(_f[0]->modelName().compare("TensileMulti") == 0 &&
79 96 : _f[1]->modelName().compare("MohrCoulombMulti") == 0))
80 0 : mooseError("Choosing specialIC=rock, you must have plasticity models of type 'TensileMulti "
81 : "MohrCoulombMulti'\n");
82 : }
83 3122 : if (_specialIC == "joint")
84 : {
85 12 : if (_num_models != 2)
86 0 : mooseError("Choosing specialIC=joint, you must have plasticity models of type "
87 : "'WeakPlaneTensile WeakPlaneShear'\n");
88 36 : if (!(_f[0]->modelName().compare("WeakPlaneTensile") == 0 &&
89 24 : _f[1]->modelName().compare("WeakPlaneShear") == 0))
90 0 : mooseError("Choosing specialIC=joint, you must have plasticity models of type "
91 : "'WeakPlaneTensile WeakPlaneShear'\n");
92 : }
93 3122 : }
94 :
95 : void
96 1420128 : MultiPlasticityRawComponentAssembler::yieldFunction(const RankTwoTensor & stress,
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 1420128 : 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 3848864 : for (unsigned model = 0; model < _num_models; ++model)
109 : {
110 2428736 : activeModelSurfaces(model, active, active_surfaces_of_model);
111 2428736 : if (active_surfaces_of_model.size() > 0)
112 : {
113 2013388 : _f[model]->yieldFunctionV(stress, intnl[model], model_f);
114 2013388 : for (active_surface = active_surfaces_of_model.begin();
115 4295240 : active_surface != active_surfaces_of_model.end();
116 : ++active_surface)
117 2281852 : f.push_back(model_f[*active_surface]);
118 : }
119 : }
120 1420128 : }
121 :
122 : void
123 706816 : MultiPlasticityRawComponentAssembler::dyieldFunction_dstress(
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 706816 : 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 1936176 : for (unsigned model = 0; model < _num_models; ++model)
137 : {
138 1229360 : activeModelSurfaces(model, active, active_surfaces_of_model);
139 1229360 : if (active_surfaces_of_model.size() > 0)
140 : {
141 927892 : _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
142 927892 : for (active_surface = active_surfaces_of_model.begin();
143 1891128 : active_surface != active_surfaces_of_model.end();
144 : ++active_surface)
145 963236 : df_dstress.push_back(model_df_dstress[*active_surface]);
146 : }
147 : }
148 706816 : }
149 :
150 : void
151 685824 : MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl(const RankTwoTensor & stress,
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 685824 : 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 1787312 : for (unsigned model = 0; model < _num_models; ++model)
164 : {
165 1101488 : activeModelSurfaces(model, active, active_surfaces_of_model);
166 1101488 : if (active_surfaces_of_model.size() > 0)
167 : {
168 824640 : _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
169 824640 : for (active_surface = active_surfaces_of_model.begin();
170 1679696 : active_surface != active_surfaces_of_model.end();
171 : ++active_surface)
172 855056 : df_dintnl.push_back(model_df_dintnl[*active_surface]);
173 : }
174 : }
175 685824 : }
176 :
177 : void
178 1994724 : MultiPlasticityRawComponentAssembler::flowPotential(const RankTwoTensor & stress,
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 1994724 : 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 5216836 : for (unsigned model = 0; model < _num_models; ++model)
191 : {
192 3222112 : activeModelSurfaces(model, active, active_surfaces_of_model);
193 3222112 : if (active_surfaces_of_model.size() > 0)
194 : {
195 2602052 : _f[model]->flowPotentialV(stress, intnl[model], model_r);
196 2602052 : for (active_surface = active_surfaces_of_model.begin();
197 5266936 : active_surface != active_surfaces_of_model.end();
198 : ++active_surface)
199 2664884 : r.push_back(model_r[*active_surface]);
200 : }
201 : }
202 1994724 : }
203 :
204 : void
205 685824 : MultiPlasticityRawComponentAssembler::dflowPotential_dstress(
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 685824 : 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 1787312 : for (unsigned model = 0; model < _num_models; ++model)
219 : {
220 1101488 : activeModelSurfaces(model, active, active_surfaces_of_model);
221 1101488 : if (active_surfaces_of_model.size() > 0)
222 : {
223 885908 : _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
224 885908 : for (active_surface = active_surfaces_of_model.begin();
225 1804696 : active_surface != active_surfaces_of_model.end();
226 : ++active_surface)
227 918788 : dr_dstress.push_back(model_dr_dstress[*active_surface]);
228 : }
229 : }
230 685824 : }
231 :
232 : void
233 685824 : MultiPlasticityRawComponentAssembler::dflowPotential_dintnl(const RankTwoTensor & stress,
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 685824 : 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 1787312 : for (unsigned model = 0; model < _num_models; ++model)
246 : {
247 1101488 : activeModelSurfaces(model, active, active_surfaces_of_model);
248 1101488 : if (active_surfaces_of_model.size() > 0)
249 : {
250 885908 : _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
251 885908 : for (active_surface = active_surfaces_of_model.begin();
252 1804696 : active_surface != active_surfaces_of_model.end();
253 : ++active_surface)
254 918788 : dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
255 : }
256 : }
257 685824 : }
258 :
259 : void
260 1994724 : MultiPlasticityRawComponentAssembler::hardPotential(const RankTwoTensor & stress,
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 1994724 : 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 5216836 : for (unsigned model = 0; model < _num_models; ++model)
273 : {
274 3222112 : activeModelSurfaces(model, active, active_surfaces_of_model);
275 3222112 : if (active_surfaces_of_model.size() > 0)
276 : {
277 2602052 : _f[model]->hardPotentialV(stress, intnl[model], model_h);
278 2602052 : for (active_surface = active_surfaces_of_model.begin();
279 5266936 : active_surface != active_surfaces_of_model.end();
280 : ++active_surface)
281 2664884 : h.push_back(model_h[*active_surface]);
282 : }
283 : }
284 1994724 : }
285 :
286 : void
287 544304 : MultiPlasticityRawComponentAssembler::dhardPotential_dstress(
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 544304 : 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 1433424 : for (unsigned model = 0; model < _num_models; ++model)
301 : {
302 889120 : activeModelSurfaces(model, active, active_surfaces_of_model);
303 889120 : if (active_surfaces_of_model.size() > 0)
304 : {
305 724772 : _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
306 724772 : for (active_surface = active_surfaces_of_model.begin();
307 1464520 : active_surface != active_surfaces_of_model.end();
308 : ++active_surface)
309 739748 : dh_dstress.push_back(model_dh_dstress[*active_surface]);
310 : }
311 : }
312 544304 : }
313 :
314 : void
315 544304 : MultiPlasticityRawComponentAssembler::dhardPotential_dintnl(const RankTwoTensor & stress,
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 544304 : 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 1433424 : for (unsigned model = 0; model < _num_models; ++model)
328 : {
329 889120 : activeModelSurfaces(model, active, active_surfaces_of_model);
330 889120 : if (active_surfaces_of_model.size() > 0)
331 : {
332 724772 : _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
333 724772 : for (active_surface = active_surfaces_of_model.begin();
334 1464520 : active_surface != active_surfaces_of_model.end();
335 : ++active_surface)
336 739748 : dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
337 : }
338 : }
339 544304 : }
340 :
341 : void
342 142656 : MultiPlasticityRawComponentAssembler::buildActiveConstraints(const std::vector<Real> & f,
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 142656 : if (_specialIC == "rock")
356 16960 : buildActiveConstraintsRock(f, stress, intnl, Eijkl, act);
357 125696 : else if (_specialIC == "joint")
358 14816 : buildActiveConstraintsJoint(f, stress, intnl, Eijkl, act);
359 : else // no specialIC
360 : {
361 110880 : act.resize(0);
362 : unsigned ind = 0;
363 273936 : for (unsigned model = 0; model < _num_models; ++model)
364 : {
365 163056 : std::vector<Real> model_f(0);
366 326112 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
367 163056 : model_f.push_back(f[ind++]);
368 : std::vector<bool> model_act;
369 163056 : RankTwoTensor returned_stress;
370 163056 : _f[model]->activeConstraints(
371 : model_f, stress, intnl[model], Eijkl, model_act, returned_stress);
372 326112 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
373 163056 : act.push_back(model_act[model_surface]);
374 : }
375 : }
376 142656 : }
377 :
378 : void
379 14816 : MultiPlasticityRawComponentAssembler::buildActiveConstraintsJoint(const std::vector<Real> & f,
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 14816 : 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 14816 : f_single.assign(1, 0);
394 14816 : f_single[0] = f[0];
395 14816 : _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
396 14816 : _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
397 14816 : if (f_single[0] <= _f[1]->_f_tol)
398 : {
399 0 : act[0] = active_tensile[0];
400 0 : return;
401 : }
402 :
403 : // next try shear alone
404 14816 : f_single.assign(1, 0);
405 14816 : f_single[0] = f[1];
406 14816 : _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_shear, returned_stress);
407 14816 : _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
408 14816 : if (f_single[0] <= _f[0]->_f_tol)
409 : {
410 7424 : act[1] = active_shear[0];
411 7424 : return;
412 : }
413 :
414 : // must be mixed
415 7392 : act[0] = act[1] = true;
416 7392 : return;
417 : }
418 :
419 : void
420 16960 : MultiPlasticityRawComponentAssembler::buildActiveConstraintsRock(const std::vector<Real> & f,
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 16960 : 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 16960 : f_single.assign(3, 0);
435 16960 : f_single[0] = f[0];
436 16960 : f_single[1] = f[1];
437 16960 : f_single[2] = f[2];
438 16960 : _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
439 16960 : _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
440 16960 : if (f_single[0] <= _f[1]->_f_tol && f_single[1] <= _f[1]->_f_tol &&
441 13104 : f_single[2] <= _f[1]->_f_tol && f_single[3] <= _f[1]->_f_tol &&
442 21632 : f_single[4] <= _f[1]->_f_tol && f_single[5] <= _f[1]->_f_tol)
443 : {
444 4672 : act[0] = active_tensile[0];
445 4672 : act[1] = active_tensile[1];
446 4672 : act[2] = active_tensile[2];
447 4672 : return;
448 : }
449 :
450 : // next try MC alone
451 12288 : f_single.assign(6, 0);
452 12288 : f_single[0] = f[3];
453 12288 : f_single[1] = f[4];
454 12288 : f_single[2] = f[5];
455 12288 : f_single[3] = f[6];
456 12288 : f_single[4] = f[7];
457 12288 : f_single[5] = f[8];
458 12288 : _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_MC, returned_stress);
459 12288 : _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
460 12288 : if (f_single[0] <= _f[0]->_f_tol && f_single[1] <= _f[0]->_f_tol && f_single[2] <= _f[0]->_f_tol)
461 : {
462 6944 : act[3] = active_MC[0];
463 6944 : act[4] = active_MC[1];
464 6944 : act[5] = active_MC[2];
465 6944 : act[6] = active_MC[3];
466 6944 : act[7] = active_MC[4];
467 6944 : act[8] = active_MC[5];
468 6944 : return;
469 : }
470 :
471 : // must be a mix.
472 : // The possibilities are enumerated below.
473 :
474 : // tensile=edge, MC=tip (two possibilities)
475 5344 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
476 400 : active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
477 5744 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
478 : {
479 400 : act[1] = act[2] = act[6] = true;
480 : act[4] = true;
481 400 : return;
482 : }
483 4944 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
484 784 : active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
485 5728 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
486 : {
487 0 : act[1] = act[2] = act[6] = true; // i don't think act[4] is necessary, is it?!
488 0 : return;
489 : }
490 :
491 : // tensile = edge, MC=edge (two possibilities)
492 4944 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
493 784 : active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
494 5728 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
495 : {
496 784 : act[1] = act[2] = act[4] = act[6] = true;
497 784 : return;
498 : }
499 4160 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
500 0 : active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
501 4160 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
502 : {
503 0 : act[1] = act[2] = act[4] = act[6] = true;
504 0 : return;
505 : }
506 :
507 : // tensile = edge, MC=face
508 4160 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
509 0 : active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
510 4160 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
511 : {
512 0 : act[1] = act[2] = act[6] = true;
513 0 : return;
514 : }
515 :
516 : // tensile = face, MC=tip (two possibilities)
517 4160 : if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
518 1280 : active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
519 5440 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
520 : {
521 1280 : act[2] = act[6] = true;
522 : act[4] = true;
523 : act[8] = true;
524 1280 : return;
525 : }
526 2880 : if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
527 2880 : active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
528 3408 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
529 : {
530 128 : act[2] = act[6] = true;
531 : act[8] = true;
532 128 : return;
533 : }
534 :
535 : // tensile = face, MC=face
536 2752 : if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
537 2752 : active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
538 5104 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
539 : {
540 1232 : act[1] = act[2] = act[6] = true;
541 1232 : return;
542 : }
543 :
544 : // tensile = face, MC=edge (two possibilites).
545 : act[2] = true; // tensile face
546 1520 : act[3] = active_MC[0];
547 1520 : act[4] = active_MC[1];
548 1520 : act[5] = active_MC[2];
549 1520 : act[6] = active_MC[3];
550 1520 : act[7] = active_MC[4];
551 1520 : act[8] = active_MC[5];
552 1520 : return;
553 : }
554 :
555 : /**
556 : * Performs a returnMap for each plastic model.
557 : *
558 : * If all models actually signal "elastic" by
559 : * returning true from their returnMap, and
560 : * by returning model_plastically_active=0, then
561 : * yf will contain the yield function values
562 : * num_successful_plastic_returns will be zero
563 : * intnl = intnl_old
564 : * delta_dp will be unchanged from its input value
565 : * stress will be set to trial_stress
566 : * pm will be zero
567 : * cumulative_pm will be unchanged
568 : * return value will be true
569 : * num_successful_plastic_returns = 0
570 : *
571 : * If only one model signals "plastically active"
572 : * by returning true from its returnMap,
573 : * and by returning model_plastically_active=1, then
574 : * yf will contain the yield function values
575 : * num_successful_plastic_returns will be one
576 : * intnl will be set by the returnMap algorithm
577 : * delta_dp will be set by the returnMap algorithm
578 : * stress will be set by the returnMap algorithm
579 : * pm will be nonzero for the single model, and zero for other models
580 : * cumulative_pm will be updated
581 : * return value will be true
582 : * num_successful_plastic_returns = 1
583 : *
584 : * If >1 model signals "plastically active"
585 : * or if >=1 model's returnMap fails, then
586 : * yf will contain the yield function values
587 : * num_successful_plastic_returns will be set appropriately
588 : * intnl = intnl_old
589 : * delta_dp will be unchanged from its input value
590 : * stress will be set to trial_stress
591 : * pm will be zero
592 : * cumulative_pm will be unchanged
593 : * return value will be true if all returnMap functions returned true, otherwise it will be false
594 : * num_successful_plastic_returns is set appropriately
595 : */
596 : bool
597 6801584 : 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 6801584 : num_successful_plastic_returns = 0;
616 6801584 : yf.resize(0);
617 6801584 : pm.assign(_num_surfaces, 0.0);
618 :
619 6801584 : 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 6801584 : 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 7254180 : for (unsigned model = 0; model < _num_models; ++model)
636 : {
637 6931456 : if (using_custom_return_map)
638 : {
639 6802816 : model_pm.assign(_f[model]->numberSurfaces(), 0.0);
640 6802816 : 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 6802816 : 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 109224 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
655 : ++model_surface)
656 62180 : yf.push_back(model_f[model_surface]);
657 : }
658 6755772 : 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 619584 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
666 : ++model_surface)
667 342672 : 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 6478860 : num_successful_plastic_returns++;
681 : the_single_plastic_model = model;
682 6478860 : 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 6478860 : 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 128640 : _f[model]->yieldFunctionV(trial_stress, intnl_old[model], model_f);
700 398880 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
701 270240 : yf.push_back(model_f[model_surface]);
702 : }
703 : }
704 :
705 6801584 : 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 322724 : stress = trial_stress;
712 775320 : for (unsigned model = 0; model < _num_models; ++model)
713 452596 : 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 6478860 : std::vector<Real> yf_at_returned_stress(0);
722 : bool all_admissible = true;
723 12959480 : for (unsigned model = 0; model < _num_models; ++model)
724 : {
725 6484396 : if (model == the_single_plastic_model)
726 : {
727 : // no need to spend time calculating the yield function: we know its admissible
728 13027928 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
729 6549068 : yf_at_returned_stress.push_back(model_f[model_surface]);
730 6478860 : continue;
731 6478860 : }
732 5536 : _f[model]->yieldFunctionV(stress, intnl_old[model], model_f);
733 38752 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
734 : {
735 33216 : 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 33216 : yf_at_returned_stress.push_back(model_f[model_surface]);
739 : }
740 5536 : if (!all_admissible)
741 : // no point in continuing computing yield functions
742 : break;
743 : }
744 :
745 6478860 : 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 3776 : stress = trial_stress;
751 11328 : for (unsigned model = 0; model < _num_models; ++model)
752 7552 : intnl[model] = intnl_old[model];
753 : // and calculate the remainder of the yield functions at trial_stress
754 11328 : for (unsigned model = the_single_plastic_model; model < _num_models; ++model)
755 : {
756 7552 : _f[model]->yieldFunctionV(trial_stress, intnl[model], model_f);
757 41536 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
758 33984 : yf.push_back(model_f[model_surface]);
759 : }
760 3776 : num_successful_plastic_returns = 0;
761 3776 : 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 6475084 : yf.resize(0);
768 13023384 : for (unsigned surface = 0; surface < yf_at_returned_stress.size(); ++surface)
769 6548300 : yf.push_back(yf_at_returned_stress[surface]);
770 6475084 : delta_dp = model_delta_dp;
771 13012824 : for (unsigned model_surface = 0; model_surface < _f[the_single_plastic_model]->numberSurfaces();
772 : ++model_surface)
773 : {
774 6537740 : cumulative_pm[_surfaces_given_model[the_single_plastic_model][model_surface]] +=
775 6537740 : model_pm[model_surface];
776 6537740 : pm[_surfaces_given_model[the_single_plastic_model][model_surface]] = model_pm[model_surface];
777 : }
778 6475084 : custom_model = the_single_plastic_model;
779 6475084 : return true;
780 : }
781 :
782 : unsigned int
783 6920448 : MultiPlasticityRawComponentAssembler::modelNumber(unsigned int surface)
784 : {
785 6920448 : return _model_given_surface[surface];
786 : }
787 :
788 : bool
789 5557068 : MultiPlasticityRawComponentAssembler::anyActiveSurfaces(int model, const std::vector<bool> & active)
790 : {
791 7353576 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
792 6131308 : if (active[_surfaces_given_model[model][model_surface]])
793 : return true;
794 : return false;
795 : }
796 :
797 : void
798 2120624 : MultiPlasticityRawComponentAssembler::activeSurfaces(int model,
799 : const std::vector<bool> & active,
800 : std::vector<unsigned int> & active_surfaces)
801 : {
802 2120624 : active_surfaces.resize(0);
803 4536256 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
804 2415632 : if (active[_surfaces_given_model[model][model_surface]])
805 1746096 : active_surfaces.push_back(_surfaces_given_model[model][model_surface]);
806 2120624 : }
807 :
808 : void
809 15185024 : MultiPlasticityRawComponentAssembler::activeModelSurfaces(
810 : int model,
811 : const std::vector<bool> & active,
812 : std::vector<unsigned int> & active_surfaces_of_model)
813 : {
814 15185024 : active_surfaces_of_model.resize(0);
815 33702576 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
816 18517552 : if (active[_surfaces_given_model[model][model_surface]])
817 12746984 : active_surfaces_of_model.push_back(model_surface);
818 15185024 : }
|