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),
75 int row_side,
int col_side, EntityType row_type, EntityType col_type,
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);
100 MatrixDouble K1 = prod(row_data.
getDiffN(gg, nb_row), val);
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;
152 int row_side,
int col_side, EntityType row_type, EntityType col_type,
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_NULLPTR) {
241 int row_side,
int col_side, EntityType row_type, EntityType col_type,
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);
392 int row_side,
int col_side, EntityType row_type, EntityType col_type,
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) {}
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 ->evalFEAtThePoint<3>(point.data(), 1e-12,
"DMTHERMAL",
"THERMAL_FE",
510 fe_ptr->getOpPtrVector().pop_back();
516 if (!evalPoints.empty()) {
518 for (
auto p : evalPoints)
519 CHKERR post_proc_at_points(p, num++);
529 const std::string mesh_nodals_positions) {
538 mesh_nodals_positions);
548 ierr = it->getAttributeDataStructure(temp_data);
550 setOfBlocks[it->getMeshsetId()].cOnductivity_mat.resize(
552 setOfBlocks[it->getMeshsetId()].cOnductivity_mat.clear();
553 setOfBlocks[it->getMeshsetId()].cOnductivity_mat(0, 0) =
554 temp_data.
data.Conductivity;
555 setOfBlocks[it->getMeshsetId()].cOnductivity_mat(1, 1) =
556 temp_data.
data.Conductivity;
557 setOfBlocks[it->getMeshsetId()].cOnductivity_mat(2, 2) =
558 temp_data.
data.Conductivity;
561 setOfBlocks[it->getMeshsetId()].cApacity = temp_data.
data.HeatCapacity;
562 if (temp_data.
data.User2 != 0) {
566 it->meshset, MBTET,
setOfBlocks[it->getMeshsetId()].tEts,
true);
568 setOfBlocks[it->getMeshsetId()].tEts, MBTET,
"THERMAL_FE");
576 const std::string mesh_nodals_positions) {
588 mesh_nodals_positions);
595 it->meshset, MBTRI,
setOfFluxes[it->getMeshsetId()].tRis,
true);
597 setOfFluxes[it->getMeshsetId()].tRis, MBTRI,
"THERMAL_FLUX_FE");
603 if (std::regex_match(it->getName(), std::regex(
"(.*)HEAT_FLUX(.*)"))) {
604 std::vector<double> data;
605 CHKERR it->getAttributes(data);
606 if (data.size() != 1) {
607 SETERRQ(PETSC_COMM_SELF, 1,
"Data inconsistency");
609 strcpy(
setOfFluxes[it->getMeshsetId()].dAta.data.name,
"HeatFlu");
610 setOfFluxes[it->getMeshsetId()].dAta.data.flag1 = 1;
611 setOfFluxes[it->getMeshsetId()].dAta.data.value1 = data[0];
613 it->meshset, MBTRI,
setOfFluxes[it->getMeshsetId()].tRis,
true);
615 setOfFluxes[it->getMeshsetId()].tRis, MBTRI,
"THERMAL_FLUX_FE");
623 const std::string
field_name,
const std::string mesh_nodals_positions) {
635 mesh_nodals_positions);
641 if (std::regex_match(it->getName(), std::regex(
"(.*)CONVECTION(.*)"))) {
642 std::vector<double> data;
643 CHKERR it->getAttributes(data);
644 if (data.size() != 2) {
645 SETERRQ(PETSC_COMM_SELF, 1,
"Data inconsistency");
653 "THERMAL_CONVECTION_FE");
661 const std::string
field_name,
const std::string mesh_nodals_positions) {
673 mesh_nodals_positions);
679 if (std::regex_match(it->getName(), std::regex(
"(.*)RADIATION(.*)"))) {
680 std::vector<double> data;
681 ierr = it->getAttributes(data);
682 if (data.size() != 3) {
683 SETERRQ(PETSC_COMM_SELF, 1,
"Data inconsistency");
689 it->meshset, MBTRI,
setOfRadiation[it->getMeshsetId()].tRis,
true);
692 "THERMAL_RADIATION_FE");
702 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
716 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
726 string field_name, Vec &
F,
const std::string mesh_nodals_positions) {
728 bool hoGeometry =
false;
732 std::map<int, FluxData>::iterator sit =
setOfFluxes.begin();
742 string field_name, Vec &
F,
const std::string mesh_nodals_positions) {
744 bool hoGeometry =
false;
760 string field_name, Mat A,
const std::string mesh_nodals_positions) {
762 bool hoGeometry =
false;
777 const std::string mesh_nodals_positions) {
780 bool hoGeometry =
false;
786 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
812 std::map<int, FluxData>::iterator sit =
setOfFluxes.begin();
839 std::map<int, RadiationData>::iterator sit =
setOfRadiation.begin();
848 std::map<int, RadiationData>::iterator sit =
setOfRadiation.begin();
862 const std::string mesh_nodals_positions) {
871 loops_to_do_Rhs.push_back(
874 loops_to_do_Rhs.push_back(
877 loops_to_do_Rhs.push_back(
885 loops_to_do_Mat.push_back(
888 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_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
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 reference 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(TangentMatrix) for transient Thermal Problem
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