13static char help[] =
"...\n\n";
53const double v =
c /
n;
119 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE],
false);
120 doEntities[MBTRI] = doEntities[MBQUAD] =
true;
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 = [&]() {
286 Range camera_surface;
287 const std::string block_name =
"CAM";
291 if (
bit->getName().compare(0, block_name.size(), block_name) == 0) {
292 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Found CAM block";
294 bit->getMeshset(), 2, camera_surface,
true);
299 MOFEM_LOG(
"PHOTON", Sev::noisy) <<
"CAM block entities:\n"
305 "PHOTON_FLUENCE_RATE");
312 auto my_simple_set_up = [&]() {
328 CHKERR set_camera_skin_fe();
329 CHKERR my_simple_set_up();
360 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"EXT",
361 "PHOTON_FLUENCE_RATE", 0, 0,
false);
366 std::string entity_name = it->getName();
367 if (entity_name.compare(0, 3,
"INT") == 0) {
369 boundary_ents,
true);
373 Range boundary_verts;
376 boundary_ents.merge(boundary_verts);
380 simple->getProblemName(),
"PHOTON_FLUENCE_RATE", boundary_ents);
391 auto add_domain_base_ops = [&](
auto &pipeline) {
392 auto jac_ptr = boost::make_shared<MatrixDouble>();
393 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
394 auto det_ptr = boost::make_shared<VectorDouble>();
401 auto add_domain_lhs_ops = [&](
auto &pipeline) {
403 "PHOTON_FLUENCE_RATE",
"PHOTON_FLUENCE_RATE",
404 [](
double,
double,
double) ->
double {
return D; }));
409 "PHOTON_FLUENCE_RATE",
"PHOTON_FLUENCE_RATE", get_mass_coefficient));
412 auto add_domain_rhs_ops = [&](
auto &pipeline) {
413 auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
414 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
415 auto dot_u_at_gauss_pts = boost::make_shared<VectorDouble>();
417 "PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts));
421 "PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts));
423 "PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts,
424 [](
double,
double,
double) ->
double {
return D; }));
426 "PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts,
427 [](
const double,
const double,
const double) {
return inv_v; }));
429 "PHOTON_FLUENCE_RATE", u_at_gauss_pts,
430 [](
const double,
const double,
const double) {
return mu_a; }));
433 auto add_boundary_base_ops = [&](
auto &pipeline) {
437 auto add_boundary_lhs_ops = [&](
auto &pipeline) {
438 for (
auto b : bc_map) {
439 if (std::regex_match(b.first, std::regex(
"(.*)EXT(.*)"))) {
441 "PHOTON_FLUENCE_RATE",
"PHOTON_FLUENCE_RATE",
443 [](
const double,
const double,
const double) {
return h; },
445 b.second->getBcEntsPtr()));
450 auto add_boundary_rhs_ops = [&](
auto &pipeline) {
451 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
454 for (
auto b : bc_map) {
455 if (std::regex_match(b.first, std::regex(
"(.*)EXT(.*)"))) {
457 "PHOTON_FLUENCE_RATE", u_at_gauss_pts,
459 [](
const double,
const double,
const double) {
return h; },
461 b.second->getBcEntsPtr()));
467 add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
468 add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
469 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
470 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
472 add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
473 add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
474 add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
475 add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
490 auto create_post_process_element = [&]() {
491 auto post_froc_fe = boost::make_shared<PostProcEle>(
mField);
492 auto u_ptr = boost::make_shared<VectorDouble>();
493 auto grad_ptr = boost::make_shared<MatrixDouble>();
494 post_froc_fe->getOpPtrVector().push_back(
496 post_froc_fe->getOpPtrVector().push_back(
499 post_froc_fe->getOpPtrVector().push_back(
new OpPPMap(
500 post_froc_fe->getPostProcMesh(), post_froc_fe->getMapGaussPts(),
501 {{
"PHOTON_FLUENCE_RATE", u_ptr}},
502 {{
"GRAD_PHOTON_FLUENCE_RATE", grad_ptr}}, {}, {}));
506 auto create_post_process_camera_element = [&]() {
509 auto post_proc_skin = boost::make_shared<PostProcFaceEle>(
mField);
511 auto u_ptr = boost::make_shared<VectorDouble>();
512 auto grad_ptr = boost::make_shared<MatrixDouble>();
517 auto jac_ptr = boost::make_shared<MatrixDouble>();
518 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
519 auto det_ptr = boost::make_shared<VectorDouble>();
522 op_loop_side->getOpPtrVector().push_back(
524 op_loop_side->getOpPtrVector().push_back(
526 op_loop_side->getOpPtrVector().push_back(
528 op_loop_side->getOpPtrVector().push_back(
530 op_loop_side->getOpPtrVector().push_back(
533 post_proc_skin->getOpPtrVector().push_back(op_loop_side);
535 post_proc_skin->getOpPtrVector().push_back(
new OpPPMap(
536 post_proc_skin->getPostProcMesh(), post_proc_skin->getMapGaussPts(),
537 {{
"PHOTON_FLUENCE_RATE", u_ptr}},
538 {{
"GRAD_PHOTON_FLUENCE_RATE", grad_ptr}}, {}, {}));
540 return post_proc_skin;
542 return boost::shared_ptr<PostProcFaceEle>();
546 auto create_post_process_integ_camera_element = [&]() {
549 auto post_proc_integ_skin = boost::make_shared<BoundaryEle>(
mField);
550 post_proc_integ_skin->getOpPtrVector().push_back(
552 post_proc_integ_skin->getOpPtrVector().push_back(
555 post_proc_integ_skin->getOpPtrVector().push_back(
558 return post_proc_integ_skin;
560 return boost::shared_ptr<BoundaryEle>();
564 auto set_time_monitor = [&](
auto dm,
auto solver) {
566 boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(
567 dm, create_post_process_element(), create_post_process_camera_element(),
569 boost::shared_ptr<ForcesAndSourcesCore> null;
571 monitor_ptr, null, null);
575 auto dm =
simple->getDM();
580 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"reading vector in binary from file "
590 auto solver = pipeline_mng->createTS();
592 CHKERR TSSetSolution(solver, X);
593 CHKERR set_time_monitor(dm, solver);
594 CHKERR TSSetSolution(solver, X);
595 CHKERR TSSetFromOptions(solver);
597 CHKERR TSSolve(solver, NULL);
599 CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
600 CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
633 const int nb_integration_pts = getGaussPts().size2();
634 const double area = getMeasure();
635 auto t_w = getFTensor0IntegrationWeight();
638 double values_integ = 0;
640 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
641 const double alpha = t_w * area;
643 values_integ += alpha * t_val;
650 std::array<double, 1> values;
651 values[0] = values_integ;
657int main(
int argc,
char *argv[]) {
660 const char param_file[] =
"param_file.petsc";
664 auto core_log = logging::core::get();
673 DMType dm_name =
"DMMOFEM";
677 moab::Core mb_instance;
678 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 >::BoundaryEle BoundaryEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#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_field)=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< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
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.
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.
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
double mu_sp
scattering coefficient (cm^-1)
constexpr int SPACE_DIM
[Define dimension]
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
char init_data_file_name[255]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
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
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.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
structure for User Loop Methods on finite elements
friend class UserDataOperator
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.
Get field gradients at integration pts for scalar field 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.
Modify integration weights on face to take in account higher-order geometry.
Set inverse jacobian to base functions.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Base volume element used to integrate on skeleton.
Volume finite element base.
[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, EntitiesFieldData::EntData &data)
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