v0.14.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>
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
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, 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
133 DataForcesAndSourcesCore::EntData &data) {
135 if (type != MBVERTEX)
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>();
216 CHKERR 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_NULL, "", "-initial_file",
247 init_data_file_name, 255, PETSC_NULL);
248
249 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-from_initial", &from_initial,
250 PETSC_NULL);
251 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-output_volume", &output_volume,
252 PETSC_NULL);
253 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-output_camera", &output_camera,
254 PETSC_NULL);
255 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-mu_a", &mu_a, PETSC_NULL);
256 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-mu_sp", &mu_sp, PETSC_NULL);
257 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-coef_A", &A, PETSC_NULL);
258
259 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
260 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_step", &save_every_nth_step,
261 PETSC_NULL);
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 auto meshset_mng = mField.getInterface<MeshsetsManager>();
286
287 Range camera_surface;
288 const std::string block_name = "CAM";
289 bool add_fe = false;
290
292 if (bit->getName().compare(0, block_name.size(), block_name) == 0) {
293 MOFEM_LOG("PHOTON", Sev::inform) << "Found CAM block";
294 CHKERR mField.get_moab().get_entities_by_dimension(
295 bit->getMeshset(), 2, camera_surface, true);
296 add_fe = true;
297 }
298 }
299
300 MOFEM_LOG("PHOTON", Sev::noisy) << "CAM block entities:\n"
301 << camera_surface;
302
303 if (add_fe) {
304 CHKERR mField.add_finite_element("CAMERA_FE");
306 "PHOTON_FLUENCE_RATE");
308 "CAMERA_FE");
309 }
311 };
312
313 auto my_simple_set_up = [&]() {
315 CHKERR simple->defineFiniteElements();
316 CHKERR simple->defineProblem(PETSC_TRUE);
317 CHKERR simple->buildFields();
318 CHKERR simple->buildFiniteElements();
319
320 if (mField.check_finite_element("CAMERA_FE")) {
322 CHKERR DMMoFEMAddElement(simple->getDM(), "CAMERA_FE");
323 }
324
325 CHKERR simple->buildProblem();
327 };
328
329 CHKERR set_camera_skin_fe();
330 CHKERR my_simple_set_up();
331
333}
334
337
338 auto integration_rule = [](int o_row, int o_col, int approx_order) {
339 return 2 * approx_order;
340 };
341
342 auto *pipeline_mng = mField.getInterface<PipelineManager>();
344 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
345 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
346 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
347
349}
350
353
355}
356
359 auto bc_mng = mField.getInterface<BcManager>();
360 auto *simple = mField.getInterface<Simple>();
361 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "EXT",
362 "PHOTON_FLUENCE_RATE", 0, 0, false);
363
364 // Get boundary edges marked in block named "BOUNDARY_CONDITION"
365 Range boundary_ents;
367 std::string entity_name = it->getName();
368 if (entity_name.compare(0, 3, "INT") == 0) {
369 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
370 boundary_ents, true);
371 }
372 }
373 // Add vertices to boundary entities
374 Range boundary_verts;
375 CHKERR mField.get_moab().get_connectivity(boundary_ents, boundary_verts,
376 true);
377 boundary_ents.merge(boundary_verts);
378
379 // Remove DOFs as homogeneous boundary condition is used
380 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
381 simple->getProblemName(), "PHOTON_FLUENCE_RATE", boundary_ents);
382
384}
385
388
389 auto bc_mng = mField.getInterface<BcManager>();
390 auto &bc_map = bc_mng->getBcMapByBlockName();
391
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>();
396 pipeline.push_back(new OpCalculateHOJac<3>(jac_ptr));
397 pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
398 pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(H1, inv_jac_ptr));
399 pipeline.push_back(new OpSetHOWeights(det_ptr));
400 };
401
402 auto add_domain_lhs_ops = [&](auto &pipeline) {
403 pipeline.push_back(new OpDomainGradGrad(
404 "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
405 [](double, double, double) -> double { return D; }));
406 auto get_mass_coefficient = [&](const double, const double, const double) {
407 return inv_v * domainLhsFEPtr->ts_a + mu_a;
408 };
409 pipeline.push_back(new OpDomainMass(
410 "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE", get_mass_coefficient));
411 };
412
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>();
417 pipeline.push_back(new OpCalculateScalarFieldGradient<SPACE_DIM>(
418 "PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts));
419 pipeline.push_back(new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
420 u_at_gauss_pts));
421 pipeline.push_back(new OpCalculateScalarFieldValuesDot(
422 "PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts));
423 pipeline.push_back(new OpDomainGradTimesVec(
424 "PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts,
425 [](double, double, double) -> double { return D; }));
426 pipeline.push_back(new OpDomainTimesScalarField(
427 "PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts,
428 [](const double, const double, const double) { return inv_v; }));
429 pipeline.push_back(new OpDomainTimesScalarField(
430 "PHOTON_FLUENCE_RATE", u_at_gauss_pts,
431 [](const double, const double, const double) { return mu_a; }));
432 };
433
434 auto add_boundary_base_ops = [&](auto &pipeline) {
435 pipeline.push_back(new OpSetHOWeightsOnFace());
436 };
437
438 auto add_boundary_lhs_ops = [&](auto &pipeline) {
439 for (auto b : bc_map) {
440 if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
441 pipeline.push_back(new OpBoundaryMass(
442 "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
443
444 [](const double, const double, const double) { return h; },
445
446 b.second->getBcEntsPtr()));
447 }
448 }
449 };
450
451 auto add_boundary_rhs_ops = [&](auto &pipeline) {
452 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
453 pipeline.push_back(new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
454 u_at_gauss_pts));
455 for (auto b : bc_map) {
456 if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
457 pipeline.push_back(new OpBoundaryTimeScalarField(
458 "PHOTON_FLUENCE_RATE", u_at_gauss_pts,
459
460 [](const double, const double, const double) { return h; },
461
462 b.second->getBcEntsPtr()));
463 }
464 }
465 };
466
467 auto pipeline_mng = mField.getInterface<PipelineManager>();
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());
472
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());
477
478 domainLhsFEPtr = pipeline_mng->getDomainLhsFE();
479 boundaryLhsFEPtr = pipeline_mng->getBoundaryLhsFE();
480 boundaryRhsFEPtr = pipeline_mng->getBoundaryRhsFE();
481
483}
484
487
488 auto *simple = mField.getInterface<Simple>();
489 auto *pipeline_mng = mField.getInterface<PipelineManager>();
490
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(
496 new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE", u_ptr));
497 post_froc_fe->getOpPtrVector().push_back(
498 new OpCalculateScalarFieldGradient<SPACE_DIM>("PHOTON_FLUENCE_RATE",
499 grad_ptr));
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}}, {}, {}));
504 return post_froc_fe;
505 };
506
507 auto create_post_process_camera_element = [&]() {
508 if (mField.check_finite_element("CAMERA_FE")) {
509
510 auto post_proc_skin = boost::make_shared<PostProcFaceEle>(mField);
511
512 auto u_ptr = boost::make_shared<VectorDouble>();
513 auto grad_ptr = boost::make_shared<MatrixDouble>();
514
516 mField, simple->getDomainFEName(), SPACE_DIM);
517
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>();
521
522 // push operators to side element
523 op_loop_side->getOpPtrVector().push_back(
524 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
525 op_loop_side->getOpPtrVector().push_back(
526 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
527 op_loop_side->getOpPtrVector().push_back(
528 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
529 op_loop_side->getOpPtrVector().push_back(
530 new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE", u_ptr));
531 op_loop_side->getOpPtrVector().push_back(
532 new OpCalculateScalarFieldGradient<SPACE_DIM>("PHOTON_FLUENCE_RATE", grad_ptr));
533 // push op to boundary element
534 post_proc_skin->getOpPtrVector().push_back(op_loop_side);
535
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}}, {}, {}));
540
541 return post_proc_skin;
542 } else {
543 return boost::shared_ptr<PostProcFaceEle>();
544 }
545 };
546
547 auto create_post_process_integ_camera_element = [&]() {
548 if (mField.check_finite_element("CAMERA_FE")) {
549
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(
554 new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
555 commonDataPtr->approxVals));
556 post_proc_integ_skin->getOpPtrVector().push_back(
557 new OpCameraInteg(commonDataPtr));
558
559 return post_proc_integ_skin;
560 } else {
561 return boost::shared_ptr<BoundaryEle>();
562 }
563 };
564
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(),
569 create_post_process_integ_camera_element(), commonDataPtr));
570 boost::shared_ptr<ForcesAndSourcesCore> null;
571 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
572 monitor_ptr, null, null);
574 };
575
576 auto dm = simple->getDM();
577 auto X = createDMVector(dm);
578
579 if (from_initial) {
580
581 MOFEM_LOG("PHOTON", Sev::inform) << "reading vector in binary from file "
582 << init_data_file_name << " ...";
583 PetscViewer viewer;
584 PetscViewerBinaryOpen(PETSC_COMM_WORLD, init_data_file_name, FILE_MODE_READ,
585 &viewer);
586 VecLoad(X, viewer);
587
588 CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_REVERSE);
589 }
590
591 auto solver = pipeline_mng->createTS();
592
593 CHKERR TSSetSolution(solver, X);
594 CHKERR set_time_monitor(dm, solver);
595 CHKERR TSSetSolution(solver, X);
596 CHKERR TSSetFromOptions(solver);
597 CHKERR TSSetUp(solver);
598 CHKERR TSSolve(solver, NULL);
599
600 CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
601 CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
602 CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_REVERSE);
603
605}
606
609
610 // Processes to set output results are integrated in solveSystem()
611
613}
614
617
627
629}
630
632 EntData &data) {
634 const int nb_integration_pts = getGaussPts().size2();
635 const double area = getMeasure();
636 auto t_w = getFTensor0IntegrationWeight();
637 auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
638
639 double values_integ = 0;
640
641 for (int gg = 0; gg != nb_integration_pts; ++gg) {
642 const double alpha = t_w * area;
643
644 values_integ += alpha * t_val;
645
646 ++t_w;
647 ++t_val;
648 }
649
650 constexpr std::array<int, 1> indices = {CommonData::VALUES_INTEG};
651 std::array<double, 1> values;
652 values[0] = values_integ;
653 CHKERR VecSetValues(commonDataPtr->petscVec, 1, indices.data(), values.data(),
654 ADD_VALUES);
656}
657
658int main(int argc, char *argv[]) {
659
660 // Initialisation of MoFEM/PETSc and MOAB data structures
661 const char param_file[] = "param_file.petsc";
663
664 // Add logging channel for example
665 auto core_log = logging::core::get();
666 core_log->add_sink(
668 LogManager::setLog("PHOTON");
669 MOFEM_LOG_TAG("PHOTON", "photon_diffusion")
670
671 // Error handling
672 try {
673 // Register MoFEM discrete manager in PETSc
674 DMType dm_name = "DMMOFEM";
675 CHKERR DMRegister_MoFEM(dm_name);
676
677 // Create MOAB instance
678 moab::Core mb_instance; // mesh database
679 moab::Interface &moab = mb_instance; // mesh database interface
680
681 // Create MoFEM instance
682 MoFEM::Core core(moab); // finite element database
683 MoFEM::Interface &m_field = core; // finite element interface
684
685 // Run the main analysis
686 PhotonDiffusion heat_problem(m_field);
687 CHKERR heat_problem.runProgram();
688 }
690
691 // Finish work: cleaning memory, getting statistics, etc.
693
694 return 0;
695}
std::string param_file
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
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.
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
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
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.
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
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)
structure for User Loop Methods on finite elements
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.
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