This program was contributed by Maien Hamed <[email protected]>.
It comes without any warranty or support by its authors or the authors of deal.II.

This program is part of the deal.II code gallery and consists of the following files (click to inspect):

Annotated version of README.md

PlasticityLab: ALE finite-strain thermoplasticity (axisymmetric)

Overview

This code-gallery entry provides a deal.II-based implementation of a large-deformation thermomechanical solver for finite-strain thermoplasticity, developed to support numerical simulation of severe-deformation metal forming processes and related problems involving strong thermomechanical coupling and viscoplastic flow.

The solver implements a finite-strain associative coupled thermoplasticity model.

Key idea: ALE via incremental reference motion

The solver incorporates an Arbitrary Lagrangian–Eulerian (ALE) formulation for coupled finite-strain thermoplasticity in which the motion of the reference configuration is represented incrementally through a reference velocity field. This avoids the need to explicitly track the deformation from the initial material configuration, either as a deformation field or through storing and updating the full deformation gradient history.

The method targets regimes where large accumulated strains can cause excessive mesh distortion in purely Lagrangian finite element simulations. The ALE formulation reduces sensitivity to mesh distortion and enables stable simulation without requiring prohibitively small time steps.

Physics and models

Finite-strain thermoplasticity

The constitutive model is a finite-strain, associative, coupled thermoplasticity formulation based on a multiplicative decomposition of the deformation gradient and a J2 (von Mises) flow theory. The implementation supports viscoplastic behavior via rate-dependent flow stress models (including Johnson-Cook).

Thermomechanical coupling

The thermal problem is coupled to the mechanical response through plastic dissipation and heat conduction.

Discretization and solution strategy (high level)

Mixed finite element formulation to mitigate volumetric locking (Jacobian/pressure treated as additional unknowns).

Newton-Raphson nonlinear solve with consistent tangent moduli (targeting second-order convergence behavior).

Mechanical-thermal operator splitting per time step (mechanical, then thermal, then mechanical sub-step).

Axisymmetric reduction and benchmark problems

Although the formulation is derived for general 3D settings, the code uses an axisymmetric approximation for benchmark problems and representative manufacturing-process simulations.

The entry validates and illustrates the approach using benchmark problems including:

thermally triggered necking of a circular bar (thermoplasticity benchmark),

Taylor anvil impact of a circular bar (dynamic high-rate deformation benchmark),

To run

# in a build directory:
@f$ cmake -DDEAL_II_DIR=<path-to-deal-ii> <path-to-entry>
@f$ make release
@f$ make -j 8 &&  mpirun -n 18 ./PlasticityLab

Notes on configuration

Geometry / triangulation: configured in PlasticityLabProgDrivers.cpp in run().

Material model selection and parameters: configured in main.cpp.

Time step settings: configured in PlasticityLabProg.h.

References

@article{HamedMcBrideReddy2023_ALE_Thermoplasticity_FrictionWelding,
  author  = {Hamed, M. M. O. and McBride, A. T. and Reddy, B. D.},
  title   = {An {ALE} approach for large-deformation thermoplasticity with application to friction welding},
  journal = {Computational Mechanics},
  volume  = {72},
  pages   = {803--826},
  year    = {2023},
  doi     = {10.1007/s00466-023-02303-0}
}

Annotated version of src/BodyForceApplier.h

  /*
   * BodyForceApplier.h
   *
   *  Created on: 14 Jan 2015
   *      Author: maien
   */
  
  #ifndef BODYFORCEAPPLIER_H_
  #define BODYFORCEAPPLIER_H_
  
  namespace PlasticityLab {
  
    template <int dim, typename Number = double>
    class BodyForceApplier {
     public:
      BodyForceApplier();
      BodyForceApplier(int direction, Number bodyForceMagnitude = 0);
      virtual ~BodyForceApplier();
      inline Number apply(const unsigned int direction,
                          const Number shapeFunctionValue,
                          const Number JxW) const;
     private:
      const unsigned int direction;
      const Number bodyForceMagnitude;
    };
  
  
    template <int dim, typename Number>
    BodyForceApplier<dim, Number>::
    BodyForceApplier(int direction, Number bodyForceMagnitude)
      : direction(direction), bodyForceMagnitude(bodyForceMagnitude) {
    }
  
  
    template <int dim, typename Number>
    BodyForceApplier<dim, Number>::~BodyForceApplier() {
    }
  
    template <int dim, typename Number>
    Number BodyForceApplier<dim, Number>::
    apply(const unsigned int direction,
          const Number shapeFunctionValue,
          const Number JxW) const {
      if (this->direction == direction)
        return -shapeFunctionValue * this->bodyForceMagnitude * JxW;
      return 0.0;
    }
  
  } /* namespace PlasticityLab */
  
  #endif /* BODYFORCEAPPLIER_H_ */

Annotated version of src/BoundaryUnidirectionalPenaltySpec.h

  /*
   * BoundaryUnidirectionalPenaltySpec.h
   *
   *  Created on: 04 May 2021
   *      Author: maien
   */
  
  #ifndef BOUNDARYUNIDIRECTIONALPENALTYSPEC_H_
  #define BOUNDARYUNIDIRECTIONALPENALTYSPEC_H_
  
  namespace PlasticityLab {
  
    template<typename Number = double>
    class BoundaryUnidirectionalPenaltySpec {
    public:
      BoundaryUnidirectionalPenaltySpec(
        unsigned int boundary_id,
        Number reference_displacement_increment,
        Number residual_force,
        Number quadratic_spring_factor) :
          boundary_id(boundary_id),
          reference_displacement_increment(reference_displacement_increment),
          residual_force(residual_force),
          quadratic_spring_factor(quadratic_spring_factor) {}
  
      unsigned int get_boundary_id() const { return boundary_id; }
      Number get_reference_displacement_increment() const { return reference_displacement_increment; }
      Number get_residual_force() const { return residual_force; }
      Number get_quadratic_spring_factor() const { return quadratic_spring_factor; }
  
    private:
      const unsigned int boundary_id;
      const Number reference_displacement_increment;
      const Number residual_force;
      const Number quadratic_spring_factor;
    };
  
  } /* namespace PlasticityLab */
  
  #endif /* BOUNDARYUNIDIRECTIONALPENALTYSPEC_H_ */

Annotated version of src/Constants.h

  /*
   * Constants.h
   *
   *  Created on: 10 Feb 2015
   *      Author: maien
   */
  
  #ifndef CONSTANTS_H_
  #define CONSTANTS_H_
  
  namespace PlasticityLab {
  
  
    template <int dim, typename Number>
    class Constants {
     public:
      inline static const Number one_third() {
        return static_cast<Number>(0.33333333333333333333333333333333333333333333333333);
      }
  
      inline static const Number sqrt2thirds() {
        return static_cast<Number>(0.81649658092772603273242802490196379732198249355222);
      }
  
      inline static const Number two_thirds() {
        return static_cast<Number>(0.66666666666666666666666666666666666666666666666666);
      }
  
      inline static const Number sqrt_half() {
        return static_cast<Number>(0.70710678118654752440084436210484903928483593768847);
      }
  
      inline static const Number sqrt_2() {
        return static_cast<Number>(1.41421356237309504880168872420969807856967187537694);
      }
  
  
      inline static const Number one_over_dim() {
        if(3==dim)
          return static_cast<Number>(0.33333333333333333333333333333333333333333333333333);
        else if (2==dim)
          return static_cast<Number>(0.5);
        else
          return static_cast<Number>(1./static_cast<Number>(dim));
      }
  
      inline static const Number two_over_dim() {
        if(3==dim)
          return static_cast<Number>(0.66666666666666666666666666666666666666666666666666);
        else if (2==dim)
          return static_cast<Number>(1.0);
        else
          return static_cast<Number>(2./static_cast<Number>(dim));
      }
  
    };
  
    inline void get_generalized_alpha_method_params(
          double *alpha_m,
          double *alpha_f,
          double *gamma,
          double *beta,
          double rho_infty
        ) {
      *alpha_m = (2. * rho_infty - 1.)/(rho_infty + 1.);
      *alpha_f = rho_infty / (rho_infty + 1.);
      *gamma = 0.5 - *alpha_m + *alpha_f;
      *beta = 0.25 * (1. - *alpha_m + *alpha_f) * (1. - *alpha_m + *alpha_f);
    }
  
  } /* namespace PlasticityLab */
  
  #endif /* CONSTANTS_H_ */

Annotated version of src/ConstitModelUpdateFlags.h

  /*
   * ConstitModelUpdateFlags.h
   *
   *  Created on: 04 Jan 2015
   *      Author: maien
   */
  
  #ifndef CONSTITMODELUPDATEFLAGS_H_
  #define CONSTITMODELUPDATEFLAGS_H_
  
  namespace PlasticityLab {
  
    enum ConstitutiveModelUpdateFlags {
      update_default = 0x0000,
      update_pressure = 0x0001,
      update_pressure_tangent = 0x0002,
      update_stress_deviator = 0x0004,
      update_stress_deviator_tangent = 0x0008,
      update_heat_flux = 0x0010,
      update_heat_flux_tangent = 0x0020,
      update_elastic_entropy = 0x0040,
      update_elastic_entropy_tangent = 0x0080,
      update_mechanical_dissipation = 0x0100,
      update_mechanical_dissipation_tangent = 0x0200,
      update_thermoelastic_heating = 0x0400,
      update_thermoelastic_heating_tangent = 0x0800,
      update_stored_heat = 0x1000,
      update_stored_heat_tangent = 0x2000,
      update_material_point_history = 0x4000
    };
  
    inline
    ConstitutiveModelUpdateFlags
    operator | (ConstitutiveModelUpdateFlags f1, ConstitutiveModelUpdateFlags f2) {
      return static_cast<ConstitutiveModelUpdateFlags> (
               static_cast<unsigned int> (f1) |
               static_cast<unsigned int> (f2));
    }
  
    inline
    const ConstitutiveModelUpdateFlags &
    operator |= (ConstitutiveModelUpdateFlags &f1, ConstitutiveModelUpdateFlags f2) {
      f1 = f1 | f2;
      return f1;
    }
  
    inline
    ConstitutiveModelUpdateFlags
    operator & (ConstitutiveModelUpdateFlags f1, ConstitutiveModelUpdateFlags f2) {
      return static_cast<ConstitutiveModelUpdateFlags> (
               static_cast<unsigned int> (f1) &
               static_cast<unsigned int> (f2));
    }
  
    inline
    const ConstitutiveModelUpdateFlags &
    operator &= (ConstitutiveModelUpdateFlags &f1, ConstitutiveModelUpdateFlags f2) {
      f1 = f1 & f2;
      return f1;
    }
  
  } /* namespace PlasticityLab */
  
  #endif /*CONSTITMODELUPDATEFLAGS_H_*/

Annotated version of src/ConstitutiveModelRequest.h

  /*
   * ConstitutiveModelRequest.h
   *
   *  Created on: 03 Feb 2015
   *      Author: maien
   */
  
  #ifndef CONSTITUTIVEMODELREQUEST_H_
  #define CONSTITUTIVEMODELREQUEST_H_
  
  #include <deal.II/base/tensor.h>
  #include <deal.II/base/symmetric_tensor.h>
  
  #include "Constants.h"
  #include "ConstitModelUpdateFlags.h"
  #include "TensorUtilities.h"
  
  namespace PlasticityLab {
  
    template <int dim, typename Number>
    class ConstitutiveModelRequest {
     public:
      ConstitutiveModelRequest(ConstitutiveModelUpdateFlags);
      virtual ~ConstitutiveModelRequest();
  

Interface to be used by request client (FE system assembler) –request configuration stage–

      void set_deformation_gradient(const Tensor<2, dim, Number> &deformation_gradient);
      void set_deformation_Jacobian(const Number deformation_Jacobian);
      void set_unprojected_deformation_Jacobian(const Number unprojected_deformation_Jacobian);
      void set_previous_deformation_Jacobian(const Number previous_deformation_Jacobian);
      void set_deformation_Jacobian_time_rate(const Number deformation_Jacobian_time_rate);
      void set_temperature(const Number temperature);
      void set_previous_temperature(const Number previous_temperature);
      void set_temperature_time_rate(const Number temperature_time_rate);
      void set_thermal_gradient(const Tensor<1, dim, Number> &thermalGradient);
      void set_time_increment(const Number timeIncrement);
  

Interface to be used by request client (FE system assembler) –request response retrieval and interrogation stage–

      Number get_pressure();
      Number get_pressure_tangent(const Number volume_change_increment);
      SymmetricTensor<2, dim, Number> get_stress_deviator() const;
      SymmetricTensor<2, dim, Number> get_stress_deviator_tangent(const Tensor<2, dim, Number> &strain_increment) const;
      Tensor<1, dim, Number> get_heat_flux() const;
      Tensor<1, dim, Number> get_heat_flux_tangent(const Tensor<1, dim, Number> &thermal_gradient_increment) const;
      Number  get_stored_heat_rate() const;
      Number  get_stored_heat_rate_tangent(const Number temperature_increment) const;
      Number  get_elastic_entropy() const;
      bool    get_is_plastic() const;
      Number  get_elastic_entropy_tangent(const Number temperature_increment) const;
      Number  get_mechanical_dissipation() const;
      Number  get_mechanical_dissipation_tangent(const Number temperature_increment) const;
      Number  get_thermo_elastic_heating() const;
      Number  get_thermo_elastic_heating_tangent(const Number temperature_increment) const;
  

interface used by constitutive model object to perform computation TODO consider hiding this interface and exposing it through adapter

      ConstitutiveModelUpdateFlags get_update_flags() const;
      Tensor<2, dim, Number> get_deformation_gradient() const;
      Number get_deformation_Jacobian() const;
      Number get_unprojected_deformation_Jacobian() const;
      Number get_previous_deformation_Jacobian() const;
      Number get_deformation_Jacobian_time_rate() const;
      Number get_temperature() const;
      Number get_previous_temperature() const;
      Number get_temperature_time_rate() const;
      Tensor<1, dim, Number> get_thermal_gradient() const;
      Number get_time_increment() const;
  
      void set_pressure(Number pressure);
      void set_stress_deviator(const SymmetricTensor<2, dim, Number> &stress_deviator);
      void set_heat_flux(const Tensor<1, dim, Number> &heat_flux);
      void set_stored_heat_rate(const Number stored_heat_rate);
      void set_elastic_entropy(const Number elastic_entropy);
      void set_mechanical_dissipation(const Number mechanical_dissipation);
      void set_thermo_elastic_heating(const Number thermo_elastic_heating);
  

TODO this can be changed so that smaller objects can be set and used to construct the tangents than the full moduli tensors

      void set_pressure_tangent_modulus(const Number pressure_tangent_modulus);
      void set_b_e_bar(const SymmetricTensor<2, dim, Number> &b_e_bar);
      void set_mu(const Number mu);
      void set_is_plastic(const bool is_plastic);
      void set_delta_gamma(const Number delta_gamma);
      void set_dK(const Number dK);
      void set_dH(const Number dH);
      void set_heat_flux_tangent_moduli(const SymmetricTensor<2, dim, Number> &heat_flux_tangent_modului);
      void set_stored_heat_rate_tangent_modulus(const Number stored_heat_rate_tangent_modulus);
      void set_elastic_entropy_tangent_modulus(const Number elastic_entropy_tangent_modulus);
      void set_mechanical_dissipation_tangent_modulus(const Number mechanicalDissipationTangentModulus);
      void set_thermo_elastic_heating_tangent_modulus(const Number thermo_elastic_heating_tangent_modulus);
  
     protected:
      ConstitutiveModelUpdateFlags update_flags;
  
      Tensor<2, dim, Number> deformation_gradient;
      Number deformation_Jacobian, previous_deformation_Jacobian, deformation_Jacobian_time_rate;
      Number unprojected_deformation_Jacobian;
      Number temperature, previous_temperature, temperature_time_rate;
      Tensor<1, dim, Number> thermal_gradient;
  
      Number pressure;
      SymmetricTensor<2, dim, Number> stress_deviator;
      Tensor<1, dim, Number> heat_flux;
      Number stored_heat_rate;
      Number  elastic_entropy;
      Number  mechanical_dissipation;
      Number  thermo_elastic_heating;
  
      Number time_increment;
  
      bool is_plastic;
      Number pressure_tangent_modulus;
      Number dK, dH, mu, delta_gamma;
      SymmetricTensor<2, dim, Number> b_e_bar;
      SymmetricTensor<2, dim, Number> heat_flux_tangent_moduli;
      Number  stored_heat_rate_tangent_modulus;
      Number  elastic_entropy_tangent_modulus;
      Number  mechanical_dissipation_tangent_modulus;
      Number  thermo_elastic_heating_tangent_modulus;
    };
  
    template <int dim, typename Number>
    ConstitutiveModelRequest<dim, Number>::
    ConstitutiveModelRequest(ConstitutiveModelUpdateFlags update_flags):
      update_flags(update_flags) {
        is_plastic = true;  // not necessarily elastic
    }
  
    template <int dim, typename Number>
    bool ConstitutiveModelRequest<dim, Number>::get_is_plastic() const {
      return is_plastic;
    }
  
    template <int dim, typename Number>
    ConstitutiveModelRequest<dim, Number>::~ConstitutiveModelRequest() { }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_deformation_gradient(const Tensor<2, dim, Number> &deformation_gradient) {
      this->deformation_gradient = deformation_gradient;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_deformation_Jacobian(Number deformation_Jacobian) {
      this->deformation_Jacobian = deformation_Jacobian;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_unprojected_deformation_Jacobian(Number unprojected_deformation_Jacobian) {
      this->unprojected_deformation_Jacobian = unprojected_deformation_Jacobian;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_previous_deformation_Jacobian(Number previous_deformation_Jacobian) {
      this->previous_deformation_Jacobian = previous_deformation_Jacobian;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_deformation_Jacobian_time_rate(Number deformation_Jacobian_time_rate) {
      this->deformation_Jacobian_time_rate = deformation_Jacobian_time_rate;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_temperature(const Number temperature) {
      this->temperature = temperature;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_previous_temperature(const Number previous_temperature) {
      this->previous_temperature = previous_temperature;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_temperature_time_rate(const Number temperature_time_rate) {
      this->temperature_time_rate = temperature_time_rate;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_thermal_gradient(const Tensor<1, dim, Number> &thermal_gradient) {
      this->thermal_gradient = thermal_gradient;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_time_increment(const Number time_increment) {
      this->time_increment = time_increment;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::get_pressure() {
      return pressure;
    }
  
    template <int dim, typename Number>
    Number
    ConstitutiveModelRequest<dim, Number>::
    get_pressure_tangent(const Number volume_change_increment) {
      return pressure_tangent_modulus * volume_change_increment;
    }
  
    template <int dim, typename Number>
    SymmetricTensor<2, dim, Number> ConstitutiveModelRequest<dim, Number>::
    get_stress_deviator() const {
      return stress_deviator;
    }
  
    template <int dim, typename Number>
    SymmetricTensor<2, dim, Number> ConstitutiveModelRequest<dim, Number>::
    get_stress_deviator_tangent(const Tensor<2, dim, Number> &strain_increment) const {
      const Number twothirds = Constants<dim, Number>::two_thirds();
      const auto tensor_b_e_bar = static_cast<Tensor<2,dim,Number> >(b_e_bar);
      SymmetricTensor<2, dim, Number> d_b_e_bar = symmetrize(2 * strain_increment * tensor_b_e_bar);
      SymmetricTensor<2, dim, Number> d_dev_b_e_bar = get_log_of_tensor_variation(b_e_bar, d_b_e_bar);
      SymmetricTensor<2, dim, Number> d_trial_stress_dev = mu * d_dev_b_e_bar;
  

TODO ensure that all the debugging tests were removed

if (is_plastic) {

Number mu_bar = Constants<dim, Number>::one_third() * mu * trace(b_e_bar); Number d_mu_bar = Constants<dim, Number>::one_third() * mu * trace(d_b_e_bar); Number norm_dev_b_e_bar = (deviator(b_e_bar)).norm(); SymmetricTensor<2, dim, Number> dev_b_e_direction = deviator(b_e_bar) / norm_dev_b_e_bar;

        Number mu_bar = mu;
        Number d_mu_bar = 0;
  
        const auto epsilon_e_bar = get_log_of_tensor(b_e_bar);
        Number norm_dev_b_e_bar = (epsilon_e_bar).norm();
        SymmetricTensor<2, dim, Number> dev_b_e_direction = epsilon_e_bar / norm_dev_b_e_bar;
  
        SymmetricTensor<2, dim, Number> d_dev_b_e_direction =
          (1.0 / norm_dev_b_e_bar) * (d_dev_b_e_bar - dev_b_e_direction * (dev_b_e_direction * d_dev_b_e_bar));
        Number d_delta_gamma = (dev_b_e_direction * d_trial_stress_dev - 2 * d_mu_bar * delta_gamma) / (2 * mu_bar + twothirds * (dK + dH));
  
        return deviator(d_trial_stress_dev
                        - (  2 * mu_bar * delta_gamma * d_dev_b_e_direction
                           + 2 * mu_bar * d_delta_gamma * dev_b_e_direction
                           + 2 * d_mu_bar * delta_gamma * dev_b_e_direction));
      }
  
      return deviator(d_trial_stress_dev);
    }
  
    template <int dim, typename Number>
    Tensor<1, dim, Number>
    ConstitutiveModelRequest<dim, Number>::get_heat_flux() const {
      return heat_flux;
    }
  
    template <int dim, typename Number>
    Tensor<1, dim, Number>
    ConstitutiveModelRequest<dim, Number>::
    get_heat_flux_tangent(const Tensor<1, dim, Number> &thermal_gradient_increment) const {
      return heat_flux_tangent_moduli * thermal_gradient_increment;
    }
  
    template <int dim, typename Number>
    Number
    ConstitutiveModelRequest<dim, Number>::get_stored_heat_rate() const {
      return stored_heat_rate;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_stored_heat_rate_tangent(const Number temperature_increment) const {
      return stored_heat_rate_tangent_modulus * temperature_increment;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_elastic_entropy() const {
      return elastic_entropy;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_elastic_entropy_tangent(const Number temperature_increment) const {
      return elastic_entropy_tangent_modulus * time_increment;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_mechanical_dissipation() const {
      return mechanical_dissipation;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_mechanical_dissipation_tangent(const Number temperature_increment) const {
      return mechanical_dissipation_tangent_modulus * temperature_increment;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_thermo_elastic_heating() const {
      return thermo_elastic_heating;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_thermo_elastic_heating_tangent(const Number temperature_increment) const {
      return thermo_elastic_heating_tangent_modulus * temperature_increment;
    }
  
    template <int dim, typename Number>
    ConstitutiveModelUpdateFlags ConstitutiveModelRequest<dim, Number>::
    get_update_flags() const {
      return update_flags;
    }
  
    template <int dim, typename Number>
    Tensor<2, dim, Number> ConstitutiveModelRequest<dim, Number>::
    get_deformation_gradient() const {
      return deformation_gradient;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_deformation_Jacobian() const {
      return deformation_Jacobian;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_unprojected_deformation_Jacobian() const {
      return unprojected_deformation_Jacobian;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_previous_deformation_Jacobian() const {
      return previous_deformation_Jacobian;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_deformation_Jacobian_time_rate() const {
      return deformation_Jacobian_time_rate;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_temperature() const {
      return temperature;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_previous_temperature() const {
      return previous_temperature;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_temperature_time_rate() const {
      return temperature_time_rate;
    }
  
    template <int dim, typename Number>
    Tensor<1, dim, Number> ConstitutiveModelRequest<dim, Number>::
    get_thermal_gradient() const {
      return thermal_gradient;
    }
  
    template <int dim, typename Number>
    Number ConstitutiveModelRequest<dim, Number>::
    get_time_increment() const {
      return time_increment;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_pressure(Number pressure) {
      this->pressure = pressure;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_stress_deviator(const SymmetricTensor<2, dim, Number> &stress_deviator) {
      this->stress_deviator = stress_deviator;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_b_e_bar(const SymmetricTensor<2, dim, Number> &b_e_bar) {
      this->b_e_bar = b_e_bar;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_mu(const Number mu) {
      this->mu = mu;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_is_plastic(const bool is_plastic) {
      this->is_plastic = is_plastic;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_delta_gamma(const Number delta_gamma) {
      this->delta_gamma = delta_gamma;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_dK(const Number dK) {
      this->dK = dK;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_dH(const Number dH) {
      this->dH = dH;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_heat_flux(const Tensor<1, dim, Number> &heat_flux) {
      this->heat_flux = heat_flux;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_stored_heat_rate(const Number stored_heat_rate) {
      this->stored_heat_rate = stored_heat_rate;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_elastic_entropy(const Number elastic_entropy) {
      this->elastic_entropy = elastic_entropy;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_mechanical_dissipation(const Number mechanical_dissipation) {
      this->mechanical_dissipation = mechanical_dissipation;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_thermo_elastic_heating(const Number thermo_elastic_heating) {
      this->thermo_elastic_heating = thermo_elastic_heating;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_pressure_tangent_modulus(const Number pressure_tangent_modulus) {
      this->pressure_tangent_modulus = pressure_tangent_modulus;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_heat_flux_tangent_moduli(const SymmetricTensor<2, dim, Number> &heat_flux_tangent_modului) {
      this->heat_flux_tangent_moduli = heat_flux_tangent_modului;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_stored_heat_rate_tangent_modulus(const Number stored_heat_rate_tangent_modulus) {
      this->stored_heat_rate_tangent_modulus = stored_heat_rate_tangent_modulus;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_elastic_entropy_tangent_modulus(const Number elastic_entropy_tangent_modulus) {
      this->elastic_entropy_tangent_modulus = elastic_entropy_tangent_modulus;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_mechanical_dissipation_tangent_modulus(const Number mechanicalDissipationTangentModulus) {
      this->mechanical_dissipation_tangent_modulus = mechanicalDissipationTangentModulus;
    }
  
    template <int dim, typename Number>
    void ConstitutiveModelRequest<dim, Number>::
    set_thermo_elastic_heating_tangent_modulus(const Number thermo_elastic_heating_tangent_modulus) {
      this->thermo_elastic_heating_tangent_modulus = thermo_elastic_heating_tangent_modulus;
    }
  
  } /* namespace PlasticityLab */
  
  #endif /* CONSTITUTIVEMODELREQUEST_H_ */

Annotated version of src/ConvectionBoundaryConditionApplier.h

  /*
   * ConvectionBoundaryConditionApplier.h
   *
   *  Created on: 07 Oct 2017
   *      Author: maien
   */
  
  #ifndef CONVECTIONBOUNDARYCONDITIONAPPLIER_H_
  #define CONVECTIONBOUNDARYCONDITIONAPPLIER_H_
  
  namespace PlasticityLab {
  
    template <int dim, typename Number = double>
    class ConvectionBoundaryConditionApplier {
     public:
      ConvectionBoundaryConditionApplier();
      ConvectionBoundaryConditionApplier(
        int direction,
        Number convection_coefficient = 1.0,
        Number ambient_field_value = 0.0);
      virtual ~ConvectionBoundaryConditionApplier();
      inline Number apply(const unsigned int direction,
                          const Number test_function_value,
                          const Number field_value,
                          const Number JxW) const;
      inline Number apply_gradient(
        const unsigned int direction,
        const Number &test_gradient,
        const Number &field_gradient,
        const Number JxW) const;
     private:
      const unsigned int direction;
      const Number convection_coefficient;
      const Number ambient_field_value;
    };
  
  
    template <int dim, typename Number>
    ConvectionBoundaryConditionApplier<dim, Number>::
    ConvectionBoundaryConditionApplier(
      int direction,
      Number convection_coefficient,
      Number ambient_field_value)
      : direction(direction),
        convection_coefficient(convection_coefficient),
        ambient_field_value(ambient_field_value) {
    }
  
  
    template <int dim, typename Number>
    ConvectionBoundaryConditionApplier<dim, Number>::~ConvectionBoundaryConditionApplier() {
    }
  
    template <int dim, typename Number>
    Number ConvectionBoundaryConditionApplier<dim, Number>::
    apply(const unsigned int direction,
          const Number test_function_value,
          const Number field_value,
          const Number JxW) const {
      if (this->direction == direction)
        return test_function_value * convection_coefficient * (field_value - ambient_field_value) * JxW;
      return 0.0;
    }
  
    template <int dim, typename Number>
    Number ConvectionBoundaryConditionApplier<dim, Number>::apply_gradient(
      const unsigned int direction,
      const Number &test_gradient,
      const Number &field_gradient,
      const Number JxW) const {
      if (this->direction == direction) {
        return convection_coefficient * (test_gradient * field_gradient) * JxW;
      }
      return 0;
    }
  
  
  } /* namespace PlasticityLab */
  
  #endif /* CONVECTIONBOUNDARYCONDITIONAPPLIER_H_ */

Annotated version of src/DoFSystem.h

  /*
   * DoFSystem.h
   *
   *  Created on: 05 May 2015
   *      Author: maien
   */
  
  #ifndef DOFSYSTEM_H_
  #define DOFSYSTEM_H_
  
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/base/conditional_ostream.h>
  
  #include "InterpolatoryConstraintApplier.h"
  #include "BodyForceApplier.h"
  #include "ConvectionBoundaryConditionApplier.h"
  #include "mpi.h"
  #include "utilities.h"
  
  using namespace dealii;
  
  namespace PlasticityLab {
  
    template <int dim, typename Number=double>
    class DoFSystem {
     public:
      DoFSystem (const ::Triangulation<dim> &triangulation,
                 const ::Mapping<dim> &mapping);
  
      DoFHandler<dim>     dof_handler;
  
      AffineConstraints<Number>    nodal_constraints;
  
      IndexSet           locally_owned_dofs;
      IndexSet           locally_relevant_dofs;
  
      void setup_dof_system (const FiniteElement<dim> &fe);
      const ::Mapping<dim> &mapping;
  
    };
  
    template<int dim, typename Number>
    DoFSystem <dim, Number> :: DoFSystem(const ::Triangulation<dim> &triangulation,
                                         const ::Mapping<dim> &mapping) :
      dof_handler(triangulation),
      mapping(mapping) {
    }
  
    template <int dim, typename Number>
    void DoFSystem<dim, Number>::setup_dof_system (const FiniteElement<dim> &fe) {
      dof_handler.distribute_dofs(fe);
      locally_owned_dofs = dof_handler.locally_owned_dofs();
      locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler);
  
      nodal_constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
      DoFTools::make_hanging_node_constraints (dof_handler, nodal_constraints);
    }
  
  } /* namespace PlasticityLab */
  
  #endif /* DOFSYSTEM_H_ */

Annotated version of src/ExponentialHardeningElastoplasticMaterial.cpp

  /*
   * ExponentialHardeningElastoplasticMaterial.cpp
   *
   *  Created on: 10 Jul 2014
   *      Author: cerecam
   */
  
  
  #include <math.h>
  #include <deal.II/base/symmetric_tensor.h>
  
  #include "ExponentialHardeningElastoplasticMaterial.h"
  #include "utilities.h"
  
  namespace PlasticityLab {
  
    template <int dim, typename Number>
    ExponentialHardeningElastoplasticMaterial<dim, Number>::
    ExponentialHardeningElastoplasticMaterial
    (const Number kappa,
     const Number mu,
     const Number K_0,
     const Number K_infty,
     const Number delta,
     const Number H_bar,
     const Number beta)    :
      kappa (kappa),
      mu (mu),
      K_0(K_0),
      K_infty(K_infty),
      delta(delta),
      H_bar(H_bar),
      beta(beta),
  
      stress_strain_tensor_kappa (kappa
                                  * outer_product(unit_symmetric_tensor<dim>(),
                                                  unit_symmetric_tensor<dim>())),
      stress_strain_tensor_mu (2 * mu
                               * (identity_tensor<dim>()
                                  - outer_product(unit_symmetric_tensor<dim>(),
                                                  unit_symmetric_tensor<dim>()) / 3.0)) {
  
    }
  
    template <int dim, typename Number>
    ExponentialHardeningElastoplasticMaterial<dim, Number>::
    ~ExponentialHardeningElastoplasticMaterial() {
    }
  
    template <int dim, typename Number>
    std::vector<Number> ExponentialHardeningElastoplasticMaterial<dim, Number>::get_state_parameters(
        const point_index_t &,
        const Tensor<2, dim, Number> &) const {
      throw NotImplementedException();
    }
  
    template <int dim, typename Number>
    void ExponentialHardeningElastoplasticMaterial<dim, Number>::
    set_state_parameters(
          const point_index_t &,
          const std::vector<Number> &,
          const Tensor<2, dim, Number> &) {
      throw NotImplementedException();
    }
  
    template <int dim, typename Number>
    size_t ExponentialHardeningElastoplasticMaterial<dim, Number>::
    get_material_parameter_count() const {
      throw NotImplementedException();
    }
  
  
    template <int dim, typename Number>
    Number ExponentialHardeningElastoplasticMaterial<dim, Number>::get_material_Jacobian(const point_index_t &) const {
      throw NotImplementedException();
    }
  
    template <int dim, typename Number>
    ::SymmetricTensor<2, dim, Number> ExponentialHardeningElastoplasticMaterial<dim, Number>::get_plastic_strain(const point_index_t &) const {
      throw NotImplementedException();
    }
  
  
    template <int dim, typename Number>
    void ExponentialHardeningElastoplasticMaterial<dim, Number>::
    compute_constitutive_request(ConstitutiveModelRequest<dim, Number> &constitutive_request,
                                 const point_index_t &point_index) {
      SymmetricTensor<2, dim, Number> plastic_strain = material_point_history[point_index].plastic_strain;
      typename PointHistory<dim, Number>::HardeningParameters
      hardening_parameters = material_point_history[point_index].hardening_parameters;
  
      auto deformation_gradient = static_cast<SymmetricTensor<2, dim, Number> >(constitutive_request.get_deformation_gradient());
      SymmetricTensor<4, dim, Number> elastoplastic_tangent_moduli;
      SymmetricTensor<2, dim, Number> deviator_strain_tensor = deviator(deformation_gradient);
  
      SymmetricTensor<2, dim, Number> dev_stress_trial = 2 * mu * (deviator_strain_tensor - plastic_strain);
      SymmetricTensor<2, dim, Number>
      ksi_trial = dev_stress_trial - hardening_parameters.kinematic_hardening;
  
      Number norm_ksi_trial = ksi_trial.norm();
      if (trial_yield_criterion( norm_ksi_trial, hardening_parameters.equivalent_plastic_strain ) > 0) {
        Number delta_gamma, alpha_n_plus_1;
        determine_delta_gamma(delta_gamma, alpha_n_plus_1,
                              norm_ksi_trial,
                              hardening_parameters.equivalent_plastic_strain,
                              1e-10, 200);
        SymmetricTensor<2, dim, Number> stress_flow_direction = ksi_trial / norm_ksi_trial;
  

Update back stress, plastic strain and stress

        const Number sqrt2thirds = sqrt((Number)2 / (Number)3);
        Number H_alpha_n_plus_1, H_alpha_n, K_alpha_n_plus_1, K_alpha_n;
        Number DH_alpha_n_plus_1, DK_alpha_n_plus_1;
        exponential_hardening_values(K_alpha_n, H_alpha_n, hardening_parameters.equivalent_plastic_strain);
        exponential_hardening_values(K_alpha_n_plus_1, H_alpha_n_plus_1, alpha_n_plus_1);
        exponential_hardening_derivatives(DK_alpha_n_plus_1, DH_alpha_n_plus_1, alpha_n_plus_1);
        if (update_material_point_history & constitutive_request.get_update_flags()) {
          material_point_history[point_index].hardening_parameters.equivalent_plastic_strain = alpha_n_plus_1;
          material_point_history[point_index].hardening_parameters.kinematic_hardening =
            hardening_parameters.kinematic_hardening
            + sqrt2thirds
            * (H_alpha_n_plus_1 - H_alpha_n)
            * stress_flow_direction;
          material_point_history[point_index].plastic_strain =
            plastic_strain + delta_gamma * stress_flow_direction;
        }
  
        SymmetricTensor<2, dim, Number>
        stress = kappa * trace(deformation_gradient) * unit_symmetric_tensor<dim, Number>()
                 + dev_stress_trial
                 - 2 * mu * delta_gamma * stress_flow_direction;
  
        Number theta_n_plus_1 = 1 - 2 * mu * delta_gamma / norm_ksi_trial;
        Number theta_bar_n_plus_1 = 1 / (1 + (DK_alpha_n_plus_1 + DH_alpha_n_plus_1) / (3 * mu))
                                    - (1 - theta_n_plus_1);
        const SymmetricTensor<4, dim, Number> one_prod_one =
          outer_product(unit_symmetric_tensor<dim, Number>(), unit_symmetric_tensor<dim, Number>());
        elastoplastic_tangent_moduli = kappa * one_prod_one
                                       + 2 * mu * theta_n_plus_1 * (identity_tensor<dim, Number>() - 1 / 3 * one_prod_one)
                                       - 2 * mu * theta_bar_n_plus_1 * outer_product(stress_flow_direction, stress_flow_direction);
  

TODO change code such that request update flags are respected

constitutive_request.set_stress_deviator(stress);

constitutiveRequest.setStressDeviatorTangentModuli(elastoplastic_tangent_moduli);

      } /*if ( trial yield criterion test )*/
      else {
        elastoplastic_tangent_moduli = stress_strain_tensor_kappa + stress_strain_tensor_mu;
        SymmetricTensor<2, dim, Number> stress = (stress_strain_tensor_kappa + stress_strain_tensor_mu) * deformation_gradient;
        constitutive_request.set_stress_deviator(stress);

constitutiveRequest.setStressDeviatorTangentModuli(elastoplastic_tangent_moduli);

      }
    }
  
    template <int dim, typename Number>
    void
    ExponentialHardeningElastoplasticMaterial<dim, Number>::
    setup_point_history (const point_index_t point_count) {
      {
        std::vector< PointHistory<dim, Number> > tmp;
        tmp.swap (material_point_history);
      }
      material_point_history.resize (point_count);
    }
  
    template <int dim, typename Number>
    inline void
    ExponentialHardeningElastoplasticMaterial<dim, Number>::
    determine_delta_gamma(Number &delta_gamma, Number &alpha_n_plus_1,
                          const Number norm_ksi_trial,
                          const Number alpha_n,
                          Number tol, unsigned int max_iter) const {
      unsigned int k = 0;
      const Number sqrt2thirds = sqrt((Number)2 / (Number)3);
      Number g_of_gamma_k, Dg_of_gamma_k;
      Number K_alpha_n, K_alpha_n_plus_1, H_alpha_n, H_alpha_n_plus_1;
      Number DH_alpha_n_plus_1, DK_alpha_n_plus_1;
  
      delta_gamma = 0;
      alpha_n_plus_1 = alpha_n;
  
      exponential_hardening_values(K_alpha_n, H_alpha_n, alpha_n);
  
      do {
        k++;
        exponential_hardening_values(K_alpha_n_plus_1, H_alpha_n_plus_1, alpha_n_plus_1);
        g_of_gamma_k = -sqrt2thirds * K_alpha_n_plus_1 + norm_ksi_trial
                       - (2 * mu * delta_gamma + sqrt2thirds * (H_alpha_n_plus_1 - H_alpha_n));
  
        exponential_hardening_derivatives(DK_alpha_n_plus_1, DH_alpha_n_plus_1, alpha_n_plus_1);
        Dg_of_gamma_k = -2 * mu * (1 + (DH_alpha_n_plus_1 + DK_alpha_n_plus_1) / (3 * mu));
  
        delta_gamma = delta_gamma - g_of_gamma_k / Dg_of_gamma_k;
        alpha_n_plus_1 = alpha_n + sqrt2thirds * delta_gamma;
  
      } while (std::fabs(g_of_gamma_k) > tol && k < max_iter);
    }
  
    template <int dim, typename Number>
    inline void
    ExponentialHardeningElastoplasticMaterial<dim, Number>::
    exponential_hardening_values(Number &kinematic_hardening,
                                 Number &isotropic_hardening,
                                 const Number alpha) const {
      Number h = K_infty - (K_infty - K_0) * exp(-delta * alpha) + H_bar * alpha;
      kinematic_hardening = beta * h;
      isotropic_hardening = (1 - beta) * h;
    }
  
    template <int dim, typename Number>
    inline void
    ExponentialHardeningElastoplasticMaterial<dim, Number>::
    exponential_hardening_derivatives(Number &D_kinematic_hardening,
                                      Number &D_isotropic_hardening,
                                      const Number alpha) const {
      Number Dh = delta * (K_infty - K_0) * exp(-delta * alpha) + H_bar;
      D_kinematic_hardening = beta * Dh;
      D_isotropic_hardening = (1 - beta) * Dh;
    }
  
    template <int dim, typename Number>
    inline Number
    ExponentialHardeningElastoplasticMaterial<dim, Number>::
    trial_yield_criterion(const Number norm_ksi_trial,
                          const Number alpha) const {
      Number H_alpha, K_alpha;
      exponential_hardening_values(K_alpha, H_alpha, alpha);
      const Number sqrt2thirds = sqrt((Number)2 / (Number)3);
      Number trial_yield_criterion = norm_ksi_trial - sqrt2thirds * K_alpha;
      return trial_yield_criterion;
    }
  
    template class ExponentialHardeningElastoplasticMaterial<3, double>;
    template class ExponentialHardeningElastoplasticMaterial<2, double>;
  
  } /* namespace PlasticityLab */

Annotated version of src/ExponentialHardeningElastoplasticMaterial.h

  /*
   * ExponentialHardeningElastoplasticMaterial.h
   *
   *  Created on: 10 Jul 2014
   *      Author: cerecam
   */
  
  #ifndef EXPONENTIALHARDENINGMATERIAL_H_
  #define EXPONENTIALHARDENINGMATERIAL_H_
  
  #include "PointHistory.h"
  #include "Material.h"
  #include "ConstitutiveModelRequest.h"
  
  using namespace dealii;
  
  namespace PlasticityLab {
  
    template <int dim, typename Number = double>
    class ExponentialHardeningElastoplasticMaterial : public Material<dim, Number> {
     public:
      ExponentialHardeningElastoplasticMaterial(const Number E,
                                                const Number nu,
                                                const Number K_0,
                                                const Number K_infty,
                                                const Number delta,
                                                const Number H_bar,
                                                const Number beta);
  
      virtual ~ExponentialHardeningElastoplasticMaterial();
  
      void compute_constitutive_request(
        ConstitutiveModelRequest<dim, Number> &constitutive_request,
        const point_index_t &point_index) override;
  
      Number get_material_Jacobian(const point_index_t &point_index) const override;
      ::SymmetricTensor<2, dim, Number> get_plastic_strain(const point_index_t &point_index) const override;
  
      void setup_point_history (const point_index_t point_count) override;
  
      std::vector<Number> get_state_parameters(
                const point_index_t &point_index,
                const Tensor<2, dim, Number> &reference_transformation=unit_symmetric_tensor<dim>()) const override;
  
      void  set_state_parameters(
                const point_index_t &point_index,
                const std::vector<Number> &state_parameters,
                const Tensor<2, dim, Number> &reference_transformation) override;
      size_t get_material_parameter_count() const override;
  
     private:
      const Number kappa;
      const Number mu;
  
      const Number K_0, K_infty, delta, H_bar; // hardening parameters (Simo & Hughes pp185)
      const Number beta;               // isotropic/kinematic hardening parameter
  
  
      const SymmetricTensor<4, dim, Number> stress_strain_tensor_kappa;
      const SymmetricTensor<4, dim, Number> stress_strain_tensor_mu;
  
      std::vector< PointHistory< dim, Number> > material_point_history;
  
      inline void
      determine_delta_gamma(Number &delta_gamma, Number &alpha_n_plus_1,
                            const Number norm_ksi_trial,
                            const Number alpha_n,
                            Number tol, unsigned int max_iter) const;
  
      inline void
      exponential_hardening_values(Number &kinematic_hardening,
                                   Number &isotropic_hardening,
                                   const Number alpha) const;
  
      inline void
      exponential_hardening_derivatives(Number &D_kinematic_hardening,
                                        Number &D_isotropic_hardening,
                                        const Number alpha) const;
  
      inline Number
      trial_yield_criterion(const Number norm_ksi_trial,
                            const Number alpha) const;
  
    };
  
  } /* namespace PlasticityLab */
  
  #endif /* EXPONENTIALHARDENINGMATERIAL_H_ */

Annotated version of src/ExponentialHardeningThermoviscoplasticYieldLaw.h

  /*
   * ExponentialHardeningThermoviscoplasticYieldLaw.h
   *
   *  Created on: 22 Nov 2019
   *      Author: maien
   */
  
  #ifndef EXPONENTIALHARDENINGTHERMOVISCOPLASTICYIELDLAW_H_
  #define EXPONENTIALHARDENINGTHERMOVISCOPLASTICYIELDLAW_H_
  
  #include "Constants.h"
  
  namespace PlasticityLab {
  
    template<typename Number>
    class ExponentialHardeningThermoviscoplasticYieldLaw {
    public:
      ExponentialHardeningThermoviscoplasticYieldLaw(
        const Number K_0,
        const Number K_infty,
        const Number delta,
        const Number H_bar,
        const Number beta,
        const Number flow_stress_softening,
        const Number hardening_softening,
        const Number reference_temperature=293.0) :
      K_0(K_0),
      K_infty(K_infty),
      delta(delta),
      H_bar(H_bar),
      beta(beta),
      flow_stress_softening(flow_stress_softening),
      hardening_softening(hardening_softening),
      reference_temperature(reference_temperature),
      viscous_hardening_factor(0.0),
      sqrt2thirds(Constants<3, Number>::sqrt2thirds()) {}
  
    Number hardening_values(Number &isotropic_hardening,
                                 Number &kinematic_hardening,
                                 const Number alpha,
                                 const Number gamma,
                                 const Number time_increment,
                                 const Number temperature) const {
      Number h = K_0 * (1 - std::min(softening_threshold, flow_stress_softening * (temperature - reference_temperature)))
                 + ((K_infty - K_0) * (1 - exp(-delta * alpha)) + H_bar * alpha) * (1 - std::min(softening_threshold, hardening_softening * (temperature - reference_temperature)))
                 + viscous_hardening_factor * sqrt2thirds * gamma / time_increment;
      isotropic_hardening = beta * h;
      kinematic_hardening = (1 - beta) * h;
      return h;
    }
  
    Number hardening_alpha_derivatives(Number &D_isotropic_hardening,
                                            Number &D_kinematic_hardening,
                                            const Number alpha,
                                            [[maybe_unused]] const Number gamma,
                                            const Number time_increment,
                                            const Number temperature) const {
      Number Dh = (delta * (K_infty - K_0) * exp(-delta * alpha) + H_bar) * (1 - std::min(softening_threshold, hardening_softening * (temperature - reference_temperature)))
                  + viscous_hardening_factor / time_increment;
      D_isotropic_hardening = beta * Dh;
      D_kinematic_hardening = (1 - beta) * Dh;
      return Dh;
    }
  
    Number hardening_temperature_derivatives(Number &D_isotropic_hardening,
                                                  Number &D_kinematic_hardening,
                                                  const Number alpha,
                                                  [[maybe_unused]] const Number gamma,
                                                  [[maybe_unused]] const Number time_increment,
                                                  const Number temperature) const {
      Number Dh = (hardening_softening * (temperature - reference_temperature) < softening_threshold)?
                  -flow_stress_softening * K_0
                  - hardening_softening * ((K_infty - K_0) * (1 - exp(-delta * alpha)) + H_bar * alpha)
                  : 0.0;
      D_isotropic_hardening = beta * Dh;
      D_kinematic_hardening = (1 - beta) * Dh;
      return Dh;
    }
  
    Number trial_yield_criterion(const Number norm_ksi_trial,
                          const Number alpha,
                          const Number gamma,
                          const Number time_increment,
                          const Number temperature) const {
      Number K_alpha, H_alpha;
      hardening_values(H_alpha, K_alpha, alpha, gamma, time_increment, temperature);
      const Number trial_yield_criterion = norm_ksi_trial - sqrt2thirds * (H_alpha);
      return trial_yield_criterion;
    }
    private:
      const Number K_0, K_infty, delta, H_bar; // hardening parameters (Simo & Hughes pp185)
      const Number beta; // isotropic/kinematic hardening parameter (1.0 for isotropic)
      const Number flow_stress_softening, hardening_softening; // thermal softening parameters (Simo & Miehe 1992 pp74)
      const Number reference_temperature;
      const Number viscous_hardening_factor;
      const Number sqrt2thirds;
      const Number softening_threshold = 0.98;
    };
  
  } /* namespace PlasticityLab */
  
  #endif /* EXPONENTIALHARDENINGTHERMOVISCOPLASTICYIELDLAW_H_ */ * 

Annotated version of src/IncrementInterpolationHandler.h

  /*
   * IncrementInterpolationHandler.h
   *
   *  Created on: 25 Oct 2019
   *      Author: maien
   */
  
  #ifndef INCREMENTINTERPOLATIONHANDLER_H_
  #define INCREMENTINTERPOLATIONHANDLER_H_
  
  
  using namespace dealii;
  
  namespace PlasticityLab {
  template <int dim, typename Number=double, int components=dim>
  class IncrementInterpolationHandler {
  public:
      IncrementInterpolationHandler(
              Function<dim, Number> *increment_interpolation_function,
              bool do_interpolate,
              ComponentMask interpolation_component_mask,
              bool do_constrain,
              ComponentMask constrain_component_mask,
              types::boundary_id constrain_boundary_id,
              Mapping<dim> &mapping)
      : increment_interpolation_function(increment_interpolation_function),
        do_interpolate(do_interpolate),
        interpolation_component_mask(interpolation_component_mask),
        do_constrain(do_constrain),
        constrain_component_mask(constrain_component_mask),
        constrain_boundary_id(constrain_boundary_id),
        mapping(mapping) { }
  
      ~IncrementInterpolationHandler(){
          delete increment_interpolation_function;
      }
  
      void advance_time(const Number delta_t);
  
      template<typename VectorType>
      void distribute_step_constraints(VectorType &increment) const {
          if(do_constrain) {
              function_constraint.distribute(increment);
          }
      }
  
      template<typename VectorType>
      void interpolate(VectorType &increment, const DoFSystem<dim, Number> &dof_system) const {
          if(do_interpolate) {
              VectorTools::interpolate(
                mapping,
                dof_system.dof_handler,
                *increment_interpolation_function,
                increment,
                interpolation_component_mask);
          }
      }
  
      void reinit_constraint_matrix(const DoFSystem<dim, Number> &dof_system) {
          if(do_constrain) {
              function_constraint.reinit(dof_system.locally_relevant_dofs);
              DoFTools::make_hanging_node_constraints(dof_system.dof_handler, function_constraint);
              std::map< types::boundary_id, const Function< dim, Number > * > constraint_function_map;
              constraint_function_map.insert(std::make_pair(constrain_boundary_id, increment_interpolation_function));
              ::VectorTools::interpolate_boundary_values(dof_system.mapping, dof_system.dof_handler, constraint_function_map, function_constraint, constrain_component_mask);
              function_constraint.close();
          }
      }
  
  private:
      Function<dim, Number>* increment_interpolation_function;
      const bool do_interpolate;
      ComponentMask interpolation_component_mask;
      const bool do_constrain;
      ComponentMask constrain_component_mask;
      types::boundary_id constrain_boundary_id;
      AffineConstraints<Number> function_constraint;
      const Mapping<dim> &mapping;
  };
  
  
  template<int dim, typename Number, int components>
  void IncrementInterpolationHandler<dim, Number, components>::advance_time(const Number delta_t) {
      increment_interpolation_function->advance_time(delta_t);
  }
  
  } /* namespace PlasticityLab */
  
  #endif /* INCREMENTINTERPOLATIONHANDLER_H_ */

Annotated version of src/InterpolatoryConstraintApplier.h

  /*
   * InterpolatoryConstraintApplier.h
   *
   *  Created on: 15 Jan 2015
   *      Author: maien
   */
  
  #ifndef INTERPOLATORYCONSTRAINTAPPLIER_H_
  #define INTERPOLATORYCONSTRAINTAPPLIER_H_
  
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/numerics/vector_tools.h>
  
  namespace PlasticityLab {
  
    template <int dim, typename Number = double>
    class InterpolatoryConstraintApplier {
     public:
      InterpolatoryConstraintApplier();
      InterpolatoryConstraintApplier(const std::map< ::types::boundary_id, const ::Function< dim, Number > * > &constraintFunctionMap,
                                     ::ComponentMask componentMask);
      virtual ~InterpolatoryConstraintApplier();
  
      void configure(std::map< ::types::boundary_id, const ::Function< dim, Number > * > constraintFunctionMap,
                     ::ComponentMask componentMask);
  
      void apply(const ::Mapping<dim> &mapping,
                 ::DoFHandler<dim> &doFHandler,
                 ::AffineConstraints<Number> &constraintMatrix,
                 bool useComponentMask = true) const;
  
     private:
      std::map< ::types::boundary_id, const ::Function< dim, Number > * > constraintFunctionMap;
      ::ComponentMask componentMask;
    };
  
    template <int dim, typename Number>
    InterpolatoryConstraintApplier<dim, Number>::InterpolatoryConstraintApplier() {
    }
  
    template <int dim, typename Number>
    InterpolatoryConstraintApplier<dim, Number>::InterpolatoryConstraintApplier
    (const std::map< ::types::boundary_id, const ::Function< dim, Number > * > &constraintFunctionMap,
     ::ComponentMask componentMask):
      constraintFunctionMap(constraintFunctionMap),
      componentMask(componentMask) {
    }
  
    template <int dim, typename Number>
    InterpolatoryConstraintApplier<dim, Number>::~InterpolatoryConstraintApplier() {
    }
  
    template <int dim, typename Number>
    void InterpolatoryConstraintApplier<dim, Number>::configure
    (std::map< ::types::boundary_id, const ::Function< dim, Number > * > constraintFunctionMap,
     ::ComponentMask componentMask) {
      this->constraintFunctionMap = std::map< ::types::boundary_id, const ::Function< dim, Number > * >(constraintFunctionMap);
      this->componentMask = ::ComponentMask(componentMask);
    }
  
    template <int dim, typename Number>
    void InterpolatoryConstraintApplier<dim, Number>::apply(const ::Mapping<dim> &mapping,
                                                            ::DoFHandler<dim> &doFHandler,
                                                            ::AffineConstraints<Number> &constraintMatrix,
                                                            bool useComponentMask) const {
      if (useComponentMask)
        ::VectorTools::interpolate_boundary_values(mapping,
                                                         doFHandler,
                                                         constraintFunctionMap,
                                                         constraintMatrix,
                                                         componentMask);
      else
        ::VectorTools::interpolate_boundary_values(mapping,
                                                         doFHandler,
                                                         constraintFunctionMap,
                                                         constraintMatrix);
    }
  
  } /* namespace PlasticityLab */
  
  #endif /* INTERPOLATORYCONSTRAINTAPPLIER_H_ */

Annotated version of src/JohnsonCookThermoviscoplasticYieldLaw.h

  /*
   * JohnsonCookThermoviscoplasticYieldLaw.h
   *
   *  Created on: 22 Nov 2019
   *      Author: maien
   */
  
  #ifndef JOHNSONCOOKTHERMOVISCOPLASTICYIELDLAW_H_
  #define JOHNSONCOOKTHERMOVISCOPLASTICYIELDLAW_H_
  
  #include <math.h>
  
  #include "Constants.h"
  
  
  namespace PlasticityLab {
  
    template<typename Number>
    class JohnsonCookThermoviscoplasticYieldLaw {
    public:
      JohnsonCookThermoviscoplasticYieldLaw(
        const Number mu,
        const Number A,
        const Number B,
        const Number C,
        const Number m,
        const Number n,
        const Number melting_temperature,
        const Number reference_strain_rate=1.0,
        const Number reference_temperature=293.0,
        const Number beta=1.0) :
      mu(mu),
      A(A),
      B(B),
      C(C),
      m(m),
      n(n),
      melting_temperature(melting_temperature),
      reference_strain_rate(reference_strain_rate),
      reference_temperature(reference_temperature),
      beta(beta),
      sqrt2thirds(Constants<2, Number>::sqrt2thirds()),
      exp_one_half(std::exp(0.5)) {}
  
    Number hardening_values(Number &isotropic_hardening,
                                 Number &kinematic_hardening,
                                 const Number alpha,
                                 const Number gamma,
                                 const Number time_increment,
                                 const Number temperature) const {
      if(use_Carreau_viscous_law) {
        const Number creep_strain_rate_factor = std::pow(1-std::max(0., std::min(1., (temperature - reference_temperature)/(melting_temperature - reference_temperature))), 1.5);
        const Number minimum_strain_rate = creep_strain_rate_factor * epsilon_dot_0;
        const Number strain_rate = std::max(sqrt2thirds*gamma/time_increment, minimum_strain_rate);
        const Number softened_quasistatic_elastoplastic_stress = get_elastoplastic_factor(alpha) * get_softening_factor(temperature);
        const Number Carreau_viscocity = get_Carreau_viscocity(strain_rate, softened_quasistatic_elastoplastic_stress);
        const Number h = 3 * Carreau_viscocity * strain_rate;
        isotropic_hardening = beta * h;
        kinematic_hardening = (1 - beta) * h;
        return h;
      } else {
        const Number h =
            get_elastoplastic_factor(alpha)
                * get_viscosity_factor(gamma, time_increment)
                * get_softening_factor(temperature)
            + viscosity_regularization_factor * sqrt2thirds * gamma/time_increment;
        isotropic_hardening = beta * h;
        kinematic_hardening = (1 - beta) * h;
        return h;
      }
    }
  
    Number hardening_alpha_derivatives(Number &D_isotropic_hardening,
                                            Number &D_kinematic_hardening,
                                            const Number alpha,
                                            const Number gamma,
                                            const Number time_increment,
                                            const Number temperature) const {
      if(use_Carreau_viscous_law) {
        const Number creep_strain_rate_factor = std::pow(1-std::max(0., std::min(1., (temperature - reference_temperature)/(melting_temperature - reference_temperature))), 1.5);
        const Number minimum_strain_rate = creep_strain_rate_factor * epsilon_dot_0;
        const Number strain_rate = std::max(sqrt2thirds*gamma/time_increment, minimum_strain_rate);
        const Number softened_quasistatic_elastoplastic_stress = get_elastoplastic_factor(alpha) * get_softening_factor(temperature);
        const Number Carreau_viscocity = get_Carreau_viscocity(strain_rate, softened_quasistatic_elastoplastic_stress);
        const Number strain_rate_tangent = strain_rate > minimum_strain_rate? 1.0/time_increment : 0;
        const Number stress_tangent =
          get_elastoplastic_factor_tangent(alpha) * get_softening_factor(temperature);
        Number Carreau_viscocity_strain_rate_tangent, Carreau_viscocity_stress_tangent;
        get_Carreau_viscocity_tangents(
            strain_rate, softened_quasistatic_elastoplastic_stress,
            Carreau_viscocity_strain_rate_tangent,
            Carreau_viscocity_stress_tangent);

const Number Dh = 3 * Carreau_viscocity * strain_rate > softened_quasistatic_elastoplastic_stress? 3 * Carreau_viscocity * strain_rate_tangent

3 * Carreau_viscocity_stress_tangent * stress_tangent * strain_rate
3 * Carreau_viscocity_strain_rate_tangent * strain_rate_tangent * strain_rate : 0;

        const Number Dh =
          3 * Carreau_viscocity * strain_rate_tangent
          + 3 * Carreau_viscocity_stress_tangent * stress_tangent * strain_rate
          + 3 * Carreau_viscocity_strain_rate_tangent * strain_rate_tangent * strain_rate;
        D_isotropic_hardening = beta * Dh;
        D_kinematic_hardening = (1.0 - beta) * Dh;
        return Dh;
      } else {
        const Number Dh =
            get_elastoplastic_factor_tangent(alpha)
                * get_viscosity_factor(gamma, time_increment)
                * get_softening_factor(temperature)
            + get_elastoplastic_factor(alpha)
                * get_viscosity_factor_tangent(gamma, time_increment)
                * get_softening_factor(temperature)
            + viscosity_regularization_factor/time_increment;
        D_isotropic_hardening = beta * Dh;
        D_kinematic_hardening = (1.0 - beta) * Dh;
        return Dh;
      }
    }
  
    Number hardening_temperature_derivatives(Number &D_isotropic_hardening,
                                                  Number &D_kinematic_hardening,
                                                  const Number alpha,
                                                  const Number gamma,
                                                  const Number time_increment,
                                                  const Number temperature) const {
      if(use_Carreau_viscous_law) {
        const Number creep_strain_rate_factor = std::pow(1-std::max(0., std::min(1., (temperature - reference_temperature)/(melting_temperature - reference_temperature))), 1.5);
        const Number minimum_strain_rate = creep_strain_rate_factor * epsilon_dot_0;
        const Number strain_rate = std::max(sqrt2thirds*gamma/time_increment, minimum_strain_rate);
        const Number softened_quasistatic_elastoplastic_stress = get_elastoplastic_factor(alpha) * get_softening_factor(temperature);
        [[maybe_unused]] const Number Carreau_viscocity = get_Carreau_viscocity(strain_rate, softened_quasistatic_elastoplastic_stress);
        const Number stress_temperature_tangent =
                get_elastoplastic_factor(alpha) * get_softening_factor_tangent(temperature);
        Number Carreau_viscocity_strain_rate_tangent, Carreau_viscocity_stress_tangent;
        get_Carreau_viscocity_tangents(
            strain_rate, softened_quasistatic_elastoplastic_stress,
            Carreau_viscocity_strain_rate_tangent,
            Carreau_viscocity_stress_tangent);
        const Number Dh = 3 * Carreau_viscocity_stress_tangent * stress_temperature_tangent * strain_rate;
        D_isotropic_hardening = beta * Dh;
        D_kinematic_hardening = (1 - beta) * Dh;
        return Dh;
      } else {
        Number Dh =
            get_elastoplastic_factor(alpha)
                * get_viscosity_factor(gamma, time_increment)
                * get_softening_factor_tangent(temperature);
        D_isotropic_hardening = beta * Dh;
        D_kinematic_hardening = (1 - beta) * Dh;
        return Dh;
      }
    }
  
    Number trial_yield_criterion(const Number norm_ksi_trial,
                          const Number alpha,
                          const Number gamma,
                          const Number time_increment,
                          const Number temperature) const {
      Number K_alpha, H_alpha;
      hardening_values(H_alpha, K_alpha, alpha, gamma, time_increment, temperature);
      const Number trial_yield_criterion = norm_ksi_trial - sqrt2thirds * (H_alpha);
      return trial_yield_criterion;
    }
  
    Number get_elastoplastic_factor(const Number alpha) const {
      if(alpha > max_strain) {
        return A + B * std::pow(max_strain, n);
      }
      if(alpha >= eps) {
          return A + B * std::pow(alpha, n);
      }
      return A + B * alpha/eps * std::pow(eps, n);
    }
  
    Number get_viscosity_factor(const Number gamma, const Number time_increment) const {
      Number slope, intercept;
      get_small_hardening_fit(slope, intercept, time_increment);
      if(gamma >= intercept) {
          return 1.0 + C * std::log(sqrt2thirds*gamma/(time_increment*reference_strain_rate));
      } else if (gamma < 0.0) {
          return 1.0 - C * slope * sqrt2thirds * gamma * gamma;
      }
      return 1.0 + C * slope * sqrt2thirds * gamma * gamma;
    }
  
    Number get_softening_factor(const Number temperature) const {
      if(temperature > reference_temperature) {
        if(temperature < melting_temperature) {
          return (1.0 + softening_threshold - std::pow((temperature - reference_temperature)/(melting_temperature - reference_temperature), m));
        } else {
          return 0.0 + softening_threshold;
        }
      }
      return 1.0 + softening_threshold;
    }
  
    Number get_elastoplastic_factor_tangent(const Number alpha) const {
      if(alpha > max_strain) {
        return 0;
      }
      if(alpha >= eps) {
          return B * n * std::pow(alpha, n-1.0);
      }
      return B * 1.0/eps * std::pow(eps, n);
    }
  
    Number get_viscosity_factor_tangent(const Number gamma, const Number time_increment) const {
      Number slope, intercept;
      get_small_hardening_fit(slope, intercept, time_increment);
      if(gamma >= intercept) {
          return C / (sqrt2thirds * gamma);
      } else if (gamma < 0.0) {
          return -2 * C * slope * gamma;
      }
      return 2 * C * slope * gamma;
    }
  
    Number get_softening_factor_tangent(const Number temperature) const {
      if(temperature > reference_temperature) {
        if(temperature < melting_temperature) {
          return (-m/(melting_temperature - reference_temperature))
                  * std::pow((temperature - reference_temperature)/(melting_temperature - reference_temperature), m-1.0);
        } else {
          return 0.0;
        }
      }
      return 0.0;
    }
  
    void get_small_hardening_fit(Number &slope, Number &intercept, const Number time_increment) const {

the log factor is annoying when below 1.0. Replace it by a parabula till it behaves.

      const Number log_factor = 1./(time_increment*reference_strain_rate);
      intercept = exp_one_half/(sqrt2thirds * log_factor);
      slope = 1./(2*intercept*intercept*sqrt2thirds);
    }
  
    Number get_Carreau_viscocity(Number strain_rate, Number sigma_0_theta) const {
      if(sigma_0_theta <= 0) {
        return mu_infty;
      }
      const Number g_sigma_epsilon_dot = std::pow(sigma_0_theta/(3*epsilon_dot_0*mu_0), (n_C/(1-n_C))) * (strain_rate/epsilon_dot_0);
      return std::pow(1 + std::pow(g_sigma_epsilon_dot, 2), ((1-n_C)/(2*n_C))) * (mu_0 - mu_infty) + mu_infty;
    }
  
    void get_Carreau_viscocity_tangents(
            Number strain_rate,
            Number sigma_0_theta,
            Number &Carreau_viscocity_strain_rate_tangent,
            Number &Carreau_viscocity_stress_tangent) const {
      if(sigma_0_theta <= 0) {
        Carreau_viscocity_strain_rate_tangent = 0;
        Carreau_viscocity_stress_tangent = 0;
        return;
      }
      const Number g_sigma_epsilon_dot = std::pow(sigma_0_theta/(3*epsilon_dot_0*mu_0), (n_C/(1-n_C))) * (strain_rate/epsilon_dot_0);
      const Number g_sigma_epsilon_dot_strain_rate_tangent = std::pow(sigma_0_theta/(3*epsilon_dot_0*mu_0), (n_C/(1-n_C))) * (1/epsilon_dot_0);
      const Number g_sigma_epsilon_dot_stress_tangent =
        (n_C / (1-n_C)) * std::pow(sigma_0_theta/(3*epsilon_dot_0*mu_0), ((2*n_C-1)/(1-n_C))) * (strain_rate/epsilon_dot_0) * (1/(3 * epsilon_dot_0 * mu_0));
      const Number Carreau_viscocity_g_tangent =
        ((1-n_C)/(2*n_C)) * std::pow(1 + std::pow(g_sigma_epsilon_dot, 2), ((1-3*n_C)/(2*n_C))) * (2*g_sigma_epsilon_dot) * (mu_0 - mu_infty);
  
      Carreau_viscocity_strain_rate_tangent = Carreau_viscocity_g_tangent * g_sigma_epsilon_dot_strain_rate_tangent;
      Carreau_viscocity_stress_tangent = Carreau_viscocity_g_tangent * g_sigma_epsilon_dot_stress_tangent;
    }
  
    private:
      const Number mu;
      const Number A;
      const Number B;
      const Number C;
      const Number m;
      const Number n;
      const Number melting_temperature;
      const Number reference_strain_rate;
      const Number reference_temperature;
      const Number beta;
      const Number sqrt2thirds;
      const Number exp_one_half;
  
      const Number eps = std::pow(B/mu, 1.0/(1.0-n));
      const Number softening_threshold = 0.0;
      const Number viscosity_regularization_factor = 0;
      const Number max_strain = std::numeric_limits<Number>::max();
  

Carreau fluid parameters

      const bool use_Carreau_viscous_law = false;
      const Number epsilon_dot_0 = 1;
      const Number n_C = 3;
      const Number mu_0 = 1e18;
      const Number mu_infty = 1e-4;
  
    };
  
  } /* namespace PlasticityLab */
  
  #endif /* JOHNSONCOOKTHERMOVISCOPLASTICYIELDLAW_H_ */ * 

Annotated version of src/LBCSystem.h

  /*
   * LBCSystem.h
   *
   *  Created on: 08 Oct 2019
   *      Author: maien
   */
  
  #ifndef LBCSYSTEM_H_
  #define LBCSYSTEM_H_
  
  #include <deal.II/base/function.h>
  
  #include "InterpolatoryConstraintApplier.h"
  #include "BodyForceApplier.h"
  #include "ConvectionBoundaryConditionApplier.h"
  #include "IncrementInterpolationHandler.h"
  #include "utilities.h"
  #include "DoFSystem.h"
  #include "BoundaryUnidirectionalPenaltySpec.h"
  
  using namespace dealii;
  using namespace Functions;
  
  namespace PlasticityLab {
  
  template <int dim, typename Number=double, int components=dim>
  class LBCSystem {
  public:
      LBCSystem(): zero_function(components){}
      ~LBCSystem() {
          for(auto increment_interpolation_handler: increment_interpolation_handlers) {
              delete increment_interpolation_handler;
          }
          for(auto initial_velocity_interpolation_handler: initial_velocity_interpolation_handlers) {
              delete initial_velocity_interpolation_handler;
          }
          for(auto initial_deformation_interpolation_handler: initial_deformation_interpolation_handlers) {
              delete initial_deformation_interpolation_handler;
          }
          for(auto boundary_unidirectional_penalty_spec: boundary_unidirectional_penalty_specs) {
              delete boundary_unidirectional_penalty_spec;
          }
      }
  
      void apply_constraints(DoFSystem<dim, Number> &dof_system) const;
  
      void clear();
  
      std::vector< BodyForceApplier<dim, Number> > bodyLoadAppliers;
      std::vector< std::pair<int, BodyForceApplier<dim, Number> > > boundaryLoadAppliers;
      std::vector< std::pair<int, ConvectionBoundaryConditionApplier<dim, Number> > > convection_BC_appliers;
      std::vector< InterpolatoryConstraintApplier<dim, Number> > interpolatoryConstraintAppliers;
      std::vector<std::pair<unsigned int, std::set<types::boundary_id>>> no_normal_flux_constraints;
      std::vector<IncrementInterpolationHandler<dim, Number, components>*> increment_interpolation_handlers;
      std::vector<IncrementInterpolationHandler<dim, Number, components>*> initial_velocity_interpolation_handlers;
      std::vector<IncrementInterpolationHandler<dim, Number, components>*> initial_deformation_interpolation_handlers;
      std::vector<BoundaryUnidirectionalPenaltySpec<Number>*> boundary_unidirectional_penalty_specs;
  
      ZeroFunction<dim, Number>     zero_function;
  };
  
  
  template<int dim, typename Number, int components>
  void LBCSystem<dim, Number, components>::apply_constraints(DoFSystem<dim, Number> &dof_system) const {
      for (auto constraintApplier = interpolatoryConstraintAppliers.cbegin();
            constraintApplier != interpolatoryConstraintAppliers.end();
            ++constraintApplier) {
        constraintApplier->apply(dof_system.mapping, dof_system.dof_handler, dof_system.nodal_constraints);
      }
  
      for (auto no_normal_flux_constraint : no_normal_flux_constraints) {
        VectorTools::compute_no_normal_flux_constraints(
          dof_system.dof_handler,
          no_normal_flux_constraint.first,
          no_normal_flux_constraint.second,
          dof_system.nodal_constraints,
          dof_system.mapping);
      }
  
      dof_system.nodal_constraints.close();
  }
  
  
  template<int dim, typename Number, int components>
  void LBCSystem<dim, Number, components>::clear() {
      for(auto increment_interpolation_handler: increment_interpolation_handlers) {
          delete increment_interpolation_handler;
      }
      for(auto initial_velocity_interpolation_handler: initial_velocity_interpolation_handlers) {
          delete initial_velocity_interpolation_handler;
      }
      for(auto initial_deformation_interpolation_handler: initial_deformation_interpolation_handlers) {
          delete initial_deformation_interpolation_handler;
      }
      for(auto boundary_unidirectional_penalty_spec: boundary_unidirectional_penalty_specs) {
          delete boundary_unidirectional_penalty_spec;
      }
  
      bodyLoadAppliers.clear();
      boundaryLoadAppliers.clear();
      convection_BC_appliers.clear();
      interpolatoryConstraintAppliers.clear();
      no_normal_flux_constraints.clear();
      increment_interpolation_handlers.clear();
      initial_velocity_interpolation_handlers.clear();
      initial_deformation_interpolation_handlers.clear();
      boundary_unidirectional_penalty_specs.clear();
  
  }
  
  
  } /* namespace PlasticityLab */
  
  #endif /* LBCSYSTEM_H_ */

Annotated version of src/Material.h

  /*
   * Material.h
   *
   *  Created on: 05 Jan 2015
   *      Author: maien
   */
  
  #ifndef MATERIAL_H_
  #define MATERIAL_H_
  
  using namespace dealii;
  
  #include "ConstitutiveModelRequest.h"
  #include <stdexcept>
  
  namespace PlasticityLab {
  
    typedef size_t point_index_t;
  
    template <int dim, typename Number = double>
    class Material {
     public:
      virtual ~Material() = 0;
  
      virtual void compute_constitutive_request(
        ConstitutiveModelRequest <dim, Number> &constitutive_request,
        const point_index_t &point_index) = 0;
  
      virtual Number get_material_Jacobian(const point_index_t &point_index) const = 0;
      virtual ::SymmetricTensor<2, dim, Number> get_plastic_strain(const point_index_t &point_index) const = 0;
  
      virtual void setup_point_history(const point_index_t point_count) = 0;
  
      virtual std::vector<Number> get_state_parameters(
                  const point_index_t &point_index,
                  const Tensor<2, dim, Number> &reference_transformation=unit_symmetric_tensor<dim>()) const = 0;
  
      virtual void set_state_parameters(
                  const point_index_t &point_index,
                  const std::vector<Number> &state_parameters,
                  const Tensor<2, dim, Number> &reference_transformation) = 0;
      virtual size_t get_material_parameter_count() const = 0;
  
    };
  
    template <int dim, typename Number>
    inline Material<dim, Number>::~Material() {
    }
  
    class MaterialDomainException: public std::runtime_error {
     public:
      MaterialDomainException();
      MaterialDomainException(std::string s): std::runtime_error(s) {};
    };
  
  } /* namespace PlasticityLab */
  
  #endif /* MATERIAL_H_ */

Annotated version of src/MixedFEProjector.h

  /*
   * MixedFEProjector.h
   *
   *  Created on: 21 Jul 2015
   *      Author: maien
   */
  
  #ifndef MIXEDFEPROJECTOR_H_
  #define MIXEDFEPROJECTOR_H_
  
  #include <vector>
  
  #include <deal.II/base/tensor.h>
  #include <deal.II/fe/fe_values.h>
  
  #include <deal.II/lac/vector.h>
  
  #include "utilities.h"
  
  template <int dim, typename Number = double>
  class MixedFEProjector {
   public:
    MixedFEProjector();
    MixedFEProjector(
      const unsigned int mixed_dofs_per_cell,
      const ::FEValues<dim> &mixed_fe_values);
    virtual ~MixedFEProjector();
  
    template <typename T>
    void project(
      std::vector<T> *coefficients_of_mixed_dofs,
      const std::vector<T> &values_at_q_points) const;
  
   private:
    unsigned int n_q_points;
    unsigned int mixed_dofs_per_cell;
    std::vector<std::vector<Number > > M_inv_ksi;
  };
  
  template <int dim, typename Number>
  MixedFEProjector<dim, Number>::MixedFEProjector():
    n_q_points(0),
    mixed_dofs_per_cell(0),
    M_inv_ksi(0) {
  
  }
  
  template <int dim, typename Number>
  MixedFEProjector<dim, Number>::MixedFEProjector(
    const unsigned int mixed_dofs_per_cell,
    const ::FEValues<dim> &mixed_fe_values)
    : n_q_points (mixed_fe_values.get_quadrature().size()),
      mixed_dofs_per_cell (mixed_dofs_per_cell),
      M_inv_ksi (n_q_points, std::vector<Number>(mixed_dofs_per_cell)) {
    ::FullMatrix<Number>   M_matrix(mixed_dofs_per_cell, mixed_dofs_per_cell),
           M_inv(mixed_dofs_per_cell, mixed_dofs_per_cell);
    std::vector<::Vector<Number> > ksi(n_q_points, ::Vector<Number>(mixed_dofs_per_cell));
  
    M_matrix = 0;
    for (unsigned int q_point = 0; q_point < n_q_points;
         ++q_point) {

Prep to compute mixed primary variables (Simo & Miehe 1992)

      for (unsigned int i = 0; i < mixed_dofs_per_cell; ++i) {
        const Number i_value = mixed_fe_values.shape_value (i, q_point);
        for (unsigned int j = 0; j < mixed_dofs_per_cell; ++j) {
          const Number j_value = mixed_fe_values.shape_value (j, q_point);
          M_matrix(i, j) += i_value * j_value * mixed_fe_values.quadrature_point(q_point)[0] * mixed_fe_values.JxW(q_point);
        }
        ksi.at(q_point)[i] = i_value * mixed_fe_values.quadrature_point(q_point)[0] * mixed_fe_values.JxW(q_point);
      }
    }
  
    M_inv.invert(M_matrix);
  
    for (unsigned int q_point = 0; q_point < n_q_points;
         ++q_point) {
      ::Vector<Number> M_inv_ksi_at_q_point(mixed_dofs_per_cell);
      M_inv.vmult(M_inv_ksi_at_q_point, ksi.at(q_point), false);
      for (unsigned int i = 0; i < mixed_dofs_per_cell; ++i)
        M_inv_ksi[q_point][i] = M_inv_ksi_at_q_point(i);
    }
  }
  
  template <int dim, typename Number>
  MixedFEProjector<dim, Number>::~MixedFEProjector() {}
  
  template <int dim, typename Number>
  template <typename T>
  void MixedFEProjector<dim, Number>::project(
    std::vector<T> *coefficients_of_mixed_dofs,
    const std::vector<T> &values_at_q_points) const {
    for (unsigned int i = 0; i < mixed_dofs_per_cell; ++i) {
      coefficients_of_mixed_dofs->at(i) = M_inv_ksi[0][i] * values_at_q_points[0];
      for (unsigned int q_point = 1; q_point < n_q_points; ++q_point)
        coefficients_of_mixed_dofs->at(i) += M_inv_ksi[q_point][i] * values_at_q_points[q_point];
    }
  }
  
  #endif /* MIXEDFEPROJECTOR_H_ */

Annotated version of src/NewtonStepSystem.h

  /*
   * NewtonStepSystem.h
   *
   *  Created on: 05 May 2015
   *      Author: maien
   */
  
  #ifndef NEWTONSTEPSYSTEM_H_
  #define NEWTONSTEPSYSTEM_H_
  
  
  #include <deal.II/lac/trilinos_sparse_matrix.h>
  #include <deal.II/lac/trilinos_vector.h>
  #include <deal.II/lac/sparsity_tools.h>
  
  #include "DoFSystem.h"
  #include "Constants.h"
  
  using namespace dealii;
  
  namespace PlasticityLab {
  
    class NewtonStepSystem {
     public:
  
      template<class DoFSystemType>
      void setup(const DoFSystemType &dof_system) {
        TrilinosWrappers::SparsityPattern sparsity_pattern(
          dof_system.locally_owned_dofs,
          mpi_communicator);
  
        DoFTools::make_sparsity_pattern(
          dof_system.dof_handler, sparsity_pattern,
          dof_system.nodal_constraints,
          false,
          Utilities::MPI::this_mpi_process(mpi_communicator));
  
        sparsity_pattern.compress();
        Newton_step_matrix.reinit(sparsity_pattern);
        Newton_step_solution.reinit(
          dof_system.locally_owned_dofs,
          mpi_communicator);
        current_increment.reinit(
          dof_system.locally_owned_dofs,
          dof_system.locally_relevant_dofs,
          mpi_communicator);
        Newton_step_residual.reinit(
          dof_system.locally_owned_dofs,
          mpi_communicator);
        previous_deformation.reinit(
          dof_system.locally_owned_dofs,
          dof_system.locally_relevant_dofs,
          mpi_communicator);
  
        previous_time_derivative.reinit(
          dof_system.locally_owned_dofs,
          dof_system.locally_relevant_dofs,
          mpi_communicator);
        previous_second_time_derivative.reinit(
          dof_system.locally_owned_dofs,
          dof_system.locally_relevant_dofs,
          mpi_communicator);
  
        _locally_owned_current_increment.reinit(
          dof_system.locally_owned_dofs,
          mpi_communicator);
        _locally_owned_previous_deformation.reinit(
          dof_system.locally_owned_dofs,
          mpi_communicator);
  
        _locally_owned_previous_time_derivative.reinit(
          dof_system.locally_owned_dofs,
          mpi_communicator);
        _locally_owned_previous_second_time_derivative.reinit(
          dof_system.locally_owned_dofs,
          mpi_communicator);
      }
  
      template<class DoFSystemType>
      void update_matrix_constraints(const DoFSystemType &dof_system) {
        TrilinosWrappers::SparsityPattern sparsity_pattern(
          dof_system.locally_owned_dofs,
          mpi_communicator);
  
        DoFTools::make_sparsity_pattern(
          dof_system.dof_handler, sparsity_pattern,
          dof_system.nodal_constraints,
          false,
          Utilities::MPI::this_mpi_process(mpi_communicator));
  
        sparsity_pattern.compress();
        Newton_step_matrix.reinit(sparsity_pattern);
      }
  
      void advance_time(double delta_t, double rho_infty, const bool reset_increment=true) {
        {
          double alpha_m, alpha_f, gamma, beta;
          get_generalized_alpha_method_params(
              &alpha_m, &alpha_f, &gamma, &beta, rho_infty);
  
          _locally_owned_previous_time_derivative = previous_time_derivative;
          _locally_owned_previous_second_time_derivative = previous_second_time_derivative;
  
          const TrilinosWrappers::MPI::Vector previous_second_time_derivative_backup(_locally_owned_previous_second_time_derivative);
  
          _locally_owned_previous_second_time_derivative = current_increment;
          _locally_owned_previous_second_time_derivative.add(
                -delta_t,
                _locally_owned_previous_time_derivative,
                -delta_t*delta_t*(0.5-beta),
                previous_second_time_derivative_backup);
          _locally_owned_previous_second_time_derivative *= (1./(beta*delta_t*delta_t));
          _locally_owned_previous_time_derivative.add(
                delta_t*(1.-gamma),
                previous_second_time_derivative_backup,
                delta_t*gamma,
                _locally_owned_previous_second_time_derivative);
  
  
          previous_time_derivative = _locally_owned_previous_time_derivative;
          previous_second_time_derivative = _locally_owned_previous_second_time_derivative;
  
        }
        

Update the deformation vector with the computed increment.

        add_current_increment_to_previous_deformation();
  
        if(reset_increment) {
          current_increment = 0;
        }
      }
  

Add the current Newton increment into the deformation vector, i.e. compute previous_deformation += current_increment.

'previous_deformation' is a vector with ghost entries and therefore read-only: we are not allowed to write into it (with the exception of setting it to zero). We therefore carry out the arithmetic in fully-distributed (locally-owned) temporary vectors and only assign the result back into the ghosted vector at the very end; that assignment performs the necessary ghost-value communication.

      void add_current_increment_to_previous_deformation() {
        _locally_owned_previous_deformation = previous_deformation;
        _locally_owned_current_increment    = current_increment;
        _locally_owned_previous_deformation += _locally_owned_current_increment;
        previous_deformation = _locally_owned_previous_deformation;
      }
  

Set the deformation vector to the negative of the current increment, i.e. compute previous_deformation = -current_increment.

As above, 'previous_deformation' is a ghosted, read-only vector, so the negation is performed in a fully-distributed temporary and only the result is assigned back into the ghosted vector.

      void set_previous_deformation_to_negative_current_increment() {
        _locally_owned_current_increment = current_increment;
        _locally_owned_current_increment *= -1;
        previous_deformation = _locally_owned_current_increment;
      }
  
      TrilinosWrappers::SparseMatrix  Newton_step_matrix;
      TrilinosWrappers::MPI::Vector   previous_deformation;
      TrilinosWrappers::MPI::Vector   current_increment;
      TrilinosWrappers::MPI::Vector   Newton_step_solution;
      TrilinosWrappers::MPI::Vector   Newton_step_residual;
  
      TrilinosWrappers::MPI::Vector   previous_time_derivative;
      TrilinosWrappers::MPI::Vector   previous_second_time_derivative;
  
      TrilinosWrappers::MPI::Vector   _locally_owned_previous_deformation;
      TrilinosWrappers::MPI::Vector   _locally_owned_current_increment;
      TrilinosWrappers::MPI::Vector   _locally_owned_previous_time_derivative;
      TrilinosWrappers::MPI::Vector   _locally_owned_previous_second_time_derivative;
    };
  
  } /* namespace PlasticityLab */
  
  #endif /* NEWTONSTEPSYSTEM_H_ */

Annotated version of src/PlasticityLabProg.cpp

  /*
   * PlasticityLabProg.cpp
   *
   *  Created on: 09 Jul 2014
   *      Author: cerecam
   */
  
  #include <deal.II/lac/sparsity_tools.h>
  #include <deal.II/lac/solver_bicgstab.h>
  #include <deal.II/lac/solver_gmres.h>
  #include <deal.II/lac/precondition.h>
  #include <deal.II/lac/trilinos_block_sparse_matrix.h>
  #include <deal.II/lac/trilinos_parallel_block_vector.h>
  #include <deal.II/lac/trilinos_precondition.h>
  #include <deal.II/lac/trilinos_solver.h>
  
  #include <deal.II/lac/affine_constraints.h>
  
  #include "TimeRateUpdateFlags.h"
  #include "TimeRateRequest.h"
  
  #include "PlasticityLabProg.h"
  #include "PlasticityLabProgDrivers.cpp"
  
  #include "ReferencePoint.h"
  #include "RemappedPoint.h"
  
  using namespace dealii;
  using std::endl;
  
  namespace PlasticityLab {
    template <int dim, typename Number>
    PlasticityLabProg<dim, Number>::PlasticityLabProg(
        Material<dim+1, Number> &material) :
      mpi_communicator(MPI_COMM_WORLD),
      pcout(std::cout,
            (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
      order(1),
      mech_fe(FE_Q<dim>(order), dim, FE_Q<dim>(order), 1),
      therm_fe(order),
      mixed_var_fe(order-1),
      mesh_motion_fe(FE_Q<dim>(order), dim),
      mapping (order),
      triangulation(mpi_communicator),
      mech_dof_system(triangulation, mapping),
      therm_dof_system (triangulation, mapping),
      mixed_fe_dof_system(triangulation, mapping),
      mesh_motion_dof_system(triangulation, mapping),
      quadrature_formula(order+1),
      face_quadrature_formula(order+1),
      material(material) { }
  
  
    template <int dim, typename Number>
    PlasticityLabProg<dim, Number>::~PlasticityLabProg() {
  
    }
  
    template <int dim, typename Number>
    template <typename TriangulationType, typename MaterialType>
    void PlasticityLabProg<dim, Number>::setup_material_data(
      TriangulationType   &triangulation,
      MaterialType &material) {
      const unsigned int num_cells = triangulation.n_active_cells();
      triangulation.clear_user_data();
      material.setup_point_history(num_cells * quadrature_formula.size());
      unsigned int history_index = 0;
      for (typename Triangulation<dim>::active_cell_iterator
           cell = triangulation.begin_active();
           cell != triangulation.end(); ++cell) {
        cell->set_user_index (history_index);
        history_index += quadrature_formula.size();
      }
    }
  
    template <int dim, typename Number>
    void PlasticityLabProg<dim, Number>::setup_material_area_factors(
        const DoFSystem<dim, Number> &mesh_motion_dof_system,
        std::unordered_map<size_t, Tensor<1, dim+1, Number>> &material_area_factors) {
  
      FEFaceValues<dim> fe_face_values (
        mapping,
        mesh_motion_fe,
        face_quadrature_formula,
        update_normal_vectors);
      const unsigned int n_face_q_points = face_quadrature_formula.size();
  
      for (auto cell = mesh_motion_dof_system.dof_handler.begin_active();
           cell != mesh_motion_dof_system.dof_handler.end();
           ++cell) {
        if (cell->is_locally_owned()) {
          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face) {
            if(cell->at_boundary(face)) {
              fe_face_values.reinit(cell, face);
  
              for (unsigned int q_point = 0; q_point < n_face_q_points; q_point++) {
                const unsigned int cell_index = cell->user_index() / quadrature_formula.size();
                const unsigned int surface_point_key =
                  cell_index * GeometryInfo<dim>::faces_per_cell * n_face_q_points
                  + face * n_face_q_points
                  + q_point;
                material_area_factors[surface_point_key] = postprocess_tensor_dimension(fe_face_values.normal_vector(q_point), 0);
              }
            }
          }
        }
      }
    }
  
    template <int dim, typename Number>
    void PlasticityLabProg<dim, Number>::update_material_area_factors(
        const NewtonStepSystem &mesh_motion_nonlinear_system,
        const DoFSystem<dim, Number> &mesh_motion_dof_system,
        std::unordered_map<size_t, Tensor<1, dim+1, Number>> &material_area_factors) {
  
      FEFaceValues<dim> fe_face_values (
        mapping,
        mesh_motion_fe,
        face_quadrature_formula,
        update_values | update_gradients | update_quadrature_points);
  
      const unsigned int n_face_q_points = face_quadrature_formula.size();
  
      const FEValuesExtractors::Vector displacements (0);
  
      std::vector< Tensor<1, dim, Number> > mesh_motion_value_increments(n_face_q_points);
      std::vector< Tensor<2, dim, Number> > mesh_motion_gradient_increments(n_face_q_points);
  
      for (auto cell = mesh_motion_dof_system.dof_handler.begin_active();
           cell != mesh_motion_dof_system.dof_handler.end();
           ++cell) {
        if (cell->is_locally_owned()) {
          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face) {
            if(cell->at_boundary(face)) {
              fe_face_values.reinit(cell, face);
  
              fe_face_values[displacements].get_function_gradients(
                mesh_motion_nonlinear_system.current_increment,
                mesh_motion_gradient_increments);
  
              fe_face_values[displacements].get_function_values(
                mesh_motion_nonlinear_system.current_increment,
                mesh_motion_value_increments);
  
              for (unsigned int q_point = 0; q_point < n_face_q_points; q_point++) {
                const unsigned int cell_index = cell->user_index() / quadrature_formula.size();
                const unsigned int surface_point_key =
                  cell_index * GeometryInfo<dim>::faces_per_cell * n_face_q_points
                  + face * n_face_q_points
                  + q_point;
  
                const auto mesh_motion_gradient = get_deformation_gradient(
                                                    -mesh_motion_gradient_increments[q_point],
                                                    -mesh_motion_value_increments[q_point][0]
                                                      /fe_face_values.quadrature_point(q_point)[0]);
  
                material_area_factors[surface_point_key] = std::pow(determinant(mesh_motion_gradient), -1)
                                                           * transpose(mesh_motion_gradient)
                                                           * material_area_factors[surface_point_key];
              }
            }
          }
        }
      }
    }
  
    template <int dim, typename Number>
    template <typename TriangulationType>
    void PlasticityLabProg<dim, Number>::setup_mixed_fe_projection_data(
        const TriangulationType                      &triangulation,
        std::vector< MixedFEProjector<dim, Number> > &MixedFeProjectors,
        const FiniteElement<dim>                     &MixedFE,
        const Quadrature<dim>                        &quadrature_formula) {
  
      FEValues<dim> mixed_fe_values(
        mapping,
        MixedFE,
        quadrature_formula,
        update_values | update_quadrature_points | update_JxW_values);
  
      const unsigned int num_cells = triangulation.n_active_cells(),
                         n_q_points = quadrature_formula.size(),
                         mixed_dofs_per_cell = MixedFE.dofs_per_cell;
      MixedFeProjectors.clear();
      MixedFeProjectors.resize(num_cells);
      for (typename Triangulation<dim>::active_cell_iterator
           cell = triangulation.begin_active();
           cell != triangulation.end(); ++cell) {
        mixed_fe_values.reinit(cell);
        MixedFeProjectors.at(cell->user_index() / n_q_points) =
          MixedFEProjector<dim, Number>(mixed_dofs_per_cell,
                                        mixed_fe_values);
      }
    }
  
    template<int dim, typename Number>
    void PlasticityLabProg<dim, Number>::remap_material_state_variables(
        const NewtonStepSystem &mesh_motion_nonlinear_system,
        const DoFSystem<dim, Number> &mesh_motion_dof_system,
        const NewtonStepSystem &mechanical_nonlinear_system,
        const DoFSystem<dim, Number> &mechanical_dof_system,
        const DoFSystem<dim, Number> &mixed_fe_dof_system,
        const std::vector< MixedFEProjector<dim, Number> > &mixed_fe_projector,
        Material<dim+1, Number> &material,
        std::unordered_map<point_index_t, Tensor<2, dim+1, Number>> &remapped_deformation_gradients) {
  
      const unsigned int n_q_points = quadrature_formula.size();
      [[maybe_unused]] const unsigned int dofs_per_cell = mesh_motion_fe.dofs_per_cell;
      const unsigned int mixed_dofs_per_cell = mixed_fe_dof_system.dof_handler.get_fe().dofs_per_cell;
  
      FEValues<dim> mesh_motion_fe_values(
        mapping,
        mesh_motion_fe,
        quadrature_formula,
        update_values  | update_gradients | update_quadrature_points | update_jacobians | update_JxW_values);
  
      FEValues<dim> mechanical_fe_values(
        mapping,
        mech_fe,
        quadrature_formula,
        update_values  | update_gradients);
  
      FEValues<dim> mixed_fe_values(
        mapping,
        mixed_fe_dof_system.dof_handler.get_fe(),
        quadrature_formula,
        update_values);
  
      std::vector< Tensor<1, dim, Number> > mesh_motion_value_increments(n_q_points);
      std::vector< Tensor<2, dim, Number> > mesh_motion_gradient_increments(n_q_points);
  
      std::vector< Tensor<1, dim, Number> > previous_deformation_values(n_q_points);
      std::vector< Tensor<2, dim, Number> > previous_deformation_gradients(n_q_points);
      std::vector< Tensor<1, dim, Number> > previous_deformation_value_at_remapped_point(1);
      std::vector< Tensor<2, dim, Number> > previous_deformation_gradient_at_remapped_point(1);
  
      const unsigned int material_parameter_count = material.get_material_parameter_count();
  
      const FEValuesExtractors::Vector displacements(0);
      const FEValuesExtractors::Scalar temperature (0);
  
      std::vector<ReferencePoint<dim, Number>> reference_points;
      std::vector<Point<dim, Number>> remapped_point_positions;
  
      for (auto cell = mesh_motion_dof_system.dof_handler.begin_active();
                cell != mesh_motion_dof_system.dof_handler.end();
                ++cell) {
        if (cell->is_locally_owned()) {
          mesh_motion_fe_values.reinit(cell);
  
          mesh_motion_fe_values[displacements].get_function_values(
            mesh_motion_nonlinear_system.current_increment,
            mesh_motion_value_increments);
  
          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
            [[maybe_unused]] const point_index_t quadrature_point_index = cell->user_index() + q_point;
            const Point<dim, Number> reference_point_position = mesh_motion_fe_values.quadrature_point(q_point);
            const Point<dim, Number> remapped_point_position = reference_point_position - mesh_motion_value_increments[q_point];
  
            ReferencePoint<dim, Number> reference_point;
            reference_point.mesh_motion_cell = cell;
            reference_point.q_point = q_point;
            reference_point.reference_point = reference_point_position;
            reference_point.remapped_point = remapped_point_position;
  
            reference_points.push_back(reference_point);
            remapped_point_positions.push_back(remapped_point_position);
  
          }
  
        }
      }
  
      MPI_Barrier(mpi_communicator);
  
      MPI_Datatype PointType;
      MPI_Type_contiguous(dim, MPI_DOUBLE, &PointType);
      MPI_Type_commit(&PointType);
      std::vector<int> displs;
      std::vector< Point<dim, Number>> received_remapped_positions;
      std::vector<RemappedPoint<dim, Number>> remapped_points;
      int nprocesses, this_process;
      int num_reference_points = reference_points.size();
      MPI_Comm_size(mpi_communicator, &nprocesses);
      std::vector<int> reference_point_counts(nprocesses);
      MPI_Allgather(
        &num_reference_points,
        1, MPI_INT,
        &reference_point_counts[0],
        1, MPI_INT,
        mpi_communicator);
      MPI_Comm_rank(mpi_communicator, &this_process);
      displs.resize(nprocesses);
      displs[0] = 0;
      for (int i = 1; i < nprocesses; ++i) {
        displs[i] = displs[i - 1] + reference_point_counts[i - 1];
      }
  
      const unsigned int count_received_reference_points = displs[nprocesses - 1] + reference_point_counts[nprocesses - 1];
      received_remapped_positions.resize(count_received_reference_points);
      remapped_points.resize(count_received_reference_points);
      std::vector<char> this_process_owns_remapped_point(count_received_reference_points);
  
      MPI_Allgatherv(
        &remapped_point_positions[0],
        num_reference_points,
        PointType,
        &received_remapped_positions[0],
        &reference_point_counts[0],
        &displs[0],
        PointType,
        mpi_communicator);
  
      MPI_Barrier(mpi_communicator);
  
      for(unsigned int received_point_id=0; received_point_id < count_received_reference_points; received_point_id++) {
        auto point = received_remapped_positions[received_point_id];
        auto cell_and_point = GridTools::find_active_cell_around_point(
                            mapping,
                            mesh_motion_dof_system.dof_handler,
                            point);
        auto cell = cell_and_point.first;
        auto unit_cell_point = cell_and_point.second;

here, we're assuming that find_active_cell_around_point returns the same cell and point when called with different dof_handlers

        auto mechanical_cell_and_point = GridTools::find_active_cell_around_point(
                            mapping,
                            mechanical_dof_system.dof_handler,
                            point);
        auto mechanical_cell = mechanical_cell_and_point.first;
        auto mixed_fe_cell_and_point = GridTools::find_active_cell_around_point(
                            mapping,
                            mixed_fe_dof_system.dof_handler,
                            point);
        auto mixed_fe_cell = mixed_fe_cell_and_point.first;
  
        remapped_points[received_point_id].remapped_point = point;
        remapped_points[received_point_id].mesh_motion_cell = cell;
        remapped_points[received_point_id].unit_cell_point = unit_cell_point;
        remapped_points[received_point_id].field_cell = mechanical_cell;
        remapped_points[received_point_id].mixed_fe_cell = mixed_fe_cell;
        this_process_owns_remapped_point[received_point_id] = cell.state() == IteratorState::valid && cell->is_locally_owned()? 1 : 0;
      }
  
      MPI_Barrier(mpi_communicator);
  
      std::vector<char> remapped_point_candidates(num_reference_points * nprocesses);
      for (int process = 0; process < nprocesses; ++process) {
        if (reference_point_counts[process] > 0)
          MPI_Gather(
            &this_process_owns_remapped_point[displs[process]],
            reference_point_counts[process],
            MPI_C_BOOL,
            &remapped_point_candidates[0],
            reference_point_counts[process],
            MPI_C_BOOL,
            process,
            mpi_communicator);
      }
  
      std::vector<unsigned int> remapped_point_owning_process(num_reference_points);
  
      for (int i = 0; i < num_reference_points; ++i) {
        remapped_point_owning_process[i] = 0;
        for (int j = 1; j < nprocesses; ++j) {
          if (remapped_point_candidates[j * num_reference_points + i] == 1) {
            remapped_point_owning_process[i] = j;
            continue;
          }
        }
      }
  
      std::vector<RemappedPoint<dim, Number>> mapping_remapped_points;
      std::vector<unsigned int> mapping_reference_point_owning_process;
      std::vector<unsigned int> mapping_reference_point_index_at_remote_process;
  

This should really be a vector<bool>, but addresses of individual elements of vector<bool> cannot be taken. It's a template specialization to save space

      std::vector<char> remote_remapped_point_is_accepted(num_reference_points * nprocesses);
      std::vector<char> local_remapped_point_is_accepted(count_received_reference_points);
      for (int i = 0; i < num_reference_points; ++i) {
        for (unsigned int j = 0; j < static_cast<unsigned int>(nprocesses); ++j) {
          remote_remapped_point_is_accepted[j * num_reference_points + i] = (remapped_point_owning_process[i] == j) ? 1 : 0;
        }
      }
  
      MPI_Barrier(mpi_communicator);
  
      for (int process = 0; process < nprocesses; ++process) {
        if (reference_point_counts[process] > 0) {
          MPI_Scatter(
            &remote_remapped_point_is_accepted[0],
            reference_point_counts[process],
            MPI_C_BOOL,
            &local_remapped_point_is_accepted[displs[process]],
            reference_point_counts[process],
            MPI_C_BOOL,
            process,
            mpi_communicator);
        }
      }
  
      std::vector<unsigned int> mapping_remote_reference_point_counts(nprocesses, 0);
      for (int process = 0; process < nprocesses; ++process) {
        for (int i = 0; i < reference_point_counts[process]; ++i) {
          if (local_remapped_point_is_accepted[displs[process] + i]) {
            RemappedPoint<dim, Number> accepted_remapped_point = remapped_points[displs[process] + i];
  
            mapping_remapped_points.push_back(accepted_remapped_point);
            mapping_reference_point_owning_process.push_back(process);
            mapping_reference_point_index_at_remote_process.push_back(i);
            ++mapping_remote_reference_point_counts[process];
          }
        }
      }
  
      MPI_Barrier(mpi_communicator);
  
      std::vector<unsigned int> mapping_remote_remapped_point_counts(nprocesses);
      for (int i = 0; i < nprocesses; ++i) {
        MPI_Scatter(
          &mapping_remote_reference_point_counts[0],
          1, MPI_UNSIGNED,
          &mapping_remote_remapped_point_counts[i],
          1, MPI_UNSIGNED,
          i, mpi_communicator);
      }
  
      {
        std::vector<std::vector<Number> > local_state_parameter_groups(nprocesses);
        std::vector<std::vector<Number> > remote_state_parameter_groups(nprocesses);
  
        std::vector<std::vector<Number> > local_deformation_gradient_groups(nprocesses);
        std::vector<std::vector<Number> > remote_deformation_gradient_groups(nprocesses);
  
        for (int i = 0; i < nprocesses; ++i) {
          local_state_parameter_groups.at(i).resize(material_parameter_count * mapping_remote_reference_point_counts.at(i));
          remote_state_parameter_groups.at(i).resize(material_parameter_count * mapping_remote_remapped_point_counts.at(i));
  
          local_deformation_gradient_groups.at(i).resize((dim+1) * (dim+1) * mapping_remote_reference_point_counts.at(i));
          remote_deformation_gradient_groups.at(i).resize((dim+1) * (dim+1) * mapping_remote_remapped_point_counts.at(i));
        }
  
        std::vector<unsigned int> next_to_process(nprocesses, 0);
        for (unsigned int i = 0; i < mapping_remapped_points.size(); ++i) {
          const unsigned int group = mapping_reference_point_owning_process.at(i);
          const RemappedPoint<dim, Number> remapped_point = mapping_remapped_points.at(i);
  
          Quadrature<dim> remapped_point_quadrature(
            std::vector<Point<dim, Number>> (1, remapped_point.unit_cell_point));
  
          FEValues<dim> remapped_point_fe_values(
            mapping,
            mesh_motion_fe,
            remapped_point_quadrature,
            update_values | update_gradients | update_quadrature_points);
  
          FEValues<dim> remapped_point_mechanical_fe_values(
            mapping,
            mech_fe,
            remapped_point_quadrature,
            update_values | update_gradients);
  
          FEValues<dim> remapped_point_mixed_fe_values(
            mapping,
            mixed_fe_dof_system.dof_handler.get_fe(),
            remapped_point_quadrature,
            update_values);
  
          remapped_point_fe_values.reinit(remapped_point.mesh_motion_cell);
          mesh_motion_fe_values.reinit(remapped_point.mesh_motion_cell);
  
          remapped_point_mechanical_fe_values.reinit(remapped_point.field_cell);
          mechanical_fe_values.reinit(remapped_point.field_cell);
  
          remapped_point_mixed_fe_values.reinit(remapped_point.mixed_fe_cell);
  
          mesh_motion_fe_values[displacements].get_function_gradients(
            mesh_motion_nonlinear_system.current_increment,
            mesh_motion_gradient_increments);
  
          mesh_motion_fe_values[displacements].get_function_values(
            mesh_motion_nonlinear_system.current_increment,
            mesh_motion_value_increments);
  
          mechanical_fe_values[displacements].get_function_gradients(
            mechanical_nonlinear_system.previous_deformation,
            previous_deformation_gradients);
  
          mechanical_fe_values[displacements].get_function_values(
            mechanical_nonlinear_system.previous_deformation,
            previous_deformation_values);
  
          remapped_point_mechanical_fe_values[displacements].get_function_gradients(
            mechanical_nonlinear_system.previous_deformation,
            previous_deformation_gradient_at_remapped_point);
  
          remapped_point_mechanical_fe_values[displacements].get_function_values(
            mechanical_nonlinear_system.previous_deformation,
            previous_deformation_value_at_remapped_point);
  
          std::vector< std::vector<Number> > material_parameters_at_q_points(
                material_parameter_count,
                std::vector<Number>(n_q_points));
  
          for(unsigned int q_point=0; q_point<n_q_points; q_point++) {
            const point_index_t quadrature_point_index = remapped_point.mesh_motion_cell->user_index() + q_point;
  
            [[maybe_unused]] const auto deformation_gradient =
              get_deformation_gradient(
                previous_deformation_gradients[q_point],
                previous_deformation_values[q_point][0]/mesh_motion_fe_values.quadrature_point(q_point)[0]);
  
            std::vector<Number> state_parameters = material.get_state_parameters(quadrature_point_index, /*deformation_gradient*/unit_symmetric_tensor<dim+1, Number>());
  
            for(unsigned int parameter_index=0; parameter_index<material_parameter_count; parameter_index++) {
                material_parameters_at_q_points[parameter_index][q_point] = state_parameters[parameter_index];
            }
          }
  
          std::vector<std::vector<Number> > projected_material_parameters_coefficients(
                material_parameter_count,
                std::vector<Number>(mixed_dofs_per_cell));
  
          const unsigned int cell_index = remapped_point.mesh_motion_cell->user_index() / n_q_points;
          for(unsigned int parameter_index=0; parameter_index<material_parameter_count; parameter_index++) {
            mixed_fe_projector[cell_index].project(
              &projected_material_parameters_coefficients[parameter_index],
              material_parameters_at_q_points[parameter_index]);
          }
  
          Vector<Number> mixed_values (mixed_dofs_per_cell);
          for (unsigned int mixed_dof = 0; mixed_dof < mixed_dofs_per_cell; ++mixed_dof) {
            mixed_values(mixed_dof) = remapped_point_mixed_fe_values.shape_value(mixed_dof, 0);
          }
  
          std::vector<Number> projected_state_parameters(material_parameter_count, 0);
          for(unsigned int parameter_index=0; parameter_index<material_parameter_count; parameter_index++) {
            for (unsigned int mixed_dof = 0; mixed_dof < mixed_dofs_per_cell; mixed_dof++) {
              projected_state_parameters[parameter_index] += mixed_values(mixed_dof) * projected_material_parameters_coefficients[parameter_index][mixed_dof];
            }
          }
  
  
          for(unsigned int parameter_index=0; parameter_index<material_parameter_count; parameter_index++) {
            local_state_parameter_groups.at(group).at(next_to_process.at(group)*material_parameter_count + parameter_index) = projected_state_parameters[parameter_index];
          }
  
          const Tensor<2, dim+1, Number> previous_deformation_gradient =
            get_deformation_gradient(
              previous_deformation_gradient_at_remapped_point[0],
              previous_deformation_value_at_remapped_point[0][0]/remapped_point_fe_values.quadrature_point(0)[0]);
  
          for(unsigned int dim_i=0; dim_i<dim+1; dim_i++) {
            const unsigned int current_i_index = next_to_process.at(group) * (dim+1) + dim_i;
            for(unsigned int dim_j=0; dim_j<dim+1; dim_j++) {
              local_deformation_gradient_groups.at(group).at(current_i_index * (dim+1) + dim_j) = previous_deformation_gradient[dim_i][dim_j];
            }
          }
  
          ++next_to_process.at(group);
        }
  
  
  
        enum MessageFlag {
          MATERIAL_STATE_PARAMETER,
          REMAPPED_DEFORMATION_GRADIENT
        };
  
        const unsigned int
          material_state_parameter_requests_offset = 0,
          remapped_deformation_gradient_requests_offset = 1,
          request_array_size = 2;
  
        std::vector<MPI_Request> requests_vector(2 * nprocesses * request_array_size);
  
        for (int i = 0; i < nprocesses; ++i) {
          MPI_Isend(
            local_state_parameter_groups.at(i).data(),
            material_parameter_count * mapping_remote_reference_point_counts.at(i),
            MPI_DOUBLE, i, MATERIAL_STATE_PARAMETER,
            mpi_communicator,
            &requests_vector[i + nprocesses * material_state_parameter_requests_offset]);
        }
  
        for (int i = 0; i < nprocesses; ++i) {
          MPI_Isend(
            local_deformation_gradient_groups.at(i).data(),
            (dim+1) * (dim+1) * mapping_remote_reference_point_counts.at(i),
            MPI_DOUBLE, i, REMAPPED_DEFORMATION_GRADIENT,
            mpi_communicator,
            &requests_vector[i + nprocesses * remapped_deformation_gradient_requests_offset]);
        }
  
        for (int i = 0; i < nprocesses; ++i) {
          const unsigned int row_start = i + nprocesses * request_array_size;
  
          MPI_Irecv(
            remote_state_parameter_groups.at(i).data(),
            material_parameter_count * mapping_remote_remapped_point_counts.at(i),
            MPI_DOUBLE, i, MATERIAL_STATE_PARAMETER,
            mpi_communicator,
            &requests_vector[row_start + nprocesses * material_state_parameter_requests_offset]);
  
          MPI_Irecv(
            remote_deformation_gradient_groups.at(i).data(),
            (dim+1) * (dim+1) * mapping_remote_remapped_point_counts.at(i),
            MPI_DOUBLE, i, REMAPPED_DEFORMATION_GRADIENT,
            mpi_communicator,
            &requests_vector[row_start + nprocesses * remapped_deformation_gradient_requests_offset]);
  
        }
  
        std::vector<MPI_Status> statuses_vector(2 * request_array_size * nprocesses);
        MPI_Waitall(
          2 * request_array_size * nprocesses,
          &requests_vector[0],
          &statuses_vector[0]);
  
        next_to_process.clear();
        next_to_process.resize(nprocesses, 0);
        for (unsigned int i = 0; i < reference_points.size(); ++i) {
          unsigned int group = remapped_point_owning_process.at(i);
  
          ReferencePoint<dim, Number> &reference_point = reference_points[i];
  
          mesh_motion_fe_values.reinit(reference_point.mesh_motion_cell);
  
          std::vector<Number> remapped_state_parameters(material_parameter_count);
          for(unsigned int state_index=0; state_index<material_parameter_count; state_index++) {
            remapped_state_parameters[state_index] = remote_state_parameter_groups.at(group).at(next_to_process.at(group) * material_parameter_count + state_index);
          }
  
          Tensor<2, dim+1, Number> previous_deformation_gradient;
          for(unsigned int dim_i=0; dim_i<dim+1; dim_i++) {
            const unsigned int current_i_index = next_to_process.at(group) * (dim+1) + dim_i;
            for(unsigned int dim_j=0; dim_j<dim+1; dim_j++) {
              previous_deformation_gradient[dim_i][dim_j] = remote_deformation_gradient_groups.at(group).at(current_i_index * (dim+1) + dim_j);
            }
          }
  
          mesh_motion_fe_values[displacements].get_function_gradients(
            mesh_motion_nonlinear_system.current_increment,
            mesh_motion_gradient_increments);
  
          mesh_motion_fe_values[displacements].get_function_values(
            mesh_motion_nonlinear_system.current_increment,
            mesh_motion_value_increments);
  
          [[maybe_unused]] const auto mesh_motion_gradient = get_deformation_gradient(
                                              -mesh_motion_gradient_increments[reference_point.q_point],
                                              -mesh_motion_value_increments[reference_point.q_point][0]/mesh_motion_fe_values.quadrature_point(reference_point.q_point)[0]);
  
          const point_index_t quadrature_point_index = reference_point.mesh_motion_cell->user_index() + reference_point.q_point;
          material.set_state_parameters(quadrature_point_index, remapped_state_parameters, /*previous_deformation_gradient **/ /*mesh_motion_gradient*/unit_symmetric_tensor<dim+1, Number>());
          remapped_deformation_gradients[quadrature_point_index] = previous_deformation_gradient;
  
          ++next_to_process.at(group);
        }
  
      }
    }
  
    template<int dim, typename Number>
    void PlasticityLabProg<dim, Number>::remap_thermal_field(
        NewtonStepSystem &thermal_nonlinear_system,
        const DoFSystem<dim, Number> &thermal_dof_system,
        const NewtonStepSystem &mesh_motion_nonlinear_system,
        const DoFSystem<dim, Number> &mesh_motion_dof_system) {
  
      const Quadrature<dim> thermal_fe_support_point_quadrature(therm_fe.get_unit_support_points());
  
      const unsigned int n_q_points = thermal_fe_support_point_quadrature.size();
      [[maybe_unused]] const unsigned int dofs_per_cell = mesh_motion_fe.dofs_per_cell;
      const unsigned int thermal_dofs_per_cell = therm_fe.dofs_per_cell;
  
      FEValues<dim> mesh_motion_fe_values(
        mapping,
        mesh_motion_fe,
        thermal_fe_support_point_quadrature,
        update_values  | update_quadrature_points);
  
      FEValues<dim> thermal_fe_values(
        mapping,
        therm_fe,
        thermal_fe_support_point_quadrature,
        update_values | update_quadrature_points);
  
      std::vector< Number > previous_temperatures(1);
      std::vector< Tensor<1, dim, Number> > mesh_motion_value_increments(n_q_points);
      std::vector< Tensor<2, dim, Number> > mesh_motion_gradient_increments(n_q_points);
  
      const FEValuesExtractors::Vector displacements(0);
      const FEValuesExtractors::Scalar temperature (0);
  
      std::vector<ReferencePoint<dim, Number>> reference_points;
      std::vector<Point<dim, Number>> remapped_point_positions;
  
      auto cell = mesh_motion_dof_system.dof_handler.begin_active();
      auto thermal_cell = thermal_dof_system.dof_handler.begin_active();
      for (; cell != mesh_motion_dof_system.dof_handler.end(); ++cell, ++thermal_cell) {
        if (cell->is_locally_owned()) {
          mesh_motion_fe_values.reinit(cell);
  
          mesh_motion_fe_values[displacements].get_function_values(
            mesh_motion_nonlinear_system.current_increment,
            mesh_motion_value_increments);
  
          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
            const Point<dim, Number> reference_point_position = mesh_motion_fe_values.quadrature_point(q_point);
            const Point<dim, Number> remapped_point_position = reference_point_position - mesh_motion_value_increments[q_point];
  
            ReferencePoint<dim, Number> reference_point;
            reference_point.mesh_motion_cell = cell;
            reference_point.field_cell = thermal_cell;
            reference_point.q_point = q_point;
            reference_point.reference_point = reference_point_position;
            reference_point.remapped_point = remapped_point_position;
  
            reference_points.push_back(reference_point);
            remapped_point_positions.push_back(remapped_point_position);
  
          }
  
        }
      }
  
      MPI_Datatype PointType;
      MPI_Type_contiguous(dim, MPI_DOUBLE, &PointType);
      MPI_Type_commit(&PointType);
      std::vector<int> displs;
      std::vector< Point<dim, Number>> received_remapped_positions;
      std::vector<RemappedPoint<dim, Number>> thermal_points;
      int nprocesses, this_process;
      int num_reference_points = reference_points.size();
      MPI_Comm_size(mpi_communicator, &nprocesses);
      std::vector<int> reference_point_counts(nprocesses);
      MPI_Allgather(
        &num_reference_points,
        1, MPI_INT,
        &reference_point_counts[0],
        1, MPI_INT,
        mpi_communicator);
      MPI_Comm_rank(mpi_communicator, &this_process);
      displs.resize(nprocesses);
      displs[0] = 0;
      for (int i = 1; i < nprocesses; ++i) {
        displs[i] = displs[i - 1] + reference_point_counts[i - 1];
      }
  
      const unsigned int count_received_reference_points = displs[nprocesses - 1] + reference_point_counts[nprocesses - 1];
      received_remapped_positions.resize(count_received_reference_points);
      thermal_points.resize(count_received_reference_points);
      std::vector<char> this_process_owns_remapped_point(count_received_reference_points);
      std::vector<char> this_process_owns_thermal_point(count_received_reference_points);
  
      MPI_Allgatherv(
        &remapped_point_positions[0],
        num_reference_points,
        PointType,
        &received_remapped_positions[0],
        &reference_point_counts[0],
        &displs[0],
        PointType,
        mpi_communicator);
  
  
      for(unsigned int received_point_id=0; received_point_id < count_received_reference_points; received_point_id++) {
        auto point = received_remapped_positions[received_point_id];
        auto thermal_cell_and_point = GridTools::find_active_cell_around_point(
                            mapping,
                            thermal_dof_system.dof_handler,
                            point);
        auto thermal_cell = thermal_cell_and_point.first;
        auto thermal_unit_cell_point = thermal_cell_and_point.second;
        thermal_points[received_point_id].field_cell = thermal_cell;
        thermal_points[received_point_id].unit_cell_point = thermal_unit_cell_point;
        thermal_points[received_point_id].remapped_point = point;
        this_process_owns_thermal_point[received_point_id] = thermal_cell.state() == IteratorState::valid && thermal_cell->is_locally_owned()? 1 : 0;
      }
  
      std::vector<char> thermal_point_candidates(num_reference_points * nprocesses);
      for (int process = 0; process < nprocesses; ++process) {
        if (reference_point_counts[process] > 0)
          MPI_Gather(
            &this_process_owns_thermal_point[displs[process]],
            reference_point_counts[process],
            MPI_C_BOOL,
            &thermal_point_candidates[0],
            reference_point_counts[process],
            MPI_C_BOOL,
            process,
            mpi_communicator);
      }
  
      std::vector<unsigned int> thermal_point_owning_process(num_reference_points);
      for (int i = 0; i < num_reference_points; ++i) {
        thermal_point_owning_process[i] = 0;
        for (int j = 1; j < nprocesses; ++j) {
          if (thermal_point_candidates[j * num_reference_points + i] == 1) {
            thermal_point_owning_process[i] = j;
            continue;
          }
        }
      }
  
      std::vector<RemappedPoint<dim, Number>> mapping_thermal_points;
      std::vector<unsigned int> thermal_reference_point_owning_process;
      std::vector<unsigned int> thermal_reference_point_index_at_remote_process;
  

This should really be a vector<bool>, but addresses of individual elements of vector<bool> cannot be taken. It's a template specialization to save space

      std::vector<char> remote_thermal_point_is_accepted(num_reference_points * nprocesses);
      std::vector<char> local_thermal_point_is_accepted(count_received_reference_points);
      for (int i = 0; i < num_reference_points; ++i) {
        for (unsigned int j = 0; j < static_cast<unsigned int>(nprocesses); ++j) {
          remote_thermal_point_is_accepted[j * num_reference_points + i] = (thermal_point_owning_process[i] == j) ? 1 : 0;
        }
      }
  
      for (int process = 0; process < nprocesses; ++process) {
        if (reference_point_counts[process] > 0) {
          MPI_Scatter(
            &remote_thermal_point_is_accepted[0],
            reference_point_counts[process],
            MPI_C_BOOL,
            &local_thermal_point_is_accepted[displs[process]],
            reference_point_counts[process],
            MPI_C_BOOL,
            process,
            mpi_communicator);
        }
      }
  
      std::vector<unsigned int> thermal_remote_reference_point_counts(nprocesses, 0);
      for (int process = 0; process < nprocesses; ++process) {
        for (int i = 0; i < reference_point_counts[process]; ++i) {
          if (local_thermal_point_is_accepted[displs[process] + i]) {
            RemappedPoint<dim, Number> accepted_thermal_point = thermal_points[displs[process] + i];
  
            mapping_thermal_points.push_back(accepted_thermal_point);
            thermal_reference_point_owning_process.push_back(process);
            thermal_reference_point_index_at_remote_process.push_back(i);
            ++thermal_remote_reference_point_counts[process];
          }
        }
      }
  
      MPI_Barrier(mpi_communicator);
  
      std::vector<unsigned int> thermal_remote_remapped_point_counts(nprocesses);
      for (int i = 0; i < nprocesses; ++i) {
        MPI_Scatter(
          &thermal_remote_reference_point_counts[0],
          1, MPI_UNSIGNED,
          &thermal_remote_remapped_point_counts[i],
          1, MPI_UNSIGNED,
          i, mpi_communicator);
      }
  
  
  
      {
        std::vector<std::vector<Number> > local_previous_temperature_groups(nprocesses);
        std::vector<std::vector<Number> > remote_previous_temperature_groups(nprocesses);
  
        for (int i = 0; i < nprocesses; ++i) {
          local_previous_temperature_groups.at(i).resize(thermal_remote_reference_point_counts.at(i));
          remote_previous_temperature_groups.at(i).resize(thermal_remote_remapped_point_counts.at(i));
        }
  
        std::vector<unsigned int> next_to_process(nprocesses, 0);
        for (unsigned int i = 0; i < mapping_thermal_points.size(); ++i) {
          const unsigned int group = thermal_reference_point_owning_process.at(i);
          const RemappedPoint<dim, Number> thermal_point = mapping_thermal_points.at(i);
  
          Quadrature<dim> thermal_point_quadrature(
            std::vector<Point<dim, Number>> (1, thermal_point.unit_cell_point));
  
          FEValues<dim> thermal_point_fe_values(
            mapping,
            therm_fe,
            thermal_point_quadrature,
            update_values);
  
          thermal_point_fe_values.reinit(thermal_point.field_cell);
  
          thermal_point_fe_values[temperature].get_function_values(
            thermal_nonlinear_system.previous_deformation,
            previous_temperatures);
  
          local_previous_temperature_groups.at(group).at(next_to_process.at(group)) = previous_temperatures[0];
  
          ++next_to_process.at(group);
        }
  
  
  
        enum MessageFlag {
          PREVIOUS_TEMPERATURE
        };
  
        const unsigned int
          previous_temperature_requests_offset = 0,
          request_array_size = 1;
  
        std::vector<MPI_Request> requests_vector(2 * nprocesses * request_array_size);
  
        for (int i = 0; i < nprocesses; ++i) {
          MPI_Isend(
            local_previous_temperature_groups.at(i).data(),
            thermal_remote_reference_point_counts.at(i),
            MPI_DOUBLE, i, PREVIOUS_TEMPERATURE,
            mpi_communicator,
            &requests_vector[i + nprocesses * previous_temperature_requests_offset]);
        }
  
        for (int i = 0; i < nprocesses; ++i) {
          const unsigned int row_start = i + nprocesses * request_array_size;
  
          MPI_Irecv(
            remote_previous_temperature_groups.at(i).data(),
            thermal_remote_remapped_point_counts.at(i),
            MPI_DOUBLE, i, PREVIOUS_TEMPERATURE,
            mpi_communicator,
            &requests_vector[row_start + nprocesses * previous_temperature_requests_offset]);
  
        }
  
        std::vector<MPI_Status> statuses_vector(2 * request_array_size * nprocesses);
        MPI_Waitall(
          2 * request_array_size * nprocesses,
          &requests_vector[0],
          &statuses_vector[0]);
  
        TrilinosWrappers::SparsityPattern sparsity_pattern(
          thermal_dof_system.locally_owned_dofs,
          mpi_communicator);
  
        DoFTools::make_sparsity_pattern(
          thermal_dof_system.dof_handler, sparsity_pattern,
          AffineConstraints<Number>(),
          false,
          Utilities::MPI::this_mpi_process(mpi_communicator));
        sparsity_pattern.compress();
  
        TrilinosWrappers::SparseMatrix  projection_matrix(sparsity_pattern);
        TrilinosWrappers::MPI::Vector   projection_residual(thermal_dof_system.locally_owned_dofs, mpi_communicator);
  
        projection_matrix = 0;
        projection_residual = 0;
  
        projection_matrix = 0;
        projection_residual = 0;
        FullMatrix<Number> cell_matrix(thermal_dofs_per_cell, thermal_dofs_per_cell);
        Vector<Number> cell_residual(thermal_dofs_per_cell);
        next_to_process.clear();
        next_to_process.resize(nprocesses, 0);
        for (unsigned int i = 0; i < reference_points.size(); ++i) {
          unsigned int group = thermal_point_owning_process.at(i);
  
          cell_matrix = 0;
          cell_residual = 0;
  
          ReferencePoint<dim, Number> &reference_point = reference_points[i];
  
          thermal_fe_values.reinit(reference_point.field_cell);
  
          const Number remapped_previous_temperature = remote_previous_temperature_groups.at(group).at(next_to_process.at(group));
          for(unsigned int dof_i=0; dof_i<thermal_dofs_per_cell; dof_i++) {
            cell_residual(dof_i) +=
              remapped_previous_temperature
              * thermal_fe_values[temperature].value(dof_i, reference_point.q_point);
            for(unsigned int dof_j=0; dof_j<thermal_dofs_per_cell; dof_j++) {
              cell_matrix(dof_i, dof_j) +=
                thermal_fe_values[temperature].value(dof_i, reference_point.q_point)
                * thermal_fe_values[temperature].value(dof_j, reference_point.q_point);
            }
          }
  
          std::vector<types::global_dof_index> local_dof_indices(thermal_dofs_per_cell);
          reference_point.field_cell->get_dof_indices(local_dof_indices);
  
          projection_residual.add(local_dof_indices, cell_residual);
          for(unsigned int dof_i=0; dof_i<thermal_dofs_per_cell; dof_i++) {
            for(unsigned int dof_j=0; dof_j<thermal_dofs_per_cell; dof_j++) {
              projection_matrix.add(local_dof_indices[dof_i], local_dof_indices[dof_j], cell_matrix(dof_i, dof_j));
            }
          }
  
          ++next_to_process.at(group);
        }
  
        projection_matrix.compress(
          VectorOperation::add);
        projection_residual.compress(
          VectorOperation::add);
  

solve the thermal projection system

        TrilinosWrappers::PreconditionAMG preconditioner;
  
        const std::vector<std::vector<bool> > constant_modes
          = DoFTools::extract_constant_modes(thermal_dof_system.dof_handler,
                                             ComponentMask());
  
        TrilinosWrappers::PreconditionAMG::AdditionalData additional_data;
        additional_data.constant_modes = constant_modes;
        additional_data.elliptic = true;
        additional_data.n_cycles = 1;
        additional_data.w_cycle = false;
        additional_data.output_details = false;
        additional_data.smoother_sweeps = 2;
        additional_data.aggregation_threshold = 1e-2;
        preconditioner.initialize(projection_matrix, additional_data);
  
        TrilinosWrappers::MPI::Vector tmp(thermal_dof_system.locally_owned_dofs, mpi_communicator);
        const Number relative_accuracy = 1e-08;
        const Number solver_tolerance  = relative_accuracy
                                         * projection_matrix.residual(tmp, thermal_nonlinear_system.Newton_step_solution,
                                             projection_residual);
        SolverControl solver_control(projection_matrix.m(),
                                     solver_tolerance);
  
        SolverBicgstab<TrilinosWrappers::MPI::Vector> solver(solver_control);
  
        thermal_nonlinear_system.Newton_step_solution = 0;
  
        solver.solve(projection_matrix, thermal_nonlinear_system.Newton_step_solution,
                     projection_residual, preconditioner);
  
        thermal_nonlinear_system.previous_deformation = thermal_nonlinear_system.Newton_step_solution;
  
      }
  
    }
  
  
    template<int dim, typename Number>
    void PlasticityLabProg<dim, Number>::remap_mechanical_fields(
        NewtonStepSystem &mechanical_nonlinear_system,
        const DoFSystem<dim, Number> &mechanical_dof_system,
        const NewtonStepSystem &mesh_motion_nonlinear_system,
        const DoFSystem<dim, Number> &mesh_motion_dof_system) {
  
      const Quadrature<dim> mechanical_fe_support_point_quadrature(mech_fe.base_element(0).get_unit_support_points());
  
      const unsigned int n_q_points = mechanical_fe_support_point_quadrature.size();
      [[maybe_unused]] const unsigned int dofs_per_cell = mesh_motion_fe.dofs_per_cell;
      const unsigned int mechanical_dofs_per_cell = mech_fe.dofs_per_cell;
  
      FEValues<dim> mesh_motion_fe_values(
        mapping,
        mesh_motion_fe,
        mechanical_fe_support_point_quadrature,
        update_values  | update_quadrature_points);
  
      FEValues<dim> mechanical_fe_values(
        mapping,
        mech_fe,
        mechanical_fe_support_point_quadrature,
        update_values | update_quadrature_points);
  
      std::vector< Tensor<1, dim, Number> > previous_remapped_deformations(1);
      std::vector< Tensor<1, dim, Number> > previous_remapped_velocity(1);
      std::vector< Tensor<1, dim, Number> > previous_remapped_second_time_rate(1);
      std::vector< Number > previous_remapped_twist_deformations(1);
      std::vector< Number > previous_remapped_twist_velocity(1);
      std::vector< Number > previous_remapped_twist_second_time_rate(1);
      std::vector< Tensor<1, dim, Number> > mesh_motion_value_increments(n_q_points);
      std::vector< Tensor<2, dim, Number> > mesh_motion_gradient_increments(n_q_points);
  
      const FEValuesExtractors::Vector displacements(0);
      const FEValuesExtractors::Scalar angular_velocities(dim);
  
      std::vector<ReferencePoint<dim, Number>> reference_points;
      std::vector<Point<dim, Number>> remapped_point_positions;
      std::unordered_map<point_index_t, unsigned int> quadrature_point_reference_point_id;
      auto cell = mesh_motion_dof_system.dof_handler.begin_active();
      auto mechanical_cell = mechanical_dof_system.dof_handler.begin_active();
      for (; cell != mesh_motion_dof_system.dof_handler.end(); ++cell, ++mechanical_cell) {
        if (cell->is_locally_owned()) {
          mesh_motion_fe_values.reinit(cell);
  
          mesh_motion_fe_values[displacements].get_function_values(
            mesh_motion_nonlinear_system.current_increment,
            mesh_motion_value_increments);
  
          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
            const Point<dim, Number> reference_point_position = mesh_motion_fe_values.quadrature_point(q_point);
            const Point<dim, Number> remapped_point_position = reference_point_position - mesh_motion_value_increments[q_point];
  
            ReferencePoint<dim, Number> reference_point;
            reference_point.mesh_motion_cell = cell;
            reference_point.field_cell = mechanical_cell;
            reference_point.q_point = q_point;
            reference_point.reference_point = reference_point_position;
            reference_point.remapped_point = remapped_point_position;
  
            reference_points.push_back(reference_point);
            remapped_point_positions.push_back(remapped_point_position);
  
            const point_index_t quadrature_point_index = cell->user_index() + q_point;
            quadrature_point_reference_point_id[quadrature_point_index] = reference_points.size() - 1;
  
          }
  
        }
      }
  
      MPI_Datatype PointType;
      MPI_Type_contiguous(dim, MPI_DOUBLE, &PointType);
      MPI_Type_commit(&PointType);
      std::vector<int> displs;
      std::vector< Point<dim, Number>> received_remapped_positions;
      std::vector<RemappedPoint<dim, Number>> remapped_points;
      int nprocesses, this_process;
      int num_reference_points = reference_points.size();
      MPI_Comm_size(mpi_communicator, &nprocesses);
      std::vector<int> reference_point_counts(nprocesses);
      MPI_Allgather(
        &num_reference_points,
        1, MPI_INT,
        &reference_point_counts[0],
        1, MPI_INT,
        mpi_communicator);
      MPI_Comm_rank(mpi_communicator, &this_process);
      displs.resize(nprocesses);
      displs[0] = 0;
      for (int i = 1; i < nprocesses; ++i) {
        displs[i] = displs[i - 1] + reference_point_counts[i - 1];
      }
  
      const unsigned int count_received_reference_points = displs[nprocesses - 1] + reference_point_counts[nprocesses - 1];
      received_remapped_positions.resize(count_received_reference_points);
      remapped_points.resize(count_received_reference_points);
      std::vector<char> this_process_owns_remapped_point(count_received_reference_points);
  
      MPI_Allgatherv(
        &remapped_point_positions[0],
        num_reference_points,
        PointType,
        &received_remapped_positions[0],
        &reference_point_counts[0],
        &displs[0],
        PointType,
        mpi_communicator);
  
  
      for(unsigned int received_point_id=0; received_point_id < count_received_reference_points; received_point_id++) {
        auto point = received_remapped_positions[received_point_id];
        auto remapped_cell_and_point = GridTools::find_active_cell_around_point(
                            mapping,
                            mechanical_dof_system.dof_handler,
                            point);
        auto mechanical_cell = remapped_cell_and_point.first;
        auto mechanical_unit_cell_point = remapped_cell_and_point.second;
        remapped_points[received_point_id].field_cell = mechanical_cell;
        remapped_points[received_point_id].unit_cell_point = mechanical_unit_cell_point;
        remapped_points[received_point_id].remapped_point = point;
        this_process_owns_remapped_point[received_point_id] = mechanical_cell.state() == IteratorState::valid && mechanical_cell->is_locally_owned()? 1 : 0;
      }
  
      std::vector<char> remapped_point_candidates(num_reference_points * nprocesses);
      for (int process = 0; process < nprocesses; ++process) {
        if (reference_point_counts[process] > 0)
          MPI_Gather(
            &this_process_owns_remapped_point[displs[process]],
            reference_point_counts[process],
            MPI_C_BOOL,
            &remapped_point_candidates[0],
            reference_point_counts[process],
            MPI_C_BOOL,
            process,
            mpi_communicator);
      }
  
      std::vector<unsigned int> remapped_point_owning_process(num_reference_points);
      for (int i = 0; i < num_reference_points; ++i) {
        remapped_point_owning_process[i] = 0;
        for (int j = 1; j < nprocesses; ++j) {
          if (remapped_point_candidates[j * num_reference_points + i] == 1) {
            remapped_point_owning_process[i] = j;
            continue;
          }
        }
      }
  
      std::vector<RemappedPoint<dim, Number>> mapping_remapped_points;
      std::vector<unsigned int> reference_point_owning_process;
      std::vector<unsigned int> reference_point_index_at_remote_process;
  

This should really be a vector<bool>, but addresses of individual elements of vector<bool> cannot be taken. It's a template specialization to save space

      std::vector<char> remote_remapped_point_is_accepted(num_reference_points * nprocesses);
      std::vector<char> local_remapped_point_is_accepted(count_received_reference_points);
      for (int i = 0; i < num_reference_points; ++i) {
        for (unsigned int j = 0; j < static_cast<unsigned int>(nprocesses); ++j) {
          remote_remapped_point_is_accepted[j * num_reference_points + i] = (remapped_point_owning_process[i] == j) ? 1 : 0;
        }
      }
  
      for (int process = 0; process < nprocesses; ++process) {
        if (reference_point_counts[process] > 0) {
          MPI_Scatter(
            &remote_remapped_point_is_accepted[0],
            reference_point_counts[process],
            MPI_C_BOOL,
            &local_remapped_point_is_accepted[displs[process]],
            reference_point_counts[process],
            MPI_C_BOOL,
            process,
            mpi_communicator);
        }
      }
  
      std::vector<unsigned int> remote_reference_point_counts(nprocesses, 0);
      for (int process = 0; process < nprocesses; ++process) {
        for (int i = 0; i < reference_point_counts[process]; ++i) {
          if (local_remapped_point_is_accepted[displs[process] + i]) {
            RemappedPoint<dim, Number> accepted_remapped_point = remapped_points[displs[process] + i];
  
            mapping_remapped_points.push_back(accepted_remapped_point);
            reference_point_owning_process.push_back(process);
            reference_point_index_at_remote_process.push_back(i);
            ++remote_reference_point_counts[process];
          }
        }
      }
  
      MPI_Barrier(mpi_communicator);
  
      std::vector<unsigned int> remote_remapped_point_counts(nprocesses);
      for (int i = 0; i < nprocesses; ++i) {
        MPI_Scatter(
          &remote_reference_point_counts[0],
          1, MPI_UNSIGNED,
          &remote_remapped_point_counts[i],
          1, MPI_UNSIGNED,
          i, mpi_communicator);
      }
  
      std::vector<std::vector<Number> > local_previous_deformation_groups(nprocesses);
      std::vector<std::vector<Number> > local_previous_velocity_groups(nprocesses);
      std::vector<std::vector<Number> > local_previous_second_time_rate_groups(nprocesses);
      std::vector<std::vector<Number> > remote_previous_deformation_groups(nprocesses);
      std::vector<std::vector<Number> > remote_previous_velocity_groups(nprocesses);
      std::vector<std::vector<Number> > remote_previous_second_time_rate_groups(nprocesses);
  
      for (int i = 0; i < nprocesses; ++i) {
        local_previous_deformation_groups.at(i).resize((dim+1) * remote_reference_point_counts.at(i));
        local_previous_velocity_groups.at(i).resize((dim+1) * remote_reference_point_counts.at(i));
        local_previous_second_time_rate_groups.at(i).resize((dim+1) * remote_reference_point_counts.at(i));
        remote_previous_deformation_groups.at(i).resize((dim+1) * remote_remapped_point_counts.at(i));
        remote_previous_velocity_groups.at(i).resize((dim+1) * remote_remapped_point_counts.at(i));
        remote_previous_second_time_rate_groups.at(i).resize((dim+1) * remote_remapped_point_counts.at(i));
      }
  
      std::vector<unsigned int> next_to_process(nprocesses, 0);
      for (unsigned int i = 0; i < mapping_remapped_points.size(); ++i) {
        const unsigned int group = reference_point_owning_process.at(i);
        const RemappedPoint<dim, Number> remapped_point = mapping_remapped_points.at(i);
  
        Quadrature<dim> remapped_point_quadrature(
          std::vector<Point<dim, Number>> (1, remapped_point.unit_cell_point));
  
        FEValues<dim> remapped_point_fe_values(
          mapping,
          mech_fe,
          remapped_point_quadrature,
          update_values);
  
        remapped_point_fe_values.reinit(remapped_point.field_cell);
  
        remapped_point_fe_values[displacements].get_function_values(
          mechanical_nonlinear_system.previous_deformation,
          previous_remapped_deformations);
  
        remapped_point_fe_values[displacements].get_function_values(
          mechanical_nonlinear_system.previous_time_derivative,
          previous_remapped_velocity);
  
        remapped_point_fe_values[displacements].get_function_values(
          mechanical_nonlinear_system.previous_second_time_derivative,
          previous_remapped_second_time_rate);
  
        remapped_point_fe_values[angular_velocities].get_function_values(
          mechanical_nonlinear_system.previous_deformation,
          previous_remapped_twist_deformations);
  
        remapped_point_fe_values[angular_velocities].get_function_values(
          mechanical_nonlinear_system.previous_time_derivative,
          previous_remapped_twist_velocity);
  
        remapped_point_fe_values[angular_velocities].get_function_values(
          mechanical_nonlinear_system.previous_second_time_derivative,
          previous_remapped_twist_second_time_rate);
  
        for(unsigned int dim_i=0; dim_i<dim; dim_i++) {
          local_previous_deformation_groups.at(group).at(next_to_process.at(group) * (dim+1) + dim_i) =
              previous_remapped_deformations[0][dim_i];
          local_previous_velocity_groups.at(group).at(next_to_process.at(group) * (dim+1) + dim_i) =
              previous_remapped_velocity[0][dim_i];
          local_previous_second_time_rate_groups.at(group).at(next_to_process.at(group) * (dim+1) + dim_i) =
              previous_remapped_second_time_rate[0][dim_i];
        }
  
        local_previous_deformation_groups.at(group).at(next_to_process.at(group) * (dim+1) + dim) =
            previous_remapped_twist_deformations[0];
        local_previous_velocity_groups.at(group).at(next_to_process.at(group) * (dim+1) + dim) =
            previous_remapped_twist_velocity[0];
        local_previous_second_time_rate_groups.at(group).at(next_to_process.at(group) * (dim+1) + dim) =
            previous_remapped_twist_second_time_rate[0];
  
        ++next_to_process.at(group);
      }
  
      enum MessageFlag {
        PREVIOUS_DEFORMATION,
        PREVIOUS_VELOCITY,
        PREVIOUS_SECOND_TIME_RATE
      };
  
      const unsigned int
        previous_deformation_requests_offset = 0,
        previous_velocity_requests_offset = 1,
        previous_second_time_rate_requests_offset = 2,
        request_array_size = 3;
  
      std::vector<MPI_Request> requests_vector(2 * nprocesses * request_array_size);
  
      for (int i = 0; i < nprocesses; ++i) {
        MPI_Isend(
          local_previous_deformation_groups.at(i).data(),
          (dim+1) * remote_reference_point_counts.at(i),
          MPI_DOUBLE, i, PREVIOUS_DEFORMATION,
          mpi_communicator,
          &requests_vector[i + nprocesses * previous_deformation_requests_offset]);
  
        MPI_Isend(
          local_previous_velocity_groups.at(i).data(),
          (dim+1) * remote_reference_point_counts.at(i),
          MPI_DOUBLE, i, PREVIOUS_VELOCITY,
          mpi_communicator,
          &requests_vector[i + nprocesses * previous_velocity_requests_offset]);
  
        MPI_Isend(
          local_previous_second_time_rate_groups.at(i).data(),
          (dim+1) * remote_reference_point_counts.at(i),
          MPI_DOUBLE, i, PREVIOUS_SECOND_TIME_RATE,
          mpi_communicator,
          &requests_vector[i + nprocesses * previous_second_time_rate_requests_offset]);
      }
  
      for (int i = 0; i < nprocesses; ++i) {
        const unsigned int row_start = i + nprocesses * request_array_size;
  
        MPI_Irecv(
          remote_previous_deformation_groups.at(i).data(),
          (dim+1) * remote_remapped_point_counts.at(i),
          MPI_DOUBLE, i, PREVIOUS_DEFORMATION,
          mpi_communicator,
          &requests_vector[row_start + nprocesses * previous_deformation_requests_offset]);
  
        MPI_Irecv(
          remote_previous_velocity_groups.at(i).data(),
          (dim+1) * remote_remapped_point_counts.at(i),
          MPI_DOUBLE, i, PREVIOUS_VELOCITY,
          mpi_communicator,
          &requests_vector[row_start + nprocesses * previous_velocity_requests_offset]);
  
        MPI_Irecv(
          remote_previous_second_time_rate_groups.at(i).data(),
          (dim+1) * remote_remapped_point_counts.at(i),
          MPI_DOUBLE, i, PREVIOUS_SECOND_TIME_RATE,
          mpi_communicator,
          &requests_vector[row_start + nprocesses * previous_second_time_rate_requests_offset]);
  
      }
  
      std::vector<MPI_Status> statuses_vector(2 * request_array_size * nprocesses);
      MPI_Waitall(
        2 * request_array_size * nprocesses,
        &requests_vector[0],
        &statuses_vector[0]);
  
      TrilinosWrappers::SparsityPattern sparsity_pattern(
        mechanical_dof_system.locally_owned_dofs,
        mpi_communicator);
  
      DoFTools::make_sparsity_pattern(
        mechanical_dof_system.dof_handler, sparsity_pattern,
        AffineConstraints<Number>(),
        false,
        Utilities::MPI::this_mpi_process(mpi_communicator));
      sparsity_pattern.compress();
  
      TrilinosWrappers::SparseMatrix  projection_matrix(sparsity_pattern);
      TrilinosWrappers::MPI::Vector   projection_residual(mechanical_dof_system.locally_owned_dofs, mpi_communicator);
      TrilinosWrappers::MPI::Vector   projection_velocity_residual(mechanical_dof_system.locally_owned_dofs, mpi_communicator);
      TrilinosWrappers::MPI::Vector   projection_second_time_rate_residual(mechanical_dof_system.locally_owned_dofs, mpi_communicator);
      TrilinosWrappers::MPI::Vector   projection_solution(mechanical_dof_system.locally_owned_dofs, mpi_communicator);
  
      projection_matrix = 0;
      projection_residual = 0;
      projection_velocity_residual = 0;
      projection_second_time_rate_residual = 0;
  
      FullMatrix<Number> cell_matrix(mechanical_dofs_per_cell, mechanical_dofs_per_cell);
      Vector<Number> cell_residual(mechanical_dofs_per_cell);
      Vector<Number> cell_velocity_residual(mechanical_dofs_per_cell);
      Vector<Number> cell_second_time_rate_residual(mechanical_dofs_per_cell);
  
      next_to_process.clear();
      next_to_process.resize(nprocesses, 0);
      for (unsigned int i = 0; i < reference_points.size(); ++i) {
        unsigned int group = remapped_point_owning_process.at(i);
  
        cell_matrix = 0;
        cell_residual = 0;
        cell_velocity_residual = 0;
        cell_second_time_rate_residual = 0;
  
        ReferencePoint<dim, Number> &reference_point = reference_points[i];
  
        mechanical_fe_values.reinit(reference_point.field_cell);
  
        Tensor<1, dim, Number> remapped_previous_deformation;
        Tensor<1, dim, Number> remapped_previous_velocity;
        Tensor<1, dim, Number> remapped_previous_second_time_rate;
  
        for(unsigned int dim_i=0; dim_i<dim; dim_i++) {
          remapped_previous_deformation[dim_i] =
              reference_point.remapped_point[dim_i]
              - reference_point.reference_point[dim_i]
              + remote_previous_deformation_groups.at(group).at(next_to_process.at(group) * (dim+1) + dim_i);
          remapped_previous_velocity[dim_i] = remote_previous_velocity_groups.at(group).at(next_to_process.at(group) * (dim+1) + dim_i);
          remapped_previous_second_time_rate[dim_i] = remote_previous_second_time_rate_groups.at(group).at(next_to_process.at(group) * (dim+1) + dim_i);
        }
  
        const Number remapped_previous_twist_deformation = remote_previous_deformation_groups.at(group).at(next_to_process.at(group) * (dim+1) + dim);
        const Number remapped_previous_twist_velocity = remote_previous_velocity_groups.at(group).at(next_to_process.at(group) * (dim+1) + dim);
        const Number remapped_previous_twist_second_time_rate = remote_previous_second_time_rate_groups.at(group).at(next_to_process.at(group) * (dim+1) + dim);
  
        for(unsigned int dof_i=0; dof_i<mechanical_dofs_per_cell; dof_i++) {
          const auto shape_value_i = postprocess_tensor_dimension(
            mechanical_fe_values[displacements].value(dof_i, reference_point.q_point),
            mechanical_fe_values[angular_velocities].value(dof_i, reference_point.q_point));
  
          cell_residual(dof_i) +=
            shape_value_i * postprocess_tensor_dimension(remapped_previous_deformation, remapped_previous_twist_deformation);
  
          cell_velocity_residual(dof_i) +=
            shape_value_i * postprocess_tensor_dimension(remapped_previous_velocity, remapped_previous_twist_velocity);
  
          cell_second_time_rate_residual(dof_i) +=
            shape_value_i * postprocess_tensor_dimension(remapped_previous_second_time_rate, remapped_previous_twist_second_time_rate);
  
          for(unsigned int dof_j=0; dof_j<mechanical_dofs_per_cell; dof_j++) {
            const auto shape_value_j = postprocess_tensor_dimension(
              mechanical_fe_values[displacements].value(dof_j, reference_point.q_point),
              mechanical_fe_values[angular_velocities].value(dof_j, reference_point.q_point));
  
            cell_matrix(dof_i, dof_j) += shape_value_i * shape_value_j;
  
          }
        }
  
        std::vector<types::global_dof_index> local_dof_indices(mechanical_dofs_per_cell);
        reference_point.field_cell->get_dof_indices(local_dof_indices);
        projection_residual.add(local_dof_indices, cell_residual);
        projection_velocity_residual.add(local_dof_indices, cell_velocity_residual);
        projection_second_time_rate_residual.add(local_dof_indices, cell_second_time_rate_residual);
        for(unsigned int dof_i=0; dof_i<mechanical_dofs_per_cell; dof_i++) {
          for(unsigned int dof_j=0; dof_j<mechanical_dofs_per_cell; dof_j++) {
            projection_matrix.add(local_dof_indices[dof_i], local_dof_indices[dof_j], cell_matrix(dof_i, dof_j));
          }
        }
  
        ++next_to_process.at(group);
      }
  
      projection_matrix.compress(VectorOperation::add);
      projection_residual.compress(VectorOperation::add);
      projection_velocity_residual.compress(VectorOperation::add);
      projection_second_time_rate_residual.compress(VectorOperation::add);
  

solve the projection system

      TrilinosWrappers::PreconditionAMG preconditioner;
  
      const std::vector<std::vector<bool> > constant_modes
        = DoFTools::extract_constant_modes(mechanical_dof_system.dof_handler,
                                           ComponentMask());
  
      TrilinosWrappers::PreconditionAMG::AdditionalData additional_data;
      additional_data.constant_modes = constant_modes;
      additional_data.elliptic = true;
      additional_data.n_cycles = 1;
      additional_data.w_cycle = false;
      additional_data.output_details = false;
      additional_data.smoother_sweeps = 2;
      additional_data.aggregation_threshold = 1e-2;
      preconditioner.initialize(projection_matrix, additional_data);
  
      TrilinosWrappers::MPI::Vector tmp(mechanical_dof_system.locally_owned_dofs, mpi_communicator);
      const Number relative_accuracy = 1e-08;
      const Number solver_tolerance  = relative_accuracy
                                       * projection_matrix.residual(tmp, projection_solution,
                                           projection_residual);
      SolverControl solver_control(projection_matrix.m(),
                                   solver_tolerance);
  
      SolverBicgstab<TrilinosWrappers::MPI::Vector> solver(solver_control);
  
      projection_solution = 0;
  
      solver.solve(projection_matrix, projection_solution,
                   projection_residual, preconditioner);
  
      mechanical_nonlinear_system.previous_deformation = projection_solution;
  
      projection_solution = 0;
      solver.solve(projection_matrix, projection_solution,
                   projection_velocity_residual, preconditioner);
      mechanical_nonlinear_system.previous_time_derivative = projection_solution;
  
      projection_solution = 0;
      solver.solve(projection_matrix, projection_solution,
                   projection_second_time_rate_residual, preconditioner);
      mechanical_nonlinear_system.previous_second_time_derivative = projection_solution;
  
    }
  
  
    template <int dim, typename Number>
    void PlasticityLabProg<dim, Number>::assemble_mechanical_system(
      NewtonStepSystem &Newton_system,
      const DoFSystem<dim, Number> &mechanical_dof_system,
      const LBCSystem<dim, Number, dim+1> &mechanical_lbc_system,
      const NewtonStepSystem &mesh_motion_nonlinear_system,
      const DoFSystem<dim, Number> &mesh_motion_dof_system,
      const NewtonStepSystem &thermal_Newton_system,
      const DoFSystem<dim, Number> &thermal_dof_system,
      const DoFSystem<dim, Number> &mixed_fe_dof_system,
      Material<dim+1, Number> &material,
      const std::vector< MixedFEProjector<dim, Number> > &mixed_fe_projector,
      const bool fill_system_matrix,
      const bool update_material_state) {
      FEValues<dim> fe_values(
        mapping,
        mech_fe,
        quadrature_formula,
        update_values | update_quadrature_points | update_gradients | update_JxW_values);
  
      FEFaceValues<dim> fe_face_values(
        mapping,
        mech_fe,
        face_quadrature_formula,
        update_values | update_quadrature_points | update_normal_vectors | update_JxW_values);
  
      FEValues<dim> fe_therm_values(
        mapping,
        therm_fe,
        quadrature_formula,
        update_values);
  
      FEValues<dim> mesh_motion_fe_values(
        mapping,
        mesh_motion_fe,
        quadrature_formula,
        update_values | update_gradients);
  
      FEValues<dim> mixed_fe_values(
        mapping,
        mixed_var_fe,
        quadrature_formula,
        update_values);
  
      const unsigned int dofs_per_cell = mech_fe.dofs_per_cell;
      const unsigned int n_q_points = quadrature_formula.size();
      const unsigned int mixed_dofs_per_cell = mixed_var_fe.dofs_per_cell;
      const unsigned int n_face_q_points = face_quadrature_formula.size();
  
      std::vector< Tensor<2, dim, Number> > current_displacement_gradients(n_q_points);
      std::vector< Tensor<2, dim, Number> > displacement_gradient_increments(n_q_points);
  
      std::vector< Tensor<1, dim, Number> > current_displacement_values(n_q_points);
      std::vector< Tensor<1, dim, Number> > displacement_value_increments(n_q_points);
  
      std::vector< Tensor<2, dim, Number> > displacement_gradient_previous_time_rates(n_q_points);
      std::vector< Tensor<2, dim, Number> > displacement_gradient_previous_second_time_rates(n_q_points);
  
      std::vector< Tensor<1, dim, Number> > displacement_increments(n_q_points);
      std::vector< Tensor<1, dim, Number> > displacement_previous_time_rates(n_q_points);
      std::vector< Tensor<1, dim, Number> > displacement_previous_second_time_rates(n_q_points);
  
      std::vector< Tensor<1, dim, Number> > angular_velocity_gradient_increments(n_q_points);
      std::vector< Tensor<1, dim, Number> > angular_velocity_gradient_previous_time_rates(n_q_points);
      std::vector< Tensor<1, dim, Number> > angular_velocity_gradient_previous_second_time_rates(n_q_points);
  
      std::vector< Tensor<1, dim, Number> > mesh_motion_value_increments(n_q_points);
      std::vector< Tensor<2, dim, Number> > mesh_motion_gradient_increments(n_q_points);
  
      std::vector< Number > angular_velocity_increments(n_q_points);
      std::vector< Number > angular_velocity_previous_time_rates(n_q_points);
      std::vector< Number > angular_velocity_previous_second_time_rates(n_q_points);
  
      std::vector< Number > current_temperature_values(n_q_points);
      std::vector< Number > updated_temperature_increments(n_q_points);
      std::vector< Number > updated_temperature_values(n_q_points);
  
      std::vector< Tensor<1, dim, Number> > face_displacement_value_increments(n_face_q_points);
  
      std::vector< Number > deformation_jacobians(n_q_points);
      std::vector< Number > previous_deformation_jacobian(n_q_points);
      std::vector< std::vector<Number> > strain_divergences(
        dofs_per_cell,
        std::vector<Number>(n_q_points));
      std::vector< std::vector<Number> > jacobian_tangents(
        dofs_per_cell,
        std::vector<Number>(n_q_points));
  
      std::vector< std::vector< std::vector < Number> > > strain_divergence_tangents(
        dofs_per_cell,
        std::vector< std::vector< Number> >(
          dofs_per_cell,
          std::vector<Number>(n_q_points)));
  
      std::vector<Number> projected_temperature_coefficients(mixed_dofs_per_cell);
      std::vector<Number> projected_Jacobian_coefficients(mixed_dofs_per_cell);
      std::vector<Number> projected_previous_temperature_coefficients(mixed_dofs_per_cell);
      std::vector<Number> projected_previous_Jacobian_coefficients(mixed_dofs_per_cell);
      std::vector<std::vector<Number> > projected_strain_divergence_coefficients(
        dofs_per_cell,
        std::vector<Number>(mixed_dofs_per_cell));
      std::vector<std::vector<Number> > projected_jacobian_tangent_coefficients(
        dofs_per_cell,
        std::vector<Number>(mixed_dofs_per_cell));
      std::vector<std::vector<std::vector<Number> > > projected_strain_divergence_tangent_coefficients(
        dofs_per_cell, std::vector< std::vector< Number> >(
          dofs_per_cell,
          std::vector<Number>(mixed_dofs_per_cell)));
  
      std::vector<Number> projected_strain_divergence(dofs_per_cell);
      std::vector<Number> projected_jacobian_tangent(dofs_per_cell);
      std::vector<std::vector<Number> > projected_strain_divergence_tangent(
        dofs_per_cell,
        std::vector<Number>(dofs_per_cell));
  
      FullMatrix<Number> cell_matrix(dofs_per_cell, dofs_per_cell);
      Vector<Number> cell_residual(dofs_per_cell);
  
      Vector<Number> mixed_values (mixed_dofs_per_cell);
  
      const FEValuesExtractors::Vector displacements (0);
      const FEValuesExtractors::Scalar angular_velocity(dim);
      const FEValuesExtractors::Scalar temperature (0);
  
      Newton_system.Newton_step_matrix = 0;
      Newton_system.Newton_step_residual = 0;
  
      double alpha_m, alpha_f, gamma, beta;
      get_generalized_alpha_method_params(
          &alpha_m, &alpha_f, &gamma, &beta, rho_infty);
  
      const Number d_second_time_rate_d_increment = (1./(beta*time_increment*time_increment));
      const Number d_time_rate_d_increment = gamma/(beta*time_increment);
  
      bool kinematic_domains_are_valid = true;  // innocent until proven guilty
  
      auto cell = mechanical_dof_system.dof_handler.begin_active();
      auto endc = mechanical_dof_system.dof_handler.end();
      auto thermal_cell = thermal_dof_system.dof_handler.begin_active();
      auto mesh_motion_cell = mesh_motion_dof_system.dof_handler.begin_active();
      auto mixed_fe_cell = mixed_fe_dof_system.dof_handler.begin_active();
      for (; cell != endc; ++cell, ++thermal_cell, ++mesh_motion_cell, ++mixed_fe_cell) {
        if (cell->is_locally_owned()) {
          cell_matrix = 0;
          cell_residual = 0;
  
          fe_values.reinit (cell);
          fe_therm_values.reinit (thermal_cell);
          mixed_fe_values.reinit (mixed_fe_cell);
          mesh_motion_fe_values.reinit (mesh_motion_cell);
  
          fe_values[displacements].get_function_gradients(
            Newton_system.current_increment,
            displacement_gradient_increments);
  
          fe_values[displacements].get_function_gradients(
            Newton_system.previous_deformation,
            current_displacement_gradients);
  
          fe_values[displacements].get_function_values(
            Newton_system.current_increment,
            displacement_value_increments);
  
          fe_values[displacements].get_function_values(
            Newton_system.previous_deformation,
            current_displacement_values);
  
          fe_values[displacements].get_function_gradients(
            Newton_system.previous_time_derivative,
            displacement_gradient_previous_time_rates);
  
          fe_values[displacements].get_function_gradients(
            Newton_system.previous_second_time_derivative,
            displacement_gradient_previous_second_time_rates);
  
          fe_values[displacements].get_function_values(
            Newton_system.current_increment,
            displacement_increments);
  
          fe_values[displacements].get_function_values(
            Newton_system.previous_time_derivative,
            displacement_previous_time_rates);
  
          fe_values[displacements].get_function_values(
            Newton_system.previous_second_time_derivative,
            displacement_previous_second_time_rates);
  

Angular velocity

          fe_values[angular_velocity].get_function_gradients(
            Newton_system.current_increment,
            angular_velocity_gradient_increments);
  
          fe_values[angular_velocity].get_function_gradients(
            Newton_system.previous_time_derivative,
            angular_velocity_gradient_previous_time_rates);
  
          fe_values[angular_velocity].get_function_gradients(
            Newton_system.previous_second_time_derivative,
            angular_velocity_gradient_previous_second_time_rates);
  
          fe_values[angular_velocity].get_function_values(
            Newton_system.current_increment,
            angular_velocity_increments);
  
          fe_values[angular_velocity].get_function_values(
            Newton_system.previous_time_derivative,
            angular_velocity_previous_time_rates);
  
          fe_values[angular_velocity].get_function_values(
            Newton_system.previous_second_time_derivative,
            angular_velocity_previous_second_time_rates);
  

mesh motion

          mesh_motion_fe_values[displacements].get_function_gradients(
            mesh_motion_nonlinear_system.current_increment,
            mesh_motion_gradient_increments);
  
          mesh_motion_fe_values[displacements].get_function_gradients(
            mesh_motion_nonlinear_system.current_increment,
            mesh_motion_gradient_increments);
  

temperature

          fe_therm_values[temperature].get_function_values (
            thermal_Newton_system.previous_deformation,
            current_temperature_values);
          fe_therm_values[temperature].get_function_values (
            thermal_Newton_system.current_increment,
            updated_temperature_increments);
  

get vectors for projection onto mixed fe values

          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
            const point_index_t quadrature_point_index = cell->user_index() + q_point;
            updated_temperature_values.at(q_point) =
              current_temperature_values.at(q_point)
              + updated_temperature_increments.at(q_point);
  
            const auto current_F = get_deformation_gradient(
              current_displacement_gradients[q_point],
              current_displacement_values[q_point][0]/fe_values.quadrature_point(q_point)[0]
              );
            const Number material_Jacobian = material.get_material_Jacobian(quadrature_point_index) / determinant(current_F);
  
            [[maybe_unused]] const auto mesh_motion_gradient = get_deformation_gradient(
              -mesh_motion_gradient_increments[q_point],
              -mesh_motion_value_increments[q_point][0]/fe_values.quadrature_point(q_point)[0]);
  
            const Tensor<2, dim+1, Number> updated_F = get_deformation_gradient(
              current_displacement_gradients[q_point] + displacement_gradient_increments[q_point],
              (current_displacement_values[q_point][0] + displacement_value_increments[q_point][0])
                /fe_values.quadrature_point(q_point)[0]
              );
  
  
            const Number Jacobian = material_Jacobian * determinant(updated_F);
            deformation_jacobians.at(q_point) = Jacobian;
            previous_deformation_jacobian.at(q_point) = material_Jacobian * determinant(current_F);
  
            const auto inv_updated_F = invert(updated_F);
            std::vector<Tensor<2, dim+1, Number>> rate_gradients(dofs_per_cell);
            for (unsigned int i = 0; i < dofs_per_cell; ++i) {
              rate_gradients[i] = postprocess_tensor_dimension(
                fe_values[displacements].gradient(i, q_point),
                fe_values[displacements].value(i, q_point)[0]/fe_values.quadrature_point(q_point)[0]) * inv_updated_F;
            }
  
            for (unsigned int i = 0; i < dofs_per_cell; ++i) {
              const Number strain_divergence_i = trace(rate_gradients[i]);
              strain_divergences[i].at(q_point) = Jacobian * strain_divergence_i;
              if (fill_system_matrix) {
                jacobian_tangents[i].at(q_point) = Jacobian * strain_divergence_i;
                for (unsigned int j = 0; j < dofs_per_cell; ++j) {
                  strain_divergence_tangents[i][j].at(q_point) =
                    Jacobian
                    * (trace(rate_gradients[i]) * trace(rate_gradients[j])
                        - trace(rate_gradients[i] * rate_gradients[j]));
                }
              }
            }
          }
  
          const unsigned int cell_index = cell->user_index() / n_q_points;
  
          mixed_fe_projector[cell_index].project(
            &projected_temperature_coefficients,
            updated_temperature_values);
          mixed_fe_projector[cell_index].project(
            &projected_Jacobian_coefficients,
            deformation_jacobians);
          mixed_fe_projector[cell_index].project(
            &projected_previous_Jacobian_coefficients,
            previous_deformation_jacobian);
          mixed_fe_projector[cell_index].project(
            &projected_previous_temperature_coefficients,
            current_temperature_values);
          for (unsigned int i = 0; i < dofs_per_cell; ++i) {
            mixed_fe_projector[cell_index].project(
              &projected_strain_divergence_coefficients[i],
              strain_divergences[i]);
            if (fill_system_matrix) {
              mixed_fe_projector[cell_index].project(
                &projected_jacobian_tangent_coefficients[i],
                jacobian_tangents[i]);
              for (unsigned int j = 0; j < dofs_per_cell; ++j) {
                mixed_fe_projector[cell_index].project(
                  &projected_strain_divergence_tangent_coefficients[i][j],
                  strain_divergence_tangents[i][j]);
              }
            }
          }
  
          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
            ConstitutiveModelUpdateFlags materialUpdateFlags =
              (update_pressure | update_stress_deviator);
  
            if (fill_system_matrix) {
              materialUpdateFlags |=
                (update_pressure_tangent | update_stress_deviator_tangent);
            }
  
            ConstitutiveModelRequest<dim+1, Number> previous_constitutive_request(materialUpdateFlags);
  
            if (update_material_state) {
              materialUpdateFlags |= update_material_point_history;
            }
  
            const point_index_t quadrature_point_index = cell->user_index() + q_point;
            ConstitutiveModelRequest<dim+1, Number> constitutive_request(materialUpdateFlags);
  
            for (unsigned int i = 0; i < mixed_dofs_per_cell; ++i) {
              mixed_values(i) = mixed_fe_values.shape_value(i, q_point);
            }
  
            const Number radius = fe_values.quadrature_point(q_point)[0];
            const auto current_F = get_deformation_gradient(
              current_displacement_gradients[q_point],
              current_displacement_values[q_point][0]/radius
              );
  
            const auto mesh_motion_gradient = get_deformation_gradient(
              -mesh_motion_gradient_increments[q_point],
              -mesh_motion_value_increments[q_point][0]/radius);
  
            const Tensor<2, dim+1, Number> updated_F = get_deformation_gradient(
              current_displacement_gradients[q_point] + displacement_gradient_increments[q_point],
              (current_displacement_values[q_point][0] + displacement_value_increments[q_point][0])/radius);
  
            const Number material_Jacobian = material.get_material_Jacobian(quadrature_point_index) / determinant(current_F);
            [[maybe_unused]] const Number mesh_motion_Jacobian = determinant(mesh_motion_gradient);
  
            const auto inv_updated_F = invert(updated_F);
            const Number Jacobian = determinant(updated_F) * material_Jacobian;
            const Number previous_Jacobian = determinant(current_F) * material_Jacobian;
  
  
            const Number DENSITY = 8.96e-9;
  
            const Tensor<1, dim+1, Number> displacement_increment = postprocess_tensor_dimension(displacement_increments[q_point]);
            const Tensor<1, dim+1, Number> displacement_previous_time_rate = postprocess_tensor_dimension(displacement_previous_time_rates[q_point]);
            const Tensor<1, dim+1, Number> displacement_previous_second_time_rate = postprocess_tensor_dimension(displacement_previous_second_time_rates[q_point]);
            const Tensor<2, dim+1, Number> displacement_gradient_previous_time_rate = postprocess_tensor_dimension(
              displacement_gradient_previous_time_rates[q_point], displacement_previous_time_rates[q_point][0]/radius);
            const Tensor<2, dim+1, Number> displacement_gradient_previous_second_time_rate = postprocess_tensor_dimension(
              displacement_gradient_previous_second_time_rates[q_point], displacement_previous_second_time_rates[q_point][0]/radius);
  
            const Tensor<1, dim+1, Number> uc_increment = scalar_to_angular_tensor(angular_velocity_increments[q_point]);
            const Tensor<1, dim+1, Number> vc_n = scalar_to_angular_tensor(angular_velocity_previous_time_rates[q_point]);
            const Tensor<1, dim+1, Number> d_vc_d_t_n = scalar_to_angular_tensor(angular_velocity_previous_second_time_rates[q_point]);
  
            const Tensor<1, dim+1, Number> d2_x_dt_2_n_plus_1 =
              (1./(beta*time_increment*time_increment))
              * (displacement_increment
                 - time_increment * displacement_previous_time_rate
                 - time_increment * time_increment * (0.5-beta) * displacement_previous_second_time_rate);
  
            const Tensor<1, dim+1, Number> d_x_dt_n_plus_1 =
              displacement_previous_time_rate + time_increment * ((1-gamma) * displacement_previous_second_time_rate + gamma * d2_x_dt_2_n_plus_1);
  
            const Tensor<2, dim+1, Number> Grad_d_2_x_d_t_2_n_plus_1 =
              (1./(beta*time_increment*time_increment))
              * (postprocess_tensor_dimension(
                  displacement_gradient_increments[q_point],
                  displacement_value_increments[q_point][0]/radius)
                 - time_increment * displacement_gradient_previous_time_rate
                 - time_increment * time_increment * (0.5-beta) * displacement_gradient_previous_second_time_rate);
  
            const Tensor<2, dim+1, Number> Grad_d_x_d_t_n_plus_1 =
              displacement_gradient_previous_time_rate
              + time_increment * (
                  (1.-gamma) * displacement_gradient_previous_second_time_rate
                   + gamma * Grad_d_2_x_d_t_2_n_plus_1);
  
            const Tensor<1, dim+1, Number> d_vc_d_t_n_plus_1 =
              (1./(beta*time_increment*time_increment))
              * (uc_increment
                 - time_increment * vc_n
                 - time_increment * time_increment * (0.5-beta) * d_vc_d_t_n);
  
            const Tensor<1, dim+1, Number> vc_n_plus_1 =
              vc_n + time_increment * ((1-gamma) * d_vc_d_t_n + gamma * d_vc_d_t_n_plus_1);
  
            const Number thR_increment = angular_velocity_increments[q_point];
            const Number d_thR_d_t_n = angular_velocity_previous_time_rates[q_point];
            const Number d2_thR_d_t2_n = angular_velocity_previous_second_time_rates[q_point];
  
            const Number d2_thR_d_t2_n_plus_1 =
              (1./(beta*time_increment*time_increment))
              * (thR_increment
                 - time_increment * d_thR_d_t_n
                 - time_increment * time_increment * (0.5-beta) * d2_thR_d_t2_n);
  
            const Number d_thR_d_t_n_plus_1 =
              d_thR_d_t_n + time_increment * ((1-gamma) * d2_thR_d_t2_n + gamma * d2_thR_d_t2_n_plus_1);
  
            const Number one_plus_r_over_R_n = 1.0 + current_displacement_values[q_point][0]/radius;
            const Number one_plus_r_over_R_n_plus_1 = 1.0 + (current_displacement_values[q_point][0] + displacement_value_increments[q_point][0])/radius;
            const Number d_r_d_t_n_over_R = displacement_previous_time_rate[0] / radius;
            const Number d_r_d_t_n_plus_one_over_R = d_x_dt_n_plus_1[0] / radius;
  
            Tensor<1, dim+1, Number> e_hat_R;
            e_hat_R[0] = 1;
  
            const Tensor<1, dim+1, Number> acceleration_n_plus_1 =
              d2_x_dt_2_n_plus_1
              + scalar_to_angular_tensor(
                  d_r_d_t_n_plus_one_over_R * d_thR_d_t_n_plus_1
                  + one_plus_r_over_R_n_plus_1 * d2_thR_d_t2_n_plus_1)
              - (1.0/radius) * one_plus_r_over_R_n_plus_1 * std::pow(d_thR_d_t_n_plus_1, 2) * e_hat_R;
  
            const Tensor<1, dim+1, Number> acceleration_n =
              displacement_previous_second_time_rate
              + scalar_to_angular_tensor(
                  d_r_d_t_n_over_R * d_thR_d_t_n
                  + one_plus_r_over_R_n * d2_thR_d_t2_n)
              - (1.0/radius) * one_plus_r_over_R_n * std::pow(d_thR_d_t_n, 2) * e_hat_R;
  
            const Tensor<1, dim+1, Number> acceleration_n_plus_1_minus_alpha_m =
              alpha_m * acceleration_n + (1-alpha_m) * acceleration_n_plus_1;
  
            Tensor<2, dim+1, Number> acceleration_n_plus_1_minus_alpha_m_tangent_modulus;
            for(unsigned int i=0; i<dim; i++) {
              acceleration_n_plus_1_minus_alpha_m_tangent_modulus[i][i] += (1-alpha_m) * d_second_time_rate_d_increment;
            }
            acceleration_n_plus_1_minus_alpha_m_tangent_modulus[dim][0] +=
              (1-alpha_m) * (d_time_rate_d_increment/radius * d_thR_d_t_n_plus_1 + 1.0/radius * d2_thR_d_t2_n_plus_1);
            acceleration_n_plus_1_minus_alpha_m_tangent_modulus[0][0] +=
              (1-alpha_m) * (-(1.0/radius) * (1.0/radius) * std::pow(d_thR_d_t_n_plus_1, 2));
  
            acceleration_n_plus_1_minus_alpha_m_tangent_modulus[dim][dim] +=
              (1-alpha_m)
              * (d_r_d_t_n_plus_one_over_R * d_time_rate_d_increment
                 + one_plus_r_over_R_n_plus_1 * d_second_time_rate_d_increment);
            acceleration_n_plus_1_minus_alpha_m_tangent_modulus[0][dim] +=
              (1-alpha_m)
              * (-(1.0/radius) * one_plus_r_over_R_n * 2 * d_thR_d_t_n * d_time_rate_d_increment);
  
  
            [[maybe_unused]] const auto d2_x_dt_2_n_plus_1_minus_alpha_m =
              alpha_m * displacement_previous_second_time_rate + (1-alpha_m)*d2_x_dt_2_n_plus_1;
            [[maybe_unused]] const Number d2_x_dt_2_tangent_1_minus_alpha_m = (1-alpha_m)*d_second_time_rate_d_increment;
  
            [[maybe_unused]] const auto d_Fc_dt_vc_n_plus_1_minus_alpha_f =
                          alpha_m * displacement_gradient_previous_time_rate * vc_n
                          + (1-alpha_m) * Grad_d_x_d_t_n_plus_1 * vc_n_plus_1;
  
            [[maybe_unused]] const auto d_Fc_dt_vc_n_plus_1_minus_alpha_f_F_tangent = (1-alpha_m) * d_time_rate_d_increment * vc_n_plus_1;
            [[maybe_unused]] const auto d_Fc_dt_vc_n_plus_1_minus_alpha_f_v_tangent = (1-alpha_m) * Grad_d_x_d_t_n_plus_1;
  
            [[maybe_unused]] const auto F_c_d_vc_d_t_n_plus_1_minus_alpha_m =
              alpha_m * current_F * d_vc_d_t_n
              + (1-alpha_m) * updated_F * d_vc_d_t_n_plus_1;
  
            [[maybe_unused]] const auto F_c_d_vc_d_t_n_plus_1_minus_alpha_m_F_tangent = (1-alpha_m) * d_vc_d_t_n_plus_1;
            [[maybe_unused]] const auto F_c_d_vc_d_t_n_plus_1_minus_alpha_m_V_tangent = (1-alpha_m) * updated_F * d_second_time_rate_d_increment;
  
            std::vector<Tensor<2, dim+1, Number>> rate_gradients(dofs_per_cell);
            std::vector<Tensor<2, dim+1, Number>> angular_rate_gradients(dofs_per_cell);
            for (unsigned int i = 0; i < dofs_per_cell; ++i) {
              rate_gradients[i] =
                  postprocess_tensor_dimension(
                      fe_values[displacements].gradient(i, q_point),
                      fe_values[displacements].value(i, q_point)[0]/radius) * inv_updated_F;
              angular_rate_gradients[i] =
                  order_1_tensor_to_angular_gradient(
                      fe_values[angular_velocity].gradient(i, q_point),
                      -fe_values[angular_velocity].value(i, q_point)/radius);
            }
  
            const Tensor<2, dim+1, Number> d_X_prime_d_X = deformation_gradient_from_angular_displacement_gradient(
                      -angular_velocity_increments[q_point],
                      -angular_velocity_gradient_increments[q_point],
                      radius
                    );
  
            const Tensor<2, dim+1, Number> rotation_to_X_prime_frame = rotation_tensor_to_transform_B_e(
                      -angular_velocity_increments[q_point],
                      radius
                    );
  
            const Tensor<2, dim+1, Number> inv_d_X_prime_d_X = invert(d_X_prime_d_X);
  
            const auto f_m_n_plus_1 = inv_d_X_prime_d_X * rotation_to_X_prime_frame;
  

std::cout << "f_r: " << inv_d_X_prime_d_X << std::endl; std::cout << "R: " << previous_elastic_deformation_transformation_tensor << std::endl; std::cout << "f_m_n+1: " << f_m_n_plus_1 << std::endl;

            Number projected_jacobian = 0;
            Number projected_previous_jacobian = 0;
            Number projected_temperature = 0;
            Number projected_previous_temperature = 0;
            for (unsigned int i = 0; i < mixed_dofs_per_cell; ++i) {
              projected_jacobian +=
                mixed_values(i) * projected_Jacobian_coefficients.at(i);
              projected_previous_jacobian +=
                mixed_values(i) * projected_previous_Jacobian_coefficients.at(i);
              projected_temperature +=
                mixed_values(i) * projected_temperature_coefficients.at(i);
              projected_previous_temperature +=
                mixed_values(i) * projected_previous_temperature_coefficients.at(i);
            }
  
            const auto unnormalized_deformation_gradient_increment = updated_F * f_m_n_plus_1 * invert(current_F);
  
            const auto deformation_gradient_increment =
              std::pow(determinant(unnormalized_deformation_gradient_increment),
                       -Constants<dim, Number>::one_third()) * unnormalized_deformation_gradient_increment;
  
            if(std::abs(determinant(deformation_gradient_increment)-1.0) > 1e-8)
              std::cout << "determinant(deformation_gradient_increment): " << determinant(deformation_gradient_increment) << std::endl;
  
            if(false && std::isnan(deformation_gradient_increment.norm())) {
              std::cout << "deformation_gradient_increment is nan: " << deformation_gradient_increment << std::endl;
              std::cout << "Jacobian: " << Jacobian
                        << "\nf_m_n_plus_1: " << f_m_n_plus_1
                        << "\ndet(f_m_n_plus_1): " << determinant(f_m_n_plus_1)
                        << "\nprevious_Jacobian: " << previous_Jacobian
                        << "\nupdated_F: " << updated_F
                        << "\ncurrent_F: " << current_F
                        << "\ninvert(current_F): " << invert(current_F)
                        << std::endl;
            }
  
            constitutive_request.set_deformation_Jacobian(projected_jacobian);
            constitutive_request.set_unprojected_deformation_Jacobian(determinant(updated_F) * material_Jacobian);
            constitutive_request.set_temperature(projected_temperature);
            constitutive_request.set_deformation_gradient(deformation_gradient_increment);
            constitutive_request.set_time_increment(time_increment);
  
            try {
              material.compute_constitutive_request(constitutive_request,
                                                    quadrature_point_index);
            } catch (const MaterialDomainException &exc) {

std::cerr << "projected_jacobian: " << projected_jacobian << "\ndeformation_gradient_increment: " << deformation_gradient_increment << "\nupdated_F: " << updated_F << "\nf_m_n_plus_1: " << f_m_n_plus_1 << "\ncurrent_F: " << current_F << "\ninvert(current_F): " << invert(current_F) << "\nJacobian: " << Jacobian << "\ndeterminant(f_m_n_plus_1): " << determinant(f_m_n_plus_1) << "\nprevious_Jacobian: " << previous_Jacobian << "\n-------------------\n" << std::endl; std::cerr << exc.what() << std::endl;

              kinematic_domains_are_valid = false;
              continue;
            }
  
            previous_constitutive_request.set_deformation_Jacobian(projected_previous_jacobian);
            previous_constitutive_request.set_temperature(projected_previous_temperature);
            previous_constitutive_request.set_deformation_gradient(unit_symmetric_tensor<dim+1, Number>());
            previous_constitutive_request.set_time_increment(time_increment);
            previous_constitutive_request.set_is_plastic(false); // the elastic strain is known, no need for predictor-corrector procedure
  
            material.compute_constitutive_request(previous_constitutive_request,
                                                  quadrature_point_index);
  
            for (unsigned int j = 0; j < dofs_per_cell; ++j) {
              projected_strain_divergence[j] =
                mixed_values(0) * projected_strain_divergence_coefficients[j][0];
              if (fill_system_matrix) {
                projected_jacobian_tangent[j] =
                  mixed_values(0) * projected_jacobian_tangent_coefficients[j][0];
                for (unsigned int k = 0; k < dofs_per_cell; ++k) {
                  projected_strain_divergence_tangent[j][k] =
                    mixed_values(0)
                    * projected_strain_divergence_tangent_coefficients[j][k][0];
                }
              }
            }
  
            for (unsigned int i = 1; i < mixed_dofs_per_cell; ++i) {
              for (unsigned int j = 0; j < dofs_per_cell; ++j) {
                projected_strain_divergence[j] +=
                  mixed_values(i) * projected_strain_divergence_coefficients[j][i];
                if (fill_system_matrix) {
                  projected_jacobian_tangent[j] +=
                    mixed_values(i) * projected_jacobian_tangent_coefficients[j][i];
                  for (unsigned int k = 0; k < dofs_per_cell; ++k) {
                    projected_strain_divergence_tangent[j][k] +=
                      mixed_values(i)
                      * projected_strain_divergence_tangent_coefficients[j][k][i];
                  }
                }
              }
            }
  
            const Number RJxW = radius / material_Jacobian * fe_values.JxW(q_point);
  
            for (unsigned int i = 0; i < dofs_per_cell; ++i) {
              const auto strain_i = symmetrize(rate_gradients[i] + angular_rate_gradients[i] * inv_updated_F);
  

stress deviator term

              const SymmetricTensor<2, dim+1, Number> stress_deviator =
                alpha_f * previous_constitutive_request.get_stress_deviator()
                + (1-alpha_f) * constitutive_request.get_stress_deviator();
              const Number pressure =
                alpha_f * previous_constitutive_request.get_pressure()
                + (1-alpha_f) * constitutive_request.get_pressure();
              cell_residual(i) += strain_i * stress_deviator * RJxW;
  

pressure term

              cell_residual(i) +=
                (projected_strain_divergence.at(i)) * pressure * RJxW;
  

body force term

              const unsigned int
              component_i = mech_fe.system_to_component_index(i).first;
              for ( typename std::vector<BodyForceApplier<dim, Number> >::const_iterator
                    bodyForceApplier = mechanical_lbc_system.bodyLoadAppliers.begin();
                    bodyForceApplier != mechanical_lbc_system.bodyLoadAppliers.end();
                    ++bodyForceApplier) {
                cell_residual(i) +=
                  bodyForceApplier->apply(
                    component_i,
                    fe_values.shape_value (i, q_point),
                    RJxW);
              }
  

inertial term

              cell_residual(i) +=
                (postprocess_tensor_dimension(fe_values[displacements].value(i, q_point)) + scalar_to_angular_tensor(fe_values[angular_velocity].value(i, q_point)))
                * DENSITY
                * acceleration_n_plus_1_minus_alpha_m * RJxW;
            }
  
            if (fill_system_matrix) {
              std::vector<SymmetricTensor<2, dim+1, Number> > stress_deviator_tangents(dofs_per_cell);
              std::vector<Number> pressure_tangents(dofs_per_cell);
              for (unsigned int i = 0; i < dofs_per_cell; ++i) {
  
                const Tensor<2, dim+1, Number> d_X_prime_d_X_variation =
                    deformation_gradient_from_angular_displacement_gradient_variations(
                          -angular_velocity_increments[q_point],
                          -angular_velocity_gradient_increments[q_point],
                          radius,
                          -fe_values[angular_velocity].value(i, q_point),
                          -fe_values[angular_velocity].gradient(i, q_point)
                        );
  
                const Tensor<2, dim+1, Number> rotation_to_X_prime_frame_variation = rotation_tensor_variation_to_transform_B_e(
                          -angular_velocity_increments[q_point],
                          -fe_values[angular_velocity].value(i, q_point),
                          radius
                        );
  
                const auto f_m_n_plus_1_variation_inv_f_m_n_plus_1 =
                  (- inv_d_X_prime_d_X * d_X_prime_d_X_variation * inv_d_X_prime_d_X * rotation_to_X_prime_frame
                    + inv_d_X_prime_d_X * rotation_to_X_prime_frame_variation) * invert(f_m_n_plus_1);
  
                stress_deviator_tangents[i] = (1-alpha_f) * constitutive_request.get_stress_deviator_tangent(
                                                rate_gradients[i]
                                                - Constants<dim+1, Number>::one_third()
                                                  * trace(rate_gradients[i])
                                                  * unit_symmetric_tensor<dim+1, Number>()
                                                + updated_F * f_m_n_plus_1_variation_inv_f_m_n_plus_1 * inv_updated_F
                                                - Constants<dim+1, Number>::one_third()
                                                  * trace(f_m_n_plus_1_variation_inv_f_m_n_plus_1)
                                                  * unit_symmetric_tensor<dim+1, Number>());
  
                pressure_tangents[i] = (1-alpha_f) * constitutive_request.get_pressure_tangent(
                                         projected_jacobian_tangent[i]);
              }
              for (unsigned int i = 0; i < dofs_per_cell; ++i) {
                for (unsigned int j = 0; j < dofs_per_cell; ++j) {

stress tangent

const Number f_int_dev_tau = symmetrize(rate_gradients[i] + angular_rate_gradients[i] * inv_updated_F) * stress_deviator_tangents[j];

cell_matrix(i, j) += (f_int_dev_tau) * RJxW;

pressure_tangent

                  const Number f_int_pressure =
                    projected_strain_divergence.at(i)
                    * pressure_tangents[j];
                  cell_matrix(i, j) += (f_int_pressure) * RJxW;

geometric_tangent

                  const Tensor<2, dim+1, Number> grad_ui_grad_uj = (rate_gradients[i] + angular_rate_gradients[i] * inv_updated_F) * rate_gradients[j];
                  const SymmetricTensor<2, dim+1, Number> sym_grad_ui_grad_uj = symmetrize(grad_ui_grad_uj);
                  const Number f_int_geom =
                    projected_strain_divergence_tangent[i][j]
                    * constitutive_request.get_pressure()
                    - sym_grad_ui_grad_uj
                    * constitutive_request.get_stress_deviator();
                  cell_matrix(i, j) += (f_int_geom) * RJxW;
  

inertial tangent

                  cell_matrix(i, j) +=
                    (postprocess_tensor_dimension(fe_values[displacements].value(i, q_point))
                     + scalar_to_angular_tensor(fe_values[angular_velocity].value(i, q_point)))
                    * DENSITY
                    * (acceleration_n_plus_1_minus_alpha_m_tangent_modulus
                        * (postprocess_tensor_dimension(fe_values[displacements].value(j, q_point))
                           + scalar_to_angular_tensor(fe_values[angular_velocity].value(j, q_point))))
                    * RJxW;
                }
              }
            }
          }
  
          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face) {
            for (auto boundaryForceSpec: mechanical_lbc_system.boundaryLoadAppliers) {
              if (cell->face(face)->boundary_id() == static_cast<types::boundary_id>(boundaryForceSpec.first)) {
                fe_face_values.reinit(cell, face);
                for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) {
                  for (unsigned int i = 0; i < dofs_per_cell; ++i) {
                    const unsigned int component_i = mech_fe.system_to_component_index(i).first;
                    cell_residual(i) +=
                      boundaryForceSpec.second.apply(
                        component_i,
                        fe_face_values.shape_value(i, q_point),
                        fe_face_values.quadrature_point(q_point)[0] * fe_face_values.JxW(q_point));
                  }
                }
              }
            }
  
            for(const auto boundary_unidirectional_penalty_spec: mechanical_lbc_system.boundary_unidirectional_penalty_specs) {
              if (cell->face(face)->boundary_id() == boundary_unidirectional_penalty_spec->get_boundary_id()) {
                const Number reference_displacement_increment = boundary_unidirectional_penalty_spec->get_reference_displacement_increment();
                const Number residual_force = boundary_unidirectional_penalty_spec->get_residual_force();
                const Number quadratic_spring_factor = boundary_unidirectional_penalty_spec->get_quadratic_spring_factor();
  
                fe_face_values.reinit(cell, face);
                fe_face_values[displacements].get_function_values(
                  Newton_system.current_increment,
                  face_displacement_value_increments);
                for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) {
                  const Tensor<1, dim, Number> surface_normal = fe_face_values.normal_vector(q_point);
                  const Number surface_displacement = face_displacement_value_increments[q_point] * surface_normal;
                  const Number surface_force =
                    -residual_force
                    + surface_displacement < -reference_displacement_increment?
                      0 :
                      -0.5 * quadratic_spring_factor * std::pow(surface_displacement + reference_displacement_increment, 2);
                  const Number surface_force_tangent_modulus =
                    surface_displacement < -reference_displacement_increment?
                    0 :
                    -quadratic_spring_factor * (surface_displacement + reference_displacement_increment);
                  for (unsigned int i = 0; i < dofs_per_cell; ++i) {
                    cell_residual(i) -=
                      fe_face_values[displacements].value(i, q_point)
                      * surface_force * surface_normal
                      * fe_face_values.quadrature_point(q_point)[0] * fe_face_values.JxW(q_point);
                    for (unsigned int j = 0; j < dofs_per_cell; ++j) {
                      cell_matrix(i, j) -=
                        fe_face_values[displacements].value(i, q_point)
                        * surface_force_tangent_modulus * (fe_face_values[displacements].value(j, q_point) * surface_normal) * surface_normal
                        * fe_face_values.quadrature_point(q_point)[0] * fe_face_values.JxW(q_point);
                    }
                  }
                }
              }
            }
          }
  

const Number relative_symmetry_norm2 = cell_matrix.relative_symmetry_norm2(); if(relative_symmetry_norm2 > 1e-8) std::cout << "relative_symmetry_norm2: " << cell_matrix.relative_symmetry_norm2() << std::endl;

          std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
          cell->get_dof_indices (local_dof_indices);
          if (fill_system_matrix) {
            mechanical_dof_system.nodal_constraints.distribute_local_to_global(
              cell_matrix,
              cell_residual,
              local_dof_indices,
              Newton_system.Newton_step_matrix,
              Newton_system.Newton_step_residual,
              true);
          } else {
            mechanical_dof_system.nodal_constraints.distribute_local_to_global(
              cell_residual, local_dof_indices,
              Newton_system.Newton_step_residual);
          }
        } /* if cell is locally owned */
      } /*for (; cell!=endc; ++cell)*/
  
      const unsigned short local_domain_is_valid = kinematic_domains_are_valid ? 1 : 0;
      unsigned short all_kinematic_domains_are_valid;
  

did any of the processes fail to assemble?

      MPI_Allreduce(
        &local_domain_is_valid,
        &all_kinematic_domains_are_valid, 1, MPI_UNSIGNED_SHORT,
        MPI_MIN, mpi_communicator);
  
      if (all_kinematic_domains_are_valid < 1) {
        throw std::runtime_error("The domain is not valid...");
      }
  
      if (fill_system_matrix) {
        Newton_system.Newton_step_matrix.compress(VectorOperation::add);
      }
      Newton_system.Newton_step_residual.compress(VectorOperation::add);
  
    } /*PlasticityLabProg<dim,Number>::assemble_mechanical_system()*/
  
  
    template <int dim, typename Number>
    void PlasticityLabProg<dim, Number>::assemble_thermal_system(
      NewtonStepSystem &Newton_system,
      NewtonStepSystem &mechanical_nonlinear_system,
      const DoFSystem<dim, Number> &thermal_dof_system,
      const LBCSystem<dim, Number, 1> &thermal_lbc_system,
      const DoFSystem<dim, Number> &mechanical_dof_system,
      const DoFSystem<dim, Number> &mixed_fe_dof_system,
      Material<dim+1, Number> &material,
      const std::vector< MixedFEProjector<dim, Number> > &mixed_fe_projector,
      const std::unordered_map<size_t, Tensor<1, dim+1, Number>> &material_area_factors,
      const bool fill_system_matrix) {
  
      FEValues<dim> fe_values(
        mapping,
        therm_fe,
        quadrature_formula,
        update_values | update_quadrature_points |  update_gradients | update_JxW_values);
  
      FEFaceValues<dim> fe_face_values(
        mapping,
        therm_fe,
        face_quadrature_formula,
        update_values
        | update_quadrature_points
        | update_JxW_values);
  
      FEValues<dim> fe_mech_values(
        mapping,
        mech_fe,
        quadrature_formula,
        update_values | update_quadrature_points | update_gradients);
  
      FEFaceValues<dim> mech_fe_face_values(
        mapping,
        mech_fe,
        face_quadrature_formula,
        update_values
        | update_gradients
        | update_quadrature_points
        | update_normal_vectors
        | update_JxW_values);
  
      FEValues<dim> mixed_fe_values(
        mapping,
        mixed_var_fe,
        quadrature_formula,
        update_values);
  
      const unsigned int dofs_per_cell = therm_fe.dofs_per_cell;
      const unsigned int n_q_points = quadrature_formula.size();
      const unsigned int n_face_q_points = face_quadrature_formula.size();
      const unsigned int mixed_dofs_per_cell = mixed_var_fe.dofs_per_cell;
  
      FullMatrix<Number> cell_matrix (dofs_per_cell, dofs_per_cell);
      Vector<Number> cell_residual (dofs_per_cell);
  
      Vector<Number> mixed_values (mixed_dofs_per_cell);
      std::vector<Number> weighted_updated_J_vec(mixed_dofs_per_cell),
          weighted_current_J_vec(mixed_dofs_per_cell),
          weighted_previous_J_vec(mixed_dofs_per_cell),
          weighted_updated_theta_vec(mixed_dofs_per_cell),
          weighted_previous_theta_vec(mixed_dofs_per_cell),
          weighted_J_time_rate_vec(mixed_dofs_per_cell);
      std::vector< std::vector<Number> > weighted_shape_values(
        dofs_per_cell,
        std::vector<Number>(mixed_dofs_per_cell) );
  
      std::vector<Number> qp_updated_J_values(n_q_points),
          qp_previous_J_values(n_q_points),
          qp_updated_theta_values(n_q_points),
          qp_previous_theta_values(n_q_points),
          qp_J_time_rates(n_q_points);
      std::vector<std::vector<Number> > qp_shape_values(
        dofs_per_cell,
        std::vector<Number>(n_q_points));
  
      std::vector< Tensor<1, dim, Number> > thermal_gradient_increment(n_q_points),
          current_thermal_gradient(n_q_points);
  
      std::vector< Number > current_temperature_values(n_q_points),
          temperature_values_increment(n_q_points);
  
      std::vector< Number > current_face_temperature_values(n_face_q_points),
          face_temperature_values_increment(n_face_q_points);
  
      std::vector< Tensor<2, dim, Number> > current_displacement_gradients(n_q_points),
          displacement_gradient_increments(n_q_points);
      std::vector< Tensor<1, dim, Number> > current_displacement_values(n_q_points),
          displacement_value_increments(n_q_points);
  
      std::vector< Tensor<1, dim, Number> > current_angular_velocity_gradients(n_q_points),
          angular_velocity_gradient_increments(n_q_points),
          angular_velocity_gradient_previous_time_rates(n_q_points);
  
      std::vector< Number > angular_velocity_increments(n_q_points),
          current_angular_velocities(n_q_points),
          angular_velocity_previous_time_rates(n_q_points);
  
      std::vector< Tensor<2, dim, Number> > current_face_displacement_gradients(n_face_q_points);
      std::vector< Tensor<2, dim, Number> > face_displacement_gradient_increments(n_face_q_points);
  
      std::vector< Tensor<1, dim, Number> > current_face_displacement_values(n_face_q_points);
      std::vector< Tensor<1, dim, Number> > face_displacement_value_increments(n_face_q_points);
  
      const FEValuesExtractors::Scalar temperature (0);
      const FEValuesExtractors::Vector displacements (0);
      const FEValuesExtractors::Scalar angular_velocity(dim);
      Newton_system.Newton_step_matrix = 0;
      Newton_system.Newton_step_residual = 0;
      Newton_system.Newton_step_matrix.compress(VectorOperation::insert);
      Newton_system.Newton_step_residual.compress(VectorOperation::insert);
  
      auto cell = thermal_dof_system.dof_handler.begin_active();
      auto endc = thermal_dof_system.dof_handler.end();
      auto mechanical_cell = mechanical_dof_system.dof_handler.begin_active();
      auto mixed_fe_cell = mixed_fe_dof_system.dof_handler.begin_active();
      for (; cell != endc; ++cell, ++mechanical_cell, ++mixed_fe_cell) {
        if (cell->is_locally_owned()) {
          cell_matrix = 0;
          cell_residual = 0;
  
          fe_values.reinit (cell);
          fe_mech_values.reinit (mechanical_cell);
          mixed_fe_values.reinit (mixed_fe_cell);
  
          fe_values[temperature].get_function_gradients(
            Newton_system.current_increment,
            thermal_gradient_increment);
  
          fe_values[temperature].get_function_gradients(
            Newton_system.previous_deformation,
            current_thermal_gradient);
  
          fe_values[temperature].get_function_values(
            Newton_system.current_increment,
            temperature_values_increment);
  
          fe_values[temperature].get_function_values(
            Newton_system.previous_deformation,
            current_temperature_values);
  
          fe_mech_values[displacements].get_function_gradients(
            mechanical_nonlinear_system.current_increment,
            displacement_gradient_increments);
  
          fe_mech_values[displacements].get_function_gradients(
            mechanical_nonlinear_system.previous_deformation,
            current_displacement_gradients);
  
          fe_mech_values[displacements].get_function_values(
            mechanical_nonlinear_system.current_increment,
            displacement_value_increments);
  
          fe_mech_values[displacements].get_function_values(
            mechanical_nonlinear_system.previous_deformation,
            current_displacement_values);
  

Angular velocity

          fe_mech_values[angular_velocity].get_function_gradients(
            mechanical_nonlinear_system.current_increment,
            angular_velocity_gradient_increments);
  
          fe_mech_values[angular_velocity].get_function_values(
            mechanical_nonlinear_system.current_increment,
            angular_velocity_increments);
  

get vectors for projection onto mixed fe values

          for (unsigned int q_point = 0; q_point < n_q_points;
               ++q_point) {
  
            const Number radius = fe_mech_values.quadrature_point(q_point)[0];
  
            const auto previous_F = get_deformation_gradient(
              current_displacement_gradients[q_point],
              current_displacement_values[q_point][0]/radius);
  
            const auto updated_F = get_deformation_gradient(
              current_displacement_gradients[q_point] + displacement_gradient_increments[q_point],
              (current_displacement_values[q_point][0] + displacement_value_increments[q_point][0])/radius);
  
            const auto deformation_gradient_increment = postprocess_tensor_dimension(
              displacement_gradient_increments[q_point],
              displacement_value_increments[q_point][0]/radius);
  
            const point_index_t quadrature_point_index = cell->user_index() + q_point;
            const Number material_Jacobian = material.get_material_Jacobian(quadrature_point_index) / determinant(previous_F);
  
            const Number Jacobian = material_Jacobian * determinant(updated_F);
            const Number previous_Jacobian = material_Jacobian * determinant(previous_F);
  
            qp_previous_J_values.at(q_point) = previous_Jacobian;
            qp_updated_J_values.at(q_point) = Jacobian;
            qp_J_time_rates.at(q_point) = Jacobian * trace(deformation_gradient_increment * invert(updated_F)) / time_increment;
  
            qp_previous_theta_values.at(q_point) = current_temperature_values.at(q_point);
            qp_updated_theta_values.at(q_point) = current_temperature_values.at(q_point) + temperature_values_increment.at(q_point);
            for (unsigned int i = 0; i < dofs_per_cell; ++i) {
              qp_shape_values.at(i).at(q_point) = fe_values[temperature].value(i, q_point);
            }
          }
  
          const unsigned int cell_index = cell->user_index() / n_q_points;
  
          mixed_fe_projector[cell_index].project(
            &weighted_updated_J_vec,
            qp_updated_J_values);
          mixed_fe_projector[cell_index].project(
            &weighted_J_time_rate_vec,
            qp_J_time_rates);
          mixed_fe_projector[cell_index].project(
            &weighted_updated_theta_vec,
            qp_updated_theta_values);
          mixed_fe_projector[cell_index].project(
            &weighted_previous_J_vec,
            qp_previous_J_values);
          mixed_fe_projector[cell_index].project(
            &weighted_previous_theta_vec,
            qp_previous_theta_values);
          for (unsigned int i = 0; i < dofs_per_cell; ++i) {
            mixed_fe_projector[cell_index].project(
              &weighted_shape_values.at(i),
              qp_shape_values.at(i));
          }
  
          for (unsigned int q_point = 0; q_point < n_q_points;
               ++q_point) {
            const ConstitutiveModelUpdateFlags material_update_flags =
              fill_system_matrix ?
              (update_heat_flux | update_heat_flux_tangent
               | update_mechanical_dissipation
               | update_mechanical_dissipation_tangent
               | update_stored_heat | update_stored_heat_tangent)
              :
              (update_heat_flux | update_mechanical_dissipation
               | update_stored_heat);
  
            const ConstitutiveModelUpdateFlags heating_update_flags =
              fill_system_matrix ?
              (update_thermoelastic_heating
               | update_thermoelastic_heating_tangent)
              :
              (update_thermoelastic_heating);
  
            point_index_t quadrature_point_index = cell->user_index() + q_point;
            ConstitutiveModelRequest<dim+1, Number> constitutive_request(material_update_flags);
            ConstitutiveModelRequest<dim+1, Number> heating_request(heating_update_flags);
  
            for (unsigned int i = 0; i < mixed_dofs_per_cell; ++i) {
              mixed_values(i) = mixed_fe_values.shape_value(i, q_point);
            }
  
            Number projected_updated_J = 0;
            Number projected_previous_J = 0;
            Number projected_updated_theta = 0;
            Number projected_previous_theta = 0;
            Number projected_J_time_rate = 0;
            Vector<Number> projected_shape_values(dofs_per_cell);
            for (unsigned int i = 0; i < mixed_dofs_per_cell; ++i) {
              projected_updated_J +=
                mixed_values(i) * weighted_updated_J_vec.at(i);
              projected_J_time_rate +=
                mixed_values(i) * weighted_J_time_rate_vec.at(i);
              projected_updated_theta +=
                mixed_values(i) * weighted_updated_theta_vec.at(i);
              projected_previous_J +=
                mixed_values(i) * weighted_previous_J_vec.at(i);
              projected_previous_theta +=
                mixed_values(i) * weighted_previous_theta_vec.at(i);
              for (unsigned int j = 0; j < dofs_per_cell; ++j) {
                projected_shape_values(j) +=
                  mixed_values(i) * weighted_shape_values.at(j).at(i);
              }
            }
  
            const Number d_theta_dt_n_plus_1 = temperature_values_increment[q_point] / time_increment;
  
            const Number d_theta_dt_tangent = 1.0 / time_increment;
  
            heating_request.set_deformation_Jacobian(projected_updated_J);
            heating_request.set_deformation_Jacobian_time_rate(projected_J_time_rate);
            heating_request.set_temperature(projected_updated_theta);
            heating_request.set_previous_deformation_Jacobian(projected_previous_J);
            heating_request.set_previous_temperature(projected_previous_theta);
            heating_request.set_time_increment(time_increment);
  
            const auto current_F = get_deformation_gradient(
              current_displacement_gradients[q_point],
              current_displacement_values[q_point][0]/fe_mech_values.quadrature_point(q_point)[0]
              );
            const auto updated_F = get_deformation_gradient(
              current_displacement_gradients[q_point]
              + displacement_gradient_increments[q_point],
              (current_displacement_values[q_point][0]
                + displacement_value_increments[q_point][0])/fe_mech_values.quadrature_point(q_point)[0]);
  
            const auto inv_updated_F = invert(updated_F);
            const Number material_Jacobian = material.get_material_Jacobian(quadrature_point_index) / determinant(current_F);
            [[maybe_unused]] const Number previous_Jacobian = determinant(current_F) * material_Jacobian;
            [[maybe_unused]] const Number Jacobian = determinant(updated_F) * material_Jacobian;
  
            const Number radius = fe_values.quadrature_point(q_point)[0];
  
            const Tensor<2, dim+1, Number> d_X_prime_d_X = deformation_gradient_from_angular_displacement_gradient(
                      -angular_velocity_increments[q_point],
                      -angular_velocity_gradient_increments[q_point],
                      radius
                    );
  
            const Tensor<2, dim+1, Number> rotation_to_X_prime_frame = rotation_tensor_to_transform_B_e(
                      -angular_velocity_increments[q_point],
                      radius
                    );
  
            const Tensor<2, dim+1, Number> inv_d_X_prime_d_X = invert(d_X_prime_d_X);
  
            const auto f_m_n_plus_1 = inv_d_X_prime_d_X * rotation_to_X_prime_frame;
  

std::cout << "f_r: " << inv_d_X_prime_d_X << std::endl; std::cout << "R: " << previous_elastic_deformation_transformation_tensor << std::endl; std::cout << "f_m_n+1: " << f_m_n_plus_1 << std::endl;

            const auto unnormalized_deformation_gradient_increment = updated_F * f_m_n_plus_1 * invert(current_F);
  
            const auto deformation_gradient_increment =
              std::pow(determinant(unnormalized_deformation_gradient_increment),
                        -Constants<dim, Number>::one_third()) * unnormalized_deformation_gradient_increment;
  
            const Number previous_temperature = current_temperature_values[q_point];
            const Number updated_temperature = previous_temperature + temperature_values_increment[q_point];
            const auto thermal_gradient =
              postprocess_tensor_dimension(current_thermal_gradient[q_point] + thermal_gradient_increment[q_point]) * inv_updated_F;
  
            constitutive_request.set_deformation_gradient(deformation_gradient_increment);
            constitutive_request.set_temperature_time_rate(d_theta_dt_n_plus_1);
            constitutive_request.set_temperature(updated_temperature);
            constitutive_request.set_thermal_gradient(thermal_gradient);
            constitutive_request.set_time_increment(time_increment);
  
            material.compute_constitutive_request(
              constitutive_request,
              quadrature_point_index);
            material.compute_constitutive_request(
              heating_request,
              quadrature_point_index);
  
            std::vector<Tensor<1, dim+1, Number>> rate_gradients(dofs_per_cell);
            std::vector<Number> rate_temperatures(dofs_per_cell);
            for (unsigned int i = 0; i < dofs_per_cell; ++i) {
              rate_gradients[i] = postprocess_tensor_dimension(fe_values[temperature].gradient(i, q_point)) * inv_updated_F;
              rate_temperatures[i] =fe_values[temperature].value(i, q_point);
            }
  
            const Number RJxW = fe_values.quadrature_point(q_point)[0] / material_Jacobian * fe_values.JxW(q_point);
  
            for (unsigned int i = 0; i < dofs_per_cell; ++i) {
  

heat flux term

              const auto heat_flux = constitutive_request.get_heat_flux();
              cell_residual(i) += rate_gradients[i]
                                  * heat_flux
                                  * RJxW;
  

stored heat term

              cell_residual(i) += rate_temperatures[i]
                                  * constitutive_request.get_stored_heat_rate()
                                  * RJxW;
  

mechanical dissipation term

              cell_residual(i) -= rate_temperatures[i]
                                  * constitutive_request.get_mechanical_dissipation()
                                  * RJxW;
  

elastoplastic heating term

              cell_residual(i) += (projected_shape_values(i))
                                  * heating_request.get_thermo_elastic_heating()
                                  * RJxW;
  
              for ( typename std::vector<BodyForceApplier<dim, Number> >::const_iterator
                    bodyHeatSourceApplier = thermal_lbc_system.bodyLoadAppliers.cbegin();
                    bodyHeatSourceApplier != thermal_lbc_system.bodyLoadAppliers.cend();
                    ++bodyHeatSourceApplier) {
                cell_residual(i) += bodyHeatSourceApplier->apply(
                                      0, rate_temperatures[i],
                                      fe_values.quadrature_point(q_point)[0] * fe_values.JxW(q_point));
              }
            }
  
            if (fill_system_matrix) {
              for (unsigned int i = 0; i < dofs_per_cell; ++i) {
                for (unsigned int j = 0; j < dofs_per_cell; ++j) {
  

heat flux tangent

                  const Number f_int_q =
                    rate_gradients[i]
                    * constitutive_request.get_heat_flux_tangent(rate_gradients[j]);
                  cell_matrix(i, j) += f_int_q * RJxW;
  

stored heat rate tangent

                  const Number f_int_cThetaDot =
                    rate_temperatures[i]
                    * constitutive_request.get_stored_heat_rate_tangent(d_theta_dt_tangent*rate_temperatures[j]);
                  cell_matrix(i, j) += f_int_cThetaDot * RJxW;
  

mechanical dissipation tangent

                  const Number f_int_mech_dissipation =
                    rate_temperatures[i]
                    * constitutive_request.get_mechanical_dissipation_tangent(rate_temperatures[j]);
                  cell_matrix(i, j) -= f_int_mech_dissipation * RJxW;
  

elastoplastic heating tangent

                  const Number f_int_elastoplastic_heating =
                    (projected_shape_values(i))
                    * heating_request.get_thermo_elastic_heating_tangent(projected_shape_values(j));
                  cell_matrix(i, j) += f_int_elastoplastic_heating * RJxW;
                }
              }
            }
          }
  
          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face) {
            fe_face_values.reinit(cell, face);
            mech_fe_face_values.reinit(mechanical_cell, face);
  
            fe_face_values[temperature].get_function_values(
              Newton_system.previous_deformation,
              current_face_temperature_values);
            fe_face_values[temperature].get_function_values(
              Newton_system.current_increment,
              face_temperature_values_increment);
  
            mech_fe_face_values[displacements].get_function_gradients(
              mechanical_nonlinear_system.previous_deformation,
              current_face_displacement_gradients);
            mech_fe_face_values[displacements].get_function_gradients(
              mechanical_nonlinear_system.current_increment,
              face_displacement_gradient_increments);
  
            mech_fe_face_values[displacements].get_function_values(
              mechanical_nonlinear_system.previous_deformation,
              current_face_displacement_values);
            mech_fe_face_values[displacements].get_function_values(
              mechanical_nonlinear_system.current_increment,
              face_displacement_value_increments);
  
  
            for ( typename std::vector<std::pair<int, BodyForceApplier<dim, Number> > >::const_iterator
                  boundaryHeatSource = thermal_lbc_system.boundaryLoadAppliers.cbegin();
                  boundaryHeatSource != thermal_lbc_system.boundaryLoadAppliers.cend();
                  ++boundaryHeatSource) {
              if (cell->face(face)->boundary_id() == static_cast<types::boundary_id>(boundaryHeatSource->first)) {
                for (unsigned int i = 0; i < dofs_per_cell; ++i) {
                  for (unsigned int q_point = 0;
                       q_point < face_quadrature_formula.size();
                       ++q_point) {
  
  
                    const auto updated_F = get_deformation_gradient(
                      current_face_displacement_gradients[q_point]
                      + face_displacement_gradient_increments[q_point],
                      (current_face_displacement_values[q_point][0]
                        + face_displacement_value_increments[q_point][0])/fe_mech_values.quadrature_point(q_point)[0]);
  
  
                    const Number J = determinant(updated_F);
                    const Tensor<2, dim+1, Number> inv_deformation_gradient = invert(updated_F);
                    const Tensor<1, dim+1, Number> reference_normal = postprocess_tensor_dimension(mech_fe_face_values.normal_vector(q_point), 0);
                    const Tensor<1, dim+1, Number> F_inv_transpose_N = transpose(inv_deformation_gradient) * reference_normal;
                    const Number norm_F_inv_transpose_N = (F_inv_transpose_N).norm();
  
                    cell_residual(i) +=
                      boundaryHeatSource->second.apply(
                        0,
                        fe_face_values.shape_value(i, q_point),
                        norm_F_inv_transpose_N * J *
                        fe_face_values.quadrature_point(q_point)[0] * fe_face_values.JxW(q_point));
                  }
                }
              }
            }
  
            for ( typename std::vector<std::pair<int, ConvectionBoundaryConditionApplier<dim, Number> > >::const_iterator
                  convectionBC = thermal_lbc_system.convection_BC_appliers.cbegin();
                  convectionBC != thermal_lbc_system.convection_BC_appliers.cend();
                  ++convectionBC) {
              if (cell->face(face)->boundary_id() == static_cast<types::boundary_id>(convectionBC->first)) {
                for (unsigned int q_point = 0;
                     q_point < face_quadrature_formula.size();
                     ++q_point) {
  
                  const unsigned int cell_index = cell->user_index() / quadrature_formula.size();
                  const unsigned int surface_point_key =
                    cell_index * GeometryInfo<dim>::faces_per_cell * n_face_q_points
                    + face * n_face_q_points
                    + q_point;
  
                  const auto updated_F = get_deformation_gradient(
                    current_face_displacement_gradients[q_point]
                    + face_displacement_gradient_increments[q_point],
                    (current_face_displacement_values[q_point][0]
                      + face_displacement_value_increments[q_point][0])/fe_mech_values.quadrature_point(q_point)[0]);
  
                  [[maybe_unused]] const Number J = determinant(updated_F);
                  const Tensor<2, dim+1, Number> inv_deformation_gradient = invert(updated_F);
                  const Tensor<1, dim+1, Number> reference_normal = postprocess_tensor_dimension(mech_fe_face_values.normal_vector(q_point), 0);
                  const Tensor<1, dim+1, Number> F_inv_transpose_N = transpose(inv_deformation_gradient) * reference_normal;
                  [[maybe_unused]] const Number norm_F_inv_transpose_N = (F_inv_transpose_N).norm();
  
                  const Number RJxW = fe_face_values.quadrature_point(q_point)[0] / material_area_factors.at(surface_point_key).norm() * fe_face_values.JxW(q_point);
  
                  const Number updated_face_temperature_value =
                    current_face_temperature_values[q_point] +
                    face_temperature_values_increment[q_point];
                  for (unsigned int i = 0; i < dofs_per_cell; ++i) {
                    cell_residual(i) +=
                      convectionBC->second.apply(
                        0,
                        fe_face_values.shape_value(i, q_point),
                        updated_face_temperature_value,
                        RJxW);
                    if (fill_system_matrix) {
                      for (unsigned int j = 0; j < dofs_per_cell; ++j) {
                        cell_matrix(i, j) +=
                          convectionBC->second.apply_gradient(
                            0,
                            fe_face_values.shape_value(i, q_point),
                            fe_face_values.shape_value(j, q_point),
                            RJxW);
                      }
                    }
                  }
                }
              }
            }
          }
  
  
          std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
          cell->get_dof_indices (local_dof_indices);
          if (fill_system_matrix) {
            thermal_dof_system.nodal_constraints.distribute_local_to_global(
              cell_matrix,
              cell_residual,
              local_dof_indices,
              Newton_system.Newton_step_matrix,
              Newton_system.Newton_step_residual,
              true);
          } else {
            thermal_dof_system.nodal_constraints.distribute_local_to_global(
              cell_residual, local_dof_indices,
              Newton_system.Newton_step_residual);
          }
        } /*for (; cell!=endc; ++cell) if(cell->is_locally_owned())*/
      }
  
      if (fill_system_matrix) Newton_system.Newton_step_matrix.compress(VectorOperation::add);
      Newton_system.Newton_step_residual.compress(VectorOperation::add);
    }
  
    template <int dim, typename Number>
    void PlasticityLabProg<dim, Number>::assemble_mesh_motion_system(
          NewtonStepSystem &mesh_motion_nonlinear_system,
          const DoFSystem<dim, Number> &mesh_motion_dof_system,
          const LBCSystem<dim, Number, dim> &,
          const NewtonStepSystem &deformation_nonlinear_system,
          const DoFSystem<dim, Number> &deformation_dof_system,
          const DoFSystem<dim, Number> &mixed_fe_dof_system,
          const std::vector< MixedFEProjector<dim, Number> > &mixed_fe_projector,
          const bool fill_system_matrix) {
  
  
      const Number mesh_motion_mu = 1.0;
      const Number mesh_motion_kappa = 5.0;
      const Number cell_jacobian_exponent = -0.0;
  
  
      FEValues<dim> deformation_fe_values(
        mapping,
        mech_fe,
        quadrature_formula,
        update_values | update_gradients);
  
      FEValues<dim> mesh_motion_fe_values(
        mapping,
        mesh_motion_fe,
        quadrature_formula,
        update_values  | update_gradients | update_quadrature_points |
        update_jacobians | update_JxW_values);
  
      FEValues<dim> mixed_fe_values(
        mapping,
        mixed_var_fe,
        quadrature_formula,
        update_values);
  
      const unsigned int dofs_per_cell = mesh_motion_fe.dofs_per_cell;
      const unsigned int n_q_points = quadrature_formula.size();
      const unsigned int mixed_dofs_per_cell = mixed_var_fe.dofs_per_cell;
  
      std::vector< Tensor<2, dim, Number> > mesh_motion_gradient_increments(n_q_points);
      std::vector< Tensor<1, dim, Number> > mesh_motion_value_increments(n_q_points);
  
      std::vector< Tensor<2, dim, Number> > current_deformation_gradients(n_q_points);
      std::vector< Tensor<2, dim, Number> > deformation_gradient_increments(n_q_points);
  
      std::vector< Tensor<1, dim, Number> > current_deformation_values(n_q_points);
      std::vector< Tensor<1, dim, Number> > deformation_value_increments(n_q_points);
  
      std::vector< Number > mesh_motion_jacobians(n_q_points);
      std::vector< std::vector<Number> > strain_divergences(dofs_per_cell, std::vector<Number>(n_q_points));
      std::vector< std::vector<Number> > jacobian_tangents(dofs_per_cell, std::vector<Number>(n_q_points));
  
      std::vector< std::vector< std::vector < Number> > > strain_divergence_tangents(
        dofs_per_cell,
        std::vector< std::vector< Number> >(dofs_per_cell,std::vector<Number>(n_q_points)));
  
      std::vector<Number> projected_Jacobian_coefficients(mixed_dofs_per_cell);
      std::vector<std::vector<Number> > projected_strain_divergence_coefficients(
        dofs_per_cell,
        std::vector<Number>(mixed_dofs_per_cell));
      std::vector<std::vector<Number> > projected_jacobian_tangent_coefficients(
        dofs_per_cell,
        std::vector<Number>(mixed_dofs_per_cell));
      std::vector<std::vector<std::vector<Number> > > projected_strain_divergence_tangent_coefficients(
        dofs_per_cell, std::vector< std::vector< Number> >(
          dofs_per_cell,
          std::vector<Number>(mixed_dofs_per_cell)));
  
      std::vector<Number> projected_strain_divergence(dofs_per_cell);
      std::vector<Number> projected_jacobian_tangent(dofs_per_cell);
      std::vector<std::vector<Number> > projected_strain_divergence_tangent(dofs_per_cell, std::vector<Number>(dofs_per_cell));
  
      FullMatrix<Number> cell_matrix(dofs_per_cell, dofs_per_cell);
      Vector<Number> cell_residual(dofs_per_cell);
  
      Vector<Number> mixed_values (mixed_dofs_per_cell);
  
      const FEValuesExtractors::Vector displacements (0);
      const FEValuesExtractors::Scalar angular_velocity(dim);
      const FEValuesExtractors::Scalar temperature (0);
  
      mesh_motion_nonlinear_system.Newton_step_matrix = 0;
      mesh_motion_nonlinear_system.Newton_step_residual = 0;
  
      bool kinematic_domains_are_valid = true;  // innocent until proven guilty
  
      auto cell = mesh_motion_dof_system.dof_handler.begin_active();
      auto endc = mesh_motion_dof_system.dof_handler.end();
      auto deformation_cell = deformation_dof_system.dof_handler.begin_active();
      auto mixed_fe_cell = mixed_fe_dof_system.dof_handler.begin_active();
  
      for (; cell != endc; ++cell, ++deformation_cell, ++mixed_fe_cell) {
        if (cell->is_locally_owned()) {
  
          cell_matrix = 0;
          cell_residual = 0;
  
          mesh_motion_fe_values.reinit (cell);
          deformation_fe_values.reinit (deformation_cell);
          mixed_fe_values.reinit (mixed_fe_cell);
  
          deformation_fe_values[displacements].get_function_gradients(
            deformation_nonlinear_system.previous_deformation,
            current_deformation_gradients);
  
          deformation_fe_values[displacements].get_function_gradients(
            deformation_nonlinear_system.current_increment,
            deformation_gradient_increments);
  
          deformation_fe_values[displacements].get_function_values(
            deformation_nonlinear_system.previous_deformation,
            current_deformation_values);
  
          deformation_fe_values[displacements].get_function_values(
            deformation_nonlinear_system.current_increment,
            deformation_value_increments);
  
          mesh_motion_fe_values[displacements].get_function_gradients(
            mesh_motion_nonlinear_system.current_increment,
            mesh_motion_gradient_increments);
  
          mesh_motion_fe_values[displacements].get_function_values(
            mesh_motion_nonlinear_system.current_increment,
            mesh_motion_value_increments);
  

get vectors for projection onto mixed fe values

          for (unsigned int q_point = 0; q_point < n_q_points;
               ++q_point) {
  
            const auto mesh_motion_gradient = get_deformation_gradient(
                                                -mesh_motion_gradient_increments[q_point],
                                                -mesh_motion_value_increments[q_point][0]/mesh_motion_fe_values.quadrature_point(q_point)[0]);
  
            const Number Jacobian = determinant(mesh_motion_gradient);
            mesh_motion_jacobians.at(q_point) = Jacobian;
  
            const auto inv_mesh_motion_gradient = invert(mesh_motion_gradient);
            std::vector<Tensor<2, dim+1, Number>> rate_gradients(dofs_per_cell);
            for (unsigned int i = 0; i < dofs_per_cell; ++i) {
              rate_gradients[i] = postprocess_tensor_dimension(
                                    -mesh_motion_fe_values[displacements].gradient(i, q_point),
                                    -mesh_motion_fe_values[displacements].value(i, q_point)[0]
                                      /mesh_motion_fe_values.quadrature_point(q_point)[0]) * inv_mesh_motion_gradient;
            }
  
            for (unsigned int i = 0; i < dofs_per_cell; ++i) {
              const Number strain_divergence_i = trace(rate_gradients[i]);
              strain_divergences[i].at(q_point) = Jacobian * strain_divergence_i;
              if (fill_system_matrix) {
                jacobian_tangents[i].at(q_point) = Jacobian * strain_divergence_i;
                for (unsigned int j = 0; j < dofs_per_cell; ++j) {
                  strain_divergence_tangents[i][j].at(q_point) = Jacobian * (trace(rate_gradients[i]) * trace(rate_gradients[j]) - trace(rate_gradients[i] * rate_gradients[j]));
                }
              }
            }
          }
  
          const unsigned int cell_index = cell->user_index() / n_q_points;
  
          mixed_fe_projector[cell_index].project(
            &projected_Jacobian_coefficients,
            mesh_motion_jacobians);
          for (unsigned int i = 0; i < dofs_per_cell; ++i) {
            mixed_fe_projector[cell_index].project(
              &projected_strain_divergence_coefficients[i],
              strain_divergences[i]);
            if (fill_system_matrix) {
              mixed_fe_projector[cell_index].project(
                &projected_jacobian_tangent_coefficients[i],
                jacobian_tangents[i]);
              for (unsigned int j = 0; j < dofs_per_cell; ++j) {
                mixed_fe_projector[cell_index].project(
                  &projected_strain_divergence_tangent_coefficients[i][j],
                  strain_divergence_tangents[i][j]);
              }
            }
          }
  
          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
  
            [[maybe_unused]] const point_index_t quadrature_point_index = cell->user_index() + q_point;
  
            const Number cell_jacobian = determinant(static_cast<Tensor <2, dim, Number>>(mesh_motion_fe_values.jacobian(q_point)));
  
            for (unsigned int i = 0; i < mixed_dofs_per_cell; ++i) {
              mixed_values(i) = mixed_fe_values.shape_value(i, q_point);
            }
  
            const auto mesh_motion_gradient = get_deformation_gradient(
                                                -mesh_motion_gradient_increments[q_point],
                                                -mesh_motion_value_increments[q_point][0]/mesh_motion_fe_values.quadrature_point(q_point)[0]);
  
            const auto deformation_gradient = get_deformation_gradient(
                                                current_deformation_gradients[q_point] + deformation_gradient_increments[q_point],
                                                (current_deformation_values[q_point][0] + deformation_value_increments[q_point][0])
                                                / mesh_motion_fe_values.quadrature_point(q_point)[0]);
  
            const auto inv_mesh_motion_gradient = invert(mesh_motion_gradient);
            const auto inv_deformation_gradient = invert(deformation_gradient);
  
            [[maybe_unused]] const Number deformation_Jacobian = determinant(deformation_gradient);
            const Number mesh_motion_Jacobian = determinant(mesh_motion_gradient);
  
            Number projected_jacobian = 0;
            for (unsigned int i = 0; i < mixed_dofs_per_cell; ++i) {
              projected_jacobian +=
                mixed_values(i) * projected_Jacobian_coefficients.at(i);
            }
  
            for (unsigned int j = 0; j < dofs_per_cell; ++j) {
              projected_strain_divergence[j] =
                mixed_values(0) * projected_strain_divergence_coefficients[j][0];
              if (fill_system_matrix) {
                projected_jacobian_tangent[j] =
                  mixed_values(0) * projected_jacobian_tangent_coefficients[j][0];
                for (unsigned int k = 0; k < dofs_per_cell; ++k) {
                  projected_strain_divergence_tangent[j][k] =
                    mixed_values(0)
                    * projected_strain_divergence_tangent_coefficients[j][k][0];
                }
              }
            }
  
            for (unsigned int i = 1; i < mixed_dofs_per_cell; ++i) {
              for (unsigned int j = 0; j < dofs_per_cell; ++j) {
                projected_strain_divergence[j] +=
                  mixed_values(i) * projected_strain_divergence_coefficients[j][i];
                if (fill_system_matrix) {
                  projected_jacobian_tangent[j] +=
                    mixed_values(i) * projected_jacobian_tangent_coefficients[j][i];
                  for (unsigned int k = 0; k < dofs_per_cell; ++k) {
                    projected_strain_divergence_tangent[j][k] +=
                      mixed_values(i)
                      * projected_strain_divergence_tangent_coefficients[j][k][i];
                  }
                }
              }
            }
  
            std::vector<Tensor<2, dim+1, Number>> rate_gradients(dofs_per_cell);
            for (unsigned int i = 0; i < dofs_per_cell; ++i) {
              rate_gradients[i] = postprocess_tensor_dimension(
                                    -mesh_motion_fe_values[displacements].gradient(i, q_point),
                                    -mesh_motion_fe_values[displacements].value(i, q_point)[0]
                                      /mesh_motion_fe_values.quadrature_point(q_point)[0]) * inv_mesh_motion_gradient;
            }
  
            const SymmetricTensor<2, dim+1, Number> stress_deviator =
                    std::pow(cell_jacobian, cell_jacobian_exponent) * mesh_motion_mu
                    * std::pow(mesh_motion_Jacobian, -Constants<dim, Number>::two_thirds())
                    * deviator(
                        symmetrize(
                          (deformation_gradient * mesh_motion_gradient)
                          * transpose(deformation_gradient * mesh_motion_gradient)));
  
            const Number pressure = mesh_motion_kappa * std::log(projected_jacobian);
  
            for (unsigned int i = 0; i < dofs_per_cell; ++i) {
              const auto strain_i = deformation_gradient * rate_gradients[i] * inv_deformation_gradient;
  

stress deviator term

              cell_residual(i) += symmetrize(strain_i) * stress_deviator
                                  * mesh_motion_fe_values.quadrature_point(q_point)[0] * mesh_motion_fe_values.JxW(q_point);
  

pressure term

              cell_residual(i) += (projected_strain_divergence.at(i) * pressure)
                                    * mesh_motion_fe_values.quadrature_point(q_point)[0] * mesh_motion_fe_values.JxW(q_point);
  
              if (fill_system_matrix) {
  
                for (unsigned int j = 0; j < dofs_per_cell; ++j) {
                  const auto strain_j = deformation_gradient * rate_gradients[j] * inv_deformation_gradient;
  

stress tangent

                  const SymmetricTensor<2, dim+1, Number> stress_deviator_tangent_j =
                      std::pow(cell_jacobian, cell_jacobian_exponent) * mesh_motion_mu
                      * std::pow(mesh_motion_Jacobian, -Constants<dim, Number>::two_thirds())
                      * deviator(
                          symmetrize(
                            2 * (strain_j - Constants<dim, Number>::one_third() * trace(strain_j) * unit_symmetric_tensor<dim+1, Number>())
                            * (deformation_gradient * mesh_motion_gradient)
                            * transpose(deformation_gradient * mesh_motion_gradient)));
  
                  cell_matrix(i, j) += (symmetrize(strain_i) * stress_deviator_tangent_j - symmetrize(strain_i * strain_j) * stress_deviator)
                                       * mesh_motion_fe_values.quadrature_point(q_point)[0] * mesh_motion_fe_values.JxW(q_point);
  

pressure_tangent

                  const Number pressure_tangent_j = mesh_motion_kappa * (1.0 / projected_jacobian) * projected_jacobian_tangent[j];
                  cell_matrix(i, j) += (projected_strain_divergence.at(i) * pressure_tangent_j + projected_strain_divergence_tangent[i][j] * pressure)
                                        * mesh_motion_fe_values.quadrature_point(q_point)[0] * mesh_motion_fe_values.JxW(q_point);
  
                }
              }
            }
          }
  

const Number relative_symmetry_norm2 = cell_matrix.relative_symmetry_norm2(); if(relative_symmetry_norm2 > 1e-8) std::cout << "relative_symmetry_norm2: " << cell_matrix.relative_symmetry_norm2() << std::endl;

          std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
          cell->get_dof_indices (local_dof_indices);
          if (fill_system_matrix) {
            mesh_motion_dof_system.nodal_constraints.distribute_local_to_global(
              cell_matrix,
              cell_residual,
              local_dof_indices,
              mesh_motion_nonlinear_system.Newton_step_matrix,
              mesh_motion_nonlinear_system.Newton_step_residual,
              true);
          } else {
            mesh_motion_dof_system.nodal_constraints.distribute_local_to_global(
              cell_residual, local_dof_indices,
              mesh_motion_nonlinear_system.Newton_step_residual);
          }
        } /* if cell is locally owned */
      } /*for (; cell!=endc; ++cell)*/
  
      const unsigned short local_domain_is_valid = kinematic_domains_are_valid ? 1 : 0;
      unsigned short all_kinematic_domains_are_valid;
  

did any of the processes fail to assemble?

      MPI_Allreduce(
        &local_domain_is_valid,
        &all_kinematic_domains_are_valid, 1, MPI_UNSIGNED_SHORT,
        MPI_MIN, mpi_communicator);
  
      if (all_kinematic_domains_are_valid < 1) {
        throw std::runtime_error("The domain is not valid...");
      }
  
      if (fill_system_matrix) {
        mesh_motion_nonlinear_system.Newton_step_matrix.compress(VectorOperation::add);
      }
      mesh_motion_nonlinear_system.Newton_step_residual.compress(VectorOperation::add);
  
    }
  
    template<typename BlockType>
    class SumOfMatrices : public EnableObserverPointer {
    public:
      SumOfMatrices(
          const BlockType &m1,
          const BlockType &m2):
        m1(m1),
        m2(m2) {
      }
  
      template<typename VectorType>
      void vmult(VectorType &dst, const VectorType &src) const {
        m1.vmult(dst, src);
        m2.vmult_add(dst, src);
      }
  
    private:
      const BlockType &m1;
      const BlockType &m2;
    };
  

TODO encorporate mechanical and thermal subsystems into structs and include functions in them

    template <int dim, typename Number>
    void PlasticityLabProg<dim, Number>::solve_system(
      const DoFSystem<dim, Number> &dof_system,
      NewtonStepSystem &nonlinear_system,
      const bool reset_solution) {
      TrilinosWrappers::PreconditionAMG preconditioner;
  
      const std::vector<std::vector<bool> > constant_modes
        = DoFTools::extract_constant_modes(dof_system.dof_handler,
                                           ComponentMask());
      
      TrilinosWrappers::PreconditionAMG::AdditionalData additional_data;
      additional_data.constant_modes = constant_modes;
      additional_data.elliptic = true;
      additional_data.n_cycles = 1;
      additional_data.w_cycle = false;
      additional_data.output_details = false;
      additional_data.smoother_sweeps = 2;
      additional_data.aggregation_threshold = 1e-2;
      preconditioner.initialize(nonlinear_system.Newton_step_matrix, additional_data);
  
      TrilinosWrappers::MPI::Vector tmp(dof_system.locally_owned_dofs, mpi_communicator);
      const Number relative_accuracy = 1e-08;
      const Number solver_tolerance  = relative_accuracy
                                       * nonlinear_system.Newton_step_matrix.residual(tmp, nonlinear_system.Newton_step_solution,
                                           nonlinear_system.Newton_step_residual);
      SolverControl solver_control(nonlinear_system.Newton_step_matrix.m(),
                                   solver_tolerance);
  
      SolverBicgstab<TrilinosWrappers::MPI::Vector> solver(solver_control);
  
      if (reset_solution) {
        nonlinear_system.Newton_step_solution = 0;
        nonlinear_system.Newton_step_solution.compress(VectorOperation::insert);
      }
  
      solver.solve(nonlinear_system.Newton_step_matrix, nonlinear_system.Newton_step_solution,
                   nonlinear_system.Newton_step_residual, preconditioner);
  
      pcout << "solved in " << solver_control.last_step() << " steps to residual value of " << solver_control.last_value() << endl;
      pcout << "solution norm is: " << nonlinear_system.Newton_step_solution.l2_norm() << endl;
  
      dof_system.nodal_constraints.distribute (nonlinear_system.Newton_step_solution);
    } /*ElasticProblem<dim,Number>::solve_system*/
  
  
    template <int dim, typename Number>
    void PlasticityLabProg<dim, Number>::get_plastic_strain(
      TrilinosWrappers::MPI::Vector &plastic_strain,
      const DoFHandler<dim> &discontinuous_dof_handler,
      const Material<dim+1, Number> &material,
      const std::vector< MixedFEProjector<dim, Number> > &qp_values_projectors) {
  
      const FiniteElement<dim> &fe = discontinuous_dof_handler.get_fe();
  
      const unsigned int n_q_points = quadrature_formula.size();
      const unsigned int disc_dofs_per_cell = fe.dofs_per_cell;
  
      std::vector<Number> plastic_strain_qp_values(n_q_points);
      std::vector<Number> projected_plastic_strains(disc_dofs_per_cell);
  
      auto cell = discontinuous_dof_handler.begin_active();
      auto endc = discontinuous_dof_handler.end();
      for (; cell != endc; ++cell) {
        if (cell->is_locally_owned()) {
          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
            const point_index_t quadrature_point_index = cell->user_index() + q_point;
            plastic_strain_qp_values[q_point] = std::exp(material.get_state_parameters(quadrature_point_index).at(0)) - 1;
          }
  
          const unsigned int cell_index = cell->user_index() / n_q_points;
          qp_values_projectors[cell_index].project(
            &projected_plastic_strains,
            plastic_strain_qp_values);
  
          std::vector<types::global_dof_index> local_dof_indices (disc_dofs_per_cell);
          cell->get_dof_indices (local_dof_indices);
          for (unsigned int dof_i = 0; dof_i < disc_dofs_per_cell; ++dof_i) {
            plastic_strain(local_dof_indices[dof_i]) = projected_plastic_strains.at(dof_i);
          }
        }
      }
    }
  
  
    template <int dim, typename Number>
    void PlasticityLabProg<dim, Number>::get_pressure(
      TrilinosWrappers::MPI::Vector &pressure,
      const DoFHandler<dim> &mixed_fe_dof_handler,
      TrilinosWrappers::MPI::Vector &von_mises_stress,
      const DoFHandler<dim> &discontinuous_dof_handler,
      NewtonStepSystem &Newton_system,
      const DoFSystem<dim, Number> &mechanical_dof_system,
      const NewtonStepSystem &thermal_Newton_system,
      const DoFSystem<dim, Number> &thermal_dof_system,
      Material<dim+1, Number> &material,
      const std::vector< MixedFEProjector<dim, Number> > &mixed_fe_projector,
      const std::vector< MixedFEProjector<dim, Number> > &qp_values_projectors) {
  
      FEValues<dim> fe_values(
        mapping,
        mech_fe,
        quadrature_formula,
        update_values | update_quadrature_points |  update_gradients);
  
      FEValues<dim> fe_therm_values(
        mapping,
        therm_fe,
        quadrature_formula,
        update_values);
  
      FEValues<dim> mixed_fe_values(
        mapping,
        mixed_var_fe,
        quadrature_formula,
        update_values);
  
      const unsigned int n_q_points = quadrature_formula.size();
      const unsigned int mixed_dofs_per_cell = mixed_var_fe.dofs_per_cell;
      const unsigned int disc_dofs_per_cell = discontinuous_dof_handler.get_fe().dofs_per_cell;
  
  
      std::vector< Tensor<2, dim, Number> > current_displacement_gradients(n_q_points);
      std::vector< Tensor<2, dim, Number> > displacement_gradient_increments(n_q_points);
      std::vector< Tensor<1, dim, Number> > current_displacement_values(n_q_points);
      std::vector< Tensor<1, dim, Number> > displacement_value_increments(n_q_points);
  
      std::vector< Number > current_temperature_values(n_q_points);
      std::vector< Number > updated_temperature_increments(n_q_points);
      std::vector< Number > deformation_jacobians(n_q_points);
      std::vector< Number > pressure_values(n_q_points);
      std::vector< Number > von_mises_stress_values(n_q_points);
  
      std::vector<Number> projected_temperature_coefficients(mixed_dofs_per_cell);
      std::vector<Number> projected_Jacobian_coefficients(mixed_dofs_per_cell);
      std::vector<Number> projected_pressure_coefficients(mixed_dofs_per_cell);
      std::vector<Number> projected_von_mises_stress_coefficients(disc_dofs_per_cell);
  
      Vector<Number> mixed_values(mixed_dofs_per_cell);
  
      const FEValuesExtractors::Vector displacements (0);
      const FEValuesExtractors::Scalar temperature (0);
  
      auto cell = mechanical_dof_system.dof_handler.begin_active();
      auto endc = mechanical_dof_system.dof_handler.end();
      auto thermal_cell = thermal_dof_system.dof_handler.begin_active();
      auto mixed_fe_cell = mixed_fe_dof_handler.begin_active();
      auto discontinuous_fe_cell = discontinuous_dof_handler.begin_active();
      for (; cell != endc; ++cell, ++thermal_cell, ++mixed_fe_cell, ++discontinuous_fe_cell) {
        if (cell->is_locally_owned()) {
  
          fe_values.reinit (cell);
          fe_therm_values.reinit (thermal_cell);
          mixed_fe_values.reinit (mixed_fe_cell);
  
          fe_values[displacements].get_function_gradients(
            Newton_system.current_increment,
            displacement_gradient_increments);
  
          fe_values[displacements].get_function_gradients(
            Newton_system.previous_deformation,
            current_displacement_gradients);
  
          fe_values[displacements].get_function_values(
            Newton_system.current_increment,
            displacement_value_increments);
  
          fe_values[displacements].get_function_values(
            Newton_system.previous_deformation,
            current_displacement_values);
  
          fe_therm_values[temperature].get_function_values (
            thermal_Newton_system.previous_deformation,
            current_temperature_values);
  
          fe_therm_values[temperature].get_function_values (
            thermal_Newton_system.current_increment,
            updated_temperature_increments);
  

get vectors for projection onto mixed fe values

          for (unsigned int q_point = 0; q_point < n_q_points;
               ++q_point) {
            updated_temperature_increments.at(q_point) += current_temperature_values.at(q_point);
  
            const auto current_F = get_deformation_gradient(
              current_displacement_gradients[q_point],
              current_displacement_values[q_point][0]/fe_values.quadrature_point(q_point)[0]
              );
            const auto updated_F = get_deformation_gradient(
              current_displacement_gradients[q_point]
              + displacement_gradient_increments[q_point],
              (current_displacement_values[q_point][0]
                + displacement_value_increments[q_point][0])/fe_values.quadrature_point(q_point)[0]);
  
            const point_index_t quadrature_point_index = cell->user_index() + q_point;
            const Number material_Jacobian = material.get_material_Jacobian(quadrature_point_index) / determinant(current_F);
  
            deformation_jacobians.at(q_point) = material_Jacobian * determinant(updated_F);
          }
  
          const unsigned int cell_index = cell->user_index() / n_q_points;
  
          mixed_fe_projector[cell_index].project(
            &projected_temperature_coefficients,
            updated_temperature_increments);
          mixed_fe_projector[cell_index].project(
            &projected_Jacobian_coefficients,
            deformation_jacobians);
  
          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
  
            const point_index_t quadrature_point_index = cell->user_index() + q_point;
            ConstitutiveModelRequest<dim+1, Number> constitutive_request(update_pressure | update_stress_deviator);
  
            for (unsigned int i = 0; i < mixed_dofs_per_cell; ++i) {
              mixed_values(i) = mixed_fe_values.shape_value(i, q_point);
            }
  
            const auto current_F = get_deformation_gradient(
              current_displacement_gradients[q_point],
              current_displacement_values[q_point][0]/fe_values.quadrature_point(q_point)[0]
              );
            const auto updated_F = get_deformation_gradient(
              current_displacement_gradients[q_point]
              + displacement_gradient_increments[q_point],
              (current_displacement_values[q_point][0]
                + displacement_value_increments[q_point][0])/fe_values.quadrature_point(q_point)[0]);
  
            [[maybe_unused]] const auto inv_updated_F = invert(updated_F);
            const Number Jacobian = determinant(updated_F);
            const Number previous_Jacobian = determinant(current_F);
  
            const auto deformation_gradient_increment = std::pow(Jacobian / previous_Jacobian, -Constants<dim, Number>::one_third()) * updated_F * invert(current_F);
  
            Number projected_jacobian = 0;
            Number projected_temperature = 0;
            for (unsigned int i = 0; i < mixed_dofs_per_cell; ++i) {
              projected_jacobian +=
                mixed_values(i) * projected_Jacobian_coefficients.at(i);
              projected_temperature +=
                mixed_values(i) * projected_temperature_coefficients.at(i);
            }
  
            constitutive_request.set_deformation_Jacobian(projected_jacobian);
            constitutive_request.set_temperature(projected_temperature);
            constitutive_request.set_deformation_gradient(deformation_gradient_increment);
            constitutive_request.set_time_increment(time_increment);
  
            material.compute_constitutive_request(constitutive_request,
                                                  quadrature_point_index);
  

pressure term

            pressure_values.at(q_point) = constitutive_request.get_pressure();
            von_mises_stress_values.at(q_point) =
              constitutive_request.get_stress_deviator().norm() / Constants<dim, Number>::sqrt2thirds();
          }
  
          mixed_fe_projector[cell_index].project(
            &projected_pressure_coefficients,
            pressure_values);
  
          qp_values_projectors[cell_index].project(
            &projected_von_mises_stress_coefficients,
            von_mises_stress_values);
  
          std::vector<types::global_dof_index> local_dof_indices (mixed_dofs_per_cell);
          mixed_fe_cell->get_dof_indices (local_dof_indices);
          for (unsigned int dof_i = 0; dof_i < mixed_dofs_per_cell; ++dof_i) {
            pressure(local_dof_indices[dof_i]) = projected_pressure_coefficients.at(dof_i);
          }
  
          std::vector<types::global_dof_index> local_discontinuous_dof_indices(disc_dofs_per_cell);
          discontinuous_fe_cell->get_dof_indices (local_discontinuous_dof_indices);
          for (unsigned int dof_i = 0; dof_i < disc_dofs_per_cell; ++dof_i) {
            von_mises_stress(local_discontinuous_dof_indices[dof_i]) =
              projected_von_mises_stress_coefficients.at(dof_i);
          }
  
        } /* if cell is locally owned */
      } /*for (; cell!=endc; ++cell)*/
    }
  
    template <int dim, typename Number>
    void PlasticityLabProg<dim, Number>::prepare_output_results(
      DataOut<dim> &data_out,
      const DoFSystem<dim, Number> &mechanical_dof_system,
      const NewtonStepSystem &mechanical_nonlinear_system,
      const DoFSystem<dim, Number> &thermal_dof_system,
      const NewtonStepSystem &thermal_nonlinear_system,
      const DoFSystem<dim, Number> &mesh_motion_dof_system,
      const NewtonStepSystem &mesh_motion_nonlinear_system) const {
  
      std::vector<std::string> displacement_names(dim, "displacement");
      displacement_names.emplace_back("angular_displacement");
      std::vector<std::string> velocity_names(dim, "displacement_time_rate");
      velocity_names.emplace_back("angular_velocity");
      std::vector<DataComponentInterpretation::DataComponentInterpretation>
        data_component_interpretation(
          dim, DataComponentInterpretation::component_is_part_of_vector);
      data_component_interpretation.push_back(
        DataComponentInterpretation::component_is_scalar);
  
      std::vector<DataComponentInterpretation::DataComponentInterpretation>
        mesh_motion_data_component_interpretation(
          dim, DataComponentInterpretation::component_is_part_of_vector);
  
      data_out.add_data_vector(mechanical_dof_system.dof_handler,
                               mechanical_nonlinear_system.previous_deformation,
                               displacement_names,
                               data_component_interpretation);
      data_out.add_data_vector(mechanical_dof_system.dof_handler,
                               mechanical_nonlinear_system.previous_time_derivative,
                               velocity_names,
                               data_component_interpretation);
      data_out.add_data_vector(thermal_dof_system.dof_handler,
                               thermal_nonlinear_system.previous_deformation,
                               "Temperature");
      data_out.add_data_vector(mesh_motion_dof_system.dof_handler,
                               mesh_motion_nonlinear_system.previous_deformation,
                               std::vector<std::string>(dim, "mesh_motion"),
                               mesh_motion_data_component_interpretation);
      data_out.add_data_vector(mesh_motion_dof_system.dof_handler,
                               mesh_motion_nonlinear_system.previous_time_derivative,
                               std::vector<std::string>(dim, "mesh_velocity"),
                               mesh_motion_data_component_interpretation);
      data_out.build_patches(mapping, 2);
    }
  
    template <int dim, typename Number>
    template <typename TriangulationType>
    void PlasticityLabProg<dim, Number>::write_output_results(
      DataOut<dim> &data_out,
      const TriangulationType &tria,
      const std::string &filename_base) const {
  
      const std::string filename =
        (filename_base + "-"
         + Utilities::int_to_string(tria.locally_owned_subdomain(), 4));
  
      std::ofstream output_vtu((filename + ".vtu").c_str());
      data_out.write_vtu(output_vtu);
  
      if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) {
        std::vector<std::string> filenames;
        for (unsigned int i = 0;
             i < Utilities::MPI::n_mpi_processes(mpi_communicator); ++i)
          filenames.push_back(filename_base + "-"
                              + Utilities::int_to_string(i, 4)
                              + ".vtu");
        std::ofstream pvtu_master_output((filename_base + ".pvtu").c_str());
        data_out.write_pvtu_record(pvtu_master_output, filenames);
        std::ofstream visit_master_output((filename_base + ".visit").c_str());
        DataOutBase::write_visit_record(visit_master_output, filenames);
      }
    } /* output_results */
  
    template class PlasticityLabProg<2, double>;
  
  } /*namespace PlasticityLab*/

Annotated version of src/PlasticityLabProg.h

  /*
   * PlasticityLabProg.h
   *
   *  Created on: 09 Jul 2014
   *      Author: cerecam
   */
  
  #ifndef PLASTICITYLABPROG_H_
  #define PLASTICITYLABPROG_H_
  
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/fe/fe_dgp.h>
  #include <deal.II/fe/fe_system.h>
  #include <deal.II/fe/mapping_q.h>
  
  #include <deal.II/distributed/tria.h>
  
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/base/function.h>
  #include <stdexcept>
  
  #include "DoFSystem.h"
  #include "LBCSystem.h"
  #include "MixedFEProjector.h"
  #include "Material.h"
  #include "NewtonStepSystem.h"
  
  namespace PlasticityLab {
    using namespace dealii;
  
    template <int dim, typename Number = double>
    class PlasticityLabProg {
     public:
      PlasticityLabProg(Material<dim+1, Number> &);
      virtual ~PlasticityLabProg();
      void run();
  
     private:
      void make_grid (int);
      void make_cylindrical_grid(
        Triangulation<dim> &triangulation,
        LBCSystem<dim, Number, dim+1> &mech_lbc_system,
        LBCSystem<dim, Number, 1> &therm_lbc_system,
        int n_initial_global_refinements);
      void make_grid_();
      void make_necking_grid(
        Triangulation<dim> &triangulation,
        LBCSystem<dim, Number, dim+1> &mech_lbc_system,
        LBCSystem<dim, Number, 1> &therm_lbc_system,
        int n_initial_global_refinements);
      void make_interference_cylinder_grid(
        Triangulation<dim> &triangulation,
        LBCSystem<dim, Number, dim+1> &mech_lbc_system,
        LBCSystem<dim, Number, 1> &therm_lbc_system,
        int n_initial_global_refinements);
      void make_cylindrical_impact_grid(
        Triangulation<dim> &triangulation,
        LBCSystem<dim, Number, dim+1> &mech_lbc_system,
        LBCSystem<dim, Number, 1> &therm_lbc_system,
        int n_initial_global_refinements);
      void make_ball_in_hypershell_grid(
        Triangulation<dim> &triangulation,
        LBCSystem<dim, Number, dim+1> &mech_lbc_system,
        LBCSystem<dim, Number, 1> &therm_lbc_system,
        int n_initial_global_refinements);
      void make_hook_membrane_grid(int);
  
      void set_mesh_motion_LBCs(
        Triangulation<dim> &triangulation,
        LBCSystem<dim, Number, dim> &mesh_motion_lbc_system);
  
  
      Tensor<2, dim+1, Number> get_rotation_tensor(const Tensor<2, dim+1, Number> &skew_symmetric_rotation) const;
      Tensor<2, dim+1, Number> get_rotation_tensor_variation(
        const Tensor<2, dim+1, Number> &skew_symmetric_rotation,
        const Tensor<2, dim+1, Number> &skew_symmetric_rotation_variation) const;
  
  
      void remap_material_state_variables(
        const NewtonStepSystem &mesh_motion_nonlinear_system,
        const DoFSystem<dim, Number> &mesh_motion_dof_system,
        const NewtonStepSystem &mechanical_nonlinear_system,
        const DoFSystem<dim, Number> &mechanical_dof_system,
        const DoFSystem<dim, Number> &mixed_fe_dof_system,
        const std::vector< MixedFEProjector<dim, Number> > &mixed_fe_projector,
        Material<dim+1, Number> &material,
        std::unordered_map<point_index_t, Tensor<2, dim+1, Number>> &remapped_deformation_gradients);
  
  
      void remap_thermal_field(
        NewtonStepSystem &thermal_nonlinear_system,
        const DoFSystem<dim, Number> &thermal_dof_system,
        const NewtonStepSystem &mesh_motion_nonlinear_system,
        const DoFSystem<dim, Number> &mesh_motion_dof_system);
  
  
      void remap_mechanical_fields(
        NewtonStepSystem &mechanical_nonlinear_system,
        const DoFSystem<dim, Number> &mechanical_dof_system,
        const NewtonStepSystem &mesh_motion_nonlinear_system,
        const DoFSystem<dim, Number> &mesh_motion_dof_system);
  
  
      template <typename TriangulationType, typename MaterialType>
      void setup_material_data(TriangulationType &triangulation,
                               MaterialType &material);
  
      void setup_material_area_factors(
        const DoFSystem<dim, Number> &mesh_motion_dof_system,
        std::unordered_map<size_t, Tensor<1, dim+1, Number>> &material_area_factors);
  
      void update_material_area_factors(
          const NewtonStepSystem &mesh_motion_nonlinear_system,
          const DoFSystem<dim, Number> &mesh_motion_dof_system,
          std::unordered_map<size_t, Tensor<1, dim+1, Number>> &material_area_factors);
  
      template <typename TriangulationType>
      void setup_mixed_fe_projection_data(
        const TriangulationType &triangulation,
        std::vector< MixedFEProjector<dim, Number> > &MixedFeProjectors,
        const FiniteElement<dim> &MixedFE,
        const Quadrature<dim> &quadrature_formula);
  
      void assemble_mechanical_system(
        NewtonStepSystem &Newton_system,
        const DoFSystem<dim, Number> &mechanical_dof_system,
        const LBCSystem<dim, Number, dim+1> &mechanical_lbc_system,
        const NewtonStepSystem &mesh_motion_nonlinear_system,
        const DoFSystem<dim, Number> &mesh_motion_dof_system,
        const NewtonStepSystem &thermal_Newton_system,
        const DoFSystem<dim, Number> &thermal_dof_system,
        const DoFSystem<dim, Number> &mixed_fe_dof_system,
        Material<dim+1, Number> &material,
        const std::vector< MixedFEProjector<dim, Number> > &mixed_fe_projector,
        const bool fill_system_matrix = true,
        const bool update_material_state = false);
  
      void assemble_thermal_system(
        NewtonStepSystem &Newton_system,
        NewtonStepSystem &mechanical_nonlinear_system,
        const DoFSystem<dim, Number> &thermal_dof_system,
        const LBCSystem<dim, Number, 1> &thermal_lbc_system,
        const DoFSystem<dim, Number> &mechanical_dof_system,
        const DoFSystem<dim, Number> &mixed_fe_dof_system,
        Material<dim+1, Number> &material,
        const std::vector< MixedFEProjector<dim, Number> > &mixed_fe_projector,
        const std::unordered_map<size_t, Tensor<1, dim+1, Number>> &material_area_factors,
        const bool fill_system_matrix = true);
  
      void assemble_mesh_motion_system(
            NewtonStepSystem &mesh_motion_nonlinear_system,
            const DoFSystem<dim, Number> &mesh_motion_dof_system,
            const LBCSystem<dim, Number, dim> &mesh_motion_lbc_system,
            const NewtonStepSystem &deformation_nonlinear_system,
            const DoFSystem<dim, Number> &deformation_dof_system,
            const DoFSystem<dim, Number> &mixed_fe_dof_system,
            const std::vector< MixedFEProjector<dim, Number> > &mixed_fe_projector,
            const bool fill_system_matrix = true);
  
      void solve_system(const DoFSystem<dim, Number> &dof_system,
                        NewtonStepSystem &nonlinear_system,
                        const bool reset_solution=true);
  
      void prepare_output_results(DataOut<dim> &data_out,
                                  const DoFSystem<dim, Number> &dof_system,
                                  const NewtonStepSystem &nonlinear_system,
                                  const DoFSystem<dim, Number> &thermal_dof_system,
                                  const NewtonStepSystem &thermal_nonlinear_system,
                                  const DoFSystem<dim, Number> &mesh_motion_dof_system,
                                  const NewtonStepSystem &mesh_motion_nonlinear_system) const;
  
      template <typename TriangulationType>
      void write_output_results(DataOut<dim> &data_out,
                                const TriangulationType &tria,
                                const std::string &filename_base) const;
  
      void get_plastic_strain(
        TrilinosWrappers::MPI::Vector &plastic_strain,
        const DoFHandler<dim> &discontinuous_dof_handler,
        const Material<dim+1, Number> &material,
        const std::vector< MixedFEProjector<dim, Number> > &qp_values_projectors);
  
      void get_pressure(
        TrilinosWrappers::MPI::Vector &pressure,
        const DoFHandler<dim> &mixed_fe_dof_handler,
        TrilinosWrappers::MPI::Vector &von_mises_stress,
        const DoFHandler<dim> &discontinuous_dof_handler,
        NewtonStepSystem &Newton_system,
        const DoFSystem<dim, Number> &mechanical_dof_system,
        const NewtonStepSystem &thermal_Newton_system,
        const DoFSystem<dim, Number> &thermal_dof_system,
        Material<dim+1, Number> &material,
        const std::vector< MixedFEProjector<dim, Number> > &mixed_fe_projector,
        const std::vector< MixedFEProjector<dim, Number> > &qp_values_projectors);
  
      void solve_mechanical_step(int time_step);
      void solve_thermal_step(int time_step);
  
      void solve_mesh_motion_step(
            NewtonStepSystem &mesh_motion_nonlinear_system,
            const DoFSystem<dim, Number> &mesh_motion_dof_system,
            const LBCSystem<dim, Number, dim> &mesh_motion_lbc_system,
            const NewtonStepSystem &deformation_nonlinear_system,
            const DoFSystem<dim, Number> &deformation_dof_system,
            const DoFSystem<dim, Number> &mixed_fe_dof_system,
            const std::vector< MixedFEProjector<dim, Number> > &mixed_fe_projector,
            const int time_step);
  
      Tensor<2, dim+1, Number> get_deformation_gradient(
          const Tensor<2, dim, Number> &increment_gradient,
          const Number increment_0_over_rho) {
        Tensor<2, dim+1, Number> deformation_gradient = unit_symmetric_tensor<dim+1, Number>();
        for(unsigned int i=0; i<dim; ++i){
          for(unsigned int j=0; j<dim; ++j) {
            deformation_gradient[i][j] += increment_gradient[i][j];
          }
        }
        deformation_gradient[dim][dim] += increment_0_over_rho;
        return deformation_gradient;
      }
  
      Tensor<2, dim+1, Number> postprocess_tensor_dimension(
          const Tensor<2, dim, Number> &dimension_short_tensor,
          const Number entry_0_over_rho) {
        Tensor<2, dim+1, Number> postprocessed_tensor;
        for(unsigned int i=0; i<dim; ++i){
          for(unsigned int j=0; j<dim; ++j) {
            postprocessed_tensor[i][j] = dimension_short_tensor[i][j];
          }
        }
        postprocessed_tensor[dim][dim] = entry_0_over_rho;
        return postprocessed_tensor;
      }
  
      Tensor<1, dim+1, Number> postprocess_tensor_dimension(
            const Tensor<1, dim, Number> &dimension_short_tensor,
            const Number entry_at_dim=static_cast<Number>(0.0)) {
        Tensor<1, dim+1, Number> postprocessed_tensor;
        for(unsigned int i=0; i<dim; ++i){
          postprocessed_tensor[i] += dimension_short_tensor[i];
        }
        postprocessed_tensor[dim] = entry_at_dim;
        return postprocessed_tensor;
      }
  
      Tensor<1, dim+1, Number> scalar_to_angular_tensor(const Number angular_value) {
        Tensor<1, dim+1, Number> postprocessed_tensor = Tensor<1, dim+1, Number>();
        postprocessed_tensor[dim] = angular_value;
        return postprocessed_tensor;
      }
  
      Tensor<2, dim+1, Number> order_1_tensor_to_angular_gradient(
          const Tensor<1, dim, Number> &in_plane_gradient,
          const Number minus_entry_over_rho) {
        Tensor<2, dim+1, Number> postprocessed_tensor;
        for(unsigned int i=0; i<dim; ++i) {
          postprocessed_tensor[dim][i] = in_plane_gradient[i];
        }
        postprocessed_tensor[0][dim] = minus_entry_over_rho;
        return postprocessed_tensor;
      }
  
      Tensor<2, dim+1, Number> deformation_gradient_from_angular_displacement_gradient(
          const Number angular_displacement,
          const Tensor<1, dim, Number> &angular_displacement_gradient,
          const Number radius
        ) {
        if(2 != dim) {
          throw std::logic_error("Radial deformation gradient not implemented for dim!=2");
        }
  
        Tensor<2, dim+1, Number> result;
        const Number theta = angular_displacement / radius;
        Tensor<1, dim, Number> angular_displacmenet_over_r_squared;
        angular_displacmenet_over_r_squared[0] = angular_displacement / (radius * radius);
        const Tensor<1, dim, Number> theta_gradient = angular_displacement_gradient / radius - angular_displacmenet_over_r_squared;
  
        result[0][0] = std::cos(theta) - radius * std::sin(theta) * theta_gradient[0];
        result[0][1] = - radius * std::sin(theta) * theta_gradient[1];
        result[0][2] = -std::sin(theta);
  
        result[1][0] = 0;
        result[1][1] = 1;
        result[1][2] = 0;
  
        result[2][0] = std::sin(theta) + radius * std::cos(theta) * theta_gradient[0];
        result[2][1] = radius * std::cos(theta) * theta_gradient[1];
        result[2][2] = std::cos(theta);
  
        return result;
      }
  
      Tensor<2, dim+1, Number> deformation_gradient_from_angular_displacement_gradient_variations(
          const Number angular_displacement,
          const Tensor<1, dim, Number> &angular_displacement_gradient,
          const Number radius,
          const Number angular_displacement_variation,
          const Tensor<1, dim, Number> &angular_displacement_gradient_variation
        ) {
        if(2 != dim) {
          throw std::logic_error("Radial deformation gradient not implemented for dim!=2");
        }
  
        Tensor<2, dim+1, Number> result;
        const Number theta = angular_displacement / radius;
        const Number theta_variation = angular_displacement_variation / radius;
  
        Tensor<1, dim, Number> angular_displacmenet_over_r_squared;
        angular_displacmenet_over_r_squared[0] = angular_displacement / (radius * radius);
        const Tensor<1, dim, Number> theta_gradient = angular_displacement_gradient / radius - angular_displacmenet_over_r_squared;
  
        Tensor<1, dim, Number> angular_displacmenet_variation_over_r_squared;
        angular_displacmenet_variation_over_r_squared[0] = angular_displacement_variation / (radius * radius);
        const Tensor<1, dim, Number> theta_gradient_variation = angular_displacement_gradient_variation / radius - angular_displacmenet_variation_over_r_squared;
  
        result[0][0] = - std::sin(theta) * theta_variation
                       - radius * std::cos(theta) * theta_variation * theta_gradient[0]
                       - radius * std::sin(theta) * theta_gradient_variation[0];
  
        result[0][1] = - radius * std::cos(theta) * theta_variation * theta_gradient[1]
                       - radius * std::sin(theta) * theta_gradient_variation[1];
  
        result[0][2] = -std::cos(theta) * theta_variation;
  
        result[1][0] = 0;
        result[1][1] = 0;
        result[1][2] = 0;
  
        result[2][0] = std::cos(theta) * theta_variation
                       - radius * std::sin(theta) * theta_variation * theta_gradient[0]
                       + radius * std::cos(theta) * theta_gradient_variation[0];
  
        result[2][1] = - radius * std::sin(theta) * theta_variation * theta_gradient[1]
                       + radius * std::cos(theta) * theta_gradient_variation[1];
  
        result[2][2] = - std::sin(theta) * theta_variation;
  
        return result;
      }
  
      Tensor<2, dim+1, Number> rotation_tensor_to_transform_B_e(
          const Number angular_displacement,
          const Number radius
        ) {
        if(2 != dim) {
          throw std::logic_error("Radial deformation gradient not implemented for dim!=2");
        }
  
        Tensor<2, dim+1, Number> result;
        const Number theta = angular_displacement / radius;
  
        result[0][0] = std::cos(theta);
        result[0][1] = 0;
        result[0][2] = -std::sin(theta);
  
        result[1][0] = 0;
        result[1][1] = 1;
        result[1][2] = 0;
  
        result[2][0] = std::sin(theta);
        result[2][1] = 0;
        result[2][2] = std::cos(theta);
  
        return result;
      }
  
      Tensor<2, dim+1, Number> rotation_tensor_variation_to_transform_B_e(
          const Number angular_displacement,
          const Number angular_displacement_variation,
          const Number radius
        ) {
        if(2 != dim) {
          throw std::logic_error("Radial deformation gradient not implemented for dim!=2");
        }
  
        Tensor<2, dim+1, Number> result;
        const Number theta = angular_displacement / radius;
        const Number theta_variation = angular_displacement_variation / radius;
  
        result[0][0] = -std::sin(theta) * theta_variation;
        result[0][1] = 0;
        result[0][2] = -std::cos(theta) * theta_variation;
  
        result[1][0] = 0;
        result[1][1] = 0;
        result[1][2] = 0;
  
        result[2][0] = std::cos(theta) * theta_variation;
        result[2][1] = 0;
        result[2][2] = -std::sin(theta) * theta_variation;
  
        return result;
      }
  
      MPI_Comm mpi_communicator;
      ConditionalOStream pcout;
  
      const Number      order;
      FESystem<dim>     mech_fe;
      FE_Q<dim>         therm_fe;
      FE_DGP<dim>       mixed_var_fe;
      FESystem<dim>     mesh_motion_fe;
      MappingQ<dim>     mapping;
  
      parallel::distributed::Triangulation<dim> triangulation;
  
      DoFSystem<dim, Number> mech_dof_system;
      DoFSystem<dim, Number> therm_dof_system;
      DoFSystem<dim, Number> mixed_fe_dof_system;
  
      LBCSystem<dim, Number, dim+1> mech_lbc_system;
      LBCSystem<dim, Number, 1> therm_lbc_system;
  
      NewtonStepSystem mech_nonlinear_system;
      NewtonStepSystem therm_nonlinear_system;
  
      NewtonStepSystem mesh_motion_nonlinear_system;
      DoFSystem<dim, Number> mesh_motion_dof_system;
      LBCSystem<dim, Number, dim> mesh_motion_lbc_system;
  
      NewtonStepSystem deformation_remapping_nonlinear_system;
  
      QGauss<dim>      quadrature_formula;
      QGauss<dim-1>    face_quadrature_formula;
  
      Material<dim+1, Number>  &material;
  
      std::vector< MixedFEProjector<dim, Number> > mixed_FE_projectors;
  
      std::unordered_map<size_t, Tensor<1, dim+1, Number>> material_area_factors;
  
      Number time_increment = 1.0e-01; /*0.5e-6;*/ // [s]
      unsigned int output_rate = 1;
      Number time_since_start = 0;
  
      const Number ambient_temperature = 293.0; // [K]
      const Number rho_infty = 0.0;
  
      const unsigned int surface_boundary_id = 2;
  
      const bool COMPUTE_FORCES_PER_UNIT_AREA_IN_CURRENT_CONFIGURATION = false;
      const bool use_sigmoid_friction_law = true;
  
      Number global_lagrangian_penalty_factor = 1.0;
  
    };
  
    struct NewtonIterationDivergenceException : std::exception {
      const char *what() const _GLIBCXX_USE_NOEXCEPT override {
        return "Newton step solution diverged!\n";
      }
    };
  
  } /*namespace PlasticityLab*/
  
  #endif /* PLASTICITYLABPROG_H_ */

Annotated version of src/PlasticityLabProgDrivers.cpp

  #include <sstream>
  
  #include <deal.II/grid/tria.h>
  #include <deal.II/grid/grid_generator.h>
  #include <deal.II/grid/grid_in.h>
  #include <deal.II/grid/manifold_lib.h>
  #include <deal.II/grid/grid_tools.h>
  
  
  #include <deal.II/fe/fe_dgq.h>
  
  #include "RotationFunction.h"
  #include "ScaleZFunction.h"
  #include "ScaleComponentFunction.h"
  #include "PlasticityLabProg.h"
  
  using namespace dealii;
  using std::endl;
  
  namespace PlasticityLab {
  
    template <int dim, typename Number>
    void PlasticityLabProg<dim, Number>::run() {

make_grid_(); make_ball_in_hypershell_grid( make_cylindrical_grid( make_cylindrical_impact_grid(

      make_necking_grid(
        triangulation,
        mech_lbc_system,
        therm_lbc_system,
        3);

triangulation, mech_lbc_system, therm_lbc_system, 3); make_hook_membrane_grid(1);

      set_mesh_motion_LBCs(triangulation, mesh_motion_lbc_system);
  
      const Number total_elongation = 8.0; // mm
      const Number elongation_rate = 1.0; // [mm/s]
      const unsigned int n_steps = ceil(total_elongation/(elongation_rate*time_increment));
  
      mech_dof_system.setup_dof_system(mech_fe);
      mech_lbc_system.apply_constraints(mech_dof_system);
      mech_nonlinear_system.setup(mech_dof_system);
  
      therm_dof_system.setup_dof_system(therm_fe);
      therm_lbc_system.apply_constraints(therm_dof_system);
      therm_nonlinear_system.setup(therm_dof_system);
  

Initialize the temperature solution vector. Because we are initializing with a possibly non-zero value, we can't just assign that value to the vector in a parallel setting because the solution vector has ghost entries and so is read-only. Rather, we create a completely distributed vector, assign the value to it, and then copy that into the solution vector.

      {
        TrilinosWrappers::MPI::Vector tmp (therm_dof_system.locally_owned_dofs,
                                           mpi_communicator);
        tmp = ambient_temperature;
        therm_nonlinear_system.previous_deformation = tmp;
      }
  
      mixed_fe_dof_system.setup_dof_system(mixed_var_fe);
  
      mesh_motion_dof_system.setup_dof_system(mesh_motion_fe);
      mesh_motion_lbc_system.apply_constraints(mesh_motion_dof_system);
      mesh_motion_nonlinear_system.setup(mesh_motion_dof_system);
  
      deformation_remapping_nonlinear_system.setup(mech_dof_system);
  
      setup_material_data(triangulation, material);
      setup_material_area_factors(mesh_motion_dof_system, material_area_factors);
  
      setup_mixed_fe_projection_data(
        triangulation, mixed_FE_projectors,
        mixed_var_fe, quadrature_formula);
  
      std::vector< MixedFEProjector<dim, Number> > discontinuous_projectors;
      FE_DGQ<dim> discontinuous_fe(1);
  
      setup_mixed_fe_projection_data(
        triangulation, discontinuous_projectors,
        discontinuous_fe, quadrature_formula);
  
      TrilinosWrappers::MPI::Vector pressure(
        mixed_fe_dof_system.locally_owned_dofs,
        mpi_communicator);
  
      DoFSystem<dim, Number> discontinuous_dof_system(triangulation, mapping);
      discontinuous_dof_system.setup_dof_system(discontinuous_fe);
      TrilinosWrappers::MPI::Vector plastic_strain(
        discontinuous_dof_system.locally_owned_dofs,
        mpi_communicator);
      TrilinosWrappers::MPI::Vector von_mises_stress(
        discontinuous_dof_system.locally_owned_dofs,
        mpi_communicator);
  
      {
        TrilinosWrappers::MPI::Vector initial_velocity(
          mech_dof_system.locally_owned_dofs,
          mpi_communicator);
  
        MPI_Barrier(mpi_communicator);
  
        for(const auto initial_velocity_interpolation_handler: mech_lbc_system.initial_velocity_interpolation_handlers) {
          initial_velocity_interpolation_handler->interpolate(initial_velocity, mech_dof_system);
        }
  

Assigning the locally-owned vector into the ghosted vector performs the necessary ghost import; compress() must not be called on a vector that has ghost elements (it is read-only).

        mech_nonlinear_system.previous_time_derivative = initial_velocity;
      }
  
  
      {
        TrilinosWrappers::MPI::Vector initial_deformation(
          mech_dof_system.locally_owned_dofs,
          mpi_communicator);
  
        MPI_Barrier(mpi_communicator);
  
        for(const auto initial_deformation_interpolation_handler: mech_lbc_system.initial_deformation_interpolation_handlers) {
          initial_deformation_interpolation_handler->interpolate(initial_deformation, mech_dof_system);
        }
  
        mech_nonlinear_system.previous_deformation = initial_deformation;
      }
  
  
      for (unsigned int timeStep = 0; timeStep < n_steps + 1; ++timeStep) {
        for(const auto increment_interpolation_handler: mech_lbc_system.increment_interpolation_handlers) {
          increment_interpolation_handler->advance_time(time_increment);
        }
  
        get_plastic_strain(
          plastic_strain,
          discontinuous_dof_system.dof_handler,
          material,
          discontinuous_projectors);
  
        get_pressure(
          pressure,
          mixed_fe_dof_system.dof_handler,
          von_mises_stress,
          discontinuous_dof_system.dof_handler,
          mech_nonlinear_system,
          mech_dof_system,
          therm_nonlinear_system,
          therm_dof_system,
          material,
          mixed_FE_projectors,
          discontinuous_projectors);
  
        if(0==timeStep % output_rate) {
          pcout << "\nOutputting results..." << endl;
          DataOut<dim> data_out;
  
          data_out.add_data_vector(
            discontinuous_dof_system.dof_handler,
            plastic_strain,
            "plastic_strain");
          data_out.build_patches();
  
          data_out.add_data_vector(
            mixed_fe_dof_system.dof_handler,
            pressure,
            "pressure");
          data_out.build_patches();
  
          data_out.add_data_vector(
            discontinuous_dof_system.dof_handler,
            von_mises_stress,
            "von_mises_stress");
  
          prepare_output_results(
            data_out,
            mech_dof_system,
            mech_nonlinear_system,
            therm_dof_system,
            therm_nonlinear_system,
            mesh_motion_dof_system,
            mesh_motion_nonlinear_system);
  
          std::ostringstream oss;
          oss << "step_" << timeStep;
          const std::string output_name = oss.str();
          write_output_results(data_out, triangulation, output_name);
        }
  
        if (timeStep == n_steps) break;
  
        pcout << "\n\nStarting time step " << timeStep << ":\n\n" << endl;
  
        {
          TrilinosWrappers::MPI::Vector   step_increment(
            mech_dof_system.locally_owned_dofs,
            mpi_communicator);
  
          MPI_Barrier(mpi_communicator);
  
          for(const auto increment_interpolation_handler: mech_lbc_system.increment_interpolation_handlers) {
            increment_interpolation_handler->interpolate(step_increment, mech_dof_system);
          }
  
          mech_nonlinear_system.current_increment = step_increment;
        }
  
        solve_mesh_motion_step(
          mesh_motion_nonlinear_system,
          mesh_motion_dof_system,
          mesh_motion_lbc_system,
          mech_nonlinear_system,
          mech_dof_system,
          mixed_fe_dof_system,
          mixed_FE_projectors,
          timeStep);
  
        mesh_motion_nonlinear_system.advance_time(time_increment, rho_infty, false);

The mesh-motion deformation is the negative of the just-computed increment. 'previous_deformation' is a ghosted (read-only) vector, so this negation is delegated to the nonlinear system, which performs the arithmetic in fully-distributed temporaries.

        mesh_motion_nonlinear_system.set_previous_deformation_to_negative_current_increment();
  
        std::unordered_map<point_index_t, Tensor<2, dim+1, Number>> remapped_deformation_gradients;
        remap_material_state_variables(
          mesh_motion_nonlinear_system,
          mesh_motion_dof_system,
          mech_nonlinear_system,
          mech_dof_system,
          mixed_fe_dof_system,
          mixed_FE_projectors,
          material,
          remapped_deformation_gradients);
  
        remap_thermal_field(
          therm_nonlinear_system,
          therm_dof_system,
          mesh_motion_nonlinear_system,
          mesh_motion_dof_system);
  
        remap_mechanical_fields(
          mech_nonlinear_system,
          mech_dof_system,
          mesh_motion_nonlinear_system,
          mesh_motion_dof_system);
  
        update_material_area_factors(
          mesh_motion_nonlinear_system,
          mesh_motion_dof_system,
          material_area_factors);
  
        solve_mechanical_step(timeStep);
  
        solve_thermal_step(timeStep);
  
        solve_mechanical_step(timeStep);
  

udpate material state

        pcout << "\n\t\tassembling mechanical system updating material state..." << endl;
        assemble_mechanical_system(
          mech_nonlinear_system,
          mech_dof_system,
          mech_lbc_system,
          mesh_motion_nonlinear_system,
          mesh_motion_dof_system,
          therm_nonlinear_system,
          therm_dof_system,
          mixed_fe_dof_system,
          material,
          mixed_FE_projectors,
          false,
          true);
  
        mech_nonlinear_system.advance_time(time_increment, rho_infty, true);

Accumulate the thermal increment into the (ghosted, read-only) temperature field. The nonlinear system performs the addition in fully-distributed temporaries and assigns the result back.

        therm_nonlinear_system.add_current_increment_to_previous_deformation();
  
        therm_nonlinear_system.current_increment = 0;
  
        pcout << "Next timestep..." << std::endl;
  
      } /*for(timeStep)*/
    }
  
    template<int dim, typename Number>
    void PlasticityLabProg<dim, Number>::solve_mechanical_step(int time_step) {
  
      TrilinosWrappers::MPI::Vector total_residual;
  
      for (unsigned int NewtonStep = 0; true; NewtonStep++) {
  
        pcout << "\n\ttime step " << time_step
              << ", Newton step " << NewtonStep << "..."
              << "\n\t\tassembling mechanical system with tangents..." << endl;
  
        assemble_mechanical_system(
          mech_nonlinear_system,
          mech_dof_system,
          mech_lbc_system,
          mesh_motion_nonlinear_system,
          mesh_motion_dof_system,
          therm_nonlinear_system,
          therm_dof_system,
          mixed_fe_dof_system,
          material,
          mixed_FE_projectors,
          true);
  
        total_residual = mech_nonlinear_system.Newton_step_residual;
        total_residual.compress(VectorOperation::insert);
  
        pcout << "-------------------------------------------------------------------" << endl;
  
        pcout << "Normalized system residual: "
              << std::sqrt(total_residual.norm_sqr())
              << " ..." << endl;
        pcout << "-------------------------------------------------------------------" << endl;
  
        if (std::sqrt(total_residual.norm_sqr()) <= 1e-5) {
          break;
        }
  
        const Number old_residual = total_residual.norm_sqr();
        Number previous_residual = old_residual;
  
        pcout << "solving system..." << endl;
        try {
          solve_system(mech_dof_system, mech_nonlinear_system);
        } catch (...) {
          if (std::isnan(mech_nonlinear_system.Newton_step_solution.norm_sqr())) {
            throw;
          }
          pcout << "System solution falied. Continuing with partial solution..." << endl;
        }
        mech_dof_system.nodal_constraints.distribute(
          mech_nonlinear_system.Newton_step_solution);
        const Number solution_norm = std::sqrt(
                                       mech_nonlinear_system.Newton_step_solution.norm_sqr());
        TrilinosWrappers::MPI::Vector full_step_increment(mech_nonlinear_system.current_increment);
        TrilinosWrappers::MPI::Vector temp_locally_owned_increment(mech_dof_system.locally_owned_dofs);
  
        Number clip_factor =
          (solution_norm <= std::sqrt(old_residual)) ?
          1.0 : sqrt(old_residual) / solution_norm;
        if (clip_factor < 1.0) pcout << "clip factor: " << clip_factor << endl;
  
        pcout << "doing line search..." << endl;
        [[maybe_unused]] bool hit_line_search_limit = false;
        for (unsigned int i = 0; true/*i < 18*/; ++i) {
          const Number alpha = std::pow(0.5, static_cast<Number>(i));
          if (i > 5 && clip_factor * alpha * solution_norm < 1e-1) {
            hit_line_search_limit = true;
            break;
          }
          if (i > 0) pcout << "\tline search step " << i << "..." << endl;
  
          while (true) {
            temp_locally_owned_increment = full_step_increment;
            temp_locally_owned_increment.sadd(1, -alpha * clip_factor, mech_nonlinear_system.Newton_step_solution);
            temp_locally_owned_increment.compress(VectorOperation::insert);
            mech_nonlinear_system.current_increment = temp_locally_owned_increment;
  
            try {
              assemble_mechanical_system(
                mech_nonlinear_system,
                mech_dof_system,
                mech_lbc_system,
                mesh_motion_nonlinear_system,
                mesh_motion_dof_system,
                therm_nonlinear_system,
                therm_dof_system,
                mixed_fe_dof_system,
                material,
                mixed_FE_projectors,
                false);
            } catch (const std::runtime_error &) {
              clip_factor *= 0.125;  // using a power of 0.5; better for binary arithmetic
              pcout << "\t-------------------------------------------------------------------" << endl;
              pcout << "\tDeformation too large, causing degenerate mesh..." << endl;
              pcout << "\tupdated clip factor: " << clip_factor << endl;
              pcout << "\t-------------------------------------------------------------------" << endl;
              continue;
            }
            break;
          }
  
          total_residual = mech_nonlinear_system.Newton_step_residual;
          total_residual.compress(VectorOperation::insert);
  
          const Number current_residual = total_residual.norm_sqr();
  
          pcout << "\t-------------------------------------------------------------------" << endl;
  
          pcout << "\tNormalized system residual: "
                << std::sqrt(current_residual)
                << " ..." << endl;
          pcout << "\t-------------------------------------------------------------------" << endl;
  
          if (previous_residual < old_residual and current_residual >= previous_residual) {
            pcout << "\t---Accepting previous residual: " << std::sqrt(previous_residual)
                  << " ..." << endl;
  
            temp_locally_owned_increment = full_step_increment;
            temp_locally_owned_increment.sadd(1, -2 * alpha * clip_factor, mech_nonlinear_system.Newton_step_solution);
            temp_locally_owned_increment.compress(VectorOperation::insert);
            mech_nonlinear_system.current_increment = temp_locally_owned_increment;
            break;
          }
          previous_residual = current_residual;
        }
  
      }
    }
  
  
    template<int dim, typename Number>
    void PlasticityLabProg<dim, Number>::solve_thermal_step(int time_step) {
  
        TrilinosWrappers::MPI::Vector total_therm_residual;
        const Number starting_thermal_residual_squared_norm = 1.0;
  
        for (unsigned int NewtonStep = 0; true; NewtonStep++) {
  
          pcout << "\n\ttime step " << time_step
                << ", Newton step " << NewtonStep << "..."
                << "\n\t\tassembling thermal system with tangents..." << endl;
  
          assemble_thermal_system(
            therm_nonlinear_system,
            mech_nonlinear_system,
            therm_dof_system,
            therm_lbc_system,
            mech_dof_system,
            mixed_fe_dof_system,
            material,
            mixed_FE_projectors,
            material_area_factors,
            true);
  
          total_therm_residual = therm_nonlinear_system.Newton_step_residual;
          total_therm_residual.compress(VectorOperation::insert);
  
          pcout << "-------------------------------------------------------------------" << endl;
          pcout << "Normalized system residual (contactor): "
                << std::sqrt(total_therm_residual.norm_sqr()
                             / starting_thermal_residual_squared_norm)
                << " ..." << endl;
          pcout << "-------------------------------------------------------------------" << endl;
  
          if (std::sqrt(total_therm_residual.norm_sqr()
                        / starting_thermal_residual_squared_norm) <= 1e-6) {
            break;
          }
  
          const Number old_residual = total_therm_residual.norm_sqr();
  
          pcout << "solving system..." << endl;
  
          solve_system(therm_dof_system, therm_nonlinear_system);
  
          therm_dof_system.nodal_constraints.distribute(therm_nonlinear_system.Newton_step_solution);
          TrilinosWrappers::MPI::Vector full_step_increment(therm_nonlinear_system.current_increment);
  
          TrilinosWrappers::MPI::Vector temp_locally_owned_increment(therm_dof_system.locally_owned_dofs);
  
  
          pcout << "doing line search..." << endl;
          for (unsigned int i = 0; i < (NewtonStep > 0 ? 6 : 1); ++i) {
            const Number alpha = std::pow(0.5, static_cast<Number>(i));
  
            temp_locally_owned_increment = full_step_increment;
            temp_locally_owned_increment.sadd(1, -alpha, therm_nonlinear_system.Newton_step_solution);
            therm_dof_system.nodal_constraints.distribute(temp_locally_owned_increment);
            temp_locally_owned_increment.compress(VectorOperation::insert);
            therm_nonlinear_system.current_increment = temp_locally_owned_increment;
  
            assemble_thermal_system(
              therm_nonlinear_system,
              mech_nonlinear_system,
              therm_dof_system,
              therm_lbc_system,
              mech_dof_system,
              mixed_fe_dof_system,
              material,
              mixed_FE_projectors,
              material_area_factors,
              false);
  
            total_therm_residual = therm_nonlinear_system.Newton_step_residual;
            total_therm_residual.compress(VectorOperation::insert);
  
            const Number current_residual = total_therm_residual.norm_sqr();
            if (current_residual < old_residual)
              break;
          }
        }
        pcout << endl;
    }
  
  
    template<int dim, typename Number>
    void PlasticityLabProg<dim, Number>::solve_mesh_motion_step(
            NewtonStepSystem &mesh_motion_nonlinear_system,
            const DoFSystem<dim, Number> &mesh_motion_dof_system,
            const LBCSystem<dim, Number, dim> &mesh_motion_lbc_system,
            const NewtonStepSystem &deformation_nonlinear_system,
            const DoFSystem<dim, Number> &deformation_dof_system,
            const DoFSystem<dim, Number> &mixed_fe_dof_system,
            const std::vector< MixedFEProjector<dim, Number> > &mixed_fe_projector,
            const int time_step) {
  
      TrilinosWrappers::MPI::Vector total_residual;
  
      for (unsigned int NewtonStep = 0; true; NewtonStep++) {
  
        pcout << "\n\ttime step " << time_step
              << ", Newton step " << NewtonStep << "..."
              << "\n\t\tassembling mesh motion system with tangents..." << endl;
  
        assemble_mesh_motion_system(
          mesh_motion_nonlinear_system,
          mesh_motion_dof_system,
          mesh_motion_lbc_system,
          deformation_nonlinear_system,
          deformation_dof_system,
          mixed_fe_dof_system,
          mixed_fe_projector,
          true);
  
        total_residual = mesh_motion_nonlinear_system.Newton_step_residual;
        total_residual.compress(VectorOperation::insert);
  
        pcout << "-------------------------------------------------------------------" << endl;
  
        pcout << "Normalized system residual: "
              << std::sqrt(total_residual.norm_sqr())
              << " ..." << endl;
        pcout << "-------------------------------------------------------------------" << endl;
  
        if (std::sqrt(total_residual.norm_sqr()) <= 1e-4) {
          break;
        }
  
        const Number old_residual = total_residual.norm_sqr();
        Number previous_residual = old_residual;
  
        pcout << "solving system..." << endl;
        try {
          solve_system(mesh_motion_dof_system, mesh_motion_nonlinear_system);
        } catch (...) {
          if (std::isnan(mesh_motion_nonlinear_system.Newton_step_solution.norm_sqr())) {
            throw;
          }
          pcout << "System solution falied. Continuing with partial solution..." << endl;
        }
        mesh_motion_dof_system.nodal_constraints.distribute(mesh_motion_nonlinear_system.Newton_step_solution);
        const Number solution_norm = std::sqrt(mesh_motion_nonlinear_system.Newton_step_solution.norm_sqr());
        TrilinosWrappers::MPI::Vector full_step_increment(mesh_motion_nonlinear_system.current_increment);
        TrilinosWrappers::MPI::Vector temp_locally_owned_increment(mesh_motion_dof_system.locally_owned_dofs);
  
        Number clip_factor =
          (solution_norm <= std::sqrt(old_residual)) ? 1.0 : sqrt(old_residual) / solution_norm;
        if (clip_factor < 1.0) pcout << "clip factor: " << clip_factor << endl;
  
        pcout << "doing line search..." << endl;
        [[maybe_unused]] bool hit_line_search_limit = false;
        for (unsigned int i = 0; true/*i < 18*/; ++i) {
          const Number alpha = std::pow(0.5, static_cast<Number>(i));
          if (i > 5 && clip_factor * alpha * solution_norm < 1e-1) {
            hit_line_search_limit = true;
            break;
          }
          if (i > 0) pcout << "\tline search step " << i << "..." << endl;
  
          while (true) {
            temp_locally_owned_increment = full_step_increment;
            temp_locally_owned_increment.sadd(1, -alpha * clip_factor, mesh_motion_nonlinear_system.Newton_step_solution);
            temp_locally_owned_increment.compress(VectorOperation::insert);
            mesh_motion_nonlinear_system.current_increment = temp_locally_owned_increment;
  
            try {
              assemble_mesh_motion_system(
                mesh_motion_nonlinear_system,
                mesh_motion_dof_system,
                mesh_motion_lbc_system,
                deformation_nonlinear_system,
                deformation_dof_system,
                mixed_fe_dof_system,
                mixed_fe_projector,
                false);
            } catch (const std::runtime_error &) {
              clip_factor *= 0.125;  // using a power of 0.5; better for binary arithmetic
              pcout << "\t-------------------------------------------------------------------" << endl;
              pcout << "\tDeformation too large, causing degenerate mesh..." << endl;
              pcout << "\tupdated clip factor: " << clip_factor << endl;
              pcout << "\t-------------------------------------------------------------------" << endl;
              continue;
            }
            break;
          }
  
          total_residual = mesh_motion_nonlinear_system.Newton_step_residual;
          total_residual.compress(VectorOperation::insert);
  
          const Number current_residual = total_residual.norm_sqr();
  
          pcout << "\t-------------------------------------------------------------------" << endl;
  
          pcout << "\tNormalized system residual: "
                << std::sqrt(current_residual)
                << " ..." << endl;
          pcout << "\t-------------------------------------------------------------------" << endl;
  
          if (previous_residual < old_residual and current_residual >= previous_residual) {
            pcout << "\t---Accepting previous residual: " << std::sqrt(previous_residual)
                  << " ..." << endl;
  
            temp_locally_owned_increment = full_step_increment;
            temp_locally_owned_increment.sadd(1, -2 * alpha * clip_factor, mesh_motion_nonlinear_system.Newton_step_solution);
            temp_locally_owned_increment.compress(VectorOperation::insert);
            mesh_motion_nonlinear_system.current_increment = temp_locally_owned_increment;
            break;
          }
          previous_residual = current_residual;
        }
      }
    }
  
  
    template <int dim>
    struct RefiningTransform
    {
      RefiningTransform(
        double height,
        double refining_fraction,
        double base=0,
        size_t dimension=1) :
          height(height),
          refining_fraction(refining_fraction),
          base(base),
          dimension(dimension) {}
  
      Point<dim> operator()(const Point<dim> &p) const
      {
        Point<dim> q = p;
        if ((p[dimension]-base)/(height-base) <= 0.5) {
          q[dimension] = base + refining_fraction/0.5 * (p[dimension]-base);
        } else if ((p[dimension]-base)/(height-base) > 0.5) {
          q[dimension] = base + refining_fraction * (height - base) + (1.0 - refining_fraction) / 0.5 * (p[dimension] - 0.5 * (height + base));
        }
        return q;
      }
  
      double height;
      double refining_fraction;
      double base;
      size_t dimension;
    };
  
  
    template <int dim, typename Number>
    void PlasticityLabProg<dim, Number>::set_mesh_motion_LBCs(
        Triangulation<dim> &triangulation,
        LBCSystem<dim, Number, dim> &mesh_motion_lbc_system) {
  
      std::set<types::boundary_id> all_boundary_ids;
      for(types::boundary_id id: triangulation.get_boundary_ids()) {
        all_boundary_ids.insert(id);
      }
  
      mesh_motion_lbc_system.no_normal_flux_constraints.push_back(std::make_pair(0, all_boundary_ids));
    }
  
  
    template <int dim, typename Number>
    void PlasticityLabProg<dim, Number>::make_cylindrical_grid(
      Triangulation<dim> &triangulation,
      LBCSystem<dim, Number, dim+1> &mech_lbc_system,
      LBCSystem<dim, Number, 1> &therm_lbc_system,
      int n_initial_global_refinements) {
  
      [[maybe_unused]] const Number initial_velocity = 1.9e5; // [mm/s]
      const Number height = 2.5*25.4; // [mm]
      const Number inner_radius = 12.5; // [mm]
      const Number radius = 12.5; // [mm]
  
      const Number top_coordinate = height/2;
      const Number base_coordinate = 0.0;
  
      const unsigned int base_repetitions = std::pow(2, n_initial_global_refinements);
      const unsigned int aspect_ratio = std::ceil(0.25 * height / radius);
  
      GridGenerator::subdivided_hyper_rectangle(
        triangulation,
        std::vector<unsigned int> {base_repetitions, aspect_ratio * base_repetitions},
        Point<dim>(inner_radius, base_coordinate),
        Point<dim>(inner_radius + radius, top_coordinate),
        true
      );
  
      for (auto &cell: triangulation.active_cell_iterators()) {
        for (const auto &face : cell->face_iterators()) {
          if(face->boundary_id() == 0 || face->boundary_id() == 1) {
            if(face->center()[1] > 0.95 * height/2) {
              face->set_boundary_id(4);
            }
          }
        }
      }
  

GridTools::transform(RefiningTransform<dim>(top_coordinate, 0.3, base_coordinate), triangulation); GridTools::transform(RefiningTransform<dim>(top_coordinate, 0.4, base_coordinate), triangulation); GridTools::transform(RefiningTransform<dim>/(inner_radius, 0.35, inner_radius + radius, 0), triangulation);

for(unsigned int i=0; i<2; i++) { for (auto &cell : triangulation.active_cell_iterators()) { for (const auto &face : cell->face_iterators()) { if (face->boundary_id() == 2) { cell->set_refine_flag(); break; } } } triangulation.execute_coarsening_and_refinement(); }

      ComponentMask x_component_mask(dim+1, false);
      x_component_mask.set(0, true);
      ComponentMask y_component_mask(dim+1, false);
      y_component_mask.set(1, true);
      ComponentMask x_and_y_component_mask(dim+1, false);
      x_and_y_component_mask.set(0, true);
      x_and_y_component_mask.set(1, true);
      ComponentMask z_component_mask(dim+1, false);
      z_component_mask.set(dim-1, true);
      ComponentMask rho_component_mask(dim+1, false);
      rho_component_mask.set(dim, true);
  

std::map< types::boundary_id, const Function< dim, Number > * > base_constraint_function_map; base_constraint_function_map.insert( std::pair<types::boundary_id, Function<dim, Number>*>(2, &mech_lbc_system.zero_function)); mech_lbc_system.interpolatoryConstraintAppliers.push_back( InterpolatoryConstraintApplier<dim, Number>( base_constraint_function_map, y_component_mask));

      std::map< types::boundary_id, const Function< dim, Number > * > top_constraint_function_map;
      top_constraint_function_map.insert(
        std::pair<types::boundary_id, Function<dim, Number>*>(3, &mech_lbc_system.zero_function));
      mech_lbc_system.interpolatoryConstraintAppliers.push_back(
        InterpolatoryConstraintApplier<dim, Number>(
          top_constraint_function_map,
          x_component_mask));
  
      std::map< types::boundary_id, const Function< dim, Number > * > clamp_constraint_function_map;
      clamp_constraint_function_map.insert(
        std::pair<types::boundary_id, Function<dim, Number>*>(4, &mech_lbc_system.zero_function));
      mech_lbc_system.interpolatoryConstraintAppliers.push_back(
        InterpolatoryConstraintApplier<dim, Number>(
          clamp_constraint_function_map,
          x_component_mask));
  
      if(std::abs(inner_radius) < 1e-16) {
        std::map< types::boundary_id, const Function< dim, Number > * > axial_constraint_function_map;
        axial_constraint_function_map.insert(
          std::pair<types::boundary_id, Function<dim, Number>*>(0, &mech_lbc_system.zero_function));
        mech_lbc_system.interpolatoryConstraintAppliers.push_back(
          InterpolatoryConstraintApplier<dim, Number>(
            axial_constraint_function_map,
            x_component_mask));
  
        std::map< types::boundary_id, const Function< dim, Number > * > axial_rotation_constraint_function_map;
        axial_rotation_constraint_function_map.insert(
          std::pair<types::boundary_id, Function<dim, Number>*>(0, &mech_lbc_system.zero_function));
        mech_lbc_system.interpolatoryConstraintAppliers.push_back(
          InterpolatoryConstraintApplier<dim, Number>(
            axial_rotation_constraint_function_map,
            rho_component_mask));
      }
  

therm_lbc_system.boundaryLoadAppliers.push_back( std::pair<int,BodyForceApplier<dim,Number> >( 2, BodyForceApplier<dim,Number>(0, 22e0)));

      std::map< types::boundary_id, const Function< dim, Number > * > top_rotation_constraint_function_map;
      top_rotation_constraint_function_map.insert(
        std::pair<types::boundary_id, Function<dim, Number>*>(3, &mech_lbc_system.zero_function));
      mech_lbc_system.interpolatoryConstraintAppliers.push_back(
        InterpolatoryConstraintApplier<dim, Number>(
          top_rotation_constraint_function_map,
          rho_component_mask));
  

std::map< types::boundary_id, const Function< dim, Number > * > base_rotation_constraint_function_map; base_rotation_constraint_function_map.insert( std::pair<types::boundary_id, Function<dim, Number>*>(2, &mech_lbc_system.zero_function)); mech_lbc_system.interpolatoryConstraintAppliers.push_back( InterpolatoryConstraintApplier<dim, Number>( base_rotation_constraint_function_map, rho_component_mask));

Thermal constraints

      const Number convection_coefficient = /*17.5e-6*/ 100e-6; // [J.mm^-2.s^-1.K^-1]
      therm_lbc_system.convection_BC_appliers.push_back(
        std::pair<int, ConvectionBoundaryConditionApplier<dim, Number> >(
          0,
          ConvectionBoundaryConditionApplier<dim, Number>(
            0, convection_coefficient, ambient_temperature)));
      therm_lbc_system.convection_BC_appliers.push_back(
        std::pair<int, ConvectionBoundaryConditionApplier<dim, Number> >(
          1,
          ConvectionBoundaryConditionApplier<dim, Number>(
            0, convection_coefficient, ambient_temperature)));
      therm_lbc_system.convection_BC_appliers.push_back(
        std::pair<int, ConvectionBoundaryConditionApplier<dim, Number> >(
          2,
          ConvectionBoundaryConditionApplier<dim, Number>(
            0, convection_coefficient, ambient_temperature)));
      therm_lbc_system.convection_BC_appliers.push_back(
        std::pair<int, ConvectionBoundaryConditionApplier<dim, Number> >(
          3,
          ConvectionBoundaryConditionApplier<dim, Number>(
            0, convection_coefficient, ambient_temperature)));
  

therm_lbc_system.convection_BC_appliers.push_back( std::pair<int, ConvectionBoundaryConditionApplier<dim, Number> >( 2, ConvectionBoundaryConditionApplier<dim, Number>( 0, 3000*convection_coefficient, 1350.0 /*a little less than melting