mofem/html/tutorial_plastic_problem.html

/**

 * \file plastic.cpp

 * \example plastic.cpp

 *

 * Plasticity in 2d and 3d

 *

 */


/* The above code is a preprocessor directive in C++ that checks if the macro

"EXECUTABLE_DIMENSION" has been defined. If it has not been defined, it replaces

the " */

#ifndef EXECUTABLE_DIMENSION

#define EXECUTABLE_DIMENSION 3

#endif


// #define ADD_CONTACT


#include <MoFEM.hpp>

#include <MatrixFunction.hpp>

#include <IntegrationRules.hpp>


using namespace MoFEM;


template <int DIM> struct ElementsAndOps;


template <> struct ElementsAndOps<2> {

  using DomainEle = PipelineManager::FaceEle;

  using BoundaryEle = PipelineManager::EdgeEle;

  using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;

  static constexpr FieldSpace CONTACT_SPACE = HCURL;

};


template <> struct ElementsAndOps<3> {

  using DomainEle = PipelineManager::VolEle;

  using BoundaryEle = PipelineManager::FaceEle;

  using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;

  static constexpr FieldSpace CONTACT_SPACE = HDIV;

};


constexpr int SPACE_DIM =

    EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh

constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;


constexpr AssemblyType AT =

    (SCHUR_ASSEMBLE) ? AssemblyType::SCHUR

                     : AssemblyType::PETSC; //< selected assembly type

constexpr IntegrationType IT =

    IntegrationType::GAUSS; //< selected integration type


constexpr FieldSpace ElementsAndOps<2>::CONTACT_SPACE;

constexpr FieldSpace ElementsAndOps<3>::CONTACT_SPACE;

constexpr FieldSpace CONTACT_SPACE = ElementsAndOps<SPACE_DIM>::CONTACT_SPACE;


using EntData = EntitiesFieldData::EntData;

using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;

using DomainEleOp = DomainEle::UserDataOperator;

using BoundaryEle = ElementsAndOps<SPACE_DIM>::BoundaryEle;

using BoundaryEleOp = BoundaryEle::UserDataOperator;

using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;

using SkinPostProcEle = PostProcBrokenMeshInMoab<BoundaryEle>;

using SideEle = ElementsAndOps<SPACE_DIM>::SideEle;


#ifdef ADD_CONTACT

//! [Specialisation for assembly]


// Assemble to A matrix, by default, however, some terms are assembled only to

// preconditioning.


template <>

typename MoFEM::OpBaseImpl<AT, DomainEleOp>::MatSetValuesHook

    MoFEM::OpBaseImpl<AT, DomainEleOp>::matSetValuesHook =

        [](ForcesAndSourcesCore::UserDataOperator *op_ptr,

           const EntitiesFieldData::EntData &row_data,

           const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {

          return MatSetValues<AssemblyTypeSelector<AT>>(

              op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);

        };


template <>

typename MoFEM::OpBaseImpl<AT, BoundaryEleOp>::MatSetValuesHook

    MoFEM::OpBaseImpl<AT, BoundaryEleOp>::matSetValuesHook =

        [](ForcesAndSourcesCore::UserDataOperator *op_ptr,

           const EntitiesFieldData::EntData &row_data,

           const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {

          return MatSetValues<AssemblyTypeSelector<AT>>(

              op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);

        };


/**

 * @brief Element used to specialise assembly

 *

 */

struct BoundaryEleOpStab : public BoundaryEleOp {

  using BoundaryEleOp::BoundaryEleOp;

};


/**

 * @brief Specialise assembly for Stabilised matrix

 *

 * @tparam

 */

template <>

typename MoFEM::OpBaseImpl<AT, BoundaryEleOpStab>::MatSetValuesHook

    MoFEM::OpBaseImpl<AT, BoundaryEleOpStab>::matSetValuesHook =

        [](ForcesAndSourcesCore::UserDataOperator *op_ptr,

           const EntitiesFieldData::EntData &row_data,

           const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {

          return MatSetValues<AssemblyTypeSelector<AT>>(

              op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);

        };

//! [Specialisation for assembly]

#endif // ADD_CONTACT


inline double iso_hardening_exp(double tau, double b_iso) {

  return std::exp(

      std::max(static_cast<double>(std::numeric_limits<float>::min_exponent10),

               -b_iso * tau));

}


/**

 * Isotropic hardening

 */

inline double iso_hardening(double tau, double H, double Qinf, double b_iso,

                            double sigmaY) {

  return H * tau + Qinf * (1. - iso_hardening_exp(tau, b_iso)) + sigmaY;

}


inline double iso_hardening_dtau(double tau, double H, double Qinf,

                                 double b_iso) {

  auto r = [&](auto tau) {

    return H + Qinf * b_iso * iso_hardening_exp(tau, b_iso);

  };

  constexpr double eps = 1e-12;

  return std::max(r(tau), eps * r(0));

}


/**

 * Kinematic hardening

 */

template <typename T, int DIM>

inline auto

kinematic_hardening(FTensor::Tensor2_symmetric<T, DIM> &t_plastic_strain,

                    double C1_k) {

  FTensor::Index<'i', DIM> i;

  FTensor::Index<'j', DIM> j;

  FTensor::Tensor2_symmetric<double, DIM> t_alpha;

  if (C1_k < std::numeric_limits<double>::epsilon()) {

    t_alpha(i, j) = 0;

    return t_alpha;

  }

  t_alpha(i, j) = C1_k * t_plastic_strain(i, j);

  return t_alpha;

}


template <int DIM>

inline auto kinematic_hardening_dplastic_strain(double C1_k) {

  FTensor::Index<'i', DIM> i;

  FTensor::Index<'j', DIM> j;

  FTensor::Index<'k', DIM> k;

  FTensor::Index<'l', DIM> l;

  FTensor::Ddg<double, DIM, DIM> t_diff;

  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();

  t_diff(i, j, k, l) = C1_k * (t_kd(i, k) ^ t_kd(j, l)) / 4.;

  return t_diff;

}


PetscBool is_large_strains = PETSC_TRUE; ///< Large strains

PetscBool set_timer = PETSC_FALSE;       ///< Set timer


double scale = 1.;


double young_modulus = 206913; ///< Young modulus

double poisson_ratio = 0.29;   ///< Poisson ratio

double sigmaY = 450;           ///< Yield stress

double H = 129;                ///< Hardening

double visH = 0;               ///< Viscous hardening

double zeta = 5e-2;            ///< Viscous hardening

double Qinf = 265;             ///< Saturation yield stress

double b_iso = 16.93;          ///< Saturation exponent

double C1_k = 0;               ///< Kinematic hardening


double cn0 = 1;

double cn1 = 1;


int order = 2;             ///< Order displacement

int tau_order = order - 2; ///< Order of tau files

int ep_order = order - 1;  ///< Order of ep files

int geom_order = 2;        ///< Order if fixed.


PetscBool is_quasi_static = PETSC_TRUE;

double rho = 0.0;

double alpha_damping = 0;


#include <HenckyOps.hpp>

#include <PlasticOps.hpp>

#include <PlasticNaturalBCs.hpp>


#ifdef ADD_CONTACT

#ifdef PYTHON_SFD

#include <boost/python.hpp>

#include <boost/python/def.hpp>

namespace bp = boost::python;

#endif


namespace ContactOps {


double cn_contact = 0.1;


}; // namespace ContactOps


#include <ContactOps.hpp>

#endif // ADD_CONTACT


using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<AT>::LinearForm<IT>;

using OpDomainRhsBCs =

    DomainRhsBCs::OpFlux<PlasticOps::DomainBCs, 1, SPACE_DIM>;

using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::LinearForm<IT>;

using OpBoundaryRhsBCs =

    BoundaryRhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;

using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::BiLinearForm<IT>;

using OpBoundaryLhsBCs =

    BoundaryLhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;


using namespace PlasticOps;

using namespace HenckyOps;

struct Example {


  Example(MoFEM::Interface &m_field) : mField(m_field) {}


  MoFEMErrorCode runProblem();


private:

  MoFEM::Interface &mField;


  MoFEMErrorCode setupProblem();

  MoFEMErrorCode createCommonData();

  MoFEMErrorCode bC();

  MoFEMErrorCode OPs();

  MoFEMErrorCode tsSolve();


  boost::shared_ptr<DomainEle> reactionFe;


  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;

  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;

  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;


  struct ScaledTimeScale : public MoFEM::TimeScale {

    using MoFEM::TimeScale::TimeScale;

    double getScale(const double time) {

      return scale * MoFEM::TimeScale::getScale(time);

    };

  };


#ifdef ADD_CONTACT

#ifdef PYTHON_SFD

  boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;

#endif

#endif // ADD_CONTACT

};


//! [Run problem]

MoFEMErrorCode Example::runProblem() {

  MoFEMFunctionBegin;

  CHKERR createCommonData();

  CHKERR setupProblem();

  CHKERR bC();

  CHKERR OPs();

  CHKERR tsSolve();

  MoFEMFunctionReturn(0);

}

//! [Run problem]


//! [Set up problem]

MoFEMErrorCode Example::setupProblem() {

  MoFEMFunctionBegin;

  Simple *simple = mField.getInterface<Simple>();


  Range domain_ents;

  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,

                                                     true);

  auto get_ents_by_dim = [&](const auto dim) {

    if (dim == SPACE_DIM) {

      return domain_ents;

    } else {

      Range ents;

      if (dim == 0)

        CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);

      else

        CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);

      return ents;

    }

  };


  auto get_base = [&]() {

    auto domain_ents = get_ents_by_dim(SPACE_DIM);

    if (domain_ents.empty())

      CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");

    const auto type = type_from_handle(domain_ents[0]);

    switch (type) {

    case MBQUAD:

      return DEMKOWICZ_JACOBI_BASE;

    case MBHEX:

      return DEMKOWICZ_JACOBI_BASE;

    case MBTRI:

      return AINSWORTH_LEGENDRE_BASE;

    case MBTET:

      return AINSWORTH_LEGENDRE_BASE;

    default:

      CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");

    }

    return NOBASE;

  };


  const auto base = get_base();

  MOFEM_LOG("PLASTICITY", Sev::inform)

      << "Base " << ApproximationBaseNames[base];


  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);

  CHKERR simple->addDomainField("EP", L2, base, size_symm);

  CHKERR simple->addDomainField("TAU", L2, base, 1);

  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);


  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);


  PetscBool order_edge = PETSC_FALSE;

  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_edge", &order_edge,

                             PETSC_NULL);

  PetscBool order_face = PETSC_FALSE;

  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_face", &order_face,

                             PETSC_NULL);

  PetscBool order_volume = PETSC_FALSE;

  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_volume", &order_volume,

                             PETSC_NULL);


  if (order_edge || order_face || order_volume) {


    MOFEM_LOG("PLASTICITY", Sev::inform) << "Order edge " << order_edge

        ? "true"

        : "false";

    MOFEM_LOG("PLASTICITY", Sev::inform) << "Order face " << order_face

        ? "true"

        : "false";

    MOFEM_LOG("PLASTICITY", Sev::inform) << "Order volume " << order_volume

        ? "true"

        : "false";


    auto ents = get_ents_by_dim(0);

    if (order_edge)

      ents.merge(get_ents_by_dim(1));

    if (order_face)

      ents.merge(get_ents_by_dim(2));

    if (order_volume)

      ents.merge(get_ents_by_dim(3));

    CHKERR simple->setFieldOrder("U", order, &ents);

  } else {

    CHKERR simple->setFieldOrder("U", order);

  }

  CHKERR simple->setFieldOrder("EP", ep_order);

  CHKERR simple->setFieldOrder("TAU", tau_order);


  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);


#ifdef ADD_CONTACT

  CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,

                                SPACE_DIM);

  CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,

                                  SPACE_DIM);


  auto get_skin = [&]() {

    Range body_ents;

    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);

    Skinner skin(&mField.get_moab());

    Range skin_ents;

    CHKERR skin.find_skin(0, body_ents, false, skin_ents);

    return skin_ents;

  };


  auto filter_blocks = [&](auto skin) {

    bool is_contact_block = false;

    Range contact_range;

    for (auto m :

         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(


             (boost::format("%s(.*)") % "CONTACT").str()


                 ))


    ) {

      is_contact_block =

          true; ///< bloks interation is collectibe, so that is set irrespective

                ///< if there are enerities in given rank or not in the block

      MOFEM_LOG("CONTACT", Sev::inform)

          << "Find contact block set:  " << m->getName();

      auto meshset = m->getMeshset();

      Range contact_meshset_range;

      CHKERR mField.get_moab().get_entities_by_dimension(

          meshset, SPACE_DIM - 1, contact_meshset_range, true);


      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(

          contact_meshset_range);

      contact_range.merge(contact_meshset_range);

    }

    if (is_contact_block) {

      MOFEM_LOG("SYNC", Sev::inform)

          << "Nb entities in contact surface: " << contact_range.size();

      MOFEM_LOG_SYNCHRONISE(mField.get_comm());

      skin = intersect(skin, contact_range);

    }

    return skin;

  };


  auto filter_true_skin = [&](auto skin) {

    Range boundary_ents;

    ParallelComm *pcomm =

        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);

    CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,

                                 PSTATUS_NOT, -1, &boundary_ents);

    return boundary_ents;

  };


  auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));

  CHKERR simple->setFieldOrder("SIGMA", 0);

  CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);

#endif


  CHKERR simple->setUp();

  CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");


  auto project_ho_geometry = [&]() {

    Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");

    return mField.loop_dofs("GEOMETRY", ent_method);

  };

  CHKERR project_ho_geometry();


  MoFEMFunctionReturn(0);

}

//! [Set up problem]


//! [Create common data]

MoFEMErrorCode Example::createCommonData() {

  MoFEMFunctionBegin;


  auto get_command_line_parameters = [&]() {

    MoFEMFunctionBegin;


    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",

                                 &young_modulus, PETSC_NULL);

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",

                                 &poisson_ratio, PETSC_NULL);

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,

                                 PETSC_NULL);

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,

                                 PETSC_NULL);

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn0", &cn0, PETSC_NULL);

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn1", &cn1, PETSC_NULL);

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-zeta", &zeta, PETSC_NULL);

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-C1_k", &C1_k, PETSC_NULL);

    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strains",

                               &is_large_strains, PETSC_NULL);

    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-set_timer", &set_timer,

                               PETSC_NULL);


    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);

    PetscBool tau_order_is_set; ///< true if tau order is set

    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-tau_order", &tau_order,

                              &tau_order_is_set);

    PetscBool ep_order_is_set; ///< true if tau order is set

    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-ep_order", &ep_order,

                              &ep_order_is_set);

    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,

                              PETSC_NULL);


    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_damping",

                                 &alpha_damping, PETSC_NULL);


    MOFEM_LOG("PLASTICITY", Sev::inform) << "Young modulus " << young_modulus;

    MOFEM_LOG("PLASTICITY", Sev::inform) << "Poisson ratio " << poisson_ratio;

    MOFEM_LOG("PLASTICITY", Sev::inform) << "Yield stress " << sigmaY;

    MOFEM_LOG("PLASTICITY", Sev::inform) << "Hardening " << H;

    MOFEM_LOG("PLASTICITY", Sev::inform) << "Viscous hardening " << visH;

    MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation yield stress " << Qinf;

    MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation exponent " << b_iso;

    MOFEM_LOG("PLASTICITY", Sev::inform) << "Kinematic hardening " << C1_k;

    MOFEM_LOG("PLASTICITY", Sev::inform) << "cn0 " << cn0;

    MOFEM_LOG("PLASTICITY", Sev::inform) << "cn1 " << cn1;

    MOFEM_LOG("PLASTICITY", Sev::inform) << "zeta " << zeta;


    if (tau_order_is_set == PETSC_FALSE)

      tau_order = order - 2;

    if (ep_order_is_set == PETSC_FALSE)

      ep_order = order - 1;


    MOFEM_LOG("PLASTICITY", Sev::inform) << "Approximation order " << order;

    MOFEM_LOG("PLASTICITY", Sev::inform)

        << "Ep approximation order " << ep_order;

    MOFEM_LOG("PLASTICITY", Sev::inform)

        << "Tau approximation order " << tau_order;

    MOFEM_LOG("PLASTICITY", Sev::inform)

        << "Geometry approximation order " << geom_order;


    MOFEM_LOG("PLASTICITY", Sev::inform) << "Density " << rho;

    MOFEM_LOG("PLASTICITY", Sev::inform) << "alpha_damping " << alpha_damping;


    PetscBool is_scale = PETSC_TRUE;

    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,

                               PETSC_NULL);

    if (is_scale) {

      scale /= young_modulus;

    }


    MOFEM_LOG("PLASTICITY", Sev::inform) << "Scale " << scale;


#ifdef ADD_CONTACT

    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn_contact",

                                 &ContactOps::cn_contact, PETSC_NULL);

    MOFEM_LOG("CONTACT", Sev::inform)

        << "cn_contact " << ContactOps::cn_contact;

#endif // ADD_CONTACT


    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-quasi_static",

                               &is_quasi_static, PETSC_NULL);

    MOFEM_LOG("PLASTICITY", Sev::inform)

        << "Is quasi static: " << (is_quasi_static ? "true" : "false");


    MoFEMFunctionReturn(0);

  };


  CHKERR get_command_line_parameters();


#ifdef ADD_CONTACT

#ifdef PYTHON_SFD

  sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();

  CHKERR sdfPythonPtr->sdfInit("sdf.py");

  ContactOps::sdfPythonWeakPtr = sdfPythonPtr;

#endif

#endif // ADD_CONTACT


  MoFEMFunctionReturn(0);

}

//! [Create common data]


//! [Boundary condition]

MoFEMErrorCode Example::bC() {

  MoFEMFunctionBegin;


  auto simple = mField.getInterface<Simple>();

  auto bc_mng = mField.getInterface<BcManager>();

  auto prb_mng = mField.getInterface<ProblemsManager>();


  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",

                                           "U", 0, 0);

  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",

                                           "U", 1, 1);

  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",

                                           "U", 2, 2);

  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),

                                           "REMOVE_ALL", "U", 0, 3);


#ifdef ADD_CONTACT

  for (auto b : {"FIX_X", "REMOVE_X"})

    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,

                                             "SIGMA", 0, 0, false, true);

  for (auto b : {"FIX_Y", "REMOVE_Y"})

    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,

                                             "SIGMA", 1, 1, false, true);

  for (auto b : {"FIX_Z", "REMOVE_Z"})

    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,

                                             "SIGMA", 2, 2, false, true);

  for (auto b : {"FIX_ALL", "REMOVE_ALL"})

    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,

                                             "SIGMA", 0, 3, false, true);

  CHKERR bc_mng->removeBlockDOFsOnEntities(

      simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);

#endif


  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(

      simple->getProblemName(), "U");


  auto &bc_map = bc_mng->getBcMapByBlockName();

  for (auto bc : bc_map)

    MOFEM_LOG("PLASTICITY", Sev::verbose) << "Marker " << bc.first;


  MoFEMFunctionReturn(0);

}

//! [Boundary condition]


//! [Push operators to pipeline]

MoFEMErrorCode Example::OPs() {

  MoFEMFunctionBegin;

  auto pip_mng = mField.getInterface<PipelineManager>();

  auto simple = mField.getInterface<Simple>();

  auto bc_mng = mField.getInterface<BcManager>();


  auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };


  auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order - 1; };


  auto add_boundary_ops_lhs_mechanical = [&](auto &pip) {

    MoFEMFunctionBegin;


    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},

                                                              "GEOMETRY");

    pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());


    // Add Natural BCs to LHS

    CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(

        pip, mField, "U", Sev::inform);


#ifdef ADD_CONTACT

    CHKERR

    ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(

        pip, "SIGMA", "U");

    CHKERR

    ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(

        mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",

        vol_rule);

#endif // ADD_CONTACT


    MoFEMFunctionReturn(0);

  };


  auto add_boundary_ops_rhs_mechanical = [&](auto &pip) {

    MoFEMFunctionBegin;


    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},

                                                              "GEOMETRY");

    pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());


    // Add Natural BCs to RHS

    CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(

        pip, mField, "U", {boost::make_shared<ScaledTimeScale>()}, Sev::inform);


#ifdef ADD_CONTACT

    CHKERR ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, IT, BoundaryEleOp>(

        pip, "SIGMA", "U");

#endif // ADD_CONTACT


    MoFEMFunctionReturn(0);

  };


  auto add_domain_ops_lhs = [this](auto &pip) {

    MoFEMFunctionBegin;

    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},

                                                          "GEOMETRY");


    if (is_quasi_static == PETSC_FALSE) {


      //! [Only used for dynamics]

      using OpMass = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<

          GAUSS>::OpMass<1, SPACE_DIM>;

      //! [Only used for dynamics]


      auto get_inertia_and_mass_damping = [this](const double, const double,

                                                 const double) {

        auto *pip = mField.getInterface<PipelineManager>();

        auto &fe_domain_lhs = pip->getDomainLhsFE();

        return (rho / scale) * fe_domain_lhs->ts_aa +

               (alpha_damping / scale) * fe_domain_lhs->ts_a;

      };

      pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));

    }


    CHKERR PlasticOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(

        mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");


    MoFEMFunctionReturn(0);

  };


  auto add_domain_ops_rhs = [this](auto &pip) {

    MoFEMFunctionBegin;


    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},

                                                          "GEOMETRY");


    CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(

        pip, mField, "U",

        {boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},

        Sev::inform);


    // only in case of dynamics

    if (is_quasi_static == PETSC_FALSE) {


      //! [Only used for dynamics]

      using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<

          AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;

      //! [Only used for dynamics]


      auto mat_acceleration = boost::make_shared<MatrixDouble>();

      pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(

          "U", mat_acceleration));

      pip.push_back(

          new OpInertiaForce("U", mat_acceleration, [](double, double, double) {

            return rho / scale;

          }));

      if (alpha_damping > 0) {

        auto mat_velocity = boost::make_shared<MatrixDouble>();

        pip.push_back(

            new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));

        pip.push_back(

            new OpInertiaForce("U", mat_velocity, [](double, double, double) {

              return alpha_damping / scale;

            }));

      }

    }


    CHKERR PlasticOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(

        mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");


#ifdef ADD_CONTACT

    CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(

        pip, "SIGMA", "U");

#endif // ADD_CONTACT


    MoFEMFunctionReturn(0);

  };


  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());

  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());


  // Boundary

  CHKERR add_boundary_ops_lhs_mechanical(pip_mng->getOpBoundaryLhsPipeline());

  CHKERR add_boundary_ops_rhs_mechanical(pip_mng->getOpBoundaryRhsPipeline());


  CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);

  CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);


  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);

  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);


  auto create_reaction_pipeline = [&](auto &pip) {

    MoFEMFunctionBegin;

    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},

                                                          "GEOMETRY");

    CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(

        mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");

    MoFEMFunctionReturn(0);

  };


  reactionFe = boost::make_shared<DomainEle>(mField);

  reactionFe->getRuleHook = vol_rule;

  CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());

  reactionFe->postProcessHook =

      EssentialPreProcReaction<DisplacementCubitBcData>(mField, reactionFe);


  MoFEMFunctionReturn(0);

}

//! [Push operators to pipeline]


//! [Solve]

struct SetUpSchur {


  /**

   * @brief Create data structure for handling Schur complement

   *

   * @param m_field

   * @param sub_dm  Schur complement sub dm

   * @param field_split_it IS of Schur block

   * @param ao_map AO map from sub dm to main problem

   * @return boost::shared_ptr<SetUpSchur>

   */

  static boost::shared_ptr<SetUpSchur> createSetUpSchur(


      MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,

      SmartPetscObj<IS> field_split_it, SmartPetscObj<AO> ao_map


  );

  virtual MoFEMErrorCode setUp(TS solver) = 0;


protected:

  SetUpSchur() = default;

};


MoFEMErrorCode Example::tsSolve() {

  MoFEMFunctionBegin;


  Simple *simple = mField.getInterface<Simple>();

  PipelineManager *pip_mng = mField.getInterface<PipelineManager>();

  ISManager *is_manager = mField.getInterface<ISManager>();


  auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());


  auto set_section_monitor = [&](auto solver) {

    MoFEMFunctionBegin;

    SNES snes;

    CHKERR TSGetSNES(solver, &snes);

    CHKERR SNESMonitorSet(snes,

                          (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,

                                             void *))MoFEMSNESMonitorFields,

                          (void *)(snes_ctx_ptr.get()), nullptr);

    MoFEMFunctionReturn(0);

  };


  auto create_post_process_elements = [&]() {

    auto pp_fe = boost::make_shared<PostProcEle>(mField);

    auto &pip = pp_fe->getOpPtrVector();


    auto push_vol_ops = [this](auto &pip) {

      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},

                                                            "GEOMETRY");


      auto [common_plastic_ptr, common_henky_ptr] =

          PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(

              mField, "MAT_PLASTIC", pip, "U", "EP", "TAU", 1., Sev::inform);


      if (common_henky_ptr) {

        if (common_plastic_ptr->mGradPtr != common_henky_ptr->matGradPtr)

          CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");

      }


      return std::make_pair(common_plastic_ptr, common_henky_ptr);

    };


    auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {

      MoFEMFunctionBegin;


      auto &pip = pp_fe->getOpPtrVector();


      auto [common_plastic_ptr, common_henky_ptr] = p;


      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;


      auto x_ptr = boost::make_shared<MatrixDouble>();

      pip.push_back(

          new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));

      auto u_ptr = boost::make_shared<MatrixDouble>();

      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));


      if (is_large_strains) {


        pip.push_back(


            new OpPPMap(


                pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),


                {{"PLASTIC_SURFACE",

                  common_plastic_ptr->getPlasticSurfacePtr()},

                 {"PLASTIC_MULTIPLIER",

                  common_plastic_ptr->getPlasticTauPtr()}},


                {{"U", u_ptr}, {"GEOMETRY", x_ptr}},


                {{"GRAD", common_henky_ptr->matGradPtr},

                 {"FIRST_PIOLA", common_henky_ptr->getMatFirstPiolaStress()}},


                {{"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},

                 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}


                )


        );


      } else {


        pip.push_back(


            new OpPPMap(


                pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),


                {{"PLASTIC_SURFACE",

                  common_plastic_ptr->getPlasticSurfacePtr()},

                 {"PLASTIC_MULTIPLIER",

                  common_plastic_ptr->getPlasticTauPtr()}},


                {{"U", u_ptr}, {"GEOMETRY", x_ptr}},


                {},


                {{"STRAIN", common_plastic_ptr->mStrainPtr},

                 {"STRESS", common_plastic_ptr->mStressPtr},

                 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},

                 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}


                )


        );

      }


      MoFEMFunctionReturn(0);

    };


    auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {

      PetscBool post_proc_vol = PETSC_FALSE;

      CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol",

                                 &post_proc_vol, PETSC_NULL);

      if (post_proc_vol == PETSC_FALSE)

        return boost::shared_ptr<PostProcEle>();

      auto pp_fe = boost::make_shared<PostProcEle>(mField);

      CHK_MOAB_THROW(

          push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),

          "push_vol_post_proc_ops");

      return pp_fe;

    };


    auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {

      PetscBool post_proc_skin = PETSC_TRUE;

      CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",

                                 &post_proc_skin, PETSC_NULL);

      if (post_proc_skin == PETSC_FALSE)

        return boost::shared_ptr<SkinPostProcEle>();


      auto simple = mField.getInterface<Simple>();

      auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);

      auto op_side = new OpLoopSide<SideEle>(

          mField, simple->getDomainFEName(), SPACE_DIM, Sev::verbose);

      pp_fe->getOpPtrVector().push_back(op_side);

      CHK_MOAB_THROW(push_vol_post_proc_ops(

                         pp_fe, push_vol_ops(op_side->getOpPtrVector())),

                     "push_vol_post_proc_ops");

      return pp_fe;

    };


    return std::make_pair(vol_post_proc(), skin_post_proc());

  };


  auto scatter_create = [&](auto D, auto coeff) {

    SmartPetscObj<IS> is;

    CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),

                                                   ROW, "U", coeff, coeff, is);

    int loc_size;

    CHKERR ISGetLocalSize(is, &loc_size);

    Vec v;

    CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);

    VecScatter scatter;

    CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);

    return std::make_tuple(SmartPetscObj<Vec>(v),

                           SmartPetscObj<VecScatter>(scatter));

  };


  auto set_time_monitor = [&](auto dm, auto solver) {

    MoFEMFunctionBegin;

    boost::shared_ptr<Monitor> monitor_ptr(

        new Monitor(dm, create_post_process_elements(), reactionFe, uXScatter,

                    uYScatter, uZScatter));

    boost::shared_ptr<ForcesAndSourcesCore> null;

    CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),

                               monitor_ptr, null, null);

    MoFEMFunctionReturn(0);

  };


  auto set_schur_pc = [&](auto solver,

                          boost::shared_ptr<SetUpSchur> &schur_ptr) {

    MoFEMFunctionBeginHot;


    auto bc_mng = mField.getInterface<BcManager>();

    auto name_prb = simple->getProblemName();


    // create sub dm for Schur complement

    auto create_sub_u_dm = [&](SmartPetscObj<DM> base_dm,

                               SmartPetscObj<DM> &dm_sub) {

      MoFEMFunctionBegin;

      dm_sub = createDM(mField.get_comm(), "DMMOFEM");

      CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SUB_U");

      CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);

      CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());

      CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());

      for (auto f : {"U"}) {

        CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);

        CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);

      }

      CHKERR DMSetUp(dm_sub);


      MoFEMFunctionReturn(0);

    };


    // Create nested (sub BC) Schur DM

    if constexpr (AT == AssemblyType::SCHUR) {

      SmartPetscObj<IS> is_epp;

      CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(

          simple->getProblemName(), ROW, "EP", 0, MAX_DOFS_ON_ENTITY, is_epp);

      SmartPetscObj<IS> is_tau;

      CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(

          simple->getProblemName(), ROW, "TAU", 0, MAX_DOFS_ON_ENTITY, is_tau);


      IS is_union_raw;

      CHKERR ISExpand(is_epp, is_tau, &is_union_raw);

      SmartPetscObj<IS> is_union(is_union_raw);


#ifdef ADD_CONTACT

      auto add_sigma_to_is = [&](auto is_union) {

        SmartPetscObj<IS> is_union_sigma;

        auto add_sigma_to_is_impl = [&]() {

          MoFEMFunctionBegin;

          SmartPetscObj<IS> is_sigma;

          CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(

              simple->getProblemName(), ROW, "SIGMA", 0, MAX_DOFS_ON_ENTITY,

              is_sigma);

          IS is_union_raw_sigma;

          CHKERR ISExpand(is_union, is_sigma, &is_union_raw_sigma);

          is_union_sigma = SmartPetscObj<IS>(is_union_raw_sigma);

          MoFEMFunctionReturn(0);

        };

        CHK_THROW_MESSAGE(add_sigma_to_is_impl(), "Can not add sigma to IS");

        return is_union_sigma;

      };

      is_union = add_sigma_to_is(is_union);

#endif // ADD_CONTACT


      SmartPetscObj<DM> dm_u_sub;

      CHKERR create_sub_u_dm(simple->getDM(), dm_u_sub);


      // Indices has to be map fro very to level, while assembling Schur

      // complement.

      auto is_up = getDMSubData(dm_u_sub)->getSmartRowIs();

      auto ao_up = createAOMappingIS(is_up, PETSC_NULL);

      schur_ptr =

          SetUpSchur::createSetUpSchur(mField, dm_u_sub, is_union, ao_up);

      CHKERR schur_ptr->setUp(solver);

    }


    MoFEMFunctionReturnHot(0);

  };


  auto dm = simple->getDM();

  auto D = createDMVector(dm);

  auto DD = vectorDuplicate(D);

  uXScatter = scatter_create(D, 0);

  uYScatter = scatter_create(D, 1);

  if constexpr (SPACE_DIM == 3)

    uZScatter = scatter_create(D, 2);


  auto create_solver = [pip_mng]() {

    if (is_quasi_static == PETSC_TRUE)

      return pip_mng->createTSIM();

    else

      return pip_mng->createTSIM2();

  };


  auto solver = create_solver();


  auto active_pre_lhs = []() {

    MoFEMFunctionBegin;

    std::fill(PlasticOps::CommonData::activityData.begin(),

              PlasticOps::CommonData::activityData.end(), 0);

    MoFEMFunctionReturn(0);

  };


  auto active_post_lhs = [&]() {

    MoFEMFunctionBegin;

    auto get_iter = [&]() {

      SNES snes;

      CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");

      int iter;

      CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),

                        "Can not get iter");

      return iter;

    };


    auto iter = get_iter();

    if (iter >= 0) {


      std::array<int, 5> activity_data;

      std::fill(activity_data.begin(), activity_data.end(), 0);

      MPI_Allreduce(PlasticOps::CommonData::activityData.data(),

                    activity_data.data(), activity_data.size(), MPI_INT,

                    MPI_SUM, mField.get_comm());


      int &active_points = activity_data[0];

      int &avtive_full_elems = activity_data[1];

      int &avtive_elems = activity_data[2];

      int &nb_points = activity_data[3];

      int &nb_elements = activity_data[4];


      if (nb_points) {


        double proc_nb_points =

            100 * static_cast<double>(active_points) / nb_points;

        double proc_nb_active =

            100 * static_cast<double>(avtive_elems) / nb_elements;

        double proc_nb_full_active = 100;

        if (avtive_elems)

          proc_nb_full_active =

              100 * static_cast<double>(avtive_full_elems) / avtive_elems;


        MOFEM_LOG_C("PLASTICITY", Sev::inform,

                    "Iter %d nb pts %d nb avtive pts %d (%3.3f\%) nb active "

                    "elements %d "

                    "(%3.3f\%) nb full active elems %d (%3.3f\%)",

                    iter, nb_points, active_points, proc_nb_points,

                    avtive_elems, proc_nb_active, avtive_full_elems,

                    proc_nb_full_active, iter);

      }

    }


    MoFEMFunctionReturn(0);

  };


  auto add_active_dofs_elem = [&](auto dm) {

    MoFEMFunctionBegin;

    auto fe_pre_proc = boost::make_shared<FEMethod>();

    fe_pre_proc->preProcessHook = active_pre_lhs;

    auto fe_post_proc = boost::make_shared<FEMethod>();

    fe_post_proc->postProcessHook = active_post_lhs;

    auto ts_ctx_ptr = getDMTsCtx(dm);

    ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);

    ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);

    MoFEMFunctionReturn(0);

  };


  auto set_essential_bc = [&](auto dm, auto solver) {

    MoFEMFunctionBegin;

    // This is low level pushing finite elements (pipelines) to solver


    auto pre_proc_ptr = boost::make_shared<FEMethod>();

    auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();

    auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();


    // Add boundary condition scaling

    auto disp_time_scale = boost::make_shared<TimeScale>();


    auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {

      EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,

                                                     {disp_time_scale}, false);

      return hook;

    };

    pre_proc_ptr->preProcessHook = get_bc_hook_rhs();


    auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {

      MoFEMFunctionBegin;

      CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(

          mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();

      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(

          mField, post_proc_rhs_ptr, 1.)();

      MoFEMFunctionReturn(0);

    };

    auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {

      return EssentialPostProcLhs<DisplacementCubitBcData>(

          mField, post_proc_lhs_ptr, 1.);

    };

    post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;


    auto ts_ctx_ptr = getDMTsCtx(dm);

    ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);

    ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);

    ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);


    SNES snes;

    CHKERR TSGetSNES(solver, &snes);

    KSP ksp;

    CHKERR SNESGetKSP(snes, &ksp);

    PC pc;

    CHKERR KSPGetPC(ksp, &pc);

    PetscBool is_pcfs = PETSC_FALSE;

    PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);


    if (is_pcfs == PETSC_FALSE) {

      post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();

      ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);

    }

    MoFEMFunctionReturn(0);

  };


  if (is_quasi_static == PETSC_TRUE) {

    CHKERR TSSetSolution(solver, D);

  } else {

    CHKERR TS2SetSolution(solver, D, DD);

  }


  CHKERR set_section_monitor(solver);

  CHKERR set_time_monitor(dm, solver);

  CHKERR TSSetFromOptions(solver);

  CHKERR TSSetUp(solver);


  CHKERR add_active_dofs_elem(dm);

  boost::shared_ptr<SetUpSchur> schur_ptr;

  CHKERR set_schur_pc(solver, schur_ptr);

  CHKERR set_essential_bc(dm, solver);


  MOFEM_LOG_CHANNEL("TIMER");

  MOFEM_LOG_TAG("TIMER", "timer");

  if (set_timer)

    BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());

  MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp";

  CHKERR TSSetUp(solver);

  MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp <= done";

  MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve";

  CHKERR TSSolve(solver, NULL);

  MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve <= done";


  MoFEMFunctionReturn(0);

}

//! [Solve]


static char help[] = "...\n\n";


int main(int argc, char *argv[]) {


#ifdef ADD_CONTACT

#ifdef PYTHON_SFD

  Py_Initialize();

#endif

#endif // ADD_CONTACT


  // Initialisation of MoFEM/PETSc and MOAB data structures

  const char param_file[] = "param_file.petsc";

  MoFEM::Core::Initialize(&argc, &argv, param_file, help);


  // Add logging channel for example

  auto core_log = logging::core::get();

  core_log->add_sink(

      LogManager::createSink(LogManager::getStrmWorld(), "PLASTICITY"));

  core_log->add_sink(

      LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));

  LogManager::setLog("PLASTICITY");

  MOFEM_LOG_TAG("PLASTICITY", "Plasticity");


#ifdef ADD_CONTACT

  core_log->add_sink(

      LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));

  LogManager::setLog("CONTACT");

  MOFEM_LOG_TAG("CONTACT", "Contact");

#endif // ADD_CONTACT


  try {


    //! [Register MoFEM discrete manager in PETSc]

    DMType dm_name = "DMMOFEM";

    CHKERR DMRegister_MoFEM(dm_name);

    //! [Register MoFEM discrete manager in PETSc


    //! [Create MoAB]

    moab::Core mb_instance;              ///< mesh database

    moab::Interface &moab = mb_instance; ///< mesh database interface

    //! [Create MoAB]


    //! [Create MoFEM]

    MoFEM::Core core(moab);           ///< finite element database

    MoFEM::Interface &m_field = core; ///< finite element database interface

    //! [Create MoFEM]


    //! [Load mesh]

    Simple *simple = m_field.getInterface<Simple>();

    CHKERR simple->getOptions();

    CHKERR simple->loadFile();

    //! [Load mesh]


    //! [Example]

    Example ex(m_field);

    CHKERR ex.runProblem();

    //! [Example]

  }

  CATCH_ERRORS;


  CHKERR MoFEM::Core::Finalize();


#ifdef ADD_CONTACT

#ifdef PYTHON_SFD

  if (Py_FinalizeEx() < 0) {

    exit(120);

  }

#endif

#endif // ADD_CONTACT


  return 0;

}


struct SetUpSchurImpl : public SetUpSchur {


  SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,

                 SmartPetscObj<IS> field_split_is, SmartPetscObj<AO> ao_up)

      : SetUpSchur(), mField(m_field), subDM(sub_dm),

        fieldSplitIS(field_split_is), aoUp(ao_up) {

    if (S) {

      CHK_THROW_MESSAGE(

          MOFEM_DATA_INCONSISTENCY,

          "Is expected that schur matrix is not allocated. This is "

          "possible only is PC is set up twice");

    }

  }

  virtual ~SetUpSchurImpl() {

#ifdef ADD_CONTACT

    A.reset();

    P.reset();

#endif // ADD_CONTACT

    S.reset();

  }


  MoFEMErrorCode setUp(TS solver);

  MoFEMErrorCode preProc();

  MoFEMErrorCode postProc();


private:

#ifdef ADD_CONTACT

  SmartPetscObj<Mat> A;

  SmartPetscObj<Mat> P;

#endif // ADD_CONTACT

  SmartPetscObj<Mat> S;


  MoFEM::Interface &mField;

  SmartPetscObj<DM> subDM;        ///< field split sub dm

  SmartPetscObj<IS> fieldSplitIS; ///< IS for split Schur block

  SmartPetscObj<AO> aoUp;         ///> main DM to subDM

};


MoFEMErrorCode SetUpSchurImpl::setUp(TS solver) {

  MoFEMFunctionBegin;

  auto simple = mField.getInterface<Simple>();

  auto pip_mng = mField.getInterface<PipelineManager>();


  SNES snes;

  CHKERR TSGetSNES(solver, &snes);

  KSP ksp;

  CHKERR SNESGetKSP(snes, &ksp);

  CHKERR KSPSetFromOptions(ksp);


  PC pc;

  CHKERR KSPSetFromOptions(ksp);

  CHKERR KSPGetPC(ksp, &pc);

  PetscBool is_pcfs = PETSC_FALSE;

  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);

  if (is_pcfs) {

    if (S) {

      CHK_THROW_MESSAGE(

          MOFEM_DATA_INCONSISTENCY,

          "Is expected that schur matrix is not allocated. This is "

          "possible only is PC is set up twice");

    }


#ifdef ADD_CONTACT

    auto ts_ctx_ptr = getDMTsCtx(simple->getDM());

    A = createDMMatrix(simple->getDM());

    P = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);

    CHKERR TSSetIJacobian(solver, A, P, TsSetIJacobian, ts_ctx_ptr.get());

#endif // ADD_CONTACT

    S = createDMMatrix(subDM);

    CHKERR MatSetBlockSize(S, SPACE_DIM);


    auto set_ops = [&]() {

      MoFEMFunctionBegin;

      auto pip_mng = mField.getInterface<PipelineManager>();


#ifndef ADD_CONTACT

      // Boundary

      pip_mng->getOpBoundaryLhsPipeline().push_front(

          new OpSchurAssembleBegin());

      pip_mng->getOpBoundaryLhsPipeline().push_back(

          new OpSchurAssembleEnd<SCHUR_DGESV>(


              {"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), aoUp},

              {SmartPetscObj<Mat>(), S}, {false, false}


              ));

      // Domain

      pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());

      pip_mng->getOpDomainLhsPipeline().push_back(

          new OpSchurAssembleEnd<SCHUR_DGESV>(


              {"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), aoUp},

              {SmartPetscObj<Mat>(), S}, {false, false}


              ));

#else


      double eps_stab = 1e-4;

      CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,

                                   PETSC_NULL);


      using B = FormsIntegrators<BoundaryEleOpStab>::Assembly<

          SCHUR>::BiLinearForm<IT>;

      using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;


      // Boundary

      pip_mng->getOpBoundaryLhsPipeline().push_front(

          new OpSchurAssembleBegin());

      pip_mng->getOpBoundaryLhsPipeline().push_back(

          new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {

            return eps_stab;

          }));

      pip_mng->getOpBoundaryLhsPipeline().push_back(

          new OpSchurAssembleEnd<SCHUR_DGESV>(


              {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},

              {SmartPetscObj<AO>(), SmartPetscObj<AO>(), aoUp},

              {SmartPetscObj<Mat>(), SmartPetscObj<Mat>(), S},

              {false, false, false}


              ));

      // Domain

      pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());

      pip_mng->getOpDomainLhsPipeline().push_back(

          new OpSchurAssembleEnd<SCHUR_DGESV>(


              {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},

              {SmartPetscObj<AO>(), SmartPetscObj<AO>(), aoUp},

              {SmartPetscObj<Mat>(), SmartPetscObj<Mat>(), S},

              {false, false, false}


              ));

#endif // ADD_CONTACT

      MoFEMFunctionReturn(0);

    };


    auto set_assemble_elems = [&]() {

      MoFEMFunctionBegin;

      auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();

      schur_asmb_pre_proc->preProcessHook = [this]() {

        MoFEMFunctionBegin;

#ifdef ADD_CONTACT

        CHKERR MatZeroEntries(A);

        CHKERR MatZeroEntries(P);

#endif // ADD_CONTACT

        CHKERR MatZeroEntries(S);

        MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";

        MoFEMFunctionReturn(0);

      };

      auto schur_asmb_post_proc = boost::make_shared<FEMethod>();


      schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {

        MoFEMFunctionBegin;

        MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";


#ifndef ADD_CONTACT

        CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(

            mField, schur_asmb_post_proc, 1)();

#else  // ADD_CONTACT

        CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);

        CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

        // Apply essential constrains to A matrix

        CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(

            mField, schur_asmb_post_proc, 1, A)();

        CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);

        CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);

        CHKERR MatAXPY(P, 1, A, SAME_NONZERO_PATTERN);

#endif // ADD_CONTACT


        // Apply essential constrains to Schur complement

        CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);

        CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);

        CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(

            mField, schur_asmb_post_proc, 1, S, aoUp)();


        MoFEMFunctionReturn(0);

      };

      auto ts_ctx_ptr = getDMTsCtx(simple->getDM());

      ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);

      ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);

      MoFEMFunctionReturn(0);

    };


    auto set_pc = [&]() {

      MoFEMFunctionBegin;

      CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);

      CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);

      MoFEMFunctionReturn(0);

    };


    CHKERR set_ops();

    CHKERR set_pc();

    CHKERR set_assemble_elems();


  } else {

    pip_mng->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());

    pip_mng->getOpBoundaryLhsPipeline().push_back(

        new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));

    pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());

    pip_mng->getOpDomainLhsPipeline().push_back(

        new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));

  }


  // we do not those anymore

  subDM.reset();

  fieldSplitIS.reset();

  // aoUp.reset();

  MoFEMFunctionReturn(0);

}


boost::shared_ptr<SetUpSchur>

SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,

                             SmartPetscObj<DM> sub_dm, SmartPetscObj<IS> is_sub,

                             SmartPetscObj<AO> ao_up) {

  return boost::shared_ptr<SetUpSchur>(

      new SetUpSchurImpl(m_field, sub_dm, is_sub, ao_up));

}

p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90

param_file
std::string param_file
Definition: DefaultParams.hpp:113

HenckyOps.hpp

IntegrationRules.hpp

MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345

MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311

MatrixFunction.hpp

MoFEM.hpp

PlasticNaturalBCs.hpp

PlasticOps.hpp

simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69

help
static char help[]
Definition: activate_deactivate_dofs.cpp:13

main
int main()
Definition: adol-c_atom.cpp:46

eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11

SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16

DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34

BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39

B

BiLinearForm

FTensor::Ddg
Definition: Ddg_value.hpp:8

FTensor::Index
Definition: Index.hpp:24

FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49

FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:14

LinearForm

Range

UBlasMatrix< double >

double

ROW
@ ROW
Definition: definitions.h:123

CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372

MAX_DOFS_ON_ENTITY
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236

AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60

NOBASE
@ NOBASE
Definition: definitions.h:59

DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66

CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595

MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447

FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82

L2
@ L2
field with C-1 continuity
Definition: definitions.h:88

H1
@ H1
continuous field
Definition: definitions.h:85

HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86

HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87

MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215

MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346

CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576

MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33

MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31

ApproximationBaseNames
static const char *const ApproximationBaseNames[]
Definition: definitions.h:72

MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416

CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535

MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440

order
constexpr int order
Definition: dg_projection.cpp:18

SPACE_DIM
constexpr int SPACE_DIM
[Define dimension]
Definition: elastic.cpp:18

dim
const int dim
Definition: elec_phys_2D.cpp:12

OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
Definition: fluid_structure_eigenproblem.cpp:59

poisson_ratio
auto poisson_ratio
Definition: fluid_structure_eigenproblem.cpp:89

m
FTensor::Index< 'm', SPACE_DIM > m
Definition: fluid_structure_eigenproblem.cpp:65

t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137

MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:219

MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483

MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:442

MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242

MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47

MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003

MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:275

MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:988

MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:128

MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition: FormsIntegrators.hpp:104

MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:128

MoFEM::SCHUR
@ SCHUR
Definition: FormsIntegrators.hpp:104

MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308

MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339

MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284

MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.

MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:563

i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27

D
double D
Definition: initial_diffusion.cpp:44

v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40

l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21

j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19

k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20

OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
Definition: multifield_plasticity.cpp:73

ContactOps
Definition: EshelbianContact.hpp:10

ContactOps::cn_contact
double cn_contact
Definition: EshelbianContact.hpp:19

EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66

ElectroPhysiology::b
const double b
Definition: ElecPhysOperators.hpp:30

HenckyOps
Definition: HenckyOps.hpp:11

MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56

MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10

MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1918

MoFEM::DMMoFEMTSSetMonitor
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1042

MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1045

MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142

MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182

MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162

MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:217

MoFEM::getDMSubData
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1061

MoFEM::createAOMappingIS
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
Definition: PetscSmartObj.hpp:302

MoFEM::matDuplicate
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
Definition: PetscSmartObj.hpp:230

MoFEM::getDMSnesCtx
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1031

MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:137

PlasticOps
Definition: PlasticNaturalBCs.hpp:13

convert.int
int
Definition: convert.py:64

convert.type
type
Definition: convert.py:64

sdf.r
int r
Definition: sdf.py:8

OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29

young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:172

AT
constexpr AssemblyType AT
Definition: plastic.cpp:44

C1_k
double C1_k
Kinematic hardening.
Definition: plastic.cpp:180

Qinf
double Qinf
Saturation yield stress.
Definition: plastic.cpp:178

IT
constexpr IntegrationType IT
Definition: plastic.cpp:47

rho
double rho
Definition: plastic.cpp:191

EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13

is_quasi_static
PetscBool is_quasi_static
Definition: plastic.cpp:190

alpha_damping
double alpha_damping
Definition: plastic.cpp:192

visH
double visH
Viscous hardening.
Definition: plastic.cpp:176

kinematic_hardening
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
Definition: plastic.cpp:142

set_timer
PetscBool set_timer
Set timer.
Definition: plastic.cpp:168

iso_hardening_dtau
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
Definition: plastic.cpp:128

scale
double scale
Definition: plastic.cpp:170

size_symm
constexpr auto size_symm
Definition: plastic.cpp:42

zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:177

H
double H
Hardening.
Definition: plastic.cpp:175

tau_order
int tau_order
Order of tau files.
Definition: plastic.cpp:186

iso_hardening_exp
double iso_hardening_exp(double tau, double b_iso)
Definition: plastic.cpp:114

cn0
double cn0
Definition: plastic.cpp:182

b_iso
double b_iso
Saturation exponent.
Definition: plastic.cpp:179

is_large_strains
PetscBool is_large_strains
Large strains.
Definition: plastic.cpp:167

geom_order
int geom_order
Order if fixed.
Definition: plastic.cpp:188

sigmaY
double sigmaY
Yield stress.
Definition: plastic.cpp:174

iso_hardening
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
Definition: plastic.cpp:123

kinematic_hardening_dplastic_strain
auto kinematic_hardening_dplastic_strain(double C1_k)
Definition: plastic.cpp:156

SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61

ep_order
int ep_order
Order of ep files.
Definition: plastic.cpp:187

cn1
double cn1
Definition: plastic.cpp:183

CONTACT_SPACE
constexpr FieldSpace CONTACT_SPACE
Definition: plastic.cpp:52

m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80

BoundaryEleOpStab
Element used to specialise assembly.
Definition: contact.cpp:86

ElementsAndOps
Definition: child_and_parent.cpp:18

Example::ScaledTimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: plastic.cpp:249

Example
[Example]
Definition: plastic.cpp:226

Example::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:245

Example::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:243

Example::tsSolve
MoFEMErrorCode tsSolve()
Definition: plastic.cpp:778

Example::base
FieldApproximationBase base
Definition: plot_base.cpp:68

Example::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:440

Example::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:244

Example::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:593

Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:262

Example::mField
MoFEM::Interface & mField
Definition: plastic.cpp:233

Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:274

Example::bC
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:548

Example::reactionFe
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:241

MoFEM::AddFluxToLhsPipelineImpl
Definition: Natural.hpp:49

MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46

MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503

MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25

MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21

MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0

MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0

MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82

MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72

MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112

MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16

MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:72

MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30

MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127

MoFEM::EssentialPostProcLhs
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition: Essential.hpp:47

MoFEM::EssentialPostProcRhs
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition: Essential.hpp:55

MoFEM::EssentialPreProc
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:39

MoFEM::EssentialPreProcReaction
Class (Function) to calculate residual side diagonal.
Definition: Essential.hpp:63

MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:96

MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:25

MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:546

MoFEM::ForcesAndSourcesCore::UserDataOperator::getKSPA
Mat getKSPA() const
Definition: ForcesAndSourcesCore.hpp:1092

MoFEM::ForcesAndSourcesCore::UserDataOperator::getKSPB
Mat getKSPB() const
Definition: ForcesAndSourcesCore.hpp:1100

MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:296

MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23

MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104

MoFEM::NaturalBC::Assembly::BiLinearForm
Definition: Natural.hpp:74

MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67

MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65

MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:170

MoFEM::OpBaseImpl::MatSetValuesHook
boost::function< MoFEMErrorCode(ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
Definition: FormsIntegrators.hpp:210

MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:596

MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:468

MoFEM::OpFluxLhsImpl
Definition: Natural.hpp:43

MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39

MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1286

MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698

MoFEM::OpSchurAssembleBegin
Clear Schur complement internal data.
Definition: Schur.hpp:22

MoFEM::OpSchurAssembleEnd
Assemble Schur complement.
Definition: Schur.hpp:104

MoFEM::OpSetHOWeightsOnSubDim
Definition: HODataOperators.hpp:145

MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38

MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24

MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401

MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667

MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21

MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24

MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27

MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition: PetscSmartObj.hpp:79

MoFEM::TimeScale
Force scale operator for reading two columns.
Definition: ScalingMethod.hpp:32

MoFEM::TimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: ScalingMethod.cpp:122

MoFEM::TimeScale::TimeScale
TimeScale(std::string file_name="", bool error_if_file_not_given=false)
TimeScale constructor.
Definition: ScalingMethod.cpp:22

MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93

MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:109

MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26

Monitor
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:783

PlasticOps::CommonData::activityData
static std::array< int, 5 > activityData
Definition: PlasticOps.hpp:107

SetUpSchur
[Push operators to pipeline]
Definition: plastic.cpp:755

SetUpSchur::SetUpSchur
SetUpSchur()=default

SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
Create data structure for handling Schur complement.
Definition: plastic.cpp:1474

SetUpSchur::setUp
virtual MoFEMErrorCode setUp(TS solver)=0

SetUpSchurImpl
Definition: SetUpSchurImpl.cpp:9

SetUpSchurImpl::A
SmartPetscObj< Mat > A
Definition: contact.cpp:854

SetUpSchurImpl::subDM
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1296

SetUpSchurImpl::setUp
MoFEMErrorCode setUp(KSP solver)
Definition: SetUpSchurImpl.cpp:39

SetUpSchurImpl::S
SmartPetscObj< Mat > S
Definition: SetUpSchurImpl.cpp:30

SetUpSchurImpl::P
SmartPetscObj< Mat > P
Definition: SetUpSchurImpl.cpp:29

SetUpSchurImpl::fieldSplitIS
SmartPetscObj< IS > fieldSplitIS
IS for split Schur block.
Definition: plastic.cpp:1297

SetUpSchurImpl::~SetUpSchurImpl
virtual ~SetUpSchurImpl()
Definition: SetUpSchurImpl.cpp:21

SetUpSchurImpl::postProc
MoFEMErrorCode postProc()
Definition: SetUpSchurImpl.cpp:204

SetUpSchurImpl::aoUp
SmartPetscObj< AO > aoUp
Definition: SetUpSchurImpl.cpp:31

SetUpSchurImpl::preProc
MoFEMErrorCode preProc()
Definition: SetUpSchurImpl.cpp:193

SetUpSchurImpl::mField
MoFEM::Interface & mField
Definition: SetUpSchurImpl.cpp:35