Line data Source code
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 :
10 : #include "MultiPlasticityRawComponentAssembler.h"
11 : #include "RankFourTensor.h"
12 :
13 : InputParameters
14 2087 : MultiPlasticityRawComponentAssembler::validParams()
15 : {
16 2087 : InputParameters params = emptyInputParameters();
17 4174 : MooseEnum specialIC("none rock joint", "none");
18 4174 : 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 4174 : 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 2087 : params.addClassDescription("RawComponentAssembler class to calculate yield functions, etc, used "
32 : "in multi-surface finite-strain plasticity");
33 2087 : return params;
34 2087 : }
35 :
36 1561 : MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(
37 1561 : const MooseObject * moose_object)
38 : : UserObjectInterface(moose_object),
39 1561 : _params(moose_object->parameters()),
40 1561 : _num_models(_params.get<std::vector<UserObjectName>>("plastic_models").size()),
41 1561 : _num_surfaces(0),
42 1561 : _specialIC(_params.get<MooseEnum>("specialIC"))
43 : {
44 1561 : _f.resize(_num_models);
45 3971 : for (unsigned model = 0; model < _num_models; ++model)
46 2410 : _f[model] = &getUserObjectByName<TensorMechanicsPlasticModel>(
47 2410 : _params.get<std::vector<UserObjectName>>("plastic_models")[model]);
48 :
49 3971 : for (unsigned model = 0; model < _num_models; ++model)
50 2410 : _num_surfaces += _f[model]->numberSurfaces();
51 :
52 1561 : _model_given_surface.resize(_num_surfaces);
53 1561 : _model_surface_given_surface.resize(_num_surfaces);
54 : unsigned int surface = 0;
55 3971 : for (unsigned model = 0; model < _num_models; ++model)
56 5600 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
57 : {
58 3190 : _model_given_surface[surface] = model;
59 3190 : _model_surface_given_surface[surface] = model_surface;
60 3190 : surface++;
61 : }
62 :
63 1561 : _surfaces_given_model.resize(_num_models);
64 : surface = 0;
65 3971 : for (unsigned model = 0; model < _num_models; ++model)
66 : {
67 2410 : _surfaces_given_model[model].resize(_f[model]->numberSurfaces());
68 5600 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
69 3190 : _surfaces_given_model[model][model_surface] = surface++;
70 : }
71 :
72 : // check the plastic_models for specialIC
73 1561 : if (_specialIC == "rock")
74 : {
75 24 : if (_num_models != 2)
76 0 : mooseError("Choosing specialIC=rock, you must have plasticity models of type 'TensileMulti "
77 : "MohrCoulombMulti'\n");
78 72 : if (!(_f[0]->modelName().compare("TensileMulti") == 0 &&
79 48 : _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 1561 : if (_specialIC == "joint")
84 : {
85 6 : if (_num_models != 2)
86 0 : mooseError("Choosing specialIC=joint, you must have plasticity models of type "
87 : "'WeakPlaneTensile WeakPlaneShear'\n");
88 18 : if (!(_f[0]->modelName().compare("WeakPlaneTensile") == 0 &&
89 12 : _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 1561 : }
94 :
95 : void
96 753696 : 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 753696 : 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 2058152 : for (unsigned model = 0; model < _num_models; ++model)
109 : {
110 1304456 : activeModelSurfaces(model, active, active_surfaces_of_model);
111 1304456 : if (active_surfaces_of_model.size() > 0)
112 : {
113 1094648 : _f[model]->yieldFunctionV(stress, intnl[model], model_f);
114 1094648 : for (active_surface = active_surfaces_of_model.begin();
115 2324040 : active_surface != active_surfaces_of_model.end();
116 : ++active_surface)
117 1229392 : f.push_back(model_f[*active_surface]);
118 : }
119 : }
120 753696 : }
121 :
122 : void
123 363540 : 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 363540 : 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 1000748 : for (unsigned model = 0; model < _num_models; ++model)
137 : {
138 637208 : activeModelSurfaces(model, active, active_surfaces_of_model);
139 637208 : if (active_surfaces_of_model.size() > 0)
140 : {
141 484066 : _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
142 484066 : for (active_surface = active_surfaces_of_model.begin();
143 985868 : active_surface != active_surfaces_of_model.end();
144 : ++active_surface)
145 501802 : df_dstress.push_back(model_df_dstress[*active_surface]);
146 : }
147 : }
148 363540 : }
149 :
150 : void
151 352810 : 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 352810 : 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 924514 : for (unsigned model = 0; model < _num_models; ++model)
164 : {
165 571704 : activeModelSurfaces(model, active, active_surfaces_of_model);
166 571704 : if (active_surfaces_of_model.size() > 0)
167 : {
168 431140 : _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
169 431140 : for (active_surface = active_surfaces_of_model.begin();
170 877552 : active_surface != active_surfaces_of_model.end();
171 : ++active_surface)
172 446412 : df_dintnl.push_back(model_df_dintnl[*active_surface]);
173 : }
174 : }
175 352810 : }
176 :
177 : void
178 1050224 : 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 1050224 : 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 2770192 : for (unsigned model = 0; model < _num_models; ++model)
191 : {
192 1719968 : activeModelSurfaces(model, active, active_surfaces_of_model);
193 1719968 : if (active_surfaces_of_model.size() > 0)
194 : {
195 1406718 : _f[model]->flowPotentialV(stress, intnl[model], model_r);
196 1406718 : for (active_surface = active_surfaces_of_model.begin();
197 2844980 : active_surface != active_surfaces_of_model.end();
198 : ++active_surface)
199 1438262 : r.push_back(model_r[*active_surface]);
200 : }
201 : }
202 1050224 : }
203 :
204 : void
205 352810 : 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 352810 : 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 924514 : for (unsigned model = 0; model < _num_models; ++model)
219 : {
220 571704 : activeModelSurfaces(model, active, active_surfaces_of_model);
221 571704 : if (active_surfaces_of_model.size() > 0)
222 : {
223 462606 : _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
224 462606 : for (active_surface = active_surfaces_of_model.begin();
225 941716 : active_surface != active_surfaces_of_model.end();
226 : ++active_surface)
227 479110 : dr_dstress.push_back(model_dr_dstress[*active_surface]);
228 : }
229 : }
230 352810 : }
231 :
232 : void
233 352810 : 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 352810 : 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 924514 : for (unsigned model = 0; model < _num_models; ++model)
246 : {
247 571704 : activeModelSurfaces(model, active, active_surfaces_of_model);
248 571704 : if (active_surfaces_of_model.size() > 0)
249 : {
250 462606 : _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
251 462606 : for (active_surface = active_surfaces_of_model.begin();
252 941716 : active_surface != active_surfaces_of_model.end();
253 : ++active_surface)
254 479110 : dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
255 : }
256 : }
257 352810 : }
258 :
259 : void
260 1050224 : 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 1050224 : 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 2770192 : for (unsigned model = 0; model < _num_models; ++model)
273 : {
274 1719968 : activeModelSurfaces(model, active, active_surfaces_of_model);
275 1719968 : if (active_surfaces_of_model.size() > 0)
276 : {
277 1406718 : _f[model]->hardPotentialV(stress, intnl[model], model_h);
278 1406718 : for (active_surface = active_surfaces_of_model.begin();
279 2844980 : active_surface != active_surfaces_of_model.end();
280 : ++active_surface)
281 1438262 : h.push_back(model_h[*active_surface]);
282 : }
283 : }
284 1050224 : }
285 :
286 : void
287 281794 : 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 281794 : 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 746610 : for (unsigned model = 0; model < _num_models; ++model)
301 : {
302 464816 : activeModelSurfaces(model, active, active_surfaces_of_model);
303 464816 : if (active_surfaces_of_model.size() > 0)
304 : {
305 381718 : _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
306 381718 : for (active_surface = active_surfaces_of_model.begin();
307 770956 : active_surface != active_surfaces_of_model.end();
308 : ++active_surface)
309 389238 : dh_dstress.push_back(model_dh_dstress[*active_surface]);
310 : }
311 : }
312 281794 : }
313 :
314 : void
315 281794 : 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 281794 : 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 746610 : for (unsigned model = 0; model < _num_models; ++model)
328 : {
329 464816 : activeModelSurfaces(model, active, active_surfaces_of_model);
330 464816 : if (active_surfaces_of_model.size() > 0)
331 : {
332 381718 : _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
333 381718 : for (active_surface = active_surfaces_of_model.begin();
334 770956 : active_surface != active_surfaces_of_model.end();
335 : ++active_surface)
336 389238 : dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
337 : }
338 : }
339 281794 : }
340 :
341 : void
342 71712 : 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 71712 : if (_specialIC == "rock")
356 8512 : buildActiveConstraintsRock(f, stress, intnl, Eijkl, act);
357 63200 : else if (_specialIC == "joint")
358 7408 : buildActiveConstraintsJoint(f, stress, intnl, Eijkl, act);
359 : else // no specialIC
360 : {
361 55792 : act.resize(0);
362 : unsigned ind = 0;
363 138216 : for (unsigned model = 0; model < _num_models; ++model)
364 : {
365 82424 : std::vector<Real> model_f(0);
366 164848 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
367 82424 : model_f.push_back(f[ind++]);
368 : std::vector<bool> model_act;
369 82424 : RankTwoTensor returned_stress;
370 82424 : _f[model]->activeConstraints(
371 : model_f, stress, intnl[model], Eijkl, model_act, returned_stress);
372 164848 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
373 82424 : act.push_back(model_act[model_surface]);
374 : }
375 : }
376 71712 : }
377 :
378 : void
379 7408 : 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 7408 : 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 7408 : f_single.assign(1, 0);
394 7408 : f_single[0] = f[0];
395 7408 : _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
396 7408 : _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
397 7408 : 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 7408 : f_single.assign(1, 0);
405 7408 : f_single[0] = f[1];
406 7408 : _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_shear, returned_stress);
407 7408 : _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
408 7408 : if (f_single[0] <= _f[0]->_f_tol)
409 : {
410 3712 : act[1] = active_shear[0];
411 3712 : return;
412 : }
413 :
414 : // must be mixed
415 3696 : act[0] = act[1] = true;
416 3696 : return;
417 : }
418 :
419 : void
420 8512 : 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 8512 : 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 8512 : f_single.assign(3, 0);
435 8512 : f_single[0] = f[0];
436 8512 : f_single[1] = f[1];
437 8512 : f_single[2] = f[2];
438 8512 : _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
439 8512 : _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
440 8512 : if (f_single[0] <= _f[1]->_f_tol && f_single[1] <= _f[1]->_f_tol &&
441 6584 : f_single[2] <= _f[1]->_f_tol && f_single[3] <= _f[1]->_f_tol &&
442 10848 : f_single[4] <= _f[1]->_f_tol && f_single[5] <= _f[1]->_f_tol)
443 : {
444 2336 : act[0] = active_tensile[0];
445 2336 : act[1] = active_tensile[1];
446 2336 : act[2] = active_tensile[2];
447 2336 : return;
448 : }
449 :
450 : // next try MC alone
451 6176 : f_single.assign(6, 0);
452 6176 : f_single[0] = f[3];
453 6176 : f_single[1] = f[4];
454 6176 : f_single[2] = f[5];
455 6176 : f_single[3] = f[6];
456 6176 : f_single[4] = f[7];
457 6176 : f_single[5] = f[8];
458 6176 : _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_MC, returned_stress);
459 6176 : _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
460 6176 : if (f_single[0] <= _f[0]->_f_tol && f_single[1] <= _f[0]->_f_tol && f_single[2] <= _f[0]->_f_tol)
461 : {
462 3504 : act[3] = active_MC[0];
463 3504 : act[4] = active_MC[1];
464 3504 : act[5] = active_MC[2];
465 3504 : act[6] = active_MC[3];
466 3504 : act[7] = active_MC[4];
467 3504 : act[8] = active_MC[5];
468 3504 : return;
469 : }
470 :
471 : // must be a mix.
472 : // The possibilities are enumerated below.
473 :
474 : // tensile=edge, MC=tip (two possibilities)
475 2672 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
476 200 : active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
477 2872 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
478 : {
479 200 : act[1] = act[2] = act[6] = true;
480 : act[4] = true;
481 200 : return;
482 : }
483 2472 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
484 392 : active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
485 2864 : 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 2472 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
493 392 : active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
494 2864 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
495 : {
496 392 : act[1] = act[2] = act[4] = act[6] = true;
497 392 : return;
498 : }
499 2080 : 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 2080 : 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 2080 : 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 2080 : 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 2080 : if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
518 640 : active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
519 2720 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
520 : {
521 640 : act[2] = act[6] = true;
522 : act[4] = true;
523 : act[8] = true;
524 640 : return;
525 : }
526 1440 : if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
527 1440 : active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
528 1704 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
529 : {
530 64 : act[2] = act[6] = true;
531 : act[8] = true;
532 64 : return;
533 : }
534 :
535 : // tensile = face, MC=face
536 1376 : if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
537 1376 : active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
538 2552 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
539 : {
540 616 : act[1] = act[2] = act[6] = true;
541 616 : return;
542 : }
543 :
544 : // tensile = face, MC=edge (two possibilites).
545 : act[2] = true; // tensile face
546 760 : act[3] = active_MC[0];
547 760 : act[4] = active_MC[1];
548 760 : act[5] = active_MC[2];
549 760 : act[6] = active_MC[3];
550 760 : act[7] = active_MC[4];
551 760 : act[8] = active_MC[5];
552 760 : 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 3478248 : 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 3478248 : num_successful_plastic_returns = 0;
616 3478248 : yf.resize(0);
617 3478248 : pm.assign(_num_surfaces, 0.0);
618 :
619 3478248 : 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 3478248 : 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 3784082 : for (unsigned model = 0; model < _num_models; ++model)
636 : {
637 3553920 : if (using_custom_return_map)
638 : {
639 3488704 : model_pm.assign(_f[model]->numberSurfaces(), 0.0);
640 3488704 : 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 3488704 : 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 253716 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
655 : ++model_surface)
656 152130 : yf.push_back(model_f[model_surface]);
657 : }
658 3387118 : 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 310944 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
666 : ++model_surface)
667 171912 : 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 3248086 : num_successful_plastic_returns++;
681 : the_single_plastic_model = model;
682 3248086 : 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 65216 : _f[model]->yieldFunctionV(trial_stress, intnl_old[model], model_f);
700 201232 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
701 136016 : yf.push_back(model_f[model_surface]);
702 : }
703 : }
704 :
705 3478248 : 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 230162 : stress = trial_stress;
712 535996 : for (unsigned model = 0; model < _num_models; ++model)
713 305834 : 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 3248086 : std::vector<Real> yf_at_returned_stress(0);
722 : bool all_admissible = true;
723 6497084 : for (unsigned model = 0; model < _num_models; ++model)
724 : {
725 3250950 : if (model == the_single_plastic_model)
726 : {
727 : // no need to spend time calculating the yield function: we know its admissible
728 6531660 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
729 3283574 : yf_at_returned_stress.push_back(model_f[model_surface]);
730 3248086 : continue;
731 3248086 : }
732 2864 : _f[model]->yieldFunctionV(stress, intnl_old[model], model_f);
733 20048 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
734 : {
735 17184 : 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 17184 : yf_at_returned_stress.push_back(model_f[model_surface]);
739 : }
740 2864 : if (!all_admissible)
741 : // no point in continuing computing yield functions
742 : break;
743 : }
744 :
745 3248086 : 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 1952 : stress = trial_stress;
751 5856 : for (unsigned model = 0; model < _num_models; ++model)
752 3904 : intnl[model] = intnl_old[model];
753 : // and calculate the remainder of the yield functions at trial_stress
754 5856 : for (unsigned model = the_single_plastic_model; model < _num_models; ++model)
755 : {
756 3904 : _f[model]->yieldFunctionV(trial_stress, intnl[model], model_f);
757 21472 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
758 17568 : yf.push_back(model_f[model_surface]);
759 : }
760 1952 : num_successful_plastic_returns = 0;
761 1952 : 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 3246134 : yf.resize(0);
768 6529324 : for (unsigned surface = 0; surface < yf_at_returned_stress.size(); ++surface)
769 3283190 : yf.push_back(yf_at_returned_stress[surface]);
770 3246134 : delta_dp = model_delta_dp;
771 6523852 : for (unsigned model_surface = 0; model_surface < _f[the_single_plastic_model]->numberSurfaces();
772 : ++model_surface)
773 : {
774 3277718 : cumulative_pm[_surfaces_given_model[the_single_plastic_model][model_surface]] +=
775 3277718 : model_pm[model_surface];
776 3277718 : pm[_surfaces_given_model[the_single_plastic_model][model_surface]] = model_pm[model_surface];
777 : }
778 3246134 : custom_model = the_single_plastic_model;
779 3246134 : return true;
780 : }
781 :
782 : unsigned int
783 3647434 : MultiPlasticityRawComponentAssembler::modelNumber(unsigned int surface)
784 : {
785 3647434 : return _model_given_surface[surface];
786 : }
787 :
788 : bool
789 2947178 : MultiPlasticityRawComponentAssembler::anyActiveSurfaces(int model, const std::vector<bool> & active)
790 : {
791 3855424 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
792 3235194 : if (active[_surfaces_given_model[model][model_surface]])
793 : return true;
794 : return false;
795 : }
796 :
797 : void
798 1148264 : MultiPlasticityRawComponentAssembler::activeSurfaces(int model,
799 : const std::vector<bool> & active,
800 : std::vector<unsigned int> & active_surfaces)
801 : {
802 1148264 : active_surfaces.resize(0);
803 2444480 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
804 1296216 : if (active[_surfaces_given_model[model][model_surface]])
805 959152 : active_surfaces.push_back(_surfaces_given_model[model][model_surface]);
806 1148264 : }
807 :
808 : void
809 8026344 : MultiPlasticityRawComponentAssembler::activeModelSurfaces(
810 : int model,
811 : const std::vector<bool> & active,
812 : std::vector<unsigned int> & active_surfaces_of_model)
813 : {
814 8026344 : active_surfaces_of_model.resize(0);
815 17723880 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
816 9697536 : if (active[_surfaces_given_model[model][model_surface]])
817 6790826 : active_surfaces_of_model.push_back(model_surface);
818 8026344 : }
|
