v0.13.2
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
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
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_NULL) {
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
452}
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,
486 DataForcesAndSourcesCore::EntData &data) {
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));
505 CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint3D(
506 point.data(), 1e-12, "DMTHERMAL", "THERMAL_FE", dataFieldEval,
508 QUIET);
509 fe_ptr->getOpPtrVector().pop_back();
510 }
511
513 };
514
515 if (!evalPoints.empty()) {
516 int num = 0;
517 for (auto p : evalPoints)
518 CHKERR post_proc_at_points(p, num++);
519 }
520
522
524}
525
528 const std::string mesh_nodals_positions) {
530
535 if (mField.check_field(mesh_nodals_positions)) {
537 mesh_nodals_positions);
538 }
539
540 // takes skin of block of entities
541 // Skinner skin(&mField.get_moab());
542 // loop over all blocksets and get data which name is FluidPressure
545
546 Mat_Thermal temp_data;
547 ierr = it->getAttributeDataStructure(temp_data);
548
549 setOfBlocks[it->getMeshsetId()].cOnductivity_mat.resize(
550 3, 3); //(3X3) conductivity matrix
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;
558 // setOfBlocks[it->getMeshsetId()].cOnductivity =
559 // temp_data.data.Conductivity;
560 setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
561 if (temp_data.data.User2 != 0) {
562 setOfBlocks[it->getMeshsetId()].initTemp = temp_data.data.User2;
563 }
564 CHKERR mField.get_moab().get_entities_by_type(
565 it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
567 setOfBlocks[it->getMeshsetId()].tEts, MBTET, "THERMAL_FE");
568 }
569
571}
572
575 const std::string mesh_nodals_positions) {
577
578 CHKERR mField.add_finite_element("THERMAL_FLUX_FE", MF_ZERO);
580 field_name);
582 field_name);
584 field_name);
585 if (mField.check_field(mesh_nodals_positions)) {
587 mesh_nodals_positions);
588 }
589
591 it)) {
592 CHKERR it->getBcDataStructure(setOfFluxes[it->getMeshsetId()].dAta);
593 CHKERR mField.get_moab().get_entities_by_type(
594 it->meshset, MBTRI, setOfFluxes[it->getMeshsetId()].tRis, true);
596 setOfFluxes[it->getMeshsetId()].tRis, MBTRI, "THERMAL_FLUX_FE");
597 }
598
599 // this is alternative method for setting boundary conditions, to bypass bu
600 // in cubit file reader. not elegant, but good enough
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");
607 }
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];
611 CHKERR mField.get_moab().get_entities_by_type(
612 it->meshset, MBTRI, setOfFluxes[it->getMeshsetId()].tRis, true);
614 setOfFluxes[it->getMeshsetId()].tRis, MBTRI, "THERMAL_FLUX_FE");
615 }
616 }
617
619}
620
622 const std::string field_name, const std::string mesh_nodals_positions) {
624
625 CHKERR mField.add_finite_element("THERMAL_CONVECTION_FE", MF_ZERO);
626 CHKERR mField.modify_finite_element_add_field_row("THERMAL_CONVECTION_FE",
627 field_name);
628 CHKERR mField.modify_finite_element_add_field_col("THERMAL_CONVECTION_FE",
629 field_name);
630 CHKERR mField.modify_finite_element_add_field_data("THERMAL_CONVECTION_FE",
631 field_name);
632 if (mField.check_field(mesh_nodals_positions)) {
633 CHKERR mField.modify_finite_element_add_field_data("THERMAL_CONVECTION_FE",
634 mesh_nodals_positions);
635 }
636
637 // this is alternative method for setting boundary conditions, to bypass bu
638 // in cubit file reader. not elegant, but good enough
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");
645 }
646 setOfConvection[it->getMeshsetId()].cOnvection = data[0];
647 setOfConvection[it->getMeshsetId()].tEmperature = data[1];
648 CHKERR mField.get_moab().get_entities_by_type(
649 it->meshset, MBTRI, setOfConvection[it->getMeshsetId()].tRis, true);
651 setOfConvection[it->getMeshsetId()].tRis, MBTRI,
652 "THERMAL_CONVECTION_FE");
653 }
654 }
655
657}
658
660 const std::string field_name, const std::string mesh_nodals_positions) {
662
663 CHKERR mField.add_finite_element("THERMAL_RADIATION_FE", MF_ZERO);
665 field_name);
667 field_name);
669 field_name);
670 if (mField.check_field(mesh_nodals_positions)) {
672 mesh_nodals_positions);
673 }
674
675 // this is alternative method for setting boundary conditions, to bypass bu
676 // in cubit file reader. not elegant, but good enough
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");
683 }
684 setOfRadiation[it->getMeshsetId()].sIgma = data[0];
685 setOfRadiation[it->getMeshsetId()].eMissivity = data[1];
686 setOfRadiation[it->getMeshsetId()].aMbienttEmp = data[2];
687 CHKERR mField.get_moab().get_entities_by_type(
688 it->meshset, MBTRI, setOfRadiation[it->getMeshsetId()].tRis, true);
690 setOfRadiation[it->getMeshsetId()].tRis, MBTRI,
691 "THERMAL_RADIATION_FE");
692 }
693 }
694
696}
697
701 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
702 feRhs.getOpPtrVector().push_back(
704 for (; sit != setOfBlocks.end(); sit++) {
705 // add finite element
706 feRhs.getOpPtrVector().push_back(
707 new OpThermalRhs(field_name, F, sit->second, commonData));
708 }
710}
711
715 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
716 for (; sit != setOfBlocks.end(); sit++) {
717 // add finite elemen
718 feLhs.getOpPtrVector().push_back(
719 new OpThermalLhs(field_name, A, sit->second, commonData));
720 }
722}
723
725 string field_name, Vec &F, const std::string mesh_nodals_positions) {
727 bool hoGeometry = false;
728 if (mField.check_field(mesh_nodals_positions)) {
729 hoGeometry = true;
730 }
731 std::map<int, FluxData>::iterator sit = setOfFluxes.begin();
732 for (; sit != setOfFluxes.end(); sit++) {
733 // add finite element
734 feFlux.getOpPtrVector().push_back(
735 new OpHeatFlux(field_name, F, sit->second, hoGeometry));
736 }
738}
739
741 string field_name, Vec &F, const std::string mesh_nodals_positions) {
743 bool hoGeometry = false;
744 if (mField.check_field(mesh_nodals_positions)) {
745 hoGeometry = true;
746 }
747 std::map<int, ConvectionData>::iterator sit = setOfConvection.begin();
748 for (; sit != setOfConvection.end(); sit++) {
749 // add finite element
753 field_name, F, sit->second, commonData, hoGeometry));
754 }
756}
757
759 string field_name, Mat A, const std::string mesh_nodals_positions) {
761 bool hoGeometry = false;
762 if (mField.check_field(mesh_nodals_positions)) {
763 hoGeometry = true;
764 }
765 std::map<int, ConvectionData>::iterator sit = setOfConvection.begin();
766 for (; sit != setOfConvection.end(); sit++) {
767 // add finite element
769 new OpConvectionLhs(field_name, A, sit->second, hoGeometry));
770 }
772}
773
775 string field_name, string rate_name,
776 const std::string mesh_nodals_positions) {
778
779 bool hoGeometry = false;
780 if (mField.check_field(mesh_nodals_positions)) {
781 hoGeometry = true;
782 }
783
784 {
785 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
786 for (; sit != setOfBlocks.end(); sit++) {
787 // add finite element
788 // those methods are to calculate matrices on Lhs
789 // feLhs.getOpPtrVector().push_back(new
790 // OpGetTetTemperatureAtGaussPts(field_name,commonData));
791 feLhs.getOpPtrVector().push_back(
792 new OpThermalLhs(field_name, sit->second, commonData));
793 feLhs.getOpPtrVector().push_back(
794 new OpHeatCapacityLhs(field_name, sit->second, commonData));
795 // those methods are to calculate vectors on Rhs
796 feRhs.getOpPtrVector().push_back(
798 feRhs.getOpPtrVector().push_back(
799 new OpGetTetRateAtGaussPts(rate_name, commonData));
800 feRhs.getOpPtrVector().push_back(
802 feRhs.getOpPtrVector().push_back(
803 new OpThermalRhs(field_name, sit->second, commonData));
804 feRhs.getOpPtrVector().push_back(
805 new OpHeatCapacityRhs(field_name, sit->second, commonData));
806 }
807 }
808
809 // Flux
810 {
811 std::map<int, FluxData>::iterator sit = setOfFluxes.begin();
812 for (; sit != setOfFluxes.end(); sit++) {
813 feFlux.getOpPtrVector().push_back(
814 new OpHeatFlux(field_name, sit->second, hoGeometry));
815 }
816 }
817
818 // Convection
819 {
820 std::map<int, ConvectionData>::iterator sit = setOfConvection.begin();
821 for (; sit != setOfConvection.end(); sit++) {
825 new OpConvectionRhs(field_name, sit->second, commonData, hoGeometry));
826 }
827 }
828 {
829 std::map<int, ConvectionData>::iterator sit = setOfConvection.begin();
830 for (; sit != setOfConvection.end(); sit++) {
832 new OpConvectionLhs(field_name, sit->second, hoGeometry));
833 }
834 }
835
836 // Radiation
837 {
838 std::map<int, RadiationData>::iterator sit = setOfRadiation.begin();
839 for (; sit != setOfRadiation.end(); sit++) {
840 feRadiationRhs.getOpPtrVector().push_back(
842 feRadiationRhs.getOpPtrVector().push_back(
843 new OpRadiationRhs(field_name, sit->second, commonData, hoGeometry));
844 }
845 }
846 {
847 std::map<int, RadiationData>::iterator sit = setOfRadiation.begin();
848 for (; sit != setOfRadiation.end(); sit++) {
849 feRadiationLhs.getOpPtrVector().push_back(
851 feRadiationLhs.getOpPtrVector().push_back(
852 new OpRadiationLhs(field_name, sit->second, commonData, hoGeometry));
853 }
854 }
855
857}
858
860 TsCtx &ts_ctx, string field_name, string rate_name,
861 const std::string mesh_nodals_positions) {
863
864 CHKERR setTimeSteppingProblem(field_name, rate_name, mesh_nodals_positions);
865
866 // rhs
867 TsCtx::FEMethodsSequence &loops_to_do_Rhs =
869 loops_to_do_Rhs.push_back(TsCtx::PairNameFEMethodPtr("THERMAL_FE", &feRhs));
870 loops_to_do_Rhs.push_back(
871 TsCtx::PairNameFEMethodPtr("THERMAL_FLUX_FE", &feFlux));
872 if (mField.check_finite_element("THERMAL_CONVECTION_FE"))
873 loops_to_do_Rhs.push_back(
874 TsCtx::PairNameFEMethodPtr("THERMAL_CONVECTION_FE", &feConvectionRhs));
875 if (mField.check_finite_element("THERMAL_RADIATION_FE"))
876 loops_to_do_Rhs.push_back(
877 TsCtx::PairNameFEMethodPtr("THERMAL_RADIATION_FE", &feRadiationRhs));
878
879 // lhs
880 TsCtx::FEMethodsSequence &loops_to_do_Mat =
882 loops_to_do_Mat.push_back(TsCtx::PairNameFEMethodPtr("THERMAL_FE", &feLhs));
883 if (mField.check_finite_element("THERMAL_CONVECTION_FE"))
884 loops_to_do_Mat.push_back(
885 TsCtx::PairNameFEMethodPtr("THERMAL_CONVECTION_FE", &feConvectionLhs));
886 if (mField.check_finite_element("THERMAL_RADIATION_FE"))
887 loops_to_do_Mat.push_back(
888 TsCtx::PairNameFEMethodPtr("THERMAL_RADIATION_FE", &feRadiationLhs));
889
891}
static Index< 'M', 3 > M
static Index< 'p', 3 > p
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
Operators and data structures for thermal analysis.
@ QUIET
Definition: definitions.h:208
@ ROW
Definition: definitions.h:123
@ MF_ZERO
Definition: definitions.h:98
@ MF_EXIST
Definition: definitions.h:100
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ HEATFLUXSET
Definition: definitions.h:156
@ SIDESET
Definition: definitions.h:147
@ MAT_THERMALSET
block name is "MAT_THERMAL"
Definition: definitions.h:161
@ BLOCKSET
Definition: definitions.h:148
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
@ 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_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.
Definition: LogManager.hpp:308
#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
Definition: level_set.cpp:1884
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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.
Definition: Templates.hpp:135
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr AssemblyType A
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:15
MoFEM::FEMethodsSequence FEMethodsSequence
Definition: TsCtx.hpp:24
FEMethodsSequence & getLoopsIFunction()
Get the loops to do IFunction object.
Definition: TsCtx.hpp:69
FEMethodsSequence & getLoopsIJacobian()
Get the loops to do IJacobian object.
Definition: TsCtx.hpp:89
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
VectorDouble temperatureRateAtGaussPts
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(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
std::map< int, FluxData > setOfFluxes
maps side set id with appropriate FluxData
MyTriFE feRadiationLhs
MyTriFE feConvectionLhs
MyVolumeFE feRhs
cauclate right hand side for tetrahedral elements
MyTriFE feConvectionRhs
MoFEM::Interface & mField
CommonData commonData
MoFEMErrorCode setThermalConvectionFiniteElementLhsOperators(string field_name, Mat A, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
MyVolumeFE feLhs
MyTriFE feRadiationRhs
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