LCOV - code coverage report
Current view: top level - include/utils - MultiPlasticityRawComponentAssembler.h (source / functions) Hit Total Coverage
Test: idaholab/moose solid_mechanics: f45d79 Lines: 1 1 100.0 %
Date: 2025-07-25 05:00:39 Functions: 1 2 50.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             : #pragma once
      11             : /* MOOSE - Multiphysics Object Oriented Simulation Environment  */
      12             : /*                                                              */
      13             : /*          All contents are licensed under LGPL V2.1           */
      14             : /*             See LICENSE for full restrictions                */
      15             : /****************************************************************/
      16             : 
      17             : #include "SolidMechanicsPlasticModel.h"
      18             : #include "UserObjectInterface.h"
      19             : 
      20             : /**
      21             :  * MultiPlasticityRawComponentAssembler holds and computes yield functions,
      22             :  * flow directions, etc, for use in FiniteStrainMultiPlasticity
      23             :  *
      24             :  * Here there are a number of numbering systems.
      25             :  *
      26             :  * There are _num_models of plastic models, read from the input file.
      27             :  *    Each of these models can, in principal, contain multiple
      28             :  *    'internal surfaces'.
      29             :  *    Models are numbered from 0 to _num_models - 1.
      30             :  *
      31             :  * There are _num_surfaces surfaces.  This
      32             :  *    = sum_{plastic_models} (numberSurfaces in model)
      33             :  *    Evidently _num_surface >= _num_models
      34             :  *    Surfaces are numbered from 0 to _num_surfaces - 1.
      35             :  *
      36             :  * The std::vectors _model_given_surface, _model_surface_given_surface
      37             :  * and _surfaces_given_model allow translation between these
      38             :  */
      39             : class MultiPlasticityRawComponentAssembler : public UserObjectInterface
      40             : {
      41             : public:
      42             :   static InputParameters validParams();
      43             : 
      44             :   MultiPlasticityRawComponentAssembler(const MooseObject * moose_object);
      45             : 
      46        6192 :   virtual ~MultiPlasticityRawComponentAssembler() {}
      47             : 
      48             : protected:
      49             :   const InputParameters & _params;
      50             : 
      51             :   /// Number of plastic models for this material
      52             :   unsigned int _num_models;
      53             : 
      54             :   /**
      55             :    * Number of surfaces within the plastic models.
      56             :    * For many situations this will be = _num_models
      57             :    * since each model will contain just one surface.
      58             :    * More generally it is >= _num_models.  For instance,
      59             :    * Mohr-Coulomb is a single model with 6 surfaces
      60             :    */
      61             :   unsigned int _num_surfaces;
      62             : 
      63             :   /// _surfaces_given_model[model_number] = vector of surface numbers for this model
      64             :   std::vector<std::vector<unsigned int>> _surfaces_given_model;
      65             : 
      66             :   /// Allows initial set of active constraints to be chosen optimally
      67             :   MooseEnum _specialIC;
      68             : 
      69             :   /// User objects that define the yield functions, flow potentials, etc
      70             :   std::vector<const SolidMechanicsPlasticModel *> _f;
      71             : 
      72             :   /**
      73             :    * The active yield function(s)
      74             :    * @param stress the stress at which to calculate the yield function
      75             :    * @param intnl vector of internal parameters
      76             :    * @param active set of active constraints - only the active yield functions are put into "f"
      77             :    * @param[out] f the yield function (or functions in the case of multisurface plasticity)
      78             :    */
      79             :   virtual void yieldFunction(const RankTwoTensor & stress,
      80             :                              const std::vector<Real> & intnl,
      81             :                              const std::vector<bool> & active,
      82             :                              std::vector<Real> & f);
      83             : 
      84             :   /**
      85             :    * The derivative of the active yield function(s) with respect to stress
      86             :    * @param stress the stress at which to calculate the yield function
      87             :    * @param intnl vector of internal parameters
      88             :    * @param active set of active constraints - only the active derivatives are put into "df_dstress"
      89             :    * @param[out] df_dstress the derivative (or derivatives in the case of multisurface plasticity).
      90             :    * df_dstress[alpha](i, j) = dyieldFunction[alpha]/dstress(i, j)
      91             :    */
      92             :   virtual void dyieldFunction_dstress(const RankTwoTensor & stress,
      93             :                                       const std::vector<Real> & intnl,
      94             :                                       const std::vector<bool> & active,
      95             :                                       std::vector<RankTwoTensor> & df_dstress);
      96             : 
      97             :   /**
      98             :    * The derivative of active yield function(s) with respect to their internal parameters (the user
      99             :    * objects assume there is exactly one internal param per yield function)
     100             :    * @param stress the stress at which to calculate the yield function
     101             :    * @param intnl vector of internal parameters
     102             :    * @param active set of active constraints - only the active derivatives are put into "df_dintnl"
     103             :    * @param[out] df_dintnl the derivatives.  df_dstress[alpha] = dyieldFunction[alpha]/dintnl[alpha]
     104             :    */
     105             :   virtual void dyieldFunction_dintnl(const RankTwoTensor & stress,
     106             :                                      const std::vector<Real> & intnl,
     107             :                                      const std::vector<bool> & active,
     108             :                                      std::vector<Real> & df_dintnl);
     109             : 
     110             :   /**
     111             :    * The active flow potential(s) - one for each yield function
     112             :    * @param stress the stress at which to calculate the flow potential
     113             :    * @param intnl vector of internal parameters
     114             :    * @param active set of active constraints - only the active flow potentials are put into "r"
     115             :    * @param[out] r the flow potential (flow potentials in the multi-surface case)
     116             :    */
     117             :   virtual void flowPotential(const RankTwoTensor & stress,
     118             :                              const std::vector<Real> & intnl,
     119             :                              const std::vector<bool> & active,
     120             :                              std::vector<RankTwoTensor> & r);
     121             : 
     122             :   /**
     123             :    * The derivative of the active flow potential(s) with respect to stress
     124             :    * @param stress the stress at which to calculate the flow potential
     125             :    * @param intnl vector of internal parameters
     126             :    * @param active set of active constraints - only the active derivatives are put into "dr_dstress"
     127             :    * @param[out] dr_dstress the derivative.  dr_dstress[alpha](i, j, k, l) = dr[alpha](i,
     128             :    * j)/dstress(k, l)
     129             :    */
     130             :   virtual void dflowPotential_dstress(const RankTwoTensor & stress,
     131             :                                       const std::vector<Real> & intnl,
     132             :                                       const std::vector<bool> & active,
     133             :                                       std::vector<RankFourTensor> & dr_dstress);
     134             : 
     135             :   /**
     136             :    * The derivative of the active flow potentials with respect to the active internal parameters
     137             :    * The UserObjects explicitly assume that r[alpha] is only dependent on intnl[alpha]
     138             :    * @param stress the stress at which to calculate the flow potential
     139             :    * @param intnl vector of internal parameters
     140             :    * @param active set of active constraints - only the active derivatives are put into "dr_dintnl"
     141             :    * @param[out] dr_dintnl the derivatives.  dr_dintnl[alpha](i, j) = dr[alpha](i, j)/dintnl[alpha]
     142             :    */
     143             :   virtual void dflowPotential_dintnl(const RankTwoTensor & stress,
     144             :                                      const std::vector<Real> & intnl,
     145             :                                      const std::vector<bool> & active,
     146             :                                      std::vector<RankTwoTensor> & dr_dintnl);
     147             : 
     148             :   /**
     149             :    * The active hardening potentials (one for each internal parameter and for each yield function)
     150             :    * by assumption in the Userobjects, the h[a][alpha] is nonzero only if the surface alpha is part
     151             :    * of model a, so we only calculate those here
     152             :    * @param stress the stress at which to calculate the hardening potential
     153             :    * @param intnl vector of internal parameters
     154             :    * @param active set of active constraints - only the active hardening potentials are put into "h"
     155             :    * @param[out] h the hardening potentials.  h[alpha] = hardening potential for yield fcn alpha
     156             :    * (and, by the above assumption we know which hardening parameter, a, this belongs to)
     157             :    */
     158             :   virtual void hardPotential(const RankTwoTensor & stress,
     159             :                              const std::vector<Real> & intnl,
     160             :                              const std::vector<bool> & active,
     161             :                              std::vector<Real> & h);
     162             : 
     163             :   /**
     164             :    * The derivative of the active hardening potentials with respect to stress
     165             :    * By assumption in the Userobjects, the h[a][alpha] is nonzero only for a = alpha, so we only
     166             :    * calculate those here
     167             :    * @param stress the stress at which to calculate the hardening potentials
     168             :    * @param intnl vector of internal parameters
     169             :    * @param active set of active constraints - only the active derivatives are put into "dh_dstress"
     170             :    * @param[out] dh_dstress the derivative.  dh_dstress[a](i, j) = dh[a]/dstress(k, l)
     171             :    */
     172             :   virtual void dhardPotential_dstress(const RankTwoTensor & stress,
     173             :                                       const std::vector<Real> & intnl,
     174             :                                       const std::vector<bool> & active,
     175             :                                       std::vector<RankTwoTensor> & dh_dstress);
     176             : 
     177             :   /**
     178             :    * The derivative of the active hardening potentials with respect to the active internal
     179             :    * parameters
     180             :    * @param stress the stress at which to calculate the hardening potentials
     181             :    * @param intnl vector of internal parameters
     182             :    * @param active set of active constraints - only the active derivatives are put into "dh_dintnl"
     183             :    * @param[out] dh_dintnl the derivatives.  dh_dintnl[a][alpha][b] = dh[a][alpha]/dintnl[b].  Note
     184             :    * that the userobjects assume that there is exactly one internal parameter per yield function, so
     185             :    * the derivative is only nonzero for a=alpha=b, so that is all we calculate
     186             :    */
     187             :   virtual void dhardPotential_dintnl(const RankTwoTensor & stress,
     188             :                                      const std::vector<Real> & intnl,
     189             :                                      const std::vector<bool> & active,
     190             :                                      std::vector<Real> & dh_dintnl);
     191             : 
     192             :   /**
     193             :    * Constructs a set of active constraints, given the yield functions, f.
     194             :    * This uses SolidMechanicsPlasticModel::activeConstraints to identify the active
     195             :    * constraints for each model.
     196             :    * @param f yield functions (should be _num_surfaces of these)
     197             :    * @param stress stress tensor
     198             :    * @param intnl internal parameters
     199             :    * @param Eijkl elasticity tensor (stress = Eijkl*strain)
     200             :    * @param[out] act the set of active constraints (will be resized to _num_surfaces)
     201             :    */
     202             :   virtual void buildActiveConstraints(const std::vector<Real> & f,
     203             :                                       const RankTwoTensor & stress,
     204             :                                       const std::vector<Real> & intnl,
     205             :                                       const RankFourTensor & Eijkl,
     206             :                                       std::vector<bool> & act);
     207             : 
     208             :   /// returns the model number, given the surface number
     209             :   unsigned int modelNumber(unsigned int surface);
     210             : 
     211             :   /// returns true if any internal surfaces of the given model are active according to 'active'
     212             :   bool anyActiveSurfaces(int model, const std::vector<bool> & active);
     213             : 
     214             :   /**
     215             :    * Returns the internal surface number(s) of the active surfaces of the given model
     216             :    * This may be of size=0 if there are no active surfaces of the given model
     217             :    * @param model the model number
     218             :    * @param active array with entries being 'true' if the surface is active
     219             :    * @param[out] active_surfaces_of_model the output
     220             :    */
     221             :   void activeModelSurfaces(int model,
     222             :                            const std::vector<bool> & active,
     223             :                            std::vector<unsigned int> & active_surfaces_of_model);
     224             : 
     225             :   /**
     226             :    * Returns the external surface number(s) of the active surfaces of the given model
     227             :    * This may be of size=0 if there are no active surfaces of the given model
     228             :    * @param model the model number
     229             :    * @param active array with entries being 'true' if the surface is active
     230             :    * @param[out] active_surfaces the output
     231             :    */
     232             :   void activeSurfaces(int model,
     233             :                       const std::vector<bool> & active,
     234             :                       std::vector<unsigned int> & active_surfaces);
     235             : 
     236             :   /**
     237             :    * Performs a returnMap for each plastic model using
     238             :    * their inbuilt returnMap functions.  This may be used
     239             :    * to quickly ascertain whether a (trial_stress, intnl_old) configuration
     240             :    * is admissible, or whether a single model's customized returnMap
     241             :    * function can provide a solution to the return-map problem,
     242             :    * or whether a full Newton-Raphson approach such as implemented
     243             :    * in ComputeMultiPlasticityStress is needed.
     244             :    *
     245             :    * There are three cases mentioned below:
     246             :    * (A) The (trial_stress, intnl_old) configuration is admissible
     247             :    *     according to all plastic models
     248             :    * (B) The (trial_stress, intnl_old) configuration is inadmissible
     249             :    *     to exactly one plastic model, and that model can successfully
     250             :    *     use its customized returnMap function to provide a returned
     251             :    *     (stress, intnl) configuration, and that configuration is
     252             :    *     admissible according to all plastic models
     253             :    * (C) All other cases.  This includes customized returnMap
     254             :    *     functions failing, or more than one plastic_model being
     255             :    *     inadmissible, etc
     256             :    *
     257             :    * @param trial_stress the trial stress
     258             :    * @param intnl_old the old values of the internal parameters
     259             :    * @param E_ijkl the elasticity tensor
     260             :    * @param ep_plastic_tolerance the tolerance on the plastic strain
     261             :    * @param[out] stress is set to trial_stress in case (A) or (C), and the returned value of stress
     262             :    * in case (B).
     263             :    * @param[out] intnl is set to intnl_old in case (A) or (C), and the returned value of intnl in
     264             :    * case (B)
     265             :    * @param[out] pm  Zero in case (A) or (C), otherwise the plastic multipliers needed to bring
     266             :    * about the returnMap in case (B)
     267             :    * @param[in/out] cumulative_pm   cumulative plastic multipliers, updated in case (B), otherwise
     268             :    * left untouched
     269             :    * @param[out] delta_dp is unchanged in case (A) or (C), and is set to the change in plastic
     270             :    * strain in case(B)
     271             :    * @param[out] yf will contain the yield function values at (stress, intnl)
     272             :    * @param[out] num_successful_plastic_returns will be 0 for (A) and (C), and 1 for (B)
     273             :    * @return true in case (A) and (B), and false in case (C)
     274             :    */
     275             :   bool returnMapAll(const RankTwoTensor & trial_stress,
     276             :                     const std::vector<Real> & intnl_old,
     277             :                     const RankFourTensor & E_ijkl,
     278             :                     Real ep_plastic_tolerance,
     279             :                     RankTwoTensor & stress,
     280             :                     std::vector<Real> & intnl,
     281             :                     std::vector<Real> & pm,
     282             :                     std::vector<Real> & cumulative_pm,
     283             :                     RankTwoTensor & delta_dp,
     284             :                     std::vector<Real> & yf,
     285             :                     unsigned & num_successful_plastic_returns,
     286             :                     unsigned & custom_model);
     287             : 
     288             : private:
     289             :   /// given a surface number, this returns the model number
     290             :   std::vector<unsigned int> _model_given_surface;
     291             : 
     292             :   /// given a surface number, this returns the corresponding-model's internal surface number
     293             :   std::vector<unsigned int> _model_surface_given_surface;
     294             : 
     295             :   /**
     296             :    * "Rock" version
     297             :    * Constructs a set of active constraints, given the yield functions, f.
     298             :    * This uses SolidMechanicsPlasticModel::activeConstraints to identify the active
     299             :    * constraints for each model.
     300             :    * @param f yield functions (should be _num_surfaces of these)
     301             :    * @param stress stress tensor
     302             :    * @param intnl internal parameters
     303             :    * @param Eijkl elasticity tensor (stress = Eijkl*strain)
     304             :    * @param[out] act the set of active constraints (will be resized to _num_surfaces)
     305             :    */
     306             :   void buildActiveConstraintsRock(const std::vector<Real> & f,
     307             :                                   const RankTwoTensor & stress,
     308             :                                   const std::vector<Real> & intnl,
     309             :                                   const RankFourTensor & Eijkl,
     310             :                                   std::vector<bool> & act);
     311             : 
     312             :   /**
     313             :    * "Joint" version
     314             :    * Constructs a set of active constraints, given the yield functions, f.
     315             :    * This uses SolidMechanicsPlasticModel::activeConstraints to identify the active
     316             :    * constraints for each model.
     317             :    * @param f yield functions (should be _num_surfaces of these)
     318             :    * @param stress stress tensor
     319             :    * @param intnl internal parameters
     320             :    * @param Eijkl elasticity tensor (stress = Eijkl*strain)
     321             :    * @param[out] act the set of active constraints (will be resized to _num_surfaces)
     322             :    */
     323             :   void buildActiveConstraintsJoint(const std::vector<Real> & f,
     324             :                                    const RankTwoTensor & stress,
     325             :                                    const std::vector<Real> & intnl,
     326             :                                    const RankFourTensor & Eijkl,
     327             :                                    std::vector<bool> & act);
     328             : };

Generated by: LCOV version 1.14