23 #include <BasicFiniteElements.hpp>
25 using namespace MoFEM;
27 static char help[] =
"...\n\n";
63 #include <boost/math/quadrature/gauss_kronrod.hpp>
64 using namespace boost::math::quadrature;
67 OpError(boost::shared_ptr<VectorDouble> u_at_pts_ptr,
double &l2_error)
68 :
DomainEleOp(
"U", OPROW), uAtPtsPtr(u_at_pts_ptr), l2Error(l2_error) {
70 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
109 auto phi_pulse = [&](
const double r_s,
const double phi_s) {
110 const double xs = r_s * cos(phi_s);
111 const double ys = r_s * sin(phi_s);
112 const double xp = x - xs;
113 const double yp = y - ys;
116 const double P1 = xp * xp + yp * yp + zp1 * zp1;
117 const double P2 = xp * xp + yp * yp + zp2 * zp2;
118 return r_s * (exp(-P1 /
A) - exp(-P2 /
A));
121 auto f = [&](
const double r_s) {
122 auto g = [&](
const double phi_s) {
return phi_pulse(r_s, phi_s); };
123 return gauss_kronrod<double, 42>::integrate(
124 g, 0, 2 * M_PI, 0, std::numeric_limits<float>::epsilon());
128 gauss_kronrod<double, 32>::integrate(
129 f, 0,
spot_radius, 0, std::numeric_limits<float>::epsilon());
180 MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Refractive index: " <<
n;
181 MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Speed of light (cm/ns): " <<
c;
183 <<
"Phase velocity in medium (cm/ns): " <<
v;
185 <<
"Absorption coefficient (cm^-1): " <<
mu_a;
187 <<
"Scattering coefficient (cm^-1): " <<
mu_sp;
194 MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Approximation order: " <<
order;
235 auto add_domain_base_ops = [&](
auto &pipeline) {
236 auto jac_ptr = boost::make_shared<MatrixDouble>();
237 auto det_ptr = boost::make_shared<VectorDouble>();
243 auto add_domain_lhs_ops = [&](
auto &pipeline) {
245 "U",
"U", [](
const double,
const double,
const double) {
return 1; }));
248 auto add_domain_rhs_ops = [&](
auto &pipeline) {
253 add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
254 add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
255 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
256 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
265 auto solver = pipeline_mng->createKSP();
267 CHKERR KSPSetFromOptions(solver);
270 auto dm =
simple->getDM();
274 MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Solver start";
276 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
277 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
281 <<
"writing vector in binary to vector.dat ...";
283 PetscViewerBinaryOpen(PETSC_COMM_WORLD,
"initial_vector.dat", FILE_MODE_WRITE,
286 PetscViewerDestroy(&viewer);
288 MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Solver done";
298 auto u_vals_at_gauss_pts = boost::make_shared<VectorDouble>();
303 new OpError(u_vals_at_gauss_pts, l2_error));
305 MOFEM_LOG(
"INITIAL", Sev::inform) <<
"L2 error " << l2_error;
313 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
314 post_proc_fe->generateReferenceElementMesh();
315 post_proc_fe->addFieldValuesPostProc(
"U");
316 pipeline_mng->getDomainRhsFE() = post_proc_fe;
317 CHKERR pipeline_mng->loopFiniteElements();
318 CHKERR post_proc_fe->writeFile(
"out_initial.h5m");
338 int main(
int argc,
char *argv[]) {
345 auto core_log = logging::core::get();
347 LogManager::createSink(LogManager::getStrmWorld(),
"INITIAL"));
348 LogManager::setLog(
"INITIAL");
354 DMType dm_name =
"DMMOFEM";
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.
auto smartCreateDMVector
Get smart vector from 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.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 2 > OpDomainGradGrad
int main(int argc, char *argv[])
double mu_sp
scattering coefficient (cm^-1)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
double flux_magnitude
impulse magnitude
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpDomainTimesScalarField
double c
speed of light (cm/ns)
double v
phase velocity of light in medium (cm/ns)
constexpr int SPACE_DIM
[Define dimension]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
double n
refractive index of diffusive medium
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
double mu_a
absorption coefficient (cm^-1)
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 PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
static constexpr int approx_order
Simple interface for fast problem set-up.
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
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.
Calculate jacobian on Hex or other volume which is not simplex.
Get value at integration points for scalar field.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
boost::shared_ptr< VectorDouble > uAtPtsPtr
OpError(boost::shared_ptr< VectorDouble > u_at_pts_ptr, double &l2_error)
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
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()