v0.14.0
Loading...
Searching...
No Matches
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
*/
/// \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
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;
bool bBar;
bBar(true) {
}
};
vector<VectorDouble > &valuesAtGaussPts;
vector<MatrixDouble > &gradientAtGaussPts;
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);
};
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);
};
string field_name,
CommonData &common_data
);
PetscErrorCode doWork(
int side,EntityType type,DataForcesAndSourcesCore::EntData &data
);
};
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 {
double deltaGamma;
double tOl;
int gG;
ublas::symmetric_matrix<double,ublas::lower> C,D;
ublas::symmetric_matrix<double,ublas::lower> partial2HSigma;
ublas::symmetric_matrix<double,ublas::lower> partial2HFlux;
vector<int> tAgs;
virtual PetscErrorCode setActiveVariablesW();
virtual PetscErrorCode setActiveVariablesYH();
double w,y,h;
ublas::vector<adouble> a_sTress,a_internalFluxes;
virtual PetscErrorCode rEcordW();
virtual PetscErrorCode rEcordY();
virtual PetscErrorCode rEcordH();
virtual PetscErrorCode pLayW();
virtual PetscErrorCode pLayW_NoHessian();
virtual PetscErrorCode pLayY();
virtual PetscErrorCode pLayY_NoGradient();
virtual PetscErrorCode pLayH();
Mat A;
ublas::matrix<double,ublas::column_major> dataA;
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();
PetscErrorCode consistentTangent();
virtual PetscErrorCode freeHemholtzEnergy() {
PetscFunctionBegin;
SETERRQ(
PETSC_COMM_SELF,
"Not implemented (overload class and implement this function)"
);
PetscFunctionReturn(0);
}
virtual PetscErrorCode yieldFunction() {
PetscFunctionBegin;
SETERRQ(
PETSC_COMM_SELF,
"Not implemented (overload class and implement this function)"
);
PetscFunctionReturn(0);
}
virtual PetscErrorCode flowPotential() {
PetscFunctionBegin;
SETERRQ(
PETSC_COMM_SELF,
"Not implemented (overload class and implement this function)"
);
PetscFunctionReturn(0);
}
};
ClosestPointProjection &cP;
bool initCp;
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),
}
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
******************************************************************************/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
PetscErrorCode SmallStrainPlasticityfRes(SNES snes, Vec chi, Vec r, void *ctx)
PetscErrorCode SmallStrainPlasticityfJac(SNES snes, Vec chi, Mat A, Mat, void *ctx)
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
constexpr int order
const FTensor::Tensor2< T, Dim, Dim > Vec
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
constexpr auto field_name
Deprecated interface functions.
@ OPROW
operator doWork function is executed on FE rows
ublas::matrix< double, ublas::column_major > dataA
ublas::symmetric_matrix< double, ublas::lower > C
ublas::symmetric_matrix< double, ublas::lower > partial2HSigma
friend PetscErrorCode SmallStrainPlasticityfJac(SNES, Vec, Mat, Mat, void *ctx)
ublas::symmetric_matrix< double, ublas::lower > partial2HFlux
ublas::symmetric_matrix< double, ublas::lower > D
friend PetscErrorCode SmallStrainPlasticityfRes(SNES, Vec, Vec, void *ctx)
map< string, vector< MatrixDouble > > gradAtGaussPts
map< string, vector< VectorDouble > > dataAtGaussPts
vector< MatrixDouble > materialTangentOperator
PetscErrorCode makeB(const MatrixAdaptor &diffN, MatrixDouble &B)
PetscErrorCode addVolumetricB(const MatrixAdaptor &diffN, MatrixDouble &volB, double alpha)
int getRule(int order)
it is used to calculate nb. of Gauss integration points
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssembleLhs(string field_name, CommonData &common_data)
OpAssembleRhs(string field_name, CommonData &common_data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode setTagsData(const EntityHandle tet, const int nb_gauss_pts, const int nb_internal_variables)
OpGetCommonDataAtGaussPts(const string field_name, CommonData &common_data)
OpGetDataAtGaussPts(const string field_name, vector< VectorDouble > &values_at_gauss_pts, vector< MatrixDouble > &gardient_at_gauss_pts)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
operator calculating field data and gradients
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpUpdate(string field_name, CommonData &common_data)
Small strain finite element implementation.
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
PetscErrorCode createInternalVariablesTag()