v0.15.0
Loading...
Searching...
No Matches
ThermalElement.cpp
Go to the documentation of this file.
1/** \file ThermalElement.cpp
2 \ingroup mofem_thermal_elem
3*/
4
5
6
7#include <MoFEM.hpp>
8using namespace MoFEM;
9#include <ThermalElement.hpp>
10
11using namespace boost::numeric;
12
14 int side, EntityType type, EntitiesFieldData::EntData &data) {
16
17 if (data.getIndices().size() == 0)
19 int nb_dofs = data.getFieldData().size();
20 int nb_gauss_pts = data.getN().size1();
21
22 // initialize
23 commonData.gradAtGaussPts.resize(nb_gauss_pts, 3);
24 if (type == MBVERTEX) {
25 std::fill(commonData.gradAtGaussPts.data().begin(),
26 commonData.gradAtGaussPts.data().end(), 0);
27 }
28
29 for (int gg = 0; gg < nb_gauss_pts; gg++) {
30 ublas::noalias(commonData.getGradAtGaussPts(gg)) +=
31 prod(trans(data.getDiffN(gg, nb_dofs)), data.getFieldData());
32 }
33
35}
36
38ThermalElement::OpThermalRhs::doWork(int side, EntityType type,
41
42 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
43 dAta.tEts.end()) {
45 }
46 if (data.getIndices().size() == 0)
48
49 int nb_row_dofs = data.getIndices().size();
50 Nf.resize(nb_row_dofs);
51 Nf.clear();
52
53 for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
54
55 MatrixDouble val =
56 dAta.cOnductivity_mat * getVolume() * getGaussPts()(3, gg);
57
58 // ublas
59 ublas::noalias(Nf) += prod(prod(data.getDiffN(gg, nb_row_dofs), val),
61 }
62
63 if (useTsF) {
64 CHKERR VecSetValues(getFEMethod()->ts_F, data.getIndices().size(),
65 &data.getIndices()[0], &Nf[0], ADD_VALUES);
66 } else {
67 CHKERR VecSetValues(F, data.getIndices().size(), &data.getIndices()[0],
68 &Nf[0], ADD_VALUES);
69 }
70
72}
73
75 int row_side, int col_side, EntityType row_type, EntityType col_type,
79
80 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
81 dAta.tEts.end()) {
83 }
84
85 if (row_data.getIndices().size() == 0)
87 if (col_data.getIndices().size() == 0)
89
90 int nb_row = row_data.getN().size2();
91 int nb_col = col_data.getN().size2();
92 K.resize(nb_row, nb_col);
93 K.clear();
94 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
95
96 MatrixDouble val =
97 dAta.cOnductivity_mat * getVolume() * getGaussPts()(3, gg);
98
99 // ublas
100 MatrixDouble K1 = prod(row_data.getDiffN(gg, nb_row), val);
101 noalias(K) += prod(K1, trans(col_data.getDiffN(gg, nb_col)));
102 }
103
104 if (!useTsB) {
105 const_cast<FEMethod *>(getFEMethod())->ts_B = A;
106 }
107 CHKERR MatSetValues((getFEMethod()->ts_B), nb_row, &row_data.getIndices()[0],
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);
112 CHKERR MatSetValues((getFEMethod()->ts_B), nb_col,
113 &col_data.getIndices()[0], nb_row,
114 &row_data.getIndices()[0], &transK(0, 0), ADD_VALUES);
115 }
116
118}
119
121 int side, EntityType type, EntitiesFieldData::EntData &data) {
123
124 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
125 dAta.tEts.end()) {
127 }
128 if (data.getIndices().size() == 0)
130
131 int nb_row = data.getN().size2();
132 Nf.resize(nb_row);
133 Nf.clear();
134 for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
135 double val = getGaussPts()(3, gg);
137 ////////////
138 // cblas
139 // cblas_daxpy(nb_row,val,&data.getN()(gg,0),1,&*Nf.data().begin(),1);
140 // ublas
141 ublas::noalias(Nf) += val * data.getN(gg);
142 }
143 Nf *= getVolume() * dAta.cApacity;
144
145 CHKERR VecSetValues(getFEMethod()->ts_F, data.getIndices().size(),
146 &data.getIndices()[0], &Nf[0], ADD_VALUES);
147
149}
150
152 int row_side, int col_side, EntityType row_type, EntityType col_type,
154 EntitiesFieldData::EntData &col_data) {
156
157 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
158 dAta.tEts.end()) {
160 }
161
162 if (row_data.getIndices().size() == 0)
164 if (col_data.getIndices().size() == 0)
166
167 int nb_row = row_data.getN().size2();
168 int nb_col = col_data.getN().size2();
169 M.resize(nb_row, nb_col);
170 M.clear();
171
172 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
173
174 double val = getGaussPts()(3, gg);
175
176 // ublas
177 noalias(M) +=
178 val * outer_prod(row_data.getN(gg, nb_row), col_data.getN(gg, nb_col));
179 }
180
181 M *= getVolume() * dAta.cApacity * getFEMethod()->ts_a;
182
183 CHKERR MatSetValues((getFEMethod()->ts_B), nb_row, &row_data.getIndices()[0],
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);
188 CHKERR MatSetValues((getFEMethod()->ts_B), nb_col,
189 &col_data.getIndices()[0], nb_row,
190 &row_data.getIndices()[0], &transM(0, 0), ADD_VALUES);
191 }
192
194}
195
197ThermalElement::OpHeatFlux::doWork(int side, EntityType type,
200
201 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
202 dAta.tRis.end()) {
204 }
205 if (data.getIndices().size() == 0)
207
208 const auto &dof_ptr = data.getFieldDofs()[0];
209 int rank = dof_ptr->getNbOfCoeffs();
210
211 int nb_dofs = data.getIndices().size() / rank;
212
213 Nf.resize(data.getIndices().size());
214 Nf.clear();
215
216 for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
217
218 double val = getGaussPts()(2, gg);
219 double flux;
220 if (hoGeometry) {
221 const double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
222 flux = dAta.dAta.data.value1 * area;
223 } else {
224 flux = dAta.dAta.data.value1 * getArea();
225 }
226 ublas::noalias(Nf) += val * flux * data.getN(gg, nb_dofs);
227 }
228
229 if (useTsF || F == PETSC_NULLPTR) {
230 CHKERR VecSetValues(getFEMethod()->ts_F, data.getIndices().size(),
231 &data.getIndices()[0], &Nf[0], ADD_VALUES);
232 } else {
233 CHKERR VecSetValues(F, data.getIndices().size(), &data.getIndices()[0],
234 &Nf[0], ADD_VALUES);
235 }
236
238}
239
241 int row_side, int col_side, EntityType row_type, EntityType col_type,
243 EntitiesFieldData::EntData &col_data) {
245
246 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
247 dAta.tRis.end()) {
249 }
250
251 if (row_data.getIndices().size() == 0)
253 if (col_data.getIndices().size() == 0)
255
256 int nb_row = row_data.getN().size2();
257 int nb_col = col_data.getN().size2();
258
259 N.resize(nb_row, nb_col);
260 N.clear();
261
262 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
263 double T3_at_Gauss_pt = pow(commonData.temperatureAtGaussPts[gg], 3.0);
264
265 double radiationConst;
266 if (hoGeometry) {
267 double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
268 radiationConst = dAta.sIgma * dAta.eMissivity * area;
269 } else {
270 radiationConst = dAta.sIgma * dAta.eMissivity * getArea();
271 }
272 const double fOur = 4.0;
273 double val = fOur * getGaussPts()(2, gg) * radiationConst * T3_at_Gauss_pt;
274 noalias(N) +=
275 val * outer_prod(row_data.getN(gg, nb_row), col_data.getN(gg, nb_col));
276 }
277
278 if (!useTsB) {
279 const_cast<FEMethod *>(getFEMethod())->ts_B = A;
280 }
281 CHKERR MatSetValues((getFEMethod()->ts_B), nb_row, &row_data.getIndices()[0],
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);
286 CHKERR MatSetValues((getFEMethod()->ts_B), nb_col,
287 &col_data.getIndices()[0], nb_row,
288 &row_data.getIndices()[0], &transN(0, 0), ADD_VALUES);
289 }
290
292}
293
295 int side, EntityType type, EntitiesFieldData::EntData &data) {
297
298 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
299 dAta.tRis.end()) {
301 }
302 if (data.getIndices().size() == 0)
304
305 const auto &dof_ptr = data.getFieldDofs()[0];
306 int rank = dof_ptr->getNbOfCoeffs();
307 int nb_row_dofs = data.getIndices().size() / rank;
308
309 Nf.resize(data.getIndices().size());
310 Nf.clear();
311
312 for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
313
314 double T4_at_Gauss_pt = pow(commonData.temperatureAtGaussPts[gg], 4.0);
315 double ambientTemp = pow(dAta.aMbienttEmp, 4.0);
316 double tEmp = 0;
317
318 if (ambientTemp > 0) {
319 tEmp = -ambientTemp + T4_at_Gauss_pt;
320 }
321
322 double val = getGaussPts()(2, gg);
323 double radiationConst;
324
325 if (hoGeometry) {
326 double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
327 radiationConst = dAta.sIgma * dAta.eMissivity * tEmp * area;
328 } else {
329 radiationConst = dAta.sIgma * dAta.eMissivity * tEmp * getArea();
330 }
331 ublas::noalias(Nf) += val * radiationConst * data.getN(gg, nb_row_dofs);
332 }
333
334 if (useTsF) {
335 CHKERR VecSetValues(getFEMethod()->ts_F, data.getIndices().size(),
336 &data.getIndices()[0], &Nf[0], ADD_VALUES);
337 } else {
338 CHKERR VecSetValues(F, data.getIndices().size(), &data.getIndices()[0],
339 &Nf[0], ADD_VALUES);
340 }
341
343}
344
346 int side, EntityType type, EntitiesFieldData::EntData &data) {
348
349 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
350 dAta.tRis.end()) {
352 }
353 if (data.getIndices().size() == 0)
355
356 const auto &dof_ptr = data.getFieldDofs()[0];
357 int rank = dof_ptr->getNbOfCoeffs();
358
359 int nb_row_dofs = data.getIndices().size() / rank;
360
361 Nf.resize(data.getIndices().size());
362 Nf.clear();
363
364 for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
365
366 double T_at_Gauss_pt = commonData.temperatureAtGaussPts[gg];
367 double convectionConst;
368 if (hoGeometry) {
369 double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
370 convectionConst =
371 dAta.cOnvection * area * (T_at_Gauss_pt - dAta.tEmperature);
372 } else {
373 convectionConst =
374 dAta.cOnvection * getArea() * (T_at_Gauss_pt - dAta.tEmperature);
375 }
376 double val = getGaussPts()(2, gg) * convectionConst;
377 ublas::noalias(Nf) += val * data.getN(gg, nb_row_dofs);
378 }
379
380 if (useTsF) {
381 CHKERR VecSetValues(getFEMethod()->ts_F, data.getIndices().size(),
382 &data.getIndices()[0], &Nf[0], ADD_VALUES);
383 } else {
384 CHKERR VecSetValues(F, data.getIndices().size(), &data.getIndices()[0],
385 &Nf[0], ADD_VALUES);
386 }
387
389}
390
392 int row_side, int col_side, EntityType row_type, EntityType col_type,
394 EntitiesFieldData::EntData &col_data) {
396
397 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
398 dAta.tRis.end()) {
400 }
401 if (row_data.getIndices().size() == 0)
403 if (col_data.getIndices().size() == 0)
405
406 int nb_row = row_data.getN().size2();
407 int nb_col = col_data.getN().size2();
408 K.resize(nb_row, nb_col);
409 K.clear();
410
411 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
412
413 double convectionConst;
414 if (hoGeometry) {
415 double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
416 convectionConst = dAta.cOnvection * area;
417 } else {
418 convectionConst = dAta.cOnvection * getArea();
419 }
420 double val = getGaussPts()(2, gg) * convectionConst;
421 noalias(K) +=
422 val * outer_prod(row_data.getN(gg, nb_row), col_data.getN(gg, nb_col));
423 }
424
425 if (!useTsB) {
426 const_cast<FEMethod *>(getFEMethod())->ts_B = A;
427 }
428 CHKERR MatSetValues((getFEMethod()->ts_B), nb_row, &row_data.getIndices()[0],
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);
433 CHKERR MatSetValues((getFEMethod()->ts_B), nb_col,
434 &col_data.getIndices()[0], nb_row,
435 &row_data.getIndices()[0], &transK(0, 0), ADD_VALUES);
436 }
437
439}
440
443 CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
444 problemPtr, tempName, rateName, ROW, ts_u_t, INSERT_VALUES,
445 SCATTER_REVERSE);
447}
448
453
456
457 CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
458 problemPtr, ROW, ts_u, INSERT_VALUES, SCATTER_REVERSE);
459
460 BitRefLevel proble_bit_level = problemPtr->getBitRefLevel();
461
462 SeriesRecorder *recorder_ptr = NULL;
463 CHKERR mField.getInterface(recorder_ptr);
464 CHKERR recorder_ptr->record_begin(seriesName);
465 CHKERR recorder_ptr->record_field(seriesName, tempName, proble_bit_level,
466 mask);
467 CHKERR recorder_ptr->record_end(seriesName, ts_t);
468
469 auto post_proc_at_points = [&](std::array<double, 3> point, int num) {
471
472 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
473
474 struct OpPrint : public VolOp {
475
476 std::array<double, 3> pointCoords;
477 int pointNum;
478 boost::shared_ptr<VectorDouble> tempPtr;
479
480 OpPrint(boost::shared_ptr<VectorDouble> temp_ptr,
481 std::array<double, 3> &point_coords, int point_num)
482 : VolOp("TEMP", VolOp::OPROW), tempPtr(temp_ptr),
483 pointCoords(point_coords), pointNum(point_num) {}
484
485 MoFEMErrorCode doWork(int side, EntityType type,
488 if (type == MBVERTEX) {
489 if (getGaussPts().size2()) {
490
491 auto t_p = getFTensor0FromVec(*tempPtr);
492
493 MOFEM_LOG("THERMALSYNC", Sev::inform)
494 << "Pnt: " << std::to_string(pointNum)
495 << " Crd: " << getVectorAdaptor(pointCoords.data(), 3)
496 << " Tmp: " << t_p;
497 }
498 }
500 }
501 };
502
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",
507 dataFieldEval, mField.get_comm_rank(),
508 mField.get_comm_rank(), nullptr, MF_EXIST,
509 QUIET);
510 fe_ptr->getOpPtrVector().pop_back();
511 }
512
514 };
515
516 if (!evalPoints.empty()) {
517 int num = 0;
518 for (auto p : evalPoints)
519 CHKERR post_proc_at_points(p, num++);
520 }
521
523
525}
526
529 const std::string mesh_nodals_positions) {
531
536 if (mField.check_field(mesh_nodals_positions)) {
538 mesh_nodals_positions);
539 }
540
541 // takes skin of block of entities
542 // Skinner skin(&mField.get_moab());
543 // loop over all blocksets and get data which name is FluidPressure
546
547 Mat_Thermal temp_data;
548 ierr = it->getAttributeDataStructure(temp_data);
549
550 setOfBlocks[it->getMeshsetId()].cOnductivity_mat.resize(
551 3, 3); //(3X3) conductivity matrix
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;
559 // setOfBlocks[it->getMeshsetId()].cOnductivity =
560 // temp_data.data.Conductivity;
561 setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
562 if (temp_data.data.User2 != 0) {
563 setOfBlocks[it->getMeshsetId()].initTemp = temp_data.data.User2;
564 }
565 CHKERR mField.get_moab().get_entities_by_type(
566 it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
568 setOfBlocks[it->getMeshsetId()].tEts, MBTET, "THERMAL_FE");
569 }
570
572}
573
576 const std::string mesh_nodals_positions) {
578
579 CHKERR mField.add_finite_element("THERMAL_FLUX_FE", MF_ZERO);
581 field_name);
583 field_name);
585 field_name);
586 if (mField.check_field(mesh_nodals_positions)) {
588 mesh_nodals_positions);
589 }
590
592 it)) {
593 CHKERR it->getBcDataStructure(setOfFluxes[it->getMeshsetId()].dAta);
594 CHKERR mField.get_moab().get_entities_by_type(
595 it->meshset, MBTRI, setOfFluxes[it->getMeshsetId()].tRis, true);
597 setOfFluxes[it->getMeshsetId()].tRis, MBTRI, "THERMAL_FLUX_FE");
598 }
599
600 // this is alternative method for setting boundary conditions, to bypass bu
601 // in cubit file reader. not elegant, but good enough
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");
608 }
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];
612 CHKERR mField.get_moab().get_entities_by_type(
613 it->meshset, MBTRI, setOfFluxes[it->getMeshsetId()].tRis, true);
615 setOfFluxes[it->getMeshsetId()].tRis, MBTRI, "THERMAL_FLUX_FE");
616 }
617 }
618
620}
621
623 const std::string field_name, const std::string mesh_nodals_positions) {
625
626 CHKERR mField.add_finite_element("THERMAL_CONVECTION_FE", MF_ZERO);
627 CHKERR mField.modify_finite_element_add_field_row("THERMAL_CONVECTION_FE",
628 field_name);
629 CHKERR mField.modify_finite_element_add_field_col("THERMAL_CONVECTION_FE",
630 field_name);
631 CHKERR mField.modify_finite_element_add_field_data("THERMAL_CONVECTION_FE",
632 field_name);
633 if (mField.check_field(mesh_nodals_positions)) {
634 CHKERR mField.modify_finite_element_add_field_data("THERMAL_CONVECTION_FE",
635 mesh_nodals_positions);
636 }
637
638 // this is alternative method for setting boundary conditions, to bypass bu
639 // in cubit file reader. not elegant, but good enough
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");
646 }
647 setOfConvection[it->getMeshsetId()].cOnvection = data[0];
648 setOfConvection[it->getMeshsetId()].tEmperature = data[1];
649 CHKERR mField.get_moab().get_entities_by_type(
650 it->meshset, MBTRI, setOfConvection[it->getMeshsetId()].tRis, true);
652 setOfConvection[it->getMeshsetId()].tRis, MBTRI,
653 "THERMAL_CONVECTION_FE");
654 }
655 }
656
658}
659
661 const std::string field_name, const std::string mesh_nodals_positions) {
663
664 CHKERR mField.add_finite_element("THERMAL_RADIATION_FE", MF_ZERO);
666 field_name);
668 field_name);
670 field_name);
671 if (mField.check_field(mesh_nodals_positions)) {
673 mesh_nodals_positions);
674 }
675
676 // this is alternative method for setting boundary conditions, to bypass bu
677 // in cubit file reader. not elegant, but good enough
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");
684 }
685 setOfRadiation[it->getMeshsetId()].sIgma = data[0];
686 setOfRadiation[it->getMeshsetId()].eMissivity = data[1];
687 setOfRadiation[it->getMeshsetId()].aMbienttEmp = data[2];
688 CHKERR mField.get_moab().get_entities_by_type(
689 it->meshset, MBTRI, setOfRadiation[it->getMeshsetId()].tRis, true);
691 setOfRadiation[it->getMeshsetId()].tRis, MBTRI,
692 "THERMAL_RADIATION_FE");
693 }
694 }
695
697}
698
702 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
703 feRhs.getOpPtrVector().push_back(
705 for (; sit != setOfBlocks.end(); sit++) {
706 // add finite element
707 feRhs.getOpPtrVector().push_back(
708 new OpThermalRhs(field_name, F, sit->second, commonData));
709 }
711}
712
716 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
717 for (; sit != setOfBlocks.end(); sit++) {
718 // add finite elemen
719 feLhs.getOpPtrVector().push_back(
720 new OpThermalLhs(field_name, A, sit->second, commonData));
721 }
723}
724
726 string field_name, Vec &F, const std::string mesh_nodals_positions) {
728 bool hoGeometry = false;
729 if (mField.check_field(mesh_nodals_positions)) {
730 hoGeometry = true;
731 }
732 std::map<int, FluxData>::iterator sit = setOfFluxes.begin();
733 for (; sit != setOfFluxes.end(); sit++) {
734 // add finite element
735 feFlux.getOpPtrVector().push_back(
736 new OpHeatFlux(field_name, F, sit->second, hoGeometry));
737 }
739}
740
742 string field_name, Vec &F, const std::string mesh_nodals_positions) {
744 bool hoGeometry = false;
745 if (mField.check_field(mesh_nodals_positions)) {
746 hoGeometry = true;
747 }
748 std::map<int, ConvectionData>::iterator sit = setOfConvection.begin();
749 for (; sit != setOfConvection.end(); sit++) {
750 // add finite element
754 field_name, F, sit->second, commonData, hoGeometry));
755 }
757}
758
760 string field_name, Mat A, const std::string mesh_nodals_positions) {
762 bool hoGeometry = false;
763 if (mField.check_field(mesh_nodals_positions)) {
764 hoGeometry = true;
765 }
766 std::map<int, ConvectionData>::iterator sit = setOfConvection.begin();
767 for (; sit != setOfConvection.end(); sit++) {
768 // add finite element
770 new OpConvectionLhs(field_name, A, sit->second, hoGeometry));
771 }
773}
774
776 string field_name, string rate_name,
777 const std::string mesh_nodals_positions) {
779
780 bool hoGeometry = false;
781 if (mField.check_field(mesh_nodals_positions)) {
782 hoGeometry = true;
783 }
784
785 {
786 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
787 for (; sit != setOfBlocks.end(); sit++) {
788 // add finite element
789 // those methods are to calculate matrices on Lhs
790 // feLhs.getOpPtrVector().push_back(new
791 // OpGetTetTemperatureAtGaussPts(field_name,commonData));
792 feLhs.getOpPtrVector().push_back(
793 new OpThermalLhs(field_name, sit->second, commonData));
794 feLhs.getOpPtrVector().push_back(
795 new OpHeatCapacityLhs(field_name, sit->second, commonData));
796 // those methods are to calculate vectors on Rhs
797 feRhs.getOpPtrVector().push_back(
799 feRhs.getOpPtrVector().push_back(
800 new OpGetTetRateAtGaussPts(rate_name, commonData));
801 feRhs.getOpPtrVector().push_back(
803 feRhs.getOpPtrVector().push_back(
804 new OpThermalRhs(field_name, sit->second, commonData));
805 feRhs.getOpPtrVector().push_back(
806 new OpHeatCapacityRhs(field_name, sit->second, commonData));
807 }
808 }
809
810 // Flux
811 {
812 std::map<int, FluxData>::iterator sit = setOfFluxes.begin();
813 for (; sit != setOfFluxes.end(); sit++) {
814 feFlux.getOpPtrVector().push_back(
815 new OpHeatFlux(field_name, sit->second, hoGeometry));
816 }
817 }
818
819 // Convection
820 {
821 std::map<int, ConvectionData>::iterator sit = setOfConvection.begin();
822 for (; sit != setOfConvection.end(); sit++) {
826 new OpConvectionRhs(field_name, sit->second, commonData, hoGeometry));
827 }
828 }
829 {
830 std::map<int, ConvectionData>::iterator sit = setOfConvection.begin();
831 for (; sit != setOfConvection.end(); sit++) {
833 new OpConvectionLhs(field_name, sit->second, hoGeometry));
834 }
835 }
836
837 // Radiation
838 {
839 std::map<int, RadiationData>::iterator sit = setOfRadiation.begin();
840 for (; sit != setOfRadiation.end(); sit++) {
841 feRadiationRhs.getOpPtrVector().push_back(
843 feRadiationRhs.getOpPtrVector().push_back(
844 new OpRadiationRhs(field_name, sit->second, commonData, hoGeometry));
845 }
846 }
847 {
848 std::map<int, RadiationData>::iterator sit = setOfRadiation.begin();
849 for (; sit != setOfRadiation.end(); sit++) {
850 feRadiationLhs.getOpPtrVector().push_back(
852 feRadiationLhs.getOpPtrVector().push_back(
853 new OpRadiationLhs(field_name, sit->second, commonData, hoGeometry));
854 }
855 }
856
858}
859
861 TsCtx &ts_ctx, string field_name, string rate_name,
862 const std::string mesh_nodals_positions) {
864
865 CHKERR setTimeSteppingProblem(field_name, rate_name, mesh_nodals_positions);
866
867 // rhs
868 TsCtx::FEMethodsSequence &loops_to_do_Rhs =
870 loops_to_do_Rhs.push_back(TsCtx::PairNameFEMethodPtr("THERMAL_FE", &feRhs));
871 loops_to_do_Rhs.push_back(
872 TsCtx::PairNameFEMethodPtr("THERMAL_FLUX_FE", &feFlux));
873 if (mField.check_finite_element("THERMAL_CONVECTION_FE"))
874 loops_to_do_Rhs.push_back(
875 TsCtx::PairNameFEMethodPtr("THERMAL_CONVECTION_FE", &feConvectionRhs));
876 if (mField.check_finite_element("THERMAL_RADIATION_FE"))
877 loops_to_do_Rhs.push_back(
878 TsCtx::PairNameFEMethodPtr("THERMAL_RADIATION_FE", &feRadiationRhs));
879
880 // lhs
881 TsCtx::FEMethodsSequence &loops_to_do_Mat =
883 loops_to_do_Mat.push_back(TsCtx::PairNameFEMethodPtr("THERMAL_FE", &feLhs));
884 if (mField.check_finite_element("THERMAL_CONVECTION_FE"))
885 loops_to_do_Mat.push_back(
886 TsCtx::PairNameFEMethodPtr("THERMAL_CONVECTION_FE", &feConvectionLhs));
887 if (mField.check_finite_element("THERMAL_RADIATION_FE"))
888 loops_to_do_Mat.push_back(
889 TsCtx::PairNameFEMethodPtr("THERMAL_RADIATION_FE", &feRadiationLhs));
890
892}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Operators and data structures for thermal analysis.
@ QUIET
@ ROW
@ MF_ZERO
@ MF_EXIST
#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 ...
@ HEATFLUXSET
@ SIDESET
@ MAT_THERMALSET
block name is "MAT_THERMAL"
@ BLOCKSET
#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 ...
@ F
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
MoFEM::TsCtx * ts_ctx
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.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition Templates.hpp:31
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
const int N
Definition speed_test.cpp:3
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.
Definition TsCtx.hpp:17
MoFEM::FEMethodsSequence FEMethodsSequence
Definition TsCtx.hpp:26
FEMethodsSequence & getLoopsIFunction()
Get the loops to do IFunction object.
Definition TsCtx.hpp:63
FEMethodsSequence & getLoopsIJacobian()
Get the loops to do IJacobian object.
Definition TsCtx.hpp:83
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
ublas::matrix_row< MatrixDouble > getGradAtGaussPts(const int gg)
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
std::map< int, FluxData > setOfFluxes
maps side set id with appropriate FluxData
MyVolumeFE feRhs
cauclate right hand side for tetrahedral elements
MoFEM::Interface & mField
CommonData commonData
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