16static char help[] =
"...\n\n";
111 boost::shared_ptr<VectorDouble> &approx_vals,
double &sum_temp,
140 constexpr int order = 3;
150 approxVals = boost::make_shared<VectorDouble>();
160 auto set_initial_temperature = [&](
VectorAdaptor &&field_data,
double *xcoord,
161 double *ycoord,
double *zcoord) {
181 auto det_ptr = boost::make_shared<VectorDouble>();
182 auto jac_ptr = boost::make_shared<MatrixDouble>();
183 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
237 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
238 CHKERR TSSetFromOptions(ts);
239 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
246 CHKERR TSGetTime(ts, &ftime);
248 PetscInt steps, snesfails, rejects, nonlinits, linits;
249 CHKERR TSGetTimeStepNumber(ts, &steps);
250 CHKERR TSGetSNESFailures(ts, &snesfails);
251 CHKERR TSGetStepRejections(ts, &rejects);
252 CHKERR TSGetSNESIterations(ts, &nonlinits);
253 CHKERR TSGetKSPIterations(ts, &linits);
255 "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
257 steps, rejects, snesfails, ftime, nonlinits, linits);
271 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
273 auto t_ptr = boost::make_shared<VectorDouble>();
274 post_proc_fe->getOpPtrVector().push_back(
279 post_proc_fe->getOpPtrVector().push_back(
281 new OpPPMap(post_proc_fe->getPostProcMesh(),
282 post_proc_fe->getMapGaussPts(),
297 double sum_temperature;
300 new OpCalcSurfaceAverageTemperature(approxVals, sum_temperature,
307 CHKERR post_proc_fe->writeFile(
"out_radiation.h5m");
309 MOFEM_LOG_C(
"EXAMPLE", Sev::inform,
"Surface area %3.4e [km]", surface_area);
311 "Average subsurface temperatute %3.4e [K]",
312 sum_temperature / surface_area);
322int main(
int argc,
char *argv[]) {
325 const char param_file[] =
"param_file.petsc";
329 auto core_log = logging::core::get();
338 DMType dm_name =
"DMMOFEM";
343 moab::Core mb_instance;
344 moab::Interface &moab = mb_instance;
373 const double vol = this->getMeasure();
375 auto t_w = getFTensor0IntegrationWeight();
381 auto t_coords = getFTensor1CoordsAtGaussPts();
387 const double r_cylinder = t_coords(0);
390 const double alpha = t_w * vol *
Beta * (2 * M_PI * r_cylinder);
392 for (
int rr = 0; rr !=
nbRows; ++rr) {
395 for (
int cc = 0; cc !=
nbCols; cc++) {
396 if (std::abs(t_coords(0)) > std::numeric_limits<double>::epsilon()) {
397 locMat(rr, cc) += alpha * t_row_base * t_col_base * 4 *
398 std::pow(
static_cast<double>(t_val), 3);
417 const double vol = getMeasure();
419 auto t_w = getFTensor0IntegrationWeight();
425 auto t_coords = getFTensor1CoordsAtGaussPts();
428 for (
int gg = 0; gg != nbIntegrationPts; gg++) {
431 const double r_cylinder = t_coords(0);
434 const double alpha = t_w * vol *
Beta * (2 * M_PI * r_cylinder);
436 for (
int rr = 0; rr != nbRows; ++rr) {
437 if (std::abs(t_coords(0)) > std::numeric_limits<double>::epsilon()) {
438 locF[rr] += alpha * t_row_base *
439 (std::pow(
static_cast<double>(t_val), 4) -
457 const double vol = getMeasure();
459 auto t_w = getFTensor0IntegrationWeight();
463 auto t_coords = getFTensor1CoordsAtGaussPts();
465 const double time = getFEMethod()->ts_t;
468 constexpr double flux_p = -0.03e6;
469 constexpr double flux_c = -0.23e6;
472 for (
int gg = 0; gg != nbIntegrationPts; gg++) {
475 const double r_cylinder = t_coords(0);
477 const double r = std::sqrt(t_coords(
i) * t_coords(
i));
478 const double s = std::abs(t_coords(1)) / r;
481 const double alpha = t_w * vol * (2 * M_PI * r_cylinder);
483 for (
int rr = 0; rr != nbRows; ++rr) {
484 locF[rr] += alpha * t_row_base * (s * flux_p + flux_c) * time;
500 if (type == MBVERTEX) {
502 const double vol = getMeasure();
504 auto t_w = getFTensor0IntegrationWeight();
508 auto t_coords = getFTensor1CoordsAtGaussPts();
510 size_t nb_integration_pts = getGaussPts().size2();
513 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
516 const double r_cylinder = t_coords(0);
519 const double alpha = t_w * vol * (2 * M_PI * r_cylinder);
521 sumTemperature += alpha * t_val;
522 surfaceArea += alpha;
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#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.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< TS > createTSIM(SmartPetscObj< DM > dm=nullptr)
Create TS (time) implicit solver.
boost::ptr_deque< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorShallowArrayAdaptor< double > VectorAdaptor
implementation of Data Operators for Forces and Sources
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr double boltzmann_constant
constexpr double T_ambient
constexpr double emissivity
constexpr double heat_conductivity
OpCalcSurfaceAverageTemperature(boost::shared_ptr< VectorDouble > &approx_vals, double &sum_temp, double &surf)
boost::shared_ptr< VectorDouble > approxVals
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
[Flux Rhs]
FTensor::Index< 'i', 2 > i
summit Index
MoFEMErrorCode iNtegrate(EntData &row_data)
[Radiation Lhs]
boost::shared_ptr< VectorDouble > approxVals
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
[Radiation Lhs]
OpRadiationLhs(boost::shared_ptr< VectorDouble > &approx_vals)
OpRadiationRhs(boost::shared_ptr< VectorDouble > &approx_vals)
boost::shared_ptr< VectorDouble > approxVals
MoFEMErrorCode iNtegrate(EntData &row_data)
[Radiation Lhs]
static int integrationRule(int, int, int p_data)
boost::shared_ptr< VectorDouble > approxVals
MoFEMErrorCode kspSolve()
[Push operators to pipeline]
MoFEMErrorCode checkResults()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
boost::shared_ptr< MatrixDouble > approxGradVals
MoFEM::Interface & mField
MoFEMErrorCode postProcess()
MoFEMErrorCode setupProblem()
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 EDGE element
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user function.
@ OPROW
operator doWork function is executed on FE rows
friend class UserDataOperator
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.
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
boost::shared_ptr< FEMethod > & getBoundaryLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.