#include <stdlib.h>
#include <cmath>
#include <BasicFiniteElements.hpp>
static char help[] =
"...\n\n";
using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc,
boost::shared_ptr<PostProcFaceOnRefinedMesh> skin_post_proc)
: dM(dm), postProc(post_proc), skinPostProc(skin_post_proc){};
"out_volume_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
}
if (skinPostProc) {
CHKERR skinPostProc->writeFile(
"out_camera_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
}
}
}
private:
SmartPetscObj<DM> dM;
boost::shared_ptr<PostProcEle> postProc;
boost::shared_ptr<PostProcFaceOnRefinedMesh> skinPostProc;
};
public:
private:
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
boost::shared_ptr<FEMethod> domianLhsFEPtr;
boost::shared_ptr<FEMethod> boundaryLhsFEPtr;
boost::shared_ptr<FEMethod> boundaryRhsFEPtr;
};
}
PETSC_NULL);
PETSC_NULL);
PETSC_NULL);
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;
<<
"Absorption coefficient (cm^-1): " <<
mu_a;
<<
"Scattering coefficient (cm^-1): " <<
mu_sp;
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Impulse magnitude: " <<
flux;
PETSC_NULL);
auto set_camera_skin_fe = [&]() {
Range camera_surface;
const std::string block_name = "CAM";
bool add_fe = false;
if (
bit->getName().compare(0, block_name.size(), block_name) == 0) {
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Found CAM block";
bit->getMeshset(), 2, camera_surface,
true);
add_fe = true;
}
}
MOFEM_LOG(
"PHOTON", Sev::noisy) <<
"CAM block entities:\n"
<< camera_surface;
if (add_fe) {
"CAMERA_FE");
}
};
auto my_simple_set_up = [&]() {
}
};
}
};
}
}
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"MIX",
"U", 0,
0, false);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"SPOT",
"U",
0, 0, false);
}
auto &bc_map = bc_mng->getBcMapByBlockName();
auto add_domain_base_ops = [&](auto &pipeline) {
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_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, inv_jac_ptr));
pipeline.push_back(
new OpSetHOInvJacToScalarBases(
H1, inv_jac_ptr));
pipeline.push_back(new OpSetHOWeights(det_ptr));
};
auto add_domain_lhs_ops = [&](auto &pipeline) {
"U",
"U", [](
double,
double,
double) ->
double {
return D; }));
auto get_mass_coefficient = [&](const double, const double, const double) {
};
pipeline.push_back(
new OpDomainMass(
"U",
"U", get_mass_coefficient));
};
auto add_domain_rhs_ops = [&](auto &pipeline) {
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.push_back(new OpCalculateScalarFieldGradient<SPACE_DIM>(
"U", grad_u_at_gauss_pts));
pipeline.push_back(new OpCalculateScalarFieldValues("U", u_at_gauss_pts));
pipeline.push_back(
"U", grad_u_at_gauss_pts,
[](
double,
double,
double) ->
double {
return D; }));
"U", dot_u_at_gauss_pts,
[](
const double,
const double,
const double) {
return inv_v; }));
"U", u_at_gauss_pts,
[](
const double,
const double,
const double) {
return mu_a; }));
auto source_term = [&](const double, const double, const double) {
return 0;
};
};
auto add_boundary_base_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetHOWeightsOnFace());
};
auto add_lhs_base_ops = [&](auto &pipeline) {
if (std::regex_match(
b.first, std::regex(
"(.*)EXT(.*)"))) {
"U", "U",
[](
const double,
const double,
const double) {
return h; },
b.second->getBcEntsPtr()));
}
}
};
auto add_rhs_base_ops = [&](auto &pipeline) {
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
pipeline.push_back(new OpCalculateScalarFieldValues("U", u_at_gauss_pts));
if (std::regex_match(
b.first, std::regex(
"(.*)EXT(.*)"))) {
"U", u_at_gauss_pts,
[](
const double,
const double,
const double) {
return h; },
b.second->getBcEntsPtr()));
}
}
};
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());
add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
add_lhs_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_rhs_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
}
auto create_post_process_element = [&]() {
auto post_froc_fe = boost::make_shared<PostProcEle>(
mField);
post_froc_fe->generateReferenceElementMesh();
post_froc_fe->addFieldValuesPostProc("U");
post_froc_fe->addFieldValuesGradientPostProc("U");
return post_froc_fe;
};
auto create_post_process_camera_element = [&]() {
auto post_proc_skin =
boost::make_shared<PostProcFaceOnRefinedMesh>(
mField);
post_proc_skin->generateReferenceElementMesh();
CHKERR post_proc_skin->addFieldValuesPostProc(
"U");
CHKERR post_proc_skin->addFieldValuesGradientPostProcOnSkin(
"U",
simple->getDomainFEName());
return post_proc_skin;
} else {
return boost::shared_ptr<PostProcFaceOnRefinedMesh>();
}
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(dm, create_post_process_element(),
create_post_process_camera_element()));
boost::shared_ptr<ForcesAndSourcesCore> null;
monitor_ptr, null, null);
};
<< "reading vector in binary from vector.dat ...";
PetscViewer viewer;
PetscViewerBinaryOpen(PETSC_COMM_WORLD, "initial_vector.dat",
FILE_MODE_READ, &viewer);
}
auto solver = pipeline_mng->createTS();
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetFromOptions(solver);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
}
}
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "PHOTON"));
LogManager::setLog("PHOTON");
try {
DMType dm_name = "DMMOFEM";
CHKERR heat_problem.runProgram();
}
return 0;
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
int main(int argc, char *argv[])
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.
auto smartCreateDMVector
Get smart vector from DM.
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
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.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
#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, 2 > OpDomainGradGrad
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
double mu_sp
scattering coefficient (cm^-1)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpDomainTimesScalarField
double v
phase velocity of light in medium (cm/ns)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
double mu_a
absorption coefficient (cm^-1)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
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)
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
const double D
diffusivity
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBoundaryTimeScalarField
double flux
impulse magnitude
double h
convective heat coefficient
double duration
impulse duration (ns)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
static constexpr int approx_order
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
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.
Deprecated interface functions.
default operator for TRI element
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]
MoFEMErrorCode assembleSystem()
boost::shared_ptr< FEMethod > boundaryRhsFEPtr
MoFEMErrorCode solveSystem()
MoFEM::Interface & mField
MoFEMErrorCode readMesh()
MoFEMErrorCode outputResults()
MoFEMErrorCode initialCondition()
boost::shared_ptr< FEMethod > domianLhsFEPtr
PhotonDiffusion(MoFEM::Interface &m_field)
MoFEMErrorCode runProgram()
boost::shared_ptr< FEMethod > boundaryLhsFEPtr
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode setupProblem()