Incompressible Struct Reference

Collaboration diagram for Incompressible:

Public Member Functions
	Incompressible (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	[Run problem] More...

Private Member Functions
MoFEMErrorCode	setupProblem ()
	[Run problem] More...
MoFEMErrorCode	createCommonData ()
	[Set up problem] More...
MoFEMErrorCode	bC ()
	[Create common data] More...
MoFEMErrorCode	OPs ()
	[Boundary condition] More...
MoFEMErrorCode	tsSolve ()
MoFEMErrorCode	checkResults ()
	[Solve] More...

Private Attributes
MoFEM::Interface &	mField
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > >	uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > >	uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > >	uZScatter

Detailed Description

Definition at line 203 of file incompressible_elasticity.cpp.

Constructor & Destructor Documentation

Incompressible()

Incompressible::Incompressible ( MoFEM::Interface &  m_field )

inline

Definition at line 205 of file incompressible_elasticity.cpp.

: mField(m_field) {}

Member Function Documentation

bC()

MoFEMErrorCode Incompressible::bC ( )

private

[Create common data]

[Boundary condition]

[Natural boundary condition]

[Define gravity vector]

Definition at line 391 of file incompressible_elasticity.cpp.

                                  {
  MoFEMFunctionBegin;
 
  auto pipeline_mng = mField.getInterface<PipelineManager>();
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
  auto time_scale = boost::make_shared<TimeScale>();
 
  auto integration_rule = [](int, int, int approx_order) {
    return 2 * (approx_order - 1);
  };
 
  using DomainNaturalBC =
      NaturalBC<DomainEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
  using OpBodyForce =
      DomainNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
  using BoundaryNaturalBC =
      NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
  using OpForce = BoundaryNaturalBC::OpFlux<NaturalForceMeshsets, 1, SPACE_DIM>;
  
  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pipeline_mng->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
  //! [Natural boundary condition]
  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
      pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
      "FORCE", Sev::inform);
 
  //! [Define gravity vector]
  CHKERR DomainNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
      pipeline_mng->getOpDomainRhsPipeline(), mField, "U", {time_scale},
      "BODY_FORCE", Sev::inform);
 
  // Essential BC
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
                                           "U", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
                                           "U", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
                                           "U", 2, 2);
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      simple->getProblemName(), "U");
 
  MoFEMFunctionReturn(0);
}

checkResults()

MoFEMErrorCode Incompressible::checkResults ( )

private

[Solve]

[Check]

Definition at line 880 of file incompressible_elasticity.cpp.

{ return 0; }

createCommonData()

MoFEMErrorCode Incompressible::createCommonData ( )

private

[Set up problem]

[Create common data]

Definition at line 361 of file incompressible_elasticity.cpp.

                                                {
  MoFEMFunctionBegin;
 
  auto get_options = [&]() {
    MoFEMFunctionBegin;
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
                                 &young_modulus, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
                                 &poisson_ratio, PETSC_NULL);
 
    MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
        << "Young modulus " << young_modulus;
    MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
        << "Poisson_ratio " << poisson_ratio;
 
    mu = young_modulus / (2. + 2. * poisson_ratio);
    const double lambda_denom =
        (1. + poisson_ratio) * (1. - 2. * poisson_ratio);
    lambda = young_modulus * poisson_ratio / lambda_denom;
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR get_options();
 
  MoFEMFunctionReturn(0);
}

OPs()

MoFEMErrorCode Incompressible::OPs ( )

private

[Boundary condition]

[Push operators to pip]

[Operators used for RHS incompressible elasticity]

[Operators used for incompressible elasticity]

[Operators used for incompressible elasticity]

[Operators used for RHS incompressible elasticity]

[Operators used for RHS incompressible elasticity]

Definition at line 439 of file incompressible_elasticity.cpp.

                                   {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pip_mng = mField.getInterface<PipelineManager>();
  auto bc_mng = mField.getInterface<BcManager>();
 
  auto integration_rule_vol = [](int, int, int approx_order) {
    return 2 * approx_order + geom_order - 1;
  };
 
  auto add_domain_base_ops = [&](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, L2},
                                                          "GEOMETRY");
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    // This assemble A-matrix
    using OpMassPressure = FormsIntegrators<DomainEleOp>::Assembly<
        AT>::BiLinearForm<GAUSS>::OpMass<1, 1>;
    // This assemble B-matrix (preconditioned)
    using OpMassPressureStab = FormsIntegrators<DomainEleOpStab>::Assembly<
        AT>::BiLinearForm<GAUSS>::OpMass<1, 1>;
    //! [Operators used for RHS incompressible elasticity]
 
    //! [Operators used for incompressible elasticity]
    using OpGradSymTensorGrad = FormsIntegrators<DomainEleOp>::Assembly<
        AT>::BiLinearForm<IT>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
    using OpMixScalarTimesDiv = FormsIntegrators<DomainEleOp>::Assembly<
        AT>::BiLinearForm<IT>::OpMixScalarTimesDiv<SPACE_DIM, coord_type>;
    //! [Operators used for incompressible elasticity]
 
    auto mat_D_ptr = boost::make_shared<MatrixDouble>();
    constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
    mat_D_ptr->resize(size_symm * size_symm, 1);
 
    FTensor::Index<'i', SPACE_DIM> i;
    FTensor::Index<'j', SPACE_DIM> j;
    FTensor::Index<'k', SPACE_DIM> k;
    FTensor::Index<'l', SPACE_DIM> l;
 
    constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
    auto t_mat = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
    t_mat(i, j, k, l) = -2. * mu * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
 
    pip.push_back(new OpMixScalarTimesDiv(
        "P", "U",
        [](const double, const double, const double) constexpr { return -1.; },
        true, false));
    pip.push_back(new OpGradSymTensorGrad("U", "U", mat_D_ptr));
 
    auto get_lambda_reciprocal = [](const double, const double, const double) {
      return 1. / lambda;
    };
    if (poisson_ratio < 0.5)
      pip.push_back(new OpMassPressure("P", "P", get_lambda_reciprocal));
 
    double eps_stab = 0;
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
                                 PETSC_NULL);
    if (eps_stab > 0)
      pip.push_back(new OpMassPressureStab(
          "P", "P", [eps_stab](double, double, double) { return eps_stab; }));
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    //! [Operators used for RHS incompressible elasticity]
    using OpDomainGradTimesTensor = FormsIntegrators<DomainEleOp>::Assembly<
        AT>::LinearForm<GAUSS>::OpGradTimesSymTensor<1, SPACE_DIM, SPACE_DIM>;
 
    using OpDivDeltaUTimesP = FormsIntegrators<DomainEleOp>::Assembly<
        AT>::LinearForm<GAUSS>::OpMixDivTimesU<1, SPACE_DIM, SPACE_DIM>;
 
    using OpBaseTimesScalarValues = FormsIntegrators<DomainEleOp>::Assembly<
        AT>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
 
    //! [Operators used for RHS incompressible elasticity]
 
    auto pressure_ptr = boost::make_shared<VectorDouble>();
    pip.push_back(new OpCalculateScalarFieldValues("P", pressure_ptr));
 
    auto div_u_ptr = boost::make_shared<VectorDouble>();
    pip.push_back(
        new OpCalculateDivergenceVectorFieldValues<SPACE_DIM>("U", div_u_ptr));
 
    auto grad_u_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "U", grad_u_ptr));
 
    auto strain_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(grad_u_ptr, strain_ptr));
 
    auto get_four_mu = [](const double, const double, const double) {
      return -2. * mu;
    };
    auto minus_one = [](const double, const double, const double) constexpr {
      return -1.;
    };
 
    pip.push_back(new OpDomainGradTimesTensor("U", strain_ptr, get_four_mu));
 
    pip.push_back(new OpDivDeltaUTimesP("U", pressure_ptr, minus_one));
 
    pip.push_back(new OpBaseTimesScalarValues("P", div_u_ptr, minus_one));
 
    auto get_lambda_reciprocal = [](const double, const double, const double) {
      return 1. / lambda;
    };
    if (poisson_ratio < 0.5) {
      pip.push_back(new OpBaseTimesScalarValues("P", pressure_ptr,
                                                get_lambda_reciprocal));
    }
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
  CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
 
  CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
  CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
 
  MoFEMFunctionReturn(0);
}

runProblem()

MoFEMErrorCode Incompressible::runProblem

(

)

[Run problem]

Definition at line 273 of file incompressible_elasticity.cpp.

                                          {
  MoFEMFunctionBegin;
  CHKERR setupProblem();
  CHKERR createCommonData();
  CHKERR bC();
  CHKERR OPs();
  CHKERR tsSolve();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}

setupProblem()

MoFEMErrorCode Incompressible::setupProblem ( )

private

[Run problem]

[Set up problem]

Definition at line 286 of file incompressible_elasticity.cpp.

                                            {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
 
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
                            PETSC_NULL);
  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_discontinuous_pressure",
                             &isDiscontinuousPressure, PETSC_NULL);
 
  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Order " << order;
  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Geom order " << geom_order;
 
  // Select base
  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("INCOMP_ELASTICITY", Sev::inform)
        << "Set AINSWORTH_LEGENDRE_BASE for displacements";
    break;
  case DEMKOWICZ:
    base = DEMKOWICZ_JACOBI_BASE;
    MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
        << "Set DEMKOWICZ_JACOBI_BASE for displacements";
    break;
  default:
    base = LASTBASE;
    break;
  }
 
  // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
  // Demkowicz base. We need to implement Demkowicz H1 base on tet.
  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
 
  // Adding fields related to incompressible elasticity
  // Add displacement domain and boundary fields
  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
  CHKERR simple->setFieldOrder("U", order);
 
  // Add pressure domain and boundary fields
  // Choose either Crouzeix-Raviart element:
  if (isDiscontinuousPressure) {
    CHKERR simple->addDomainField("P", L2, base, 1);
    CHKERR simple->setFieldOrder("P", order - 2);
  } else {
    // ... or Taylor-Hood element:
    CHKERR simple->addDomainField("P", H1, base, 1);
    CHKERR simple->setFieldOrder("P", order - 1);
  }
 
  // Add geometry data field
  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
 
  CHKERR simple->setUp();
 
  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]

tsSolve()

MoFEMErrorCode Incompressible::tsSolve ( )

private

Definition at line 584 of file incompressible_elasticity.cpp.

                                       {
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
  auto pip_mng = mField.getInterface<PipelineManager>();
  auto is_manager = mField.getInterface<ISManager>();
 
  auto set_section_monitor = [&](auto solver) {
    MoFEMFunctionBegin;
    SNES snes;
    CHKERR TSGetSNES(solver, &snes);
    PetscViewerAndFormat *vf;
    CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
                                      PETSC_VIEWER_DEFAULT, &vf);
    CHKERR SNESMonitorSet(
        snes,
        (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
        vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
    MoFEMFunctionReturn(0);
  };
 
  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 create_post_process_elements = [&]() {
    auto pp_fe = boost::make_shared<PostProcEle>(mField);
    auto &pip = pp_fe->getOpPtrVector();
 
    auto push_vol_ops = [this](auto &pip) {
      MoFEMFunctionBegin;
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                            "GEOMETRY");
 
      MoFEMFunctionReturn(0);
    };
 
    auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
      MoFEMFunctionBegin;
 
      auto &pip = pp_fe->getOpPtrVector();
 
      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));
 
      auto pressure_ptr = boost::make_shared<VectorDouble>();
      pip.push_back(new OpCalculateScalarFieldValues("P", pressure_ptr));
 
      auto div_u_ptr = boost::make_shared<VectorDouble>();
      pip.push_back(new OpCalculateDivergenceVectorFieldValues<SPACE_DIM>(
          "U", div_u_ptr));
 
      auto grad_u_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
          "U", grad_u_ptr));
 
      auto strain_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(grad_u_ptr, strain_ptr));
 
      auto stress_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateLameStress<SPACE_DIM>(
          mu, stress_ptr, strain_ptr, pressure_ptr));
 
      pip.push_back(
 
          new OpPPMap(
 
              pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
 
              {{"P", pressure_ptr}},
 
              {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
 
              {},
 
              {{"STRAIN", strain_ptr}, {"STRESS", stress_ptr}}
 
              )
 
      );
 
      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<DomainEle>(
          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());
  };
 
  boost::shared_ptr<SetPtsData> field_eval_data;
  boost::shared_ptr<MatrixDouble> vector_field_ptr;
 
  std::array<double, SPACE_DIM> field_eval_coords;
  int coords_dim = SPACE_DIM;
  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
                                  field_eval_coords.data(), &coords_dim,
                                  &doEvalField);
 
  if (doEvalField) {
    vector_field_ptr = boost::make_shared<MatrixDouble>();
    field_eval_data =
        mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
    if constexpr (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 field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
    field_eval_fe_ptr->getRuleHook = no_rule;
    field_eval_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("U", vector_field_ptr));
  }
 
  auto set_time_monitor = [&](auto dm, auto solver) {
    MoFEMFunctionBegin;
    boost::shared_ptr<MonitorIncompressible> monitor_ptr(
        new MonitorIncompressible(SmartPetscObj<DM>(dm, true),
                                  create_post_process_elements(), uXScatter,
                                  uYScatter, uZScatter, field_eval_coords,
                                  field_eval_data, vector_field_ptr));
    boost::shared_ptr<ForcesAndSourcesCore> null;
    CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
                               monitor_ptr, null, null);
    MoFEMFunctionReturn(0);
  };
 
  auto set_essential_bc = [&]() {
    MoFEMFunctionBegin;
    // This is low level pushing finite elements (pipelines) to solver
    auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
    auto pre_proc_ptr = boost::make_shared<FEMethod>();
    auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
 
    // Add boundary condition scaling
    auto time_scale = boost::make_shared<TimeScale>();
 
    pre_proc_ptr->preProcessHook = EssentialPreProc<DisplacementCubitBcData>(
        mField, pre_proc_ptr, {time_scale}, false);
    post_proc_rhs_ptr->postProcessHook =
        EssentialPostProcRhs<DisplacementCubitBcData>(mField, post_proc_rhs_ptr,
                                                      1.);
    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);
    MoFEMFunctionReturn(0);
  };
 
  auto set_schur_pc = [&](auto solver) {
    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);
    boost::shared_ptr<SetUpSchur> schur_ptr;
    auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
 
    auto A = createDMMatrix(simple->getDM());
    auto B = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
 
    CHK_MOAB_THROW(
        TSSetIJacobian(solver, A, B, TsSetIJacobian, ts_ctx_ptr.get()),
        "set operators");
    auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
    auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
    pre_proc_schur_lhs_ptr->preProcessHook = [pre_proc_schur_lhs_ptr]() {
      MoFEMFunctionBegin;
      MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Zero matrices";
      CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->A);
      CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->B);
      MoFEMFunctionReturn(0);
    };
    post_proc_schur_lhs_ptr->postProcessHook = [this,
                                                post_proc_schur_lhs_ptr]() {
      MoFEMFunctionBegin;
      MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Assemble Begin";
      *(post_proc_schur_lhs_ptr->matAssembleSwitch) = false;
      CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
      CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
      CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
          mField, post_proc_schur_lhs_ptr, 1.,
          SmartPetscObj<Mat>(post_proc_schur_lhs_ptr->A, true))();
 
      CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
      CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
      CHKERR MatAXPY(post_proc_schur_lhs_ptr->B, 1, post_proc_schur_lhs_ptr->A,
                     SAME_NONZERO_PATTERN);
      MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Assemble End";
      MoFEMFunctionReturn(0);
    };
    ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
    ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
 
    if (is_pcfs == PETSC_TRUE) {
      if (AT == AssemblyType::SCHUR) {
        schur_ptr = SetUpSchur::createSetUpSchur(mField);
        CHK_MOAB_THROW(schur_ptr->setUp(solver), "setup schur preconditioner");
      } else {
        auto set_pcfieldsplit_preconditioned_ts = [&](auto solver) {
          MoFEMFunctionBegin;
          auto bc_mng = mField.getInterface<BcManager>();
          auto name_prb = simple->getProblemName();
          SmartPetscObj<IS> is_p;
          CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
              name_prb, ROW, "P", 0, 1, is_p);
          CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_p);
          MoFEMFunctionReturn(0);
        };
        CHK_MOAB_THROW(set_pcfieldsplit_preconditioned_ts(solver),
                       "set pcfieldsplit preconditioned");
      }
      return boost::make_tuple(schur_ptr, A, B);
    }
 
    return boost::make_tuple(schur_ptr, A, B);
  };
 
  auto dm = simple->getDM();
  auto D = createDMVector(dm);
 
  uXScatter = scatter_create(D, 0);
  uYScatter = scatter_create(D, 1);
  if (SPACE_DIM == 3)
    uZScatter = scatter_create(D, 2);
 
  // Add extra finite elements to SNES solver pipelines to resolve essential
  // boundary conditions
  CHKERR set_essential_bc();
 
  auto solver = pip_mng->createTSIM();
  CHKERR TSSetFromOptions(solver);
 
  CHKERR set_section_monitor(solver);
  CHKERR set_time_monitor(dm, solver);
  auto [schur_pc_ptr, A, B] = set_schur_pc(solver);
 
  CHKERR TSSetSolution(solver, D);
  CHKERR TSSetUp(solver);
  CHKERR TSSolve(solver, NULL);
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

mField

MoFEM::Interface& Incompressible::mField

private

Definition at line 210 of file incompressible_elasticity.cpp.

uXScatter

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

private

Definition at line 219 of file incompressible_elasticity.cpp.

uYScatter

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

private

Definition at line 220 of file incompressible_elasticity.cpp.

uZScatter

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

private

Definition at line 221 of file incompressible_elasticity.cpp.

---

The documentation for this struct was generated from the following file:

• users_modules/tutorials/mix-2/incompressible_elasticity.cpp