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 4754 : MultiPlasticityRawComponentAssembler::validParams()
15 : {
16 4754 : InputParameters params = emptyInputParameters();
17 9508 : MooseEnum specialIC("none rock joint", "none");
18 9508 : 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 9508 : 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 4754 : params.addClassDescription("RawComponentAssembler class to calculate yield functions, etc, used "
32 : "in multi-surface finite-strain plasticity");
33 4754 : return params;
34 4754 : }
35 :
36 3557 : MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(
37 3557 : const MooseObject * moose_object)
38 : : UserObjectInterface(moose_object),
39 3557 : _params(moose_object->parameters()),
40 3557 : _num_models(_params.get<std::vector<UserObjectName>>("plastic_models").size()),
41 3557 : _num_surfaces(0),
42 3557 : _specialIC(_params.get<MooseEnum>("specialIC"))
43 : {
44 3557 : _f.resize(_num_models);
45 9025 : for (unsigned model = 0; model < _num_models; ++model)
46 5468 : _f[model] = &getUserObjectByName<SolidMechanicsPlasticModel>(
47 5468 : _params.get<std::vector<UserObjectName>>("plastic_models")[model]);
48 :
49 9025 : for (unsigned model = 0; model < _num_models; ++model)
50 5468 : _num_surfaces += _f[model]->numberSurfaces();
51 :
52 3557 : _model_given_surface.resize(_num_surfaces);
53 3557 : _model_surface_given_surface.resize(_num_surfaces);
54 : unsigned int surface = 0;
55 9025 : for (unsigned model = 0; model < _num_models; ++model)
56 12703 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
57 : {
58 7235 : _model_given_surface[surface] = model;
59 7235 : _model_surface_given_surface[surface] = model_surface;
60 7235 : surface++;
61 : }
62 :
63 3557 : _surfaces_given_model.resize(_num_models);
64 : surface = 0;
65 9025 : for (unsigned model = 0; model < _num_models; ++model)
66 : {
67 5468 : _surfaces_given_model[model].resize(_f[model]->numberSurfaces());
68 12703 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
69 7235 : _surfaces_given_model[model][model_surface] = surface++;
70 : }
71 :
72 : // check the plastic_models for specialIC
73 3557 : if (_specialIC == "rock")
74 : {
75 51 : if (_num_models != 2)
76 0 : mooseError("Choosing specialIC=rock, you must have plasticity models of type 'TensileMulti "
77 : "MohrCoulombMulti'\n");
78 102 : if (!(_f[0]->modelName().compare("TensileMulti") == 0 &&
79 102 : _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 3557 : if (_specialIC == "joint")
84 : {
85 15 : if (_num_models != 2)
86 0 : mooseError("Choosing specialIC=joint, you must have plasticity models of type "
87 : "'WeakPlaneTensile WeakPlaneShear'\n");
88 30 : if (!(_f[0]->modelName().compare("WeakPlaneTensile") == 0 &&
89 30 : _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 3557 : }
94 :
95 : void
96 1896520 : 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 1896520 : 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 5156703 : for (unsigned model = 0; model < _num_models; ++model)
109 : {
110 3260183 : activeModelSurfaces(model, active, active_surfaces_of_model);
111 3260183 : if (active_surfaces_of_model.size() > 0)
112 : {
113 2675186 : _f[model]->yieldFunctionV(stress, intnl[model], model_f);
114 2675186 : for (active_surface = active_surfaces_of_model.begin();
115 5737836 : active_surface != active_surfaces_of_model.end();
116 : ++active_surface)
117 3062650 : f.push_back(model_f[*active_surface]);
118 : }
119 : }
120 1896520 : }
121 :
122 : void
123 930996 : 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 930996 : 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 2538614 : for (unsigned model = 0; model < _num_models; ++model)
137 : {
138 1607618 : activeModelSurfaces(model, active, active_surfaces_of_model);
139 1607618 : if (active_surfaces_of_model.size() > 0)
140 : {
141 1217522 : _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
142 1217522 : for (active_surface = active_surfaces_of_model.begin();
143 2484436 : active_surface != active_surfaces_of_model.end();
144 : ++active_surface)
145 1266914 : df_dstress.push_back(model_df_dstress[*active_surface]);
146 : }
147 : }
148 930996 : }
149 :
150 : void
151 906826 : 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 906826 : 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 2372755 : for (unsigned model = 0; model < _num_models; ++model)
164 : {
165 1465929 : activeModelSurfaces(model, active, active_surfaces_of_model);
166 1465929 : if (active_surfaces_of_model.size() > 0)
167 : {
168 1102860 : _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
169 1102860 : for (active_surface = active_surfaces_of_model.begin();
170 2247848 : active_surface != active_surfaces_of_model.end();
171 : ++active_surface)
172 1144988 : df_dintnl.push_back(model_df_dintnl[*active_surface]);
173 : }
174 : }
175 906826 : }
176 :
177 : void
178 2647630 : 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 2647630 : 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 6958366 : for (unsigned model = 0; model < _num_models; ++model)
191 : {
192 4310736 : activeModelSurfaces(model, active, active_surfaces_of_model);
193 4310736 : if (active_surfaces_of_model.size() > 0)
194 : {
195 3441644 : _f[model]->flowPotentialV(stress, intnl[model], model_r);
196 3441644 : for (active_surface = active_surfaces_of_model.begin();
197 6971960 : active_surface != active_surfaces_of_model.end();
198 : ++active_surface)
199 3530316 : r.push_back(model_r[*active_surface]);
200 : }
201 : }
202 2647630 : }
203 :
204 : void
205 906826 : 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 906826 : 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 2372755 : for (unsigned model = 0; model < _num_models; ++model)
219 : {
220 1465929 : activeModelSurfaces(model, active, active_surfaces_of_model);
221 1465929 : if (active_surfaces_of_model.size() > 0)
222 : {
223 1169182 : _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
224 1169182 : for (active_surface = active_surfaces_of_model.begin();
225 2384124 : active_surface != active_surfaces_of_model.end();
226 : ++active_surface)
227 1214942 : dr_dstress.push_back(model_dr_dstress[*active_surface]);
228 : }
229 : }
230 906826 : }
231 :
232 : void
233 906826 : 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 906826 : 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 2372755 : for (unsigned model = 0; model < _num_models; ++model)
246 : {
247 1465929 : activeModelSurfaces(model, active, active_surfaces_of_model);
248 1465929 : if (active_surfaces_of_model.size() > 0)
249 : {
250 1169182 : _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
251 1169182 : for (active_surface = active_surfaces_of_model.begin();
252 2384124 : active_surface != active_surfaces_of_model.end();
253 : ++active_surface)
254 1214942 : dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
255 : }
256 : }
257 906826 : }
258 :
259 : void
260 2647630 : 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 2647630 : 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 6958366 : for (unsigned model = 0; model < _num_models; ++model)
273 : {
274 4310736 : activeModelSurfaces(model, active, active_surfaces_of_model);
275 4310736 : if (active_surfaces_of_model.size() > 0)
276 : {
277 3441644 : _f[model]->hardPotentialV(stress, intnl[model], model_h);
278 3441644 : for (active_surface = active_surfaces_of_model.begin();
279 6971960 : active_surface != active_surfaces_of_model.end();
280 : ++active_surface)
281 3530316 : h.push_back(model_h[*active_surface]);
282 : }
283 : }
284 2647630 : }
285 :
286 : void
287 717002 : 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 717002 : 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 1898555 : for (unsigned model = 0; model < _num_models; ++model)
301 : {
302 1181553 : activeModelSurfaces(model, active, active_surfaces_of_model);
303 1181553 : if (active_surfaces_of_model.size() > 0)
304 : {
305 951790 : _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
306 951790 : for (active_surface = active_surfaces_of_model.begin();
307 1925036 : active_surface != active_surfaces_of_model.end();
308 : ++active_surface)
309 973246 : dh_dstress.push_back(model_dh_dstress[*active_surface]);
310 : }
311 : }
312 717002 : }
313 :
314 : void
315 717002 : 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 717002 : 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 1898555 : for (unsigned model = 0; model < _num_models; ++model)
328 : {
329 1181553 : activeModelSurfaces(model, active, active_surfaces_of_model);
330 1181553 : if (active_surfaces_of_model.size() > 0)
331 : {
332 951790 : _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
333 951790 : for (active_surface = active_surfaces_of_model.begin();
334 1925036 : active_surface != active_surfaces_of_model.end();
335 : ++active_surface)
336 973246 : dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
337 : }
338 : }
339 717002 : }
340 :
341 : void
342 191728 : 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 191728 : if (_specialIC == "rock")
356 24496 : buildActiveConstraintsRock(f, stress, intnl, Eijkl, act);
357 167232 : else if (_specialIC == "joint")
358 22224 : buildActiveConstraintsJoint(f, stress, intnl, Eijkl, act);
359 : else // no specialIC
360 : {
361 145008 : act.resize(0);
362 : unsigned ind = 0;
363 357504 : for (unsigned model = 0; model < _num_models; ++model)
364 : {
365 212496 : std::vector<Real> model_f(0);
366 424992 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
367 212496 : model_f.push_back(f[ind++]);
368 : std::vector<bool> model_act;
369 212496 : RankTwoTensor returned_stress;
370 212496 : _f[model]->activeConstraints(
371 : model_f, stress, intnl[model], Eijkl, model_act, returned_stress);
372 424992 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
373 212496 : act.push_back(model_act[model_surface]);
374 212496 : }
375 : }
376 191728 : }
377 :
378 : void
379 22224 : 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 22224 : 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 22224 : f_single.assign(1, 0);
394 22224 : f_single[0] = f[0];
395 22224 : _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
396 22224 : _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
397 22224 : 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 22224 : f_single.assign(1, 0);
405 22224 : f_single[0] = f[1];
406 22224 : _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_shear, returned_stress);
407 22224 : _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
408 22224 : if (f_single[0] <= _f[0]->_f_tol)
409 : {
410 11136 : act[1] = active_shear[0];
411 11136 : return;
412 : }
413 :
414 : // must be mixed
415 11088 : act[0] = act[1] = true;
416 11088 : return;
417 22224 : }
418 :
419 : void
420 24496 : 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 24496 : 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 24496 : f_single.assign(3, 0);
435 24496 : f_single[0] = f[0];
436 24496 : f_single[1] = f[1];
437 24496 : f_single[2] = f[2];
438 24496 : _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
439 24496 : _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
440 24496 : if (f_single[0] <= _f[1]->_f_tol && f_single[1] <= _f[1]->_f_tol &&
441 18712 : f_single[2] <= _f[1]->_f_tol && f_single[3] <= _f[1]->_f_tol &&
442 31504 : f_single[4] <= _f[1]->_f_tol && f_single[5] <= _f[1]->_f_tol)
443 : {
444 7008 : act[0] = active_tensile[0];
445 7008 : act[1] = active_tensile[1];
446 7008 : act[2] = active_tensile[2];
447 7008 : return;
448 : }
449 :
450 : // next try MC alone
451 17488 : f_single.assign(6, 0);
452 17488 : f_single[0] = f[3];
453 17488 : f_single[1] = f[4];
454 17488 : f_single[2] = f[5];
455 17488 : f_single[3] = f[6];
456 17488 : f_single[4] = f[7];
457 17488 : f_single[5] = f[8];
458 17488 : _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_MC, returned_stress);
459 17488 : _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
460 17488 : if (f_single[0] <= _f[0]->_f_tol && f_single[1] <= _f[0]->_f_tol && f_single[2] <= _f[0]->_f_tol)
461 : {
462 9536 : act[3] = active_MC[0];
463 9536 : act[4] = active_MC[1];
464 9536 : act[5] = active_MC[2];
465 9536 : act[6] = active_MC[3];
466 9536 : act[7] = active_MC[4];
467 9536 : act[8] = active_MC[5];
468 9536 : return;
469 : }
470 :
471 : // must be a mix.
472 : // The possibilities are enumerated below.
473 :
474 : // tensile=edge, MC=tip (two possibilities)
475 7952 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
476 600 : active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
477 8552 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
478 : {
479 600 : act[1] = act[2] = act[6] = true;
480 : act[4] = true;
481 600 : return;
482 : }
483 7352 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
484 1176 : active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
485 8528 : 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 7352 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
493 1176 : active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
494 8528 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
495 : {
496 1176 : act[1] = act[2] = act[4] = act[6] = true;
497 1176 : return;
498 : }
499 6176 : 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 6176 : 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 6176 : 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 6176 : 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 6176 : if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
518 1856 : active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
519 8032 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
520 : {
521 1856 : act[2] = act[6] = true;
522 : act[4] = true;
523 : act[8] = true;
524 1856 : return;
525 : }
526 4320 : if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
527 4320 : active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
528 5112 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
529 : {
530 192 : act[2] = act[6] = true;
531 : act[8] = true;
532 192 : return;
533 : }
534 :
535 : // tensile = face, MC=face
536 4128 : if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
537 4128 : active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
538 7656 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
539 : {
540 1848 : act[1] = act[2] = act[6] = true;
541 1848 : return;
542 : }
543 :
544 : // tensile = face, MC=edge (two possibilites).
545 : act[2] = true; // tensile face
546 2280 : act[3] = active_MC[0];
547 2280 : act[4] = active_MC[1];
548 2280 : act[5] = active_MC[2];
549 2280 : act[6] = active_MC[3];
550 2280 : act[7] = active_MC[4];
551 2280 : act[8] = active_MC[5];
552 2280 : return;
553 24496 : }
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 6956976 : 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 6956976 : num_successful_plastic_returns = 0;
616 6956976 : yf.resize(0);
617 6956976 : pm.assign(_num_surfaces, 0.0);
618 :
619 6956976 : 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 6956976 : 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 7567496 : for (unsigned model = 0; model < _num_models; ++model)
636 : {
637 7131544 : if (using_custom_return_map)
638 : {
639 6958728 : model_pm.assign(_f[model]->numberSurfaces(), 0.0);
640 6958728 : 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 6958728 : 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 151664 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
655 : ++model_surface)
656 86896 : yf.push_back(model_f[model_surface]);
657 : }
658 6893960 : 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 844512 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
666 : ++model_surface)
667 471576 : 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 6521024 : num_successful_plastic_returns++;
681 : the_single_plastic_model = model;
682 6521024 : 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 6521024 : 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 172816 : _f[model]->yieldFunctionV(trial_stress, intnl_old[model], model_f);
700 558032 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
701 385216 : yf.push_back(model_f[model_surface]);
702 : }
703 : }
704 :
705 6956976 : 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 435952 : stress = trial_stress;
712 1046472 : for (unsigned model = 0; model < _num_models; ++model)
713 610520 : 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 6521024 : std::vector<Real> yf_at_returned_stress(0);
722 : bool all_admissible = true;
723 13043808 : for (unsigned model = 0; model < _num_models; ++model)
724 : {
725 6526560 : if (model == the_single_plastic_model)
726 : {
727 : // no need to spend time calculating the yield function: we know its admissible
728 13131608 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
729 6610584 : yf_at_returned_stress.push_back(model_f[model_surface]);
730 6521024 : continue;
731 6521024 : }
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 6521024 : 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 6517248 : yf.resize(0);
768 13127064 : for (unsigned surface = 0; surface < yf_at_returned_stress.size(); ++surface)
769 6609816 : yf.push_back(yf_at_returned_stress[surface]);
770 6517248 : delta_dp = model_delta_dp;
771 13116504 : for (unsigned model_surface = 0; model_surface < _f[the_single_plastic_model]->numberSurfaces();
772 : ++model_surface)
773 : {
774 6599256 : cumulative_pm[_surfaces_given_model[the_single_plastic_model][model_surface]] +=
775 6599256 : model_pm[model_surface];
776 6599256 : pm[_surfaces_given_model[the_single_plastic_model][model_surface]] = model_pm[model_surface];
777 : }
778 6517248 : custom_model = the_single_plastic_model;
779 6517248 : return true;
780 6956976 : }
781 :
782 : unsigned int
783 9250276 : MultiPlasticityRawComponentAssembler::modelNumber(unsigned int surface)
784 : {
785 9250276 : return _model_given_surface[surface];
786 : }
787 :
788 : bool
789 7389669 : MultiPlasticityRawComponentAssembler::anyActiveSurfaces(int model, const std::vector<bool> & active)
790 : {
791 9837634 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
792 8224597 : if (active[_surfaces_given_model[model][model_surface]])
793 : return true;
794 : return false;
795 : }
796 :
797 : void
798 2844807 : MultiPlasticityRawComponentAssembler::activeSurfaces(int model,
799 : const std::vector<bool> & active,
800 : std::vector<unsigned int> & active_surfaces)
801 : {
802 2844807 : active_surfaces.resize(0);
803 6118910 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
804 3274103 : if (active[_surfaces_given_model[model][model_surface]])
805 2315374 : active_surfaces.push_back(_surfaces_given_model[model][model_surface]);
806 2844807 : }
807 :
808 : void
809 20250166 : MultiPlasticityRawComponentAssembler::activeModelSurfaces(
810 : int model,
811 : const std::vector<bool> & active,
812 : std::vector<unsigned int> & active_surfaces_of_model)
813 : {
814 20250166 : active_surfaces_of_model.resize(0);
815 45325220 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
816 25075054 : if (active[_surfaces_given_model[model][model_surface]])
817 16911560 : active_surfaces_of_model.push_back(model_surface);
818 20250166 : }
|
