Webpage/dynamic_first_order_con_law_8cpp-example.html

/**

 * \file dynamic_first_order_con_law.cpp

 * \example dynamic_first_order_con_law.cpp

 *

 * Explicit first order conservation laws for solid dynamics

 *

 */


#include <MoFEM.hpp>

#include <MatrixFunction.hpp>


using namespace MoFEM;


template <typename T> inline double trace(FTensor::Tensor2<T, 2, 2> &t_stress) {

  constexpr double third = boost::math::constants::third<double>();

  return (t_stress(0, 0) + t_stress(1, 1));

};


template <typename T> inline double trace(FTensor::Tensor2<T, 3, 3> &t_stress) {

  constexpr double third = boost::math::constants::third<double>();

  return (t_stress(0, 0) + t_stress(1, 1) + t_stress(2, 2));

};


constexpr int SPACE_DIM =

    EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh


using EntData = EntitiesFieldData::EntData;

using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;

using DomainEleOp = DomainEle::UserDataOperator;

using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;

using PostProcFaceEle =

    PostProcBrokenMeshInMoab<FaceElementForcesAndSourcesCore>;


using BoundaryEle =

    PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;

using BoundaryEleOp = BoundaryEle::UserDataOperator;

using SetPtsData = FieldEvaluatorInterface::SetPtsData;


template <int DIM> struct PostProcEleByDim;


template <> struct PostProcEleByDim<2> {

  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;

  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;

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

};


template <> struct PostProcEleByDim<3> {

  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;

  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;

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

};


using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain;

using SideEle = PostProcEleByDim<SPACE_DIM>::SideEle;

using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy;


using OpMassV = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<

    GAUSS>::OpMass<1, SPACE_DIM>;

using OpMassF = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<

    GAUSS>::OpMass<1, SPACE_DIM * SPACE_DIM>;


using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<

    PETSC>::LinearForm<GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 1>;


using OpBodyForce = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<

    GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 0>;


using DomainNaturalBC =

    NaturalBC<DomainEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;

using OpBodyForceVector =

    DomainNaturalBC::OpFlux<NaturalMeshsetTypeVectorScaling<BLOCKSET>, 1,

                            SPACE_DIM>;


using OpGradTimesTensor2 =

    FormsIntegrators<DomainEleOp>::Assembly<AssemblyType::PETSC>::LinearForm<

        IntegrationType::GAUSS>::OpGradTimesTensor<1, SPACE_DIM, SPACE_DIM>;


using OpRhsTestPiola =

    FormsIntegrators<DomainEleOp>::Assembly<AssemblyType::PETSC>::LinearForm<

        IntegrationType::GAUSS>::OpBaseTimesVector<1, SPACE_DIM * SPACE_DIM, 1>;


using OpGradTimesPiola =

    FormsIntegrators<DomainEleOp>::Assembly<AssemblyType::PETSC>::LinearForm<

        IntegrationType::GAUSS>::OpGradTimesTensor<1, SPACE_DIM, SPACE_DIM>;


using BoundaryNaturalBC =

    NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;

using OpForce = BoundaryNaturalBC::OpFlux<NaturalForceMeshsets, 1, SPACE_DIM>;


/** \brief Save field DOFS on vertices/tags

 */


constexpr double omega = 1.;

constexpr double young_modulus = 1.;

constexpr double poisson_ratio = 0.;

double bulk_modulus_K = young_modulus / (3. * (1. - 2. * poisson_ratio));

double shear_modulus_G = young_modulus / (2. * (1. + poisson_ratio));

double mu = young_modulus / (2. * (1. + poisson_ratio));

double lamme_lambda = young_modulus * poisson_ratio /

                      ((1. + poisson_ratio) * (1. - 2. * poisson_ratio));


// Operator to Calculate F

template <int DIM_0, int DIM_1>

struct OpCalculateFStab : public ForcesAndSourcesCore::UserDataOperator {

  OpCalculateFStab(boost::shared_ptr<MatrixDouble> def_grad_ptr,

                   boost::shared_ptr<MatrixDouble> def_grad_stab_ptr,

                   boost::shared_ptr<MatrixDouble> def_grad_dot_ptr,

                   double tau_F_ptr, double xi_F_ptr,

                   boost::shared_ptr<MatrixDouble> grad_x_ptr,

                   boost::shared_ptr<MatrixDouble> grad_vel_ptr)

      : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST),

        defGradPtr(def_grad_ptr), defGradStabPtr(def_grad_stab_ptr),

        defGradDotPtr(def_grad_dot_ptr), tauFPtr(tau_F_ptr), xiF(xi_F_ptr),

        gradxPtr(grad_x_ptr), gradVelPtr(grad_vel_ptr) {}


  MoFEMErrorCode doWork(int side, EntityType type,

                        DataForcesAndSourcesCore::EntData &data) {

    MoFEMFunctionBegin;

    // Define Indicies

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

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


    // Number of Gauss points

    const size_t nb_gauss_pts = getGaussPts().size2();


    defGradStabPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false);

    defGradStabPtr->clear();


    // Extract matrix from data matrix

    auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);

    auto t_Fstab = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradStabPtr);

    auto t_F_dot = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradDotPtr);


    // tau_F = alpha deltaT

    auto tau_F = tauFPtr;

    double xi_F = xiF;

    auto t_gradx = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradxPtr);

    auto t_gradVel = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradVelPtr);


    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {

      // Stabilised Deformation Gradient

      t_Fstab(i, j) = t_F(i, j) + tau_F * (t_gradVel(i, j) - t_F_dot(i, j)) +

                      xi_F * (t_gradx(i, j) - t_F(i, j));


      ++t_F;

      ++t_Fstab;

      ++t_gradVel;

      ++t_F_dot;


      ++t_gradx;

    }


    MoFEMFunctionReturn(0);

  }


private:

  double tauFPtr;

  double xiF;

  boost::shared_ptr<MatrixDouble> defGradPtr;

  boost::shared_ptr<MatrixDouble> defGradStabPtr;

  boost::shared_ptr<MatrixDouble> defGradDotPtr;

  boost::shared_ptr<MatrixDouble> gradxPtr;

  boost::shared_ptr<MatrixDouble> gradVelPtr;

};


// Operator to Calculate P

template <int DIM_0, int DIM_1>

struct OpCalculatePiola : public ForcesAndSourcesCore::UserDataOperator {

  OpCalculatePiola(double shear_modulus, double bulk_modulus, double m_u,

                   double lambda_lamme,

                   boost::shared_ptr<MatrixDouble> first_piola_ptr,

                   boost::shared_ptr<MatrixDouble> def_grad_ptr)

      : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST),

        shearModulus(shear_modulus), bulkModulus(bulk_modulus), mU(m_u),

        lammeLambda(lambda_lamme), firstPiolaPtr(first_piola_ptr),

        defGradPtr(def_grad_ptr) {}


  MoFEMErrorCode doWork(int side, EntityType type,

                        DataForcesAndSourcesCore::EntData &data) {

    MoFEMFunctionBegin;

    // Define Indicies

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

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

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


    // Define Kronecker Delta

    constexpr auto t_kd = FTensor::Kronecker_Delta<double>();


    // Number of Gauss points

    const size_t nb_gauss_pts = getGaussPts().size2();


    // Resize Piola

    firstPiolaPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false); // ignatios check

    firstPiolaPtr->clear();


    // Extract matrix from data matrix

    auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*firstPiolaPtr);

    auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);

    const double two_o_three = 2. / 3.;

    const double trace_t_dk = DIM_0;

    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {


      t_P(i, j) = shearModulus * (t_F(i, j) + t_F(j, i) - 2. * t_kd(i, j) -

                                  two_o_three * trace(t_F) * t_kd(i, j) +

                                  two_o_three * trace_t_dk * t_kd(i, j)) +

                  bulkModulus * trace(t_F) * t_kd(i, j) -

                  bulkModulus * trace_t_dk * t_kd(i, j);


      ++t_F;

      ++t_P;

    }


    MoFEMFunctionReturn(0);

  }


private:

  double shearModulus;

  double bulkModulus;

  double mU;

  double lammeLambda;

  boost::shared_ptr<MatrixDouble> firstPiolaPtr;

  boost::shared_ptr<MatrixDouble> defGradPtr;

};


template <int DIM>

struct OpCalculateDisplacement : public ForcesAndSourcesCore::UserDataOperator {

  OpCalculateDisplacement(boost::shared_ptr<MatrixDouble> spatial_pos_ptr,

                          boost::shared_ptr<MatrixDouble> reference_pos_ptr,

                          boost::shared_ptr<MatrixDouble> u_ptr)

      : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST),

        xPtr(spatial_pos_ptr), XPtr(reference_pos_ptr), uPtr(u_ptr) {}


  MoFEMErrorCode doWork(int side, EntityType type,

                        DataForcesAndSourcesCore::EntData &data) {

    MoFEMFunctionBegin;

    // Define Indicies

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


    // Number of Gauss points

    const size_t nb_gauss_pts = getGaussPts().size2();


    uPtr->resize(DIM, nb_gauss_pts, false); // ignatios check

    uPtr->clear();


    // Extract matrix from data matrix

    auto t_x = getFTensor1FromMat<DIM>(*xPtr);

    auto t_X = getFTensor1FromMat<DIM>(*XPtr);

    auto t_u = getFTensor1FromMat<DIM>(*uPtr);

    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {


      t_u(i) = t_x(i) - t_X(i);

      ++t_u;

      ++t_x;

      ++t_X;

    }


    MoFEMFunctionReturn(0);

  }


private:

  boost::shared_ptr<MatrixDouble> xPtr;

  boost::shared_ptr<MatrixDouble> XPtr;

  boost::shared_ptr<MatrixDouble> uPtr;

};


template <int DIM_0, int DIM_1>

struct OpCalculatePiolaIncompressibleNH

    : public ForcesAndSourcesCore::UserDataOperator {

  OpCalculatePiolaIncompressibleNH(

      double shear_modulus, double bulk_modulus, double m_u,

      double lambda_lamme, boost::shared_ptr<MatrixDouble> first_piola_ptr,

      boost::shared_ptr<MatrixDouble> def_grad_ptr,

      boost::shared_ptr<MatrixDouble> inv_def_grad_ptr,

      boost::shared_ptr<VectorDouble> det)

      : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST),

        shearModulus(shear_modulus), bulkModulus(bulk_modulus), mU(m_u),

        lammeLambda(lambda_lamme), firstPiolaPtr(first_piola_ptr),

        defGradPtr(def_grad_ptr), invDefGradPtr(inv_def_grad_ptr), dEt(det) {}


  MoFEMErrorCode doWork(int side, EntityType type,

                        DataForcesAndSourcesCore::EntData &data) {

    MoFEMFunctionBegin;

    // Define Indicies

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

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

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

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


    // Define Kronecker Delta

    constexpr auto t_kd = FTensor::Kronecker_Delta<double>();


    // Number of Gauss points

    const size_t nb_gauss_pts = getGaussPts().size2();


    // Resize Piola

    firstPiolaPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false); // ignatios check

    firstPiolaPtr->clear();


    // Extract matrix from data matrix

    auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*firstPiolaPtr);

    auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);

    auto t_inv_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*invDefGradPtr);

    auto t_det = getFTensor0FromVec<1>(*dEt);

    const double two_o_three = 2. / 3.;

    const double one_o_three = 1. / 3.;

    const double bulk_mod = bulkModulus;

    const double shear_mod = shearModulus;

    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {


      // Nearly incompressible NH

      // volumetric part

      t_P(i, j) = bulk_mod * (t_det - 1.) * t_det * t_inv_F(j, i);

      // deviatoric part

      t_P(i, j) +=

          shear_mod * pow(t_det, two_o_three) *

          (t_F(i, j) - one_o_three * (t_F(l, k) * t_F(l, k)) * t_inv_F(j, i));


      ++t_F;

      ++t_P;

      ++t_inv_F;

      ++t_det;

    }


    MoFEMFunctionReturn(0);

  }


private:

  double shearModulus;

  double bulkModulus;

  double mU;

  double lammeLambda;

  boost::shared_ptr<MatrixDouble> firstPiolaPtr;

  boost::shared_ptr<MatrixDouble> defGradPtr;

  boost::shared_ptr<MatrixDouble> invDefGradPtr;

  boost::shared_ptr<VectorDouble> dEt;

};


template <int DIM_0, int DIM_1>

struct OpCalculateDeformationGradient

    : public ForcesAndSourcesCore::UserDataOperator {

  OpCalculateDeformationGradient(

      boost::shared_ptr<MatrixDouble> def_grad_ptr,

      boost::shared_ptr<MatrixDouble> grad_tensor_ptr)

      : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST),

        defGradPtr(def_grad_ptr), gradTensorPtr(grad_tensor_ptr) {}


  MoFEMErrorCode doWork(int side, EntityType type,

                        DataForcesAndSourcesCore::EntData &data) {

    MoFEMFunctionBegin;

    // Define Indicies

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

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


    // Define Kronecker Delta

    constexpr auto t_kd = FTensor::Kronecker_Delta<double>();


    // Number of Gauss points

    const size_t nb_gauss_pts = getGaussPts().size2();


    // Resize Piola

    defGradPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false);

    defGradPtr->clear();


    // Extract matrix from data matrix

    auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);

    auto t_H = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradTensorPtr);

    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {


      t_F(i, j) = t_H(i, j) + t_kd(i, j);


      ++t_F;

      ++t_H;

    }


    MoFEMFunctionReturn(0);

  }


private:

  boost::shared_ptr<MatrixDouble> gradTensorPtr;

  boost::shared_ptr<MatrixDouble> defGradPtr;

};


struct Example;

struct TSPrePostProc {


  TSPrePostProc() = default;

  virtual ~TSPrePostProc() = default;


  /**

   * @brief Used to setup TS solver

   *

   * @param ts

   * @return MoFEMErrorCode

   */

  MoFEMErrorCode tsSetUp(TS ts);


  // SmartPetscObj<VecScatter> getScatter(Vec x, Vec y, enum FR fr);

  Example *fsRawPtr;

  static MoFEMErrorCode tsPostStage(TS ts, PetscReal stagetime,

                                    PetscInt stageindex, Vec *Y);

  static MoFEMErrorCode tsPostStep(TS ts);

  static MoFEMErrorCode tsPreStep(TS ts);

};


static boost::weak_ptr<TSPrePostProc> tsPrePostProc;


struct LinMomTimeScale : public MoFEM::TimeScale {

  using MoFEM::TimeScale::TimeScale;

  double getScale(const double time) {

    return sin(2. * M_PI * MoFEM::TimeScale::getScale(time));

  };

};


struct CommonData {

  SmartPetscObj<Mat> M;   ///< Mass matrix

  SmartPetscObj<KSP> ksp; ///< Linear solver

};


struct Example {


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


  MoFEMErrorCode runProblem();


private:

  MoFEM::Interface &mField;


  MoFEMErrorCode readMesh();

  MoFEMErrorCode setupProblem();

  MoFEMErrorCode boundaryCondition();

  MoFEMErrorCode assembleSystem();

  MoFEMErrorCode solveSystem();

  MoFEMErrorCode outputResults();

  MoFEMErrorCode checkResults();

  friend struct TSPrePostProc;


  struct DynamicFirstOrderConsSinusTimeScale : public MoFEM::TimeScale {

    using MoFEM::TimeScale::TimeScale;

    double getScale(const double time) { return 0.001 * sin(0.1 * time); };

  };


  struct DynamicFirstOrderConsConstantTimeScale : public MoFEM::TimeScale {

    using MoFEM::TimeScale::TimeScale;

    double getScale(const double time) { return 0.001; };

  };

};


//! [Run problem]

MoFEMErrorCode Example::runProblem() {

  MoFEMFunctionBegin;

  CHKERR readMesh();

  CHKERR setupProblem();

  CHKERR boundaryCondition();

  CHKERR assembleSystem();

  CHKERR solveSystem();

  CHKERR outputResults();

  CHKERR checkResults();

  MoFEMFunctionReturn(0);

}

//! [Run problem]


//! [Read mesh]

MoFEMErrorCode Example::readMesh() {

  MoFEMFunctionBegin;

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


  CHKERR simple->getOptions();

  CHKERR simple->loadFile();

  MoFEMFunctionReturn(0);

}

//! [Read mesh]


//! [Set up problem]

MoFEMErrorCode Example::setupProblem() {

  MoFEMFunctionBegin;

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

  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };

  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};

  PetscInt choice_base_value = AINSWORTH;

  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,

                              LASBASETOPT, &choice_base_value, PETSC_NULL);


  FieldApproximationBase base;

  switch (choice_base_value) {

  case AINSWORTH:

    base = AINSWORTH_LEGENDRE_BASE;

    MOFEM_LOG("WORLD", Sev::inform)

        << "Set AINSWORTH_LEGENDRE_BASE for displacements";

    break;

  case DEMKOWICZ:

    base = DEMKOWICZ_JACOBI_BASE;

    MOFEM_LOG("WORLD", Sev::inform)

        << "Set DEMKOWICZ_JACOBI_BASE for displacements";

    break;

  default:

    base = LASTBASE;

    break;

  }

  // Add field

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

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

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

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

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

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

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


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

  int order = 2;

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

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

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

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

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

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

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

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

  CHKERR simple->setUp();


  auto project_ho_geometry = [&]() {

    Projection10NodeCoordsOnField ent_method_x(mField, "x_1");

    CHKERR mField.loop_dofs("x_1", ent_method_x);

    Projection10NodeCoordsOnField ent_method_x_2(mField, "x_2");

    CHKERR mField.loop_dofs("x_2", ent_method_x_2);


    Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");

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

  };

  CHKERR project_ho_geometry();


  MoFEMFunctionReturn(0);

}

//! [Set up problem]


//! [Boundary condition]

MoFEMErrorCode Example::boundaryCondition() {

  MoFEMFunctionBegin;


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

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

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

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


  PetscBool sin_time_function = PETSC_FALSE;

  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-sin_time_function",

                             &sin_time_function, PETSC_NULL);


  if (sin_time_function)

    time_scale = boost::make_shared<DynamicFirstOrderConsSinusTimeScale>();

  else

    time_scale = boost::make_shared<DynamicFirstOrderConsConstantTimeScale>();


  pipeline_mng->getBoundaryExplicitRhsFE().reset();

  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(

      pipeline_mng->getOpBoundaryExplicitRhsPipeline(), {NOSPACE}, "GEOMETRY");


  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(

      pipeline_mng->getOpBoundaryExplicitRhsPipeline(), mField, "V",

      {time_scale}, "FORCE", Sev::inform);


  auto integration_rule = [](int, int, int approx_order) {

    return 2 * approx_order;

  };


  CHKERR pipeline_mng->setBoundaryExplicitRhsIntegrationRule(integration_rule);

  CHKERR pipeline_mng->setDomainExplicitRhsIntegrationRule(integration_rule);


  CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(

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


  auto get_pre_proc_hook = [&]() {

    return EssentialPreProc<DisplacementCubitBcData>(

        mField, pipeline_mng->getDomainExplicitRhsFE(), {time_scale});

  };

  pipeline_mng->getDomainExplicitRhsFE()->preProcessHook = get_pre_proc_hook();


  MoFEMFunctionReturn(0);

}

//! [Boundary condition]


MoFEMErrorCode TSPrePostProc::tsPostStage(TS ts, PetscReal stagetime,

                                          PetscInt stageindex, Vec *Y) {

  MoFEMFunctionBegin;

  // cerr << "tsPostStage " <<"\n";

  if (auto ptr = tsPrePostProc.lock()) {

    auto &m_field = ptr->fsRawPtr->mField;

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

    auto *pipeline_mng = m_field.getInterface<PipelineManager>();


    auto fb = m_field.getInterface<FieldBlas>();

    double dt;

    CHKERR TSGetTimeStep(ts, &dt);

    double time;

    CHKERR TSGetTime(ts, &time);

    PetscInt num_stages;

    Vec *stage_solutions;


    CHKERR TSGetStages(ts, &num_stages, &stage_solutions);

    PetscPrintf(PETSC_COMM_WORLD, "Check timestep %d time %e dt %e\n",

                num_stages, time, dt);


    const double inv_num_step = (double)num_stages;

    CHKERR fb->fieldCopy(1., "x_1", "x_2");

    CHKERR fb->fieldAxpy(dt, "V", "x_2");

    CHKERR fb->fieldCopy(1., "x_2", "x_1");


    CHKERR fb->fieldCopy(-inv_num_step / dt, "F_0", "F_dot");

    CHKERR fb->fieldAxpy(inv_num_step / dt, "F", "F_dot");

    CHKERR fb->fieldCopy(1., "F", "F_0");

  }

  MoFEMFunctionReturn(0);

}


MoFEMErrorCode TSPrePostProc::tsPostStep(TS ts) {

  MoFEMFunctionBegin;


  if (auto ptr = tsPrePostProc.lock()) {

    auto &m_field = ptr->fsRawPtr->mField;

    auto fb = m_field.getInterface<FieldBlas>();

    double dt;

    CHKERR TSGetTimeStep(ts, &dt);

    double time;

    CHKERR TSGetTime(ts, &time);

  }

  MoFEMFunctionReturn(0);

}


MoFEMErrorCode TSPrePostProc::tsPreStep(TS ts) {

  MoFEMFunctionBegin;


  if (auto ptr = tsPrePostProc.lock()) {

    auto &m_field = ptr->fsRawPtr->mField;

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

    auto *pipeline_mng = m_field.getInterface<PipelineManager>();


    double dt;

    CHKERR TSGetTimeStep(ts, &dt);

    double time;

    CHKERR TSGetTime(ts, &time);

    int step_num;

    CHKERR TSGetStepNumber(ts, &step_num);

  }

  MoFEMFunctionReturn(0);

}


//! [Push operators to pipeline]

MoFEMErrorCode Example::assembleSystem() {

  MoFEMFunctionBegin;

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

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


  auto get_body_force = [this](const double, const double, const double) {

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

    FTensor::Tensor1<double, SPACE_DIM> t_source;

    t_source(i) = 0.;

    t_source(0) = 0.1;

    t_source(1) = 1.;

    return t_source;

  };


  // specific time scaling

  auto get_time_scale = [this](const double time) {

    return sin(time * omega * M_PI);

  };


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

    MoFEMFunctionBegin;


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

                                                          "GEOMETRY");


    // Calculate Gradient of velocity

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

    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(

        "V", mat_v_grad_ptr));


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

    gravity_vector_ptr->resize(SPACE_DIM, 1);

    auto set_body_force = [&]() {

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

      MoFEMFunctionBegin;

      auto t_force = getFTensor1FromMat<SPACE_DIM, 0>(*gravity_vector_ptr);

      double unit_weight = 0.;

      CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-unit_weight", &unit_weight,

                                 PETSC_NULL);

      t_force(i) = 0;

      if (SPACE_DIM == 2) {

        t_force(1) = -unit_weight;

      } else if (SPACE_DIM == 3) {

        t_force(2) = unit_weight;

      }

      MoFEMFunctionReturn(0);

    };


    CHKERR set_body_force();

    pip.push_back(new OpBodyForce("V", gravity_vector_ptr,

                                  [](double, double, double) { return 1.; }));


    // Calculate unknown F

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

    pip.push_back(new OpCalculateTensor2FieldValues<SPACE_DIM, SPACE_DIM>(

        "F", mat_H_tensor_ptr));


    //  // Calculate F

    double tau = 0.2;

    CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-tau", &tau, PETSC_NULL);


    double xi = 0.;

    CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-xi", &xi, PETSC_NULL);


    // Calculate P stab

    auto one = [&](const double, const double, const double) {

      return 3. * bulk_modulus_K;

    };

    auto minus_one = [](const double, const double, const double) {

      return -1.;

    };


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

    pip.push_back(new OpCalculateTensor2FieldValues<SPACE_DIM, SPACE_DIM>(

        "F_dot", mat_dot_F_tensor_ptr));


    // Calculate Gradient of Spatial Positions

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

    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(

        "x_2", mat_x_grad_ptr));


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

    pip.push_back(new OpCalculateDeformationGradient<SPACE_DIM, SPACE_DIM>(

        mat_F_tensor_ptr, mat_H_tensor_ptr));


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

    pip.push_back(new OpCalculateFStab<SPACE_DIM, SPACE_DIM>(

        mat_F_tensor_ptr, mat_F_stab_ptr, mat_dot_F_tensor_ptr, tau, xi,

        mat_x_grad_ptr, mat_v_grad_ptr));


    PetscBool is_linear_elasticity = PETSC_TRUE;

    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_linear_elasticity",

                               &is_linear_elasticity, PETSC_NULL);


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

    if (is_linear_elasticity) {

      pip.push_back(new OpCalculatePiola<SPACE_DIM, SPACE_DIM>(

          shear_modulus_G, bulk_modulus_K, mu, lamme_lambda, mat_P_stab_ptr,

          mat_F_stab_ptr));

    } else {

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

      auto det_ptr = boost::make_shared<VectorDouble>();


      pip.push_back(

          new OpInvertMatrix<SPACE_DIM>(mat_F_stab_ptr, det_ptr, inv_F));


      // OpCalculatePiolaIncompressibleNH

      pip.push_back(new OpCalculatePiolaIncompressibleNH<SPACE_DIM, SPACE_DIM>(

          shear_modulus_G, bulk_modulus_K, mu, lamme_lambda, mat_P_stab_ptr,

          mat_F_stab_ptr, inv_F, det_ptr));

    }


    pip.push_back(new OpGradTimesTensor2("V", mat_P_stab_ptr, minus_one));

    pip.push_back(new OpRhsTestPiola("F", mat_v_grad_ptr, one));


    MoFEMFunctionReturn(0);

  };


  CHKERR apply_rhs(pipeline_mng->getOpDomainExplicitRhsPipeline());


  auto integration_rule = [](int, int, int approx_order) {

    return 2 * approx_order;

  };

  CHKERR pipeline_mng->setDomainExplicitRhsIntegrationRule(integration_rule);


  MoFEMFunctionReturn(0);

}

//! [Push operators to pipeline]


/**

 * @brief Monitor solution

 *

 * This functions is called by TS solver at the end of each step. It is used

 * to output results to the hard drive.

 */


struct Monitor : public FEMethod {

  MoFEM::Interface &mField;

  Monitor(SmartPetscObj<DM> dm, MoFEM::Interface &m_field,

          boost::shared_ptr<PostProcEle> post_proc,

          boost::shared_ptr<PostProcFaceEle> post_proc_bdry,

          boost::shared_ptr<MatrixDouble> velocity_field_ptr,

          boost::shared_ptr<MatrixDouble> x2_field_ptr,

          boost::shared_ptr<MatrixDouble> geometry_field_ptr,

          std::array<double, 3> pass_field_eval_coords,

          boost::shared_ptr<SetPtsData> pass_field_eval_data)

      : dM(dm), mField(m_field), postProc(post_proc),

        postProcBdy(post_proc_bdry), velocityFieldPtr(velocity_field_ptr),

        x2FieldPtr(x2_field_ptr), geometryFieldPtr(geometry_field_ptr),

        fieldEvalCoords(pass_field_eval_coords),

        fieldEvalData(pass_field_eval_data){};

  MoFEMErrorCode postProcess() {

    MoFEMFunctionBegin;


    // cerr << "wagawaga\n";

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


    if (SPACE_DIM == 3) {

      CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint3D(

          fieldEvalCoords.data(), 1e-12, simple->getProblemName(),

          simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),

          mField.get_comm_rank(), getCacheWeakPtr().lock(), MF_EXIST, QUIET);

    } else {

      CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint2D(

          fieldEvalCoords.data(), 1e-12, simple->getProblemName(),

          simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),

          mField.get_comm_rank(), getCacheWeakPtr().lock(), MF_EXIST, QUIET);

    }


    if (velocityFieldPtr->size1()) {

      auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velocityFieldPtr);

      auto t_x2_field = getFTensor1FromMat<SPACE_DIM>(*x2FieldPtr);

      auto t_geom = getFTensor1FromMat<SPACE_DIM>(*geometryFieldPtr);


      double u_x = t_x2_field(0) - t_geom(0);

      double u_y = t_x2_field(1) - t_geom(1);

      double u_z = t_x2_field(2) - t_geom(2);


      MOFEM_LOG("SYNC", Sev::inform)

          << "Velocities x: " << t_vel(0) << " y: " << t_vel(1)

          << " z: " << t_vel(2) << "\n";

      MOFEM_LOG("SYNC", Sev::inform) << "Displacement x: " << u_x

                                     << " y: " << u_y << " z: " << u_z << "\n";

    }


    for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(

             std::regex((boost::format("%s(.*)") % "Data_Vertex").str()))) {

      Range ents;

      mField.get_moab().get_entities_by_dimension(m->getMeshset(), 0, ents,

                                                  true);

      auto print_vets = [](boost::shared_ptr<FieldEntity> ent_ptr) {

        MoFEMFunctionBegin;

        if (!(ent_ptr->getPStatus() & PSTATUS_NOT_OWNED)) {

          MOFEM_LOG("SYNC", Sev::inform)

              << "Velocities: " << ent_ptr->getEntFieldData()[0] << " "

              << ent_ptr->getEntFieldData()[1] << " "

              << ent_ptr->getEntFieldData()[2] << "\n";

        }

        MoFEMFunctionReturn(0);

      };

      CHKERR mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(

          print_vets, "V", &ents);

    }

    MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::inform);


    PetscBool print_volume = PETSC_FALSE;

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

                               PETSC_NULL);


    PetscBool print_skin = PETSC_FALSE;

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

                               PETSC_NULL);


    int save_every_nth_step = 1;

    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_step",

                              &save_every_nth_step, PETSC_NULL);

    if (ts_step % save_every_nth_step == 0) {


      if (print_volume) {

        CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);

        CHKERR postProc->writeFile(

            "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");

      }


      if (print_skin) {

        CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcBdy);

        CHKERR postProcBdy->writeFile(

            "out_boundary_" + boost::lexical_cast<std::string>(ts_step) +

            ".h5m");

      }

    }

    MoFEMFunctionReturn(0);

  }


private:

  SmartPetscObj<DM> dM;

  boost::shared_ptr<PostProcEle> postProc;

  boost::shared_ptr<PostProcFaceEle> postProcBdy;

  boost::shared_ptr<MatrixDouble> velocityFieldPtr;

  boost::shared_ptr<MatrixDouble> x2FieldPtr;

  boost::shared_ptr<MatrixDouble> geometryFieldPtr;

  std::array<double, 3> fieldEvalCoords;

  boost::shared_ptr<SetPtsData> fieldEvalData;

};


//! [Solve]

MoFEMErrorCode Example::solveSystem() {

  MoFEMFunctionBegin;

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

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


  auto dm = simple->getDM();


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

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


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

    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("V", v_ptr));

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

    pip.push_back(

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


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

    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("x_1", x_ptr));


    // Calculate unknown F

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

    pip.push_back(new OpCalculateTensor2FieldValues<SPACE_DIM, SPACE_DIM>(

        "F", mat_H_tensor_ptr));


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

    pip.push_back(new OpCalculateDisplacement<SPACE_DIM>(x_ptr, X_ptr, u_ptr));

    // Calculate P


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

    pip.push_back(new OpCalculateDeformationGradient<SPACE_DIM, SPACE_DIM>(

        mat_F_ptr, mat_H_tensor_ptr));


    PetscBool is_linear_elasticity = PETSC_TRUE;

    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_linear_elasticity",

                               &is_linear_elasticity, PETSC_NULL);


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

    if (is_linear_elasticity) {

      pip.push_back(new OpCalculatePiola<SPACE_DIM, SPACE_DIM>(

          shear_modulus_G, bulk_modulus_K, mu, lamme_lambda, mat_P_ptr,

          mat_F_ptr));

    } else {

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

      auto det_ptr = boost::make_shared<VectorDouble>();


      pip.push_back(new OpInvertMatrix<SPACE_DIM>(mat_F_ptr, det_ptr, inv_F));


      pip.push_back(new OpCalculatePiolaIncompressibleNH<SPACE_DIM, SPACE_DIM>(

          shear_modulus_G, bulk_modulus_K, mu, lamme_lambda, mat_P_ptr,

          mat_F_ptr, inv_F, det_ptr));

    }


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

    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(

        "V", mat_v_grad_ptr));


    return boost::make_tuple(v_ptr, X_ptr, x_ptr, mat_P_ptr, mat_F_ptr, u_ptr);

  };


  auto post_proc_boundary = [&]() {

    auto boundary_post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);


    AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(

        boundary_post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");

    auto op_loop_side =

        new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);

    // push ops to side element, through op_loop_side operator

    auto [boundary_v_ptr, boundary_X_ptr, boundary_x_ptr, boundary_mat_P_ptr,

          boundary_mat_F_ptr, boundary_u_ptr] =

        calculate_stress_ops(op_loop_side->getOpPtrVector());

    boundary_post_proc_fe->getOpPtrVector().push_back(op_loop_side);


    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;


    boundary_post_proc_fe->getOpPtrVector().push_back(


        new OpPPMap(


            boundary_post_proc_fe->getPostProcMesh(),

            boundary_post_proc_fe->getMapGaussPts(),


            OpPPMap::DataMapVec{},


            OpPPMap::DataMapMat{{"V", boundary_v_ptr},

                                {"GEOMETRY", boundary_X_ptr},

                                {"x", boundary_x_ptr},

                                {"U", boundary_u_ptr}},


            OpPPMap::DataMapMat{{"FIRST_PIOLA", boundary_mat_P_ptr},

                                {"F", boundary_mat_F_ptr}},


            OpPPMap::DataMapMat{}


            )


    );

    return boundary_post_proc_fe;

  };


  // Add monitor to time solver


  double rho = 1.;

  CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-density", &rho, PETSC_NULL);

  auto get_rho = [rho](const double, const double, const double) {

    return rho;

  };


  SmartPetscObj<Mat> M;   ///< Mass matrix

  SmartPetscObj<KSP> ksp; ///< Linear solver


  auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();

  tsPrePostProc = ts_pre_post_proc;


  CHKERR DMCreateMatrix_MoFEM(dm, M);

  CHKERR MatZeroEntries(M);


  boost::shared_ptr<DomainEle> vol_mass_ele(new DomainEle(mField));


  vol_mass_ele->B = M;


  auto integration_rule = [](int, int, int approx_order) {

    return 2 * approx_order;

  };


  vol_mass_ele->getRuleHook = integration_rule;


  vol_mass_ele->getOpPtrVector().push_back(new OpMassV("V", "V", get_rho));

  vol_mass_ele->getOpPtrVector().push_back(new OpMassF("F", "F"));


  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), vol_mass_ele);

  CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);

  CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);


  auto lumpVec = createDMVector(simple->getDM());

  CHKERR MatGetRowSum(M, lumpVec);


  CHKERR MatZeroEntries(M);

  CHKERR MatDiagonalSet(M, lumpVec, INSERT_VALUES);


  // Create and septup KSP (linear solver), we need this to calculate g(t,u) =

  // M^-1G(t,u)

  ksp = createKSP(mField.get_comm());

  CHKERR KSPSetOperators(ksp, M, M);

  CHKERR KSPSetFromOptions(ksp);

  CHKERR KSPSetUp(ksp);


  auto solve_boundary_for_g = [&]() {

    MoFEMFunctionBegin;

    if (*(pipeline_mng->getBoundaryExplicitRhsFE()->vecAssembleSwitch)) {


      CHKERR VecGhostUpdateBegin(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,

                                 ADD_VALUES, SCATTER_REVERSE);

      CHKERR VecGhostUpdateEnd(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,

                               ADD_VALUES, SCATTER_REVERSE);

      CHKERR VecAssemblyBegin(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);

      CHKERR VecAssemblyEnd(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);

      *(pipeline_mng->getBoundaryExplicitRhsFE()->vecAssembleSwitch) = false;


      auto D =

          smartVectorDuplicate(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);

      CHKERR KSPSolve(ksp, pipeline_mng->getBoundaryExplicitRhsFE()->ts_F, D);

      CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);

      CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);

      CHKERR VecCopy(D, pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);

    }


    MoFEMFunctionReturn(0);

  };


  pipeline_mng->getBoundaryExplicitRhsFE()->postProcessHook =

      solve_boundary_for_g;


  MoFEM::SmartPetscObj<TS> ts;

  ts = pipeline_mng->createTSEX(dm);


  // Field eval

  PetscBool field_eval_flag = PETSC_TRUE;

  boost::shared_ptr<MatrixDouble> velocity_field_ptr;

  boost::shared_ptr<MatrixDouble> geometry_field_ptr;

  boost::shared_ptr<MatrixDouble> spatial_position_field_ptr;

  boost::shared_ptr<SetPtsData> field_eval_data;


  std::array<double, 3> field_eval_coords = {0.5, 0.5, 5.};

  int dim = 3;

  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",

                                  field_eval_coords.data(), &dim,

                                  &field_eval_flag);


  if (field_eval_flag) {

    field_eval_data =

        mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();

    if (SPACE_DIM == 3) {

      CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(

          field_eval_data, simple->getDomainFEName());

    } else {

      CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree2D(

          field_eval_data, simple->getDomainFEName());

    }


    field_eval_data->setEvalPoints(field_eval_coords.data(), 1);


    auto no_rule = [](int, int, int) { return -1; };


    auto fe_ptr = field_eval_data->feMethodPtr.lock();

    fe_ptr->getRuleHook = no_rule;

    velocity_field_ptr = boost::make_shared<MatrixDouble>();

    geometry_field_ptr = boost::make_shared<MatrixDouble>();

    spatial_position_field_ptr = boost::make_shared<MatrixDouble>();

    fe_ptr->getOpPtrVector().push_back(

        new OpCalculateVectorFieldValues<SPACE_DIM>("V", velocity_field_ptr));

    fe_ptr->getOpPtrVector().push_back(

        new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY",

                                                    geometry_field_ptr));

    fe_ptr->getOpPtrVector().push_back(

        new OpCalculateVectorFieldValues<SPACE_DIM>(

            "x_2", spatial_position_field_ptr));

  }


  auto post_proc_domain = [&]() {

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


    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;


    auto [boundary_v_ptr, boundary_X_ptr, boundary_x_ptr, boundary_mat_P_ptr,

          boundary_mat_F_ptr, boundary_u_ptr] =

        calculate_stress_ops(post_proc_fe_vol->getOpPtrVector());


    post_proc_fe_vol->getOpPtrVector().push_back(


        new OpPPMap(


            post_proc_fe_vol->getPostProcMesh(),

            post_proc_fe_vol->getMapGaussPts(),


            {},


            {{"V", boundary_v_ptr},

             {"GEOMETRY", boundary_X_ptr},

             {"x", boundary_x_ptr},

             {"U", boundary_u_ptr}},


            {{"FIRST_PIOLA", boundary_mat_P_ptr}, {"F", boundary_mat_F_ptr}},


            {}


            )


    );

    return post_proc_fe_vol;

  };


  boost::shared_ptr<FEMethod> null_fe;

  auto monitor_ptr = boost::make_shared<Monitor>(

      SmartPetscObj<DM>(dm, true), mField, post_proc_domain(),

      post_proc_boundary(), velocity_field_ptr, spatial_position_field_ptr,

      geometry_field_ptr, field_eval_coords, field_eval_data);


  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,

                             null_fe, monitor_ptr);


  double ftime = 1;

  // CHKERR TSSetMaxTime(ts, ftime);

  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);


  auto T = createDMVector(simple->getDM());

  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,

                                 SCATTER_FORWARD);

  CHKERR TSSetSolution(ts, T);

  CHKERR TSSetFromOptions(ts);


  auto fb = mField.getInterface<FieldBlas>();


  CHKERR TSSetPostStage(ts, TSPrePostProc::tsPostStage);

  CHKERR TSSetPostStep(ts, TSPrePostProc::tsPostStep);

  CHKERR TSSetPreStep(ts, TSPrePostProc::tsPreStep);


  boost::shared_ptr<ForcesAndSourcesCore> null;


  if (auto ptr = tsPrePostProc.lock()) {

    ptr->fsRawPtr = this;

    CHKERR TSSetUp(ts);

    CHKERR TSSolve(ts, NULL);

    CHKERR TSGetTime(ts, &ftime);

  }


  MoFEMFunctionReturn(0);

}

//! [Solve]


//! [Postprocess results]

MoFEMErrorCode Example::outputResults() {

  MoFEMFunctionBegin;

  PetscBool test_flg = PETSC_FALSE;

  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);

  if (test_flg) {

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

    auto T = createDMVector(simple->getDM());

    CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,

                                   SCATTER_FORWARD);

    double nrm2;

    CHKERR VecNorm(T, NORM_2, &nrm2);

    MOFEM_LOG("EXAMPLE", Sev::inform) << "Regression norm " << nrm2;

    constexpr double regression_value = 0.0194561;

    if (fabs(nrm2 - regression_value) > 1e-2)

      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,

              "Regression test failed; wrong norm value.");

  }

  MoFEMFunctionReturn(0);

}

//! [Postprocess results]


//! [Check]

MoFEMErrorCode Example::checkResults() {

  MoFEMFunctionBegin;

  MoFEMFunctionReturn(0);

}

//! [Check]


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


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


  // Initialisation of MoFEM/PETSc and MOAB data structures

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

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


  // Add logging channel for example

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

  core_log->add_sink(

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

  LogManager::setLog("EXAMPLE");

  MOFEM_LOG_TAG("EXAMPLE", "example");


  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]


    //! [Example]

    Example ex(m_field);

    CHKERR ex.runProblem();

    //! [Example]

  }

  CATCH_ERRORS;


  CHKERR MoFEM::Core::Finalize();

}

third
constexpr double third
Definition EshelbianADOL-C.cpp:17

UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition HookeElement.hpp:75

MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition LogManager.hpp:352

MatrixFunction.hpp

MoFEM.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

SPACE_DIM
constexpr int SPACE_DIM
Definition adolc_plasticity.cpp:14

SPACE_DIM
constexpr int SPACE_DIM
Definition child_and_parent.cpp:16

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

DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
Definition child_and_parent.cpp:34

BiLinearForm

BoundaryEleOp

DomainEleOp

FEMethod

FTensor::Index
Definition Index.hpp:24

FTensor::Kronecker_Delta
Kronecker Delta class.
Definition Kronecker_Delta.hpp:15

FTensor::Tensor1
Definition Tensor1_value.hpp:9

FTensor::Tensor2
Definition Tensor2_value.hpp:17

LinearForm

OpCalculateVectorFieldGradient

Range

double

QUIET
@ QUIET
Definition definitions.h:221

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

MF_EXIST
@ MF_EXIST
Definition definitions.h:113

FieldApproximationBase
FieldApproximationBase
approximation base
Definition definitions.h:58

LASTBASE
@ LASTBASE
Definition definitions.h:69

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

DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66

H1
@ H1
continuous field
Definition definitions.h:85

NOSPACE
@ NOSPACE
Definition definitions.h:83

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

MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40

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

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

order
constexpr int order
Definition dg_projection.cpp:18

PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition dynamic_first_order_con_law.cpp:53

tsPrePostProc
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
Definition dynamic_first_order_con_law.cpp:405

OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition dynamic_first_order_con_law.cpp:62

lamme_lambda
double lamme_lambda
Definition dynamic_first_order_con_law.cpp:99

OpRhsTestPiola
FormsIntegrators< DomainEleOp >::Assembly< AssemblyType::PETSC >::LinearForm< IntegrationType::GAUSS >::OpBaseTimesVector< 1, SPACE_DIM *SPACE_DIM, 1 > OpRhsTestPiola
Definition dynamic_first_order_con_law.cpp:78

OpMassV
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMassV
Definition dynamic_first_order_con_law.cpp:57

SPACE_DIM
constexpr int SPACE_DIM
Definition dynamic_first_order_con_law.cpp:24

OpGradTimesTensor2
FormsIntegrators< DomainEleOp >::Assembly< AssemblyType::PETSC >::LinearForm< IntegrationType::GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpGradTimesTensor2
Definition dynamic_first_order_con_law.cpp:74

SetPtsData
FieldEvaluatorInterface::SetPtsData SetPtsData
Definition dynamic_first_order_con_law.cpp:37

mu
double mu
Definition dynamic_first_order_con_law.cpp:98

DomainNaturalBC
NaturalBC< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS > DomainNaturalBC
Definition dynamic_first_order_con_law.cpp:68

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

OpBodyForceVector
DomainNaturalBC::OpFlux< NaturalMeshsetTypeVectorScaling< BLOCKSET >, 1, SPACE_DIM > OpBodyForceVector
Definition dynamic_first_order_con_law.cpp:70

omega
constexpr double omega
Save field DOFS on vertices/tags.
Definition dynamic_first_order_con_law.cpp:93

PostProcEleBdy
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
Definition dynamic_first_order_con_law.cpp:55

PostProcFaceEle
PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore > PostProcFaceEle
Definition dynamic_first_order_con_law.cpp:31

OpGradTimesPiola
FormsIntegrators< DomainEleOp >::Assembly< AssemblyType::PETSC >::LinearForm< IntegrationType::GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpGradTimesPiola
Definition dynamic_first_order_con_law.cpp:82

trace
double trace(FTensor::Tensor2< T, 2, 2 > &t_stress)
Definition dynamic_first_order_con_law.cpp:14

bulk_modulus_K
double bulk_modulus_K
Definition dynamic_first_order_con_law.cpp:96

shear_modulus_G
double shear_modulus_G
Definition dynamic_first_order_con_law.cpp:97

tsPrePostProc
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
Definition free_surface.cpp:467

integration_rule
auto integration_rule
Definition free_surface.cpp:185

t_kd
constexpr auto t_kd
Definition free_surface.cpp:137

MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:523

MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition DMMoFEM.cpp:1197

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

MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:586

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

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

MoFEM::PETSC
@ PETSC
Definition FormsIntegrators.hpp:105

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::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:578

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

dt
double dt
Definition heat_method.cpp:26

D
double D
Definition initial_diffusion.cpp:44

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

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::createKSP
auto createKSP(MPI_Comm comm)
Definition PetscSmartObj.hpp:261

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:1056

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

MoFEM::smartVectorDuplicate
DEPRECATED SmartPetscObj< Vec > smartVectorDuplicate(Vec vec)
Definition PetscSmartObj.hpp:230

MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition DeprecatedPetsc.hpp:152

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

MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition DeprecatedPetsc.hpp:192

MoFEM::getFTensor1FromMat
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
Definition Templates.hpp:238

MoFEM::getFTensor2FromMat
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
Definition Templates.hpp:278

MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition Templates.hpp:135

MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition DeprecatedPetsc.hpp:203

PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition PlasticOps.hpp:117

OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition operators_tests.cpp:47

save_every_nth_step
int save_every_nth_step
Definition photon_diffusion.cpp:67

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

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

rho
double rho
Definition plastic.cpp:140

EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition plastic.cpp:13

poisson_ratio
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:122

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

approx_order
static constexpr int approx_order
Definition prism_polynomial_approximation.cpp:14

OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition seepage.cpp:56

mu
constexpr double mu
Definition shallow_wave.cpp:64

omega
constexpr double omega
Definition shallow_wave.cpp:62

i
FTensor::Index< 'i', 3 > i
Definition shallow_wave.cpp:77

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

AddHOOps
Definition EshelbianOperators.hpp:597

CommonData
Definition continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22

CommonData::M
SmartPetscObj< Mat > M
Mass matrix.
Definition dynamic_first_order_con_law.cpp:415

CommonData::ksp
SmartPetscObj< KSP > ksp
Linear solver.
Definition dynamic_first_order_con_law.cpp:416

Example::DynamicFirstOrderConsConstantTimeScale::getScale
double getScale(const double time)
Get scaling at given time.
Definition dynamic_first_order_con_law.cpp:444

Example::DynamicFirstOrderConsSinusTimeScale::getScale
double getScale(const double time)
Get scaling at given time.
Definition dynamic_first_order_con_law.cpp:439

Example
[Example]
Definition plastic.cpp:213

Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition dynamic_first_order_con_law.cpp:536

Example::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition dynamic_first_order_con_law.cpp:647

Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition dynamic_first_order_con_law.cpp:463

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

Example::checkResults
MoFEMErrorCode checkResults()
[Postprocess results]
Definition dynamic_first_order_con_law.cpp:1205

Example::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition dynamic_first_order_con_law.cpp:893

Example::Example
Example(MoFEM::Interface &m_field)
Definition plastic.cpp:215

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

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

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

Example::outputResults
MoFEMErrorCode outputResults()
[Solve]
Definition dynamic_first_order_con_law.cpp:1183

LinMomTimeScale
Definition dynamic_first_order_con_law.cpp:407

LinMomTimeScale::getScale
double getScale(const double time)
Get scaling at given time.
Definition dynamic_first_order_con_law.cpp:409

MoFEM::BasicMethod::getCacheWeakPtr
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak ptr object.
Definition LoopMethods.hpp:344

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

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

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

MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() 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:76

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

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

MoFEM::FieldBlas
Basic algebra on fields.
Definition FieldBlas.hpp:21

MoFEM::FieldEvaluatorInterface::SetPtsData
Definition FieldEvaluator.hpp:29

MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition FieldEvaluator.hpp:21

MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition ForcesAndSourcesCore.hpp:22

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

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

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

MoFEM::OpCalculateTensor2FieldValues
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Definition UserDataOperators.hpp:889

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

MoFEM::OpInvertMatrix
Definition UserDataOperators.hpp:3249

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

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

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

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

MoFEM::PostProcBrokenMeshInMoabBaseContImpl
Definition PostProcBrokenMeshInMoabBase.hpp:884

MoFEM::PostProcBrokenMeshInMoab
Definition PostProcBrokenMeshInMoabBase.hpp:667

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:83

MoFEM::TSMethod::ts_step
PetscInt ts_step
time step number
Definition LoopMethods.hpp:159

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:137

MoFEM::TimeScale::TimeScale
TimeScale(std::string file_name="", bool error_if_file_not_given=false, ScalingFun def_scaling_fun=[](double time) { return time;})
TimeScale constructor.
Definition ScalingMethod.cpp:22

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

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

Monitor::postProcBdy
boost::shared_ptr< PostProcFaceEle > postProcBdy
Definition dynamic_first_order_con_law.cpp:884

Monitor::fieldEvalCoords
std::array< double, 3 > fieldEvalCoords
Definition dynamic_first_order_con_law.cpp:888

Monitor::mField
MoFEM::Interface & mField
Definition dynamic_first_order_con_law.cpp:784

Monitor::Monitor
Monitor(SmartPetscObj< DM > dm, MoFEM::Interface &m_field, boost::shared_ptr< PostProcEle > post_proc, boost::shared_ptr< PostProcFaceEle > post_proc_bdry, boost::shared_ptr< MatrixDouble > velocity_field_ptr, boost::shared_ptr< MatrixDouble > x2_field_ptr, boost::shared_ptr< MatrixDouble > geometry_field_ptr, std::array< double, 3 > pass_field_eval_coords, boost::shared_ptr< SetPtsData > pass_field_eval_data)
Definition dynamic_first_order_con_law.cpp:785

Monitor::geometryFieldPtr
boost::shared_ptr< MatrixDouble > geometryFieldPtr
Definition dynamic_first_order_con_law.cpp:887

Monitor::dM
SmartPetscObj< DM > dM
Definition dynamic_first_order_con_law.cpp:882

Monitor::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition dynamic_first_order_con_law.cpp:798

Monitor::velocityFieldPtr
boost::shared_ptr< MatrixDouble > velocityFieldPtr
Definition dynamic_first_order_con_law.cpp:885

Monitor::fieldEvalData
boost::shared_ptr< SetPtsData > fieldEvalData
Definition dynamic_first_order_con_law.cpp:889

Monitor::x2FieldPtr
boost::shared_ptr< MatrixDouble > x2FieldPtr
Definition dynamic_first_order_con_law.cpp:886

Monitor::postProc
boost::shared_ptr< PostProcEle > postProc
Definition dynamic_first_order_con_law.cpp:883

OpCalculateDeformationGradient
Definition dynamic_first_order_con_law.cpp:340

OpCalculateDeformationGradient::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition dynamic_first_order_con_law.cpp:347

OpCalculateDeformationGradient::defGradPtr
boost::shared_ptr< MatrixDouble > defGradPtr
Definition dynamic_first_order_con_law.cpp:380

OpCalculateDeformationGradient::gradTensorPtr
boost::shared_ptr< MatrixDouble > gradTensorPtr
Definition dynamic_first_order_con_law.cpp:379

OpCalculateDeformationGradient::OpCalculateDeformationGradient
OpCalculateDeformationGradient(boost::shared_ptr< MatrixDouble > def_grad_ptr, boost::shared_ptr< MatrixDouble > grad_tensor_ptr)
Definition dynamic_first_order_con_law.cpp:341

OpCalculateDisplacement
Definition dynamic_first_order_con_law.cpp:226

OpCalculateDisplacement::XPtr
boost::shared_ptr< MatrixDouble > XPtr
Definition dynamic_first_order_con_law.cpp:262

OpCalculateDisplacement::uPtr
boost::shared_ptr< MatrixDouble > uPtr
Definition dynamic_first_order_con_law.cpp:263

OpCalculateDisplacement::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition dynamic_first_order_con_law.cpp:233

OpCalculateDisplacement::OpCalculateDisplacement
OpCalculateDisplacement(boost::shared_ptr< MatrixDouble > spatial_pos_ptr, boost::shared_ptr< MatrixDouble > reference_pos_ptr, boost::shared_ptr< MatrixDouble > u_ptr)
Definition dynamic_first_order_con_law.cpp:227

OpCalculateDisplacement::xPtr
boost::shared_ptr< MatrixDouble > xPtr
Definition dynamic_first_order_con_law.cpp:261

OpCalculateFStab
Definition dynamic_first_order_con_law.cpp:104

OpCalculateFStab::gradxPtr
boost::shared_ptr< MatrixDouble > gradxPtr
Definition dynamic_first_order_con_law.cpp:162

OpCalculateFStab::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition dynamic_first_order_con_law.cpp:116

OpCalculateFStab::OpCalculateFStab
OpCalculateFStab(boost::shared_ptr< MatrixDouble > def_grad_ptr, boost::shared_ptr< MatrixDouble > def_grad_stab_ptr, boost::shared_ptr< MatrixDouble > def_grad_dot_ptr, double tau_F_ptr, double xi_F_ptr, boost::shared_ptr< MatrixDouble > grad_x_ptr, boost::shared_ptr< MatrixDouble > grad_vel_ptr)
Definition dynamic_first_order_con_law.cpp:105

OpCalculateFStab::defGradStabPtr
boost::shared_ptr< MatrixDouble > defGradStabPtr
Definition dynamic_first_order_con_law.cpp:160

OpCalculateFStab::gradVelPtr
boost::shared_ptr< MatrixDouble > gradVelPtr
Definition dynamic_first_order_con_law.cpp:163

OpCalculateFStab::tauFPtr
double tauFPtr
Definition dynamic_first_order_con_law.cpp:157

OpCalculateFStab::defGradPtr
boost::shared_ptr< MatrixDouble > defGradPtr
Definition dynamic_first_order_con_law.cpp:159

OpCalculateFStab::xiF
double xiF
Definition dynamic_first_order_con_law.cpp:158

OpCalculateFStab::defGradDotPtr
boost::shared_ptr< MatrixDouble > defGradDotPtr
Definition dynamic_first_order_con_law.cpp:161

OpCalculatePiola
Definition dynamic_first_order_con_law.cpp:168

OpCalculatePiola::mU
double mU
Definition dynamic_first_order_con_law.cpp:219

OpCalculatePiola::OpCalculatePiola
OpCalculatePiola(double shear_modulus, double bulk_modulus, double m_u, double lambda_lamme, boost::shared_ptr< MatrixDouble > first_piola_ptr, boost::shared_ptr< MatrixDouble > def_grad_ptr)
Definition dynamic_first_order_con_law.cpp:169

OpCalculatePiola::shearModulus
double shearModulus
Definition dynamic_first_order_con_law.cpp:217

OpCalculatePiola::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition dynamic_first_order_con_law.cpp:178

OpCalculatePiola::lammeLambda
double lammeLambda
Definition dynamic_first_order_con_law.cpp:220

OpCalculatePiola::defGradPtr
boost::shared_ptr< MatrixDouble > defGradPtr
Definition dynamic_first_order_con_law.cpp:222

OpCalculatePiola::firstPiolaPtr
boost::shared_ptr< MatrixDouble > firstPiolaPtr
Definition dynamic_first_order_con_law.cpp:221

OpCalculatePiola::bulkModulus
double bulkModulus
Definition dynamic_first_order_con_law.cpp:218

OpCalculatePiolaIncompressibleNH
Definition dynamic_first_order_con_law.cpp:268

OpCalculatePiolaIncompressibleNH::OpCalculatePiolaIncompressibleNH
OpCalculatePiolaIncompressibleNH(double shear_modulus, double bulk_modulus, double m_u, double lambda_lamme, boost::shared_ptr< MatrixDouble > first_piola_ptr, boost::shared_ptr< MatrixDouble > def_grad_ptr, boost::shared_ptr< MatrixDouble > inv_def_grad_ptr, boost::shared_ptr< VectorDouble > det)
Definition dynamic_first_order_con_law.cpp:269

OpCalculatePiolaIncompressibleNH::dEt
boost::shared_ptr< VectorDouble > dEt
Definition dynamic_first_order_con_law.cpp:335

OpCalculatePiolaIncompressibleNH::bulkModulus
double bulkModulus
Definition dynamic_first_order_con_law.cpp:329

OpCalculatePiolaIncompressibleNH::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition dynamic_first_order_con_law.cpp:280

OpCalculatePiolaIncompressibleNH::mU
double mU
Definition dynamic_first_order_con_law.cpp:330

OpCalculatePiolaIncompressibleNH::invDefGradPtr
boost::shared_ptr< MatrixDouble > invDefGradPtr
Definition dynamic_first_order_con_law.cpp:334

OpCalculatePiolaIncompressibleNH::shearModulus
double shearModulus
Definition dynamic_first_order_con_law.cpp:328

OpCalculatePiolaIncompressibleNH::defGradPtr
boost::shared_ptr< MatrixDouble > defGradPtr
Definition dynamic_first_order_con_law.cpp:333

OpCalculatePiolaIncompressibleNH::firstPiolaPtr
boost::shared_ptr< MatrixDouble > firstPiolaPtr
Definition dynamic_first_order_con_law.cpp:332

OpCalculatePiolaIncompressibleNH::lammeLambda
double lammeLambda
Definition dynamic_first_order_con_law.cpp:331

PostProcEleByDim
Definition dynamic_first_order_con_law.cpp:39

TSPrePostProc
Set of functions called by PETSc solver used to refine and update mesh.
Definition dynamic_first_order_con_law.cpp:384

TSPrePostProc::tsPostStep
static MoFEMErrorCode tsPostStep(TS ts)
Definition dynamic_first_order_con_law.cpp:614

TSPrePostProc::~TSPrePostProc
virtual ~TSPrePostProc()=default

TSPrePostProc::fsRawPtr
Example * fsRawPtr
Definition dynamic_first_order_con_law.cpp:398

TSPrePostProc::tsPreStep
static MoFEMErrorCode tsPreStep(TS ts)
Definition dynamic_first_order_con_law.cpp:628

TSPrePostProc::TSPrePostProc
TSPrePostProc()=default

TSPrePostProc::tsPostStage
static MoFEMErrorCode tsPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
[Boundary condition]
Definition dynamic_first_order_con_law.cpp:581

TSPrePostProc::tsSetUp
MoFEMErrorCode tsSetUp(TS ts)
Used to setup TS solver.
Definition free_surface.cpp:3142

OpForce
BoundaryNaturalBC::OpFlux< NaturalForceMeshsets, 1, SPACE_DIM > OpForce
Definition thermo_elastic.cpp:112

OpBodyForce
DomainNaturalBCRhs::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, SPACE_DIM > OpBodyForce
Definition thermo_elastic.cpp:101

BoundaryNaturalBC
NaturalBC< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS > BoundaryNaturalBC
[Body and heat source]
Definition thermo_elastic.cpp:110