* \file photon_diffusion.cpp
* \example photon_diffusion.cpp
/* 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 <stdlib.h>
#include <cmath>
#include <BasicFiniteElements.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
template <int DIM> struct ElementsAndOps {};
//! [Define dimension]
constexpr int SPACE_DIM = 3; //< Space dimension of problem, mesh
//! [Define dimension]
using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
using OpDomainGradGrad = FormsIntegrators<DomainEleOp>::Assembly<
using OpDomainTimesScalarField = FormsIntegrators<DomainEleOp>::Assembly<
using OpDomainGradTimesVec = FormsIntegrators<DomainEleOp>::Assembly<
using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
using OpBoundaryMass = FormsIntegrators<BoundaryEleOp>::Assembly<
using OpBoundaryTimeScalarField = FormsIntegrators<BoundaryEleOp>::Assembly<
using OpBoundarySource = FormsIntegrators<BoundaryEleOp>::Assembly<
double n = 1.44; ///< refractive index of diffusive medium
double c = 30.; ///< speed of light (cm/ns)
double v = c / n; ///< phase velocity of light in medium (cm/ns)
double mu_a = 0.09; ///< absorption coefficient (cm^-1)
double mu_sp = 16.5; ///< scattering coefficient (cm^-1)
double flux = 1e3; ///< impulse magnitude
double duration = 0.05; ///< impulse duration (ns)
PetscBool from_initial = PETSC_TRUE;
PetscBool output_volume = PETSC_FALSE;
int order = 3;
double A = 3.0;
double h = 0.5 / A; ///< convective heat coefficient
double D = 1. / (3. * (mu_a + mu_sp));
double inv_v = 1. / v;
* @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 {
Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc,
boost::shared_ptr<PostProcFaceOnRefinedMesh> skin_post_proc)
: dM(dm), postProc(post_proc), skinPostProc(skin_post_proc){};
MoFEMErrorCode preProcess() { return 0; }
MoFEMErrorCode operator()() { return 0; }
MoFEMErrorCode postProcess() {
if (ts_step % save_every_nth_step == 0) {
CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
CHKERR postProc->writeFile(
"out_volume_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
if (skinPostProc) {
CHKERR DMoFEMLoopFiniteElements(dM, "CAMERA_FE", skinPostProc);
CHKERR skinPostProc->writeFile(
"out_camera_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
SmartPetscObj<DM> dM;
boost::shared_ptr<PostProcEle> postProc;
boost::shared_ptr<PostProcFaceOnRefinedMesh> skinPostProc;
// Declaration of the main function to run analysis
MoFEMErrorCode runProgram();
// Declaration of other main functions called in runProgram()
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode setIntegrationRules();
MoFEMErrorCode initialCondition();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
// Main interfaces
// Object to mark boundary entities for the assembling of domain elements
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
boost::shared_ptr<FEMethod> domianLhsFEPtr;
boost::shared_ptr<FEMethod> boundaryLhsFEPtr;
boost::shared_ptr<FEMethod> boundaryRhsFEPtr;
PhotonDiffusion::PhotonDiffusion(MoFEM::Interface &m_field) : mField(m_field) {}
auto *simple = mField.getInterface<Simple>();
CHKERR simple->getOptions();
CHKERR simple->loadFile();
auto *simple = mField.getInterface<Simple>();
CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-flux", &flux, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-duration", &duration,
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-from_initial", &from_initial,
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-output_volume", &output_volume,
MOFEM_LOG("PHOTON", Sev::inform) << "Refractive index: " << n;
MOFEM_LOG("PHOTON", Sev::inform) << "Speed of light (cm/ns): " << c;
MOFEM_LOG("PHOTON", Sev::inform) << "Phase velocity in medium (cm/ns): " << v;
MOFEM_LOG("PHOTON", Sev::inform)
<< "Absorption coefficient (cm^-1): " << mu_a;
MOFEM_LOG("PHOTON", Sev::inform)
<< "Scattering coefficient (cm^-1): " << mu_sp;
MOFEM_LOG("PHOTON", Sev::inform) << "Impulse magnitude: " << flux;
MOFEM_LOG("PHOTON", Sev::inform) << "Impulse duration (ns): " << duration;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_step", &save_every_nth_step,
MOFEM_LOG("PHOTON", Sev::inform) << "Approximation order: " << order;
MOFEM_LOG("PHOTON", Sev::inform) << "Save step: " << save_every_nth_step;
CHKERR simple->setFieldOrder("U", order);
auto set_camera_skin_fe = [&]() {
auto meshset_mng = mField.getInterface<MeshsetsManager>();
Range camera_surface;
const std::string block_name = "CAM";
bool add_fe = false;
if (bit->getName().compare(0, block_name.size(), block_name) == 0) {
MOFEM_LOG("PHOTON", Sev::inform) << "Found CAM block";
CHKERR mField.get_moab().get_entities_by_dimension(
bit->getMeshset(), 2, camera_surface, true);
add_fe = true;
MOFEM_LOG("PHOTON", Sev::noisy) << "CAM block entities:\n"
<< camera_surface;
if (add_fe) {
auto my_simple_set_up = [&]() {
CHKERR simple->defineFiniteElements();
CHKERR simple->defineProblem(PETSC_TRUE);
CHKERR simple->buildFields();
CHKERR simple->buildFiniteElements();
if (mField.check_finite_element("CAMERA_FE")) {
CHKERR DMMoFEMAddElement(simple->getDM(), "CAMERA_FE");
CHKERR simple->buildProblem();
CHKERR set_camera_skin_fe();
CHKERR my_simple_set_up();
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);
auto bc_mng = mField.getInterface<BcManager>();
auto *simple = mField.getInterface<Simple>();
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "MIX", "U", 0,
0, false);
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "SPOT", "U",
0, 0, false);
auto bc_mng = mField.getInterface<BcManager>();
auto &bc_map = bc_mng->getBcMapByBlockName();
auto add_domain_base_ops = [&](auto &pipeline) {
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
pipeline.push_back(new OpCalculateHOJacVolume(jac_ptr));
pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetHOInvJacToScalarBases(H1, inv_jac_ptr));
pipeline.push_back(new OpSetHOWeights(det_ptr));
auto add_domain_lhs_ops = [&](auto &pipeline) {
pipeline.push_back(new OpDomainGradGrad(
"U", "U", [](double, double, double) -> double { return D; }));
auto get_mass_coefficient = [&](const double, const double, const double) {
return inv_v * domianLhsFEPtr->ts_a + mu_a;
pipeline.push_back(new OpDomainMass("U", "U", get_mass_coefficient));
auto add_domain_rhs_ops = [&](auto &pipeline) {
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
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 OpCalculateScalarFieldValues("U", u_at_gauss_pts));
new OpCalculateScalarFieldValuesDot("U", dot_u_at_gauss_pts));
pipeline.push_back(new OpDomainGradTimesVec(
"U", grad_u_at_gauss_pts,
[](double, double, double) -> double { return D; }));
pipeline.push_back(new OpDomainTimesScalarField(
"U", dot_u_at_gauss_pts,
[](const double, const double, const double) { return inv_v; }));
pipeline.push_back(new OpDomainTimesScalarField(
"U", u_at_gauss_pts,
[](const double, const double, const double) { return mu_a; }));
auto source_term = [&](const double, const double, const double) {
return 0;
pipeline.push_back(new OpDomainSource("U", source_term));
auto add_boundary_base_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetHOWeightsOnFace());
auto add_lhs_base_ops = [&](auto &pipeline) {
for (auto b : bc_map) {
if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
pipeline.push_back(new OpBoundaryMass(
"U", "U",
[](const double, const double, const double) { return h; },
auto add_rhs_base_ops = [&](auto &pipeline) {
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
pipeline.push_back(new OpCalculateScalarFieldValues("U", u_at_gauss_pts));
for (auto b : bc_map) {
if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
pipeline.push_back(new OpBoundaryTimeScalarField(
"U", u_at_gauss_pts,
[](const double, const double, const double) { return h; },
auto pipeline_mng = mField.getInterface<PipelineManager>();
domianLhsFEPtr = pipeline_mng->getDomainLhsFE();
boundaryLhsFEPtr = pipeline_mng->getBoundaryLhsFE();
boundaryRhsFEPtr = pipeline_mng->getBoundaryRhsFE();
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto create_post_process_element = [&]() {
auto post_froc_fe = boost::make_shared<PostProcEle>(mField);
return post_froc_fe;
auto create_post_process_camera_element = [&]() {
if (mField.check_finite_element("CAMERA_FE")) {
auto post_proc_skin =
CHKERR post_proc_skin->addFieldValuesPostProc("U");
CHKERR post_proc_skin->addFieldValuesGradientPostProcOnSkin(
"U", simple->getDomainFEName());
return post_proc_skin;
} else {
return boost::shared_ptr<PostProcFaceOnRefinedMesh>();
auto set_time_monitor = [&](auto dm, auto solver) {
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);
auto dm = simple->getDM();
auto D = smartCreateDMVector(dm);
if (from_initial) {
MOFEM_LOG("PHOTON", Sev::inform)
<< "reading vector in binary from vector.dat ...";
PetscViewer viewer;
PetscViewerBinaryOpen(PETSC_COMM_WORLD, "initial_vector.dat",
FILE_MODE_READ, &viewer);
VecLoad(D, viewer);
auto solver = pipeline_mng->createTS();
CHKERR TSSetSolution(solver, D);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetSolution(solver, D);
CHKERR TSSetFromOptions(solver);
CHKERR TSSetUp(solver);
CHKERR TSSolve(solver, NULL);
// Processes to set output results are integrated in solveSystem()
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(), "PHOTON"));
MOFEM_LOG_TAG("PHOTON", "photon_diffusion")
// Error handling
try {
// Register MoFEM discrete manager in PETSc
DMType dm_name = "DMMOFEM";
// Create MOAB instance
moab::Core mb_instance; // mesh database
moab::Interface &moab = mb_instance; // mesh database interface
// Create MoFEM instance
MoFEM::Core core(moab); // finite element database
MoFEM::Interface &m_field = core; // finite element interface
// Run the main analysis
PhotonDiffusion heat_problem(m_field);
CHKERR heat_problem.runProgram();
// Finish work: cleaning memory, getting statistics, etc.
return 0;
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
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
Definition: definitions.h:161
#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 integration_rule
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:461
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
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:544
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
#define MOFEM_LOG(channel, severity)
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
auto bit
set bit
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 2 > OpDomainGradGrad
Definition: helmholtz.cpp:37
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition: helmholtz.cpp:43
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpDomainTimesScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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: DMMMoFEM.cpp:994
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
int main(int argc, char *argv[])
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBoundaryTimeScalarField
double mu_sp
scattering coefficient (cm^-1)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
double A
static char help[]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpDomainTimesScalarField
double c
speed of light (cm/ns)
double v
phase velocity of light in medium (cm/ns)
double flux
impulse magnitude
double inv_v
constexpr int SPACE_DIM
[Define dimension]
double h
convective heat coefficient
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
double duration
impulse duration (ns)
int save_every_nth_step
double n
refractive index of diffusive medium
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
int order
PetscBool from_initial
double D
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
PetscBool output_volume
double mu_a
absorption coefficient (cm^-1)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:93
static constexpr int approx_order
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual moab::Interface & get_moab()=0
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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]
MoFEMErrorCode assembleSystem()
boost::shared_ptr< FEMethod > boundaryRhsFEPtr
MoFEMErrorCode solveSystem()
MoFEM::Interface & mField
MoFEMErrorCode readMesh()
MoFEMErrorCode outputResults()
MoFEMErrorCode initialCondition()
boost::shared_ptr< FEMethod > domianLhsFEPtr
PhotonDiffusion(MoFEM::Interface &m_field)
MoFEMErrorCode runProgram()
boost::shared_ptr< FEMethod > boundaryLhsFEPtr
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode setupProblem()
Post processing.