26 using namespace MoFEM;
28 static char help[] =
"...\n\n";
30 #include <BasicFiniteElements.hpp>
84 :
OpBase(
"T",
"T",
OpBase::OPROWCOL), approxVals(approx_vals) {
98 :
OpBase(
"T",
"T",
OpBase::OPROW), approxVals(approx_vals) {}
124 boost::shared_ptr<VectorDouble> &approx_vals,
double &sum_temp,
127 sumTemperature(sum_temp), surfaceArea(surf) {}
141 CHKERR createCommonData();
157 constexpr
int order = 3;
167 approxVals = boost::make_shared<VectorDouble>();
168 approxGradVals = boost::make_shared<MatrixDouble>();
177 auto set_initial_temperature = [&](
VectorAdaptor &&field_data,
double *xcoord,
178 double *ycoord,
double *zcoord) {
184 CHKERR mField.getInterface(field_blas);
194 auto beta = [](
const double r,
const double,
const double) {
198 auto det_ptr = boost::make_shared<VectorDouble>();
199 auto jac_ptr = boost::make_shared<MatrixDouble>();
200 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
229 new OpRadiationRhs(approxVals));
236 new OpRadiationLhs(approxVals));
250 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
251 CHKERR TSSetFromOptions(ts);
252 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
259 CHKERR TSGetTime(ts, &ftime);
261 PetscInt steps, snesfails, rejects, nonlinits, linits;
262 CHKERR TSGetTimeStepNumber(ts, &steps);
263 CHKERR TSGetSNESFailures(ts, &snesfails);
264 CHKERR TSGetStepRejections(ts, &rejects);
265 CHKERR TSGetSNESIterations(ts, &nonlinits);
266 CHKERR TSGetKSPIterations(ts, &linits);
268 "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
270 steps, rejects, snesfails, ftime, nonlinits, linits);
283 auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
284 post_proc_fe->generateReferenceElementMesh();
285 post_proc_fe->addFieldValuesPostProc(
"T");
291 double sum_temperature;
294 new OpCalcSurfaceAverageTemperature(approxVals, sum_temperature,
301 CHKERR post_proc_fe->writeFile(
"out_radiation.h5m");
303 MOFEM_LOG_C(
"EXAMPLE", Sev::inform,
"Surface area %3.4e [km]", surface_area);
305 "Average subsurface temperatute %3.4e [K]",
306 sum_temperature / surface_area);
316 int main(
int argc,
char *argv[]) {
323 auto core_log = logging::core::get();
332 DMType dm_name =
"DMMOFEM";
367 const double vol = this->getMeasure();
369 auto t_w = getFTensor0IntegrationWeight();
375 auto t_coords = getFTensor1CoordsAtGaussPts();
378 for (
int gg = 0; gg != nbIntegrationPts; gg++) {
381 const double r_cylinder = t_coords(0);
384 const double alpha = t_w * vol *
Beta * (2 * M_PI * r_cylinder);
386 for (
int rr = 0; rr != nbRows; ++rr) {
389 for (
int cc = 0; cc != nbCols; cc++) {
390 if (std::abs(t_coords(0)) > std::numeric_limits<double>::epsilon()) {
391 locMat(rr, cc) +=
alpha * t_row_base * t_col_base * 4 * pow(t_val, 3);
410 const double vol = getMeasure();
412 auto t_w = getFTensor0IntegrationWeight();
418 auto t_coords = getFTensor1CoordsAtGaussPts();
421 for (
int gg = 0; gg != nbIntegrationPts; gg++) {
424 const double r_cylinder = t_coords(0);
427 const double alpha = t_w * vol *
Beta * (2 * M_PI * r_cylinder);
429 for (
int rr = 0; rr != nbRows; ++rr) {
430 if (std::abs(t_coords(0)) > std::numeric_limits<double>::epsilon()) {
431 locF[rr] +=
alpha * t_row_base * (pow(t_val, 4) - pow(
T_ambient, 4));
448 const double vol = getMeasure();
450 auto t_w = getFTensor0IntegrationWeight();
454 auto t_coords = getFTensor1CoordsAtGaussPts();
456 const double time = getFEMethod()->ts_t;
459 constexpr
double flux_p = -0.03e6;
460 constexpr
double flux_c = -0.23e6;
463 for (
int gg = 0; gg != nbIntegrationPts; gg++) {
466 const double r_cylinder = t_coords(0);
468 const double r = std::sqrt(t_coords(
i) * t_coords(
i));
469 const double s = std::abs(t_coords(1)) /
r;
472 const double alpha = t_w * vol * (2 * M_PI * r_cylinder);
474 for (
int rr = 0; rr != nbRows; ++rr) {
475 locF[rr] +=
alpha * t_row_base * (s * flux_p + flux_c) * time;
491 if (
type == MBVERTEX) {
493 const double vol = getMeasure();
495 auto t_w = getFTensor0IntegrationWeight();
499 auto t_coords = getFTensor1CoordsAtGaussPts();
501 size_t nb_integration_pts = getGaussPts().size2();
504 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
507 const double r_cylinder = t_coords(0);
510 const double alpha = t_w * vol * (2 * M_PI * r_cylinder);
512 sumTemperature +=
alpha * t_val;
513 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 CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#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 ...
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 > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< TS > createTSIM(SmartPetscObj< DM > dm=nullptr)
Create TS (time) implicit solver.
boost::ptr_vector< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_vector< 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, 2 > 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
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
DeprecatedCoreInterface Interface
const double r
rate factor
int main(int argc, char *argv[])
[Check results]
EntitiesFieldData::EntData EntData
constexpr double boltzmann_constant
constexpr double T_ambient
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, 2 > OpDomainGradTimesVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 2 > OpDomainGradGrad
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
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.
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.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
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 getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator