PhotonDiffusion Struct Reference

Collaboration diagram for PhotonDiffusion:

Classes
struct	CommonData
struct	Monitor
struct	OpCameraInteg
struct	OpError
struct	OpGetScalarFieldGradientValuesOnSkin

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

Static Public Member Functions
static double	sourceFunction (const double x, const double y, const double z)

Private Member Functions
MoFEMErrorCode	readMesh ()
	[sourceFunction]
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	initialCondition ()
MoFEMErrorCode	boundaryCondition ()
	[interior_boundary]
MoFEMErrorCode	assembleSystem ()
	[interior_boundary]
MoFEMErrorCode	solveSystem ()
	[assembleSystem]
MoFEMErrorCode	checkResults ()
	[runProgram]
MoFEMErrorCode	outputResults ()
	[solveSystem]
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	setIntegrationRules ()
MoFEMErrorCode	initialCondition ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	outputResults ()

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

Detailed Description

Definition at line 67 of file initial_diffusion.cpp.

Constructor & Destructor Documentation

PhotonDiffusion() [1/2]

PhotonDiffusion::PhotonDiffusion

(

MoFEM::Interface &

m_field

)

Definition at line 108 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

[interior_boundary]

[assembleSystem]

Definition at line 334 of file initial_diffusion.cpp.

                                               {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto integration_rule = [](int o_row, int o_col, int approx_order) {
    return 2 * approx_order;
  };
 
  auto set_domain = [&]() {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        pipeline_mng->getOpDomainLhsPipeline(), {H1});
 
    pipeline_mng->getOpDomainLhsPipeline().push_back(new OpDomainMass(
        "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
        [](const double, const double, const double) { return 1; }));
 
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpDomainSource("PHOTON_FLUENCE_RATE", sourceFunction));
 
    CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
    CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
    MoFEMFunctionReturn(0);
  };
 
  auto set_boundary = [&]() {
    MoFEMFunctionBegin;
    CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
    CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_domain();
  CHKERR set_boundary();
 
  MoFEMFunctionReturn(0);
}

assembleSystem() [2/2]

MoFEMErrorCode PhotonDiffusion::assembleSystem ( )

private

boundaryCondition() [1/2]

MoFEMErrorCode PhotonDiffusion::boundaryCondition ( )

private

[interior_boundary]

Definition at line 282 of file initial_diffusion.cpp.

                                                  {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
 
  // 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 (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);
 
  EntityHandle meshset;
  CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
  CHKERR mField.get_moab().add_entities(meshset, boundary_faces);
 
  // Remove DOFs as homogeneous boundary condition is used
 
  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
      simple->getProblemName(), "PHOTON_FLUENCE_RATE", boundary_faces);
 
  MoFEMFunctionReturn(0);
}

boundaryCondition() [2/2]

MoFEMErrorCode PhotonDiffusion::boundaryCondition ( )

private

checkResults() [1/2]

MoFEMErrorCode PhotonDiffusion::checkResults ( )

private

[runProgram]

[Check results]

Definition at line 459 of file initial_diffusion.cpp.

                                             {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
  pipeline_mng->getDomainLhsFE().reset();
  pipeline_mng->getDomainRhsFE().reset();
  pipeline_mng->getOpDomainRhsPipeline().clear();
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
                                       commonDataPtr->uAtPtsPtr));
  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpError(commonDataPtr));
  CHKERR pipeline_mng->loopFiniteElements();
  CHKERR VecAssemblyBegin(commonDataPtr->L2Vec);
  CHKERR VecAssemblyEnd(commonDataPtr->L2Vec);
  CHKERR VecAssemblyBegin(commonDataPtr->resVec);
  CHKERR VecAssemblyEnd(commonDataPtr->resVec);
  double nrm2;
  CHKERR VecNorm(commonDataPtr->resVec, NORM_2, &nrm2);
  const double *array;
  CHKERR VecGetArrayRead(commonDataPtr->L2Vec, &array);
  if (mField.get_comm_rank() == 0)
    PetscPrintf(PETSC_COMM_SELF, "Error %6.4e Vec norm %6.4e\n",
                std::sqrt(array[0]), nrm2);
  CHKERR VecRestoreArrayRead(commonDataPtr->L2Vec, &array);
  constexpr double eps = 1e-8;
  if (nrm2 > eps)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Not converged solution");
  MoFEMFunctionReturn(0);
}

checkResults() [2/2]

MoFEMErrorCode PhotonDiffusion::checkResults ( )

private

createCommonData() [1/2]

MoFEMErrorCode PhotonDiffusion::createCommonData ( )

private

Definition at line 263 of file initial_diffusion.cpp.

                                                 {
  MoFEMFunctionBegin;
  auto *simple = mField.getInterface<Simple>();
 
  commonDataPtr = boost::make_shared<PhotonDiffusion::CommonData>();
  commonDataPtr->resVec = createDMVector(simple->getDM());
  commonDataPtr->L2Vec =
      createVectorMPI(mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
  commonDataPtr->uAtPtsPtr = boost::make_shared<VectorDouble>();
 
  MoFEMFunctionReturn(0);
}

createCommonData() [2/2]

MoFEMErrorCode PhotonDiffusion::createCommonData ( )

private

initialCondition() [1/2]

MoFEMErrorCode PhotonDiffusion::initialCondition ( )

private

Definition at line 276 of file initial_diffusion.cpp.

                                                 {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}

initialCondition() [2/2]

MoFEMErrorCode PhotonDiffusion::initialCondition ( )

private

outputResults() [1/2]

MoFEMErrorCode PhotonDiffusion::outputResults ( )

private

[solveSystem]

[timesolver] [solveSystem]

Definition at line 405 of file initial_diffusion.cpp.

                                              {
  MoFEMFunctionBegin;
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
 
  auto u_ptr = boost::make_shared<VectorDouble>();
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE", 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->getDomainPostProcFE() = post_proc_fe;
  CHKERR pipeline_mng->loopFiniteElementsPostProc();
  CHKERR post_proc_fe->writeFile("out_initial.h5m");
  MoFEMFunctionReturn(0);
}

outputResults() [2/2]

MoFEMErrorCode PhotonDiffusion::outputResults ( )

private

readMesh() [1/2]

MoFEMErrorCode PhotonDiffusion::readMesh ( )

private

[sourceFunction]

Definition at line 162 of file initial_diffusion.cpp.

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

readMesh() [2/2]

MoFEMErrorCode PhotonDiffusion::readMesh ( )

private

runProgram() [1/2]

MoFEMErrorCode PhotonDiffusion::runProgram

(

)

[runProgram]

Definition at line 439 of file initial_diffusion.cpp.

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

runProgram() [2/2]

MoFEMErrorCode PhotonDiffusion::runProgram

(

)

setIntegrationRules()

MoFEMErrorCode PhotonDiffusion::setIntegrationRules ( )

private

setupProblem() [1/2]

MoFEMErrorCode PhotonDiffusion::setupProblem ( )

private

Definition at line 173 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_NULLPTR, "", "-flux_magnitude",
                               &flux_magnitude, 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 PetscOptionsGetScalar(PETSC_NULLPTR, "", "-mu_a", &mu_a,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-mu_sp", &mu_sp,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-initial_time",
                               &initial_time, PETSC_NULLPTR);
 
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-output_file", out_file_name,
                               255, PETSC_NULLPTR);
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-output_volume",
                             &output_volume, PETSC_NULLPTR);
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-testing", &testing,
                             PETSC_NULLPTR);
 
  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_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
 
  MOFEM_LOG("INITIAL", Sev::inform) << "Approximation order: " << order;
  MOFEM_LOG("INITIAL", Sev::inform) << "Kronrod points: " << kronrod_points;
 
  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);
}

setupProblem() [2/2]

MoFEMErrorCode PhotonDiffusion::setupProblem ( )

private

solveSystem() [1/2]

MoFEMErrorCode PhotonDiffusion::solveSystem ( )

private

[assembleSystem]

[solveSystem]

[postprocess]

[postprocess]

[timesolver]

Definition at line 374 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);
 
  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);
}

solveSystem() [2/2]

MoFEMErrorCode PhotonDiffusion::solveSystem ( )

private

sourceFunction()

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

inlinestatic

Definition at line 74 of file initial_diffusion.cpp.

                                                      {
    return SourceFunction::sourceFunctionEval(
        x, y, z, beam_radius, beam_centre_x, beam_centre_y, slab_thickness,
        mu_a, mu_sp, flux_magnitude, initial_time, v, D);
  }

Member Data Documentation

boundaryLhsFEPtr

boost::shared_ptr<FEMethod> PhotonDiffusion::boundaryLhsFEPtr

private

Definition at line 131 of file photon_diffusion.cpp.

boundaryMarker

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

private

Definition at line 128 of file photon_diffusion.cpp.

boundaryRhsFEPtr

boost::shared_ptr<FEMethod> PhotonDiffusion::boundaryRhsFEPtr

private

Definition at line 132 of file photon_diffusion.cpp.

commonDataPtr

boost::shared_ptr< CommonData > PhotonDiffusion::commonDataPtr

private

Definition at line 103 of file initial_diffusion.cpp.

domainLhsFEPtr

boost::shared_ptr<FEMethod> PhotonDiffusion::domainLhsFEPtr

private

Definition at line 130 of file photon_diffusion.cpp.

mField

MoFEM::Interface & PhotonDiffusion::mField

private

Definition at line 95 of file initial_diffusion.cpp.

---

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

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