v0.15.0
Loading...
Searching...
No Matches
ThermalRadiation.hpp
Go to the documentation of this file.
1/**
2 * @file ThermalRadiation.hpp
3 * @author Andrei Shvarts (andrei.shvarts@glasgow.ac.uk)
4 * @brief Implementation of thermal radiation bc
5 * @version 0.14.0
6 * @date 2024-08-23
7 */
8
9namespace ThermoElasticOps {
10
11template <CubitBC BC> struct RadiationBcType {};
12
13template <CubitBC BC> struct GetRadiationParameters;
14
15template <> struct GetRadiationParameters<BLOCKSET> {
17
18 static MoFEMErrorCode
19 getParameters(double &stefan_boltzmann_constant, double &emissivity_parameter,
20 double &ambient_temperature, boost::shared_ptr<Range> &ents,
21 MoFEM::Interface &m_field, int ms_id, Sev sev) {
22
24
25 auto cubit_meshset_ptr =
26 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(ms_id,
27 BLOCKSET);
28 std::vector<double> block_data;
29 CHKERR cubit_meshset_ptr->getAttributes(block_data);
30 if (block_data.size() < 3) {
31 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
32 "Expected that radiation block has three attributes "
33 "(Stefan–Boltzmann "
34 "constant, emissivity parameter, and ambient temperature)");
35 }
36 stefan_boltzmann_constant = block_data[0];
37 emissivity_parameter = block_data[1];
38 ambient_temperature = block_data[2];
39
40 MOFEM_LOG_CHANNEL("WORLD");
41 MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "RadiationBc")
42 << "Stefan–Boltzmann constant " << stefan_boltzmann_constant;
43 MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "RadiationBc")
44 << "Emissivity parameter " << emissivity_parameter;
45 MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "RadiationBc")
46 << "Ambient temperature " << ambient_temperature;
47
48 ents = boost::make_shared<Range>();
49 CHKERR m_field.get_moab().get_entities_by_handle(cubit_meshset_ptr->meshset,
50 *(ents), true);
51
52 MOFEM_LOG_CHANNEL("SYNC");
53 MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "RadiationBc") << *ents;
54 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::noisy);
55
57 }
58};
59
60} // namespace ThermoElasticOps
61
62template <AssemblyType A, typename EleOp>
63struct OpFluxRhsImpl<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1, A,
64 GAUSS, EleOp> : OpBaseImpl<A, EleOp> {
65
67
68 OpFluxRhsImpl(MoFEM::Interface &m_field, int ms_id, std::string field_name,
69 boost::shared_ptr<VectorDouble> t_ptr);
70
71protected:
75 boost::shared_ptr<VectorDouble> tPtr;
76 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
77};
78
79template <AssemblyType A, typename EleOp>
80struct OpFluxLhsImpl<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1, A,
81 GAUSS, EleOp> : OpBaseImpl<A, EleOp> {
82
84
85 OpFluxLhsImpl(MoFEM::Interface &m_field, int ms_id,
86 std::string row_field_name, std::string col_field_name,
87 boost::shared_ptr<VectorDouble> t_ptr);
88
89protected:
93 boost::shared_ptr<VectorDouble>
94 tPtr; // pointer for boundary temperature field
95 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
96 EntitiesFieldData::EntData &col_data);
97};
98
99template <AssemblyType A, typename EleOp>
100OpFluxRhsImpl<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1, A, GAUSS,
101 EleOp>::OpFluxRhsImpl(MoFEM::Interface &m_field, int ms_id,
102 std::string field_name,
103 boost::shared_ptr<VectorDouble> t_ptr)
104 : OpBase(field_name, field_name, OpBase::OPROW), tPtr(t_ptr) {
107 this->stefanBoltzmanConstant, this->emmisivityParameter,
108 this->ambientTemperature, this->entsPtr, m_field, ms_id, Sev::inform),
109 "Cannot read radiation data from blockset");
110}
111
112template <AssemblyType A, typename EleOp>
113MoFEMErrorCode
115 GAUSS, EleOp>::iNtegrate(EntitiesFieldData::EntData
116 &row_data) {
118
119 // get integration weights
120 auto t_w = OpBase::getFTensor0IntegrationWeight();
121 // get base function values on rows
122 auto t_row_base = row_data.getFTensor0N();
123
124 // get temperature
125 auto t_temp = getFTensor0FromVec(*tPtr);
126 // loop over integration points
127 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
128 const auto alpha = -t_w * (std::pow(static_cast<double>(t_temp), 4) -
129 std::pow(ambientTemperature, 4));
130
131 int rr = 0;
132 for (; rr != OpBase::nbRows; ++rr) {
133 OpBase::locF[rr] += (alpha * t_row_base);
134 ++t_row_base;
135 }
136
137 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
138 ++t_row_base;
139
140 ++t_w;
141 ++t_temp;
142 }
143
144 // get normal vector
145 auto t_normal = OpBase::getFTensor1Normal();
146 OpBase::locF *= t_normal.l2() * stefanBoltzmanConstant * emmisivityParameter;
147
148 EntityType fe_type = OpBase::getNumeredEntFiniteElementPtr()->getEntType();
149 if (fe_type == MBTRI) {
150 OpBase::locF /= 2;
151 }
152
154}
155
156template <AssemblyType A, typename EleOp>
157OpFluxLhsImpl<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1, A, GAUSS,
158 EleOp>::OpFluxLhsImpl(MoFEM::Interface &m_field, int ms_id,
159 std::string row_field_name,
160 std::string col_field_name,
161 boost::shared_ptr<VectorDouble> t_ptr)
162 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL), tPtr(t_ptr) {
163
164 this->sYmm = true;
165
168 this->stefanBoltzmanConstant, this->emmisivityParameter,
169 this->ambientTemperature, this->entsPtr, m_field, ms_id,
170 Sev::verbose),
171 "Can not read radiation bc from blockset");
172}
173
174template <AssemblyType A, typename EleOp>
175MoFEMErrorCode
177 GAUSS,
178 EleOp>::iNtegrate(EntitiesFieldData::EntData &row_data,
179 EntitiesFieldData::EntData &col_data) {
180
182 // get integration weights
183 auto t_w = OpBase::getFTensor0IntegrationWeight();
184 // get base function values on rows
185 auto t_row_base = row_data.getFTensor0N();
186 // get temperature
187 auto t_temp = getFTensor0FromVec(*tPtr);
188 // loop over integration points
189 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
190 const auto alpha = t_w * 4 * std::pow(static_cast<double>(t_temp), 3);
191 ++t_w; // move to another integration weight
192 ++t_temp;
193 // loop over rows base functions
194 int rr = 0;
195 for (; rr != OpBase::nbRows; rr++) {
196 const auto beta = alpha * t_row_base;
197 // get column base functions values at gauss point gg
198 auto t_col_base = col_data.getFTensor0N(gg, 0);
199 // loop over columns
200 for (int cc = 0; cc != OpBase::nbCols; cc++) {
201 OpBase::locMat(rr, cc) += beta * t_col_base;
202 ++t_col_base;
203 }
204 ++t_row_base;
205 }
206 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
207 ++t_row_base;
208 }
209 // get normal vector
210 auto t_normal = OpBase::getFTensor1Normal();
211 // using L2 norm of normal vector for convenient area calculation for quads,
212 // tris etc.
214 -t_normal.l2() * stefanBoltzmanConstant * emmisivityParameter;
215
216 EntityType fe_type = OpBase::getNumeredEntFiniteElementPtr()->getEntType();
217 if (fe_type == MBTRI) {
218 OpBase::locMat /= 2;
219 }
220
222}
223
224template <CubitBC BCTYPE, int BASE_DIM, int FIELD_DIM, AssemblyType A,
225 IntegrationType I, typename OpBase>
226struct AddFluxToRhsPipelineImpl<
227
228 OpFluxRhsImpl<ThermoElasticOps::RadiationBcType<BCTYPE>, BASE_DIM,
229 FIELD_DIM, A, I, OpBase>,
230 A, I, OpBase
231
232 > {
233
235
236 static MoFEMErrorCode add(
237
238 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
239 MoFEM::Interface &m_field, std::string field_name,
240 boost::shared_ptr<VectorDouble> T_ptr, std::string block_name, Sev sev
241
242 ) {
244
245 using OP = OpFluxRhsImpl<ThermoElasticOps::RadiationBcType<BLOCKSET>,
247
248 auto add_op = [&](auto &&meshset_vec_ptr) {
249 for (auto m : meshset_vec_ptr) {
250 MOFEM_TAG_AND_LOG("WORLD", sev, "OpRadiationRhs") << "Add " << *m;
251 pipeline.push_back(
252 new OP(m_field, m->getMeshsetId(), field_name, T_ptr));
253 }
254 MOFEM_LOG_CHANNEL("WORLD");
255 };
256
257 switch (BCTYPE) {
258 case BLOCKSET:
259 add_op(
260
261 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
262 std::regex(
263
264 (boost::format("%s(.*)") % block_name).str()
265
266 ))
267
268 );
269
270 break;
271 default:
272 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
273 "Handling of bc type not implemented");
274 break;
275 }
277 }
278};
279
280template <CubitBC BCTYPE, int BASE_DIM, int FIELD_DIM, AssemblyType A,
281 IntegrationType I, typename OpBase>
282struct AddFluxToLhsPipelineImpl<
283
284 OpFluxLhsImpl<ThermoElasticOps::RadiationBcType<BCTYPE>, BASE_DIM,
285 FIELD_DIM, A, I, OpBase>,
286 A, I, OpBase
287
288 > {
289
291
292 static MoFEMErrorCode add(
293
294 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
295 MoFEM::Interface &m_field, std::string row_field_name,
296 std::string col_field_name, boost::shared_ptr<VectorDouble> T_ptr,
297 std::string block_name, Sev sev
298
299 ) {
301
302 using OP = OpFluxLhsImpl<ThermoElasticOps::RadiationBcType<BLOCKSET>,
304
305 auto add_op = [&](auto &&meshset_vec_ptr) {
306 for (auto m : meshset_vec_ptr) {
307 MOFEM_TAG_AND_LOG("WORLD", sev, "OpRadiationLhs") << "Add " << *m;
308 pipeline.push_back(new OP(m_field, m->getMeshsetId(), row_field_name,
309 col_field_name, T_ptr));
310 }
311 MOFEM_LOG_CHANNEL("WORLD");
312 };
313
314 switch (BCTYPE) {
315 case BLOCKSET:
316 add_op(
317
318 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
319 std::regex(
320
321 (boost::format("%s(.*)") % block_name).str()
322
323 ))
324
325 );
326
327 break;
328 default:
329 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
330 "Handling of bc type not implemented");
331 break;
332 }
334 }
335};
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
constexpr int FIELD_DIM
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
CubitBC
Types of sets and boundary conditions.
@ BLOCKSET
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int BASE_DIM
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
constexpr IntegrationType I
constexpr AssemblyType A
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static MoFEMErrorCode add(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, MoFEM::Interface &m_field, std::string row_field_name, std::string col_field_name, boost::shared_ptr< VectorDouble > T_ptr, std::string block_name, Sev sev)
static MoFEMErrorCode add(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, MoFEM::Interface &m_field, std::string field_name, boost::shared_ptr< VectorDouble > T_ptr, std::string block_name, Sev sev)
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Deprecated interface functions.
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
static MoFEMErrorCode getParameters(double &stefan_boltzmann_constant, double &emissivity_parameter, double &ambient_temperature, boost::shared_ptr< Range > &ents, MoFEM::Interface &m_field, int ms_id, Sev sev)