v0.14.0
VEC-7: Automatic Differentiation for Plasticity
Note
Prerequisites of this tutorial include FUN-1: Integration on finite element mesh, following SCL-1: Poisson's equation (homogeneous BC), SCL-4: Nonlinear Poisson's equation, VEC-0: Linear elasticity and finally VEC-2: Nonlinear elastic


Note
Intended learning outcome:
  • nonlinear elasticity
  • automatic differentiation
  • linearisation
  • tangent operators
  • discountinuous Galerkin projection
  • users modules

This code, after some refactoring and removing handling cohesive interferes and periodic boundary conditions to serve the purposes of the tutorial is based on code used in [71] (doi: https://doi.org/10.1016/j.compstruct.2016.11.059), see Figure Figure 1.

Figure 1: Three-dimensional nonlinear micro/meso-mechanical response of fiber-reinforced polymer composites. Note paraboloidal surface described bewlo.

However, some limitations of the original code are removed, i.e., the code can be used to solve complex practical problems, see tutorial code specification. Note that MoFEM tutorials are dedicated to contributors and developers keen to extend functionality, rather than to end-users.

You can jump immediately into the details of implementation here link. You also might look at ADOL-C implementation (library used for automatic differentiation), see ADOL-C manual.

Note
This tutorial utilizes automatic differentiation which is a powerful tool, however, here using it we evaluate hessian which has significant efficiency drawbacks. Despite that code has educational, value. Moreover, enables testing some ideas about plasticity models. Some efficiency drawbacks can be mitigated by parallel processing, however, plasticity is often localized to small regions, which requires load balancing. However current implementation works in parallel, load balancing is not implemented, and this tutorial will benefit from such a feature.
Todo:
Fix notation, use bold for vectors and tensors. Commas and dots at the end of equations. Other issues.

Installation

This tutorial is a separate module, located in an independent repository. Users who would like create bespoke solutions for some research or industrial problems would like to create independent repository which is their ownership, licensing and contribution policy. If necessary keep it private. This tutorial is an example of how it can be done.

You can install a module using Spack.pm, which module is an extension to generic user modules. That simplifies use and installation for end users. However for development purposes, you can install the module as follows:

  • Clone repository
    cd $HOME/mofem_install/mofem-cephas/mofem/users_modules
    git clone -b master https://bitbucket.org/mofem/adolc-plasticity.git
  • Build module
    cd $HOME/mofem_install/mofem-cephas/mofem/users_modules/um-build-Rel*
    make rebuild_cache
    make -j 4 install

Tutorial program specification

  • We provide tutorials for elastoplasticity problems in 2D and 3D on the triangle, quadrilateral, tetrahedron and hexahedron elements. In the 2D case, the code can handle plane stress and plane strain. One source code serves all cases.
  • Implementation is generic, dimension agnostic, i.e. code is implemented that dimension of the problem is unknown while implementation, but is known at compile time. Thus we exploit static polymorphism, i.e. C++ templates and partial specialization.
  • The type elements, the order of approximate is not known at compile time, is decided at run time. Thus we exploit dynamic polymorphism, i.e. C++ virtual functions.
  • Code is assumed to be run in parallel, on distributed mesh. No load balancing if performed, thus we assume that the mesh is balanced.
  • Code can handle a "zoo" of different boundary conditions, i.e. Dirichlet, Neumann, and Robin, i.e. spring boundary conditions.
  • Code can handle multiple material blocks, and if necessary can be extended to handle material property fields.
  • Users can easily add new material models.
  • Some points are different approaches between 2D and 3D while post-processing. In the 2D case, we can plot stress and strain on the interior face element, while in the 3D case, we can plot stress and strain volume of the skin, i.e. skin faces. That is a point that we perform branching, i.e. consent expression if a statement is used, at compilation time.
Note
That can/should be code extended to large deformations by following the Mieche procedure [55], which utilizes Lie algebra by applying the logarithmic mapping to the left Cauchy-Green deformation tense. This code is based on the additive decomposition of total strain into elastic and plastic parts. Logarithmic strains, also called Hencky strains, are additive, and that is why the Mieche procedure allows the use of the code below in essence in unchanged form. See VEC-2: Nonlinear elastic or ADV-0: Plastic problem for details.

Theory and problem formulation

Introduction

In the following document, we define functions \(\psi(\varepsilon^e,\alpha)\), \(f(\sigma,\overline{\alpha})\) and \(\Psi(\sigma,\overline{\alpha})\) which are free (Helmholtz) energy, yield function and flow potential, respectively. Note that free energy is defined by strains and kinematic internal state variables, whereas yield function and flow potential are defined by fluxes (stress tensor and fluxes work conjugated to internal kinematic state variables). In this formulation, we assume that all functions are convex. From the convexity of free energy emerges the ellipticity of the partial differential equation, which is solved to satisfy the conservation of the linear momentum, i.e. establish equilibrium. The convexity of the yield function guarantees the uniqueness of the return on the plastic surface, and the convexity of the flow potential guarantees positive dissipation, i.e. the second law of thermodynamics is satisfied. For details see [66]. Moreover, in the following implementation is assumed that all functions are smooth enough that appropriate derivatives can be calculated, in particular, we assume that plastic surfaces do not have corners or apexs.

In the following implementation, the implicit scheme is applied to solve a set nonlinear system of equations with unilateral constraints. The admissible stress state is enforced by the Lagrange multiplier method and (Kuhn-Tucker) loading/unloading conditions, see [66].

Following [66] can be shown that plastic loading unloading conditions can be exclusively characterized 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{\overline{\alpha}}^\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},\overline{\alpha}_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 the second condition is true, plastic strains and values of the internal state variables are calculated using Closet point projection algorithm. For consistency, we also satisfy the condition \(\Delta\gamma_{n+1}f_{n+1}\), at all converged states.

Closet point projection algorithm

Additive strain decomposition

We assume that the total strain is decomposed into elastic and plastic parts, as follows

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

where \(\varepsilon_{n+1} = \textrm{sym}[ \nabla \mathbf{u}_{n+1}]\) is the total strain, \(\varepsilon^e_{n+1}\) is elastic strain, and \(\varepsilon^p_{n+1}\) is the plastic strain. \(\mathbf{u}_{n+1}\) is displacement at time \(t_{n+1}\).

Constitutive relations

The constitutive relations are defined exclusively through the free energy potential as follows, for the stresses and internal fluxes, respectively

\[ \sigma_{n+1} := \frac{\partial \psi}{\partial \varepsilon^e}\left(\varepsilon^e_{n+1},\alpha_{n+1}\right), \quad \mathbf{\overline{\alpha}}_{n+1} := \frac{\partial \psi}{\partial \alpha}\left(\varepsilon^e_{n+1},\alpha_{n+1}\right) \]

Note that the stresses, i.e. fluxes of linear momentum, are work-conjugated to the strains and the internal fluxes \(\overline{\alpha}\) is work conjugated to the kinematic internal state variables \(\alpha\).

Constraints

The admissible inelastic domain is defined by the yield function

\[ f_{n+1}:=f(\sigma_{n+1},\mathbf{\overline{\alpha}}_{n+1}), \quad f_{n+1} \leq 0 \]

The yield stress could be understood as functional, with zero value implicitly defines the yield surface.

Evolution functions

The flow rule defines the evolution of plastic strains, such that

\[ \mathbf{n}_{n+1}:=\frac{\partial \Psi} {\partial \sigma} \left(\sigma_{n+1},\mathbf{\overline{\alpha}}_{n+1}\right), \quad \Delta \varepsilon^p_{n+1} = \Delta \gamma_{n+1} \mathbf{n}_{n+1} \]

where \(\Delta \gamma_{n+1}\) is the plastic multiplier, and has to be determined that for the plastic process, stress remains on the yield surface, i.e. \(f_{n+1} = 0\). The evolution of internal variables is
also defined consistently by the flow potential, such that

\[ \mathbf{h}_{n+1}:=-\frac{\partial \Psi}{\partial \mathbf{\overline{\alpha}}} \left(\sigma_{n+1},\mathbf{\overline{\alpha}}_{n+1}\right), \quad \Delta \alpha_{n+1} = \Delta \gamma_{n+1} \mathbf{h}_{n+1} \]

Residual

Notting that in computational setting loading increment is finite, If we load the body at some point we will exceed the yield stress, thus we reach an inadmissible state. As a consequence, we have to return to the yield surface and stay on it, while material plastically yields. That task leads to an optimization problem with unilateral constraints [66].

Once we establish that the material plastically flows, by calculating the state for trial stress, which is done individually for each integration point, we solve the optimization problem. The standard procedure is the Lagrange method, which introduces an additional variable, i.e. the Lagrange multiplier, which in this context is called the plastic multiplier, i.e. \(\Delta\gamma\), mentioned above. In some cases, we can associate some physical interpretation to the plastic multiplier, however, in general, it is just a mathematical tool.

Finally, the problem at the integration point boils down to the problem of solving a system of nonlinear equations. We have a whole spectrum of possibilities, including the Newthon-Raphson method, trust region methods, and line searchers, available through the PETSc library and its SNES solver.

To use SNES solver, we need to define the residual and the tangent operator. Here residual is defined as follows

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

where \(i\) is the iteration number. The last row of the residual vector has an interpretation of the distance to yield surface and in the process of the solution of the optimization problem, we reduce that distance to the surface, thus naming closes point projection. Residual can be defined in several alternative ways, here we choose the most straightforward naive way to define it, which works well, at least for the problems we tested. From residual, by applying chain derivatives we obtain the tangent matrix, i.e. second quantity required by SNES solver. Tangent matrix is derived as follows

\[ \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]^i \left\{ \begin{array}{c} \delta \varepsilon^p \\ \delta \alpha \\ \delta \gamma \end{array} \right\} = \mathbf{R}^{i} \]

where \(\left(\varepsilon^p_{n+1}\right)^{i+1}=(\varepsilon^p_{n+1})^i + \delta\varepsilon^p_{n+1}\), \(\left(\alpha^p_{n+1}\right)^{i+1}=(\alpha_n)^i + \delta\alpha_{n+1}\) and \(\left(\Delta\gamma_{n+1}\right)^{i+1} = (\Delta\gamma_n)^i+\delta\gamma_{n+1}\). The quantities proceeded by \(\delta\) are sub increments of the Newton-Raphson method, calculated at series of iterations. Usually, it is assumed free energy split, since is easier to construct a convex function, as follows

\[ \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]^i \left\{ \begin{array}{c} \delta \varepsilon^p \\ \delta \alpha \\ \delta \gamma \end{array} \right\} = -\mathbf{R}_n^i \]

Where sub increment of the Newton-Raphson method is calculated as follows

\[ \delta \chi = \left\{ \begin{array}{c} \delta \varepsilon^p \\ \delta \alpha \\ \delta \gamma \end{array} \right\} = - \mathbf{A}^{-1} \mathbf{R}_n^i, \quad \chi^{i+1} = \chi^i + \delta \chi \]

where \(\mathbf{A}\) is the tangent metrix defined above.

Note
If would you like to improve residual definition, let us know, and we will be happy to guide you.

Algorithmic tangent matrix

The classical approach to computational plasticity is based on a nested approach, we solve unilateral optimization-constrained problems at the integration point, and then we assemble the global tangent operator, which we use to solve the global problem, of emerging from the conservation of linear momentum.

The global problem is nonlinear, as a consequence of the nonlinear material model, i.e. elastoplasticity. The global problem is computationally expensive, and solving it efficiently using the Newton-Raphson. The Newton-Raphson method promises quadratic convergence, near the solution, however, it requires the calculation of the tangent operator, which is consistent with the residual. In our case, the residual depends on a nested solution, which is also solved by the Newton-Raphson method. To calculate a consistent tangent operator, at the global system level, we need to calculate the change of plastic strain as a function of of change of total strain, which implicitly depends on the change of internal material state variables.

To derive the consistent tangent operator, start with the constitutive equation

\[ \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}]\) is change of total strain.

Then we introduce the three operators, for completeness, however, we need only first one for the derivation of system level tangent operator. Mentioned operators have form

\[ \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. \]

The operator \(\mathbf{E}^p\) provides us the infinitesimal increment of plastic strans, as a consequence of the infinitesimal increment of total strain strain. Note that through inverse of the tangent operator \(\mathbf{A}\), the increment of plastic strain is coupled with a change of Lagrange plastic multiplier, i.e. containing infinitesimal increment of stress such that it remains plastic surface. We note that consistent tangent operator is implicitly connected to algorithm, and our choice of definition of the residual.

Taking the above, we can see that we can express an increment of the plastic stain internal fluxes, and plastic multiplier as follows

\[ \delta \varepsilon^p = \mathbf{E}^p \delta \varepsilon,\quad \delta \alpha = \mathbf{L}^p \delta \varepsilon,\quad \delta \gamma = \mathbf{Y}^p \delta \varepsilon \]

Plugin first equation from above into the constitutive equation we obtain

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

which directly gives us material consistent tangent operator

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

where \(\mathbf{C}^{ep} = \mathbf{C}^e-\mathbf{C}^e\mathbf{E}^p\). We use \(\mathbf{C}^{ep}\) to assemble the global tangent operator, i.e. the holy grail of computational plasticity.

B-bar method

When material plastic flows, additional constraints are imposed on the finite element deformation, emerging from the flow rule, in particular for J2 plasticity where plastic flow is isochoric, i.e. volume-preserving. That is a case for 3D case and plane strain. In such cases, low-order elements are locking, i.e. elements do not converge fast to the exact solution, fast as Céa's lemma states for elliptic problems, i.e. that the solution error is the best approximation error.

For low-order quadrilateral elements, there are several methods to mitigate the problem. That is reduced integration with hourglass stabilization, selective integration, and the B-bar method. Here we use the B-bar method, whose name comes from the modification finite element differential operator, i.e. \(\mathbf{B}\) operator, in the finite element expression for the stiffness matric \(\int_{V^e} \mathbf{B}^\textrm{T} \mathbf{C} \mathbf{B} \textrm{d}V\). The B-bar added on the top of the \(\overline{\mathbf{B}}\) operator, indicates that the strain axiator of the deformation is averaged on the element, which is equivalent to selective integration. Since the evaluation of constitutive equations are computationally expensive, in particular with automatic differentiation, we use the B-bar method to be more attractive.

The principle of the B-bar method is expressed by the following formula for variation of strain work on stress, i.e.

\[ \delta \mathcal{W} = \int_{V^e} \mathbf{B}^\textrm{T} \sigma \textrm{d}V = \int_{V^e} \delta \varepsilon \sigma \textrm{d}V \approx \int_{V^e} \left[ \delta \varepsilon + \frac{1}{3} \left( (\textrm{tr}[\delta\varepsilon])^\textrm{ave} - \textrm{tr}[\delta\varepsilon] \right)I \right] \sigma \textrm{d}V = \int_{V^e} \overline{\mathbf{B}}^\textrm{T} \sigma \textrm{d}V \]

, where is an average of the trace of the strain tensor, i.e.

\[ (\textrm{tr}[\delta\varepsilon])^\textrm{ave} = \frac{1}{V^e} \int_{V^e} \textrm{tr}[\delta\varepsilon] \textrm{d}V. \]

where \(V^e\) is the volume of the element.

The above formula is implemented using tensorial notation, rather than matrix notation, typical for traditional nodal-based finite element codes. See file ADOLCPlasticity.cpp and snippet below.

/**
* @brief Calculate base for B-bar method
*
* @tparam DIM Dimension of problem
* @param data access to base functions
* @param storage storage for base functions at integration points
* @param b_bar flag to sith on B-bar method
* @param nb_integration_pts number of integration points
* @param w_ptr integration weights
* @return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 6>, 3>
*/
template <int DIM>
MatrixDouble &storage, const bool b_bar,
const int nb_integration_pts, double *w_ptr,
const auto nb_dofs = data.getFieldData().size();
const auto nb_bases = data.getN().size2();
storage.resize(nb_integration_pts, 6 * nb_dofs, false);
// get tensorial base
auto get_ftensor2_symmetric = [&](const int gg, const int rr) {
return getFTensor2SymmetricFromPtr<3>(&storage(gg, 6 * rr));
};
// calculate base
auto calc_base = [&]() {
auto t_t2_diff = get_ftensor2_symmetric(0, 0);
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
auto t_diff = data.getFTensor1DiffN<DIM>(gg, 0);
for (auto b = 0; b != nb_dofs / DIM; ++b) {
t_grad(i, j) = 0;
t_grad(N0, J) = t_diff(J);
t_t2_diff(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
++t_t2_diff;
t_grad(i, j) = 0;
t_grad(N1, J) = t_diff(J);
t_t2_diff(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
++t_t2_diff;
t_grad(i, j) = 0;
if constexpr (DIM == 3) {
t_grad(N2, J) = t_diff(J);
t_t2_diff(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
++t_t2_diff;
}
++t_diff;
}
}
return get_ftensor2_symmetric(0, 0);
};
// calculate volume average of trace
auto calc_vol = [&](auto &&t_t2_dff) {
std::vector<double> vol(nb_dofs, 0);
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
for (auto b = 0; b != nb_dofs; ++b) {
vol[b] += w_ptr[gg] * t_t2_dff(i, i) / DIM;
++t_t2_dff;
}
}
double sum = 0;
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
sum += w_ptr[gg];
}
for (auto &v : vol) {
v /= sum;
}
return vol;
};
// modify base for B-bar method
auto make_b_bar = [&](auto &&vol) {
auto t_t2_diff = get_ftensor2_symmetric(0, 0);
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
for (auto b = 0; b != nb_dofs; ++b) {
const auto trace = t_t2_diff(J, J) / DIM;
t_t2_diff(I, J) += (vol[b] - trace) * t_kd(I, J);
++t_t2_diff;
}
}
return get_ftensor2_symmetric(0, 0);
};
return b_bar ? make_b_bar(calc_vol(calc_base())) : calc_base();
};

As a result of the above code, we obtain the tensorial base functions, which are used to evaluate the stiffness matrix and the right-hand side vector. For the latter case, code is shown in the snippet form ADOLCPlasticity.hpp below

template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode OpRhsImpl<DIM, GAUSS, AssemblyDomainEleOp>::iNtegrate(
double *w_ptr = &(OP::getGaussPts()(DIM, 0));
// Calulate tensorial base functions. Apply bBar method when needed
auto t_diff = MakeB::getFTensor2SymmetricDiffBase(
data, baseStorage, commonDataPtr->bBar, OP::getGaussPts().size2(), w_ptr,
auto t_stress = getFTensor2SymmetricFromMat<3>(commonDataPtr->stressMatrix);
const auto vol = OP::getMeasure();
auto t_w = OP::getFTensor0IntegrationWeight();
for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
double alpha = vol * t_w;
for (int bb = 0; bb != OP::nbRows; ++bb) {
OP::locF[bb] += alpha * (t_stress(i, j) * t_diff(i, j));
++t_diff;
}
++t_w;
++t_stress;
}
}
Note
Another method for resolving the issue with locking would be adding a bubble function on the element, which is a natural way in MoFEM, look to ADV-0: Plastic problem to see how this can be done. You can also look at VEC-0: Linear elasticity to learn how bubble could be efficiently by removing block diagonal entries from the stiffness matrix using the Schur complement and block solver.

Material models

Each material model is a class derived from the base class ADOLCPlasticity::ClosestPointProjection. Each class has three abstract methods, which must be implemented by the user. That is ADOLCPlasticity::ClosestPointProjection::freeHemholtzEnergy, ADOLCPlasticity::ClosestPointProjection::flowPotential and ADOLCPlasticity::ClosestPointProjection::yieldFunction.

If a new model is implemented, it must be registered in the factory, see and follow the pattern below from adolc_plasticity.cpp

/**
* FIXME: Here only array of material models is initalized. Each model has
* unique set of the ADOL-C tags. Pointer is attached based on block name to
* which entity belongs. That will enable heterogeneity of the model, in
* addition of heterogeneity of the material properties.
*/
enum MaterialModel {
VonMisses,
VonMissesPlaneStress,
Paraboloidal,
LastModel
};
const char *list_materials[LastModel] = {"VonMisses", "VonMissesPlaneStress",
"Paraboloidal"};
PetscInt choice_material = VonMisses;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-material", list_materials,
LastModel, &choice_material, PETSC_NULL);
switch (choice_material) {
case VonMisses:
cpPtr = createMaterial<J2Plasticity<3>>();
break;
case VonMissesPlaneStress:
if (SPACE_DIM != 2)
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"VonMissesPlaneStrain is only for 2D case");
cpPtr = createMaterial<J2Plasticity<2>>();
break;
case Paraboloidal:
cpPtr = createMaterial<ParaboloidalPlasticity>();
break;
default:
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Wrong material model");
}
if (approxBase == DEMKOWICZ_JACOBI_BASE && SPACE_DIM == 2)
cpPtr->integrationRule = [](int, int, int p) { return 2 * (p - 2); };
Todo:
For simplicity we predefine one model to the mesh, however, using MoFEM::MeshsetsManager it is possible to assign different models to mesh blocks baes on the block name. We should implement this extension. Block names can be set in Cubit, gMesh, or Salome. In particular MED format can handle block names.

Auxiliary functions

The class has two auxiliary virtual functions, ADOLCPlasticity::ClosestPointProjection::addMatBlockOps, and ADOLCPlasticity::ClosestPointProjection::setParams. Both are optional and can be used to read material properties from mesh blocks or to set material properties. Using this functionality more complex cases can be handled, for example, by reading material properties from a stochastic field.

Reading of material properties from mesh blocks is implemented as an example for von Misses material model, see ADOLCPlasticity::J2Plasticity<3>::addMatBlockOps and ADOLCPlasticity::J2Plasticity<3>::setParams.

J2 plasticity

Von Mises yield function, know also as J2 plasticity, also called HMH plasticity after Huber, Mises and Hencky, who invented pressure-independent strength criterion independently. This model is useful to describe the behavior of a large class of metals.

Free energy density

Free energy density has form

\[ \psi(\varepsilon^e, \alpha, \beta) = \frac{1}{2} \lambda \textrm{tr}[\varepsilon^e]^2 + \mu \varepsilon^e : \varepsilon^e + \sigma_y\alpha + \frac{1}{2} \phi H \alpha^2 + \frac{1}{2} (1-\phi)K \beta^2 \]

where \(\lambda\) and \(\mu\) are Lame'e parameters, \(\sigma_y\) is yield stress, \(\varepsilon^e\) is elastic strain \(\alpha\) is the internal variable associated with isotropic hardening and \(\beta\) is the internal variable associated with kinematics hardening. \(\phi\) is the fraction controlling how the amount of kinematic and isotropic hardening is useful for studying the Bauchinger effect, i.e. plastic hysteresis under cyclic loading.

In general, such a function needs to be convex, which yields elasticity of the problem, and the uniqueness of the solution. Here we restricted implementation to polyconvex form, i.e. that each term separated by addition is convex. That separates in an additive way strains and internal variables, as stated in Introduction, and exploited while assembly the nested tangent operator. Models in small strains, since a symmetric gradient of displacements is used, are invariant to infinitesimal rotations, which provides objectivity to the problem in infinitesimal cases. Also, note that the model can be easily extended to anisotropic cases.

The code for which the above function is implemented is given below from ADOLCPlasticityMaterialModels.hpp

MoFEMErrorCode freeHemholtzEnergy() {
// ADOL-C register parameters (those can varied for integration points)
auto p_lambda = mkparam(lambda); // 0
auto p_mu = mkparam(mu); // 1
auto p_H = mkparam(H); // 2
auto p_K = mkparam(K); // 3
auto p_phi = mkparam(phi); // 4
auto p_sigma_y = mkparam(sigmaY); // 5
// Operator converting Voight and tensor notation
auto t_voight_strain_op = voight_to_strain_op();
// Variable in voight notation
auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
auto t_voight_internal_tensor =
getFTensor1FromPtr<6>(&a_internalVariables[1]);
// Convert variables to tensor notation
tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
tPlasticStrain(i, j) =
t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
tInternalTensor(i, j) =
t_voight_strain_op(i, j, Z) * t_voight_internal_tensor(Z);
// Evaluate elastic strain
tElasticStrain(i, j) = tTotalStrain(i, j) - tPlasticStrain(i, j);
tR = tElasticStrain(i, i);
// Strain energy
a_w = 0.5 * p_lambda * tR * tR +
p_mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
// Isotropic strain energy
a_w +=
a_internalVariables[0] * p_sigma_y +
0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
// Kinematic strain energy
a_w += (0.5 * (1 - p_phi) * p_K) *
(tInternalTensor(i, j) * tInternalTensor(i, j));
}

Yield function and flow potential

The yield function has form

\[ f(\sigma, \overline{\alpha}) = \sqrt{\frac{3}{2} \eta:\eta} - \overline{\alpha} \]

where

\[ \eta=\textrm{dev}[\sigma]-\overline{\beta} \]

where \(\sigma\) is the stress tensor and \(\overline{\alpha}\) is the scalar internal flux variable, and \(\overline{\beta}\) is tensorial internal flux.

Note that dissipation density can be expressed as

\[ \mathcal{D} = \alpha \overline{\alpha} + \beta : \overline{\beta} \geq 0 \]

The code for which the above function is implemented is given below from ADOLCPlasticityMaterialModels.hpp

MoFEMErrorCode evalF() {
auto t_voight_stress_op = voight_to_stress_op();
auto t_vioght_stress = getFTensor1FromPtr<6>(&a_sTress[0]);
auto t_back_stress = getFTensor1FromPtr<6>(&a_internalFluxes[1]);
tStress(i, j) = t_voight_stress_op(i, j, Z) * t_vioght_stress(Z);
tBackStress(i, j) = t_voight_stress_op(i, j, Z) * t_back_stress(Z);
constexpr double third = boost::math::constants::third<double>();
tR = tStress(i, i) / 3.;
tDeviator(i, j) = tStress(i, j) - tR * t_kd(i, j) - tBackStress(i, j);
// Fix the constant to mach uniaxial test
constexpr double c = 3. / 2.;
f = c * tDeviator(i, j) * tDeviator(i, j);
}
MoFEMErrorCode yieldFunction() {
CHKERR evalF();
a_y = sqrt(f) - a_internalFluxes[0];
}

Since for J2 plasticity, we assume the associated flow rule, the flow potential is the same as the yield function, i.e. \(\Psi(\sigma, \overline{\alpha}) = f(\sigma, \overline{\alpha})\). In such cases, we can show that the problem could be expressed as a constrained minimization problem of the Lagrangian function

\[ \min_{\alpha, \beta} \max_{\dot{\lambda}} \mathcal{L}(\alpha, \beta,\overline{\alpha}, \overline{\beta}, \dot{\lambda}) := \alpha \overline{\alpha} + \beta : \overline{\beta} + \dot{\gamma} f(\sigma(\varepsilon^e), \overline{\alpha}) \]

where implicit relations between elastic strain and stress, and internal kinematic variables and fluxes are defined by free energy potential. For such cases we dissipation of energy is maximal. If for such a case model prediction matches within reason experiments for curtain range of loading, it indicates that the model captures all fundamental physics of the problem and constitutive equations can not be expressed with a smaller number of parameters that need to be estimated experimentally by independent tests.

Note
Any nonassociative flow rule requires additional parameters, and thus the model is less fundamental. For example, if we model granular material such as soil, their dissipation also takes place by friction grains sliding against each other or rotating. The letter part, associated with rotation is not captured by the kinematics of classical continuum description, and many experimental results with nonassociative flow rule are added, with the cost of additional parameters.

J2 plasticity in plane stress

Here we use the classical implementation of finite elements, within displacements as a primary variable, we use the so-called strain-driven approach. For the plane strain case such an approach is natural, and simply emerges from type of use of finite elements, i.e. 2D elements or 3D, and the model presented above remains unchanged.

This is not the case for plane stress, where the constraint that the stress, let's say, in z-direction is zero, needs to be explicitly enforced. That typically leads to some implementation difficulties, i.e. plastic surface is no longer circular, it has an elliptical shape, and a radial return mapping algorithm no longer works. That is a direct result of constraints, which yield a change in the definition of deviatoric stress defined in the 2D word of the plane stress problem.

However, this implementation is straightforward. Since the model is based on potential functions, we choose to implement integration points always in 3D. Specialization between 2D and 3D cases is at the final step while the assembly and construction of finite element base/shape functions, in particular, when the b-bar method is used, to solve the problem of volumetric locking. Note that the plane stress problem is it is not an actual 2D problem, but the reduction of the 3D problem with additional static constrain.

Thus, the implementation of plane stress requires only a change in the definition of free energy, taking into account static constraints, see a snippet from ADOLCPlasticityMaterialModels.hpp

template <> struct J2Plasticity<2> : public J2Plasticity<3> {
using J2Plasticity<3>::J2Plasticity;
MoFEMErrorCode freeHemholtzEnergy() {
auto p_lambda = mkparam(lambda); // 0
auto p_mu = mkparam(mu); // 1
auto p_H = mkparam(H); // 2
auto p_K = mkparam(K); // 3
auto p_phi = mkparam(phi); // 4
auto p_sigma_y = mkparam(sigmaY); // 5
auto t_voight_strain_op = voight_to_strain_op();
auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
auto t_voight_internal_tensor =
getFTensor1FromPtr<6>(&a_internalVariables[1]);
tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
tPlasticStrain(i, j) =
t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
tInternalTensor(i, j) =
t_voight_strain_op(i, j, Z) * t_voight_internal_tensor(Z);
tElasticStrain(i, j) = tTotalStrain(i, j) - tPlasticStrain(i, j);
//! [Plane stress]
// Calculate st rain at ezz enforcing plane strain case
tElasticStrain(2, 2) =
(-tElasticStrain(0, 0) * p_lambda - tElasticStrain(1, 1) * p_lambda) /
(p_lambda + 2 * p_mu);
//! [Plane stress]
tR = tElasticStrain(i, i);
a_w = 0.5 * p_lambda * tR * tR +
p_mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
// plastic part
// isotropic
a_w +=
a_internalVariables[0] * p_sigma_y +
0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
// kinematic
a_w += (0.5 * (1 - p_phi) * p_K) *
(tInternalTensor(i, j) * tInternalTensor(i, j));
}
};

in particular plane stress constrain is explicitly enforced in the lines

// Calculate st rain at ezz enforcing plane strain case
tElasticStrain(2, 2) =
(-tElasticStrain(0, 0) * p_lambda - tElasticStrain(1, 1) * p_lambda) /
(p_lambda + 2 * p_mu);

what is

\[ \varepsilon_{zz} = \frac{-\lambda}{\lambda + 2\mu} \left(\varepsilon_{xx} + \varepsilon_{yy}\right) \]

Paraboloidal model

That is pressure dependent model, which non-associated flow rule. For some applications see papers for reference [71] and [68].

Free energy density

Free energy density has form

\[ \psi = \frac{1}{2} \lambda \textrm{tr}[\varepsilon]^2 + \mu \varepsilon : \varepsilon + \sigma_{yt}\alpha_0 + \frac{1}{2} H_t \alpha_0^2 + \sigma_{yc}\alpha_1 + \frac{1}{2} H_c \alpha_1^2 \]

where \(\lambda\) and \(\mu\) are Lam\'e parameters, \(H_t\) is the hardening modulus in tension, \(H_c\) is the hardening modulus in compression.

The code for which the above function is implemented is given in snippet form ADOLCPlasticityMaterialModels.hpp

MoFEMErrorCode freeHemholtzEnergy() {
auto t_voight_strain_op = voight_to_strain_op();
auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
tPlasticStrain(i, j) =
t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
tElasticStrain(i, j) = tTotalStrain(i, j) - tPlasticStrain(i, j);
tR = tElasticStrain(i, i);
a_w = 0.5 * lambda * tR * tR +
mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
a_w += sigmaYt * a_internalVariables[0] +
0.5 * Ht * a_internalVariables[0] * a_internalVariables[0];
a_w += sigmaYc * a_internalVariables[1] +
0.5 * Hc * a_internalVariables[1] * a_internalVariables[1];
}

Yield function

The yield function has form

\[ y = 6J_2 + 2I_1\left(\overline{\alpha_1}-\overline{\alpha_0}\right) - 2\overline{\alpha_0} \,\overline{\alpha_1} \]

where the internal fluxes are defined by free energy functional and have from, for compression and tension, respectively, as follows

\[ \overline{\alpha_0}=\frac{\partial \psi}{\partial \alpha_0}=\sigma_{yt} + H_t \alpha_0, \quad \overline{\alpha_1}=\frac{\partial \psi}{\partial \alpha_1}=\sigma_{yc} + H_c \alpha_1 \]

where \(\sigma_{yt}\) is the yield stress in tension, \(\sigma_{yc}\) is the yield stress in compression.

The J2 and I1 stress invariants are defined as follows

\[ I_1 = \textrm{tr} (\boldsymbol{\sigma}), \quad J_2 = \frac{1}{2} \eta:\eta, \]

where \(\eta\) is the deviatoric stress tensor, defined as follows

\[ \eta=\textrm{dev}[\boldsymbol{\sigma}] = \sigma - \frac{1}{3} \textrm{tr} (\boldsymbol{\sigma}) I \]

The implementation of the yield function is given snippet from ADOLCPlasticityMaterialModels.hpp

MoFEMErrorCode evalF() {
auto t_voight_stress_op = voight_to_stress_op();
auto t_vioght_stress = getFTensor1FromPtr<6>(&a_sTress[0]);
tStress(i, j) = t_voight_stress_op(i, j, Z) * t_vioght_stress(Z);
tR = tStress(i, i);
I1 = tR;
tR /= 3.;
tDeviator(i, j) = tStress(i, j) - tR * t_kd(i, j);
J2 = 0.5 * (tDeviator(i, j) * tDeviator(i, j));
}
MoFEMErrorCode yieldFunction() {
CHKERR evalF();
a_y = 6 * J2 + 2 * I1 * (a_internalFluxes[1] - a_internalFluxes[0]) -
2 * a_internalFluxes[0] * a_internalFluxes[1];
}

Flow potential

The flow potential has form

\[ \Psi = 6J_2 + 2\alpha I_1 \left(\overline{\alpha_1}-\overline{\alpha_0}\right) - 2\overline{\alpha_0} \,\overline{\alpha_1} \]

\[ \alpha= \frac{1-2\nu_p}{1+\nu_p} \]

where \(\nu_p\) is the Poisson ratio of the plastic strain, which for \(\nu_p=0\) model reduces to the associated flow rule, for \(\nu_p=0.5\) the model reduces to the isochoric flow rule. The \(\nu_p\) is an extra material parameter, which has to be tuned to match experimental results. Note that for associated flow rule yield surface and flow potential are parameter-free.

The implementation of the flow potential is given snippet from ADOLCPlasticityMaterialModels.hpp

MoFEMErrorCode flowPotential() {
CHKERR evalF();
double alpha =
(1 - 2 * nup) /
(1 + nup); // relation between alpha and plastic poisson's ratio
a_h = 6 * J2 +
2 * alpha * I1 * (a_internalFluxes[1] - a_internalFluxes[0]) -
2 * a_internalFluxes[0] * a_internalFluxes[1];
}

Assembly and integration

The integration of ADOL-C plasticity models with core MoFEM functionality involves three elements, pushing operators for assembly of the right-hand side, i.e. internal forces, and pushing operators to the left-hand side, i.e. tangent matrix. Users can do that using so-called operators factory, which encapsulates all necessary functionality.

The principle behind the operators and operators factory is to create resilient code, that each element of the finite element system is isolated, which interferes with other elements by strictly defined rules. That enforces the modularity of the code. That involves that data are local, and each operator is a kind of black box, which does not depend on the state of other operators. That simplicity of code development, testing, and allows the composition of many components, i.e. different phenomena together.

In the case of plasticity, beyond the assembly linear system of equations, by pushing operators to the left-hand side, and the right-hand side, we need to update history model internal variables, and plastic strains, since elastic-plastic material remembers load history,

Also, the postprocessing of the elastoplastic model requires special attention, while post-processing, since stresses can not be evaluated at post-processing mesh nodes, since the history of deformation is unknown at arbitrary material points, only at integration points used while assembly of the linear system of equations.

As you will see each element, i.e. assembly, updating history, and postprocessing is a black box, the only common data structure is the material model, in which all elements share. This tutorial with little effort can be integrated with ADV-1: Contact problem. Moreover, by evaluating dissipation and assuming that part of it is converted to heat, the user can integrate code with ADV-2: Thermo-elastic example. Another option is coupling with rates of change of the volume, rendering complex physics for hydromechanical problems, i.e. ADV-5: Seepage in elsatic body example.

Note that the tutorial has a general structure that resembles the tasks presented above, see the following code snippet for adolc_plasticity.cpp

Assembly

The pushing of operators to the left-hand side and right-hand side is done by a snippet from adolc_plasticity.cpp

auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
// assemble operator to the right hand side
auto add_domain_ops_lhs = [&](auto &pip) {
// push forward finite element bases from reference to physical element
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
// create local common data
auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
// calculate displacement gradients at integration points
"U", common_data_ptr->getGardAtGaussPtsPtr()));
// assemble tangent operator
CHKERR opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr);
};
auto add_domain_ops_rhs = [&](auto &pip) {
// push forward finite element bases from reference to physical element
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
// create local common data
auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
// calculate displacement gradients at integration points
"U", common_data_ptr->getGardAtGaussPtsPtr()));
// assemble residual
CHKERR opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr, Sev::inform);
};
// Push operators to the left hand side pipeline. Indicate that domain (i.e.
// volume/interior) element is used.
CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
// Push operators to the right hand side pipeline.
CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
}

Note "ADOLCMAT" is an arbitrary name of the mesh block, with carrying information about the material model parameters.

Implementation is dimension agnostic, i.e. dimensions of the problem are specialized while code compilation and shape of the elements shape and approximation base orders are decided during the runtimes. However, the implementation of the physics, i.e. elastoplastic phenomena, is free from those considerations about dimensions, order base, or element type.

Have particular attention to ADOLCPlasticity::opFactoryDomainRhs and ADOLCPlasticity::opFactoryDomainLhs, which are directly responsible for pushing operators to the right-hand side and left-hand side, respectively. Those follow patterns that you will find in the other tutorials, i.e. templates specialized for the dimension of the problem, assembly type, integrate method, and operator type. The dimensions of the people are obvious, i.e. 1D, 2D, or 3D.

Assembly could be for example PETSC assembly, that uses simply PETSCc functionality and uses PETSC function MatSetValues. In principle, technically, MatSetValues is specialized in MoFEM, see for example MoFEM::MatSetValues<AssemblyTypeSelector<SCHUR>>. However, users can imagine more complex cases, for example, special assembly for block diagonal parts, i.e. bubbles functions, and prefer calculating Schur complement, see for example VEC-0: Linear elasticity. Other types of assembly are to be written for matrix-free approach, or to GPU. This typically can be separate from the implementation of the physics, i.e. elastoplasticity, that why templatised and specialized if necessary.

Similar considerations are for the integration method, typically GAUSS integrate is used, however, users seeking efficiency can use other methods, i.e. Duffy integrate for fast integration for high-order elements.

Finally, the operator type is responsible for handling abstractions related to specific shapes, dimensions, and element types. We exploit here powerful C++ type checking, to validate that dimension, assembly type, integration method, and operator type are compatible.

Updating history

Setting up the update of history variables is related to the solver you use, i.e. if you would use SNES solver, and program load stepping by yourself, at the end of each step, you need to call a finite element with the right pipeline operators by yourself.

However, in most tutorials, we use PETSC time solver (TS) to carry out time-stepping and inform TS to perform an update of history variables. By doing that we utilize the functionality of PETSC, tested by many users around the world. TS has generic functionality to perform such tasks, i.e. to update history variables, see TSSetPostStep. The only issue is that TSSetPostStep requires passing a pointer to pure "C" function so we have to take care by passing data about the material model
through the same static global variable.

We start with a snippet from adolc_plasticity.cpp, creating time solver

auto ts = pipeline_mng->createTSIM();
// Set time solver
double ftime = 1;
CHKERR TSSetMaxTime(ts, ftime);
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
auto D = createDMVector(simple->getDM());
CHKERR TSSetSolution(ts, D);
// Set monitor, step adaptivity, and post-step to update history variables
CHKERR set_up_monitor(ts);
CHKERR set_up_post_step(ts);
CHKERR set_up_adapt(ts);
CHKERR TSSetFromOptions(ts);
CHKERR TSSolve(ts, NULL);

Next, we call function setting TS post-step function,

auto set_up_post_step = [&](auto ts) {
// create finite element (pipeline) and populate it with operators to
// update history variables
auto create_update_ptr = [&]() {
auto fe_ptr = boost::make_shared<DomainEle>(mField);
fe_ptr->getRuleHook = cpPtr->integrationRule;
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(fe_ptr->getOpPtrVector(),
{H1});
auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
fe_ptr->getOpPtrVector().push_back(
"U", common_data_ptr->getGardAtGaussPtsPtr()));
opFactoryDomainUpdate<SPACE_DIM>(mField, fe_ptr->getOpPtrVector(),
"ADOLCMAT", common_data_ptr, cpPtr),
"push update ops");
return fe_ptr;
};
// ts_update_ptr is global static variable
createTSUpdate(simple->getDomainFEName(), create_update_ptr());
//! [TS post-step function]
// This is pure "C" function which we can to the TS solver
auto ts_step_post_proc = [](TS ts) {
CHKERR ts_update_ptr->postProcess(ts);
};
//! [TS post-step function]
// finally set up post-step
CHKERR TSSetPostStep(ts, ts_step_post_proc);
};

Note that this involves creating finite elements, creating data structures used internally by TS post-step function, and finally registering the function with TS.

Finally, we have to implement the intervention on the elements and execute operators associated with the pipeline. See snippet in ADOLCPlasticity.cpp

/**
* @brief TS update history variables
*/
struct TSUpdateImpl : public TSUpdate {
TSUpdateImpl(std::string fe_name, boost::shared_ptr<FEMethod> fe_ptr);
MoFEMErrorCode postProcess(TS ts);
private:
boost::shared_ptr<FEMethod> fePtr;
const std::string feName;
};
TSUpdateImpl::TSUpdateImpl(std::string fe_name,
boost::shared_ptr<FEMethod> fe_ptr)
: feName(fe_name), fePtr(fe_ptr) {}
MoFEMErrorCode TSUpdateImpl::postProcess(TS ts) {
DM dm;
CHKERR TSGetDM(ts, &dm);
// This is where update is preformed. Element list is iterated and on each
// element pipeline accessed throng fePtr is executed.
CHKERR DMoFEMLoopFiniteElements(dm, feName, fePtr);
}
/**
* @brief Create test update object
*
* @param fe_name
* @param fe_ptr
* @return boost::shared_ptr<TSUpdate>
*/
boost::shared_ptr<TSUpdate> createTSUpdate(std::string fe_name,
boost::shared_ptr<FEMethod> fe_ptr) {
return boost::make_shared<TSUpdateImpl>(fe_name, fe_ptr);
}

You can see again here several levels of code encapsulation, which is necessary for integrating the code with other modules. Note that the user in the future can have other models that also require some update of the state, then it will simply call several update functions inside this snippet in adolc_plasticity.cpp

// This is pure "C" function which we can to the TS solver
auto ts_step_post_proc = [](TS ts) {
CHKERR ts_update_ptr->postProcess(ts);
};

Postprocessing

Postprocessing is a bit more complicated than in the case of linear elasticity. We use Lagrangian formulation, also called material formulation, in which each integration point has an associated material point that does not change through the analysis. For each material point, i.e. integration point, we record history, since the current state, stress, and elastic strain, depend on it.

However, to post-process results, we need to know the current state, such as stress at the mesh nodes, to visualize results in post-processing software, for example Paraview.

To solve that issue we use Discontinuous Galerkin (DG) projection, see implementation in MoFEM::OpLoopThis or see code in dg_projection.cpp. The idea is to build local L2 based on the element, then project the stress and plastic strain on the local base, and then extrapolate to the mesh nodes. DG projection provides general functionality for scalar, vectorial and tensorial quantities, in this tutorial, we use it for stress and plastic strain, which are second-order symmetric tensors. See a snippet from adolc_plasticity.cpp which implements the above idea

// Post-process domain, i.e. volume elements. If 3d case creates stresses
// are evaluated at points which are trace of the volume element on boundary
auto post_proc_ele_domain = [this](auto &pip_domain, auto &fe_name) {
// Push bases on reference element to the phyiscal element
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip_domain, {H1});
// calculate displacement gradients at nodes of post processing mesh. For
// gradient DG projection is obsolete, since displacements can be
// evaluated at arbitrary points.
auto grad_ptr = boost::make_shared<MatrixDouble>();
pip_domain.push_back(
grad_ptr));
// Create operator to run (neted) pipeline in operator. InteriorElem
// element will have iteration rule and bases as domain element used to
// solve problem, and remember history variables.
using InteriorElem = DomainEle;
auto op_this = new OpLoopThis<InteriorElem>(mField, fe_name, Sev::noisy);
// add interior element to post-processing pipeline
pip_domain.push_back(op_this);
// get pointer to interior element, which is in essence pipeline which
// pipeline.
auto fe_physics = op_this->getThisFEPtr();
auto evaluate_stress_on_physical_element = [&]() {
// evaluate stress and plastic strain at Gauss points
fe_physics->getRuleHook = cpPtr->integrationRule;
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
fe_physics->getOpPtrVector(), {H1});
auto common_data_ptr =
boost::make_shared<ADOLCPlasticity::CommonData>();
// note that gradient of displacements is evaluate again, at
// physical nodes
fe_physics->getOpPtrVector().push_back(
"U", common_data_ptr->getGardAtGaussPtsPtr()));
CHKERR cpPtr->addMatBlockOps(mField, fe_physics->getOpPtrVector(),
"ADOLCMAT", Sev::noisy);
fe_physics->getOpPtrVector().push_back(
getRawPtrOpCalculateStress<SPACE_DIM>(mField, common_data_ptr,
cpPtr, false));
return common_data_ptr;
};
auto dg_projection_froward_and_back = [&](auto &common_data_ptr) {
// here we do actual projection of stress and plastic strain to DG space
int dg_order = approxOrder - 1;
auto entity_data_l2 =
boost::make_shared<EntitiesFieldData>(MBENTITYSET);
auto mass_ptr =
boost::make_shared<MatrixDouble>(); //< projection operator (mass
// matrix)
fe_physics->getOpPtrVector().push_back(new OpDGProjectionMassMatrix(
dg_order, mass_ptr, entity_data_l2, approxBase, L2));
auto coeffs_ptr_stress = boost::make_shared<MatrixDouble>();
// project stress on DG space on physical element
fe_physics->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
common_data_ptr->getStressMatrixPtr(), coeffs_ptr_stress, mass_ptr,
entity_data_l2, approxBase, L2));
// project strains plastic strains DG space on physical element
auto coeffs_ptr_plastic_strain = boost::make_shared<MatrixDouble>();
fe_physics->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
common_data_ptr->getPlasticStrainMatrixPtr(),
coeffs_ptr_plastic_strain, mass_ptr, entity_data_l2, approxBase,
L2));
// project stress and plastic strain for DG space on post-process
// element
pip_domain.push_back(new OpDGProjectionEvaluation(
common_data_ptr->getStressMatrixPtr(), coeffs_ptr_stress,
entity_data_l2, approxBase, L2));
pip_domain.push_back(new OpDGProjectionEvaluation(
common_data_ptr->getPlasticStrainMatrixPtr(),
coeffs_ptr_plastic_strain, entity_data_l2, approxBase, L2));
};
auto common_data_ptr = evaluate_stress_on_physical_element();
dg_projection_froward_and_back(common_data_ptr);
return boost::make_tuple(grad_ptr, common_data_ptr->getStressMatrixPtr(),
common_data_ptr->getPlasticStrainMatrixPtr());
};
Note
For 3D code implementation have an additional abstraction layer, in which the post-processing element is boundary elements, i.e. face. Then post-processing points at which mesh is evaluated are projected from the face element (boundary) to the volume element (domain), and then the DG projection is run on the domain element. That is part of more general MoFEM functionality, which is used in ADV-1: Contact problem, and in evaluating integrals on mesh skeleton, for methods such as Discontinuous Galerkin, see SCL-11: Discountinous Galrekin for Poisson problem. The key operator in this case is MoFEM::OpLoopSide.

Example

The best way to run and play with the code is to use the JupyterHub and MoFEM docker container. Follow the instructions in Installation with Docker - JupyterHub Then install the release version for code development, you will find a notebook in the JupterHub folder to do that. If you install the module from scratch, i.e. clone repository, build and install the module, you have to recreate Spack symlink, you simply rerun the last section of the install_from_source_release.md. The outcome of this will be that in folder um_view_release you will find the directory adolc_plasticity, and in it adolc_plasticity.md notebook. Run it and play with the parameters.

You can upload a mesh file and set up the problem in Cubit or gMesh. Also notebook creates zipped VTK files which you can download and post-process in Paraview. Any questions please write to Q&A.

The following figures show the results of the kind of simulation which you can do in Jupyter notebook, see Figure 2 and Figure 3.

Figure 2: Load history and load vs displacement response for a plate with the hole. Red color is for J2 plasticity response, blue paraboloidal model

Figure 3: On top is the solution for paraboloidal model and on the bottom solution for J2 plasticity.

How to cite?

If you find code, or this tutorial useful, and use it as a guide for your work, please cite the following paper The Journal of Open Source Software (follow the link on the badge).

Code

/**
* \file adolc_plasticity.cpp
* \ingroup adoc_plasticity
* \example adolc_plasticity.cpp
*
* \brife Implementation of plasticity problem using ADOLC-C
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
/** Select finite element type for integration on domain based on problem
* dimension
*/
/** Select finite element type for integrate on boundary based on problem
* dimension
*/
using BoundaryEle =
/** Operators used to assemble domain integrals
*/
/** Operators used to assemble boundary integrals
*/
/** Linear forms used to integrate body forces
*/
/** Select linear froms reading data from blockest (e.g. "BODY_FORCE") and
* applying body force.
*/
using OpBodyForce =
// Select natural BC boundary condition
// Use natural boundary conditions implemented with adv-1 example which includes
// spring BC
/** Use Gauss quadrature rule and PETSc assembly to integrate neural boundary
* conditions. Select linear forms operators.
*/
/* Use specialization from adv-1 integrating boundary conditions on forces and
* with springs. OP is used to integrate the right hand side
*/
/** Use Gauss quadrature rule and PETSc assembly to integrate neural boundary
* conditions. Select bi-linear forms operators.
*/
/** Use specialization from adv-1 integrating boundary conditions on forces and
* with springs
*/
// Note: We create only skin post-processing mesh.
template <int DIM>
struct PostProcEleByDim; //< Select finite elements and operators used for
// postprocessing
template <> struct PostProcEleByDim<2> {
};
template <> struct PostProcEleByDim<3> {
};
using namespace ADOLCPlasticity;
PlasticProblem(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode gettingNorms();
MoFEMErrorCode outputResults();
MoFEMErrorCode checkResults();
int approxOrder = 2;
boost::shared_ptr<ClosestPointProjection> cpPtr; ///< Closest point
///< projection
};
//! [Run problem]
CHKERR readMesh();
CHKERR setupProblem();
CHKERR boundaryCondition();
CHKERR assembleSystem();
CHKERR solveSystem();
CHKERR gettingNorms();
CHKERR outputResults();
CHKERR checkResults();
}
//! [Run problem]
//! [Read mesh]
auto simple = mField.getInterface<Simple>();
CHKERR simple->getOptions();
CHKERR simple->loadFile();
if (simple->getDim() != SPACE_DIM)
SETERRQ2(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Wrong dimension of mesh %d != %d", simple->getDim(), SPACE_DIM);
}
//! [Read mesh]
//! [Set up problem]
Simple *simple = mField.getInterface<Simple>();
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
PetscInt choice_base_value = AINSWORTH;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
LASBASETOPT, &choice_base_value, PETSC_NULL);
switch (choice_base_value) {
case AINSWORTH:
MOFEM_LOG("WORLD", Sev::inform)
<< "Set AINSWORTH_LEGENDRE_BASE for displacements";
break;
case DEMKOWICZ:
approxBase = DEMKOWICZ_JACOBI_BASE;
MOFEM_LOG("WORLD", Sev::inform)
<< "Set DEMKOWICZ_JACOBI_BASE for displacements";
break;
default:
approxBase = LASTBASE;
break;
}
// Add field
CHKERR simple->addDomainField("U", H1, approxBase, SPACE_DIM);
CHKERR simple->addBoundaryField("U", H1, approxBase, SPACE_DIM);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &approxOrder, PETSC_NULL);
CHKERR simple->setFieldOrder("U", approxOrder);
CHKERR simple->setUp();
//! [Material models selection]
/**
* FIXME: Here only array of material models is initalized. Each model has
* unique set of the ADOL-C tags. Pointer is attached based on block name to
* which entity belongs. That will enable heterogeneity of the model, in
* addition of heterogeneity of the material properties.
*/
enum MaterialModel {
VonMisses,
VonMissesPlaneStress,
Paraboloidal,
LastModel
};
const char *list_materials[LastModel] = {"VonMisses", "VonMissesPlaneStress",
"Paraboloidal"};
PetscInt choice_material = VonMisses;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-material", list_materials,
LastModel, &choice_material, PETSC_NULL);
switch (choice_material) {
case VonMisses:
cpPtr = createMaterial<J2Plasticity<3>>();
break;
case VonMissesPlaneStress:
if (SPACE_DIM != 2)
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"VonMissesPlaneStrain is only for 2D case");
cpPtr = createMaterial<J2Plasticity<2>>();
break;
case Paraboloidal:
cpPtr = createMaterial<ParaboloidalPlasticity>();
break;
default:
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Wrong material model");
}
if (approxBase == DEMKOWICZ_JACOBI_BASE && SPACE_DIM == 2)
cpPtr->integrationRule = [](int, int, int p) { return 2 * (p - 2); };
//! [Material models selection]
}
//! [Set up problem]
//! [Boundary condition]
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto simple = mField.getInterface<Simple>();
auto bc_mng = mField.getInterface<BcManager>();
auto time_scale = boost::make_shared<TimeScale>();
auto rule = [](int, int, int p) { return 2 * p; };
CHKERR pipeline_mng->setDomainRhsIntegrationRule(cpPtr->integrationRule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(cpPtr->integrationRule);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule);
// Add Natural BCs to RHS
pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
Sev::inform);
// Add Natural BCs to LHS
pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", Sev::inform);
//! Add body froces
pipeline_mng->getOpDomainRhsPipeline(), mField, "U", {time_scale},
"BODY_FORCE", Sev::inform);
// Essential BC
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
"U", 2, 2);
CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
}
//! [Boundary condition]
//! [Push operators to pipeline]
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
// assemble operator to the right hand side
auto add_domain_ops_lhs = [&](auto &pip) {
// push forward finite element bases from reference to physical element
// create local common data
auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
// calculate displacement gradients at integration points
"U", common_data_ptr->getGardAtGaussPtsPtr()));
// assemble tangent operator
CHKERR opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr);
};
auto add_domain_ops_rhs = [&](auto &pip) {
// push forward finite element bases from reference to physical element
// create local common data
auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
// calculate displacement gradients at integration points
"U", common_data_ptr->getGardAtGaussPtsPtr()));
// assemble residual
CHKERR opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr, Sev::inform);
};
// Push operators to the left hand side pipeline. Indicate that domain (i.e.
// volume/interior) element is used.
CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
// Push operators to the right hand side pipeline.
CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
}
//! [Push operators to pipeline]
/**
* @brief Monitor solution
*
* This functions is called by TS solver at the end of each step. It is used
* to output results to the hard drive.
*/
struct Monitor : public FEMethod {
std::pair<boost::shared_ptr<PostProcEleDomain>,
boost::shared_ptr<PostProcEleBdy>>
pair_post_proc_fe,
boost::shared_ptr<DomainEle> reaction_fe,
std::vector<boost::shared_ptr<ScalingMethod>> smv)
: dM(dm), reactionFE(reaction_fe), vecOfTimeScalingMethods(smv) {
fRes = createDMVector(dM);
domainPostProcFe = pair_post_proc_fe.first;
skinPostProcFe = pair_post_proc_fe.second;
}
MoFEMErrorCode postProcess() {
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every_nth_step",
&save_every_nth_step, PETSC_NULL);
auto make_vtk = [&]() {
if (domainPostProcFe) {
CHKERR DMoFEMLoopFiniteElements(dM, "dFE", domainPostProcFe,
getCacheWeakPtr());
CHKERR domainPostProcFe->writeFile(
"out_plastic_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
}
if (skinPostProcFe) {
CHKERR DMoFEMLoopFiniteElements(dM, "bFE", skinPostProcFe,
getCacheWeakPtr());
CHKERR skinPostProcFe->writeFile(
"out_skin_plastic_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
}
};
if (!(ts_step % save_every_nth_step)) {
CHKERR make_vtk();
}
if (reactionFE) {
CHKERR VecZeroEntries(fRes);
reactionFE->f = fRes;
CHKERR DMoFEMLoopFiniteElements(dM, "dFE", reactionFE, getCacheWeakPtr());
CHKERR VecAssemblyBegin(fRes);
CHKERR VecAssemblyEnd(fRes);
CHKERR VecGhostUpdateBegin(fRes, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(fRes, ADD_VALUES, SCATTER_REVERSE);
MoFEM::Interface *m_field_ptr;
CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
*m_field_ptr, reactionFE, fRes)();
double nrm;
CHKERR VecNorm(fRes, NORM_2, &nrm);
MOFEM_LOG("PlasticPrb", Sev::verbose)
<< "Residual norm " << nrm << " at time step " << ts_step;
}
auto get_min_max_displacement = [&]() {
MoFEM::Interface *m_field_ptr;
CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
std::array<double, 4> a_min = {DBL_MAX, DBL_MAX, DBL_MAX, 0};
std::array<double, 4> a_max = {-DBL_MAX, -DBL_MAX, -DBL_MAX, 0};
auto get_min_max = [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
int d = 0;
for (auto v : field_entity_ptr->getEntFieldData()) {
a_min[d] = std::min(a_min[d], v);
a_max[d] = std::max(a_max[d], v);
++d;
}
};
a_min[SPACE_DIM] = 0;
a_max[SPACE_DIM] = 0;
Range verts;
CHKERR m_field_ptr->get_field_entities_by_type("U", MBVERTEX, verts);
CHKERR m_field_ptr->getInterface<FieldBlas>()->fieldLambdaOnEntities(
get_min_max, "U", &verts);
auto mpi_reduce = [&](auto &a, auto op) {
std::array<double, 3> a_mpi = {0, 0, 0};
MPI_Allreduce(a.data(), a_mpi.data(), 3, MPI_DOUBLE, op,
m_field_ptr->get_comm());
return a_mpi;
};
auto a_min_mpi = mpi_reduce(a_min, MPI_MIN);
auto a_max_mpi = mpi_reduce(a_max, MPI_MAX);
MOFEM_LOG("PlasticPrb", Sev::inform)
<< "Min displacement " << a_min_mpi[0] << " " << a_min_mpi[1] << " "
<< a_min_mpi[2];
MOFEM_LOG("PlasticPrb", Sev::inform)
<< "Max displacement " << a_max_mpi[0] << " " << a_max_mpi[1] << " "
<< a_max_mpi[2];
};
CHKERR get_min_max_displacement();
double scale = 1;
for (auto s : vecOfTimeScalingMethods) {
scale *= s->getScale(this->ts_t);
}
MOFEM_LOG("PlasticPrb", Sev::inform)
<< "Time: " << this->ts_t << " scale: " << scale;
}
private:
boost::shared_ptr<PostProcEleDomain> domainPostProcFe;
boost::shared_ptr<PostProcEleBdy> skinPostProcFe;
boost::shared_ptr<FEMethod> reactionFE;
VecOfTimeScalingMethods vecOfTimeScalingMethods;
};
static boost::shared_ptr<TSUpdate> ts_update_ptr = nullptr;
//! [Solve]
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto dm = simple->getDM();
auto time_scale = boost::make_shared<TimeScale>();
// Setup postprocessing
auto create_post_proc_fe = [dm, this, simple]() {
//! [DG projection]
// Post-process domain, i.e. volume elements. If 3d case creates stresses
// are evaluated at points which are trace of the volume element on boundary
auto post_proc_ele_domain = [this](auto &pip_domain, auto &fe_name) {
// Push bases on reference element to the phyiscal element
// calculate displacement gradients at nodes of post processing mesh. For
// gradient DG projection is obsolete, since displacements can be
// evaluated at arbitrary points.
auto grad_ptr = boost::make_shared<MatrixDouble>();
pip_domain.push_back(
grad_ptr));
// Create operator to run (neted) pipeline in operator. InteriorElem
// element will have iteration rule and bases as domain element used to
// solve problem, and remember history variables.
using InteriorElem = DomainEle;
auto op_this = new OpLoopThis<InteriorElem>(mField, fe_name, Sev::noisy);
// add interior element to post-processing pipeline
pip_domain.push_back(op_this);
// get pointer to interior element, which is in essence pipeline which
// pipeline.
auto fe_physics = op_this->getThisFEPtr();
auto evaluate_stress_on_physical_element = [&]() {
// evaluate stress and plastic strain at Gauss points
fe_physics->getRuleHook = cpPtr->integrationRule;
fe_physics->getOpPtrVector(), {H1});
auto common_data_ptr =
boost::make_shared<ADOLCPlasticity::CommonData>();
// note that gradient of displacements is evaluate again, at
// physical nodes
fe_physics->getOpPtrVector().push_back(
"U", common_data_ptr->getGardAtGaussPtsPtr()));
CHKERR cpPtr->addMatBlockOps(mField, fe_physics->getOpPtrVector(),
"ADOLCMAT", Sev::noisy);
fe_physics->getOpPtrVector().push_back(
getRawPtrOpCalculateStress<SPACE_DIM>(mField, common_data_ptr,
cpPtr, false));
return common_data_ptr;
};
auto dg_projection_froward_and_back = [&](auto &common_data_ptr) {
// here we do actual projection of stress and plastic strain to DG space
int dg_order = approxOrder - 1;
auto entity_data_l2 =
boost::make_shared<EntitiesFieldData>(MBENTITYSET);
auto mass_ptr =
boost::make_shared<MatrixDouble>(); //< projection operator (mass
// matrix)
fe_physics->getOpPtrVector().push_back(new OpDGProjectionMassMatrix(
dg_order, mass_ptr, entity_data_l2, approxBase, L2));
auto coeffs_ptr_stress = boost::make_shared<MatrixDouble>();
// project stress on DG space on physical element
fe_physics->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
common_data_ptr->getStressMatrixPtr(), coeffs_ptr_stress, mass_ptr,
entity_data_l2, approxBase, L2));
// project strains plastic strains DG space on physical element
auto coeffs_ptr_plastic_strain = boost::make_shared<MatrixDouble>();
fe_physics->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
common_data_ptr->getPlasticStrainMatrixPtr(),
coeffs_ptr_plastic_strain, mass_ptr, entity_data_l2, approxBase,
L2));
// project stress and plastic strain for DG space on post-process
// element
pip_domain.push_back(new OpDGProjectionEvaluation(
common_data_ptr->getStressMatrixPtr(), coeffs_ptr_stress,
entity_data_l2, approxBase, L2));
pip_domain.push_back(new OpDGProjectionEvaluation(
common_data_ptr->getPlasticStrainMatrixPtr(),
coeffs_ptr_plastic_strain, entity_data_l2, approxBase, L2));
};
auto common_data_ptr = evaluate_stress_on_physical_element();
dg_projection_froward_and_back(common_data_ptr);
return boost::make_tuple(grad_ptr, common_data_ptr->getStressMatrixPtr(),
common_data_ptr->getPlasticStrainMatrixPtr());
};
//! [DG projection]
// Create tags on post-processing mesh, i.e. those tags are visible in
// Preview
auto add_post_proc_map = [&](auto post_proc_fe, auto u_ptr, auto grad_ptr,
auto stress_ptr, auto plastic_strain_ptr) {
post_proc_fe->getOpPtrVector().push_back(
new OpPPMapSPACE_DIM(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{},
{{"U", u_ptr}},
{{"GRAD", grad_ptr}},
{}
)
);
using OpPPMap3D = OpPostProcMapInMoab<3, 3>;
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap3D(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{},
{},
{},
{{"STRESS", stress_ptr}, {"PLASTIC_STRAIN", plastic_strain_ptr}}
)
);
return post_proc_fe;
};
auto vol_post_proc = [this, simple, post_proc_ele_domain,
add_post_proc_map]() {
PetscBool post_proc_vol = PETSC_FALSE;
if (SPACE_DIM == 2)
post_proc_vol = PETSC_TRUE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol",
&post_proc_vol, PETSC_NULL);
if (post_proc_vol == PETSC_FALSE)
return boost::shared_ptr<PostProcEleDomain>();
auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
auto u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
auto [grad_ptr, stress_ptr, plastic_strain_ptr] = post_proc_ele_domain(
post_proc_fe->getOpPtrVector(), simple->getDomainFEName());
return add_post_proc_map(post_proc_fe, u_ptr, grad_ptr, stress_ptr,
plastic_strain_ptr);
};
auto skin_post_proc = [this, simple, post_proc_ele_domain,
add_post_proc_map]() {
// create skin of the volume mesh for post-processing,
// i.e. boundary of the volume mesh
PetscBool post_proc_skin = PETSC_TRUE;
if (SPACE_DIM == 2)
post_proc_skin = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",
&post_proc_skin, PETSC_NULL);
if (post_proc_skin == PETSC_FALSE)
return boost::shared_ptr<PostProcEleBdy>();
auto post_proc_fe = boost::make_shared<PostProcEleBdy>(mField);
auto u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
auto op_loop_side =
new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
auto [grad_ptr, stress_ptr, plastic_strain_ptr] = post_proc_ele_domain(
op_loop_side->getOpPtrVector(), simple->getDomainFEName());
post_proc_fe->getOpPtrVector().push_back(op_loop_side);
return add_post_proc_map(post_proc_fe, u_ptr, grad_ptr, stress_ptr,
plastic_strain_ptr);
};
return std::make_pair(vol_post_proc(), skin_post_proc());
};
auto create_reaction_fe = [&]() {
auto fe_ptr = boost::make_shared<DomainEle>(mField);
fe_ptr->getRuleHook = cpPtr->integrationRule;
auto &pip = fe_ptr->getOpPtrVector();
auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
"U", common_data_ptr->getGardAtGaussPtsPtr()));
CHKERR opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr, Sev::noisy);
return fe_ptr;
};
auto add_extra_finite_elements_to_ksp_solver_pipelines = [&]() {
auto pre_proc_ptr = boost::make_shared<FEMethod>();
auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
auto get_bc_hook_rhs = [this, pre_proc_ptr, time_scale]() {
{time_scale}, false)();
};
pre_proc_ptr->preProcessHook = get_bc_hook_rhs;
auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
mField, post_proc_rhs_ptr, 1.)();
};
auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
mField, post_proc_lhs_ptr, 1.)();
};
post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs;
// This is low level pushing finite elements (pipelines) to solver
auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
};
// Add extra finite elements to SNES solver pipelines to resolve essential
// boundary conditions
CHKERR add_extra_finite_elements_to_ksp_solver_pipelines();
auto create_monitor_fe = [dm, time_scale](auto &&post_proc_fe,
auto &&reaction_fe) {
return boost::make_shared<Monitor>(
dm, post_proc_fe, reaction_fe,
std::vector<boost::shared_ptr<ScalingMethod>>{time_scale});
};
//! [Set up post-step]
auto set_up_post_step = [&](auto ts) {
// create finite element (pipeline) and populate it with operators to
// update history variables
auto create_update_ptr = [&]() {
auto fe_ptr = boost::make_shared<DomainEle>(mField);
fe_ptr->getRuleHook = cpPtr->integrationRule;
{H1});
auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
fe_ptr->getOpPtrVector().push_back(
"U", common_data_ptr->getGardAtGaussPtsPtr()));
opFactoryDomainUpdate<SPACE_DIM>(mField, fe_ptr->getOpPtrVector(),
"ADOLCMAT", common_data_ptr, cpPtr),
"push update ops");
return fe_ptr;
};
// ts_update_ptr is global static variable
createTSUpdate(simple->getDomainFEName(), create_update_ptr());
//! [TS post-step function]
// This is pure "C" function which we can to the TS solver
auto ts_step_post_proc = [](TS ts) {
CHKERR ts_update_ptr->postProcess(ts);
};
//! [TS post-step function]
// finally set up post-step
CHKERR TSSetPostStep(ts, ts_step_post_proc);
};
//! [Set up post-step]
// Set monitor which postprocessing results and saves them to the hard drive
auto set_up_monitor = [&](auto ts) {
boost::shared_ptr<FEMethod> null_fe;
auto monitor_ptr =
create_monitor_fe(create_post_proc_fe(), create_reaction_fe());
CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
null_fe, monitor_ptr);
};
auto set_up_adapt = [&](auto ts) {
TSAdapt adapt;
CHKERR TSGetAdapt(ts, &adapt);
};
//! [Create TS]
auto ts = pipeline_mng->createTSIM();
// Set time solver
double ftime = 1;
CHKERR TSSetMaxTime(ts, ftime);
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
auto D = createDMVector(simple->getDM());
CHKERR TSSetSolution(ts, D);
// Set monitor, step adaptivity, and post-step to update history variables
CHKERR set_up_monitor(ts);
CHKERR set_up_post_step(ts);
CHKERR set_up_adapt(ts);
CHKERR TSSetFromOptions(ts);
CHKERR TSSolve(ts, NULL);
//! [Create TS]
CHKERR TSGetTime(ts, &ftime);
PetscInt steps, snesfails, rejects, nonlinits, linits;
CHKERR TSGetStepNumber(ts, &steps);
CHKERR TSGetSNESFailures(ts, &snesfails);
CHKERR TSGetStepRejections(ts, &rejects);
CHKERR TSGetSNESIterations(ts, &nonlinits);
CHKERR TSGetKSPIterations(ts, &linits);
MOFEM_LOG_C("PlasticPrb", Sev::inform,
"steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
"%d, linits %d",
steps, rejects, snesfails, ftime, nonlinits, linits);
}
//! [Solve]
//! [Getting norms]
auto dm = simple->getDM();
auto T = createDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
SCATTER_FORWARD);
double nrm2;
CHKERR VecNorm(T, NORM_2, &nrm2);
MOFEM_LOG("PlasticPrb", Sev::inform) << "Solution norm " << nrm2;
auto post_proc_norm_fe = boost::make_shared<DomainEle>(mField);
post_proc_norm_fe->getRuleHook = cpPtr->integrationRule;
post_proc_norm_fe->getOpPtrVector(), {H1});
enum NORMS { U_NORM_L2 = 0, PIOLA_NORM, LAST_NORM };
auto norms_vec =
(mField.get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
CHKERR VecZeroEntries(norms_vec);
auto u_ptr = boost::make_shared<MatrixDouble>();
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
new OpCalcNormL2Tensor1<SPACE_DIM>(u_ptr, norms_vec, U_NORM_L2));
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
post_proc_norm_fe);
CHKERR VecAssemblyBegin(norms_vec);
CHKERR VecAssemblyEnd(norms_vec);
MOFEM_LOG_CHANNEL("SELF"); // Clear channel from old tags
if (mField.get_comm_rank() == 0) {
const double *norms;
CHKERR VecGetArrayRead(norms_vec, &norms);
MOFEM_TAG_AND_LOG("SELF", Sev::inform, "example")
<< "norm_u: " << std::scientific << std::sqrt(norms[U_NORM_L2]);
CHKERR VecRestoreArrayRead(norms_vec, &norms);
}
}
//! [Getting norms]
//! [Postprocessing results]
PetscInt test_nb = 0;
PetscBool test_flg = PETSC_FALSE;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-test", &test_nb, &test_flg);
if (test_flg) {
auto T = createDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
SCATTER_FORWARD);
double nrm2;
CHKERR VecNorm(T, NORM_2, &nrm2);
MOFEM_LOG("PlasticPrb", Sev::verbose) << "Regression norm " << nrm2;
double regression_value = 0;
switch (test_nb) {
default:
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong test number.");
break;
}
if (fabs(nrm2 - regression_value) > 1e-2)
SETERRQ2(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
"Regression test field; wrong norm value. %6.4e != %6.4e", nrm2,
regression_value);
}
}
//! [Postprocessing results]
//! [Check]
}
//! [Check]
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "PlasticPrb"));
LogManager::setLog("PlasticPrb");
MOFEM_LOG_TAG("PlasticPrb", "PlasticPrb");
MOFEM_LOG("PlasticPrb", Sev::inform) << "SPACE_DIM " << SPACE_DIM;
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database interface
//! [Create MoFEM]
//! [PlasticProblem]
PlasticProblem ex(m_field);
CHKERR ex.runProblem();
//! [PlasticProblem]
}
}
/** \file ADOLCPlasticityMaterialModels.hpp
* \ingroup adoc_plasticity
* \brief Matetial models for plasticity
*
* \example ADOLCPlasticityMaterialModels.hpp
*
*/
#ifndef __ADOLC_MATERIAL_MODELS_HPP
#define __ADOLC_MATERIAL_MODELS_HPP
#ifndef WITH_ADOL_C
#error "MoFEM need to be compiled with ADOL-C"
#endif
namespace ADOLCPlasticity {
/** \brief J2 plasticity (Kinematic Isotropic (Linear) Hardening)
* \ingroup adoc_plasticity
*/
/**
* @brief J2 plasticity (Kinematic Isotropic (Linear) Hardening)
*
* @tparam DIM model dimension (2 - is plane stress)
*/
template <int DIM> struct J2Plasticity;
/**
* @brief J2 (Von Misses) plasticity
* @ingroup adoc_plasticity
*
* Free energy:
*
* \f[
* \psi(\varepsilon^e, \alpha, \beta) =
* \frac{1}{2} \lambda \textrm{tr}[\varepsilon^e]^2
* + \mu \varepsilon^e : \varepsilon^e
* +
* \sigma_y\alpha
* +
* \frac{1}{2} \phi H \alpha^2
* +
* \frac{1}{2} (1-\phi)K \beta^2
* \f]
*
* Isotropic hardening variable \f$\alpha\f$ is in first index of
* internal variables vector. Kinematic hardening variable is in
* the remaining indices of internal variables vector.
*
* Yield function:
*
* Deviator of actual stress:
* \f[
* \eta=\textrm{dev}[\sigma]-\overline{\beta}
* \f]
* where \f$\overline{\beta}\f$ is back stress.
*
* Square of the deviator norm
* \f[
* f = \frac{3}{2} \eta:\eta
* \f]
*
* Yield function:
* \f[
* y = \sqrt{f} - \overline{\alpha}
* \f]
* where \f$f\f$ is defined in \ref eva
*
* This is associated potential model, so flow potential is equal to yield
*
*/
template <> struct J2Plasticity<3> : public ClosestPointProjection {
J2Plasticity() : ClosestPointProjection() {
nbInternalVariables = 7;
activeVariablesW.resize(12 + nbInternalVariables, false);
}
adouble tR;
double mu; //< Lamé parameter
double lambda; //< Lamé parameters
double H; //< Isotropic hardening
double K; //< Kinematic hardening
double phi; //< Combined isotropic/kinematic hardening
double sigmaY; //< Yield stress
enum ModelParams { LAMBDA, MU, ISOH, KINK, PHI, SIGMA_Y, LAST_PARAM };
using ParamsArray = std::array<double, LAST_PARAM>;
boost::shared_ptr<ParamsArray> defaultParamsArrayPtr = nullptr;
boost::shared_ptr<ParamsArray> paramsArrayPtr = nullptr;
boost::shared_ptr<ParamsArray> oldParamsArrayPtr = nullptr;
MoFEMErrorCode getDefaultMaterialParameters() {
double young_modulus = 1;
double poisson_ratio = 0.25;
sigmaY = 1;
H = 0.1;
K = 0;
phi = 1;
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-young_modulus", &young_modulus,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-poisson_ratio", &poisson_ratio,
0);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-sigma_y", &sigmaY, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-H", &H, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-K", &K, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-phi", &phi, PETSC_NULL);
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "Young modulus: " << young_modulus;
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "Poisson ratio: " << poisson_ratio;
defaultParamsArrayPtr = boost::make_shared<ParamsArray>();
(*defaultParamsArrayPtr)[LAMBDA] = lambda;
(*defaultParamsArrayPtr)[MU] = mu;
(*defaultParamsArrayPtr)[H] = H;
(*defaultParamsArrayPtr)[K] = K;
(*defaultParamsArrayPtr)[PHI] = phi;
(*defaultParamsArrayPtr)[SIGMA_Y] = sigmaY;
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "LAMBDA: " << (*defaultParamsArrayPtr)[LAMBDA];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "MU: " << (*defaultParamsArrayPtr)[MU];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "H: " << (*defaultParamsArrayPtr)[ISOH];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "K: " << (*defaultParamsArrayPtr)[KINK];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "PHI: " << (*defaultParamsArrayPtr)[PHI];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "SIGMA_Y: " << (*defaultParamsArrayPtr)[SIGMA_Y];
paramsArrayPtr = boost::make_shared<ParamsArray>();
std::copy(defaultParamsArrayPtr->begin(), defaultParamsArrayPtr->end(),
paramsArrayPtr->begin());
oldParamsArrayPtr = boost::make_shared<ParamsArray>();
std::fill(oldParamsArrayPtr->begin(), oldParamsArrayPtr->end(), 0);
};
/**
* @brief Get material parameters from block
*
* @param m_field
* @param pip
* @param block_name
* @param sev
* @return MoFEMErrorCode
*/
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string block_name, Sev sev) {
struct OpMatBlocks : public ForcesAndSourcesCore::UserDataOperator {
OpMatBlocks(boost::shared_ptr<ParamsArray> params_array_ptr,
boost::shared_ptr<ParamsArray> def_params_array_otr,
MoFEM::Interface &m_field, Sev sev,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
: OP(NOSPACE, OP::OPSPACE), paramsArrayPtr(params_array_ptr),
defParamsArray(def_params_array_otr) {
CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
"Can not get data from block");
}
MoFEMErrorCode doWork(int side, EntityType type,
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
std::copy(b.paramsArray.begin(), b.paramsArray.end(),
paramsArrayPtr->begin());
}
}
std::copy(defParamsArray->begin(), defParamsArray->end(),
paramsArrayPtr->begin());
}
private:
boost::shared_ptr<ParamsArray> paramsArrayPtr;
boost::shared_ptr<ParamsArray> defParamsArray;
struct BlockData {
ParamsArray paramsArray;
Range blockEnts;
};
std::vector<BlockData> blockData;
extractBlockData(MoFEM::Interface &m_field,
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
for (auto m : meshset_vec_ptr) {
MOFEM_TAG_AND_LOG("ADOLCPlasticityWold", sev, "JP2 parameters") << *m;
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() != 6) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Expected that block has two attribute");
}
auto get_block_ents = [&]() {
Range ents;
m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
return ents;
};
blockData.push_back({*defParamsArray, get_block_ents()});
std::copy(block_data.begin(), block_data.end(),
blockData.back().paramsArray.begin());
const auto young_modulus = blockData.back().paramsArray[LAMBDA];
const auto poisson_ratio = blockData.back().paramsArray[MU];
// Is assumed that uset provide young_modulus and poisson_ratio in
// fist two argiumens of the block
MOFEM_LOG("ADOLCPlasticityWold", sev)
<< "Young modulus: " << (blockData.back().paramsArray)[LAMBDA];
MOFEM_LOG("ADOLCPlasticityWold", sev)
<< "Poisson ratio: " << (blockData.back().paramsArray)[MU];
blockData.back().paramsArray[LAMBDA] =
blockData.back().paramsArray[MU] = MU(young_modulus, poisson_ratio);
MOFEM_LOG("ADOLCPlasticityWold", sev)
<< "LAMBDA: " << (blockData.back().paramsArray)[LAMBDA];
MOFEM_LOG("ADOLCPlasticityWold", sev)
<< "MU: " << (blockData.back().paramsArray)[MU];
MOFEM_LOG("ADOLCPlasticityWold", sev)
<< "H: " << (blockData.back().paramsArray)[ISOH];
MOFEM_LOG("ADOLCPlasticityWold", sev)
<< "K: " << (blockData.back().paramsArray)[KINK];
MOFEM_LOG("ADOLCPlasticityWold", sev)
<< "PHI: " << (blockData.back().paramsArray)[PHI];
MOFEM_LOG("ADOLCPlasticityWold", sev)
<< "SIGMA_Y: " << (blockData.back().paramsArray)[SIGMA_Y];
}
}
};
auto cubit_meshsets_vec_ptr =
m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
std::regex((boost::format("%s(.*)") % block_name).str()));
pip.push_back(new OpMatBlocks(paramsArrayPtr, defaultParamsArrayPtr,
m_field, sev, cubit_meshsets_vec_ptr));
}
/**
* @brief Set material parameters at integration point
*
* @param tag
* @param recalculate_elastic_tangent
* @return MoFEMErrorCode
*/
MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent) {
switch (tag) {
case TypesTags::W: {
set_param_vec(tapesTags[tag], paramsArrayPtr->size(),
paramsArrayPtr->data());
for (auto p = 0; p != LAST_PARAM; p++) {
if ((*paramsArrayPtr)[p] != (*oldParamsArrayPtr)[p]) {
recalculate_elastic_tangent = true;
std::copy(paramsArrayPtr->begin(), paramsArrayPtr->end(),
oldParamsArrayPtr->begin());
break;
}
}
} break;
case TypesTags::Y:
break;
default:
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Unknown tag for setParams");
}
}
//! [Free energy J2]
MoFEMErrorCode freeHemholtzEnergy() {
// ADOL-C register parameters (those can varied for integration points)
auto p_lambda = mkparam(lambda); // 0
auto p_mu = mkparam(mu); // 1
auto p_H = mkparam(H); // 2
auto p_K = mkparam(K); // 3
auto p_phi = mkparam(phi); // 4
auto p_sigma_y = mkparam(sigmaY); // 5
// Operator converting Voight and tensor notation
auto t_voight_strain_op = voight_to_strain_op();
// Variable in voight notation
auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
auto t_voight_internal_tensor =
getFTensor1FromPtr<6>(&a_internalVariables[1]);
// Convert variables to tensor notation
tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
tPlasticStrain(i, j) =
t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
tInternalTensor(i, j) =
t_voight_strain_op(i, j, Z) * t_voight_internal_tensor(Z);
// Evaluate elastic strain
tElasticStrain(i, j) = tTotalStrain(i, j) - tPlasticStrain(i, j);
tR = tElasticStrain(i, i);
// Strain energy
a_w = 0.5 * p_lambda * tR * tR +
p_mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
// Isotropic strain energy
a_w +=
a_internalVariables[0] * p_sigma_y +
0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
// Kinematic strain energy
a_w += (0.5 * (1 - p_phi) * p_K) *
(tInternalTensor(i, j) * tInternalTensor(i, j));
}
//! [Free energy J2]
//! [Yield function J2]
MoFEMErrorCode evalF() {
auto t_voight_stress_op = voight_to_stress_op();
auto t_vioght_stress = getFTensor1FromPtr<6>(&a_sTress[0]);
auto t_back_stress = getFTensor1FromPtr<6>(&a_internalFluxes[1]);
tStress(i, j) = t_voight_stress_op(i, j, Z) * t_vioght_stress(Z);
tBackStress(i, j) = t_voight_stress_op(i, j, Z) * t_back_stress(Z);
constexpr double third = boost::math::constants::third<double>();
tR = tStress(i, i) / 3.;
tDeviator(i, j) = tStress(i, j) - tR * t_kd(i, j) - tBackStress(i, j);
// Fix the constant to mach uniaxial test
constexpr double c = 3. / 2.;
f = c * tDeviator(i, j) * tDeviator(i, j);
}
MoFEMErrorCode yieldFunction() {
CHKERR evalF();
a_y = sqrt(f) - a_internalFluxes[0];
}
//! [Yield function J2]
MoFEMErrorCode flowPotential() {
CHKERR evalF();
a_h = sqrt(f) - a_internalFluxes[0];
}
};
//! [J2 2D]
template <> struct J2Plasticity<2> : public J2Plasticity<3> {
MoFEMErrorCode freeHemholtzEnergy() {
auto p_lambda = mkparam(lambda); // 0
auto p_mu = mkparam(mu); // 1
auto p_H = mkparam(H); // 2
auto p_K = mkparam(K); // 3
auto p_phi = mkparam(phi); // 4
auto p_sigma_y = mkparam(sigmaY); // 5
auto t_voight_strain_op = voight_to_strain_op();
auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
auto t_voight_internal_tensor =
getFTensor1FromPtr<6>(&a_internalVariables[1]);
tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
tPlasticStrain(i, j) =
t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
tInternalTensor(i, j) =
t_voight_strain_op(i, j, Z) * t_voight_internal_tensor(Z);
tElasticStrain(i, j) = tTotalStrain(i, j) - tPlasticStrain(i, j);
//! [Plane stress]
// Calculate st rain at ezz enforcing plane strain case
tElasticStrain(2, 2) =
(-tElasticStrain(0, 0) * p_lambda - tElasticStrain(1, 1) * p_lambda) /
(p_lambda + 2 * p_mu);
//! [Plane stress]
tR = tElasticStrain(i, i);
a_w = 0.5 * p_lambda * tR * tR +
p_mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
// plastic part
// isotropic
a_w +=
a_internalVariables[0] * p_sigma_y +
0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
// kinematic
a_w += (0.5 * (1 - p_phi) * p_K) *
(tInternalTensor(i, j) * tInternalTensor(i, j));
}
};
//! [J2 2D]
/** \brief Paraboloidal yield criterion
*
* See paper for reference \cite ullah2017three and \cite stier2015comparing
*
* Linear hardening
* \f[
* \psi =
* \frac{1}{2} \lambda \textrm{tr}[\varepsilon]^2 + \mu \varepsilon :
* \varepsilon
* +
* \sigma_{yt}\alpha_0
* +
* \frac{1}{2} H_t \alpha_0^2
* +
* \sigma_{yc}\alpha_1
* +
* \frac{1}{2} H_c \alpha_1^2
* \f]
*
* Exponential hardening
* \f[
* \psi =
* \frac{1}{2} \lambda \textrm{tr}[\varepsilon]^2 + \mu \varepsilon :
* \varepsilon\\
* +
* (\sigma_{yt}+H_t)\alpha_0
* +
* \frac{H_t}{n_t} \exp{(-n_t \alpha_0)}\\
* +
* (\sigma_{yc}+H_c)\alpha_0
* +
* \frac{H_c}{n_c} \exp{(-n_c \alpha_1)}
* \f]
*
* \f[
* I_1 = \textrm{tr} (\boldsymbol{\sigma})
* \f]
*
* \f[
* \eta=\textrm{dev}[\boldsymbol{\sigma}]
* \f]
*
* \f[
* J_2 = \frac{1}{2} \eta:\eta
* \f]
*
* \f[
* y = 6J_2 + 2I_1\left(\overline{\alpha_1}-\overline{\alpha_0}\right) -
* 2\overline{\alpha_0} \,\overline{\alpha_1} \f]
* where
* \f[
* \overline{\alpha_0}=\frac{\partial \psi}{\partial \alpha_0}=\sigma_{yt} + H_t
* \alpha_0
* \f]
*
* \f[
* \overline{\alpha_1}=\frac{\partial \psi}{\partial \alpha_1}=\sigma_{yc} + H_c
* \alpha_1
* \f]
*
* \f[
* \Psi = 6J_2 + 2\alpha I_1 \left(\overline{\alpha_1}-\overline{\alpha_0}\right)
* - 2\overline{\alpha_0} \,\overline{\alpha_1} \f] \f[ \alpha=
* \frac{1-2\nu_p}{1+\nu_p}
* \f]
*
* \ingroup adoc_plasticity
*/
struct ParaboloidalPlasticity : public ClosestPointProjection {
}
double lambda; //< Lamé parameters
double mu; //< Lamé parameter
double nup; //< Plastic Poisson’s ratio
double Ht, Hc; //< Isotropic hardening for tension and compression
double sigmaYt, sigmaYc; //< Yield stress in tension and compression
double nt, nc; //< Hardening parameters for tension and compression
enum ModelParams {
MU,
NUP,
HT,
HC,
NT,
NC,
};
using ParamsArray = std::array<double, LAST_PARAM>;
boost::shared_ptr<ParamsArray> defaultParamsArrayPtr = nullptr;
double young_modulus = 1;
double poisson_ratio = 0.25;
sigmaYt = 1;
sigmaYc = 1;
Ht = 0.1;
Hc = 0.1;
nc = 1.;
nt = 1.;
nup = 0.;
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-young_modulus", &young_modulus,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-poisson_ratio", &poisson_ratio,
0);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-sigma_yt", &sigmaYt, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-sigma_yc", &sigmaYc, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-Ht", &Ht, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-Hc", &Hc, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-nt", &nt, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-nc", &nc, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-nup", &nup, PETSC_NULL);
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "Young modulus: " << young_modulus;
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "Poisson ratio: " << poisson_ratio;
defaultParamsArrayPtr = boost::make_shared<ParamsArray>();
(*defaultParamsArrayPtr)[LAMBDA] = lambda;
(*defaultParamsArrayPtr)[MU] = mu;
(*defaultParamsArrayPtr)[NUP] = nup;
(*defaultParamsArrayPtr)[SIGMA_YT] = sigmaYt;
(*defaultParamsArrayPtr)[SIGMA_YC] = sigmaYc;
(*defaultParamsArrayPtr)[HT] = Ht;
(*defaultParamsArrayPtr)[HC] = Hc;
(*defaultParamsArrayPtr)[NT] = nt;
(*defaultParamsArrayPtr)[NC] = nc;
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "LAMBDA: " << (*defaultParamsArrayPtr)[LAMBDA];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "MU: " << (*defaultParamsArrayPtr)[MU];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "NUP: " << (*defaultParamsArrayPtr)[NUP];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "SIGMA_YT: " << (*defaultParamsArrayPtr)[SIGMA_YT];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "SIGMA_YC: " << (*defaultParamsArrayPtr)[SIGMA_YC];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "HT: " << (*defaultParamsArrayPtr)[HT];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "HC: " << (*defaultParamsArrayPtr)[HC];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "NT: " << (*defaultParamsArrayPtr)[NT];
MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
<< "NC: " << (*defaultParamsArrayPtr)[NC];
};
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string block_name, Sev sev) {
return 0;
}
MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent) {
recalculate_elastic_tangent = true;
return 0;
}
//! [Paraboloidal free energy]
auto t_voight_strain_op = voight_to_strain_op();
auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
a_w = 0.5 * lambda * tR * tR +
}
//! [Paraboloidal free energy]
//! [Paraboloidal yield function]
auto t_voight_stress_op = voight_to_stress_op();
auto t_vioght_stress = getFTensor1FromPtr<6>(&a_sTress[0]);
tStress(i, j) = t_voight_stress_op(i, j, Z) * t_vioght_stress(Z);
tR = tStress(i, i);
I1 = tR;
tR /= 3.;
tDeviator(i, j) = tStress(i, j) - tR * t_kd(i, j);
J2 = 0.5 * (tDeviator(i, j) * tDeviator(i, j));
}
a_y = 6 * J2 + 2 * I1 * (a_internalFluxes[1] - a_internalFluxes[0]) -
}
//! [Paraboloidal yield function]
//! [Paraboloidal flow potential]
double alpha =
(1 - 2 * nup) /
(1 + nup); // relation between alpha and plastic poisson's ratio
a_h = 6 * J2 +
2 * alpha * I1 * (a_internalFluxes[1] - a_internalFluxes[0]) -
}
//! [Paraboloidal flow potential]
};
template <typename MODEL>
inline auto createMaterial(
std::array<int, ClosestPointProjection::LAST_TAPE> tapes_tags = {0, 1, 2}) {
auto cp_ptr = boost::make_shared<MODEL>();
cp_ptr->getDefaultMaterialParameters();
cp_ptr->createMatAVecR();
cp_ptr->snesCreate();
cp_ptr->recordTapes();
cp_ptr->tapesTags = tapes_tags;
return cp_ptr;
}
} // namespace ADOLCPlasticity
#endif //__ADOLC_MATERIAL_MODELS_HPP
PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition: adolc_plasticity.cpp:97
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
ADOLCPlasticity::ClosestPointProjection::a_h
adouble a_h
Definition: ADOLCPlasticity.hpp:330
ADOLCPlasticity::getFTensor2SymmetricDiffBaseImpl
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 6 >, 3 > getFTensor2SymmetricDiffBaseImpl(DataForcesAndSourcesCore::EntData &data, MatrixDouble &storage, const bool b_bar, const int nb_integration_pts, double *w_ptr, FTensor::Number< DIM >)
[BBar method]
Definition: ADOLCPlasticity.cpp:85
PlasticProblem::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: adolc_plasticity.cpp:463
MoFEM::EssentialPreProcReaction< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:151
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
ADOLCPlasticity::ParaboloidalPlasticity::defaultParamsArrayPtr
boost::shared_ptr< ParamsArray > defaultParamsArrayPtr
Definition: ADOLCPlasticityMaterialModels.hpp:571
sdf_hertz.d
float d
Definition: sdf_hertz.py:5
LASTBASE
@ LASTBASE
Definition: definitions.h:69
MoFEM::OpLoopThis
Execute "this" element in the operator.
Definition: DGProjection.hpp:68
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
ADOLCPlasticity::createTSUpdate
boost::shared_ptr< TSUpdate > createTSUpdate(std::string fe_name, boost::shared_ptr< FEMethod > fe_ptr)
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:172
ADOLCPlasticity::ParaboloidalPlasticity::tElasticStrain
FTensor::Tensor2_symmetric< adouble, 3 > tElasticStrain
Definition: ADOLCPlasticityMaterialModels.hpp:544
ADOLCPlasticity::ParaboloidalPlasticity::SIGMA_YC
@ SIGMA_YC
Definition: ADOLCPlasticityMaterialModels.hpp:562
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
ADOLCPlasticity::ClosestPointProjection::ClosestPointProjection
ClosestPointProjection()
Definition: ClosetPointProjection.cpp:15
sigmaY
double sigmaY
Yield stress.
Definition: plastic.cpp:174
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
MoFEM::EssentialPostProcLhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:130
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::NaturalBC::Assembly::BiLinearForm
Definition: Natural.hpp:74
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
ADOLCPlasticity::ClosestPointProjection::a_plasticStrain
ublas::vector< adouble > a_plasticStrain
Definition: ADOLCPlasticity.hpp:325
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
ADOLCPlasticity::ParaboloidalPlasticity::HT
@ HT
Definition: ADOLCPlasticityMaterialModels.hpp:563
ADOLCPlasticity::ParaboloidalPlasticity::nt
double nt
Definition: ADOLCPlasticityMaterialModels.hpp:555
ADOLCPlasticity::ParaboloidalPlasticity::tR
adouble tR
Definition: ADOLCPlasticityMaterialModels.hpp:548
ADOLCPlasticity::ParaboloidalPlasticity::j
FTensor::Index< 'j', 3 > j
Definition: ADOLCPlasticityMaterialModels.hpp:539
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:29
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
ADOLCPlasticity::ParaboloidalPlasticity::Ht
double Ht
Definition: ADOLCPlasticityMaterialModels.hpp:553
ADOLCPlasticity::ParaboloidalPlasticity::sigmaYt
double sigmaYt
Definition: ADOLCPlasticityMaterialModels.hpp:554
ADOLCPlasticity::ParaboloidalPlasticity::Hc
double Hc
Definition: ADOLCPlasticityMaterialModels.hpp:553
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
ts_update_ptr
static boost::shared_ptr< TSUpdate > ts_update_ptr
Definition: adolc_plasticity.cpp:460
ADOLCPlasticity::ClosestPointProjection::a_sTrain
ublas::vector< adouble > a_sTrain
Definition: ADOLCPlasticity.hpp:324
MoFEM::EssentialPostProcRhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:111
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:104
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
ADOLCPlasticity::ParaboloidalPlasticity::NUP
@ NUP
Definition: ADOLCPlasticityMaterialModels.hpp:560
MoFEM.hpp
MoFEM::OpDGProjectionCoefficients
Definition: DGProjection.hpp:119
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:527
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
SideEle
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
Definition: adolc_plasticity.cpp:98
W
double W
Definition: free_surface.cpp:166
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::getFTensor2SymmetricFromPtr< 3 >
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 6 >, 3 > getFTensor2SymmetricFromPtr< 3 >(double *ptr)
Definition: Templates.hpp:998
ADOLCPlasticity::ParaboloidalPlasticity::flowPotential
MoFEMErrorCode flowPotential()
[Paraboloidal yield function]
Definition: ADOLCPlasticityMaterialModels.hpp:705
ADOLCPlasticity::ParaboloidalPlasticity::mu
double mu
Definition: ADOLCPlasticityMaterialModels.hpp:551
FTensor::Tensor2< double, 3, 3 >
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ADOLCPlasticity::ParaboloidalPlasticity::NT
@ NT
Definition: ADOLCPlasticityMaterialModels.hpp:565
MoFEM::OpFluxLhsImpl
Definition: Natural.hpp:43
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
PostProcEleBdy
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
Definition: adolc_plasticity.cpp:99
PlasticProblem
Definition: adolc_plasticity.cpp:105
ADOLCPlasticity::ParaboloidalPlasticity::ModelParams
ModelParams
Definition: ADOLCPlasticityMaterialModels.hpp:557
PlasticProblem::outputResults
MoFEMErrorCode outputResults()
[Getting norms]
Definition: adolc_plasticity.cpp:876
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
PostProcEleByDim
Definition: adolc_plasticity.cpp:82
MoFEM::Sev
MoFEM::LogManager::SeverityLevel Sev
Definition: CoreTemplates.hpp:17
ADOLCPlasticity::ParaboloidalPlasticity::nc
double nc
Definition: ADOLCPlasticityMaterialModels.hpp:555
PlasticProblem::gettingNorms
MoFEMErrorCode gettingNorms()
[Solve]
Definition: adolc_plasticity.cpp:823
save_every_nth_step
int save_every_nth_step
Definition: photon_diffusion.cpp:67
FTensor::Number
Definition: Number.hpp:11
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
ADOLCPlasticity::J2Plasticity< 3 >::J2Plasticity
J2Plasticity()
Definition: ADOLCPlasticityMaterialModels.hpp:75
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1018
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
ADOLCPlasticity::ParaboloidalPlasticity::getDefaultMaterialParameters
MoFEMErrorCode getDefaultMaterialParameters()
Definition: ADOLCPlasticityMaterialModels.hpp:573
BlockData
Definition: gel_analysis.cpp:74
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
PlasticProblem::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: adolc_plasticity.cpp:160
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
ADOLCPlasticity::ParaboloidalPlasticity::evalF
MoFEMErrorCode evalF()
[Paraboloidal yield function]
Definition: ADOLCPlasticityMaterialModels.hpp:681
ADOLCPlasticity::ParaboloidalPlasticity::I1
adouble I1
[Paraboloidal free energy]
Definition: ADOLCPlasticityMaterialModels.hpp:678
convert.type
type
Definition: convert.py:64
ADOLCPlasticity::ParaboloidalPlasticity::i
FTensor::Index< 'i', 3 > i
Definition: ADOLCPlasticityMaterialModels.hpp:538
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:173
ADOLCPlasticity::createMaterial
auto createMaterial(std::array< int, ClosestPointProjection::LAST_TAPE > tapes_tags={0, 1, 2})
Definition: ADOLCPlasticityMaterialModels.hpp:721
ADOLCPlasticity::voight_to_strain_op
FTensor::Dg< double, 3, 6 > voight_to_strain_op()
Op convert Vight strain vector to strain tensor.
Definition: ADOLCPlasticity.hpp:29
ADOLCPlasticity::ParaboloidalPlasticity::Z
FTensor::Index< 'Z', 6 > Z
Definition: ADOLCPlasticityMaterialModels.hpp:540
MatrixFunction.hpp
ADOLCPlasticity::ParaboloidalPlasticity::ParaboloidalPlasticity
ParaboloidalPlasticity()
Definition: ADOLCPlasticityMaterialModels.hpp:533
MoFEM::AddFluxToLhsPipelineImpl
Definition: Natural.hpp:49
ADOLCPlasticity::ParaboloidalPlasticity::tPlasticStrain
FTensor::Tensor2_symmetric< adouble, 3 > tPlasticStrain
Definition: ADOLCPlasticityMaterialModels.hpp:543
MoFEM::OpCalcNormL2Tensor1
Get norm of input MatrixDouble for Tensor1.
Definition: NormsOperators.hpp:42
TSADAPTMOFEM
#define TSADAPTMOFEM
Definition: TsCtx.hpp:10
ADOLCPlasticity::ParaboloidalPlasticity::addMatBlockOps
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, Sev sev)
Get model parameters from blocks.
Definition: ADOLCPlasticityMaterialModels.hpp:641
ADOLCPlasticity::ClosestPointProjection::a_w
adouble a_w
Definition: ADOLCPlasticity.hpp:328
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
ADOLCPlasticity
Definition: ADOLCPlasticity.hpp:24
PlasticProblem::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition: adolc_plasticity.cpp:145
ADOLCPlasticity::ParaboloidalPlasticity::nup
double nup
Definition: ADOLCPlasticityMaterialModels.hpp:552
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
ADOLCPlasticity::ClosestPointProjection::a_internalFluxes
ublas::vector< adouble > a_internalFluxes
Definition: ADOLCPlasticity.hpp:327
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
main
int main(int argc, char *argv[])
Definition: activate_deactivate_dofs.cpp:15
ADOLCPlasticity::ClosestPointProjection::nbInternalVariables
int nbInternalVariables
Definition: ADOLCPlasticity.hpp:141
MoFEM::PostProcBrokenMeshInMoabBase
Definition: PostProcBrokenMeshInMoabBase.hpp:94
ADOLCPlasticity::ClosestPointProjection::activeVariablesW
VectorAdaptor activeVariablesW
Definition: ADOLCPlasticity.hpp:173
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index
Definition: Index.hpp:23
PlasticProblem::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: adolc_plasticity.cpp:241
PlasticProblem::assembleSystem
MoFEMErrorCode assembleSystem()
[Boundary condition]
Definition: adolc_plasticity.cpp:281
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
ADOLCPlasticity::ParaboloidalPlasticity::tDeviator
FTensor::Tensor2_symmetric< adouble, 3 > tDeviator
Definition: ADOLCPlasticityMaterialModels.hpp:546
ADOLCPlasticity::ParaboloidalPlasticity::LAST_PARAM
@ LAST_PARAM
Definition: ADOLCPlasticityMaterialModels.hpp:567
PlasticProblem::checkResults
MoFEMErrorCode checkResults()
[Postprocessing results]
Definition: adolc_plasticity.cpp:906
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
DomainEleOp
adouble
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::OpDGProjectionEvaluation
Definition: DGProjection.hpp:136
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
ADOLCPlasticity::ParaboloidalPlasticity::tTotalStrain
FTensor::Tensor2_symmetric< adouble, 3 > tTotalStrain
Definition: ADOLCPlasticityMaterialModels.hpp:542
ADOLCPlasticity::ClosestPointProjection::a_internalVariables
ublas::vector< adouble > a_internalVariables
Definition: ADOLCPlasticity.hpp:326
MoFEM::DMoFEMGetInterfacePtr
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:418
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
ADOLCPlasticityMaterialModels.hpp
Matetial models for plasticity.
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
ADOLCPlasticity.hpp
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
Monitor
[Push operators to pipeline]
Definition: adolc_plasticity.cpp:333
ADOLCPlasticity::ClosestPointProjection::a_y
adouble a_y
Definition: ADOLCPlasticity.hpp:329
MoFEM::VecOfTimeScalingMethods
std::vector< boost::shared_ptr< ScalingMethod > > VecOfTimeScalingMethods
Vector of time scaling methods.
Definition: Natural.hpp:20
MoFEM::DMMoFEMTSSetMonitor
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1060
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
ADOLCPlasticity::ParaboloidalPlasticity::LAMBDA
@ LAMBDA
Definition: ADOLCPlasticityMaterialModels.hpp:558
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
ADOLCPlasticity::voight_to_stress_op
FTensor::Dg< double, 3, 6 > voight_to_stress_op()
Op convert Vight stress vector to stress tensor.
Definition: ADOLCPlasticity.hpp:49
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:198
ADOLCPlasticity::ClosestPointProjection::a_sTress
ublas::vector< adouble > a_sTress
Definition: ADOLCPlasticity.hpp:327
ADOLCPlasticity::ParaboloidalPlasticity::ParamsArray
std::array< double, LAST_PARAM > ParamsArray
Definition: ADOLCPlasticityMaterialModels.hpp:570
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
ADOLCPlasticity::ParaboloidalPlasticity::MU
@ MU
Definition: ADOLCPlasticityMaterialModels.hpp:559
PlasticProblem::cpPtr
boost::shared_ptr< ClosestPointProjection > cpPtr
Definition: adolc_plasticity.cpp:125
ADOLCPlasticity::ParaboloidalPlasticity::sigmaYc
double sigmaYc
Definition: ADOLCPlasticityMaterialModels.hpp:554
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
ContactNaturalBC.hpp
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
ADOLCPlasticity::ParaboloidalPlasticity::setParams
MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent)
Set parameters for ADOL-C of constitutive relations.
Definition: ADOLCPlasticityMaterialModels.hpp:647
ADOLCPlasticity::ParaboloidalPlasticity::yieldFunction
MoFEMErrorCode yieldFunction()
Set yield function.
Definition: ADOLCPlasticityMaterialModels.hpp:695
PlasticOps::trace
double trace(FTensor::Tensor2_symmetric< T, 2 > &t_stress)
Definition: PlasticOpsGeneric.hpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
AssemblyDomainEleOp
FormsIntegrators< DomainEleOp >::Assembly< A >::OpBase AssemblyDomainEleOp
Definition: tensor_divergence_operator.cpp:59
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
ADOLCPlasticity::ParaboloidalPlasticity::SIGMA_YT
@ SIGMA_YT
Definition: ADOLCPlasticityMaterialModels.hpp:561
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
third
constexpr double third
Definition: EshelbianADOL-C.cpp:14
ADOLCPlasticity::ParaboloidalPlasticity::NC
@ NC
Definition: ADOLCPlasticityMaterialModels.hpp:566
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::SmartPetscObj< DM >
ADOLCPlasticity::ParaboloidalPlasticity::HC
@ HC
Definition: ADOLCPlasticityMaterialModels.hpp:564
PlasticProblem::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: adolc_plasticity.cpp:130
MoFEM::CoreInterface::get_field_entities_by_type
virtual MoFEMErrorCode get_field_entities_by_type(const std::string name, EntityType type, Range &ents) const =0
get entities in the field by type
ADOLCPlasticity::ParaboloidalPlasticity::tStress
FTensor::Tensor2_symmetric< adouble, 3 > tStress
Definition: ADOLCPlasticityMaterialModels.hpp:545
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:590
ADOLCPlasticity::ParaboloidalPlasticity::freeHemholtzEnergy
MoFEMErrorCode freeHemholtzEnergy()
[Paraboloidal free energy]
Definition: ADOLCPlasticityMaterialModels.hpp:653
OP
ADOLCPlasticity::ParaboloidalPlasticity::lambda
double lambda
Definition: ADOLCPlasticityMaterialModels.hpp:550
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1060
H
double H
Hardening.
Definition: plastic.cpp:175
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1289
OpCalculateVectorFieldGradient
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
PlasticProblem::mField
MoFEM::Interface & mField
Definition: adolc_plasticity.cpp:112
MU
#define MU(E, NU)
Definition: fem_tools.h:23
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
MoFEM::TSAdaptCreateMoFEM
PetscErrorCode TSAdaptCreateMoFEM(TSAdapt adapt)
Craete MOFEM adapt.
Definition: TsCtx.cpp:797
MoFEM::OpDGProjectionMassMatrix
Definition: DGProjection.hpp:106
ADOLCPlasticity::ParaboloidalPlasticity::J2
adouble J2
Definition: ADOLCPlasticityMaterialModels.hpp:678
PlasticOps::addMatBlockOps
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip, boost::shared_ptr< MatrixDouble > mat_D_Ptr, boost::shared_ptr< CommonData::BlockParams > mat_params_ptr, double scale_value, Sev sev)
Definition: PlasticOps.hpp:248