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 2655 : MultiPlasticityRawComponentAssembler::validParams()
15 : {
16 2655 : InputParameters params = emptyInputParameters();
17 5310 : MooseEnum specialIC("none rock joint", "none");
18 5310 : 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 5310 : 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 2655 : params.addClassDescription("RawComponentAssembler class to calculate yield functions, etc, used "
32 : "in multi-surface finite-strain plasticity");
33 2655 : return params;
34 2655 : }
35 :
36 1987 : MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(
37 1987 : const MooseObject * moose_object)
38 : : UserObjectInterface(moose_object),
39 1987 : _params(moose_object->parameters()),
40 1987 : _num_models(_params.get<std::vector<UserObjectName>>("plastic_models").size()),
41 1987 : _num_surfaces(0),
42 1987 : _specialIC(_params.get<MooseEnum>("specialIC"))
43 : {
44 1987 : _f.resize(_num_models);
45 5021 : for (unsigned model = 0; model < _num_models; ++model)
46 3034 : _f[model] = &getUserObjectByName<SolidMechanicsPlasticModel>(
47 3034 : _params.get<std::vector<UserObjectName>>("plastic_models")[model]);
48 :
49 5021 : for (unsigned model = 0; model < _num_models; ++model)
50 3034 : _num_surfaces += _f[model]->numberSurfaces();
51 :
52 1987 : _model_given_surface.resize(_num_surfaces);
53 1987 : _model_surface_given_surface.resize(_num_surfaces);
54 : unsigned int surface = 0;
55 5021 : for (unsigned model = 0; model < _num_models; ++model)
56 7034 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
57 : {
58 4000 : _model_given_surface[surface] = model;
59 4000 : _model_surface_given_surface[surface] = model_surface;
60 4000 : surface++;
61 : }
62 :
63 1987 : _surfaces_given_model.resize(_num_models);
64 : surface = 0;
65 5021 : for (unsigned model = 0; model < _num_models; ++model)
66 : {
67 3034 : _surfaces_given_model[model].resize(_f[model]->numberSurfaces());
68 7034 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
69 4000 : _surfaces_given_model[model][model_surface] = surface++;
70 : }
71 :
72 : // check the plastic_models for specialIC
73 1987 : 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 48 : 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 1987 : 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 12 : 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 1987 : }
94 :
95 : void
96 1117768 : 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 1117768 : 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 2951607 : for (unsigned model = 0; model < _num_models; ++model)
109 : {
110 1833839 : activeModelSurfaces(model, active, active_surfaces_of_model);
111 1833839 : if (active_surfaces_of_model.size() > 0)
112 : {
113 1529484 : _f[model]->yieldFunctionV(stress, intnl[model], model_f);
114 1529484 : for (active_surface = active_surfaces_of_model.begin();
115 3252700 : active_surface != active_surfaces_of_model.end();
116 : ++active_surface)
117 1723216 : f.push_back(model_f[*active_surface]);
118 : }
119 : }
120 1117768 : }
121 :
122 : void
123 547452 : 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 547452 : 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 1452910 : for (unsigned model = 0; model < _num_models; ++model)
137 : {
138 905458 : activeModelSurfaces(model, active, active_surfaces_of_model);
139 905458 : if (active_surfaces_of_model.size() > 0)
140 : {
141 698236 : _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
142 698236 : for (active_surface = active_surfaces_of_model.begin();
143 1422740 : active_surface != active_surfaces_of_model.end();
144 : ++active_surface)
145 724504 : df_dstress.push_back(model_df_dstress[*active_surface]);
146 : }
147 : }
148 547452 : }
149 :
150 : void
151 534362 : 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 534362 : 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 1363235 : for (unsigned model = 0; model < _num_models; ++model)
164 : {
165 828873 : activeModelSurfaces(model, active, active_surfaces_of_model);
166 828873 : if (active_surfaces_of_model.size() > 0)
167 : {
168 636368 : _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
169 636368 : for (active_surface = active_surfaces_of_model.begin();
170 1295372 : active_surface != active_surfaces_of_model.end();
171 : ++active_surface)
172 659004 : df_dintnl.push_back(model_df_dintnl[*active_surface]);
173 : }
174 : }
175 534362 : }
176 :
177 : void
178 1568828 : 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 1568828 : 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 4009172 : for (unsigned model = 0; model < _num_models; ++model)
191 : {
192 2440344 : activeModelSurfaces(model, active, active_surfaces_of_model);
193 2440344 : if (active_surfaces_of_model.size() > 0)
194 : {
195 1986158 : _f[model]->flowPotentialV(stress, intnl[model], model_r);
196 1986158 : for (active_surface = active_surfaces_of_model.begin();
197 4018224 : active_surface != active_surfaces_of_model.end();
198 : ++active_surface)
199 2032066 : r.push_back(model_r[*active_surface]);
200 : }
201 : }
202 1568828 : }
203 :
204 : void
205 534362 : 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 534362 : 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 1363235 : for (unsigned model = 0; model < _num_models; ++model)
219 : {
220 828873 : activeModelSurfaces(model, active, active_surfaces_of_model);
221 828873 : if (active_surfaces_of_model.size() > 0)
222 : {
223 672056 : _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
224 672056 : for (active_surface = active_surfaces_of_model.begin();
225 1368564 : active_surface != active_surfaces_of_model.end();
226 : ++active_surface)
227 696508 : dr_dstress.push_back(model_dr_dstress[*active_surface]);
228 : }
229 : }
230 534362 : }
231 :
232 : void
233 534362 : 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 534362 : 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 1363235 : for (unsigned model = 0; model < _num_models; ++model)
246 : {
247 828873 : activeModelSurfaces(model, active, active_surfaces_of_model);
248 828873 : if (active_surfaces_of_model.size() > 0)
249 : {
250 672056 : _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
251 672056 : for (active_surface = active_surfaces_of_model.begin();
252 1368564 : active_surface != active_surfaces_of_model.end();
253 : ++active_surface)
254 696508 : dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
255 : }
256 : }
257 534362 : }
258 :
259 : void
260 1568828 : 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 1568828 : 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 4009172 : for (unsigned model = 0; model < _num_models; ++model)
273 : {
274 2440344 : activeModelSurfaces(model, active, active_surfaces_of_model);
275 2440344 : if (active_surfaces_of_model.size() > 0)
276 : {
277 1986158 : _f[model]->hardPotentialV(stress, intnl[model], model_h);
278 1986158 : for (active_surface = active_surfaces_of_model.begin();
279 4018224 : active_surface != active_surfaces_of_model.end();
280 : ++active_surface)
281 2032066 : h.push_back(model_h[*active_surface]);
282 : }
283 : }
284 1568828 : }
285 :
286 : void
287 423138 : 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 423138 : 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 1090235 : for (unsigned model = 0; model < _num_models; ++model)
301 : {
302 667097 : activeModelSurfaces(model, active, active_surfaces_of_model);
303 667097 : if (active_surfaces_of_model.size() > 0)
304 : {
305 546676 : _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
306 546676 : for (active_surface = active_surfaces_of_model.begin();
307 1104080 : active_surface != active_surfaces_of_model.end();
308 : ++active_surface)
309 557404 : dh_dstress.push_back(model_dh_dstress[*active_surface]);
310 : }
311 : }
312 423138 : }
313 :
314 : void
315 423138 : 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 423138 : 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 1090235 : for (unsigned model = 0; model < _num_models; ++model)
328 : {
329 667097 : activeModelSurfaces(model, active, active_surfaces_of_model);
330 667097 : if (active_surfaces_of_model.size() > 0)
331 : {
332 546676 : _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
333 546676 : for (active_surface = active_surfaces_of_model.begin();
334 1104080 : active_surface != active_surfaces_of_model.end();
335 : ++active_surface)
336 557404 : dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
337 : }
338 : }
339 423138 : }
340 :
341 : void
342 111604 : 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 111604 : if (_specialIC == "rock")
356 12248 : buildActiveConstraintsRock(f, stress, intnl, Eijkl, act);
357 99356 : else if (_specialIC == "joint")
358 11112 : buildActiveConstraintsJoint(f, stress, intnl, Eijkl, act);
359 : else // no specialIC
360 : {
361 88244 : act.resize(0);
362 : unsigned ind = 0;
363 213916 : for (unsigned model = 0; model < _num_models; ++model)
364 : {
365 125672 : std::vector<Real> model_f(0);
366 251344 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
367 125672 : model_f.push_back(f[ind++]);
368 : std::vector<bool> model_act;
369 125672 : RankTwoTensor returned_stress;
370 125672 : _f[model]->activeConstraints(
371 : model_f, stress, intnl[model], Eijkl, model_act, returned_stress);
372 251344 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
373 125672 : act.push_back(model_act[model_surface]);
374 125672 : }
375 : }
376 111604 : }
377 :
378 : void
379 11112 : 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 11112 : 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 11112 : f_single.assign(1, 0);
394 11112 : f_single[0] = f[0];
395 11112 : _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
396 11112 : _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
397 11112 : 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 11112 : f_single.assign(1, 0);
405 11112 : f_single[0] = f[1];
406 11112 : _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_shear, returned_stress);
407 11112 : _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
408 11112 : if (f_single[0] <= _f[0]->_f_tol)
409 : {
410 5568 : act[1] = active_shear[0];
411 5568 : return;
412 : }
413 :
414 : // must be mixed
415 5544 : act[0] = act[1] = true;
416 5544 : return;
417 11112 : }
418 :
419 : void
420 12248 : 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 12248 : 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 12248 : f_single.assign(3, 0);
435 12248 : f_single[0] = f[0];
436 12248 : f_single[1] = f[1];
437 12248 : f_single[2] = f[2];
438 12248 : _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
439 12248 : _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
440 12248 : if (f_single[0] <= _f[1]->_f_tol && f_single[1] <= _f[1]->_f_tol &&
441 9356 : f_single[2] <= _f[1]->_f_tol && f_single[3] <= _f[1]->_f_tol &&
442 15752 : f_single[4] <= _f[1]->_f_tol && f_single[5] <= _f[1]->_f_tol)
443 : {
444 3504 : act[0] = active_tensile[0];
445 3504 : act[1] = active_tensile[1];
446 3504 : act[2] = active_tensile[2];
447 3504 : return;
448 : }
449 :
450 : // next try MC alone
451 8744 : f_single.assign(6, 0);
452 8744 : f_single[0] = f[3];
453 8744 : f_single[1] = f[4];
454 8744 : f_single[2] = f[5];
455 8744 : f_single[3] = f[6];
456 8744 : f_single[4] = f[7];
457 8744 : f_single[5] = f[8];
458 8744 : _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_MC, returned_stress);
459 8744 : _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
460 8744 : if (f_single[0] <= _f[0]->_f_tol && f_single[1] <= _f[0]->_f_tol && f_single[2] <= _f[0]->_f_tol)
461 : {
462 4768 : act[3] = active_MC[0];
463 4768 : act[4] = active_MC[1];
464 4768 : act[5] = active_MC[2];
465 4768 : act[6] = active_MC[3];
466 4768 : act[7] = active_MC[4];
467 4768 : act[8] = active_MC[5];
468 4768 : return;
469 : }
470 :
471 : // must be a mix.
472 : // The possibilities are enumerated below.
473 :
474 : // tensile=edge, MC=tip (two possibilities)
475 3976 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
476 300 : active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
477 4276 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
478 : {
479 300 : act[1] = act[2] = act[6] = true;
480 : act[4] = true;
481 300 : return;
482 : }
483 3676 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
484 588 : active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
485 4264 : 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 3676 : if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
493 588 : active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
494 4264 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
495 : {
496 588 : act[1] = act[2] = act[4] = act[6] = true;
497 588 : return;
498 : }
499 3088 : 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 3088 : 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 3088 : 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 3088 : 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 3088 : if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
518 928 : active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
519 4016 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
520 : {
521 928 : act[2] = act[6] = true;
522 : act[4] = true;
523 : act[8] = true;
524 928 : return;
525 : }
526 2160 : if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
527 2160 : active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
528 2556 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
529 : {
530 96 : act[2] = act[6] = true;
531 : act[8] = true;
532 96 : return;
533 : }
534 :
535 : // tensile = face, MC=face
536 2064 : if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
537 2064 : active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
538 3828 : active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
539 : {
540 924 : act[1] = act[2] = act[6] = true;
541 924 : return;
542 : }
543 :
544 : // tensile = face, MC=edge (two possibilites).
545 : act[2] = true; // tensile face
546 1140 : act[3] = active_MC[0];
547 1140 : act[4] = active_MC[1];
548 1140 : act[5] = active_MC[2];
549 1140 : act[6] = active_MC[3];
550 1140 : act[7] = active_MC[4];
551 1140 : act[8] = active_MC[5];
552 1140 : return;
553 12248 : }
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 2906940 : 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 2906940 : num_successful_plastic_returns = 0;
616 2906940 : yf.resize(0);
617 2906940 : pm.assign(_num_surfaces, 0.0);
618 :
619 2906940 : 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 2906940 : 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 3253342 : for (unsigned model = 0; model < _num_models; ++model)
636 : {
637 2995976 : if (using_custom_return_map)
638 : {
639 2903672 : model_pm.assign(_f[model]->numberSurfaces(), 0.0);
640 2903672 : 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 2903672 : 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 84532 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
655 : ++model_surface)
656 48278 : yf.push_back(model_f[model_surface]);
657 : }
658 2867418 : 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 485008 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
666 : ++model_surface)
667 267164 : 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 2649574 : num_successful_plastic_returns++;
681 : the_single_plastic_model = model;
682 2649574 : 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 2649574 : 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 92304 : _f[model]->yieldFunctionV(trial_stress, intnl_old[model], model_f);
700 290808 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
701 198504 : yf.push_back(model_f[model_surface]);
702 : }
703 : }
704 :
705 2906940 : 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 257366 : stress = trial_stress;
712 603768 : for (unsigned model = 0; model < _num_models; ++model)
713 346402 : 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 2649574 : std::vector<Real> yf_at_returned_stress(0);
722 : bool all_admissible = true;
723 5300028 : for (unsigned model = 0; model < _num_models; ++model)
724 : {
725 2652342 : if (model == the_single_plastic_model)
726 : {
727 : // no need to spend time calculating the yield function: we know its admissible
728 5353604 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
729 2704030 : yf_at_returned_stress.push_back(model_f[model_surface]);
730 2649574 : continue;
731 2649574 : }
732 2768 : _f[model]->yieldFunctionV(stress, intnl_old[model], model_f);
733 19376 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
734 : {
735 16608 : 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 16608 : yf_at_returned_stress.push_back(model_f[model_surface]);
739 : }
740 2768 : if (!all_admissible)
741 : // no point in continuing computing yield functions
742 : break;
743 : }
744 :
745 2649574 : 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 1888 : stress = trial_stress;
751 5664 : for (unsigned model = 0; model < _num_models; ++model)
752 3776 : intnl[model] = intnl_old[model];
753 : // and calculate the remainder of the yield functions at trial_stress
754 5664 : for (unsigned model = the_single_plastic_model; model < _num_models; ++model)
755 : {
756 3776 : _f[model]->yieldFunctionV(trial_stress, intnl[model], model_f);
757 20768 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
758 16992 : yf.push_back(model_f[model_surface]);
759 : }
760 1888 : num_successful_plastic_returns = 0;
761 1888 : 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 2647686 : yf.resize(0);
768 5351332 : for (unsigned surface = 0; surface < yf_at_returned_stress.size(); ++surface)
769 2703646 : yf.push_back(yf_at_returned_stress[surface]);
770 2647686 : delta_dp = model_delta_dp;
771 5346052 : for (unsigned model_surface = 0; model_surface < _f[the_single_plastic_model]->numberSurfaces();
772 : ++model_surface)
773 : {
774 2698366 : cumulative_pm[_surfaces_given_model[the_single_plastic_model][model_surface]] +=
775 2698366 : model_pm[model_surface];
776 2698366 : pm[_surfaces_given_model[the_single_plastic_model][model_surface]] = model_pm[model_surface];
777 : }
778 2647686 : custom_model = the_single_plastic_model;
779 2647686 : return true;
780 2906940 : }
781 :
782 : unsigned int
783 5236972 : MultiPlasticityRawComponentAssembler::modelNumber(unsigned int surface)
784 : {
785 5236972 : return _model_given_surface[surface];
786 : }
787 :
788 : bool
789 4185695 : MultiPlasticityRawComponentAssembler::anyActiveSurfaces(int model, const std::vector<bool> & active)
790 : {
791 5450822 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
792 4603159 : if (active[_surfaces_given_model[model][model_surface]])
793 : return true;
794 : return false;
795 : }
796 :
797 : void
798 1611471 : MultiPlasticityRawComponentAssembler::activeSurfaces(int model,
799 : const std::vector<bool> & active,
800 : std::vector<unsigned int> & active_surfaces)
801 : {
802 1611471 : active_surfaces.resize(0);
803 3437590 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
804 1826119 : if (active[_surfaces_given_model[model][model_surface]])
805 1335558 : active_surfaces.push_back(_surfaces_given_model[model][model_surface]);
806 1611471 : }
807 :
808 : void
809 11440798 : MultiPlasticityRawComponentAssembler::activeModelSurfaces(
810 : int model,
811 : const std::vector<bool> & active,
812 : std::vector<unsigned int> & active_surfaces_of_model)
813 : {
814 11440798 : active_surfaces_of_model.resize(0);
815 25347920 : for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
816 13907122 : if (active[_surfaces_given_model[model][model_surface]])
817 9678680 : active_surfaces_of_model.push_back(model_surface);
818 11440798 : }
|
