11#define BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND
15static char help[] =
"...\n\n";
40const double v =
c /
n;
58#include <boost/math/quadrature/gauss_kronrod.hpp>
59using namespace boost::math::quadrature;
62 OpError(boost::shared_ptr<VectorDouble> u_at_pts_ptr,
double &l2_error)
94 static double sourceFunction(
const double x,
const double y,
const double z) {
99 auto phi_pulse = [&](
const double r_s,
const double phi_s) {
100 const double xs = r_s * cos(phi_s);
101 const double ys = r_s * sin(phi_s);
106 const double P1 = xp * xp + yp * yp + zp1 * zp1;
107 const double P2 = xp * xp + yp * yp + zp2 * zp2;
108 return r_s * (exp(-P1 /
A) - exp(-P2 /
A));
111 auto f = [&](
const double r_s) {
112 auto g = [&](
const double phi_s) {
return phi_pulse(r_s, phi_s); };
113 return gauss_kronrod<double, 15>::integrate(
114 g, 0, 2 * M_PI, 0, std::numeric_limits<float>::epsilon());
118 gauss_kronrod<double, 15>::integrate(
119 f, 0,
beam_radius, 0, std::numeric_limits<float>::epsilon());
182 MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Refractive index: " <<
n;
183 MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Speed of light (cm/ns): " <<
c;
185 <<
"Phase velocity in medium (cm/ns): " <<
v;
187 <<
"Absorption coefficient (cm^-1): " <<
mu_a;
189 <<
"Scattering coefficient (cm^-1): " <<
mu_sp;
190 MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Diffusion coefficient D : " <<
D;
198 MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Approximation order: " <<
order;
265 std::string entity_name = it->getName();
266 if (entity_name.compare(0, 3,
"INT") == 0) {
268 boundary_ents,
true);
272 Range boundary_verts;
275 boundary_ents.merge(boundary_verts);
279 simple->getProblemName(),
"PHOTON_FLUENCE_RATE", boundary_ents);
290 auto add_domain_base_ops = [&](
auto &pipeline) {
291 auto jac_ptr = boost::make_shared<MatrixDouble>();
292 auto det_ptr = boost::make_shared<VectorDouble>();
298 auto add_domain_lhs_ops = [&](
auto &pipeline) {
300 "PHOTON_FLUENCE_RATE",
"PHOTON_FLUENCE_RATE",
301 [](
const double,
const double,
const double) {
return 1; }));
304 auto add_domain_rhs_ops = [&](
auto &pipeline) {
310 add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
311 add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
312 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
313 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
322 auto solver = pipeline_mng->createKSP();
324 CHKERR KSPSetFromOptions(solver);
327 auto dm =
simple->getDM();
331 MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Solver start";
333 CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
334 CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
338 <<
"writing vector in binary to " <<
out_file_name <<
" ...";
340 PetscViewerBinaryOpen(PETSC_COMM_WORLD,
out_file_name, FILE_MODE_WRITE,
343 PetscViewerDestroy(&viewer);
345 MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Solver done";
353 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
355 auto u_ptr = boost::make_shared<VectorDouble>();
356 post_proc_fe->getOpPtrVector().push_back(
361 post_proc_fe->getOpPtrVector().push_back(
365 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
367 {{
"PHOTON_FLUENCE_RATE", u_ptr}},
377 pipeline_mng->getDomainRhsFE() = post_proc_fe;
378 CHKERR pipeline_mng->loopFiniteElements();
379 CHKERR post_proc_fe->writeFile(
"out_initial.h5m");
399int main(
int argc,
char *argv[]) {
406 auto core_log = logging::core::get();
415 DMType dm_name =
"DMMOFEM";
419 moab::Core mb_instance;
420 moab::Interface &moab = mb_instance;
void simple(double P1[], double P2[], double P3[], double c[], const int N)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
auto createDMVector(DM dm)
Get smart vector from DM.
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.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
double mu_sp
scattering coefficient (cm^-1)
double flux_magnitude
impulse magnitude
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
constexpr int SPACE_DIM
[Define dimension]
const double c
speed of light (cm/ns)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
double mu_a
absorption coefficient (cm^-1)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
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)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
OpCalculateHOJacForVolume OpCalculateHOJacVolume
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
double mu_sp
scattering coefficient (cm^-1)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
const double c
speed of light (cm/ns)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
double mu_a
absorption coefficient (cm^-1)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
static constexpr int approx_order
Simple interface for fast problem set-up.
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
virtual moab::Interface & get_moab()=0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
@ 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.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
boost::shared_ptr< VectorDouble > uAtPtsPtr
OpError(boost::shared_ptr< VectorDouble > u_at_pts_ptr, double &l2_error)
MoFEMErrorCode assembleSystem()
MoFEMErrorCode solveSystem()
MoFEM::Interface & mField
MoFEMErrorCode readMesh()
MoFEMErrorCode outputResults()
MoFEMErrorCode initialCondition()
PhotonDiffusion(MoFEM::Interface &m_field)
MoFEMErrorCode checkResults()
MoFEMErrorCode runProgram()
MoFEMErrorCode boundaryCondition()
static double sourceFunction(const double x, const double y, const double z)
Pulse is infinitely short.
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode setupProblem()