#include <stdlib.h>
#include <cmath>
#include <MoFEM.hpp>
#include <boost/math/constants/constants.hpp>
#include <SourceFunction.hpp>
#include <boost/math/quadrature/gauss_kronrod.hpp>

Classes
struct	ElementsAndOps< DIM >
	[Operators_definition] More...

struct	PhotonDiffusion

struct	PhotonDiffusion::CommonData

struct	PhotonDiffusion::OpCameraInteg

struct	PhotonDiffusion::OpGetScalarFieldGradientValuesOnSkin

struct	PhotonDiffusion::Monitor

struct	PhotonDiffusion::OpError

Typedefs
using	DomainEle = VolumeElementForcesAndSourcesCore
	[Define dimension]

using	DomainEleOp = DomainEle::UserDataOperator

using	BoundaryEle = FaceElementForcesAndSourcesCore

using	BoundaryEleOp = BoundaryEle::UserDataOperator

using	PostProcEle = PostProcBrokenMeshInMoab< VolumeElementForcesAndSourcesCore >

using	PostProcFaceEle = PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore >

using	OpPPMap = OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM >

using	VolSideFe = VolumeElementForcesAndSourcesCoreOnSide

using	OpDomainMass = FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 >

using	OpDomainGradGrad = FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM >

using	OpDomainTimesScalarField = FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 >

using	OpDomainGradTimesVec = FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM >

using	OpDomainSource = FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 >

using	OpBoundaryMass = FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 >

using	OpBoundaryTimeScalarField = FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 >

using	OpBoundarySource = FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 >

Functions
VectorDouble	interp_function (const std::string sens_image, MatrixDouble &m_ref_coords, int nb_gauss_pts, double cam_len_x, double cam_len_y, const std::string block_name)

int	main (int argc, char *argv[])
	[Integral_calc]

Variables
static char	help [] = "...\n\n"

constexpr int	SPACE_DIM = 3
	[Define dimension]

const double	n = 1.44
	refractive index of diffusive medium

const double	c = 30.
	speed of light (cm/ns)

const double	v = c / n
	phase velocity of light in medium (cm/ns)

const double	inv_v = 1. / v

double	mu_a
	absorption coefficient (cm^-1)

double	mu_sp
	scattering coefficient (cm^-1)

double	D

double	A

double	h

double	slab_thickness

double	beam_radius

double	beam_centre_x

double	beam_centre_y

double	flux_magnitude = 1e3
	impulse magnitude

const int	kronrod_points
	This has been tested and gives the same result for any number of points. Increasing the number of points will increase the compute time, so 15 is used as default.

double	cam_len_x = 5.0

double	cam_len_y = 5.0

PetscBool	from_initial = PETSC_TRUE

PetscBool	output_volume = PETSC_FALSE

PetscBool	output_camera = PETSC_FALSE

PetscBool	testing = PETSC_FALSE

int	order = 2

int	save_every_nth_step = 1

char	init_data_file_name [255] = "init_file.dat"

char	interp_file_name [255] = "interp_sensitivity.py"

char	interp_image_name [255] = "sens_image.png"

PetscBool	enable_python = PETSC_FALSE

int	numHoLevels = 1

Typedef Documentation

◆ BoundaryEle

using BoundaryEle = FaceElementForcesAndSourcesCore

Definition at line 34 of file photon_diffusion.cpp.

◆ BoundaryEleOp

using BoundaryEleOp = BoundaryEle::UserDataOperator

Definition at line 35 of file photon_diffusion.cpp.

◆ DomainEle

using DomainEle = VolumeElementForcesAndSourcesCore

[Define dimension]

Definition at line 32 of file photon_diffusion.cpp.

◆ DomainEleOp

using DomainEleOp = DomainEle::UserDataOperator

Definition at line 33 of file photon_diffusion.cpp.

◆ OpBoundaryMass

using OpBoundaryMass = FormsIntegrators<BoundaryEleOp>::Assembly< PETSC>::BiLinearForm<GAUSS>::OpMass<1, 1>

Definition at line 55 of file photon_diffusion.cpp.

◆ OpBoundarySource

using OpBoundarySource = FormsIntegrators<BoundaryEleOp>::Assembly< PETSC>::LinearForm<GAUSS>::OpSource<1, 1>

Definition at line 59 of file photon_diffusion.cpp.

◆ OpBoundaryTimeScalarField

using OpBoundaryTimeScalarField = FormsIntegrators<BoundaryEleOp>::Assembly< PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>

Definition at line 57 of file photon_diffusion.cpp.

◆ OpDomainGradGrad

using OpDomainGradGrad = FormsIntegrators<DomainEleOp>::Assembly< PETSC>::BiLinearForm<GAUSS>::OpGradGrad<1, 1, SPACE_DIM>

Definition at line 46 of file photon_diffusion.cpp.

◆ OpDomainGradTimesVec

using OpDomainGradTimesVec = FormsIntegrators<DomainEleOp>::Assembly< PETSC>::LinearForm<GAUSS>::OpGradTimesTensor<1, 1, SPACE_DIM>

Definition at line 50 of file photon_diffusion.cpp.

◆ OpDomainMass

using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly< PETSC>::BiLinearForm<GAUSS>::OpMass<1, 1>

Definition at line 44 of file photon_diffusion.cpp.

◆ OpDomainSource

using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly< PETSC>::LinearForm<GAUSS>::OpSource<1, 1>

Definition at line 52 of file photon_diffusion.cpp.

◆ OpDomainTimesScalarField

using OpDomainTimesScalarField = FormsIntegrators<DomainEleOp>::Assembly< PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>

Definition at line 48 of file photon_diffusion.cpp.

◆ OpPPMap

using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>

Examples: /home/lk58p/mofem_install/vanilla_dev_release/mofem-cephas/mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp, PlasticOpsMonitor.hpp, SolutionMapping.hpp, mofem/atom_tests/operators_tests.cpp, mofem/atom_tests/scalar_check_approximation.cpp, mofem/atom_tests/tensor_divergence_operator.cpp, mofem/atom_tests/test_broken_space.cpp, mofem/tools/mesh_smoothing.cpp, mofem/tutorials/adv-0_plasticity/plastic.cpp, mofem/tutorials/adv-0_plasticity/src/PlasticOpsMonitor.hpp, mofem/tutorials/adv-3_level_set/level_set.cpp, mofem/tutorials/adv-4_dynamic_first_order_con_law/dynamic_first_order_con_law.cpp, mofem/tutorials/adv-5_poroelasticity/seepage.cpp, mofem/tutorials/adv-6_dg_projection/between_meshes_dg_projection.cpp, mofem/tutorials/clx-0_acoustics/helmholtz.cpp, mofem/tutorials/cor-0to1_unsaturated_transport/src/UnsaturatedFlow.hpp, mofem/tutorials/cor-10_navier_stokes/navier_stokes.cpp, mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp, mofem/tutorials/cor-2to5_analytical_poisson/src/PoissonOperators.hpp, mofem/tutorials/cor-6_simple_elasticity/simple_elasticity.cpp, mofem/tutorials/cor-7_elasticity_mixed_formulation/elasticity_mixed_formulation.cpp, mofem/tutorials/cor-9_reaction_diffusion/reaction_diffusion.cpp, mofem/tutorials/fun-2_plot_base/plot_base.cpp, mofem/tutorials/max-0_magnetostatics/src/MagneticElement.hpp, mofem/tutorials/mix-0_mixed_poisson/mixed_poisson.cpp, mofem/tutorials/mix-1_light_intensity_equation/phase.cpp, mofem/tutorials/scl-0_least_squares/approximation.cpp, mofem/tutorials/scl-12_electrostatics/electrostatics.cpp, mofem/tutorials/scl-1_poisson_2d_homogeneous/poisson_2d_homogeneous.cpp, mofem/tutorials/scl-6_heat_equation/heat_equation.cpp, mofem/tutorials/scl-7_wave_equation/wave_equation.cpp, mofem/tutorials/scl-8_heat_radiation/radiation.cpp, mofem/tutorials/scl-9_heat_method/heat_method.cpp, mofem/tutorials/vec-1_eigen_elasticity/eigen_elastic.cpp, mofem/tutorials/vec-2_nonlinear_elasticity/src/NonlinearElasticExample.hpp, mofem/tutorials/vec-3_nonlinear_dynamics/nonlinear_dynamic_elastic.cpp, mofem/tutorials/vec-4_shallow_wave/approx_sphere.cpp, mofem/tutorials/vec-4_shallow_wave/shallow_wave.cpp, mofem/tutorials/vec-5_free_surface/free_surface.cpp, mofem/tutorials/vec-6_dg_kirchhoff_plate/plate.cpp, mofem/tutorials/vec-7_shape_optimisation/adjoint.cpp, mofem/tutorials/vec-9_arc_length/src/ArcLengthExample.hpp, mofem/users_modules/basic_finite_elements/elasticity/elasticity.cpp, mofem/users_modules/basic_finite_elements/nonlinear_elastic_materials/src/NonlinearElasticElementInterface.hpp, mofem/users_modules/basic_finite_elements/nonlinear_elasticity/nonlinear_dynamics.cpp, nonlinear_elastic.cpp, plastic.cpp, thermo_elastic.cpp, and thermoplastic.cpp.

Definition at line 40 of file photon_diffusion.cpp.

◆ PostProcEle

using PostProcEle = PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>

Definition at line 36 of file photon_diffusion.cpp.

◆ PostProcFaceEle

using PostProcFaceEle = PostProcBrokenMeshInMoab<FaceElementForcesAndSourcesCore>

Definition at line 37 of file photon_diffusion.cpp.

◆ VolSideFe

using VolSideFe = VolumeElementForcesAndSourcesCoreOnSide

Definition at line 42 of file photon_diffusion.cpp.

Function Documentation

◆ interp_function()

VectorDouble interp_function	(	const std::string	sens_image,
		MatrixDouble &	m_ref_coords,
		int	nb_gauss_pts,
		double	cam_len_x,
		double	cam_len_y,
		const std::string	block_name
	)

inline

Definition at line 399 of file photon_diffusion.cpp.

                                                                {
#ifdef ENABLE_PYTHON_BINDING
  if (auto interp_ptr = interpPythonWeakPtr.lock()) {
    VectorDouble v_ref_coords = m_ref_coords.data();
 
    bp::list python_coords;
 
    for (int idx = 0; idx < 3; ++idx) {
      python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
    }
 
    np::ndarray np_interp = np::empty(bp::make_tuple(nb_gauss_pts, 3),
                                      np::dtype::get_builtin<double>());
 
    auto interp_block_name = "(.*)INTERPOLATION(.*)";
    std::regex reg_interp_name(interp_block_name);
    if (std::regex_match(block_name, reg_interp_name)) {
      CHK_MOAB_THROW(interp_ptr->evalInterp(
                         sens_image, bp::extract<np::ndarray>(python_coords[0]),
                         bp::extract<np::ndarray>(python_coords[1]), cam_len_x,
                         cam_len_y, np_interp),
                     "Failed py_Interp() python call");
    } else {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Unknown block");
    }
 
    // check the shape of returned array
    if (np_interp.get_shape()[0] != nb_gauss_pts ||
        np_interp.get_shape()[1] != 1) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Wrong shape of analytical expression returned from "
                        "python, expected: (" +
                            std::to_string(nb_gauss_pts) + ", 1), got: (" +
                            std::to_string(np_interp.get_shape()[0]) + ", " +
                            std::to_string(np_interp.get_shape()[1]) + ")");
    }
    double *interp_val_ptr = reinterpret_cast<double *>(np_interp.get_data());
 
    VectorDouble v_interp;
    v_interp.resize(nb_gauss_pts, false);
    for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
      v_interp(gg) = *(interp_val_ptr + gg);
    }
    return v_interp;
  } else {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "InterpPython pointer is expired");
  }
#endif
}

◆ main()

int main	(	int	argc,
		char *	argv[]
	)

[Integral_calc]

Definition at line 1000 of file photon_diffusion.cpp.

                                 {
  // Initialisation of MoFEM/PETSc and MOAB data structures
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  // Add logging channel for example
  auto core_log = logging::core::get();
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "PHOTON"));
  LogManager::setLog("PHOTON");
  MOFEM_LOG_TAG("PHOTON", "photon_diffusion")
 
#ifdef ENABLE_PYTHON_BINDING
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-enable_python",
                             &enable_python, //this is a flag in case
                             PETSC_NULLPTR); //the user wishes to not use python
  if (enable_python) {
    Py_Initialize();
    np::initialize();
    MOFEM_LOG("PHOTON", Sev::inform) << "Python initialised";
  } else {
    MOFEM_LOG("PHOTON", Sev::inform) << "Python NOT initialised";
  }
#endif
 
  // Error handling
  try {
    // Register MoFEM discrete manager in PETSc
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    // Create MOAB instance
    moab::Core mb_instance;              // mesh database
    moab::Interface &moab = mb_instance; // mesh database interface
 
    // Create MoFEM instance
    MoFEM::Core core(moab);           // finite element database
    MoFEM::Interface &m_field = core; // finite element interface
 
    // Run the main analysis
    PhotonDiffusion heat_problem(m_field);
    CHKERR heat_problem.runProgram();
  }
  CATCH_ERRORS;
 
  // Finish work: cleaning memory, getting statistics, etc.
 
#ifdef ENABLE_PYTHON_BINDING
  if (enable_python) {
    MOFEM_LOG("PHOTON", Sev::inform) << "Finalizing Python";
    if (Py_FinalizeEx() < 0) {
      exit(120);
    }
  }
#endif
 
  MoFEM::Core::Finalize();
 
  return 0;
}

Variable Documentation

◆ A

double A

Definition at line 70 of file photon_diffusion.cpp.

◆ beam_centre_x

double beam_centre_x

Definition at line 75 of file photon_diffusion.cpp.

◆ beam_centre_y

double beam_centre_y

Definition at line 76 of file photon_diffusion.cpp.

◆ beam_radius

double beam_radius

Definition at line 74 of file photon_diffusion.cpp.

◆ c

const double c = 30.

speed of light (cm/ns)

Definition at line 63 of file photon_diffusion.cpp.

◆ cam_len_x

double cam_len_x = 5.0

Definition at line 82 of file photon_diffusion.cpp.

◆ cam_len_y

double cam_len_y = 5.0

Definition at line 83 of file photon_diffusion.cpp.

◆ D

double D

Definition at line 69 of file photon_diffusion.cpp.

◆ enable_python

PetscBool enable_python = PETSC_FALSE

Definition at line 96 of file photon_diffusion.cpp.

◆ flux_magnitude

double flux_magnitude = 1e3

impulse magnitude

Definition at line 77 of file photon_diffusion.cpp.

◆ from_initial

PetscBool from_initial = PETSC_TRUE

Definition at line 85 of file photon_diffusion.cpp.

◆ h

double h

Examples: /home/lk58p/mofem_install/vanilla_dev_release/mofem-cephas/mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp, and mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 71 of file photon_diffusion.cpp.

◆ help

char help[] = "...\n\n"

static

Definition at line 24 of file photon_diffusion.cpp.

◆ init_data_file_name

char init_data_file_name[255] = "init_file.dat"

Definition at line 93 of file photon_diffusion.cpp.

◆ interp_file_name

char interp_file_name[255] = "interp_sensitivity.py"

Definition at line 94 of file photon_diffusion.cpp.

◆ interp_image_name

char interp_image_name[255] = "sens_image.png"

Definition at line 95 of file photon_diffusion.cpp.

◆ inv_v

const double inv_v = 1. / v

Definition at line 65 of file photon_diffusion.cpp.

◆ kronrod_points

const int kronrod_points

Initial value:

=

15

This has been tested and gives the same result for any number of points. Increasing the number of points will increase the compute time, so 15 is used as default.

number of points for kronrod integration, can be 15, 31, 41, 51, or 61 (from boost library docs)

Definition at line 78 of file photon_diffusion.cpp.

◆ mu_a

double mu_a

absorption coefficient (cm^-1)

Definition at line 67 of file photon_diffusion.cpp.

◆ mu_sp

double mu_sp

scattering coefficient (cm^-1)

Definition at line 68 of file photon_diffusion.cpp.

◆ n

const double n = 1.44

refractive index of diffusive medium

Definition at line 62 of file photon_diffusion.cpp.

◆ numHoLevels

int numHoLevels = 1

Definition at line 98 of file photon_diffusion.cpp.

◆ order

int order = 2

Definition at line 90 of file photon_diffusion.cpp.

◆ output_camera

PetscBool output_camera = PETSC_FALSE

Definition at line 87 of file photon_diffusion.cpp.

◆ output_volume

PetscBool output_volume = PETSC_FALSE

Definition at line 86 of file photon_diffusion.cpp.

◆ save_every_nth_step

int save_every_nth_step = 1

Examples: mofem/tutorials/adv-4_dynamic_first_order_con_law/dynamic_first_order_con_law.cpp, mofem/tutorials/vec-3_nonlinear_dynamics/nonlinear_dynamic_elastic.cpp, mofem/tutorials/vec-4_shallow_wave/shallow_wave.cpp, mofem/tutorials/vec-5_free_surface/free_surface.cpp, mofem/users_modules/adolc-plasticity/adolc_plasticity.cpp, and nonlinear_elastic.cpp.

Definition at line 91 of file photon_diffusion.cpp.

◆ slab_thickness

double slab_thickness

Definition at line 73 of file photon_diffusion.cpp.

◆ SPACE_DIM

constexpr int SPACE_DIM = 3

constexpr