HeatEquation Struct Reference

Collaboration diagram for HeatEquation:

Public Member Functions
	HeatEquation (MoFEM::Interface &m_field)
MoFEMErrorCode	runProgram ()

Private Member Functions
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	setIntegrationRules ()
MoFEMErrorCode	initialCondition ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()

Private Attributes
MoFEM::Interface &	mField
boost::shared_ptr< std::vector< unsigned char > >	boundaryMarker

Detailed Description

Examples

heat_equation.cpp.

Definition at line 91 of file heat_equation.cpp.

Constructor & Destructor Documentation

HeatEquation()

HeatEquation::HeatEquation

(

MoFEM::Interface &

m_field

)

Examples

heat_equation.cpp.

Definition at line 116 of file heat_equation.cpp.

: mField(m_field) {}

Member Function Documentation

assembleSystem()

MoFEMErrorCode HeatEquation::assembleSystem ( )

private

Examples

heat_equation.cpp.

Definition at line 205 of file heat_equation.cpp.

                                            {
  MoFEMFunctionBegin;
 
  auto add_domain_lhs_ops = [&](auto &pipeline) {
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pipeline, {H1});
    pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
    pipeline.push_back(new OpDomainGradGrad(
        "U", "U", [](double, double, double) -> double { return k; }));
    auto get_c = [this](const double, const double, const double) {
      auto pipeline_mng = mField.getInterface<PipelineManager>();
      auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
      return c * fe_domain_lhs->ts_a;
    };
    pipeline.push_back(new OpDomainMass("U", "U", get_c));
    pipeline.push_back(new OpUnSetBc("U"));
  };
 
  auto add_domain_rhs_ops = [&](auto &pipeline) {
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pipeline, {H1});
    pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
    auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
    auto dot_u_at_gauss_pts = boost::make_shared<VectorDouble>();
    pipeline.push_back(new OpCalculateScalarFieldGradient<SPACE_DIM>(
        "U", grad_u_at_gauss_pts));
    pipeline.push_back(
        new OpCalculateScalarFieldValuesDot("U", dot_u_at_gauss_pts));
    pipeline.push_back(new OpDomainGradTimesVec(
        "U", grad_u_at_gauss_pts,
        [](double, double, double) -> double { return k; }));
    pipeline.push_back(new OpDomainTimesScalarField(
        "U", dot_u_at_gauss_pts,
        [](const double, const double, const double) { return c; }));
    auto source_term = [&](const double x, const double y, const double z) {
      auto pipeline_mng = mField.getInterface<PipelineManager>();
      auto &fe_domain_lhs = pipeline_mng->getDomainRhsFE();
      const auto t = fe_domain_lhs->ts_t;
      return 1e1 * pow(M_E, -M_PI * M_PI * t) * sin(1. * M_PI * x) *
             sin(2. * M_PI * y);
    };
    pipeline.push_back(new OpDomainSource("U", source_term));
    pipeline.push_back(new OpUnSetBc("U"));
  };
 
  auto add_boundary_lhs_ops = [&](auto &pipeline) {
    pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
    pipeline.push_back(new OpBoundaryMass(
        "U", "U", [](const double, const double, const double) { return c; }));
    pipeline.push_back(new OpUnSetBc("U"));
  };
  auto add_boundary_rhs_ops = [&](auto &pipeline) {
    pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
    auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
    auto boundary_function = [&](const double x, const double y,
                                 const double z) {
      auto pipeline_mng = mField.getInterface<PipelineManager>();
      auto &fe_rhs = pipeline_mng->getBoundaryRhsFE();
      const auto t = fe_rhs->ts_t;
      return 0;
      // abs(0.1 * pow(M_E, -M_PI * M_PI * t) * sin(2. * M_PI * x) *
      //     sin(3. * M_PI * y));
    };
    pipeline.push_back(new OpCalculateScalarFieldValues("U", u_at_gauss_pts));
    pipeline.push_back(new OpBoundaryTimeScalarField(
        "U", u_at_gauss_pts,
        [](const double, const double, const double) { return c; }));
    pipeline.push_back(new OpBoundarySource("U", boundary_function));
    pipeline.push_back(new OpUnSetBc("U"));
  };
 
  auto pipeline_mng = mField.getInterface<PipelineManager>();
  add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
  add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
  add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
  add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
 
  MoFEMFunctionReturn(0);
}

boundaryCondition()

MoFEMErrorCode HeatEquation::boundaryCondition ( )

private

Examples

heat_equation.cpp.

Definition at line 183 of file heat_equation.cpp.

                                               {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  auto *simple = mField.getInterface<Simple>();
  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ESSENTIAL",
                                        "U", 0, 0);
 
  auto &bc_map = bc_mng->getBcMapByBlockName();
  boundaryMarker = boost::make_shared<std::vector<char unsigned>>();
  for (auto b : bc_map) {
    if (std::regex_match(b.first, std::regex("(.*)ESSENTIAL(.*)"))) {
      boundaryMarker->resize(b.second->bcMarkers.size(), 0);
      for (int i = 0; i != b.second->bcMarkers.size(); ++i) {
        (*boundaryMarker)[i] |= b.second->bcMarkers[i];
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}

initialCondition()

MoFEMErrorCode HeatEquation::initialCondition ( )

private

Examples

heat_equation.cpp.

Definition at line 161 of file heat_equation.cpp.

                                              {
  MoFEMFunctionBegin;
 
  // Get surface entities form blockset, set initial values in those
  // blocksets. To keep it simple, it is assumed that inital values are on
  // blockset 1
  if (mField.getInterface<MeshsetsManager>()->checkMeshset(1, BLOCKSET)) {
    Range inner_surface;
    CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
        1, BLOCKSET, 2, inner_surface, true);
    if (!inner_surface.empty()) {
      Range inner_surface_verts;
      CHKERR mField.get_moab().get_connectivity(inner_surface,
                                                inner_surface_verts, false);
      CHKERR mField.getInterface<FieldBlas>()->setField(
          init_u, MBVERTEX, inner_surface_verts, "U");
    }
  }
 
  MoFEMFunctionReturn(0);
}

outputResults()

MoFEMErrorCode HeatEquation::outputResults ( )

private

Examples

heat_equation.cpp.

Definition at line 416 of file heat_equation.cpp.

                                           {
  MoFEMFunctionBegin;
 
  // Processes to set output results are integrated in solveSystem()
 
  MoFEMFunctionReturn(0);
}

readMesh()

MoFEMErrorCode HeatEquation::readMesh ( )

private

Examples

heat_equation.cpp.

Definition at line 118 of file heat_equation.cpp.

                                      {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
  CHKERR mField.getInterface(simple);
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
 
  MoFEMFunctionReturn(0);
}

runProgram()

MoFEMErrorCode HeatEquation::runProgram

(

)

Examples

heat_equation.cpp.

Definition at line 424 of file heat_equation.cpp.

                                        {
  MoFEMFunctionBegin;
 
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR setIntegrationRules();
  CHKERR initialCondition();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR outputResults();
 
  MoFEMFunctionReturn(0);
}

setIntegrationRules()

MoFEMErrorCode HeatEquation::setIntegrationRules ( )

private

Examples

heat_equation.cpp.

Definition at line 145 of file heat_equation.cpp.

                                                 {
  MoFEMFunctionBegin;
 
  auto integration_rule = [](int o_row, int o_col, int approx_order) {
    return 2 * approx_order;
  };
 
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
 
  MoFEMFunctionReturn(0);
}

setupProblem()

MoFEMErrorCode HeatEquation::setupProblem ( )

private

Examples

heat_equation.cpp.

Definition at line 129 of file heat_equation.cpp.

                                          {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
 
  int order = 3;
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
  CHKERR simple->setFieldOrder("U", order);
 
  CHKERR simple->setUp();
 
  MoFEMFunctionReturn(0);
}

solveSystem()

MoFEMErrorCode HeatEquation::solveSystem ( )

private

That to work, you have to create solver, as follows,

auto solver = // pipeline_mng->createTSIM( simple->getDM());

That is explicitly use use Simple DM to create solver for DM. Pipeline menage by default creat copy of DM, in case several solvers are used the same DM.

Alternatively you can get dm directly from the solver, i.e.

DM ts_dm;
CHKERR TSGetDM(solver, &ts_dm);
CHKERR DMTSSetIJacobian(
  ts_dm, CalcJacobian::set, getDMTsCtx(ts_dm).get());

Examples

heat_equation.cpp.

Definition at line 300 of file heat_equation.cpp.

                                         {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto create_post_process_element = [&]() {
    auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
    auto det_ptr = boost::make_shared<VectorDouble>();
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    post_proc_fe->getOpPtrVector().push_back(new OpCalculateHOJac<2>(jac_ptr));
    post_proc_fe->getOpPtrVector().push_back(
        new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
    post_proc_fe->getOpPtrVector().push_back(
        new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
 
    auto u_ptr = boost::make_shared<VectorDouble>();
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("U", u_ptr));
 
    using OpPPMap = OpPostProcMapInMoab<2, 2>;
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(post_proc_fe->getPostProcMesh(),
                    post_proc_fe->getMapGaussPts(),
 
                    {{"U", u_ptr}},
 
                    {}, {}, {}
 
                    )
 
    );
 
    return post_proc_fe;
  };
 
  auto set_time_monitor = [&](auto dm, auto solver) {
    MoFEMFunctionBegin;
    boost::shared_ptr<Monitor> monitor_ptr(
        new Monitor(dm, create_post_process_element()));
    boost::shared_ptr<ForcesAndSourcesCore> null;
    CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
                               monitor_ptr, null, null);
    MoFEMFunctionReturn(0);
  };
 
  auto set_fieldsplit_preconditioner = [&](auto solver) {
    MoFEMFunctionBeginHot;
    SNES snes;
    CHKERR TSGetSNES(solver, &snes);
    KSP ksp;
    CHKERR SNESGetKSP(snes, &ksp);
    PC pc;
    CHKERR KSPGetPC(ksp, &pc);
    PetscBool is_pcfs = PETSC_FALSE;
    PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
 
    if (is_pcfs == PETSC_TRUE) {
      auto bc_mng = mField.getInterface<BcManager>();
      auto name_prb = simple->getProblemName();
      auto is_all_bc = bc_mng->getBlockIS(name_prb, "ESSENTIAL", "U", 0, 0);
      int is_all_bc_size;
      CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
      MOFEM_LOG("EXAMPLE", Sev::inform)
          << "Field split block size " << is_all_bc_size;
      CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
                               is_all_bc); // boundary block
    }
    MoFEMFunctionReturnHot(0);
  };
 
  /**
   *  That to work, you have to create solver, as follows,
      \code
      auto solver = // pipeline_mng->createTSIM( simple->getDM());
      \endcode
      That is explicitly use use Simple DM to create solver for DM. Pipeline
      menage by default creat copy of DM, in case several solvers are used the
      same DM.
 
      Alternatively you can get dm directly from the solver, i.e.
      \code
      DM ts_dm;
      CHKERR TSGetDM(solver, &ts_dm);
      CHKERR DMTSSetIJacobian(
        ts_dm, CalcJacobian::set, getDMTsCtx(ts_dm).get());
      \endcode
  */
  auto set_user_ts_jacobian = [&](auto dm) {
    MoFEMFunctionBegin;
    CHKERR DMTSSetIJacobian(dm, CalcJacobian::set, getDMTsCtx(dm).get());
    MoFEMFunctionReturn(0);
  };
 
  auto dm = simple->getDM();
  auto D = createDMVector(dm);
  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
 
  auto solver = pipeline_mng->createTSIM(
      simple->getDM()); // Note DM is set as argument. If DM is not, internal
                        // copy of pipeline DM is created.
  CHKERR set_user_ts_jacobian(dm);
  CHKERR set_time_monitor(dm, solver);
  CHKERR TSSetSolution(solver, D);
  CHKERR TSSetFromOptions(solver);
  CHKERR set_fieldsplit_preconditioner(solver);
  CHKERR TSSetUp(solver);
 
  CHKERR TSSolve(solver, D);
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

boundaryMarker

boost::shared_ptr<std::vector<unsigned char> > HeatEquation::boundaryMarker

private

Examples

heat_equation.cpp.

Definition at line 113 of file heat_equation.cpp.

mField

MoFEM::Interface& HeatEquation::mField

private

Examples

heat_equation.cpp.

Definition at line 110 of file heat_equation.cpp.

---

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

• tutorials/scl-6/heat_equation.cpp