11 using namespace MoFEM;
13 static char help[] =
"...\n\n";
38 PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
42 PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
47 PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
49 PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
51 const double n = 1.44;
53 const double v =
c /
n;
118 commonDataPtr(common_data_ptr) {
119 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE],
false);
120 doEntities[MBTRI] = doEntities[MBQUAD] =
true;
130 :
BoundaryEleOp(
"PHOTON_FLUENCE_RATE", OPROW), sideOpFe(side_fe) {}
135 if (
type != MBVERTEX)
137 CHKERR loopSideVolumes(
"dFE", *sideOpFe);
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)
148 : dM(dm), postProc(post_proc), skinPostProc(skin_post_proc),
149 skinPostProcInteg(skin_post_proc_integ),
150 commonDataPtr(common_data_ptr){};
164 CHKERR VecZeroEntries(commonDataPtr->petscVec);
165 CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
167 CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
170 CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
171 CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
172 CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, ADD_VALUES,
174 CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, ADD_VALUES,
176 CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
178 CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
181 CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
182 MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Fluence rate integral: " << array[0];
187 CHKERR postProc->writeFile(
"out_volume_" +
188 boost::lexical_cast<std::string>(ts_step) +
193 CHKERR skinPostProc->writeFile(
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 = [&]() {
547 if (mField.check_finite_element(
"CAMERA_FE")) {
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(
554 commonDataPtr->approxVals));
555 post_proc_integ_skin->getOpPtrVector().push_back(
556 new OpCameraInteg(commonDataPtr));
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(),
568 create_post_process_integ_camera_element(), commonDataPtr));
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;
657 int main(
int argc,
char *argv[]) {
660 const char param_file[] =
"param_file.petsc";
664 auto core_log = logging::core::get();
666 LogManager::createSink(LogManager::getStrmWorld(),
"PHOTON"));
667 LogManager::setLog(
"PHOTON");
673 DMType dm_name =
"DMMOFEM";