v0.14.0
Loading...
Searching...
No Matches
photon_diffusion.cpp
/**
* \file photon_diffusion.cpp
* \example photon_diffusion.cpp
*
**/
#include <stdlib.h>
#include <cmath>
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 DomainEle = VolumeElementForcesAndSourcesCore;
using VolSideFe = VolumeElementForcesAndSourcesCoreOnSide;
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
private:
// Declaration of other main functions called in runProgram()
// 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 {
};
};
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, 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) {}
DataForcesAndSourcesCore::EntData &data) {
if (type != MBVERTEX)
}
};
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){};
}
}
CHKERR VecZeroEntries(commonDataPtr->petscVec);
CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
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];
CHKERR postProc->writeFile("out_volume_" +
boost::lexical_cast<std::string>(ts_step) +
".h5m");
}
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->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()
}
}
EntData &data) {
const int nb_integration_pts = getGaussPts().size2();
const double area = getMeasure();
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";
// 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;
}
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()
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
constexpr int SPACE_DIM
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ 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 ...
@ BLOCKSET
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
FTensor::Index< 'n', SPACE_DIM > n
auto integration_rule
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:483
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:509
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
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:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1003
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
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#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.
auto bit
set bit
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition helmholtz.cpp:30
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition helmholtz.cpp:24
double mu_sp
scattering coefficient (cm^-1)
int numHoLevels
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
const double c
speed of light (cm/ns)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
double D
PetscBool output_volume
double mu_a
absorption coefficient (cm^-1)
const double v
phase velocity of light in medium (cm/ns)
VolumeElementForcesAndSourcesCoreOnSide VolSideFe
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
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:1042
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)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
auto createGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr AssemblyType A
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
double mu_sp
scattering coefficient (cm^-1)
double A
static char help[]
const double inv_v
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
PetscBool output_camera
double h
const double c
speed of light (cm/ns)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
char init_data_file_name[255]
int save_every_nth_step
int order
PetscBool from_initial
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
double D
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
PetscBool output_volume
double mu_a
absorption coefficient (cm^-1)
const double v
phase velocity of light in medium (cm/ns)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
const double n
refractive index of diffusive medium
static constexpr int approx_order
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition seepage.cpp:70
Simple interface for fast problem set-up.
Definition BcManager.hpp:25
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =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
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode loopSideVolumes(const string fe_name, VolumeElementForcesAndSourcesCoreOnSide &fe_method)
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Interface for managing meshsets containing materials and boundary conditions.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
[Push operators to pipeline]
boost::shared_ptr< VectorDouble > approxVals
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< PostProcFaceEle > skinPostProc
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)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode operator()()
function is run for every finite element
boost::shared_ptr< BoundaryEle > skinPostProcInteg
boost::shared_ptr< PostProcEle > postProc
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
OpCameraInteg(boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetScalarFieldGradientValuesOnSkin(boost::shared_ptr< VolSideFe > side_fe)
MoFEMErrorCode assembleSystem()
boost::shared_ptr< FEMethod > boundaryRhsFEPtr
MoFEMErrorCode solveSystem()
MoFEM::Interface & mField
MoFEMErrorCode readMesh()
MoFEMErrorCode outputResults()
MoFEMErrorCode initialCondition()
PhotonDiffusion(MoFEM::Interface &m_field)
boost::shared_ptr< FEMethod > domainLhsFEPtr
MoFEMErrorCode runProgram()
MoFEMErrorCode createCommonData()
boost::shared_ptr< FEMethod > boundaryLhsFEPtr
MoFEMErrorCode boundaryCondition()
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode setupProblem()
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore > PostProcFaceEle