static char help[] =
"...\n\n";
private:
static int integrationRule(int, int, int p_data) { return 2 * p_data; };
boost::shared_ptr<VectorDouble> approxVals;
boost::shared_ptr<MatrixDouble> approxGradVals;
struct OpRadiationLhs :
public OpBase {
private:
boost::shared_ptr<VectorDouble> approxVals;
public:
OpRadiationLhs(boost::shared_ptr<VectorDouble> &approx_vals)
:
OpBase(
"T",
"T",
OpBase::OPROWCOL), approxVals(approx_vals) {
this->sYmm = false;
}
};
struct OpRadiationRhs :
public OpBase {
private:
boost::shared_ptr<VectorDouble> approxVals;
public:
OpRadiationRhs(boost::shared_ptr<VectorDouble> &approx_vals)
:
OpBase(
"T",
"T",
OpBase::OPROW), approxVals(approx_vals) {}
};
struct OpFluxRhs :
public OpBase {
private:
public:
};
struct OpCalcSurfaceAverageTemperature :
public EdgeEleOp {
private:
boost::shared_ptr<VectorDouble> approxVals;
double &sumTemperature;
double &surfaceArea;
public:
OpCalcSurfaceAverageTemperature(
boost::shared_ptr<VectorDouble> &approx_vals, double &sum_temp,
double &surf)
sumTemperature(sum_temp), surfaceArea(surf) {}
};
};
}
}
approxVals = boost::make_shared<VectorDouble>();
approxGradVals = boost::make_shared<MatrixDouble>();
}
auto set_initial_temperature = [&](
VectorAdaptor &&field_data,
double *xcoord,
double *ycoord, double *zcoord) {
};
}
};
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
new OpRadiationRhs(approxVals));
new OpRadiationLhs(approxVals));
}
double ftime = 1;
CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
SCATTER_FORWARD);
PetscInt steps, snesfails, rejects, nonlinits, linits;
CHKERR TSGetTimeStepNumber(ts, &steps);
CHKERR TSGetSNESFailures(ts, &snesfails);
CHKERR TSGetStepRejections(ts, &rejects);
CHKERR TSGetSNESIterations(ts, &nonlinits);
CHKERR TSGetKSPIterations(ts, &linits);
"steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
"%d, linits %d",
steps, rejects, snesfails, ftime, nonlinits, linits);
}
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
auto t_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"T", t_ptr}},
{}, {}, {}
)
);
double sum_temperature;
double surface_area;
new OpCalcSurfaceAverageTemperature(approxVals, sum_temperature,
surface_area));
sum_temperature = 0;
surface_area = 0;
CHKERR post_proc_fe->writeFile(
"out_radiation.h5m");
MOFEM_LOG_C(
"EXAMPLE", Sev::inform,
"Surface area %3.4e [km]", surface_area);
"Average subsurface temperatute %3.4e [K]",
sum_temperature / surface_area);
}
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(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
try {
DMType dm_name = "DMMOFEM";
}
}
const double vol = this->getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
const double r_cylinder = t_coords(0);
const double alpha = t_w * vol *
Beta * (2 * M_PI * r_cylinder);
for (
int rr = 0; rr !=
nbRows; ++rr) {
for (
int cc = 0; cc !=
nbCols; cc++) {
if (std::abs(t_coords(0)) > std::numeric_limits<double>::epsilon()) {
locMat(rr, cc) += alpha * t_row_base * t_col_base * 4 * pow(t_val, 3);
}
++t_col_base;
}
++t_row_base;
}
++t_val;
++t_coords;
++t_w;
}
}
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
const double r_cylinder = t_coords(0);
const double alpha = t_w * vol *
Beta * (2 * M_PI * r_cylinder);
for (int rr = 0; rr != nbRows; ++rr) {
if (std::abs(t_coords(0)) > std::numeric_limits<double>::epsilon()) {
locF[rr] += alpha * t_row_base * (pow(t_val, 4) - pow(
T_ambient, 4));
}
++t_row_base;
}
++t_coords;
++t_val;
++t_w;
}
}
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
const double time = getFEMethod()->ts_t;
constexpr double flux_p = -0.03e6;
constexpr double flux_c = -0.23e6;
for (int gg = 0; gg != nbIntegrationPts; gg++) {
const double r_cylinder = t_coords(0);
const double r = std::sqrt(t_coords(
i) * t_coords(
i));
const double s = std::abs(t_coords(1)) /
r;
const double alpha = t_w * vol * (2 * M_PI * r_cylinder);
for (int rr = 0; rr != nbRows; ++rr) {
locF[rr] += alpha * t_row_base * (s * flux_p + flux_c) * time;
++t_row_base;
}
++t_coords;
++t_w;
}
}
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
size_t nb_integration_pts = getGaussPts().size2();
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
const double r_cylinder = t_coords(0);
const double alpha = t_w * vol * (2 * M_PI * r_cylinder);
sumTemperature += alpha * t_val;
surfaceArea += alpha;
++t_coords;
++t_val;
++t_w;
}
}
}