v0.14.0
photon_diffusion.cpp
/**
* \file photon_diffusion.cpp
* \example photon_diffusion.cpp
*
**/
#include <stdlib.h>
#include <cmath>
#include <MoFEM.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]
PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
const double n = 1.44; ///< refractive index of diffusive medium
const double c = 30.; ///< speed of light (cm/ns)
const double v = c / n; ///< phase velocity of light in medium (cm/ns)
const double inv_v = 1. / v;
double mu_a; ///< absorption coefficient (cm^-1)
double mu_sp; ///< scattering coefficient (cm^-1)
double D;
double A;
double h;
PetscBool from_initial = PETSC_TRUE;
PetscBool output_volume = PETSC_FALSE;
PetscBool output_camera = PETSC_FALSE;
int order = 2;
char init_data_file_name[255] = "init_file.dat";
int numHoLevels = 1;
public:
// Declaration of the main function to run analysis
MoFEMErrorCode runProgram();
private:
// Declaration of other main functions called in runProgram()
MoFEMErrorCode readMesh();
MoFEMErrorCode createCommonData();
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> domainLhsFEPtr;
boost::shared_ptr<FEMethod> boundaryLhsFEPtr;
boost::shared_ptr<FEMethod> boundaryRhsFEPtr;
struct CommonData {
boost::shared_ptr<VectorDouble> approxVals;
enum VecElements {
VALUES_INTEG = 0,
LAST_ELEMENT
};
};
boost::shared_ptr<CommonData> commonDataPtr;
struct OpCameraInteg : public BoundaryEleOp {
boost::shared_ptr<CommonData> commonDataPtr;
OpCameraInteg(boost::shared_ptr<CommonData> common_data_ptr)
: BoundaryEleOp("PHOTON_FLUENCE_RATE", OPROW),
commonDataPtr(common_data_ptr) {
std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
doEntities[MBTRI] = doEntities[MBQUAD] = true;
}
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data);
};
struct OpGetScalarFieldGradientValuesOnSkin : public BoundaryEleOp {
boost::shared_ptr<VolSideFe> sideOpFe;
OpGetScalarFieldGradientValuesOnSkin(boost::shared_ptr<VolSideFe> side_fe)
: BoundaryEleOp("PHOTON_FLUENCE_RATE", OPROW), sideOpFe(side_fe) {}
MoFEMErrorCode doWork(int side, EntityType type,
if (type != MBVERTEX)
CHKERR loopSideVolumes("dFE", *sideOpFe);
}
};
struct Monitor : public FEMethod {
Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc,
boost::shared_ptr<PostProcFaceEle> skin_post_proc,
boost::shared_ptr<BoundaryEle> skin_post_proc_integ,
boost::shared_ptr<CommonData> common_data_ptr)
: dM(dm), postProc(post_proc), skinPostProc(skin_post_proc),
skinPostProcInteg(skin_post_proc_integ),
commonDataPtr(common_data_ptr){};
MoFEMErrorCode preProcess() {
}
MoFEMErrorCode operator()() {
}
MoFEMErrorCode postProcess() {
CHKERR VecZeroEntries(commonDataPtr->petscVec);
CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR DMoFEMLoopFiniteElements(dM, "CAMERA_FE", skinPostProcInteg);
CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, ADD_VALUES,
SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, ADD_VALUES,
SCATTER_REVERSE);
CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
const double *array;
CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
MOFEM_LOG("PHOTON", Sev::inform) << "Fluence rate integral: " << array[0];
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 (output_camera && skinPostProc) {
CHKERR DMoFEMLoopFiniteElements(dM, "CAMERA_FE", skinPostProc);
CHKERR skinPostProc->writeFile(
"out_camera_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
}
}
}
private:
boost::shared_ptr<PostProcEle> postProc;
boost::shared_ptr<PostProcFaceEle> skinPostProc;
boost::shared_ptr<BoundaryEle> skinPostProcInteg;
boost::shared_ptr<CommonData> commonDataPtr;
};
};
PhotonDiffusion::PhotonDiffusion(MoFEM::Interface &m_field) : mField(m_field) {}
CHKERR simple->getOptions();
CHKERR simple->loadFile();
}
commonDataPtr = boost::make_shared<CommonData>();
PetscInt ghosts[1] = {0};
commonDataPtr->petscVec =
createGhostVector(mField.get_comm(), 1, 1, 0, ghosts);
else
commonDataPtr->petscVec =
createGhostVector(mField.get_comm(), 0, 1, 1, ghosts);
commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
}
CHKERR simple->addDomainField("PHOTON_FLUENCE_RATE", H1,
CHKERR simple->addBoundaryField("PHOTON_FLUENCE_RATE", H1,
CHKERR PetscOptionsGetString(PETSC_NULL, "", "-initial_file",
init_data_file_name, 255, PETSC_NULL);
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-from_initial", &from_initial,
PETSC_NULL);
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-output_volume", &output_volume,
PETSC_NULL);
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-output_camera", &output_camera,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-mu_a", &mu_a, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-mu_sp", &mu_sp, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-coef_A", &A, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_step", &save_every_nth_step,
PETSC_NULL);
h = 0.5 / A;
D = 1. / (3. * (mu_a + mu_sp));
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) << "Inverse velocity : " << inv_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) << "Diffusion coefficient D : " << D;
MOFEM_LOG("PHOTON", Sev::inform) << "Coefficient A : " << A;
MOFEM_LOG("PHOTON", Sev::inform) << "Coefficient h : " << h;
MOFEM_LOG("PHOTON", Sev::inform) << "Approximation order: " << order;
MOFEM_LOG("PHOTON", Sev::inform) << "Save step: " << save_every_nth_step;
CHKERR simple->setFieldOrder("PHOTON_FLUENCE_RATE", 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) {
"PHOTON_FLUENCE_RATE");
"CAMERA_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>();
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "EXT",
"PHOTON_FLUENCE_RATE", 0, 0, false);
// Get boundary edges marked in block named "BOUNDARY_CONDITION"
Range boundary_ents;
std::string entity_name = it->getName();
if (entity_name.compare(0, 3, "INT") == 0) {
CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
boundary_ents, true);
}
}
// Add vertices to boundary entities
Range boundary_verts;
CHKERR mField.get_moab().get_connectivity(boundary_ents, boundary_verts,
true);
boundary_ents.merge(boundary_verts);
// Remove DOFs as homogeneous boundary condition is used
CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
simple->getProblemName(), "PHOTON_FLUENCE_RATE", boundary_ents);
}
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 OpCalculateHOJac<3>(jac_ptr));
pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(H1, inv_jac_ptr));
pipeline.push_back(new OpSetHOWeights(det_ptr));
};
auto add_domain_lhs_ops = [&](auto &pipeline) {
pipeline.push_back(new OpDomainGradGrad(
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
[](double, double, double) -> double { return D; }));
auto get_mass_coefficient = [&](const double, const double, const double) {
return inv_v * domainLhsFEPtr->ts_a + mu_a;
};
pipeline.push_back(new OpDomainMass(
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE", 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>();
"PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts));
pipeline.push_back(new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
u_at_gauss_pts));
pipeline.push_back(new OpCalculateScalarFieldValuesDot(
"PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts));
pipeline.push_back(new OpDomainGradTimesVec(
"PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts,
[](double, double, double) -> double { return D; }));
pipeline.push_back(new OpDomainTimesScalarField(
"PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts,
[](const double, const double, const double) { return inv_v; }));
pipeline.push_back(new OpDomainTimesScalarField(
"PHOTON_FLUENCE_RATE", u_at_gauss_pts,
[](const double, const double, const double) { return mu_a; }));
};
auto add_boundary_base_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetHOWeightsOnFace());
};
auto add_boundary_lhs_ops = [&](auto &pipeline) {
for (auto b : bc_map) {
if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
pipeline.push_back(new OpBoundaryMass(
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
[](const double, const double, const double) { return h; },
b.second->getBcEntsPtr()));
}
}
};
auto add_boundary_rhs_ops = [&](auto &pipeline) {
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
pipeline.push_back(new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
u_at_gauss_pts));
for (auto b : bc_map) {
if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
pipeline.push_back(new OpBoundaryTimeScalarField(
"PHOTON_FLUENCE_RATE", u_at_gauss_pts,
[](const double, const double, const double) { return h; },
b.second->getBcEntsPtr()));
}
}
};
auto pipeline_mng = mField.getInterface<PipelineManager>();
add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
domainLhsFEPtr = pipeline_mng->getDomainLhsFE();
boundaryLhsFEPtr = pipeline_mng->getBoundaryLhsFE();
boundaryRhsFEPtr = pipeline_mng->getBoundaryRhsFE();
}
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto create_post_process_element = [&]() {
auto post_froc_fe = boost::make_shared<PostProcEle>(mField);
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
post_froc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE", u_ptr));
post_froc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>("PHOTON_FLUENCE_RATE",
grad_ptr));
post_froc_fe->getOpPtrVector().push_back(new OpPPMap(
post_froc_fe->getPostProcMesh(), post_froc_fe->getMapGaussPts(),
{{"PHOTON_FLUENCE_RATE", u_ptr}},
{{"GRAD_PHOTON_FLUENCE_RATE", grad_ptr}}, {}, {}));
return post_froc_fe;
};
auto create_post_process_camera_element = [&]() {
if (mField.check_finite_element("CAMERA_FE")) {
auto post_proc_skin = boost::make_shared<PostProcFaceEle>(mField);
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
mField, simple->getDomainFEName(), SPACE_DIM);
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
// push operators to side element
op_loop_side->getOpPtrVector().push_back(
op_loop_side->getOpPtrVector().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
op_loop_side->getOpPtrVector().push_back(
op_loop_side->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE", u_ptr));
op_loop_side->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>("PHOTON_FLUENCE_RATE", grad_ptr));
// push op to boundary element
post_proc_skin->getOpPtrVector().push_back(op_loop_side);
post_proc_skin->getOpPtrVector().push_back(new OpPPMap(
post_proc_skin->getPostProcMesh(), post_proc_skin->getMapGaussPts(),
{{"PHOTON_FLUENCE_RATE", u_ptr}},
{{"GRAD_PHOTON_FLUENCE_RATE", grad_ptr}}, {}, {}));
return post_proc_skin;
} else {
return boost::shared_ptr<PostProcFaceEle>();
}
};
auto create_post_process_integ_camera_element = [&]() {
if (mField.check_finite_element("CAMERA_FE")) {
auto post_proc_integ_skin = boost::make_shared<BoundaryEle>(mField);
post_proc_integ_skin->getOpPtrVector().push_back(
post_proc_integ_skin->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
commonDataPtr->approxVals));
post_proc_integ_skin->getOpPtrVector().push_back(
new OpCameraInteg(commonDataPtr));
return post_proc_integ_skin;
} else {
return boost::shared_ptr<BoundaryEle>();
}
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(new Monitor(
dm, create_post_process_element(), create_post_process_camera_element(),
create_post_process_integ_camera_element(), commonDataPtr));
boost::shared_ptr<ForcesAndSourcesCore> null;
CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
monitor_ptr, null, null);
};
auto dm = simple->getDM();
auto X = createDMVector(dm);
if (from_initial) {
MOFEM_LOG("PHOTON", Sev::inform) << "reading vector in binary from file "
<< init_data_file_name << " ...";
PetscViewer viewer;
PetscViewerBinaryOpen(PETSC_COMM_WORLD, init_data_file_name, FILE_MODE_READ,
&viewer);
VecLoad(X, viewer);
CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_REVERSE);
}
auto solver = pipeline_mng->createTS();
CHKERR TSSetSolution(solver, X);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetSolution(solver, X);
CHKERR TSSetFromOptions(solver);
CHKERR TSSetUp(solver);
CHKERR TSSolve(solver, NULL);
CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_REVERSE);
}
// Processes to set output results are integrated in solveSystem()
}
}
const int nb_integration_pts = getGaussPts().size2();
const double area = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
double values_integ = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * area;
values_integ += alpha * t_val;
++t_w;
++t_val;
}
constexpr std::array<int, 1> indices = {CommonData::VALUES_INTEG};
std::array<double, 1> values;
values[0] = values_integ;
CHKERR VecSetValues(commonDataPtr->petscVec, 1, indices.data(), values.data(),
ADD_VALUES);
}
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "PHOTON"));
LogManager::setLog("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;
}
n
const double n
refractive index of diffusive medium
Definition: photon_diffusion.cpp:51
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
PhotonDiffusion::boundaryRhsFEPtr
boost::shared_ptr< FEMethod > boundaryRhsFEPtr
Definition: photon_diffusion.cpp:100
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
main
int main(int argc, char *argv[])
Definition: photon_diffusion.cpp:658
help
static char help[]
Definition: photon_diffusion.cpp:13
OpDomainGradTimesVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
Definition: initial_diffusion.cpp:34
order
int order
Definition: photon_diffusion.cpp:66
A
double A
Definition: photon_diffusion.cpp:59
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
PhotonDiffusion::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: photon_diffusion.cpp:112
PhotonDiffusion
Definition: initial_diffusion.cpp:74
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
PhotonDiffusion::boundaryLhsFEPtr
boost::shared_ptr< FEMethod > boundaryLhsFEPtr
Definition: photon_diffusion.cpp:99
OpDomainGradGrad
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition: photon_diffusion.cpp:36
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
OpBoundaryTimeScalarField
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
Definition: photon_diffusion.cpp:47
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
OpDomainGradTimesVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
Definition: photon_diffusion.cpp:40
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
PhotonDiffusion::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
Definition: initial_diffusion.cpp:236
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
PhotonDiffusion::PhotonDiffusion
PhotonDiffusion(MoFEM::Interface &m_field)
Definition: initial_diffusion.cpp:138
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
Definition: child_and_parent.cpp:53
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
PhotonDiffusion::readMesh
MoFEMErrorCode readMesh()
Definition: initial_diffusion.cpp:140
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
numHoLevels
int numHoLevels
Definition: photon_diffusion.cpp:70
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
PhotonDiffusion::domainLhsFEPtr
boost::shared_ptr< FEMethod > domainLhsFEPtr
Definition: photon_diffusion.cpp:98
OpBoundaryMass
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: seepage.cpp:71
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
MoFEM::OpCalculateHOJac
Definition: HODataOperators.hpp:267
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
PhotonDiffusion::setupProblem
MoFEMErrorCode setupProblem()
Definition: initial_diffusion.cpp:151
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
VolSideFe
VolumeElementForcesAndSourcesCoreOnSide VolSideFe
Definition: photon_diffusion.cpp:31
output_volume
PetscBool output_volume
Definition: photon_diffusion.cpp:63
save_every_nth_step
int save_every_nth_step
Definition: photon_diffusion.cpp:67
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::CoreInterface::add_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
MoFEM::createGhostVector
auto createGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
Definition: PetscSmartObj.hpp:179
from_initial
PetscBool from_initial
Definition: photon_diffusion.cpp:62
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
mu_sp
double mu_sp
scattering coefficient (cm^-1)
Definition: photon_diffusion.cpp:57
OpDomainTimesScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
Definition: initial_diffusion.cpp:32
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
BoundaryEleOp
mu_a
double mu_a
absorption coefficient (cm^-1)
Definition: photon_diffusion.cpp:56
PhotonDiffusion::mField
MoFEM::Interface & mField
Definition: initial_diffusion.cpp:135
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
PhotonDiffusion::initialCondition
MoFEMErrorCode initialCondition()
Definition: initial_diffusion.cpp:252
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::OpSetHOInvJacToScalarBases
Set inverse jacobian to base functions.
Definition: HODataOperators.hpp:73
double
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
OpDomainTimesScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
Definition: photon_diffusion.cpp:38
convert.type
type
Definition: convert.py:64
MoFEM::BcManager::getBcMapByBlockName
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:243
OpBoundaryMass
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
Definition: photon_diffusion.cpp:45
OpBoundarySource
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition: helmholtz.cpp:31
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
h
double h
Definition: photon_diffusion.cpp:60
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
PhotonDiffusion::OpCameraInteg::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: photon_diffusion.cpp:631
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
D
double D
Definition: photon_diffusion.cpp:58
MoFEM::PostProcBrokenMeshInMoab< VolumeElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:670
PhotonDiffusion::OpCameraInteg::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: photon_diffusion.cpp:115
MoFEM::OpSetHOWeightsOnFace
Modify integration weights on face to take in account higher-order geometry.
Definition: HODataOperators.hpp:122
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::OpCalculateHOJac< 3 >
Definition: HODataOperators.hpp:269
PhotonDiffusion::runProgram
MoFEMErrorCode runProgram()
Definition: initial_diffusion.cpp:383
PhotonDiffusion::createCommonData
MoFEMErrorCode createCommonData()
Definition: photon_diffusion.cpp:223
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
PhotonDiffusion::CommonData::VALUES_INTEG
@ VALUES_INTEG
Definition: photon_diffusion.cpp:107
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
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
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
PhotonDiffusion::solveSystem
MoFEMErrorCode solveSystem()
Definition: initial_diffusion.cpp:318
BiLinearForm
SPACE_DIM
constexpr int SPACE_DIM
[Define dimension]
Definition: photon_diffusion.cpp:18
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
PhotonDiffusion::boundaryCondition
MoFEMErrorCode boundaryCondition()
Definition: initial_diffusion.cpp:257
c
const double c
speed of light (cm/ns)
Definition: photon_diffusion.cpp:52
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
v
const double v
phase velocity of light in medium (cm/ns)
Definition: photon_diffusion.cpp:53
inv_v
const double inv_v
Definition: photon_diffusion.cpp:54
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
ElementsAndOps
Definition: child_and_parent.cpp:18
Range
output_camera
PetscBool output_camera
Definition: photon_diffusion.cpp:64
DomainEleOp
PhotonDiffusion::assembleSystem
MoFEMErrorCode assembleSystem()
Definition: initial_diffusion.cpp:284
MoFEM::CoreTmp< 0 >::Initialize
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
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
OpDomainGradGrad
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition: helmholtz.cpp:25
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
Monitor
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:783
approx_order
int approx_order
Definition: test_broken_space.cpp:50
CommonData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22
MoFEM::OpSetHOWeights
Set inverse jacobian to base functions.
Definition: HODataOperators.hpp:159
MoFEM::DMMoFEMTSSetMonitor
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: DMMoFEM.cpp:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3249
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
Definition: photon_diffusion.cpp:34
init_data_file_name
char init_data_file_name[255]
Definition: photon_diffusion.cpp:69
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::CoreInterface::check_finite_element
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
MoFEM::VolumeElementForcesAndSourcesCoreOnSide
Base volume element used to integrate on skeleton.
Definition: VolumeElementForcesAndSourcesCoreOnSide.hpp:22
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::SmartPetscObj< Vec >
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
Definition: child_and_parent.cpp:55
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1290
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:677
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
PhotonDiffusion::outputResults
MoFEMErrorCode outputResults()
Definition: initial_diffusion.cpp:349