#include <stdlib.h>
#include <cmath>
static char help[] =
"...\n\n";
PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
public:
private:
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<VectorDouble> approxVals;
enum VecElements {
VALUES_INTEG = 0,
LAST_ELEMENT
};
};
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<CommonData> commonDataPtr;
OpCameraInteg(boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {
std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
doEntities[MBTRI] = doEntities[MBQUAD] = true;
}
};
struct OpGetScalarFieldGradientValuesOnSkin :
public BoundaryEleOp {
boost::shared_ptr<VolSideFe> sideOpFe;
OpGetScalarFieldGradientValuesOnSkin(boost::shared_ptr<VolSideFe> side_fe)
:
BoundaryEleOp(
"PHOTON_FLUENCE_RATE", OPROW), sideOpFe(side_fe) {}
CHKERR loopSideVolumes(
"dFE", *sideOpFe);
}
};
boost::shared_ptr<PostProcFaceEle> skin_post_proc,
boost::shared_ptr<BoundaryEle> skin_post_proc_integ,
boost::shared_ptr<CommonData> common_data_ptr)
: dM(dm), postProc(post_proc), skinPostProc(skin_post_proc),
skinPostProcInteg(skin_post_proc_integ),
commonDataPtr(common_data_ptr){};
}
}
CHKERR VecZeroEntries(commonDataPtr->petscVec);
CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, ADD_VALUES,
SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, ADD_VALUES,
SCATTER_REVERSE);
CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
const double *array;
CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Fluence rate integral: " << array[0];
CHKERR postProc->writeFile(
"out_volume_" +
boost::lexical_cast<std::string>(ts_step) +
".h5m");
}
CHKERR skinPostProc->writeFile(
"out_camera_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
}
}
}
private:
boost::shared_ptr<PostProcEle> postProc;
boost::shared_ptr<PostProcFaceEle> skinPostProc;
boost::shared_ptr<BoundaryEle> skinPostProcInteg;
boost::shared_ptr<CommonData> commonDataPtr;
};
};
}
PetscInt ghosts[1] = {0};
else
}
PETSC_NULL);
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) <<
"Diffusion coefficient D : " <<
D;
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Coefficient A : " <<
A;
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Coefficient h : " <<
h;
auto set_camera_skin_fe = [&]() {
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) {
"PHOTON_FLUENCE_RATE");
"CAMERA_FE");
}
};
auto my_simple_set_up = [&]() {
}
};
}
};
}
}
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"EXT",
"PHOTON_FLUENCE_RATE", 0, 0, false);
std::string entity_name = it->getName();
if (entity_name.compare(0, 3, "INT") == 0) {
boundary_ents, true);
}
}
true);
boundary_ents.merge(boundary_verts);
simple->getProblemName(),
"PHOTON_FLUENCE_RATE", boundary_ents);
}
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>();
};
auto add_domain_lhs_ops = [&](auto &pipeline) {
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
[](
double,
double,
double) ->
double {
return D; }));
};
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE", 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>();
"PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts));
u_at_gauss_pts));
"PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts));
"PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts,
[](
double,
double,
double) ->
double {
return D; }));
"PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts,
[](
const double,
const double,
const double) {
return inv_v; }));
"PHOTON_FLUENCE_RATE", u_at_gauss_pts,
[](
const double,
const double,
const double) {
return mu_a; }));
};
auto add_boundary_base_ops = [&](auto &pipeline) {
};
auto add_boundary_lhs_ops = [&](auto &pipeline) {
for (auto b : bc_map) {
if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
[](
const double,
const double,
const double) {
return h; },
b.second->getBcEntsPtr()));
}
}
};
auto add_boundary_rhs_ops = [&](auto &pipeline) {
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
u_at_gauss_pts));
for (auto b : bc_map) {
if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
"PHOTON_FLUENCE_RATE", 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_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
}
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(
post_froc_fe->getOpPtrVector().push_back(
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 = [&]() {
auto post_proc_skin = boost::make_shared<PostProcFaceEle>(mField);
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
op_loop_side->getOpPtrVector().push_back(
op_loop_side->getOpPtrVector().push_back(
op_loop_side->getOpPtrVector().push_back(
op_loop_side->getOpPtrVector().push_back(
op_loop_side->getOpPtrVector().push_back(
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 = [&]() {
auto post_proc_integ_skin = boost::make_shared<BoundaryEle>(mField);
post_proc_integ_skin->getOpPtrVector().push_back(
post_proc_integ_skin->getOpPtrVector().push_back(
commonDataPtr->approxVals));
post_proc_integ_skin->getOpPtrVector().push_back(
new OpCameraInteg(commonDataPtr));
return post_proc_integ_skin;
} else {
return boost::shared_ptr<BoundaryEle>();
}
};
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(),
create_post_process_integ_camera_element(), commonDataPtr));
boost::shared_ptr<ForcesAndSourcesCore> null;
monitor_ptr, null, null);
};
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"reading vector in binary from file "
PetscViewer viewer;
&viewer);
VecLoad(X, viewer);
}
auto solver = pipeline_mng->createTS();
CHKERR TSSetSolution(solver, X);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetSolution(solver, X);
CHKERR TSSetFromOptions(solver);
CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
}
}
}
const int nb_integration_pts = getGaussPts().size2();
const double area = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
double values_integ = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * area;
values_integ += alpha * t_val;
++t_w;
++t_val;
}
std::array<double, 1> values;
values[0] = values_integ;
ADD_VALUES);
}
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
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;
}