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.
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* License for more details.
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
double dt = 2; // relative to edge length
#include <BasicFiniteElements.hpp>
// Use forms iterators for Grad-Grad term
using OpGradGrad = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
// Use forms for Mass term
using OpMass = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
Simple *simpleInterface;
Range pinchNodes;
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode setIntegrationRules();
MoFEMErrorCode createCommonData();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
MoFEMErrorCode checkResults();
* Use problem specific implementation for second stage of heat methid
struct OpRhs : public DomainEleOp {
OpRhs(boost::shared_ptr<MatrixDouble> u_grad_ptr);
MoFEMErrorCode doWork(int side, EntityType type,
boost::shared_ptr<MatrixDouble> uGradPtr;
//! [Run programme]
CHKERR readMesh();
CHKERR setupProblem();
CHKERR createCommonData();
CHKERR boundaryCondition();
CHKERR assembleSystem();
CHKERR setIntegrationRules();
CHKERR solveSystem();
CHKERR outputResults();
CHKERR checkResults();
//! [Run programme]
//! [Read mesh]
CHKERR mField.getInterface(simpleInterface);
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);
CHKERR simpleInterface->setUp();
//! [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.
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
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]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
// 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 simple = mField.getInterface<Simple>();
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) {
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 = [&]() {
auto simple = mField.getInterface<Simple>();
// 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) {
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
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<3>(jac_ptr, det_ptr, inv_jac_ptr));
new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_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();
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 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
int main(int argc, char *argv[])
static char help[]
constexpr double a
EntitiesFieldData::EntData EntData
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
constexpr double rho
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr double alpha
Catch errors.
Definition: definitions.h:385
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
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:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
FTensor::Index< 'i', SPACE_DIM > i
double dt
Definition: heat_method.cpp:40
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 3 > OpGradGrad
Definition: heat_method.cpp:50
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
boost::shared_ptr< MatrixDouble > uGradPtr
Definition: heat_method.cpp:89
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpRhs(boost::shared_ptr< MatrixDouble > u_grad_ptr)
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode setIntegrationRules()
[Set up problem]
Definition: plot_base.cpp:219
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode outputResults()
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
static UId getUniqueIdCalculate(const DofIdx dof, UId ent_uid)
Data on single entity (This is passed as argument to DataOperator::doWork)
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
MoFEM::FaceElementForcesAndSourcesCore FaceEle
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator