photon_diffusion.cpp

Go to the documentation of this file.

/**
 * \file photon_diffusion.cpp
 * \example mofem/tutorials/scl-10/photon_diffusion.cpp
 *
 **/
 
#include <stdlib.h>
#include <cmath>
#include <MoFEM.hpp>
#include <boost/math/constants/constants.hpp>
#include <SourceFunction.hpp>
 
#ifdef ENABLE_PYTHON_BINDING
  #include <boost/python.hpp>
  #include <boost/python/def.hpp>
  #include <boost/python/numpy.hpp>
 
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
 
using namespace MoFEM;
 
static char help[] = "...\n\n";
 
template <int DIM> struct ElementsAndOps {};
 
//! [Define dimension]
constexpr int SPACE_DIM = 3; //< Space dimension of problem, mesh
//! [Define dimension]
 
using DomainEle = VolumeElementForcesAndSourcesCore;
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>;
 
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; //< spot radius
double beam_centre_x;
double beam_centre_y;
double flux_magnitude = 1e3; ///< impulse magnitude
const int kronrod_points =
    15; ///< number of points for kronrod integration, can be 15, 31, 41, 51, or 61 (from boost library docs)
///< 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;
 
#include <boost/math/quadrature/gauss_kronrod.hpp>
using namespace boost::math::quadrature;
 
struct PhotonDiffusion {
public:
  PhotonDiffusion(MoFEM::Interface &m_field);
 
  // Declaration of the main function to run analysis
  MoFEMErrorCode runProgram();
  struct InterpPython;
 
private:
  // Declaration of other main functions called in runProgram()
  MoFEMErrorCode readMesh();
  MoFEMErrorCode createCommonData();
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode setIntegrationRules();
  MoFEMErrorCode initialCondition();
  MoFEMErrorCode boundaryCondition();
  MoFEMErrorCode assembleSystem();
  MoFEMErrorCode solveSystem();
  MoFEMErrorCode checkResults();
  MoFEMErrorCode outputResults();
 
  // Main interfaces
  MoFEM::Interface &mField;
 
  // Object to mark boundary entities for the assembling of domain elements
  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
 
  boost::shared_ptr<FEMethod> domainLhsFEPtr;
  boost::shared_ptr<FEMethod> boundaryLhsFEPtr;
  boost::shared_ptr<FEMethod> boundaryRhsFEPtr;
 
#ifdef ENABLE_PYTHON_BINDING
  boost::shared_ptr<InterpPython> interpPythonPtr;
#endif
  struct CommonData {
    boost::shared_ptr<VectorDouble> approxVals;
    boost::shared_ptr<VectorDouble> uAtPtsPtr;
    SmartPetscObj<Vec> petscVec;
    SmartPetscObj<Vec> L2Vec;
    SmartPetscObj<Vec> resVec;
 
    enum VecElements { VALUES_INTEG = 0, LAST_ELEMENT };
  };
 
  boost::shared_ptr<CommonData> commonDataPtr;
 
  struct OpError;
 
  struct OpCameraInteg : public BoundaryEleOp {
    boost::shared_ptr<CommonData> commonDataPtr;
    OpCameraInteg(boost::shared_ptr<CommonData> common_data_ptr)
        : BoundaryEleOp("PHOTON_FLUENCE_RATE", OPROW),
          commonDataPtr(common_data_ptr) {
      std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
      doEntities[MBTRI] = doEntities[MBQUAD] = true;
    }
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data);
  };
 
  struct OpGetScalarFieldGradientValuesOnSkin : public BoundaryEleOp {
 
    boost::shared_ptr<VolSideFe> sideOpFe;
 
    OpGetScalarFieldGradientValuesOnSkin(boost::shared_ptr<VolSideFe> side_fe)
        : BoundaryEleOp("PHOTON_FLUENCE_RATE", OPROW), sideOpFe(side_fe) {}
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data) {
      MoFEMFunctionBegin;
      if (type != MBVERTEX)
        MoFEMFunctionReturnHot(0);
      CHKERR loopSideVolumes("dFE", *sideOpFe);
      MoFEMFunctionReturn(0);
    }
  };
 
  struct Monitor : public FEMethod {
 
    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,
            MoFEM::Interface &m_field)
        : dM(dm), postProc(post_proc), skinPostProc(skin_post_proc),
          skinPostProcInteg(skin_post_proc_integ),
          commonDataPtr(common_data_ptr), mField(m_field) {}
 
    MoFEMErrorCode preProcess() {
      MoFEMFunctionBegin;
      MoFEMFunctionReturn(0);
    }
    MoFEMErrorCode operator()() {
      MoFEMFunctionBegin;
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode postProcess() {
      MoFEMFunctionBegin;
 
      auto vector_update = [&](auto vec) {
        MoFEMFunctionBegin;
 
        CHKERR VecZeroEntries(vec);
        CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR DMoFEMLoopFiniteElements(dM, "CAMERA_FE", skinPostProcInteg);
        CHKERR VecAssemblyBegin(vec);
        CHKERR VecAssemblyEnd(vec);
        CHKERR VecGhostUpdateBegin(vec, ADD_VALUES, SCATTER_REVERSE);
        CHKERR VecGhostUpdateEnd(vec, ADD_VALUES, SCATTER_REVERSE);
        CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
        MoFEMFunctionReturn(0);
      };
 
      if (!testing) {
        vector_update(commonDataPtr->petscVec);
        const double *array;
        CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
        MOFEM_LOG("PHOTON", Sev::inform)
            << "Fluence rate integral: " << array[0];
        CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
      } else {
        MOFEM_LOG("PHOTON", Sev::inform) << "Testing";
        CHKERR VecZeroEntries(commonDataPtr->resVec);
        vector_update(commonDataPtr->L2Vec);
        CHKERR VecAssemblyBegin(commonDataPtr->resVec);
        CHKERR VecAssemblyEnd(commonDataPtr->resVec);
        double nrm2;
        CHKERR VecNorm(commonDataPtr->resVec, NORM_2, &nrm2);
        const double *array2;
        CHKERR VecGetArrayRead(commonDataPtr->L2Vec, &array2);
        MOFEM_LOG("PHOTON", Sev::inform)
            << "Error " << array2[0] << " N of Entities " << nrm2;
 
        constexpr double eps = 0.02;
        if (array2[0] / nrm2 > eps)
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Not converged solution");
        CHKERR VecRestoreArrayRead(commonDataPtr->L2Vec, &array2);
      }
 
      if (ts_step % save_every_nth_step == 0) {
        if (output_volume) {
          CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
          CHKERR postProc->writeFile("out_volume_" +
                                     boost::lexical_cast<std::string>(ts_step) +
                                     ".h5m");
        }
        if (output_camera && skinPostProc) {
          CHKERR DMoFEMLoopFiniteElements(dM, "CAMERA_FE", skinPostProc);
          CHKERR skinPostProc->writeFile(
              "out_camera_" + boost::lexical_cast<std::string>(ts_step) +
              ".h5m");
        }
      }
      MoFEMFunctionReturn(0);
    }
 
  private:
    SmartPetscObj<DM> dM;
    boost::shared_ptr<PostProcEle> postProc;
    boost::shared_ptr<PostProcFaceEle> skinPostProc;
    boost::shared_ptr<BoundaryEle> skinPostProcInteg;
    boost::shared_ptr<CommonData> commonDataPtr;
    MoFEM::Interface &mField;
  };
};
 
struct PhotonDiffusion::OpError : public BoundaryEleOp {
  boost::shared_ptr<CommonData> commonDataPtr;
 
  static inline double sourceFunction(const double x, const double y,
                                      const double z, double time) {
    return SourceFunction::sourceFunctionEval(
        x, y, z, beam_radius, beam_centre_x, beam_centre_y, slab_thickness,
        mu_a, mu_sp, flux_magnitude, time, v, D);
  };
 
  OpError(boost::shared_ptr<CommonData> &common_data_ptr)
      : BoundaryEleOp("PHOTON_FLUENCE_RATE", OPROW),
        commonDataPtr(common_data_ptr) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
    const int nb_integration_pts = getGaussPts().size2();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_val = getFTensor0FromVec(*(commonDataPtr->uAtPtsPtr));
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    double nf = 1;
 
    FTensor::Index<'i', 3> i;
    const double volume = getMeasure();
 
    double error = 0;
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      const double alpha = t_w * volume;
      double analytical =
          sourceFunction(t_coords(0), t_coords(1), slab_thickness / 2, 0.01);
      double diff = t_val - analytical;
 
      error += alpha * pow(diff, 2) / std::abs(t_val);
 
      ++t_w;
      ++t_val;
      ++t_coords;
    }
 
    const int index = 0;
 
    CHKERR VecSetValue(commonDataPtr->L2Vec, index, error / nb_integration_pts,
                       ADD_VALUES);
    CHKERR VecSetValue(commonDataPtr->resVec, index, nf, ADD_VALUES);
 
    MoFEMFunctionReturn(0);
  }
};
 
#ifdef ENABLE_PYTHON_BINDING
struct PhotonDiffusion::InterpPython {
  InterpPython() = default;
  virtual ~InterpPython() = default;
 
  MoFEMErrorCode InterpInit(const std::string py_file);
  MoFEMErrorCode evalInterp(const std::string sens_image,
                            np::ndarray gauss_coords_x,
                            np::ndarray gauss_coords_y, double cam_len_x,
                            double cam_len_y, np::ndarray &sens_vals);
 
  template <typename T>
  inline std::vector<T>
  py_list_to_std_vector(const boost::python::object &iterable) {
    return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
                          boost::python::stl_input_iterator<T>());
  }
 
private:
  bp::object mainNamespace;
  bp::object InterpFun;
};
 
static boost::weak_ptr<PhotonDiffusion::InterpPython> interpPythonWeakPtr;
#endif
 
#ifdef ENABLE_PYTHON_BINDING
MoFEMErrorCode
PhotonDiffusion::InterpPython::InterpInit(const std::string py_file) {
  MoFEMFunctionBegin;
 
  try {
    bp::object main_module = bp::import("__main__");
    mainNamespace = main_module.attr("__dict__");
 
    bp::object ignored = bp::exec_file(py_file.c_str(), mainNamespace);
 
    InterpFun = mainNamespace
        ["py_Interpolate"]; // this is the function that will get called in the python script
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhotonDiffusion::InterpPython::evalInterp(
    const std::string sens_image, np::ndarray gauss_coords_x,
    np::ndarray gauss_coords_y, double cam_len_x, double cam_len_y,
    np::ndarray &sens_vals) {
  MoFEMFunctionBegin;
  try {
    sens_vals = bp::extract<np::ndarray>(InterpFun(
        sens_image, gauss_coords_x, gauss_coords_y, cam_len_x, cam_len_y));
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
 
  MoFEMFunctionReturn(0);
}
 
inline np::ndarray convert_to_numpy(VectorDouble &data, int nb_gauss_pts,
                                    int id) {
  auto dtype = np::dtype::get_builtin<double>();
  auto size = bp::make_tuple(nb_gauss_pts);
  auto stride = bp::make_tuple(3 * sizeof(double));
  return (np::from_data(&data[id], dtype, size, stride, bp::object()));
}
#endif
 
inline 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) {
#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
}
 
PhotonDiffusion::PhotonDiffusion(MoFEM::Interface &m_field) : mField(m_field) {}
 
MoFEMErrorCode PhotonDiffusion::readMesh() {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
  CHKERR mField.getInterface(simple);
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhotonDiffusion::createCommonData() {
  MoFEMFunctionBegin;
  commonDataPtr = boost::make_shared<CommonData>();
 
  auto *simple = mField.getInterface<Simple>();
  commonDataPtr->resVec = createDMVector(simple->getDM());
  commonDataPtr->uAtPtsPtr = boost::make_shared<VectorDouble>();
 
  PetscInt ghosts[1] = {0};
  if (!mField.get_comm_rank()) {
    commonDataPtr->petscVec =
        createGhostVector(mField.get_comm(), 1, 1, 0, ghosts);
    commonDataPtr->L2Vec =
        createGhostVector(mField.get_comm(), 1, 1, 0, ghosts);
  } else {
    commonDataPtr->petscVec =
        createGhostVector(mField.get_comm(), 0, 1, 1, ghosts);
    commonDataPtr->L2Vec =
        createGhostVector(mField.get_comm(), 0, 1, 1, ghosts);
  }
  commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhotonDiffusion::setupProblem() {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
  CHKERR simple->addDomainField("PHOTON_FLUENCE_RATE", H1,
                                AINSWORTH_LEGENDRE_BASE, 1);
  CHKERR simple->addBoundaryField("PHOTON_FLUENCE_RATE", H1,
                                  AINSWORTH_LEGENDRE_BASE, 1);
 
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-initial_file",
                               init_data_file_name, 255, PETSC_NULLPTR);
 
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-from_initial", &from_initial,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-output_volume",
                             &output_volume, PETSC_NULLPTR);
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-output_camera",
                             &output_camera, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-mu_a", &mu_a,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-mu_sp", &mu_sp,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-coef_A", &A, PETSC_NULLPTR);
 
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-slab_thickness",
                               &slab_thickness, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-beam_radius", &beam_radius,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-beam_centre_x",
                               &beam_centre_x, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-beam_centre_y",
                               &beam_centre_y, PETSC_NULLPTR);
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-testing", &testing,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cam_len_x", &cam_len_x,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cam_len_y", &cam_len_y,
                               PETSC_NULLPTR);
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-save_step",
                            &save_every_nth_step, PETSC_NULLPTR);
 
#ifdef ENABLE_PYTHON_BINDING
  if (enable_python) {
    auto file_exists = [](std::string myfile) {
      std::ifstream file(myfile.c_str());
      if (file) {
        return true;
      }
      return false;
    };
    CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-interp_file",
                                 interp_file_name, 255, PETSC_NULLPTR);
    CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-sens_image",
                                 interp_image_name, 255, PETSC_NULLPTR);
  #ifdef ENABLE_PYTHON_BINDING
    if (file_exists(interp_file_name)) {
      MOFEM_LOG("PHOTON", Sev::inform) << interp_file_name << " file found";
      interpPythonPtr = boost::make_shared<PhotonDiffusion::InterpPython>();
      CHKERR interpPythonPtr->InterpInit(interp_file_name);
      interpPythonWeakPtr = interpPythonPtr;
    } else {
      MOFEM_LOG("PHOTON", Sev::warning)
          << interp_file_name << " file NOT found";
    }
  }
  #endif
#endif
 
  h = 0.5 / A;
  D = 1. / (3. * (mu_a + mu_sp));
 
  MOFEM_LOG("PHOTON", Sev::inform) << "Refractive index: " << n;
  MOFEM_LOG("PHOTON", Sev::inform) << "Speed of light (cm/ns): " << c;
  MOFEM_LOG("PHOTON", Sev::inform) << "Phase velocity in medium (cm/ns): " << v;
  MOFEM_LOG("PHOTON", Sev::inform) << "Inverse velocity : " << inv_v;
  MOFEM_LOG("PHOTON", Sev::inform)
      << "Absorption coefficient (cm^-1): " << mu_a;
  MOFEM_LOG("PHOTON", Sev::inform)
      << "Scattering coefficient (cm^-1): " << mu_sp;
  MOFEM_LOG("PHOTON", Sev::inform) << "Diffusion coefficient D : " << D;
  MOFEM_LOG("PHOTON", Sev::inform) << "Coefficient A : " << A;
  MOFEM_LOG("PHOTON", Sev::inform) << "Coefficient h : " << h;
 
  MOFEM_LOG("PHOTON", Sev::inform) << "Approximation order: " << order;
  MOFEM_LOG("PHOTON", Sev::inform) << "Save step: " << save_every_nth_step;
 
  CHKERR simple->setFieldOrder("PHOTON_FLUENCE_RATE", order);
 
  auto set_camera_skin_fe = [&]() {
    MoFEMFunctionBegin;
 
    Range camera_surface;
    const std::string block_name = "CAM";
    bool add_fe = false;
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, bit)) {
      if (bit->getName().compare(0, block_name.size(), block_name) == 0) {
        MOFEM_LOG("PHOTON", Sev::inform) << "Found CAM block";
        CHKERR mField.get_moab().get_entities_by_dimension(
            bit->getMeshset(), 2, camera_surface, true);
        add_fe = true;
      }
    }
 
    MOFEM_LOG("PHOTON", Sev::noisy) << "CAM block entities:\n"
                                    << camera_surface;
 
    if (add_fe) {
      CHKERR mField.add_finite_element("CAMERA_FE");
      CHKERR mField.modify_finite_element_add_field_data("CAMERA_FE",
                                                         "PHOTON_FLUENCE_RATE");
      CHKERR mField.add_ents_to_finite_element_by_dim(camera_surface, 2,
                                                      "CAMERA_FE");
    }
    MoFEMFunctionReturn(0);
  };
 
  auto my_simple_set_up = [&]() {
    MoFEMFunctionBegin;
    CHKERR simple->defineFiniteElements();
    CHKERR simple->defineProblem(PETSC_TRUE);
    CHKERR simple->buildFields();
    CHKERR simple->buildFiniteElements();
 
    if (mField.check_finite_element("CAMERA_FE")) {
      CHKERR mField.build_finite_elements("CAMERA_FE");
      CHKERR DMMoFEMAddElement(simple->getDM(), "CAMERA_FE");
    }
 
    CHKERR simple->buildProblem();
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_camera_skin_fe();
  CHKERR my_simple_set_up();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhotonDiffusion::initialCondition() {
  MoFEMFunctionBegin;
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhotonDiffusion::boundaryCondition() {
  MoFEMFunctionBegin;
  auto bc_mng = mField.getInterface<BcManager>();
  auto *simple = mField.getInterface<Simple>();
  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "EXT",
                                        "PHOTON_FLUENCE_RATE", 0, 0, false);
 
  // Get boundary faces marked in block named "INT"
  Range boundary_faces;
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
    std::string entity_name = it->getName();
    if (entity_name.compare(0, 3, "INT") == 0) {
      CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2,
                                                 boundary_faces, true);
    }
  }
 
  // If the previous search has not returned anything, assume the boundary has been created as a sideset.
  if (boundary_faces.empty()) {
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, SIDESET, it)) {
      std::string entity_name = it->getName();
      CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2,
                                                 boundary_faces, true);
    }
  }
 
  // Get boundary edges in "INT"
  Range boundary_ents;
  CHKERR mField.get_moab().get_adjacencies(
      boundary_faces, 1, false, boundary_ents, moab::Interface::UNION);
  // Add vertices to boundary entities
  Range boundary_verts;
  CHKERR mField.get_moab().get_adjacencies(
      boundary_faces, 0, false, boundary_verts, moab::Interface::UNION);
 
  boundary_faces.merge(boundary_verts);
  boundary_faces.merge(boundary_ents);
 
  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
      boundary_faces);
 
  // Remove DOFs as homogeneous boundary condition is used
  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
      simple->getProblemName(), "PHOTON_FLUENCE_RATE", boundary_faces);
 
  MoFEMFunctionReturn(0);
}
 
//! [assembleSystem]
MoFEMErrorCode PhotonDiffusion::assembleSystem() {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto integration_rule = [](int o_row, int o_col, int approx_order) {
    return 2 * approx_order;
  };
 
  auto bc_mng = mField.getInterface<BcManager>();
  auto &bc_map = bc_mng->getBcMapByBlockName();
 
  auto set_domain = [&]() {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        pipeline_mng->getOpDomainLhsPipeline(), {H1});
 
    pipeline_mng->getOpDomainLhsPipeline().push_back(new OpDomainGradGrad(
        "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
        [](double, double, double) -> double { return D; }));
 
    auto get_mass_coefficient = [&](const double, const double, const double) {
      return inv_v * domainLhsFEPtr->ts_a + mu_a;
    };
    pipeline_mng->getOpDomainLhsPipeline().push_back(new OpDomainMass(
        "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE", get_mass_coefficient));
 
    auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
    auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
    auto dot_u_at_gauss_pts = boost::make_shared<VectorDouble>();
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpCalculateScalarFieldGradient<SPACE_DIM>("PHOTON_FLUENCE_RATE",
                                                      grad_u_at_gauss_pts));
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
                                         u_at_gauss_pts));
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpCalculateScalarFieldValuesDot("PHOTON_FLUENCE_RATE",
                                            dot_u_at_gauss_pts));
    pipeline_mng->getOpDomainRhsPipeline().push_back(new OpDomainGradTimesVec(
        "PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts,
        [](double, double, double) -> double { return D; }));
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpDomainTimesScalarField(
            "PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts,
            [](const double, const double, const double) { return inv_v; }));
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpDomainTimesScalarField(
            "PHOTON_FLUENCE_RATE", u_at_gauss_pts,
            [](const double, const double, const double) { return mu_a; }));
 
    CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
    CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
    MoFEMFunctionReturn(0);
  };
 
  auto set_boundary = [&]() {
    MoFEMFunctionBegin;
 
    pipeline_mng->getOpBoundaryLhsPipeline().push_back(
        new OpSetHOWeightsOnFace());
    pipeline_mng->getOpBoundaryRhsPipeline().push_back(
        new OpSetHOWeightsOnFace());
 
    auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
    pipeline_mng->getOpBoundaryRhsPipeline().push_back(
        new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
                                         u_at_gauss_pts));
    for (auto b : bc_map) {
      if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
        pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpBoundaryMass(
            "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
 
            [](const double, const double, const double) { return h; },
 
            b.second->getBcEntsPtr()));
 
        pipeline_mng->getOpBoundaryRhsPipeline().push_back(
            new OpBoundaryTimeScalarField(
                "PHOTON_FLUENCE_RATE", u_at_gauss_pts,
 
                [](const double, const double, const double) { return h; },
 
                b.second->getBcEntsPtr()));
      }
    }
 
    CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
    CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_domain();
  CHKERR set_boundary();
 
  domainLhsFEPtr = pipeline_mng->getDomainLhsFE();
  boundaryLhsFEPtr = pipeline_mng->getBoundaryLhsFE();
  boundaryRhsFEPtr = pipeline_mng->getBoundaryRhsFE();
 
  MoFEMFunctionReturn(0);
}
//! [assembleSystem]
 
MoFEMErrorCode PhotonDiffusion::solveSystem() {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  //! [postprocess]
  auto create_post_process_element = [&]() {
    auto post_froc_fe = boost::make_shared<PostProcEle>(mField);
    auto u_ptr = boost::make_shared<VectorDouble>();
    auto grad_ptr = boost::make_shared<MatrixDouble>();
    post_froc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE", u_ptr));
    post_froc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldGradient<SPACE_DIM>("PHOTON_FLUENCE_RATE",
                                                      grad_ptr));
    post_froc_fe->getOpPtrVector().push_back(new OpPPMap(
        post_froc_fe->getPostProcMesh(), post_froc_fe->getMapGaussPts(),
        {{"PHOTON_FLUENCE_RATE", u_ptr}},
        {{"GRAD_PHOTON_FLUENCE_RATE", grad_ptr}}, {}, {}));
    return post_froc_fe;
  };
 
  auto create_post_process_camera_element = [&]() {
    if (mField.check_finite_element("CAMERA_FE")) {
      auto post_proc_skin = boost::make_shared<PostProcFaceEle>(mField);
 
      auto u_ptr = boost::make_shared<VectorDouble>();
      auto grad_ptr = boost::make_shared<MatrixDouble>();
 
      auto op_loop_side = new OpLoopSide<FaceElementForcesAndSourcesCoreOnSide>(
          mField, simple->getDomainFEName(), SPACE_DIM);
 
      // push operators to side element
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_loop_side->getOpPtrVector(), {H1});
      op_loop_side->getOpPtrVector().push_back(
          new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE", u_ptr));
      op_loop_side->getOpPtrVector().push_back(
          new OpCalculateScalarFieldGradient<SPACE_DIM>("PHOTON_FLUENCE_RATE",
                                                        grad_ptr));
      // push op to boundary element
      post_proc_skin->getOpPtrVector().push_back(op_loop_side);
 
      post_proc_skin->getOpPtrVector().push_back(new OpPPMap(
          post_proc_skin->getPostProcMesh(), post_proc_skin->getMapGaussPts(),
          {{"PHOTON_FLUENCE_RATE", u_ptr}},
          {{"GRAD_PHOTON_FLUENCE_RATE", grad_ptr}}, {}, {}));
 
      return post_proc_skin;
    } else {
      return boost::shared_ptr<PostProcFaceEle>();
    }
  };
 
  auto create_post_process_integ_camera_element = [&]() {
    if (mField.check_finite_element("CAMERA_FE")) {
      auto post_proc_integ_skin = boost::make_shared<BoundaryEle>(mField);
 
      if (!testing) {
        MOFEM_LOG("PHOTON", Sev::inform)
            << "Creating post process integ camera element";
        post_proc_integ_skin->getOpPtrVector().push_back(
            new OpSetHOWeightsOnFace());
        post_proc_integ_skin->getOpPtrVector().push_back(
            new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
                                             commonDataPtr->approxVals));
        post_proc_integ_skin->getOpPtrVector().push_back(
            new OpCameraInteg(commonDataPtr));
      } else {
        MOFEM_LOG("PHOTON", Sev::inform)
            << "Creating testing camera error element";
        post_proc_integ_skin->getOpPtrVector().push_back(
            new OpSetHOWeightsOnFace());
        post_proc_integ_skin->getRuleHook = [](int, int, int approx_order) {
          return 2 * approx_order;
        };
        CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
            post_proc_integ_skin->getOpPtrVector(), {NOSPACE});
        post_proc_integ_skin->getOpPtrVector().push_back(
            new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
                                             commonDataPtr->uAtPtsPtr));
        post_proc_integ_skin->getOpPtrVector().push_back(
            new OpError(commonDataPtr));
      }
      return post_proc_integ_skin;
    } else {
      return boost::shared_ptr<BoundaryEle>();
    }
  };
  //! [postprocess]
 
  //! [timesolver]
  auto set_time_monitor = [&](auto dm, auto solver) {
    MoFEMFunctionBegin;
    boost::shared_ptr<Monitor> monitor_ptr(new Monitor(
        dm, create_post_process_element(), create_post_process_camera_element(),
        create_post_process_integ_camera_element(), commonDataPtr, mField));
    boost::shared_ptr<ForcesAndSourcesCore> null;
    CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
                               monitor_ptr, null, null);
    MoFEMFunctionReturn(0);
  };
 
  auto dm = simple->getDM();
  auto X = createDMVector(dm);
 
  if (from_initial) {
 
    MOFEM_LOG("PHOTON", Sev::inform) << "reading vector in binary from file "
                                     << init_data_file_name << " ...";
    PetscViewer viewer;
    PetscViewerBinaryOpen(PETSC_COMM_WORLD, init_data_file_name, FILE_MODE_READ,
                          &viewer);
    VecLoad(X, viewer);
 
    CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_REVERSE);
  }
 
  auto solver = pipeline_mng->createTSIM();
 
  CHKERR TSSetSolution(solver, X);
  CHKERR set_time_monitor(dm, solver);
  CHKERR TSSetSolution(solver, X);
  CHKERR TSSetFromOptions(solver);
  auto B = createDMMatrix(dm);
  CHKERR TSSetIJacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
  CHKERR TSSetUp(solver);
  CHKERR TSSolve(solver, NULL);
 
  CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_REVERSE);
 
  MoFEMFunctionReturn(0);
}
//! [timesolver]
//! [solveSystem]
 
MoFEMErrorCode PhotonDiffusion::outputResults() {
  MoFEMFunctionBegin;
 
  // Processes to set output results are integrated in solveSystem()
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhotonDiffusion::runProgram() {
  MoFEMFunctionBegin;
 
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR createCommonData();
  CHKERR initialCondition();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR outputResults();
 
  MoFEMFunctionReturn(0);
}
 
//! [Integral_calc]
MoFEMErrorCode
PhotonDiffusion::OpCameraInteg::doWork(int side, EntityType type,
                                       EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
  const int nb_integration_pts = getGaussPts().size2();
  const double area = getMeasure();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
 
  double values_integ = 0;
 
#ifdef ENABLE_PYTHON_BINDING
  VectorDouble sens_vals_vec;
  if (enable_python) {
    MatrixDouble ref_coords = getCoordsAtGaussPts();
    sens_vals_vec =
        interp_function(interp_image_name, ref_coords, nb_integration_pts,
                        cam_len_x, cam_len_y, "(.*)INTERPOLATION(.*)");
  }
#endif
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
    double sens = 1;
#ifdef ENABLE_PYTHON_BINDING
    if (enable_python) {
      sens = sens_vals_vec(gg);
    }
#endif
 
    const double alpha = t_w * area * sens;
 
    values_integ += alpha * t_val;
 
    ++t_w;
    ++t_val;
  }
 
  constexpr std::array<int, 1> indices = {CommonData::VALUES_INTEG};
  std::array<double, 1> values;
  values[0] = values_integ;
  CHKERR VecSetValues(commonDataPtr->petscVec, 1, indices.data(), values.data(),
                      ADD_VALUES);
  MoFEMFunctionReturn(0);
}
//! [Integral_calc]
 
int main(int argc, char *argv[]) {
  // 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;
}