LCOV - code coverage report
Current view: top level - src/utils - MultiPlasticityRawComponentAssembler.C (source / functions) Hit Total Coverage
Test: idaholab/moose solid_mechanics: f45d79 Lines: 308 322 95.7 %
Date: 2025-07-25 05:00:39 Functions: 19 19 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "MultiPlasticityRawComponentAssembler.h"
      11             : #include "RankFourTensor.h"
      12             : 
      13             : InputParameters
      14        4174 : MultiPlasticityRawComponentAssembler::validParams()
      15             : {
      16        4174 :   InputParameters params = emptyInputParameters();
      17        8348 :   MooseEnum specialIC("none rock joint", "none");
      18        8348 :   params.addParam<MooseEnum>("specialIC",
      19             :                              specialIC,
      20             :                              "For certain combinations of plastic models, the set of active "
      21             :                              "constraints can be initialized optimally.  'none': no special "
      22             :                              "initialization is performed.  For all other choices, the "
      23             :                              "plastic_models must be chosen to have the following types.  'rock': "
      24             :                              "'TensileMulti MohrCoulombMulti'.  'joint': 'WeakPlaneTensile "
      25             :                              "WeakPlaneShear'.");
      26        8348 :   params.addParam<std::vector<UserObjectName>>(
      27             :       "plastic_models",
      28             :       "List of names of user objects that define the plastic models that could "
      29             :       "be active for this material.  If no plastic_models are provided, only "
      30             :       "elasticity will be used.");
      31        4174 :   params.addClassDescription("RawComponentAssembler class to calculate yield functions, etc, used "
      32             :                              "in multi-surface finite-strain plasticity");
      33        4174 :   return params;
      34        4174 : }
      35             : 
      36        3122 : MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(
      37        3122 :     const MooseObject * moose_object)
      38             :   : UserObjectInterface(moose_object),
      39        3122 :     _params(moose_object->parameters()),
      40        3122 :     _num_models(_params.get<std::vector<UserObjectName>>("plastic_models").size()),
      41        3122 :     _num_surfaces(0),
      42        3122 :     _specialIC(_params.get<MooseEnum>("specialIC"))
      43             : {
      44        3122 :   _f.resize(_num_models);
      45        7942 :   for (unsigned model = 0; model < _num_models; ++model)
      46        4820 :     _f[model] = &getUserObjectByName<SolidMechanicsPlasticModel>(
      47        4820 :         _params.get<std::vector<UserObjectName>>("plastic_models")[model]);
      48             : 
      49        7942 :   for (unsigned model = 0; model < _num_models; ++model)
      50        4820 :     _num_surfaces += _f[model]->numberSurfaces();
      51             : 
      52        3122 :   _model_given_surface.resize(_num_surfaces);
      53        3122 :   _model_surface_given_surface.resize(_num_surfaces);
      54             :   unsigned int surface = 0;
      55        7942 :   for (unsigned model = 0; model < _num_models; ++model)
      56       11200 :     for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
      57             :     {
      58        6380 :       _model_given_surface[surface] = model;
      59        6380 :       _model_surface_given_surface[surface] = model_surface;
      60        6380 :       surface++;
      61             :     }
      62             : 
      63        3122 :   _surfaces_given_model.resize(_num_models);
      64             :   surface = 0;
      65        7942 :   for (unsigned model = 0; model < _num_models; ++model)
      66             :   {
      67        4820 :     _surfaces_given_model[model].resize(_f[model]->numberSurfaces());
      68       11200 :     for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
      69        6380 :       _surfaces_given_model[model][model_surface] = surface++;
      70             :   }
      71             : 
      72             :   // check the plastic_models for specialIC
      73        3122 :   if (_specialIC == "rock")
      74             :   {
      75          48 :     if (_num_models != 2)
      76           0 :       mooseError("Choosing specialIC=rock, you must have plasticity models of type 'TensileMulti "
      77             :                  "MohrCoulombMulti'\n");
      78         144 :     if (!(_f[0]->modelName().compare("TensileMulti") == 0 &&
      79          96 :           _f[1]->modelName().compare("MohrCoulombMulti") == 0))
      80           0 :       mooseError("Choosing specialIC=rock, you must have plasticity models of type 'TensileMulti "
      81             :                  "MohrCoulombMulti'\n");
      82             :   }
      83        3122 :   if (_specialIC == "joint")
      84             :   {
      85          12 :     if (_num_models != 2)
      86           0 :       mooseError("Choosing specialIC=joint, you must have plasticity models of type "
      87             :                  "'WeakPlaneTensile WeakPlaneShear'\n");
      88          36 :     if (!(_f[0]->modelName().compare("WeakPlaneTensile") == 0 &&
      89          24 :           _f[1]->modelName().compare("WeakPlaneShear") == 0))
      90           0 :       mooseError("Choosing specialIC=joint, you must have plasticity models of type "
      91             :                  "'WeakPlaneTensile WeakPlaneShear'\n");
      92             :   }
      93        3122 : }
      94             : 
      95             : void
      96     1420128 : MultiPlasticityRawComponentAssembler::yieldFunction(const RankTwoTensor & stress,
      97             :                                                     const std::vector<Real> & intnl,
      98             :                                                     const std::vector<bool> & active,
      99             :                                                     std::vector<Real> & f)
     100             : {
     101             :   mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
     102             :   mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
     103             : 
     104     1420128 :   f.resize(0);
     105             :   std::vector<unsigned int> active_surfaces_of_model;
     106             :   std::vector<unsigned int>::iterator active_surface;
     107             :   std::vector<Real> model_f;
     108     3848864 :   for (unsigned model = 0; model < _num_models; ++model)
     109             :   {
     110     2428736 :     activeModelSurfaces(model, active, active_surfaces_of_model);
     111     2428736 :     if (active_surfaces_of_model.size() > 0)
     112             :     {
     113     2013388 :       _f[model]->yieldFunctionV(stress, intnl[model], model_f);
     114     2013388 :       for (active_surface = active_surfaces_of_model.begin();
     115     4295240 :            active_surface != active_surfaces_of_model.end();
     116             :            ++active_surface)
     117     2281852 :         f.push_back(model_f[*active_surface]);
     118             :     }
     119             :   }
     120     1420128 : }
     121             : 
     122             : void
     123      706816 : MultiPlasticityRawComponentAssembler::dyieldFunction_dstress(
     124             :     const RankTwoTensor & stress,
     125             :     const std::vector<Real> & intnl,
     126             :     const std::vector<bool> & active,
     127             :     std::vector<RankTwoTensor> & df_dstress)
     128             : {
     129             :   mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
     130             :   mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
     131             : 
     132      706816 :   df_dstress.resize(0);
     133             :   std::vector<unsigned int> active_surfaces_of_model;
     134             :   std::vector<unsigned int>::iterator active_surface;
     135             :   std::vector<RankTwoTensor> model_df_dstress;
     136     1936176 :   for (unsigned model = 0; model < _num_models; ++model)
     137             :   {
     138     1229360 :     activeModelSurfaces(model, active, active_surfaces_of_model);
     139     1229360 :     if (active_surfaces_of_model.size() > 0)
     140             :     {
     141      927892 :       _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
     142      927892 :       for (active_surface = active_surfaces_of_model.begin();
     143     1891128 :            active_surface != active_surfaces_of_model.end();
     144             :            ++active_surface)
     145      963236 :         df_dstress.push_back(model_df_dstress[*active_surface]);
     146             :     }
     147             :   }
     148      706816 : }
     149             : 
     150             : void
     151      685824 : MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl(const RankTwoTensor & stress,
     152             :                                                             const std::vector<Real> & intnl,
     153             :                                                             const std::vector<bool> & active,
     154             :                                                             std::vector<Real> & df_dintnl)
     155             : {
     156             :   mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
     157             :   mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
     158             : 
     159      685824 :   df_dintnl.resize(0);
     160             :   std::vector<unsigned int> active_surfaces_of_model;
     161             :   std::vector<unsigned int>::iterator active_surface;
     162             :   std::vector<Real> model_df_dintnl;
     163     1787312 :   for (unsigned model = 0; model < _num_models; ++model)
     164             :   {
     165     1101488 :     activeModelSurfaces(model, active, active_surfaces_of_model);
     166     1101488 :     if (active_surfaces_of_model.size() > 0)
     167             :     {
     168      824640 :       _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
     169      824640 :       for (active_surface = active_surfaces_of_model.begin();
     170     1679696 :            active_surface != active_surfaces_of_model.end();
     171             :            ++active_surface)
     172      855056 :         df_dintnl.push_back(model_df_dintnl[*active_surface]);
     173             :     }
     174             :   }
     175      685824 : }
     176             : 
     177             : void
     178     1994724 : MultiPlasticityRawComponentAssembler::flowPotential(const RankTwoTensor & stress,
     179             :                                                     const std::vector<Real> & intnl,
     180             :                                                     const std::vector<bool> & active,
     181             :                                                     std::vector<RankTwoTensor> & r)
     182             : {
     183             :   mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
     184             :   mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
     185             : 
     186     1994724 :   r.resize(0);
     187             :   std::vector<unsigned int> active_surfaces_of_model;
     188             :   std::vector<unsigned int>::iterator active_surface;
     189             :   std::vector<RankTwoTensor> model_r;
     190     5216836 :   for (unsigned model = 0; model < _num_models; ++model)
     191             :   {
     192     3222112 :     activeModelSurfaces(model, active, active_surfaces_of_model);
     193     3222112 :     if (active_surfaces_of_model.size() > 0)
     194             :     {
     195     2602052 :       _f[model]->flowPotentialV(stress, intnl[model], model_r);
     196     2602052 :       for (active_surface = active_surfaces_of_model.begin();
     197     5266936 :            active_surface != active_surfaces_of_model.end();
     198             :            ++active_surface)
     199     2664884 :         r.push_back(model_r[*active_surface]);
     200             :     }
     201             :   }
     202     1994724 : }
     203             : 
     204             : void
     205      685824 : MultiPlasticityRawComponentAssembler::dflowPotential_dstress(
     206             :     const RankTwoTensor & stress,
     207             :     const std::vector<Real> & intnl,
     208             :     const std::vector<bool> & active,
     209             :     std::vector<RankFourTensor> & dr_dstress)
     210             : {
     211             :   mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
     212             :   mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
     213             : 
     214      685824 :   dr_dstress.resize(0);
     215             :   std::vector<unsigned int> active_surfaces_of_model;
     216             :   std::vector<unsigned int>::iterator active_surface;
     217             :   std::vector<RankFourTensor> model_dr_dstress;
     218     1787312 :   for (unsigned model = 0; model < _num_models; ++model)
     219             :   {
     220     1101488 :     activeModelSurfaces(model, active, active_surfaces_of_model);
     221     1101488 :     if (active_surfaces_of_model.size() > 0)
     222             :     {
     223      885908 :       _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
     224      885908 :       for (active_surface = active_surfaces_of_model.begin();
     225     1804696 :            active_surface != active_surfaces_of_model.end();
     226             :            ++active_surface)
     227      918788 :         dr_dstress.push_back(model_dr_dstress[*active_surface]);
     228             :     }
     229             :   }
     230      685824 : }
     231             : 
     232             : void
     233      685824 : MultiPlasticityRawComponentAssembler::dflowPotential_dintnl(const RankTwoTensor & stress,
     234             :                                                             const std::vector<Real> & intnl,
     235             :                                                             const std::vector<bool> & active,
     236             :                                                             std::vector<RankTwoTensor> & dr_dintnl)
     237             : {
     238             :   mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
     239             :   mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
     240             : 
     241      685824 :   dr_dintnl.resize(0);
     242             :   std::vector<unsigned int> active_surfaces_of_model;
     243             :   std::vector<unsigned int>::iterator active_surface;
     244             :   std::vector<RankTwoTensor> model_dr_dintnl;
     245     1787312 :   for (unsigned model = 0; model < _num_models; ++model)
     246             :   {
     247     1101488 :     activeModelSurfaces(model, active, active_surfaces_of_model);
     248     1101488 :     if (active_surfaces_of_model.size() > 0)
     249             :     {
     250      885908 :       _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
     251      885908 :       for (active_surface = active_surfaces_of_model.begin();
     252     1804696 :            active_surface != active_surfaces_of_model.end();
     253             :            ++active_surface)
     254      918788 :         dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
     255             :     }
     256             :   }
     257      685824 : }
     258             : 
     259             : void
     260     1994724 : MultiPlasticityRawComponentAssembler::hardPotential(const RankTwoTensor & stress,
     261             :                                                     const std::vector<Real> & intnl,
     262             :                                                     const std::vector<bool> & active,
     263             :                                                     std::vector<Real> & h)
     264             : {
     265             :   mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
     266             :   mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
     267             : 
     268     1994724 :   h.resize(0);
     269             :   std::vector<unsigned int> active_surfaces_of_model;
     270             :   std::vector<unsigned int>::iterator active_surface;
     271             :   std::vector<Real> model_h;
     272     5216836 :   for (unsigned model = 0; model < _num_models; ++model)
     273             :   {
     274     3222112 :     activeModelSurfaces(model, active, active_surfaces_of_model);
     275     3222112 :     if (active_surfaces_of_model.size() > 0)
     276             :     {
     277     2602052 :       _f[model]->hardPotentialV(stress, intnl[model], model_h);
     278     2602052 :       for (active_surface = active_surfaces_of_model.begin();
     279     5266936 :            active_surface != active_surfaces_of_model.end();
     280             :            ++active_surface)
     281     2664884 :         h.push_back(model_h[*active_surface]);
     282             :     }
     283             :   }
     284     1994724 : }
     285             : 
     286             : void
     287      544304 : MultiPlasticityRawComponentAssembler::dhardPotential_dstress(
     288             :     const RankTwoTensor & stress,
     289             :     const std::vector<Real> & intnl,
     290             :     const std::vector<bool> & active,
     291             :     std::vector<RankTwoTensor> & dh_dstress)
     292             : {
     293             :   mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
     294             :   mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
     295             : 
     296      544304 :   dh_dstress.resize(0);
     297             :   std::vector<unsigned int> active_surfaces_of_model;
     298             :   std::vector<unsigned int>::iterator active_surface;
     299             :   std::vector<RankTwoTensor> model_dh_dstress;
     300     1433424 :   for (unsigned model = 0; model < _num_models; ++model)
     301             :   {
     302      889120 :     activeModelSurfaces(model, active, active_surfaces_of_model);
     303      889120 :     if (active_surfaces_of_model.size() > 0)
     304             :     {
     305      724772 :       _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
     306      724772 :       for (active_surface = active_surfaces_of_model.begin();
     307     1464520 :            active_surface != active_surfaces_of_model.end();
     308             :            ++active_surface)
     309      739748 :         dh_dstress.push_back(model_dh_dstress[*active_surface]);
     310             :     }
     311             :   }
     312      544304 : }
     313             : 
     314             : void
     315      544304 : MultiPlasticityRawComponentAssembler::dhardPotential_dintnl(const RankTwoTensor & stress,
     316             :                                                             const std::vector<Real> & intnl,
     317             :                                                             const std::vector<bool> & active,
     318             :                                                             std::vector<Real> & dh_dintnl)
     319             : {
     320             :   mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
     321             :   mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
     322             : 
     323      544304 :   dh_dintnl.resize(0);
     324             :   std::vector<unsigned int> active_surfaces_of_model;
     325             :   std::vector<unsigned int>::iterator active_surface;
     326             :   std::vector<Real> model_dh_dintnl;
     327     1433424 :   for (unsigned model = 0; model < _num_models; ++model)
     328             :   {
     329      889120 :     activeModelSurfaces(model, active, active_surfaces_of_model);
     330      889120 :     if (active_surfaces_of_model.size() > 0)
     331             :     {
     332      724772 :       _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
     333      724772 :       for (active_surface = active_surfaces_of_model.begin();
     334     1464520 :            active_surface != active_surfaces_of_model.end();
     335             :            ++active_surface)
     336      739748 :         dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
     337             :     }
     338             :   }
     339      544304 : }
     340             : 
     341             : void
     342      142656 : MultiPlasticityRawComponentAssembler::buildActiveConstraints(const std::vector<Real> & f,
     343             :                                                              const RankTwoTensor & stress,
     344             :                                                              const std::vector<Real> & intnl,
     345             :                                                              const RankFourTensor & Eijkl,
     346             :                                                              std::vector<bool> & act)
     347             : {
     348             :   mooseAssert(f.size() == _num_surfaces,
     349             :               "buildActiveConstraints called with f.size = " << f.size() << " while there are "
     350             :                                                              << _num_surfaces << " surfaces");
     351             :   mooseAssert(intnl.size() == _num_models,
     352             :               "buildActiveConstraints called with intnl.size = "
     353             :                   << intnl.size() << " while there are " << _num_models << " models");
     354             : 
     355      142656 :   if (_specialIC == "rock")
     356       16960 :     buildActiveConstraintsRock(f, stress, intnl, Eijkl, act);
     357      125696 :   else if (_specialIC == "joint")
     358       14816 :     buildActiveConstraintsJoint(f, stress, intnl, Eijkl, act);
     359             :   else // no specialIC
     360             :   {
     361      110880 :     act.resize(0);
     362             :     unsigned ind = 0;
     363      273936 :     for (unsigned model = 0; model < _num_models; ++model)
     364             :     {
     365      163056 :       std::vector<Real> model_f(0);
     366      326112 :       for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
     367      163056 :         model_f.push_back(f[ind++]);
     368             :       std::vector<bool> model_act;
     369      163056 :       RankTwoTensor returned_stress;
     370      163056 :       _f[model]->activeConstraints(
     371             :           model_f, stress, intnl[model], Eijkl, model_act, returned_stress);
     372      326112 :       for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
     373      163056 :         act.push_back(model_act[model_surface]);
     374             :     }
     375             :   }
     376      142656 : }
     377             : 
     378             : void
     379       14816 : MultiPlasticityRawComponentAssembler::buildActiveConstraintsJoint(const std::vector<Real> & f,
     380             :                                                                   const RankTwoTensor & stress,
     381             :                                                                   const std::vector<Real> & intnl,
     382             :                                                                   const RankFourTensor & Eijkl,
     383             :                                                                   std::vector<bool> & act)
     384             : {
     385             :   act.assign(2, false);
     386             : 
     387       14816 :   RankTwoTensor returned_stress;
     388             :   std::vector<bool> active_tensile;
     389             :   std::vector<bool> active_shear;
     390             :   std::vector<Real> f_single;
     391             : 
     392             :   // first try tensile alone
     393       14816 :   f_single.assign(1, 0);
     394       14816 :   f_single[0] = f[0];
     395       14816 :   _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
     396       14816 :   _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
     397       14816 :   if (f_single[0] <= _f[1]->_f_tol)
     398             :   {
     399           0 :     act[0] = active_tensile[0];
     400           0 :     return;
     401             :   }
     402             : 
     403             :   // next try shear alone
     404       14816 :   f_single.assign(1, 0);
     405       14816 :   f_single[0] = f[1];
     406       14816 :   _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_shear, returned_stress);
     407       14816 :   _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
     408       14816 :   if (f_single[0] <= _f[0]->_f_tol)
     409             :   {
     410        7424 :     act[1] = active_shear[0];
     411        7424 :     return;
     412             :   }
     413             : 
     414             :   // must be mixed
     415        7392 :   act[0] = act[1] = true;
     416        7392 :   return;
     417             : }
     418             : 
     419             : void
     420       16960 : MultiPlasticityRawComponentAssembler::buildActiveConstraintsRock(const std::vector<Real> & f,
     421             :                                                                  const RankTwoTensor & stress,
     422             :                                                                  const std::vector<Real> & intnl,
     423             :                                                                  const RankFourTensor & Eijkl,
     424             :                                                                  std::vector<bool> & act)
     425             : {
     426             :   act.assign(9, false);
     427             : 
     428       16960 :   RankTwoTensor returned_stress;
     429             :   std::vector<bool> active_tensile;
     430             :   std::vector<bool> active_MC;
     431             :   std::vector<Real> f_single;
     432             : 
     433             :   // first try tensile alone
     434       16960 :   f_single.assign(3, 0);
     435       16960 :   f_single[0] = f[0];
     436       16960 :   f_single[1] = f[1];
     437       16960 :   f_single[2] = f[2];
     438       16960 :   _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
     439       16960 :   _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
     440       16960 :   if (f_single[0] <= _f[1]->_f_tol && f_single[1] <= _f[1]->_f_tol &&
     441       13104 :       f_single[2] <= _f[1]->_f_tol && f_single[3] <= _f[1]->_f_tol &&
     442       21632 :       f_single[4] <= _f[1]->_f_tol && f_single[5] <= _f[1]->_f_tol)
     443             :   {
     444        4672 :     act[0] = active_tensile[0];
     445        4672 :     act[1] = active_tensile[1];
     446        4672 :     act[2] = active_tensile[2];
     447        4672 :     return;
     448             :   }
     449             : 
     450             :   // next try MC alone
     451       12288 :   f_single.assign(6, 0);
     452       12288 :   f_single[0] = f[3];
     453       12288 :   f_single[1] = f[4];
     454       12288 :   f_single[2] = f[5];
     455       12288 :   f_single[3] = f[6];
     456       12288 :   f_single[4] = f[7];
     457       12288 :   f_single[5] = f[8];
     458       12288 :   _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_MC, returned_stress);
     459       12288 :   _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
     460       12288 :   if (f_single[0] <= _f[0]->_f_tol && f_single[1] <= _f[0]->_f_tol && f_single[2] <= _f[0]->_f_tol)
     461             :   {
     462        6944 :     act[3] = active_MC[0];
     463        6944 :     act[4] = active_MC[1];
     464        6944 :     act[5] = active_MC[2];
     465        6944 :     act[6] = active_MC[3];
     466        6944 :     act[7] = active_MC[4];
     467        6944 :     act[8] = active_MC[5];
     468        6944 :     return;
     469             :   }
     470             : 
     471             :   // must be a mix.
     472             :   // The possibilities are enumerated below.
     473             : 
     474             :   // tensile=edge, MC=tip (two possibilities)
     475        5344 :   if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
     476         400 :       active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
     477        5744 :       active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
     478             :   {
     479         400 :     act[1] = act[2] = act[6] = true;
     480             :     act[4] = true;
     481         400 :     return;
     482             :   }
     483        4944 :   if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
     484         784 :       active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
     485        5728 :       active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
     486             :   {
     487           0 :     act[1] = act[2] = act[6] = true; // i don't think act[4] is necessary, is it?!
     488           0 :     return;
     489             :   }
     490             : 
     491             :   // tensile = edge, MC=edge (two possibilities)
     492        4944 :   if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
     493         784 :       active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
     494        5728 :       active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
     495             :   {
     496         784 :     act[1] = act[2] = act[4] = act[6] = true;
     497         784 :     return;
     498             :   }
     499        4160 :   if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
     500           0 :       active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
     501        4160 :       active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
     502             :   {
     503           0 :     act[1] = act[2] = act[4] = act[6] = true;
     504           0 :     return;
     505             :   }
     506             : 
     507             :   // tensile = edge, MC=face
     508        4160 :   if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
     509           0 :       active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
     510        4160 :       active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
     511             :   {
     512           0 :     act[1] = act[2] = act[6] = true;
     513           0 :     return;
     514             :   }
     515             : 
     516             :   // tensile = face, MC=tip (two possibilities)
     517        4160 :   if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
     518        1280 :       active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
     519        5440 :       active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
     520             :   {
     521        1280 :     act[2] = act[6] = true;
     522             :     act[4] = true;
     523             :     act[8] = true;
     524        1280 :     return;
     525             :   }
     526        2880 :   if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
     527        2880 :       active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
     528        3408 :       active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
     529             :   {
     530         128 :     act[2] = act[6] = true;
     531             :     act[8] = true;
     532         128 :     return;
     533             :   }
     534             : 
     535             :   // tensile = face, MC=face
     536        2752 :   if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
     537        2752 :       active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
     538        5104 :       active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
     539             :   {
     540        1232 :     act[1] = act[2] = act[6] = true;
     541        1232 :     return;
     542             :   }
     543             : 
     544             :   // tensile = face, MC=edge (two possibilites).
     545             :   act[2] = true; // tensile face
     546        1520 :   act[3] = active_MC[0];
     547        1520 :   act[4] = active_MC[1];
     548        1520 :   act[5] = active_MC[2];
     549        1520 :   act[6] = active_MC[3];
     550        1520 :   act[7] = active_MC[4];
     551        1520 :   act[8] = active_MC[5];
     552        1520 :   return;
     553             : }
     554             : 
     555             : /**
     556             :  * Performs a returnMap for each plastic model.
     557             :  *
     558             :  * If all models actually signal "elastic" by
     559             :  * returning true from their returnMap, and
     560             :  * by returning model_plastically_active=0, then
     561             :  *   yf will contain the yield function values
     562             :  *   num_successful_plastic_returns will be zero
     563             :  *   intnl = intnl_old
     564             :  *   delta_dp will be unchanged from its input value
     565             :  *   stress will be set to trial_stress
     566             :  *   pm will be zero
     567             :  *   cumulative_pm will be unchanged
     568             :  *   return value will be true
     569             :  *   num_successful_plastic_returns = 0
     570             :  *
     571             :  * If only one model signals "plastically active"
     572             :  * by returning true from its returnMap,
     573             :  * and by returning model_plastically_active=1, then
     574             :  *   yf will contain the yield function values
     575             :  *   num_successful_plastic_returns will be one
     576             :  *   intnl will be set by the returnMap algorithm
     577             :  *   delta_dp will be set by the returnMap algorithm
     578             :  *   stress will be set by the returnMap algorithm
     579             :  *   pm will be nonzero for the single model, and zero for other models
     580             :  *   cumulative_pm will be updated
     581             :  *   return value will be true
     582             :  *   num_successful_plastic_returns = 1
     583             :  *
     584             :  * If >1 model signals "plastically active"
     585             :  * or if >=1 model's returnMap fails, then
     586             :  *   yf will contain the yield function values
     587             :  *   num_successful_plastic_returns will be set appropriately
     588             :  *   intnl = intnl_old
     589             :  *   delta_dp will be unchanged from its input value
     590             :  *   stress will be set to trial_stress
     591             :  *   pm will be zero
     592             :  *   cumulative_pm will be unchanged
     593             :  *   return value will be true if all returnMap functions returned true, otherwise it will be false
     594             :  *   num_successful_plastic_returns is set appropriately
     595             :  */
     596             : bool
     597     6801584 : MultiPlasticityRawComponentAssembler::returnMapAll(const RankTwoTensor & trial_stress,
     598             :                                                    const std::vector<Real> & intnl_old,
     599             :                                                    const RankFourTensor & E_ijkl,
     600             :                                                    Real ep_plastic_tolerance,
     601             :                                                    RankTwoTensor & stress,
     602             :                                                    std::vector<Real> & intnl,
     603             :                                                    std::vector<Real> & pm,
     604             :                                                    std::vector<Real> & cumulative_pm,
     605             :                                                    RankTwoTensor & delta_dp,
     606             :                                                    std::vector<Real> & yf,
     607             :                                                    unsigned & num_successful_plastic_returns,
     608             :                                                    unsigned & custom_model)
     609             : {
     610             :   mooseAssert(intnl_old.size() == _num_models,
     611             :               "returnMapAll: Incorrect size of internal parameters");
     612             :   mooseAssert(intnl.size() == _num_models, "returnMapAll: Incorrect size of internal parameters");
     613             :   mooseAssert(pm.size() == _num_surfaces, "returnMapAll: Incorrect size of pm");
     614             : 
     615     6801584 :   num_successful_plastic_returns = 0;
     616     6801584 :   yf.resize(0);
     617     6801584 :   pm.assign(_num_surfaces, 0.0);
     618             : 
     619     6801584 :   RankTwoTensor returned_stress; // each model will give a returned_stress.  if only one model is
     620             :                                  // plastically active, i set stress=returned_stress, so as to
     621             :                                  // record this returned value
     622             :   std::vector<Real> model_f;
     623     6801584 :   RankTwoTensor model_delta_dp;
     624             :   std::vector<Real> model_pm;
     625             :   bool trial_stress_inadmissible;
     626             :   bool successful_return = true;
     627             :   unsigned the_single_plastic_model = 0;
     628             :   bool using_custom_return_map = true;
     629             : 
     630             :   // run through all the plastic models, performing their
     631             :   // returnMap algorithms.
     632             :   // If one finds (trial_stress, intnl) inadmissible and
     633             :   // successfully returns, break from the loop to evaluate
     634             :   // all the others at that returned stress
     635     7254180 :   for (unsigned model = 0; model < _num_models; ++model)
     636             :   {
     637     6931456 :     if (using_custom_return_map)
     638             :     {
     639     6802816 :       model_pm.assign(_f[model]->numberSurfaces(), 0.0);
     640     6802816 :       bool model_returned = _f[model]->returnMap(trial_stress,
     641             :                                                  intnl_old[model],
     642             :                                                  E_ijkl,
     643             :                                                  ep_plastic_tolerance,
     644             :                                                  returned_stress,
     645             :                                                  intnl[model],
     646             :                                                  model_pm,
     647             :                                                  model_delta_dp,
     648             :                                                  model_f,
     649             :                                                  trial_stress_inadmissible);
     650     6802816 :       if (!trial_stress_inadmissible)
     651             :       {
     652             :         // in the elastic zone: record the yield-function values (returned_stress, intnl, model_pm
     653             :         // and model_delta_dp are undefined)
     654      109224 :         for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
     655             :              ++model_surface)
     656       62180 :           yf.push_back(model_f[model_surface]);
     657             :       }
     658     6755772 :       else if (trial_stress_inadmissible && !model_returned)
     659             :       {
     660             :         // in the plastic zone, and the customized returnMap failed
     661             :         // for some reason (or wasn't implemented).  The coder
     662             :         // should have correctly returned model_f(trial_stress, intnl_old)
     663             :         // so record them
     664             :         // (returned_stress, intnl, model_pm and model_delta_dp are undefined)
     665      619584 :         for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
     666             :              ++model_surface)
     667      342672 :           yf.push_back(model_f[model_surface]);
     668             :         // now there's almost zero point in using the custom
     669             :         // returnMap functions
     670             :         using_custom_return_map = false;
     671             :         successful_return = false;
     672             :       }
     673             :       else
     674             :       {
     675             :         // in the plastic zone, and the customized returnMap
     676             :         // succeeded.
     677             :         // record the first returned_stress and delta_dp if everything is going OK
     678             :         // as they could be the actual answer
     679             :         if (trial_stress_inadmissible)
     680     6478860 :           num_successful_plastic_returns++;
     681             :         the_single_plastic_model = model;
     682     6478860 :         stress = returned_stress;
     683             :         // note that i break here, and don't push_back
     684             :         // model_f to yf.  So now yf contains only the values of
     685             :         // model_f from previous models to the_single_plastic_model
     686             :         // also i don't set delta_dp = model_delta_dp yet, because
     687             :         // i might find problems later on
     688             :         // also, don't increment cumulative_pm for the same reason
     689             : 
     690     6478860 :         break;
     691             :       }
     692             :     }
     693             :     else
     694             :     {
     695             :       // not using custom returnMap functions because one
     696             :       // has already failed and that one said trial_stress
     697             :       // was inadmissible.  So now calculate the yield functions
     698             :       // at the trial stress
     699      128640 :       _f[model]->yieldFunctionV(trial_stress, intnl_old[model], model_f);
     700      398880 :       for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
     701      270240 :         yf.push_back(model_f[model_surface]);
     702             :     }
     703             :   }
     704             : 
     705     6801584 :   if (num_successful_plastic_returns == 0)
     706             :   {
     707             :     // here either all the models were elastic (successful_return=true),
     708             :     // or some were plastic and either the customized returnMap failed
     709             :     // or wasn't implemented (successful_return=false).
     710             :     // In either case, have to set the following:
     711      322724 :     stress = trial_stress;
     712      775320 :     for (unsigned model = 0; model < _num_models; ++model)
     713      452596 :       intnl[model] = intnl_old[model];
     714             :     return successful_return;
     715             :   }
     716             : 
     717             :   // Now we know that num_successful_plastic_returns == 1 and all the other
     718             :   // models (with model number < the_single_plastic_model) must have been
     719             :   // admissible at (trial_stress, intnl).  However, all models might
     720             :   // not be admissible at (trial_stress, intnl), so must check that
     721     6478860 :   std::vector<Real> yf_at_returned_stress(0);
     722             :   bool all_admissible = true;
     723    12959480 :   for (unsigned model = 0; model < _num_models; ++model)
     724             :   {
     725     6484396 :     if (model == the_single_plastic_model)
     726             :     {
     727             :       // no need to spend time calculating the yield function: we know its admissible
     728    13027928 :       for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
     729     6549068 :         yf_at_returned_stress.push_back(model_f[model_surface]);
     730     6478860 :       continue;
     731     6478860 :     }
     732        5536 :     _f[model]->yieldFunctionV(stress, intnl_old[model], model_f);
     733       38752 :     for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
     734             :     {
     735       33216 :       if (model_f[model_surface] > _f[model]->_f_tol)
     736             :         // bummer, this model is not admissible at the returned_stress
     737             :         all_admissible = false;
     738       33216 :       yf_at_returned_stress.push_back(model_f[model_surface]);
     739             :     }
     740        5536 :     if (!all_admissible)
     741             :       // no point in continuing computing yield functions
     742             :       break;
     743             :   }
     744             : 
     745     6478860 :   if (!all_admissible)
     746             :   {
     747             :     // we tried using the returned value of stress predicted by
     748             :     // the_single_plastic_model, but it wasn't admissible according
     749             :     // to other plastic models.  We need to set:
     750        3776 :     stress = trial_stress;
     751       11328 :     for (unsigned model = 0; model < _num_models; ++model)
     752        7552 :       intnl[model] = intnl_old[model];
     753             :     // and calculate the remainder of the yield functions at trial_stress
     754       11328 :     for (unsigned model = the_single_plastic_model; model < _num_models; ++model)
     755             :     {
     756        7552 :       _f[model]->yieldFunctionV(trial_stress, intnl[model], model_f);
     757       41536 :       for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
     758       33984 :         yf.push_back(model_f[model_surface]);
     759             :     }
     760        3776 :     num_successful_plastic_returns = 0;
     761        3776 :     return false;
     762             :   }
     763             : 
     764             :   // So the customized returnMap algorithm can provide a returned
     765             :   // (stress, intnl) configuration, and that is admissible according
     766             :   // to all plastic models
     767     6475084 :   yf.resize(0);
     768    13023384 :   for (unsigned surface = 0; surface < yf_at_returned_stress.size(); ++surface)
     769     6548300 :     yf.push_back(yf_at_returned_stress[surface]);
     770     6475084 :   delta_dp = model_delta_dp;
     771    13012824 :   for (unsigned model_surface = 0; model_surface < _f[the_single_plastic_model]->numberSurfaces();
     772             :        ++model_surface)
     773             :   {
     774     6537740 :     cumulative_pm[_surfaces_given_model[the_single_plastic_model][model_surface]] +=
     775     6537740 :         model_pm[model_surface];
     776     6537740 :     pm[_surfaces_given_model[the_single_plastic_model][model_surface]] = model_pm[model_surface];
     777             :   }
     778     6475084 :   custom_model = the_single_plastic_model;
     779     6475084 :   return true;
     780             : }
     781             : 
     782             : unsigned int
     783     6920448 : MultiPlasticityRawComponentAssembler::modelNumber(unsigned int surface)
     784             : {
     785     6920448 :   return _model_given_surface[surface];
     786             : }
     787             : 
     788             : bool
     789     5557068 : MultiPlasticityRawComponentAssembler::anyActiveSurfaces(int model, const std::vector<bool> & active)
     790             : {
     791     7353576 :   for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
     792     6131308 :     if (active[_surfaces_given_model[model][model_surface]])
     793             :       return true;
     794             :   return false;
     795             : }
     796             : 
     797             : void
     798     2120624 : MultiPlasticityRawComponentAssembler::activeSurfaces(int model,
     799             :                                                      const std::vector<bool> & active,
     800             :                                                      std::vector<unsigned int> & active_surfaces)
     801             : {
     802     2120624 :   active_surfaces.resize(0);
     803     4536256 :   for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
     804     2415632 :     if (active[_surfaces_given_model[model][model_surface]])
     805     1746096 :       active_surfaces.push_back(_surfaces_given_model[model][model_surface]);
     806     2120624 : }
     807             : 
     808             : void
     809    15185024 : MultiPlasticityRawComponentAssembler::activeModelSurfaces(
     810             :     int model,
     811             :     const std::vector<bool> & active,
     812             :     std::vector<unsigned int> & active_surfaces_of_model)
     813             : {
     814    15185024 :   active_surfaces_of_model.resize(0);
     815    33702576 :   for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
     816    18517552 :     if (active[_surfaces_given_model[model][model_surface]])
     817    12746984 :       active_surfaces_of_model.push_back(model_surface);
     818    15185024 : }

Generated by: LCOV version 1.14