13static char help[] =
"...\n\n";
53const double v =
c /
n;
133 DataForcesAndSourcesCore::EntData &data) {
135 if (type != MBVERTEX)
145 boost::shared_ptr<PostProcFaceEle> skin_post_proc,
146 boost::shared_ptr<BoundaryEle> skin_post_proc_integ,
147 boost::shared_ptr<CommonData> common_data_ptr)
182 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Fluence rate integral: " << array[0];
188 boost::lexical_cast<std::string>(
ts_step) +
194 "out_camera_" + boost::lexical_cast<std::string>(
ts_step) +
226 PetscInt ghosts[1] = {0};
233 commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
266 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Refractive index: " <<
n;
267 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Speed of light (cm/ns): " <<
c;
268 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Phase velocity in medium (cm/ns): " <<
v;
269 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Inverse velocity : " <<
inv_v;
271 <<
"Absorption coefficient (cm^-1): " <<
mu_a;
273 <<
"Scattering coefficient (cm^-1): " <<
mu_sp;
274 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Diffusion coefficient D : " <<
D;
275 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Coefficient A : " <<
A;
276 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Coefficient h : " <<
h;
278 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Approximation order: " <<
order;
283 auto set_camera_skin_fe = [&]() {
287 Range camera_surface;
288 const std::string block_name =
"CAM";
292 if (
bit->getName().compare(0, block_name.size(), block_name) == 0) {
293 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Found CAM block";
295 bit->getMeshset(), 2, camera_surface,
true);
300 MOFEM_LOG(
"PHOTON", Sev::noisy) <<
"CAM block entities:\n"
306 "PHOTON_FLUENCE_RATE");
313 auto my_simple_set_up = [&]() {
329 CHKERR set_camera_skin_fe();
330 CHKERR my_simple_set_up();
361 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"EXT",
362 "PHOTON_FLUENCE_RATE", 0, 0,
false);
367 std::string entity_name = it->getName();
368 if (entity_name.compare(0, 3,
"INT") == 0) {
370 boundary_ents,
true);
374 Range boundary_verts;
377 boundary_ents.merge(boundary_verts);
381 simple->getProblemName(),
"PHOTON_FLUENCE_RATE", boundary_ents);
392 auto add_domain_base_ops = [&](
auto &pipeline) {
393 auto jac_ptr = boost::make_shared<MatrixDouble>();
394 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
395 auto det_ptr = boost::make_shared<VectorDouble>();
402 auto add_domain_lhs_ops = [&](
auto &pipeline) {
404 "PHOTON_FLUENCE_RATE",
"PHOTON_FLUENCE_RATE",
405 [](
double,
double,
double) ->
double {
return D; }));
410 "PHOTON_FLUENCE_RATE",
"PHOTON_FLUENCE_RATE", get_mass_coefficient));
413 auto add_domain_rhs_ops = [&](
auto &pipeline) {
414 auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
415 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
416 auto dot_u_at_gauss_pts = boost::make_shared<VectorDouble>();
418 "PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts));
422 "PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts));
424 "PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts,
425 [](
double,
double,
double) ->
double {
return D; }));
427 "PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts,
428 [](
const double,
const double,
const double) {
return inv_v; }));
430 "PHOTON_FLUENCE_RATE", u_at_gauss_pts,
431 [](
const double,
const double,
const double) {
return mu_a; }));
434 auto add_boundary_base_ops = [&](
auto &pipeline) {
438 auto add_boundary_lhs_ops = [&](
auto &pipeline) {
439 for (
auto b : bc_map) {
440 if (std::regex_match(
b.first, std::regex(
"(.*)EXT(.*)"))) {
442 "PHOTON_FLUENCE_RATE",
"PHOTON_FLUENCE_RATE",
444 [](
const double,
const double,
const double) {
return h; },
446 b.second->getBcEntsPtr()));
451 auto add_boundary_rhs_ops = [&](
auto &pipeline) {
452 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
455 for (
auto b : bc_map) {
456 if (std::regex_match(
b.first, std::regex(
"(.*)EXT(.*)"))) {
458 "PHOTON_FLUENCE_RATE", u_at_gauss_pts,
460 [](
const double,
const double,
const double) {
return h; },
462 b.second->getBcEntsPtr()));
468 add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
469 add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
470 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
471 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
473 add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
474 add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
475 add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
476 add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
491 auto create_post_process_element = [&]() {
492 auto post_froc_fe = boost::make_shared<PostProcEle>(
mField);
493 auto u_ptr = boost::make_shared<VectorDouble>();
494 auto grad_ptr = boost::make_shared<MatrixDouble>();
495 post_froc_fe->getOpPtrVector().push_back(
497 post_froc_fe->getOpPtrVector().push_back(
500 post_froc_fe->getOpPtrVector().push_back(
new OpPPMap(
501 post_froc_fe->getPostProcMesh(), post_froc_fe->getMapGaussPts(),
502 {{
"PHOTON_FLUENCE_RATE", u_ptr}},
503 {{
"GRAD_PHOTON_FLUENCE_RATE", grad_ptr}}, {}, {}));
507 auto create_post_process_camera_element = [&]() {
510 auto post_proc_skin = boost::make_shared<PostProcFaceEle>(
mField);
512 auto u_ptr = boost::make_shared<VectorDouble>();
513 auto grad_ptr = boost::make_shared<MatrixDouble>();
518 auto jac_ptr = boost::make_shared<MatrixDouble>();
519 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
520 auto det_ptr = boost::make_shared<VectorDouble>();
523 op_loop_side->getOpPtrVector().push_back(
525 op_loop_side->getOpPtrVector().push_back(
527 op_loop_side->getOpPtrVector().push_back(
529 op_loop_side->getOpPtrVector().push_back(
531 op_loop_side->getOpPtrVector().push_back(
534 post_proc_skin->getOpPtrVector().push_back(op_loop_side);
536 post_proc_skin->getOpPtrVector().push_back(
new OpPPMap(
537 post_proc_skin->getPostProcMesh(), post_proc_skin->getMapGaussPts(),
538 {{
"PHOTON_FLUENCE_RATE", u_ptr}},
539 {{
"GRAD_PHOTON_FLUENCE_RATE", grad_ptr}}, {}, {}));
541 return post_proc_skin;
543 return boost::shared_ptr<PostProcFaceEle>();
547 auto create_post_process_integ_camera_element = [&]() {
550 auto post_proc_integ_skin = boost::make_shared<BoundaryEle>(
mField);
551 post_proc_integ_skin->getOpPtrVector().push_back(
553 post_proc_integ_skin->getOpPtrVector().push_back(
556 post_proc_integ_skin->getOpPtrVector().push_back(
559 return post_proc_integ_skin;
561 return boost::shared_ptr<BoundaryEle>();
565 auto set_time_monitor = [&](
auto dm,
auto solver) {
567 boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(
568 dm, create_post_process_element(), create_post_process_camera_element(),
570 boost::shared_ptr<ForcesAndSourcesCore> null;
572 monitor_ptr, null, null);
576 auto dm =
simple->getDM();
581 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"reading vector in binary from file "
591 auto solver = pipeline_mng->createTS();
593 CHKERR TSSetSolution(solver, X);
594 CHKERR set_time_monitor(dm, solver);
595 CHKERR TSSetSolution(solver, X);
596 CHKERR TSSetFromOptions(solver);
598 CHKERR TSSolve(solver, NULL);
600 CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
601 CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
634 const int nb_integration_pts = getGaussPts().size2();
635 const double area = getMeasure();
636 auto t_w = getFTensor0IntegrationWeight();
639 double values_integ = 0;
641 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
642 const double alpha = t_w * area;
644 values_integ += alpha * t_val;
651 std::array<double, 1> values;
652 values[0] = values_integ;
658int main(
int argc,
char *argv[]) {
665 auto core_log = logging::core::get();
674 DMType dm_name =
"DMMOFEM";
678 moab::Core mb_instance;
679 moab::Interface &moab = mb_instance;
void simple(double P1[], double P2[], double P3[], double c[], const int N)
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
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
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
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#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.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
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.
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.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
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)
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
double mu_sp
scattering coefficient (cm^-1)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
constexpr int SPACE_DIM
[Define dimension]
const double c
speed of light (cm/ns)
char init_data_file_name[255]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double mu_a
absorption coefficient (cm^-1)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
static constexpr int approx_order
Simple interface for fast problem set-up.
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
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
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)
structure for User Loop Methods on finite elements
default operator for TRI element
MoFEMErrorCode loopSideVolumes(const string fe_name, VolumeElementForcesAndSourcesCoreOnSide &fe_method)
@ OPROW
operator doWork function is executed on FE rows
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
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.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
Base volume element used to integrate on skeleton.
[Push operators to pipeline]
boost::shared_ptr< VectorDouble > approxVals
SmartPetscObj< Vec > petscVec
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)
boost::shared_ptr< VolSideFe > sideOpFe
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