SmallStrainPlasticity.hpp

Problem formulation

Theory

Introduction

In the following document, we define functions, \(\psi(\varepsilon^e,\alpha)\), \(f(\sigma,\mathbf{q})\) and \(\Psi(\sigma,\mathbf{q})\) which are free Helmholtz energy, yield function and flow potential, respectively. Note that free Helmholtz energy is defined by stains and kinematic internal state variables, whereas yield function and flow potential is defined by fluxes (stress tensor and fluxes work conjugated to internal kinematic state variables). In this implantation, we assume that all functions are convex. Also, is assumed that energy potential and flow potential is positive and takes zero value at zero strain and zero internal variable states. It can be shown that for such defined problem, the first law of thermodynamics and non-negative dissipation of energy is automatically satisfied, see [19] [45]. Moreover, in the following implementation is assumed that all functions are smooth enough that appropriate derivatives could be calculated.

In the following implementation, the backward-Euler difference scheme is applied to solve set nonlinear system of equations with unilateral constrains. Admissible stress state is enforced by Lagrange multiplier method and (Khun-Tucker) loading/unloading conditions, see [45]. That lead to well known Closet point projection algorithm.

Following [45] can be shown that plastic loading unloading conditions can be exclusively characterised in terms of trial state

\[ \varepsilon^{e,\textrm{trial}}_{n+1} := \varepsilon_{n+1} - \varepsilon^p_n \]

\[ \alpha^\textrm{trial}_{n+1} := \alpha_n \]

\[ \sigma^\textrm{trial}_{n+1} := \frac{\partial \psi}{\partial \varepsilon^e}(\varepsilon^{e,\textrm{trial}}_{n+1},\alpha_n) \]

\[ \mathbf{q}^\textrm{trial}_{n+1} := \frac{\partial \psi}{\partial \alpha}(\varepsilon^{e,\textrm{trial}}_{n+1},\alpha_n) \]

\[ f^\textrm{trial}_{n+1} := f(\sigma^\textrm{trial}_{n+1},\mathbf{q}_n) \]

then loading/unloading conditions are given

\[ \begin{split} f^\textrm{trial}_{n+1} < 0,\quad\textrm{elastic step}\,\to\,\Delta\gamma_{n+1} = 0 \\ f^\textrm{trial}_{n+1} \geq 0,\quad\textrm{plastic step}\,\to\,\Delta\gamma_{n+1} > 0 \end{split} \]

If first condition is true \(f^\textrm{trial}_{n+1} < 0\), trial stress is solution at integration point. If second condition is true, plastic strains and values of internal state variables are calculated using Closet point projection algorithm.

Closet point projection algorithm

Additive strain decomposition

\[ \varepsilon^e_{n+1} := \varepsilon_{n+1} - \varepsilon^p_{n+1} \]

Constitutive relations

\[ \sigma_{n+1} := \frac{\partial \psi}{\partial \varepsilon^e}(\varepsilon^e_{n+1},\alpha_{n+1}) \]

\[ \mathbf{q}_{n+1} := \frac{\partial \psi}{\partial \alpha}(\varepsilon^e_{n+1},\alpha_{n+1}) \]

Constrains

\[ f_{n+1}:=f(\sigma_{n+1},\mathbf{q}_{n+1}) \]

Evolution functions

Flow vector

\[ \mathbf{n}_{n+1}:=\frac{\partial \Psi}{\partial \sigma}(\sigma_{n+1},\mathbf{q}_{n+1}) \]

Hardening law:

\[ \mathbf{h}_{n+1}:=-\frac{\partial \Psi}{\partial \mathbf{q}}(\sigma_{n+1},\mathbf{q}_{n+1}) \]

Residual

\[ \mathbf{R}_{n+1}:= \left\{ \begin{array}{c} -\varepsilon^p_{n+1}+\varepsilon^p_n \\ -\alpha_{n+1}+\alpha_n \\ f_{n+1} \end{array} \right\} + \Delta\gamma_{n+1} \left\{ \begin{array}{c} \mathbf{n}_{n+1} \\ \mathbf{h}_{n+1} \\ 0 \end{array} \right\} \]

\[ \left[ \begin{array}{ccc} -1+\gamma \left( \frac{\partial\mathbf{n}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p} + \frac{\partial\mathbf{n}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial\alpha\varepsilon^p} \right) & \gamma \left( \frac{\partial\mathbf{n}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2} + \frac{\partial\mathbf{n}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \alpha} \right) & \mathbf{n} \\ \gamma \left( \frac{\partial\mathbf{h}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p} + \frac{\partial\mathbf{h}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial\alpha\partial\varepsilon^p} \right) & -1+\gamma \left( \frac{\partial\mathbf{h}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2} + \frac{\partial\mathbf{h}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \alpha} \right) & \mathbf{h}\\ \frac{\partial f}{\partial \sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p} + \frac{\partial f}{\partial \mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha\partial\varepsilon^p} & \frac{\partial f}{\partial \mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2} + \frac{\partial f}{\partial \sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial\alpha} & 0 \end{array} \right]_n \left\{ \begin{array}{c} \delta \varepsilon^p \\ \delta \alpha \\ \delta \gamma \end{array} \right\}_{n+1} = \mathbf{R}_n \]

where \(\varepsilon^p_{n+1}=\varepsilon^p_n + \delta\varepsilon^p_{n+1}\), \(\alpha^p_{n+1}=\alpha_n + \delta\alpha_{n+1}\) and \(\Delta\gamma_{n+1} = \Delta\gamma_n+\delta\gamma_{n+1}\). Usually it is assumed free energy split, i.e.

\[ \psi(\varepsilon^e,\alpha) = \psi^{\varepsilon}(\varepsilon^e) + \psi^\alpha(\alpha) \]

in that case linearized residual takes the form

\[ \left[ \begin{array}{ccc} -1+\Delta\gamma \frac{\partial\mathbf{n}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p} & \Delta\gamma \frac{\partial\mathbf{n}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2} & \mathbf{n} \\ \Delta\gamma \frac{\partial\mathbf{h}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p} & -1+\Delta\gamma \frac{\partial\mathbf{h}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2} & \mathbf{h}\\ \frac{\partial f}{\partial \sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p} & \frac{\partial f}{\partial \mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2} & 0 \end{array} \right]_n \left\{ \begin{array}{c} \delta \varepsilon^p \\ \delta \alpha \\ \delta \gamma \end{array} \right\}_{n+1} = -\mathbf{R}_n \]

\[ \chi_{n+1} := \left\{ \begin{array}{c} \delta \varepsilon^p \\ \delta \alpha \\ \delta \gamma \end{array} \right\}_{n+1} = - \mathbf{A}^{-1} \mathbf{R}_n \]

Algorithmic tangent matrix

\[ \delta \sigma := \mathbf{C}^e(\delta \varepsilon - \delta \varepsilon^p) \]

where \(\mathbf{C}^e = \frac{\partial^2 \psi^e}{\partial {\varepsilon^e}^2}\) and \(\delta \varepsilon = \textrm{sym}[ \nabla \delta \mathbf{u}]\)

\[ \left\{ \begin{array}{c} \mathbf{E}^p \\ \mathbf{L}^p \\ \mathbf{Y}^p \end{array} \right\} \delta \varepsilon = - \mathbf{A}^{-1} \left\{ \begin{array}{c} \Delta\gamma \frac{\partial \mathbf{n}}{\partial \sigma}\mathbf{C}^e \\ \Delta\gamma \frac{\partial \mathbf{h}}{\partial \sigma}\mathbf{C}^e \\ \frac{\partial f}{\partial \sigma}\mathbf{C}^e \end{array} \right\} \delta \varepsilon \]

\[ \delta \varepsilon^p = \mathbf{E}^p \delta \varepsilon \]

\[ \delta \alpha = \mathbf{L}^p \delta \varepsilon \]

\[ \delta \gamma = \mathbf{Y}^p \delta \varepsilon \]

\[ \delta \sigma = \mathbf{C}^e(\mathbf{1} - \mathbf{E}^p)\delta \varepsilon \]

\[ \delta \sigma = \mathbf{C}^{ep}\delta \varepsilon \]

where \(\mathbf{C}^{ep} = \mathbf{C}^e-\mathbf{C}^e\mathbf{E}^p\)

For more details and easy access to all source files related to this see group Small Strain Plasticity.

/** \file SmallStrainPlasticity.hpp
 * \ingroup small_strain_plasticity
 * \brief Operators and data structures for small strain plasticity
 * \example SmallStrainPlasticity.hpp
 
 \section problem_formulation Problem formulation
 
 \subsection theory_small_strain Theory
 
 \subsubsection introduction Introduction
 
 In the following document, we define functions,
 \f$\psi(\varepsilon^e,\alpha)\f$, \f$f(\sigma,\mathbf{q})\f$ and
 \f$\Psi(\sigma,\mathbf{q})\f$ which are free Helmholtz energy, yield function
 and flow potential, respectively. Note that free Helmholtz energy is defined
 by stains and kinematic internal state variables, whereas yield function and
 flow potential is defined by fluxes (stress tensor and fluxes work conjugated
 to internal kinematic state variables). In this implantation, we assume that all
 functions are convex. Also, is assumed that energy potential and flow
 potential is positive and takes zero value at zero strain and zero internal
 variable states. It can be shown that for such defined problem, the first law
 of thermodynamics and non-negative dissipation of energy is automatically
 satisfied, see \cite de2011computational \cite simo2006computational.
 Moreover, in the following implementation is assumed that all functions are
 smooth enough that appropriate derivatives could be calculated.
 
 In the following implementation, the backward-Euler difference scheme is
 applied to solve set nonlinear system of equations with unilateral constrains.
 Admissible stress state is enforced by Lagrange multiplier method and
 (Khun-Tucker) loading/unloading conditions, see \cite  simo2006computational.
 That lead to well known \em Closet \em point \em  projection algorithm.
 
 Following \cite simo2006computational can be shown that plastic loading
 unloading conditions can be exclusively characterised in terms of  trial
 state
 \f[
 \varepsilon^{e,\textrm{trial}}_{n+1} := \varepsilon_{n+1} - \varepsilon^p_n
 \f]
 \f[
 \alpha^\textrm{trial}_{n+1} := \alpha_n
 \f]
 \f[
 \sigma^\textrm{trial}_{n+1} := \frac{\partial \psi}{\partial \varepsilon^e}(\varepsilon^{e,\textrm{trial}}_{n+1},\alpha_n)
 \f]
 \f[
 \mathbf{q}^\textrm{trial}_{n+1} := \frac{\partial \psi}{\partial \alpha}(\varepsilon^{e,\textrm{trial}}_{n+1},\alpha_n)
 \f]
 \f[
 f^\textrm{trial}_{n+1} := f(\sigma^\textrm{trial}_{n+1},\mathbf{q}_n)
 \f]
 then loading/unloading conditions are given
 \f[
 \begin{split}
 f^\textrm{trial}_{n+1} < 0,\quad\textrm{elastic step}\,\to\,\Delta\gamma_{n+1} = 0 \\
 f^\textrm{trial}_{n+1} \geq 0,\quad\textrm{plastic step}\,\to\,\Delta\gamma_{n+1} > 0
 \end{split}
 \f]
 If first condition is true \f$f^\textrm{trial}_{n+1} < 0\f$, trial stress is solution
 at integration point. If second condition is true, plastic strains and values
 of internal state variables are calculated using  \em Closet \em point \em
 projection algorithm.
 
 \subsection cloaset_point_projection Closet point projection algorithm
 
 \subsubsection elastic_strain Additive strain decomposition
 
 \f[
 \varepsilon^e_{n+1} := \varepsilon_{n+1} - \varepsilon^p_{n+1}
 \f]
 
 \subsubsection constitutive_equation Constitutive relations
 
 \f[
 \sigma_{n+1} := \frac{\partial \psi}{\partial \varepsilon^e}(\varepsilon^e_{n+1},\alpha_{n+1})
 \f]
 
 \f[
 \mathbf{q}_{n+1} := \frac{\partial \psi}{\partial \alpha}(\varepsilon^e_{n+1},\alpha_{n+1})
 \f]
 
 \subsubsection constrains Constrains
 
 \f[
 f_{n+1}:=f(\sigma_{n+1},\mathbf{q}_{n+1})
 \f]
 
 \subsubsection evolution Evolution functions
 
 Flow vector
 \f[
 \mathbf{n}_{n+1}:=\frac{\partial \Psi}{\partial \sigma}(\sigma_{n+1},\mathbf{q}_{n+1})
 \f]
 
 Hardening law:
 \f[
 \mathbf{h}_{n+1}:=-\frac{\partial \Psi}{\partial \mathbf{q}}(\sigma_{n+1},\mathbf{q}_{n+1})
 \f]
 
 \subsubsection residual Residual
 
 \f[
 \mathbf{R}_{n+1}:=
 \left\{
 \begin{array}{c}
 -\varepsilon^p_{n+1}+\varepsilon^p_n \\
 -\alpha_{n+1}+\alpha_n \\
 f_{n+1}
 \end{array}
 \right\}
 +
 \Delta\gamma_{n+1}
 \left\{
 \begin{array}{c}
 \mathbf{n}_{n+1} \\
 \mathbf{h}_{n+1} \\
 0
 \end{array}
 \right\}
 \f]
 
 \f[
 \left[
 \begin{array}{ccc}
 -1+\gamma
 \left(
 \frac{\partial\mathbf{n}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p}
 +
 \frac{\partial\mathbf{n}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial\alpha\varepsilon^p}
 \right)
 &
 \gamma
 \left(
 \frac{\partial\mathbf{n}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2}
 +
 \frac{\partial\mathbf{n}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \alpha}
 \right)
 &
 \mathbf{n}
 \\
 \gamma
 \left(
 \frac{\partial\mathbf{h}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p}
 +
 \frac{\partial\mathbf{h}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial\alpha\partial\varepsilon^p}
 \right)
 &
 -1+\gamma
 \left(
 \frac{\partial\mathbf{h}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2}
 +
 \frac{\partial\mathbf{h}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \alpha}
 \right)
 &
 \mathbf{h}\\
 \frac{\partial f}{\partial \sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p}
 +
 \frac{\partial f}{\partial \mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha\partial\varepsilon^p}
 &
 \frac{\partial f}{\partial \mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2}
 +
 \frac{\partial f}{\partial \sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial\alpha}
 & 0
 \end{array}
 \right]_n
 \left\{
 \begin{array}{c}
 \delta \varepsilon^p \\
 \delta \alpha \\
 \delta \gamma
 \end{array}
 \right\}_{n+1}
 =
 \mathbf{R}_n
 \f]
 where
 \f$\varepsilon^p_{n+1}=\varepsilon^p_n + \delta\varepsilon^p_{n+1}\f$,
 \f$\alpha^p_{n+1}=\alpha_n + \delta\alpha_{n+1}\f$
 and
 \f$\Delta\gamma_{n+1} = \Delta\gamma_n+\delta\gamma_{n+1}\f$.
 Usually it is assumed free energy split, i.e.
 \f[
 \psi(\varepsilon^e,\alpha) = \psi^{\varepsilon}(\varepsilon^e) + \psi^\alpha(\alpha)
 \f]
 in that case linearized residual takes the form
 \f[
 \left[
 \begin{array}{ccc}
 -1+\Delta\gamma
 \frac{\partial\mathbf{n}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p}
 &
 \Delta\gamma
 \frac{\partial\mathbf{n}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2}
 &
 \mathbf{n}
 \\
 \Delta\gamma
 \frac{\partial\mathbf{h}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p}
 &
 -1+\Delta\gamma
 \frac{\partial\mathbf{h}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2}
 &
 \mathbf{h}\\
 \frac{\partial f}{\partial \sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p}
 &
 \frac{\partial f}{\partial \mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2}
 & 0
 \end{array}
 \right]_n
 \left\{
 \begin{array}{c}
 \delta \varepsilon^p \\
 \delta \alpha \\
 \delta \gamma
 \end{array}
 \right\}_{n+1}
 =
 -\mathbf{R}_n
 \f]
 
 \f[
 \chi_{n+1} :=
 \left\{
 \begin{array}{c}
 \delta \varepsilon^p \\
 \delta \alpha \\
 \delta \gamma
 \end{array}
 \right\}_{n+1}
 =
 -
 \mathbf{A}^{-1}
 \mathbf{R}_n
 \f]
 
 \subsubsection small_strain_algebraic_tangent Algorithmic tangent matrix
 
 \f[
 \delta \sigma := \mathbf{C}^e(\delta \varepsilon - \delta \varepsilon^p)
 \f]
 where \f$\mathbf{C}^e = \frac{\partial^2 \psi^e}{\partial {\varepsilon^e}^2}\f$ and
 \f$\delta \varepsilon = \textrm{sym}[ \nabla \delta \mathbf{u}]\f$
 
 \f[
 \left\{
 \begin{array}{c}
 \mathbf{E}^p \\
 \mathbf{L}^p \\
 \mathbf{Y}^p
 \end{array}
 \right\}
 \delta \varepsilon
 =
 -
 \mathbf{A}^{-1}
 \left\{
 \begin{array}{c}
 \Delta\gamma \frac{\partial \mathbf{n}}{\partial \sigma}\mathbf{C}^e \\
 \Delta\gamma \frac{\partial \mathbf{h}}{\partial \sigma}\mathbf{C}^e \\
 \frac{\partial f}{\partial \sigma}\mathbf{C}^e
 \end{array}
 \right\}
 \delta \varepsilon
 \f]
 
 \f[
 \delta \varepsilon^p = \mathbf{E}^p \delta \varepsilon
 \f]
 
 \f[
 \delta \alpha = \mathbf{L}^p \delta \varepsilon
 \f]
 
 \f[
 \delta \gamma = \mathbf{Y}^p \delta \varepsilon
 \f]
 
 \f[
 \delta \sigma = \mathbf{C}^e(\mathbf{1} - \mathbf{E}^p)\delta \varepsilon
 \f]
 
 \f[
 \delta \sigma = \mathbf{C}^{ep}\delta \varepsilon
 \f]
 where \f$\mathbf{C}^{ep} = \mathbf{C}^e-\mathbf{C}^e\mathbf{E}^p\f$
 
For more details and easy access to all source files related to this see group
\ref small_strain_plasticity.
 
 */
 
/*
 * This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
#ifndef __SMALL_STRAIN_PLASTICITY_HPP
#define __SMALL_STRAIN_PLASTICITY_HPP
 
#ifndef WITH_ADOL_C
  #error "MoFEM need to be compiled with ADOL-C"
#endif
 
PetscErrorCode SmallStrainPlasticityfRes(SNES snes,Vec chi,Vec r,void *ctx);
PetscErrorCode SmallStrainPlasticityfJac(SNES,Vec,Mat,Mat,void *ctx);
 
/** \brief Small strain finite element implementation
* \ingroup small_strain_plasticity
 
*/
struct SmallStrainPlasticity {
 
  
  
 
  MoFEM::Interface &mField;
 
  Tag thPlasticStrain;
  Tag thInternalVariables;
 
  /// \brief  definition of volume element
  struct MyVolumeFE: public MoFEM::VolumeElementForcesAndSourcesCore {
 
    int addToRule;
    MyVolumeFE(MoFEM::Interface &m_field):
    MoFEM::VolumeElementForcesAndSourcesCore(m_field),
    addToRule(0) { }
 
    /** \brief it is used to calculate nb. of Gauss integration points
     *
     * for more details pleas look
     *   Reference:
     *
     * Albert Nijenhuis, Herbert Wilf,
     * Combinatorial Algorithms for Computers and Calculators,
     * Second Edition,
     * Academic Press, 1978,
     * ISBN: 0-12-519260-6,
     * LC: QA164.N54.
     *
     * More details about algorithm
     * http://people.sc.fsu.edu/~jburkardt/cpp_src/gm_rule/gm_rule.html
    **/
    int getRule(int order) {
      // cerr << "order " << order << endl;
      return 2*(order-1)+addToRule;
    }
 
  };
 
  MyVolumeFE feRhs; ///< calculate right hand side for tetrahedral elements
  MyVolumeFE feLhs; //< calculate left hand side for tetrahedral elements
  MyVolumeFE feUpdate; //< update history variables
 
  SmallStrainPlasticity(MoFEM::Interface &m_field):
  mField(m_field),
  feRhs(m_field),
  feLhs(m_field),
  feUpdate(m_field) {
  }
 
  PetscErrorCode createInternalVariablesTag();
 
  /** \brief common data used by volume elements
  * \ingroup nonlinear_elastic_elem
  */
  struct CommonData {
    map<string,vector<VectorDouble> > dataAtGaussPts;
    map<string,vector<MatrixDouble> > gradAtGaussPts;
    vector<VectorDouble> sTress;
    vector<MatrixDouble> materialTangentOperator;
    vector<VectorDouble> plasticStrain;
    vector<VectorDouble> internalVariables;
    vector<VectorDouble> internalFluxes;
    vector<double> deltaGamma;
    double *plasticStrainPtr;
    double *internalVariablesPtr;
    int plasticStrainSize;
    int internalVariablesSize;
    bool bBar;
    CommonData():
    bBar(true) {
    }
  };
  CommonData commonData;
 
  struct OpGetDataAtGaussPts: public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
    vector<VectorDouble > &valuesAtGaussPts;
    vector<MatrixDouble > &gradientAtGaussPts;
    const EntityType zeroAtType;
 
    OpGetDataAtGaussPts(const string field_name,
      vector<VectorDouble > &values_at_gauss_pts,
      vector<MatrixDouble > &gardient_at_gauss_pts
    );
 
    /** \brief operator calculating field data and gradients
    */
    PetscErrorCode doWork(
      int side,EntityType type,DataForcesAndSourcesCore::EntData &data
    );
  };
 
  struct OpGetCommonDataAtGaussPts: public OpGetDataAtGaussPts {
    OpGetCommonDataAtGaussPts(const string field_name,CommonData &common_data);
  };
 
 
  struct OpUpdate: public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
 
    OpUpdate(
      string field_name,
      CommonData &common_data
    );
 
    PetscErrorCode doWork(
      int side,EntityType type,DataForcesAndSourcesCore::EntData &data
    );
 
  };
 
  struct MakeB {
    PetscErrorCode makeB(const MatrixAdaptor &diffN,MatrixDouble &B);
    PetscErrorCode addVolumetricB(const MatrixAdaptor &diffN,MatrixDouble &volB,double alpha);
  };
 
  struct OpAssembleRhs: public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator,MakeB {
 
    CommonData &commonData;
    VectorDouble nF;
    MatrixDouble B,volB;
    
 
    OpAssembleRhs(
      string field_name,
      CommonData &common_data
    );
 
    PetscErrorCode doWork(
      int side,EntityType type,DataForcesAndSourcesCore::EntData &data
    );
 
  };
 
  struct OpAssembleLhs: public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator,MakeB {
 
    CommonData &commonData;
    MatrixDouble K,transK;
    MatrixDouble rowB,colB,CB;
    MatrixDouble rowVolB,colVolB;
    
 
    OpAssembleLhs(
      string field_name,
      CommonData &common_data
    );
 
    PetscErrorCode doWork(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSourcesCore::EntData &row_data,
      DataForcesAndSourcesCore::EntData &col_data
    );
 
  };
 
  /** \brief Closest point projection algorithm
  * \ingroup small_strain_plasticity
 
  */
  struct ClosestPointProjection {
 
    
    VectorDouble sTrain;
    VectorDouble internalVariables;
    VectorDouble internalVariables0;
    VectorDouble plasticStrain0;
    VectorDouble plasticStrain;
    double deltaGamma;
    double tOl;
    int gG;
 
    VectorAdaptor sTress;
    VectorAdaptor internalFluxes;
    ublas::symmetric_matrix<double,ublas::lower> C,D;
    MatrixDouble partialWStrainPlasticStrain;
    VectorAdaptor partialYSigma;
    VectorAdaptor partialYFlux;
    VectorAdaptor partialHSigma;
    VectorAdaptor partialHFlux;
    ublas::symmetric_matrix<double,ublas::lower> partial2HSigma;
    ublas::symmetric_matrix<double,ublas::lower> partial2HFlux;
    MatrixDouble partial2HSigmaFlux;
 
    vector<int> tAgs;
 
    ClosestPointProjection() {}
 
    VectorDouble activeVariablesW;
    virtual PetscErrorCode setActiveVariablesW();
 
    VectorDouble activeVariablesYH;
    virtual PetscErrorCode setActiveVariablesYH();
 
    double w,y,h;
    ublas::vector<adouble> a_sTrain,a_plasticStrain,a_internalVariables;
    ublas::vector<adouble> a_sTress,a_internalFluxes;
    adouble a_w,a_y,a_h;
 
    virtual PetscErrorCode rEcordW();
 
    virtual PetscErrorCode rEcordY();
 
    virtual PetscErrorCode rEcordH();
 
    VectorDouble gradientW;
    MatrixDouble hessianW;
    virtual PetscErrorCode pLayW();
    virtual PetscErrorCode pLayW_NoHessian();
 
    VectorDouble gradientY;
    virtual PetscErrorCode pLayY();
    virtual PetscErrorCode pLayY_NoGradient();
 
    VectorDouble gradientH;
    MatrixDouble hessianH;
    virtual PetscErrorCode pLayH();
 
    Mat A;
    ublas::matrix<double,ublas::column_major> dataA;
    Vec R,Chi,dChi;
    virtual PetscErrorCode createMatAVecR();
 
    virtual PetscErrorCode destroyMatAVecR();
 
    virtual PetscErrorCode evaluatePotentials();
 
    virtual PetscErrorCode cAlculateR(Vec R);
 
    virtual PetscErrorCode cAlculateA();
 
    friend PetscErrorCode SmallStrainPlasticityfRes(SNES,Vec,Vec,void *ctx);
    friend PetscErrorCode SmallStrainPlasticityfJac(SNES,Vec,Mat,Mat,void *ctx);
 
    SNES sNes;
    KSP kSp;
    PC pC;
    SNESLineSearch lineSearch;
    PetscErrorCode snesCreate();
 
    PetscErrorCode snesDestroy();
 
    PetscErrorCode solveColasetProjection();
 
    MatrixDouble Ep,Lp,Cp,Cep;
    VectorDouble Yp;
 
    PetscErrorCode consistentTangent();
 
    virtual PetscErrorCode freeHemholtzEnergy() {
      PetscFunctionBegin;
      SETERRQ(
        PETSC_COMM_SELF,
        MOFEM_NOT_IMPLEMENTED,
        "Not implemented (overload class and implement this function)"
      );
      PetscFunctionReturn(0);
    }
 
    virtual PetscErrorCode yieldFunction() {
      PetscFunctionBegin;
      SETERRQ(
        PETSC_COMM_SELF,
        MOFEM_NOT_IMPLEMENTED,
        "Not implemented (overload class and implement this function)"
      );
      PetscFunctionReturn(0);
    }
 
    virtual PetscErrorCode flowPotential() {
      PetscFunctionBegin;
      SETERRQ(
        PETSC_COMM_SELF,
        MOFEM_NOT_IMPLEMENTED,
        "Not implemented (overload class and implement this function)"
      );
      PetscFunctionReturn(0);
    }
 
  };
 
  struct OpCalculateStress: public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
    MoFEM::Interface &mField;
    CommonData &commonData;
    ClosestPointProjection &cP;
 
    bool initCp;
    bool isLinear,hessianWCalculated;
    
 
 
    OpCalculateStress(
      MoFEM::Interface &m_field,
      string field_name,
      CommonData &common_data,
      ClosestPointProjection &cp,
      bool is_linear = true
    ):
    MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
    mField(m_field),
    commonData(common_data),
    cP(cp),
    initCp(false),
    isLinear(is_linear),
    hessianWCalculated(false) {
    }
 
    Tag thPlasticStrain;
    Tag thInternalVariables;
    PetscErrorCode getTags();
 
    PetscErrorCode setTagsData(
      const EntityHandle tet,const int nb_gauss_pts,const int nb_internal_variables
    );
 
    PetscErrorCode doWork(
      int side,EntityType type,DataForcesAndSourcesCore::EntData &data
    );
 
  };
 
};
 
#endif //__SMALL_STRAIN_PLASTICITY_HPP
 
/***************************************************************************//**
 * \defgroup small_strain_plasticity Small Strain Plasticity
 * \ingroup user_modules
 ******************************************************************************/
 
/***************************************************************************//**
 * \defgroup user_modules User modules
 ******************************************************************************/