PhotonDiffusion Struct Reference

Collaboration diagram for PhotonDiffusion:

Classes
struct	CommonData
struct	Monitor
struct	OpCameraInteg
struct	OpGetScalarFieldGradientValuesOnSkin

Public Member Functions
	PhotonDiffusion (MoFEM::Interface &m_field)
MoFEMErrorCode	runProgram ()
	PhotonDiffusion (MoFEM::Interface &m_field)
MoFEMErrorCode	runProgram ()

Static Public Member Functions
static double	sourceFunction (const double x, const double y, const double z)
	Pulse is infinitely short. More...

Private Member Functions
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	setIntegrationRules ()
MoFEMErrorCode	initialCondition ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	outputResults ()
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	setIntegrationRules ()
MoFEMErrorCode	initialCondition ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()

Private Attributes
MoFEM::Interface &	mField
boost::shared_ptr< std::vector< unsigned char > >	boundaryMarker
boost::shared_ptr< FEMethod >	domainLhsFEPtr
boost::shared_ptr< FEMethod >	boundaryLhsFEPtr
boost::shared_ptr< FEMethod >	boundaryRhsFEPtr
boost::shared_ptr< CommonData >	commonDataPtr

Detailed Description

Examples

photon_diffusion.cpp.

Definition at line 74 of file initial_diffusion.cpp.

Constructor & Destructor Documentation

PhotonDiffusion() [1/2]

PhotonDiffusion::PhotonDiffusion

(

MoFEM::Interface &

m_field

)

Examples

photon_diffusion.cpp.

Definition at line 138 of file initial_diffusion.cpp.

: mField(m_field) {}

PhotonDiffusion() [2/2]

PhotonDiffusion::PhotonDiffusion

(

MoFEM::Interface &

m_field

)

Member Function Documentation

assembleSystem() [1/2]

MoFEMErrorCode PhotonDiffusion::assembleSystem ( )

private

assembleSystem() [2/2]

MoFEMErrorCode PhotonDiffusion::assembleSystem ( )

private

Examples

photon_diffusion.cpp.

Definition at line 284 of file initial_diffusion.cpp.

                                               {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  auto &bc_map = bc_mng->getBcMapByBlockName();
 
  auto add_domain_base_ops = [&](auto &pipeline) {
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto det_ptr = boost::make_shared<VectorDouble>();
    pipeline.push_back(new OpCalculateHOJacVolume(jac_ptr));
    pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, nullptr));
    pipeline.push_back(new OpSetHOWeights(det_ptr));
  };
 
  auto add_domain_lhs_ops = [&](auto &pipeline) {
    pipeline.push_back(new OpDomainMass(
        "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
        [](const double, const double, const double) { return 1; }));
  };
 
  auto add_domain_rhs_ops = [&](auto &pipeline) {
    pipeline.push_back(
        new OpDomainSource("PHOTON_FLUENCE_RATE", sourceFunction));
  };
 
  auto pipeline_mng = mField.getInterface<PipelineManager>();
  add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
  add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
  add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
  add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
 
  MoFEMFunctionReturn(0);
}

boundaryCondition() [1/2]

MoFEMErrorCode PhotonDiffusion::boundaryCondition ( )

private

boundaryCondition() [2/2]

MoFEMErrorCode PhotonDiffusion::boundaryCondition ( )

private

Examples

photon_diffusion.cpp.

Definition at line 257 of file initial_diffusion.cpp.

                                                  {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
 
  // Get boundary edges marked in block named "BOUNDARY_CONDITION"
  Range boundary_ents;
  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(), 1,
                                                 boundary_ents, true);
    }
  }
  // Add vertices to boundary entities
  Range boundary_verts;
  CHKERR mField.get_moab().get_connectivity(boundary_ents, boundary_verts,
                                            true);
  boundary_ents.merge(boundary_verts);
 
  // Remove DOFs as homogeneous boundary condition is used
  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
      simple->getProblemName(), "PHOTON_FLUENCE_RATE", boundary_ents);
 
  MoFEMFunctionReturn(0);
}

checkResults()

MoFEMErrorCode PhotonDiffusion::checkResults ( )

private

createCommonData()

MoFEMErrorCode PhotonDiffusion::createCommonData ( )

private

Examples

photon_diffusion.cpp.

Definition at line 223 of file photon_diffusion.cpp.

                                                 {
  MoFEMFunctionBegin;
  commonDataPtr = boost::make_shared<CommonData>();
  PetscInt ghosts[1] = {0};
  if (!mField.get_comm_rank())
    commonDataPtr->petscVec =
        createGhostVector(mField.get_comm(), 1, 1, 0, ghosts);
  else
    commonDataPtr->petscVec =
        createGhostVector(mField.get_comm(), 0, 1, 1, ghosts);
  commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
  MoFEMFunctionReturn(0);
}

initialCondition() [1/2]

MoFEMErrorCode PhotonDiffusion::initialCondition ( )

private

initialCondition() [2/2]

MoFEMErrorCode PhotonDiffusion::initialCondition ( )

private

Examples

photon_diffusion.cpp.

Definition at line 252 of file initial_diffusion.cpp.

                                                 {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}

outputResults() [1/2]

MoFEMErrorCode PhotonDiffusion::outputResults ( )

private

outputResults() [2/2]

MoFEMErrorCode PhotonDiffusion::outputResults ( )

private

Examples

photon_diffusion.cpp.

Definition at line 349 of file initial_diffusion.cpp.

                                              {
  MoFEMFunctionBegin;
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
  pipeline_mng->getDomainLhsFE().reset();
  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
 
  auto u_ptr = boost::make_shared<VectorDouble>();
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues("U", u_ptr));
 
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
  post_proc_fe->getOpPtrVector().push_back(
 
      new OpPPMap(
 
          post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(), 
          
          {{"PHOTON_FLUENCE_RATE", u_ptr}},
 
          {},
 
          {},
 
          {})
 
  );     
 
  pipeline_mng->getDomainRhsFE() = post_proc_fe;
  CHKERR pipeline_mng->loopFiniteElements();
  CHKERR post_proc_fe->writeFile("out_initial.h5m");
  MoFEMFunctionReturn(0);
}

readMesh() [1/2]

MoFEMErrorCode PhotonDiffusion::readMesh ( )

private

readMesh() [2/2]

MoFEMErrorCode PhotonDiffusion::readMesh ( )

private

Examples

photon_diffusion.cpp.

Definition at line 140 of file initial_diffusion.cpp.

                                         {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
  CHKERR mField.getInterface(simple);
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
 
  MoFEMFunctionReturn(0);
}

runProgram() [1/2]

MoFEMErrorCode PhotonDiffusion::runProgram

(

)

runProgram() [2/2]

MoFEMErrorCode PhotonDiffusion::runProgram

(

)

Examples

photon_diffusion.cpp.

Definition at line 383 of file initial_diffusion.cpp.

                                           {
  MoFEMFunctionBegin;
 
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR setIntegrationRules();
  CHKERR initialCondition();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  if (output_volume)
    CHKERR outputResults();
 
  MoFEMFunctionReturn(0);
}

setIntegrationRules() [1/2]

MoFEMErrorCode PhotonDiffusion::setIntegrationRules ( )

private

setIntegrationRules() [2/2]

MoFEMErrorCode PhotonDiffusion::setIntegrationRules ( )

private

Examples

photon_diffusion.cpp.

Definition at line 236 of file initial_diffusion.cpp.

                                                    {
  MoFEMFunctionBegin;
 
  auto integration_rule = [](int o_row, int o_col, int approx_order) {
    return 2 * approx_order;
  };
 
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
 
  MoFEMFunctionReturn(0);
}

setupProblem() [1/2]

MoFEMErrorCode PhotonDiffusion::setupProblem ( )

private

setupProblem() [2/2]

MoFEMErrorCode PhotonDiffusion::setupProblem ( )

private

Examples

photon_diffusion.cpp.

Definition at line 151 of file initial_diffusion.cpp.

                                             {
  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 PetscOptionsGetScalar(PETSC_NULL, "", "-flux_magnitude",
                               &flux_magnitude, PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-slab_thickness",
                               &slab_thickness, PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-beam_radius", &beam_radius,
                               PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-beam_centre_x", &beam_centre_x,
                               PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-beam_centre_y", &beam_centre_y,
                               PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-mu_a", &mu_a, PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-mu_sp", &mu_sp, PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-initial_time", &initial_time,
                               PETSC_NULL);
 
  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-output_file", out_file_name,
                               255, PETSC_NULL);
  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-output_volume", &output_volume,
                             PETSC_NULL);
 
  D = 1. / (3. * (mu_a + mu_sp));                    
 
  MOFEM_LOG("INITIAL", Sev::inform) << "Refractive index: " << n;
  MOFEM_LOG("INITIAL", Sev::inform) << "Speed of light (cm/ns): " << c;
  MOFEM_LOG("INITIAL", Sev::inform)
      << "Phase velocity in medium (cm/ns): " << v;
  MOFEM_LOG("INITIAL", Sev::inform)
      << "Absorption coefficient (cm^-1): " << mu_a;
  MOFEM_LOG("INITIAL", Sev::inform)
      << "Scattering coefficient (cm^-1): " << mu_sp;
  MOFEM_LOG("INITIAL", Sev::inform) << "Diffusion coefficient D : " << D;
  MOFEM_LOG("INITIAL", Sev::inform) << "Impulse magnitude: " << flux_magnitude;
  MOFEM_LOG("INITIAL", Sev::inform) << "Compute time (ns): " << initial_time;
  MOFEM_LOG("INITIAL", Sev::inform) << "Slab thickness: " << slab_thickness;
 
  int order = 2;
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
 
  MOFEM_LOG("INITIAL", Sev::inform) << "Approximation order: " << order;
 
  CHKERR simple->setFieldOrder("PHOTON_FLUENCE_RATE", order);
 
  // if (numHoLevels > 0) {
 
  //   Range ho_ents;
  //   for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
  //     if (it->getName().compare(0, 3, "CAM") == 0) {
  //       CHKERR mField.get_moab().get_entities_by_dimension(it->getMeshset(), 2,
  //                                                          ho_ents, true);
  //     }
  //   }
 
  //   EntityHandle meshset;
  //   CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
  //   CHKERR mField.get_moab().add_entities(meshset, ho_ents);
  //   std::string field_name;
  //   field_name = "out_test_" +
  //                boost::lexical_cast<std::string>(mField.get_comm_rank()) +
  //                ".vtk";
  //   CHKERR mField.get_moab().write_file(field_name.c_str(), "VTK", "", &meshset,
  //                                       1);
  //   CHKERR mField.get_moab().delete_entities(&meshset, 1);
 
  //   CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ho_ents);
 
  //   CHKERR simple->setFieldOrder("PHOTON_FLUENCE_RATE", order + 1, &ho_ents);
 
  //   CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
  //       "PHOTON_FLUENCE_RATE");
  // }
 
  CHKERR simple->setUp();
 
  MoFEMFunctionReturn(0);
}

solveSystem() [1/2]

MoFEMErrorCode PhotonDiffusion::solveSystem ( )

private

solveSystem() [2/2]

MoFEMErrorCode PhotonDiffusion::solveSystem ( )

private

Examples

photon_diffusion.cpp.

Definition at line 318 of file initial_diffusion.cpp.

                                            {
  MoFEMFunctionBegin;
  auto *simple = mField.getInterface<Simple>();
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
  auto solver = pipeline_mng->createKSP();
 
  CHKERR KSPSetFromOptions(solver);
  // CHKERR KSPSetUp(solver);
 
  auto dm = simple->getDM();
  auto X = createDMVector(dm);
  auto F = vectorDuplicate(X);
 
  MOFEM_LOG("INITIAL", Sev::inform) << "Solver start";
  CHKERR KSPSolve(solver, F, X);
  CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_REVERSE);
 
  MOFEM_LOG("INITIAL", Sev::inform)
      << "writing vector in binary to " << out_file_name << " ...";
  PetscViewer viewer;
  PetscViewerBinaryOpen(PETSC_COMM_WORLD, out_file_name, FILE_MODE_WRITE,
                        &viewer);
  VecView(X, viewer);
  PetscViewerDestroy(&viewer);
 
  MOFEM_LOG("INITIAL", Sev::inform) << "Solver done";
  MoFEMFunctionReturn(0);
}

sourceFunction()

static double PhotonDiffusion::sourceFunction ( const double  x, const double  y, const double  z  )

inlinestatic

Pulse is infinitely short.

Note

It is good approximation of pulse in femtosecond scale. To make it longer one can apply third integral over time.

Note analysis in photon_diffusion is shifted in time by initial_time.

Parameters

x
y
z

Returns

double

Definition at line 94 of file initial_diffusion.cpp.

                                                                               {
    const double A = 4. * D * v * initial_time;
    const double T =
        (v / pow(M_PI * A, 3. / 2.)) * exp(-mu_a * v * initial_time);
 
    auto phi_pulse = [&](const double r_s, const double phi_s) {
      const double xs = r_s * cos(phi_s);
      const double ys = r_s * sin(phi_s);
      const double xp = x - xs - beam_centre_x;
      const double yp = y - ys - beam_centre_y;
      const double zp1 = z + slab_thickness / 2. - 1. / mu_sp;
      const double zp2 = z + slab_thickness / 2. + 1. / mu_sp;
      const double P1 = xp * xp + yp * yp + zp1 * zp1;
      const double P2 = xp * xp + yp * yp + zp2 * zp2;
      return r_s * (exp(-P1 / A) - exp(-P2 / A));
    };
 
    auto f = [&](const double r_s) {
      auto g = [&](const double phi_s) { return phi_pulse(r_s, phi_s); };
      return gauss_kronrod<double, 15>::integrate(
          g, 0, 2 * M_PI, 0, std::numeric_limits<float>::epsilon());
    };
 
    return T * flux_magnitude *
           gauss_kronrod<double, 15>::integrate(
               f, 0, beam_radius, 0, std::numeric_limits<float>::epsilon());
  };

Member Data Documentation

boundaryLhsFEPtr

boost::shared_ptr<FEMethod> PhotonDiffusion::boundaryLhsFEPtr

private

Examples

photon_diffusion.cpp.

Definition at line 99 of file photon_diffusion.cpp.

boundaryMarker

boost::shared_ptr<std::vector<unsigned char> > PhotonDiffusion::boundaryMarker

private

Definition at line 96 of file photon_diffusion.cpp.

boundaryRhsFEPtr

boost::shared_ptr<FEMethod> PhotonDiffusion::boundaryRhsFEPtr

private

Examples

photon_diffusion.cpp.

Definition at line 100 of file photon_diffusion.cpp.

commonDataPtr

boost::shared_ptr<CommonData> PhotonDiffusion::commonDataPtr

private

Examples

photon_diffusion.cpp.

Definition at line 112 of file photon_diffusion.cpp.

domainLhsFEPtr

boost::shared_ptr<FEMethod> PhotonDiffusion::domainLhsFEPtr

private

Examples

photon_diffusion.cpp.

Definition at line 98 of file photon_diffusion.cpp.

mField

MoFEM::Interface & PhotonDiffusion::mField

private

Examples

photon_diffusion.cpp.

Definition at line 135 of file initial_diffusion.cpp.

---

The documentation for this struct was generated from the following files:

• tutorials/scl-10/initial_diffusion.cpp
• tutorials/scl-10/photon_diffusion.cpp