No Matches
SCL-9: Heat method

Table of Contents

Prerequisites of this tutorial include SCL-1: Poisson's equation (homogeneous BC)

Intended learning outcome:
  • general structure of a program developed using MoFEM
  • idea of Simple Interface in MoFEM and how to use it
  • Solving scalar problem embedded in 3d space
  • How to push the developed UDOs to the Pipeline

Source code

* \file heat_method.cpp \example heat_method.cpp
* Calculating geodetic distances using heat method. For details see,
* K. Crane, C. Weischedel, M. Wardetzky, Geodesics in heat: A new approach to
* computing distance based on heat flow, ACM Transactions on Graphics , vol.
* 32, no. 5, pp. 152:1-152:11, 2013.
* Interent resources:
* http://www.cs.cmu.edu/~kmcrane/Projects/HeatMethod/
* http://www.numerical-tours.com/matlab/meshproc_7_geodesic_poisson/
* Solving two problems in sequence. Show hot to use form integrators and how to
* implement user data operator.
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
double dt = 2; // relative to edge length
constexpr int SPACE_DIM = 3;
// Use forms iterators for Grad-Grad term
// Use forms for Mass term
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
* Use problem specific implementation for second stage of heat methid
struct OpRhs : public DomainEleOp {
OpRhs(boost::shared_ptr<MatrixDouble> u_grad_ptr);
boost::shared_ptr<MatrixDouble> uGradPtr;
//! [Run programme]
//! [Run programme]
//! [Read mesh]
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
// FIXME: THis part of algorithm is not working in parallel. If you have bit
// of free time, feel free to generalise code below.
Range nodes;
CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes);
// pinchNodes could be get from BLOCKSET
Range edges;
CHKERR mField.get_moab().get_adjacencies(pinchNodes, 1, false, edges,
double l2;
for (auto e : edges) {
Range nodes;
CHKERR mField.get_moab().get_connectivity(Range(e, e), nodes, false);
MatrixDouble coords(nodes.size(), 3);
CHKERR mField.get_moab().get_coords(nodes, &coords(0, 0));
double l2e = 0;
for (int j = 0; j != 3; ++j) {
l2e += pow(coords(0, j) - coords(1, j), 2);
l2 = std::max(l2, l2e);
dt *= std::sqrt(l2);
//! [Read mesh]
//! [Set up problem]
// Only one field
CHKERR simpleInterface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
constexpr int order = 1;
CHKERR simpleInterface->setFieldOrder("U", order);
//! [Set up problem]
//! [Set integration rule]
// Set integration order. To 2p is enough to integrate mass matrix exactly.
auto rule = [](int, int, int p) -> int { return 2 * p; };
// Set integration rule for elements assembling matrix and vector. Note that
// rule for vector could be 2p-1, but that not change computation time
// significantly. So shorter code is better.
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
//! [Set integration rule]
//! [Create common data]
//! [Create common data]
//! [Boundary condition]
//! [Boundary condition]
//! [Push operators to pipeline]
// This will store gradients at integration points on element
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
// Push element from reference configuration to current configuration in 3d
// space
auto set_domain_general = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
// Operators for assembling matrix for first stage
auto set_domain_lhs_first = [&](auto &pipeline) {
auto one = [](double, double, double) -> double { return 1.; };
pipeline.push_back(new OpGradGrad("U", "U", one));
auto rho = [](double, double, double) -> double { return 1. / dt; };
pipeline.push_back(new OpMass("U", "U", rho));
// Nothing is assembled for vector for first stage of heat method. Later
// simply value of one is set to elements of right hand side vector at pinch
// nodes.
auto set_domain_rhs_first = [&](auto &pipeline) {};
// Operators for assembly of left hand side. This time only Grad-Grad
// operator.
auto set_domain_lhs_second = [&](auto &pipeline) {
auto one = [](double, double, double) { return 1.; };
pipeline.push_back(new OpGradGrad("U", "U", one));
// Now apply on the right hand side vector X = Grad/||Grad||
auto set_domain_rhs_second = [&](auto &pipeline) {
pipeline.push_back(new OpCalculateScalarFieldGradient<3>("U", grad_u_ptr));
pipeline.push_back(new OpRhs(grad_u_ptr));
// Solver first problem
auto solve_first = [&]() {
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto dm = simpleInterface->getDM();
auto D = smartCreateDMVector(dm);
// Note add one at nodes which are pinch nodes
CHKERR VecZeroEntries(F);
auto problem_ptr = mField.get_problem(simple->getProblemName());
auto bit_number = mField.get_field_bit_number("U");
auto dofs_ptr = problem_ptr->getNumeredRowDofsPtr();
for (auto v : pinchNodes) {
const auto uid = DofEntity::getUniqueIdCalculate(
0, FieldEntity::getLocalUniqueIdCalculate(bit_number, v));
auto dof = dofs_ptr->get<Unique_mi_tag>().find(uid);
if (dof != dofs_ptr->end())
VecSetValue(F, (*dof)->getPetscGlobalDofIdx(), 1, INSERT_VALUES);
CHKERR VecAssemblyBegin(F);
CHKERR VecAssemblyEnd(F);
// Solve problem
CHKERR KSPSolve(solver, F, D);
// Solve second problem
auto solve_second = [&]() {
// Note that DOFs on pinch nodes are removed from the problem
auto prb_mng = mField.getInterface<ProblemsManager>();
CHKERR prb_mng->removeDofsOnEntities(simple->getProblemName(), "U",
pinchNodes, 0, 1);
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto dm = simpleInterface->getDM();
auto D = smartCreateDMVector(dm);
CHKERR KSPSolve(solver, F, D);
// Post-process results
auto post_proc = [&](const std::string out_name) {
auto tmp_lhs_fe = pipeline_mng->getDomainLhsFE();
auto tmp_rhs_fe = pipeline_mng->getDomainRhsFE();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto post_proc_fe =
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
new OpCalculateScalarFieldValues("U", u_ptr));
post_proc_fe->getOpPtrVector().push_back(new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"U", u_ptr}},
{{"GRAD_U", grad_ptr}},
{}, {}));
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile(out_name);
tmp_lhs_fe = pipeline_mng->getDomainLhsFE() = tmp_lhs_fe;
tmp_rhs_fe = pipeline_mng->getDomainRhsFE() = tmp_rhs_fe;
CHKERR solve_first();
CHKERR post_proc("out_heat_method_first.h5m");
CHKERR solve_second();
CHKERR post_proc("out_heat_method_second.h5m");
//! [Push operators to pipeline]
//! [Solve]
//! [Solve]
//! [Postprocess results]
//! [Postprocess results]
//! [Check results]
//! [Check results]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Add logging channel for example
auto core_log = logging::core::get();
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
MOFEM_LOG_TAG("EXAMPLE", "example");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [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 insterface
//! [Create MoFEM]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
Example::OpRhs::OpRhs(boost::shared_ptr<MatrixDouble> u_grad_ptr)
: DomainEleOp("U", DomainEleOp::OPROW), uGradPtr(u_grad_ptr) {}
auto nb_dofs = data.getIndices().size();
if (nb_dofs) {
auto t_grad = getFTensor1FromMat<3>(*uGradPtr);
auto nb_base_functions = data.getN().size2();
auto nb_gauss_pts = getGaussPts().size2();
std::array<double, MAX_DOFS_ON_ENTITY> nf;
std::fill(nf.begin(), &nf[nb_dofs], 0);
auto t_diff_base = data.getFTensor1DiffN<3>();
auto t_w = getFTensor0IntegrationWeight();
auto a = getMeasure();
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
double alpha = t_w * a;
const auto l2 = t_grad(i) * t_grad(i);
if (l2 > std::numeric_limits<double>::epsilon())
t_one(i) = t_grad(i) / std::sqrt(l2);
t_one(i) = 0;
size_t bb = 0;
for (; bb != nb_dofs; ++bb) {
nf[bb] -= alpha * t_diff_base(i) * t_one(i);
for (; bb < nb_base_functions; ++bb) {
CHKERR VecSetValues<MoFEM::EssentialBcStorage>(getKSPf(), data, &nf[0],
static Index< 'p', 3 > p
std::string param_file
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
static char help[]
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Catch errors.
Definition: definitions.h:372
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
double dt
Definition: heat_method.cpp:26
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 3 > OpGradGrad
Definition: heat_method.cpp:38
double D
const double v
phase velocity of light in medium (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
double rho
Definition: plastic.cpp:98
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
boost::shared_ptr< MatrixDouble > uGradPtr
Definition: heat_method.cpp:77
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpRhs(boost::shared_ptr< MatrixDouble > u_grad_ptr)
Definition: plastic.cpp:137
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: helmholtz.cpp:106
MoFEMErrorCode assembleSystem()
[Applying essential BC]
Definition: helmholtz.cpp:154
MoFEMErrorCode readMesh()
[run problem]
Definition: helmholtz.cpp:73
MoFEMErrorCode setIntegrationRules()
[Set up problem]
Definition: plot_base.cpp:203
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: contact.cpp:659
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
Definition: helmholtz.cpp:235
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:257
Simple * simpleInterface
Definition: helmholtz.cpp:45
Range pinchNodes
Definition: heat_method.cpp:54
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:172
MoFEM::Interface & mField
Definition: plastic.cpp:144
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:184
MoFEMErrorCode outputResults()
Definition: helmholtz.cpp:255
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.