11using namespace boost::numeric;
20 int nb_gauss_pts = data.
getN().size1();
24 if (type == MBVERTEX) {
29 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
42 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
50 Nf.resize(nb_row_dofs);
53 for (
unsigned int gg = 0; gg < data.
getN().size1(); gg++) {
56 dAta.cOnductivity_mat * getVolume() * getGaussPts()(3, gg);
59 ublas::noalias(Nf) += prod(prod(data.
getDiffN(gg, nb_row_dofs), val),
80 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
90 int nb_row = row_data.
getN().size2();
91 int nb_col = col_data.
getN().size2();
92 K.resize(nb_row, nb_col);
94 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
97 dAta.cOnductivity_mat * getVolume() * getGaussPts()(3, gg);
101 noalias(
K) += prod(K1, trans(col_data.
getDiffN(gg, nb_col)));
105 const_cast<FEMethod *
>(getFEMethod())->ts_B =
A;
108 nb_col, &col_data.
getIndices()[0], &
K(0, 0), ADD_VALUES);
109 if (row_side != col_side || row_type != col_type) {
110 transK.resize(nb_col, nb_row);
111 noalias(transK) = trans(
K);
114 &row_data.
getIndices()[0], &transK(0, 0), ADD_VALUES);
124 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
131 int nb_row = data.
getN().size2();
134 for (
unsigned int gg = 0; gg < data.
getN().size1(); gg++) {
135 double val = getGaussPts()(3, gg);
141 ublas::noalias(Nf) += val * data.
getN(gg);
143 Nf *= getVolume() * dAta.cApacity;
157 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
167 int nb_row = row_data.
getN().size2();
168 int nb_col = col_data.
getN().size2();
169 M.resize(nb_row, nb_col);
172 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
174 double val = getGaussPts()(3, gg);
178 val * outer_prod(row_data.
getN(gg, nb_row), col_data.
getN(gg, nb_col));
181 M *= getVolume() * dAta.cApacity * getFEMethod()->ts_a;
184 nb_col, &col_data.
getIndices()[0], &
M(0, 0), ADD_VALUES);
185 if (row_side != col_side || row_type != col_type) {
186 transM.resize(nb_col, nb_row);
187 noalias(transM) = trans(
M);
190 &row_data.
getIndices()[0], &transM(0, 0), ADD_VALUES);
201 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
209 int rank = dof_ptr->getNbOfCoeffs();
211 int nb_dofs = data.
getIndices().size() / rank;
216 for (
unsigned int gg = 0; gg < data.
getN().size1(); gg++) {
218 double val = getGaussPts()(2, gg);
221 const double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
222 flux = dAta.dAta.data.value1 * area;
224 flux = dAta.dAta.data.value1 * getArea();
226 ublas::noalias(Nf) += val * flux * data.
getN(gg, nb_dofs);
229 if (useTsF ||
F == PETSC_NULL) {
246 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
256 int nb_row = row_data.
getN().size2();
257 int nb_col = col_data.
getN().size2();
259 N.resize(nb_row, nb_col);
262 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
265 double radiationConst;
267 double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
268 radiationConst = dAta.sIgma * dAta.eMissivity * area;
270 radiationConst = dAta.sIgma * dAta.eMissivity * getArea();
272 const double fOur = 4.0;
273 double val = fOur * getGaussPts()(2, gg) * radiationConst * T3_at_Gauss_pt;
275 val * outer_prod(row_data.
getN(gg, nb_row), col_data.
getN(gg, nb_col));
279 const_cast<FEMethod *
>(getFEMethod())->ts_B =
A;
282 nb_col, &col_data.
getIndices()[0], &
N(0, 0), ADD_VALUES);
283 if (row_side != col_side || row_type != col_type) {
284 transN.resize(nb_col, nb_row);
285 noalias(transN) = trans(
N);
288 &row_data.
getIndices()[0], &transN(0, 0), ADD_VALUES);
298 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
306 int rank = dof_ptr->getNbOfCoeffs();
307 int nb_row_dofs = data.
getIndices().size() / rank;
312 for (
unsigned int gg = 0; gg < data.
getN().size1(); gg++) {
315 double ambientTemp = pow(dAta.aMbienttEmp, 4.0);
318 if (ambientTemp > 0) {
319 tEmp = -ambientTemp + T4_at_Gauss_pt;
322 double val = getGaussPts()(2, gg);
323 double radiationConst;
326 double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
327 radiationConst = dAta.sIgma * dAta.eMissivity * tEmp * area;
329 radiationConst = dAta.sIgma * dAta.eMissivity * tEmp * getArea();
331 ublas::noalias(Nf) += val * radiationConst * data.
getN(gg, nb_row_dofs);
349 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
357 int rank = dof_ptr->getNbOfCoeffs();
359 int nb_row_dofs = data.
getIndices().size() / rank;
364 for (
unsigned int gg = 0; gg < data.
getN().size1(); gg++) {
367 double convectionConst;
369 double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
371 dAta.cOnvection * area * (T_at_Gauss_pt - dAta.tEmperature);
374 dAta.cOnvection * getArea() * (T_at_Gauss_pt - dAta.tEmperature);
376 double val = getGaussPts()(2, gg) * convectionConst;
377 ublas::noalias(Nf) += val * data.
getN(gg, nb_row_dofs);
397 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
406 int nb_row = row_data.
getN().size2();
407 int nb_col = col_data.
getN().size2();
408 K.resize(nb_row, nb_col);
411 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
413 double convectionConst;
415 double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
416 convectionConst = dAta.cOnvection * area;
418 convectionConst = dAta.cOnvection * getArea();
420 double val = getGaussPts()(2, gg) * convectionConst;
422 val * outer_prod(row_data.
getN(gg, nb_row), col_data.
getN(gg, nb_col));
426 const_cast<FEMethod *
>(getFEMethod())->ts_B =
A;
429 nb_col, &col_data.
getIndices()[0], &
K(0, 0), ADD_VALUES);
430 if (row_side != col_side || row_type != col_type) {
431 transK.resize(nb_col, nb_row);
432 noalias(transK) = trans(
K);
435 &row_data.
getIndices()[0], &transK(0, 0), ADD_VALUES);
444 problemPtr, tempName, rateName,
ROW, ts_u_t, INSERT_VALUES,
458 problemPtr,
ROW, ts_u, INSERT_VALUES, SCATTER_REVERSE);
460 BitRefLevel proble_bit_level = problemPtr->getBitRefLevel();
469 auto post_proc_at_points = [&](std::array<double, 3> point,
int num) {
472 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
474 struct OpPrint :
public VolOp {
476 std::array<double, 3> pointCoords;
478 boost::shared_ptr<VectorDouble> tempPtr;
480 OpPrint(boost::shared_ptr<VectorDouble> temp_ptr,
481 std::array<double, 3> &point_coords,
int point_num)
483 pointCoords(point_coords), pointNum(point_num) {}
486 DataForcesAndSourcesCore::EntData &data) {
488 if (type == MBVERTEX) {
489 if (getGaussPts().size2()) {
494 <<
"Pnt: " << std::to_string(pointNum)
503 if (
auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
504 fe_ptr->getOpPtrVector().push_back(
new OpPrint(tempPtr, point, num));
506 point.data(), 1e-12,
"DMTHERMAL",
"THERMAL_FE", dataFieldEval,
509 fe_ptr->getOpPtrVector().pop_back();
515 if (!evalPoints.empty()) {
517 for (
auto p : evalPoints)
518 CHKERR post_proc_at_points(
p, num++);
528 const std::string mesh_nodals_positions) {
537 mesh_nodals_positions);
547 ierr = it->getAttributeDataStructure(temp_data);
549 setOfBlocks[it->getMeshsetId()].cOnductivity_mat.resize(
551 setOfBlocks[it->getMeshsetId()].cOnductivity_mat.clear();
552 setOfBlocks[it->getMeshsetId()].cOnductivity_mat(0, 0) =
553 temp_data.
data.Conductivity;
554 setOfBlocks[it->getMeshsetId()].cOnductivity_mat(1, 1) =
555 temp_data.
data.Conductivity;
556 setOfBlocks[it->getMeshsetId()].cOnductivity_mat(2, 2) =
557 temp_data.
data.Conductivity;
560 setOfBlocks[it->getMeshsetId()].cApacity = temp_data.
data.HeatCapacity;
561 if (temp_data.
data.User2 != 0) {
565 it->meshset, MBTET,
setOfBlocks[it->getMeshsetId()].tEts,
true);
567 setOfBlocks[it->getMeshsetId()].tEts, MBTET,
"THERMAL_FE");
575 const std::string mesh_nodals_positions) {
587 mesh_nodals_positions);
594 it->meshset, MBTRI,
setOfFluxes[it->getMeshsetId()].tRis,
true);
596 setOfFluxes[it->getMeshsetId()].tRis, MBTRI,
"THERMAL_FLUX_FE");
602 if (std::regex_match(it->getName(), std::regex(
"(.*)HEAT_FLUX(.*)"))) {
603 std::vector<double> data;
604 CHKERR it->getAttributes(data);
605 if (data.size() != 1) {
606 SETERRQ(PETSC_COMM_SELF, 1,
"Data inconsistency");
608 strcpy(
setOfFluxes[it->getMeshsetId()].dAta.data.name,
"HeatFlu");
609 setOfFluxes[it->getMeshsetId()].dAta.data.flag1 = 1;
610 setOfFluxes[it->getMeshsetId()].dAta.data.value1 = data[0];
612 it->meshset, MBTRI,
setOfFluxes[it->getMeshsetId()].tRis,
true);
614 setOfFluxes[it->getMeshsetId()].tRis, MBTRI,
"THERMAL_FLUX_FE");
622 const std::string
field_name,
const std::string mesh_nodals_positions) {
634 mesh_nodals_positions);
640 if (std::regex_match(it->getName(), std::regex(
"(.*)CONVECTION(.*)"))) {
641 std::vector<double> data;
642 CHKERR it->getAttributes(data);
643 if (data.size() != 2) {
644 SETERRQ(PETSC_COMM_SELF, 1,
"Data inconsistency");
652 "THERMAL_CONVECTION_FE");
660 const std::string
field_name,
const std::string mesh_nodals_positions) {
672 mesh_nodals_positions);
678 if (std::regex_match(it->getName(), std::regex(
"(.*)RADIATION(.*)"))) {
679 std::vector<double> data;
680 ierr = it->getAttributes(data);
681 if (data.size() != 3) {
682 SETERRQ(PETSC_COMM_SELF, 1,
"Data inconsistency");
688 it->meshset, MBTRI,
setOfRadiation[it->getMeshsetId()].tRis,
true);
691 "THERMAL_RADIATION_FE");
701 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
715 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
725 string field_name, Vec &
F,
const std::string mesh_nodals_positions) {
727 bool hoGeometry =
false;
731 std::map<int, FluxData>::iterator sit =
setOfFluxes.begin();
741 string field_name, Vec &
F,
const std::string mesh_nodals_positions) {
743 bool hoGeometry =
false;
759 string field_name, Mat A,
const std::string mesh_nodals_positions) {
761 bool hoGeometry =
false;
776 const std::string mesh_nodals_positions) {
779 bool hoGeometry =
false;
785 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
811 std::map<int, FluxData>::iterator sit =
setOfFluxes.begin();
838 std::map<int, RadiationData>::iterator sit =
setOfRadiation.begin();
847 std::map<int, RadiationData>::iterator sit =
setOfRadiation.begin();
861 const std::string mesh_nodals_positions) {
870 loops_to_do_Rhs.push_back(
873 loops_to_do_Rhs.push_back(
876 loops_to_do_Rhs.push_back(
884 loops_to_do_Mat.push_back(
887 loops_to_do_Mat.push_back(
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Operators and data structures for thermal analysis.
#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 ...
@ MAT_THERMALSET
block name is "MAT_THERMAL"
#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 ...
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define MOFEM_LOG(channel, severity)
Log.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
virtual MoFEMErrorCode record_field(const std::string &serie_name, const std::string &field_name, const BitRefLevel &bit, const BitRefLevel &mask)
virtual MoFEMErrorCode record_begin(const std::string &serie_name)
virtual MoFEMErrorCode record_end(const std::string &serie_name, double time=0)
MoFEMErrorCode addThermalElements(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
add thermal element on tets
MoFEMErrorCode setThermalFluxFiniteElementRhsOperators(string field_name, Vec &F, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
this function is used in case of stationary problem for heat flux terms
MoFEMErrorCode setThermalFiniteElementRhsOperators(string field_name, Vec &F)
this function is used in case of stationary problem to set elements for rhs
MoFEMErrorCode setThermalFiniteElementLhsOperators(string field_name, Mat A)
this function is used in case of stationary heat conductivity problem for lhs
MoFEMErrorCode addThermalFluxElement(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
add heat flux element
MoFEMErrorCode addThermalConvectionElement(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
add convection element
MoFEMErrorCode addThermalRadiationElement(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
add Non-linear Radiation element
MoFEMErrorCode setTimeSteppingProblem(string field_name, string rate_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
set up operators for unsteady heat flux; convection; radiation problem
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr auto field_name
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure for User Loop Methods on finite elements
Field evaluator interface.
@ OPROW
operator doWork function is executed on FE rows
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Thermal material data structure.
Interface for Time Stepping (TS) solver.
MoFEM::FEMethodsSequence FEMethodsSequence
FEMethodsSequence & getLoopsIFunction()
Get the loops to do IFunction object.
FEMethodsSequence & getLoopsIJacobian()
Get the loops to do IJacobian object.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
VectorDouble temperatureAtGaussPts
VectorDouble temperatureRateAtGaussPts
ublas::matrix_row< MatrixDouble > getGradAtGaussPts(const int gg)
MatrixDouble gradAtGaussPts
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
calculate thermal convection term in the lhs of equations
operator to calculate convection therms on body surface and assemble to rhs of equations
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator to calculate temperature gradient at Gauss points
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating temperature gradients
operator to calculate temperature rate at Gauss pts
operator to calculate temperature at Gauss pts
operator to calculate temperature at Gauss pts
operator to calculate left hand side of heat capacity terms
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
calculate heat capacity matrix
operator to calculate right hand side of heat capacity terms
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate thermal conductivity matrix
operator for calculate heat flux and assemble to right hand side
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate heat flux
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
calculate thermal radiation term in the lhs of equations(Tangent Matrix) for transient Thermal Proble...
operator to calculate radiation therms on body surface and assemble to rhs of transient equations(Res...
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate Transient Radiation condition on the right hand side residual
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
calculate thermal conductivity matrix
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate thermal conductivity matrix
VolumeElementForcesAndSourcesCore::UserDataOperator VolOp
MoFEMErrorCode postProcess()
MoFEMErrorCode preProcess()
MoFEMErrorCode postProcess()
std::map< int, FluxData > setOfFluxes
maps side set id with appropriate FluxData
MyVolumeFE feRhs
cauclate right hand side for tetrahedral elements
MoFEM::Interface & mField
MoFEMErrorCode setThermalConvectionFiniteElementLhsOperators(string field_name, Mat A, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
std::map< int, RadiationData > setOfRadiation
std::map< int, ConvectionData > setOfConvection
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MoFEMErrorCode setThermalConvectionFiniteElementRhsOperators(string field_name, Vec &F, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
VolEle::UserDataOperator VolOp