v0.15.0
Loading...
Searching...
No Matches
ThermalConvection.hpp
Go to the documentation of this file.
1/**
2 * @file ThermalConvection.hpp
3 * @author Andrei Shvarts (andrei.shvarts@glasgow.ac.uk)
4 * @brief Implementation of thermal convection bc
5 * @version 0.14.0
6 * @date 2024-08-23
7 */
8
9namespace ThermoElasticOps {
10
11template <CubitBC BC> struct ConvectionBcType {};
12
13template <CubitBC BC> struct GetConvectionParameters;
14
15template <> struct GetConvectionParameters<BLOCKSET> {
17
18 static MoFEMErrorCode getParameters(double &heat_transfer_coefficient,
19 double &ambient_temperature,
20 boost::shared_ptr<Range> &ents,
21 MoFEM::Interface &m_field, int ms_id,
22 Sev sev) {
23
25
26 auto cubit_meshset_ptr =
27 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(ms_id,
28 BLOCKSET);
29 std::vector<double> block_data;
30 CHKERR cubit_meshset_ptr->getAttributes(block_data);
31 if (block_data.size() < 2) {
32 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
33 "Expected that convection block has two attributes (heat "
34 "transfer coefficient and ambient temperature)");
35 }
36 heat_transfer_coefficient = block_data[0];
37 ambient_temperature = block_data[1];
38
39 MOFEM_LOG_CHANNEL("WORLD");
40 MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "ConvectionBc")
41 << "Heat transfer coefficient " << heat_transfer_coefficient;
42 MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "ConvectionBc")
43 << "Ambient temperature " << ambient_temperature;
44
45 ents = boost::make_shared<Range>();
46 CHKERR
47 m_field.get_moab().get_entities_by_handle(cubit_meshset_ptr->meshset,
48 *(ents), true);
49
50 MOFEM_LOG_CHANNEL("SYNC");
51 MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "ConvectionBc") << *ents;
52 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::noisy);
53
55 }
56};
57
58} // namespace ThermoElasticOps
59
60template <AssemblyType A, typename EleOp>
61struct OpFluxRhsImpl<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1, A,
62 GAUSS, EleOp> : OpBaseImpl<A, EleOp> {
63
65
66 OpFluxRhsImpl(MoFEM::Interface &m_field, int ms_id, std::string field_name,
67 boost::shared_ptr<VectorDouble> t_ptr);
68
69protected:
72 boost::shared_ptr<VectorDouble> tPtr;
73 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
74};
75
76template <AssemblyType A, typename EleOp>
77struct OpFluxLhsImpl<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1, A,
78 GAUSS, EleOp> : OpBaseImpl<A, EleOp> {
79
81
82 OpFluxLhsImpl(MoFEM::Interface &m_field, int ms_id,
83 std::string row_field_name, std::string col_field_name);
84
85protected:
88 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
89 EntitiesFieldData::EntData &col_data);
90};
91
92template <AssemblyType A, typename EleOp>
93OpFluxRhsImpl<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1, A, GAUSS,
94 EleOp>::OpFluxRhsImpl(MoFEM::Interface &m_field, int ms_id,
95 std::string field_name,
96 boost::shared_ptr<VectorDouble> t_ptr)
97 : OpBase(field_name, field_name, OpBase::OPROW), tPtr(t_ptr) {
100 this->heatTransferCoefficient, this->ambientTemperature,
101 this->entsPtr, m_field, ms_id, Sev::inform),
102 "Cannot read convection data from blockset");
103}
104
105template <AssemblyType A, typename EleOp>
106MoFEMErrorCode
108 GAUSS, EleOp>::iNtegrate(EntitiesFieldData::EntData
109 &row_data) {
111
112 // get integration weights
113 auto t_w = OpBase::getFTensor0IntegrationWeight();
114 // get base function values on rows
115 auto t_row_base = row_data.getFTensor0N();
116
117 // get temperature
118 auto t_temp = getFTensor0FromVec(*tPtr);
119 // loop over integration points
120 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
121 const auto alpha =
122 -t_w * heatTransferCoefficient * (t_temp - ambientTemperature);
123 int rr = 0;
124 for (; rr != OpBase::nbRows; ++rr) {
125 OpBase::locF[rr] += alpha * t_row_base;
126 ++t_row_base;
127 }
128
129 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
130 ++t_row_base;
131
132 ++t_w;
133 ++t_temp;
134 }
135
136 // get normal vector
137 auto t_normal = OpBase::getFTensor1Normal();
138 OpBase::locF *= t_normal.l2();
139
140 EntityType fe_type = OpBase::getNumeredEntFiniteElementPtr()->getEntType();
141 if (fe_type == MBTRI) {
142 OpBase::locF /= 2;
143 }
144
146}
147
148template <AssemblyType A, typename EleOp>
149OpFluxLhsImpl<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1, A, GAUSS,
150 EleOp>::OpFluxLhsImpl(MoFEM::Interface &m_field, int ms_id,
151 std::string row_field_name,
152 std::string col_field_name)
153 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL) {
154
155 this->sYmm = true;
156
159 this->heatTransferCoefficient, this->ambientTemperature,
160 this->entsPtr, m_field, ms_id, Sev::verbose),
161 "Cannot read convection data from blockset");
162}
163
164template <AssemblyType A, typename EleOp>
165MoFEMErrorCode
167 GAUSS,
168 EleOp>::iNtegrate(EntitiesFieldData::EntData &row_data,
169 EntitiesFieldData::EntData &col_data) {
170
172 // get integration weights
173 auto t_w = OpBase::getFTensor0IntegrationWeight();
174 // get base function values on rows
175 auto t_row_base = row_data.getFTensor0N();
176 // loop over integration points
177 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
178 // loop over rows base functions
179 int rr = 0;
180 for (; rr != OpBase::nbRows; rr++) {
181 const auto alpha = t_w * t_row_base;
182 // get column base functions values at gauss point gg
183 auto t_col_base = col_data.getFTensor0N(gg, 0);
184 // loop over columns
185 for (int cc = 0; cc != OpBase::nbCols; cc++) {
186 OpBase::locMat(rr, cc) += alpha * t_col_base;
187 ++t_col_base;
188 }
189 ++t_row_base;
190 }
191 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
192 ++t_row_base;
193 ++t_w; // move to another integration weight
194 }
195 // get normal vector
196 auto t_normal = OpBase::getFTensor1Normal();
197 // using L2 norm of normal vector for convenient area calculation for quads,
198 // tris etc.
199 OpBase::locMat *= -t_normal.l2() * heatTransferCoefficient;
200
201 EntityType fe_type = OpBase::getNumeredEntFiniteElementPtr()->getEntType();
202 if (fe_type == MBTRI) {
203 OpBase::locMat /= 2;
204 }
205
207}
208
209template <CubitBC BCTYPE, int BASE_DIM, int FIELD_DIM, AssemblyType A,
210 IntegrationType I, typename OpBase>
211struct AddFluxToRhsPipelineImpl<
212
213 OpFluxRhsImpl<ThermoElasticOps::ConvectionBcType<BCTYPE>, BASE_DIM,
214 FIELD_DIM, A, I, OpBase>,
215 A, I, OpBase
216
217 > {
218
220
221 static MoFEMErrorCode add(
222
223 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
224 MoFEM::Interface &m_field, std::string field_name,
225 boost::shared_ptr<VectorDouble> u_ptr, std::string block_name, Sev sev
226
227 ) {
229
230 using OP = OpFluxRhsImpl<ThermoElasticOps::ConvectionBcType<BLOCKSET>,
232
233 auto add_op = [&](auto &&meshset_vec_ptr) {
234 for (auto m : meshset_vec_ptr) {
235 MOFEM_TAG_AND_LOG("WORLD", sev, "OpConvectionRhs") << "Add " << *m;
236 pipeline.push_back(
237 new OP(m_field, m->getMeshsetId(), field_name, u_ptr));
238 }
239 MOFEM_LOG_CHANNEL("WORLD");
240 };
241
242 switch (BCTYPE) {
243 case BLOCKSET:
244 add_op(
245
246 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
247 std::regex(
248
249 (boost::format("%s(.*)") % block_name).str()
250
251 ))
252
253 );
254
255 break;
256 default:
257 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
258 "Handling of bc type not implemented");
259 break;
260 }
262 }
263};
264
265template <CubitBC BCTYPE, int BASE_DIM, int FIELD_DIM, AssemblyType A,
266 IntegrationType I, typename OpBase>
267struct AddFluxToLhsPipelineImpl<
268
269 OpFluxLhsImpl<ThermoElasticOps::ConvectionBcType<BCTYPE>, BASE_DIM,
270 FIELD_DIM, A, I, OpBase>,
271 A, I, OpBase
272
273 > {
274
276
277 static MoFEMErrorCode add(
278
279 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
280 MoFEM::Interface &m_field, std::string row_field_name,
281 std::string col_field_name, std::string block_name, Sev sev
282
283 ) {
285
286 using OP = OpFluxLhsImpl<ThermoElasticOps::ConvectionBcType<BLOCKSET>,
288
289 auto add_op = [&](auto &&meshset_vec_ptr) {
290 for (auto m : meshset_vec_ptr) {
291 MOFEM_TAG_AND_LOG("WORLD", sev, "OpConvectionLhs") << "Add " << *m;
292 pipeline.push_back(
293 new OP(m_field, m->getMeshsetId(), row_field_name, col_field_name));
294 }
295 MOFEM_LOG_CHANNEL("WORLD");
296 };
297
298 switch (BCTYPE) {
299 case BLOCKSET:
300 add_op(
301
302 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
303 std::regex(
304
305 (boost::format("%s(.*)") % block_name).str()
306
307 ))
308
309 );
310
311 break;
312 default:
313 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
314 "Handling of bc type not implemented");
315 break;
316 }
318 }
319};
#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, 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 > u_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 &heat_transfer_coefficient, double &ambient_temperature, boost::shared_ptr< Range > &ents, MoFEM::Interface &m_field, int ms_id, Sev sev)