v0.15.0
Loading...
Searching...
No Matches
photon_diffusion.cpp
Go to the documentation of this file.
1/**
2 * \file photon_diffusion.cpp
3 * \example photon_diffusion.cpp
4 *
5 **/
6
7#include <stdlib.h>
8#include <cmath>
9#include <MoFEM.hpp>
10
11using namespace MoFEM;
12
13static char help[] = "...\n\n";
14
15template <int DIM> struct ElementsAndOps {};
16
17//! [Define dimension]
18constexpr int SPACE_DIM = 3; //< Space dimension of problem, mesh
19//! [Define dimension]
20
22using DomainEleOp = DomainEle::UserDataOperator;
28
30
32
43
50
51const double n = 1.44; ///< refractive index of diffusive medium
52const double c = 30.; ///< speed of light (cm/ns)
53const double v = c / n; ///< phase velocity of light in medium (cm/ns)
54const double inv_v = 1. / v;
55
56double mu_a; ///< absorption coefficient (cm^-1)
57double mu_sp; ///< scattering coefficient (cm^-1)
58double D;
59double A;
60double h;
61
62PetscBool from_initial = PETSC_TRUE;
63PetscBool output_volume = PETSC_FALSE;
64PetscBool output_camera = PETSC_FALSE;
65
66int order = 2;
68
69char init_data_file_name[255] = "init_file.dat";
71
72
73struct PhotonDiffusion {
74public:
76
77 // Declaration of the main function to run analysis
79
80private:
81 // Declaration of other main functions called in runProgram()
91
92 // Main interfaces
94
95 // Object to mark boundary entities for the assembling of domain elements
96 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
97
98 boost::shared_ptr<FEMethod> domainLhsFEPtr;
99 boost::shared_ptr<FEMethod> boundaryLhsFEPtr;
100 boost::shared_ptr<FEMethod> boundaryRhsFEPtr;
101
102 struct CommonData {
103 boost::shared_ptr<VectorDouble> approxVals;
105
110 };
111
112 boost::shared_ptr<CommonData> commonDataPtr;
113
114 struct OpCameraInteg : public BoundaryEleOp {
115 boost::shared_ptr<CommonData> commonDataPtr;
116 OpCameraInteg(boost::shared_ptr<CommonData> common_data_ptr)
117 : BoundaryEleOp("PHOTON_FLUENCE_RATE", OPROW),
118 commonDataPtr(common_data_ptr) {
119 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
120 doEntities[MBTRI] = doEntities[MBQUAD] = true;
121 }
122 MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data);
123 };
124
126
127 boost::shared_ptr<VolSideFe> sideOpFe;
128
129 OpGetScalarFieldGradientValuesOnSkin(boost::shared_ptr<VolSideFe> side_fe)
130 : BoundaryEleOp("PHOTON_FLUENCE_RATE", OPROW), sideOpFe(side_fe) {}
131
132 MoFEMErrorCode doWork(int side, EntityType type,
133 DataForcesAndSourcesCore::EntData &data) {
135 if (type != MBVERTEX)
137 CHKERR loopSideVolumes("dFE", *sideOpFe);
139 }
140 };
141
142 struct Monitor : public FEMethod {
143
144 Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc,
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){};
151
160
163
164 CHKERR VecZeroEntries(commonDataPtr->petscVec);
165 CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
166 SCATTER_FORWARD);
167 CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
168 SCATTER_FORWARD);
170 CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
171 CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
172 CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, ADD_VALUES,
173 SCATTER_REVERSE);
174 CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, ADD_VALUES,
175 SCATTER_REVERSE);
176 CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
177 SCATTER_FORWARD);
178 CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
179 SCATTER_FORWARD);
180 const double *array;
181 CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
182 MOFEM_LOG("PHOTON", Sev::inform) << "Fluence rate integral: " << array[0];
183
184 if (ts_step % save_every_nth_step == 0) {
185 if (output_volume) {
187 CHKERR postProc->writeFile("out_volume_" +
188 boost::lexical_cast<std::string>(ts_step) +
189 ".h5m");
190 }
193 CHKERR skinPostProc->writeFile(
194 "out_camera_" + boost::lexical_cast<std::string>(ts_step) +
195 ".h5m");
196 }
197 }
199 }
200
201 private:
203 boost::shared_ptr<PostProcEle> postProc;
204 boost::shared_ptr<PostProcFaceEle> skinPostProc;
205 boost::shared_ptr<BoundaryEle> skinPostProcInteg;
206 boost::shared_ptr<CommonData> commonDataPtr;
207 };
208};
209
210PhotonDiffusion::PhotonDiffusion(MoFEM::Interface &m_field) : mField(m_field) {}
211
214
215 auto *simple = mField.getInterface<Simple>();
217 CHKERR simple->getOptions();
218 CHKERR simple->loadFile();
219
221}
222
225 commonDataPtr = boost::make_shared<CommonData>();
226 PetscInt ghosts[1] = {0};
227 if (!mField.get_comm_rank())
228 commonDataPtr->petscVec =
229 createGhostVector(mField.get_comm(), 1, 1, 0, ghosts);
230 else
231 commonDataPtr->petscVec =
232 createGhostVector(mField.get_comm(), 0, 1, 1, ghosts);
233 commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
235}
236
239
240 auto *simple = mField.getInterface<Simple>();
241 CHKERR simple->addDomainField("PHOTON_FLUENCE_RATE", H1,
243 CHKERR simple->addBoundaryField("PHOTON_FLUENCE_RATE", H1,
245
246 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-initial_file",
247 init_data_file_name, 255, PETSC_NULLPTR);
248
249 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-from_initial", &from_initial,
250 PETSC_NULLPTR);
251 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-output_volume", &output_volume,
252 PETSC_NULLPTR);
253 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-output_camera", &output_camera,
254 PETSC_NULLPTR);
255 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-mu_a", &mu_a, PETSC_NULLPTR);
256 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-mu_sp", &mu_sp, PETSC_NULLPTR);
257 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-coef_A", &A, PETSC_NULLPTR);
258
259 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
260 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-save_step", &save_every_nth_step,
261 PETSC_NULLPTR);
262
263 h = 0.5 / A;
264 D = 1. / (3. * (mu_a + mu_sp));
265
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;
270 MOFEM_LOG("PHOTON", Sev::inform)
271 << "Absorption coefficient (cm^-1): " << mu_a;
272 MOFEM_LOG("PHOTON", Sev::inform)
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;
277
278 MOFEM_LOG("PHOTON", Sev::inform) << "Approximation order: " << order;
279 MOFEM_LOG("PHOTON", Sev::inform) << "Save step: " << save_every_nth_step;
280
281 CHKERR simple->setFieldOrder("PHOTON_FLUENCE_RATE", order);
282
283 auto set_camera_skin_fe = [&]() {
285
286 Range camera_surface;
287 const std::string block_name = "CAM";
288 bool add_fe = false;
289
291 if (bit->getName().compare(0, block_name.size(), block_name) == 0) {
292 MOFEM_LOG("PHOTON", Sev::inform) << "Found CAM block";
293 CHKERR mField.get_moab().get_entities_by_dimension(
294 bit->getMeshset(), 2, camera_surface, true);
295 add_fe = true;
296 }
297 }
298
299 MOFEM_LOG("PHOTON", Sev::noisy) << "CAM block entities:\n"
300 << camera_surface;
301
302 if (add_fe) {
303 CHKERR mField.add_finite_element("CAMERA_FE");
305 "PHOTON_FLUENCE_RATE");
307 "CAMERA_FE");
308 }
310 };
311
312 auto my_simple_set_up = [&]() {
314 CHKERR simple->defineFiniteElements();
315 CHKERR simple->defineProblem(PETSC_TRUE);
316 CHKERR simple->buildFields();
317 CHKERR simple->buildFiniteElements();
318
319 if (mField.check_finite_element("CAMERA_FE")) {
321 CHKERR DMMoFEMAddElement(simple->getDM(), "CAMERA_FE");
322 }
323
324 CHKERR simple->buildProblem();
326 };
327
328 CHKERR set_camera_skin_fe();
329 CHKERR my_simple_set_up();
330
332}
333
336
337 auto integration_rule = [](int o_row, int o_col, int approx_order) {
338 return 2 * approx_order;
339 };
340
341 auto *pipeline_mng = mField.getInterface<PipelineManager>();
343 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
344 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
345 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
346
348}
349
352
354}
355
358 auto bc_mng = mField.getInterface<BcManager>();
359 auto *simple = mField.getInterface<Simple>();
360 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "EXT",
361 "PHOTON_FLUENCE_RATE", 0, 0, false);
362
363 // Get boundary edges marked in block named "BOUNDARY_CONDITION"
364 Range boundary_ents;
366 std::string entity_name = it->getName();
367 if (entity_name.compare(0, 3, "INT") == 0) {
368 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
369 boundary_ents, true);
370 }
371 }
372 // Add vertices to boundary entities
373 Range boundary_verts;
374 CHKERR mField.get_moab().get_connectivity(boundary_ents, boundary_verts,
375 true);
376 boundary_ents.merge(boundary_verts);
377
378 // Remove DOFs as homogeneous boundary condition is used
379 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
380 simple->getProblemName(), "PHOTON_FLUENCE_RATE", boundary_ents);
381
383}
384
387
388 auto bc_mng = mField.getInterface<BcManager>();
389 auto &bc_map = bc_mng->getBcMapByBlockName();
390
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>();
395 pipeline.push_back(new OpCalculateHOJac<3>(jac_ptr));
396 pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
397 pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(H1, inv_jac_ptr));
398 pipeline.push_back(new OpSetHOWeights(det_ptr));
399 };
400
401 auto add_domain_lhs_ops = [&](auto &pipeline) {
402 pipeline.push_back(new OpDomainGradGrad(
403 "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
404 [](double, double, double) -> double { return D; }));
405 auto get_mass_coefficient = [&](const double, const double, const double) {
406 return inv_v * domainLhsFEPtr->ts_a + mu_a;
407 };
408 pipeline.push_back(new OpDomainMass(
409 "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE", get_mass_coefficient));
410 };
411
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>();
416 pipeline.push_back(new OpCalculateScalarFieldGradient<SPACE_DIM>(
417 "PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts));
418 pipeline.push_back(new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
419 u_at_gauss_pts));
420 pipeline.push_back(new OpCalculateScalarFieldValuesDot(
421 "PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts));
422 pipeline.push_back(new OpDomainGradTimesVec(
423 "PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts,
424 [](double, double, double) -> double { return D; }));
425 pipeline.push_back(new OpDomainTimesScalarField(
426 "PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts,
427 [](const double, const double, const double) { return inv_v; }));
428 pipeline.push_back(new OpDomainTimesScalarField(
429 "PHOTON_FLUENCE_RATE", u_at_gauss_pts,
430 [](const double, const double, const double) { return mu_a; }));
431 };
432
433 auto add_boundary_base_ops = [&](auto &pipeline) {
434 pipeline.push_back(new OpSetHOWeightsOnFace());
435 };
436
437 auto add_boundary_lhs_ops = [&](auto &pipeline) {
438 for (auto b : bc_map) {
439 if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
440 pipeline.push_back(new OpBoundaryMass(
441 "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
442
443 [](const double, const double, const double) { return h; },
444
445 b.second->getBcEntsPtr()));
446 }
447 }
448 };
449
450 auto add_boundary_rhs_ops = [&](auto &pipeline) {
451 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
452 pipeline.push_back(new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
453 u_at_gauss_pts));
454 for (auto b : bc_map) {
455 if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
456 pipeline.push_back(new OpBoundaryTimeScalarField(
457 "PHOTON_FLUENCE_RATE", u_at_gauss_pts,
458
459 [](const double, const double, const double) { return h; },
460
461 b.second->getBcEntsPtr()));
462 }
463 }
464 };
465
466 auto pipeline_mng = mField.getInterface<PipelineManager>();
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());
471
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());
476
477 domainLhsFEPtr = pipeline_mng->getDomainLhsFE();
478 boundaryLhsFEPtr = pipeline_mng->getBoundaryLhsFE();
479 boundaryRhsFEPtr = pipeline_mng->getBoundaryRhsFE();
480
482}
483
486
487 auto *simple = mField.getInterface<Simple>();
488 auto *pipeline_mng = mField.getInterface<PipelineManager>();
489
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(
495 new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE", u_ptr));
496 post_froc_fe->getOpPtrVector().push_back(
497 new OpCalculateScalarFieldGradient<SPACE_DIM>("PHOTON_FLUENCE_RATE",
498 grad_ptr));
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}}, {}, {}));
503 return post_froc_fe;
504 };
505
506 auto create_post_process_camera_element = [&]() {
507 if (mField.check_finite_element("CAMERA_FE")) {
508
509 auto post_proc_skin = boost::make_shared<PostProcFaceEle>(mField);
510
511 auto u_ptr = boost::make_shared<VectorDouble>();
512 auto grad_ptr = boost::make_shared<MatrixDouble>();
513
515 mField, simple->getDomainFEName(), SPACE_DIM);
516
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>();
520
521 // push operators to side element
522 op_loop_side->getOpPtrVector().push_back(
523 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
524 op_loop_side->getOpPtrVector().push_back(
525 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
526 op_loop_side->getOpPtrVector().push_back(
527 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
528 op_loop_side->getOpPtrVector().push_back(
529 new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE", u_ptr));
530 op_loop_side->getOpPtrVector().push_back(
531 new OpCalculateScalarFieldGradient<SPACE_DIM>("PHOTON_FLUENCE_RATE", grad_ptr));
532 // push op to boundary element
533 post_proc_skin->getOpPtrVector().push_back(op_loop_side);
534
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}}, {}, {}));
539
540 return post_proc_skin;
541 } else {
542 return boost::shared_ptr<PostProcFaceEle>();
543 }
544 };
545
546 auto create_post_process_integ_camera_element = [&]() {
547 if (mField.check_finite_element("CAMERA_FE")) {
548
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(
553 new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
554 commonDataPtr->approxVals));
555 post_proc_integ_skin->getOpPtrVector().push_back(
556 new OpCameraInteg(commonDataPtr));
557
558 return post_proc_integ_skin;
559 } else {
560 return boost::shared_ptr<BoundaryEle>();
561 }
562 };
563
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;
570 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
571 monitor_ptr, null, null);
573 };
574
575 auto dm = simple->getDM();
576 auto X = createDMVector(dm);
577
578 if (from_initial) {
579
580 MOFEM_LOG("PHOTON", Sev::inform) << "reading vector in binary from file "
581 << init_data_file_name << " ...";
582 PetscViewer viewer;
583 PetscViewerBinaryOpen(PETSC_COMM_WORLD, init_data_file_name, FILE_MODE_READ,
584 &viewer);
585 VecLoad(X, viewer);
586
587 CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_REVERSE);
588 }
589
590 auto solver = pipeline_mng->createTS();
591
592 CHKERR TSSetSolution(solver, X);
593 CHKERR set_time_monitor(dm, solver);
594 CHKERR TSSetSolution(solver, X);
595 CHKERR TSSetFromOptions(solver);
596 CHKERR TSSetUp(solver);
597 CHKERR TSSolve(solver, NULL);
598
599 CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
600 CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
601 CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_REVERSE);
602
604}
605
608
609 // Processes to set output results are integrated in solveSystem()
610
612}
613
616
626
628}
629
633 const int nb_integration_pts = getGaussPts().size2();
634 const double area = getMeasure();
635 auto t_w = getFTensor0IntegrationWeight();
636 auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
637
638 double values_integ = 0;
639
640 for (int gg = 0; gg != nb_integration_pts; ++gg) {
641 const double alpha = t_w * area;
642
643 values_integ += alpha * t_val;
644
645 ++t_w;
646 ++t_val;
647 }
648
649 constexpr std::array<int, 1> indices = {CommonData::VALUES_INTEG};
650 std::array<double, 1> values;
651 values[0] = values_integ;
652 CHKERR VecSetValues(commonDataPtr->petscVec, 1, indices.data(), values.data(),
653 ADD_VALUES);
655}
656
657int main(int argc, char *argv[]) {
658
659 // Initialisation of MoFEM/PETSc and MOAB data structures
660 const char param_file[] = "param_file.petsc";
661 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
662
663 // Add logging channel for example
664 auto core_log = logging::core::get();
665 core_log->add_sink(
667 LogManager::setLog("PHOTON");
668 MOFEM_LOG_TAG("PHOTON", "photon_diffusion")
669
670 // Error handling
671 try {
672 // Register MoFEM discrete manager in PETSc
673 DMType dm_name = "DMMOFEM";
674 CHKERR DMRegister_MoFEM(dm_name);
675
676 // Create MOAB instance
677 moab::Core mb_instance; // mesh database
678 moab::Interface &moab = mb_instance; // mesh database interface
679
680 // Create MoFEM instance
681 MoFEM::Core core(moab); // finite element database
682 MoFEM::Interface &m_field = core; // finite element interface
683
684 // Run the main analysis
685 PhotonDiffusion heat_problem(m_field);
686 CHKERR heat_problem.runProgram();
687 }
689
690 // Finish work: cleaning memory, getting statistics, etc.
692
693 return 0;
694}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
int main()
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.
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.
auto integration_rule
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
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:514
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
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:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
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.
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
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
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:1046
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)
double A
static char help[]
int numHoLevels
constexpr int SPACE_DIM
[Define dimension]
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
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
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:118
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
structure for User Loop Methods on finite elements
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.
Definition Simple.hpp:27
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.
Definition Simple.cpp:261
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference 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, EntitiesFieldData::EntData &data)
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