static char help[] =
"...\n\n";
#include <BasicFiniteElements.hpp>
using OpBase = OpBaseImpl<PETSC, EdgeEleOp>;
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) {}
};
};
}
Simple *
simple = mField.getInterface<Simple>();
}
approxVals = boost::make_shared<VectorDouble>();
approxGradVals = boost::make_shared<MatrixDouble>();
}
auto set_initial_temperature = [&](
VectorAdaptor &&field_data,
double *xcoord,
double *ycoord, double *zcoord) {
};
FieldBlas *field_blas;
CHKERR mField.getInterface(field_blas);
CHKERR field_blas->setVertexDofs(set_initial_temperature,
"T");
}
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto beta = [](
const double r,
const double,
const double) {
};
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetInvJacH1ForFace(inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integrationRule);
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSetHOWeightsOnFace());
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetInvJacH1ForFace(inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldGradient<2>("T", approxGradVals));
pipeline_mng->getOpDomainRhsPipeline().push_back(
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integrationRule);
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpCalculateScalarFieldValues("T", approxVals));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpRadiationRhs(approxVals));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpFluxRhs());
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integrationRule);
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpCalculateScalarFieldValues("T", approxVals));
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpRadiationLhs(approxVals));
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integrationRule);
}
Simple *
simple = mField.getInterface<Simple>();
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto ts = pipeline_mng->createTSIM();
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);
}
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
post_proc_fe->generateReferenceElementMesh();
post_proc_fe->addFieldValuesPostProc("T");
pipeline_mng->getDomainRhsFE() = post_proc_fe;
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpCalculateScalarFieldValues("T", approxVals));
double sum_temperature;
double surface_area;
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpCalcSurfaceAverageTemperature(approxVals, sum_temperature,
surface_area));
auto calc_surfcae_area_op = pipeline_mng->getOpBoundaryRhsPipeline().back();
sum_temperature = 0;
surface_area = 0;
CHKERR pipeline_mng->loopFiniteElements();
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[]) {
auto core_log = logging::core::get();
core_log->add_sink(
try {
DMType dm_name = "DMMOFEM";
}
}
const double vol = this->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) {
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()) {
}
++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;
++t_coords;
++t_val;
++t_w;
}
}
}
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
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.
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
constexpr double boltzmann_constant
constexpr double T_ambient
constexpr double emissivity
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
[Flux Rhs]
MoFEMErrorCode iNtegrate(EntData &row_data)
[Radiation Lhs]
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
[Radiation Lhs]
MoFEMErrorCode iNtegrate(EntData &row_data)
[Radiation Lhs]
MoFEMErrorCode kspSolve()
[Push operators to pipeline]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode postProcess()
[Solve]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode bC()
[Create common data]
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.
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.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator