b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2element_heat.H
1//------------------------------------------------------------------------
2// b2element_heat.H --
3//
4//
5// written by Mathias Doreille
6// Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
7// Thomas Blome <thomas.blome@dlr.de>
8//
9// (c) 2007-2012,2016 SMR Engineering & Development SA
10// 2502 Bienne, Switzerland
11//
12// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
13// Linder Höhe, 51147 Köln
14//
15// All Rights Reserved. Proprietary source code. The contents of
16// this file may not be disclosed to third parties, copied or
17// duplicated in any form, in whole or in part, without the prior
18// written permission of SMR.
19//------------------------------------------------------------------------
20
21#ifndef B2_ELEMENT_HEAT_H_
22#define B2_ELEMENT_HEAT_H_
23
24#include "b2element_continuum_util.H"
25#include "elements/properties/b2heat_material.H"
26#include "model/b2element.H"
27#include "model/b2solution.H"
28#include "model_imp/b2hnode.H"
29#include "utils/b2constants.H"
30#include "utils/b2exception.H"
31#include "utils/b2linear_algebra.H"
32#include "utils/b2tensor_calculus.H"
33
34// #define LUMPED_CAPACITY
35
36namespace b2000 {
37
38template <typename SHAPE, typename INTEGRATION, int NB_GAUSS>
39class GenericInterpolation {
40public:
41 static const int nb_dimensions = SHAPE::num_dimensions;
42 static const int nb_nodes = SHAPE::num_nodes;
43 static const int nb_gauss = NB_GAUSS;
44
45 GenericInterpolation() {
46 b2linalg::Vector<double, b2linalg::Vdense> v(nb_nodes);
47 b2linalg::Matrix<double, b2linalg::Mrectangle> m(nb_dimensions, nb_nodes);
48 INTEGRATION integration;
49 integration.set_num(NB_GAUSS);
50 if (NB_GAUSS != integration.get_num()) {
51 Exception() << "Cannot found an integration schemas that contain " << NB_GAUSS
52 << " integration points." << THROW;
53 }
54 IP_ID ip_id;
55 for (int l = 0; l != nb_gauss; ++l) {
56 integration.get_point(l, ip_id, array[l].weight);
57 for (int i = 0; i != nb_dimensions; ++i) { array[l].IntCoor[i] = ip_id[i]; }
58 SHAPE::eval_h(array[l].IntCoor, v);
59 std::copy(v.begin(), v.end(), array[l].N);
60 SHAPE::eval_h_derived(array[l].IntCoor, m);
61 for (size_t i = 0; i != m.size1(); ++i) {
62 for (size_t j = 0; j != m.size2(); ++j) { array[l].dN[j][i] = m(i, j); }
63 }
64 }
65 }
66
67 struct OnePoint {
68 // Integration point coordinate
69 double IntCoor[nb_dimensions];
70
71 // Integration point weight
72 double weight;
73
74 // The shape functions at the Gauss points.
75 double N[nb_nodes];
76
77 // The derivative of shape function at the Gauss points.
78 double dN[nb_nodes][nb_dimensions];
79 };
80
81 static OnePoint array[nb_gauss];
82};
83
84template <typename SHAPE, typename INTEGRATION, int NB_GAUSS>
85typename GenericInterpolation<SHAPE, INTEGRATION, NB_GAUSS>::OnePoint
86 GenericInterpolation<SHAPE, INTEGRATION, NB_GAUSS>::array[nb_gauss];
87
88template <typename SHAPE, typename INTEGRATION, int NB_GAUSS>
89const int GenericInterpolation<SHAPE, INTEGRATION, NB_GAUSS>::nb_dimensions;
90
91template <typename SHAPE, typename INTEGRATION, int NB_GAUSS>
92const int GenericInterpolation<SHAPE, INTEGRATION, NB_GAUSS>::nb_nodes;
93
94template <typename SHAPE, typename INTEGRATION, int NB_GAUSS>
95const int GenericInterpolation<SHAPE, INTEGRATION, NB_GAUSS>::nb_gauss;
96
97template <int NB_DIMMENSION>
98struct HeatMaterialType;
99
100template <>
101struct HeatMaterialType<2> {
102 typedef HeatMaterial2D heat_material_type;
103 static double get_thickness(
104 const heat_material_type* properties, const Element* element,
105 b2linalg::Vector<double, b2linalg::Vdense_constref> nodes_interpolation) {
106 return properties->get_thickness(element, nodes_interpolation);
107 }
108};
109
110template <>
111struct HeatMaterialType<3> {
112 typedef HeatMaterial3D heat_material_type;
113 static double get_thickness(
114 const heat_material_type* properties, const Element* element,
115 b2linalg::Vector<double, b2linalg::Vdense_constref> nodes_interpolation) {
116 return 0;
117 }
118};
119
120template <typename SHAPE, typename INTEGRATION, int NB_GAUSS, bool AXISYMMETRIC = false>
121class ElementHeatConduction : public TypedElement<double> {
122public:
123 using N_t = GenericInterpolation<SHAPE, INTEGRATION, NB_GAUSS>;
124
125 static const int nb_dimensions = N_t::nb_dimensions;
126
127 using node_type = node::HNode<
128 coordof::Coor<coordof::Trans<coordof::X, coordof::Y, coordof::Z>>,
129 coordof::Dof<coordof::Temperature>>;
130
131 using heat_material_type = typename HeatMaterialType<nb_dimensions>::heat_material_type;
132
133 using type_t = ObjectTypeComplete<ElementHeatConduction, typename TypedElement<double>::type_t>;
134
135 void set_nodes(std::pair<int, Node* const*> nodes_) override {
136 if (nodes_.first != N_t::nb_nodes) {
137 Exception() << "This element has " << N_t::nb_nodes << " nodes." << THROW;
138 }
139 for (int i = 0; i != N_t::nb_nodes; ++i) { nodes[i] = nodes_.second[i]; }
140 }
141
142 std::pair<int, Node* const*> get_nodes() const override {
143 return std::pair<int, Node* const*>(N_t::nb_nodes, nodes);
144 }
145
146 void set_property(ElementProperty* property_) override {
147 property = dynamic_cast<heat_material_type*>(property_);
148 if (property == 0) {
149 TypeError() << "Bad property type " << typeid(*property_)
150 << " for heat element (forgot MID or PID element attribute?)" << THROW;
151 }
152 }
153
154 const ElementProperty* get_property() const override { return property; }
155
156 void get_dof_numbering(b2linalg::Index& dof_numbering) override {
157 dof_numbering.resize(N_t::nb_nodes);
158 size_t* ptr = &dof_numbering[0];
159 for (int k = 0; k != N_t::nb_nodes; ++k) {
160 ptr = node_type::get_global_dof_numbering_s(ptr, nodes[k]);
161 }
162 }
163
164 const std::vector<VariableInfo> get_value_info() const override {
165 return std::vector<VariableInfo>(2, Element::nonlinear);
166 }
167
168 void get_value(
169 Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
170 const double time, const double delta_time,
171 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof,
172 GradientContainer* gradient_container, SolverHints* solver_hints,
173 b2linalg::Index& dof_numbering, b2linalg::Vector<double, b2linalg::Vdense>& value,
174 const std::vector<bool>& d_value_d_dof_flags,
175 std::vector<b2linalg::Matrix<double, b2linalg::Mpacked>>& d_value_d_dof,
176 b2linalg::Vector<double, b2linalg::Vdense>& d_value_d_time) override;
177
178 int edge_field_order(const int edge_id, const std::string& field_name) override {
179 return SHAPE::edge_shape_type_0::order;
180 }
181
182 bool edge_field_linear_on_dof(const int edge_id, const std::string& field_name) override {
183 return true;
184 }
185
186 int edge_field_polynomial_sub_edge(
187 const int edge_id, const std::string& field_name,
188 b2linalg::Vector<double, b2linalg::Vdense>& sub_nodes) override {
189 sub_nodes.clear();
190 sub_nodes.reserve(2);
191 sub_nodes.push_back(-1);
192 sub_nodes.push_back(1);
193 return 1;
194 }
195
196 void edge_geom(
197 const int edge_id, const double internal_coor, b2linalg::Vector<double>& geom,
198 b2linalg::Vector<double>& d_geom_d_icoor) override {
199 typedef typename SHAPE::edge_shape_type_0 edge_shape;
200 const int* edge_node = SHAPE::edge_node[edge_id - 1];
201 double node_coor[edge_shape::num_nodes][3];
202 for (int k = 0; k != edge_shape::num_nodes; ++k) {
203 node_type::get_coor_s(node_coor[k], nodes[edge_node[k]]);
204 }
205 if (!geom.is_null()) {
206 b2linalg::Vector<double> shape(edge_shape::num_nodes);
207 edge_shape::eval_h(&internal_coor, shape);
208 geom.resize(3);
209 for (size_t i = 0; i != 3; ++i) {
210 double tmp = 0;
211 for (int k = 0; k != edge_shape::num_nodes; ++k) {
212 tmp += shape[k] * node_coor[k][i];
213 }
214 geom[i] = tmp;
215 }
216 }
217 if (!d_geom_d_icoor.is_null()) {
218 b2linalg::Matrix<double> d_shape(1, edge_shape::num_nodes);
219 edge_shape::eval_h_derived(&internal_coor, d_shape);
220 d_geom_d_icoor.resize(3);
221 for (size_t i = 0; i != 3; ++i) {
222 double tmp = 0;
223 for (int k = 0; k != edge_shape::num_nodes; ++k) {
224 tmp += d_shape(0, k) * node_coor[k][i];
225 }
226 d_geom_d_icoor[i] = tmp;
227 }
228 }
229 }
230
231 void edge_field_value(
232 const int edge_id, const std::string& field_name, const double internal_coor,
233 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof, const double time,
234 b2linalg::Vector<double, b2linalg::Vdense>& value,
235 b2linalg::Vector<double, b2linalg::Vdense>& d_value_d_icoor,
236 b2linalg::Index& dof_numbering,
237 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_value_d_dof,
238 b2linalg::Index& d_value_d_dof_dep) override {
239 typedef typename SHAPE::edge_shape_type_0 edge_shape;
240 const int* edge_node = SHAPE::edge_node[edge_id - 1];
241 size_t dn[edge_shape::num_nodes];
242 {
243 size_t* ptr = &dn[0];
244 for (int k = 0; k != edge_shape::num_nodes; ++k) {
245 ptr = node_type::get_global_dof_numbering_s(ptr, nodes[edge_node[k]]);
246 }
247 }
248 if (!value.is_null()) {
249 b2linalg::Vector<double> shape(edge_shape::num_nodes);
250 edge_shape::eval_h(&internal_coor, shape);
251 value.resize(1);
252 double& tmp = value[0];
253 tmp = 0;
254 for (int k = 0; k != edge_shape::num_nodes; ++k) { tmp += shape[k] * dof(0, dn[k]); }
255 }
256 if (!d_value_d_icoor.is_null()) {
257 b2linalg::Matrix<double> d_shape(1, edge_shape::num_nodes);
258 edge_shape::eval_h_derived(&internal_coor, d_shape);
259 d_value_d_icoor.resize(1);
260 double& tmp = d_value_d_icoor[0];
261 tmp = 0;
262 for (int k = 0; k != edge_shape::num_nodes; ++k) {
263 tmp += d_shape(0, k) * dof(0, dn[k]);
264 }
265 }
266 if (!d_value_d_dof.is_null()) {
267 dof_numbering.resize(edge_shape::num_nodes);
268 std::copy(dn, dn + edge_shape::num_nodes, dof_numbering.begin());
269 b2linalg::Vector<double> shape(edge_shape::num_nodes);
270 edge_shape::eval_h(&internal_coor, shape);
271 d_value_d_dof.resize(1, edge_shape::num_nodes);
272 d_value_d_dof.set_zero();
273 for (int k = 0; k != edge_shape::num_nodes; ++k) { d_value_d_dof(0, k) = shape[k]; }
274 }
275 }
276
277 int face_field_order(const int face_id, const std::string& field_name) override {
278 switch (SHAPE::face_type_index[face_id - 1]) {
279 case 0:
280 return SHAPE::face_shape_type_0::order;
281 case 1:
282 return SHAPE::face_shape_type_1::order;
283 }
284 return 0;
285 }
286
287 bool face_field_linear_on_dof(const int face_id, const std::string& field_name) override {
288 return true;
289 }
290
291 int face_field_polynomial_sub_face(
292 const int face_id, const std::string& field_name,
293 b2linalg::Matrix<double, b2linalg::Mrectangle>& sub_nodes,
294 std::vector<Triangle>& sub_faces) override {
295 if (SHAPE::is_t_face[face_id]) {
296 sub_nodes.resize(2, 3);
297 sub_nodes(0, 0) = 0;
298 sub_nodes(1, 0) = 0;
299 sub_nodes(0, 1) = 1;
300 sub_nodes(1, 1) = 0;
301 sub_nodes(0, 2) = 0;
302 sub_nodes(1, 2) = 1;
303 sub_faces.clear();
304 sub_faces.reserve(1);
305 sub_faces.push_back(Triangle(0, 1, 2));
306 } else {
307 // we split the rectangular face in two triangles
308 sub_nodes.resize(2, 4);
309 sub_nodes(0, 0) = -1;
310 sub_nodes(1, 0) = -1;
311 sub_nodes(0, 1) = 1;
312 sub_nodes(1, 1) = -1;
313 sub_nodes(0, 2) = 1;
314 sub_nodes(1, 2) = 1;
315 sub_nodes(0, 3) = -1;
316 sub_nodes(1, 3) = 1;
317 sub_faces.clear();
318 sub_faces.reserve(2);
319 sub_faces.push_back(Triangle(0, 1, 2));
320 sub_faces.push_back(Triangle(2, 3, 0));
321 }
322 return face_field_order(face_id, field_name);
323 }
324
325 template <typename FACE_SHAPE>
326 void face_geom_for_face(
327 const int face_id,
328 const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
329 b2linalg::Vector<double>& geom,
330 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_geom_d_icoor) {
331 const int* edge_node = SHAPE::face_node[face_id - 1];
332 double node_coor[FACE_SHAPE::num_nodes][3];
333 for (int k = 0; k != FACE_SHAPE::num_nodes; ++k) {
334 node_type::get_coor_s(node_coor[k], nodes[edge_node[k]]);
335 }
336 if (!geom.is_null()) {
337 b2linalg::Vector<double> shape(FACE_SHAPE::num_nodes);
338 FACE_SHAPE::eval_h(&internal_coor[0], shape);
339 geom.resize(3);
340 for (size_t i = 0; i != 3; ++i) {
341 double tmp = 0;
342 for (int k = 0; k != FACE_SHAPE::num_nodes; ++k) {
343 tmp += shape[k] * node_coor[k][i];
344 }
345 geom[i] = tmp;
346 }
347 }
348 if (!d_geom_d_icoor.is_null()) {
349 b2linalg::Matrix<double> d_shape(2, FACE_SHAPE::num_nodes);
350 FACE_SHAPE::eval_h_derived(&internal_coor[0], d_shape);
351 d_geom_d_icoor.resize(3, 2);
352 for (size_t i = 0; i != 3; ++i) {
353 for (size_t j = 0; j != 2; ++j) {
354 double tmp = 0;
355 for (int k = 0; k != FACE_SHAPE::num_nodes; ++k) {
356 tmp += d_shape(j, k) * node_coor[k][i];
357 }
358 d_geom_d_icoor(i, j) = tmp;
359 }
360 }
361 }
362 }
363
364 void face_geom(
365 const int face_id,
366 const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
367 b2linalg::Vector<double>& geom,
368 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_geom_d_icoor) override {
369 switch (SHAPE::face_type_index[face_id - 1]) {
370 case 0:
371 face_geom_for_face<typename SHAPE::face_shape_type_0>(
372 face_id, internal_coor, geom, d_geom_d_icoor);
373 break;
374 case 1:
375 face_geom_for_face<typename SHAPE::face_shape_type_1>(
376 face_id, internal_coor, geom, d_geom_d_icoor);
377 break;
378 }
379 }
380
381 template <typename FACE_SHAPE>
382 void face_field_value_for_face(
383 const int face_id, const std::string& field_name,
384 const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
385 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof, const double time,
386 b2linalg::Vector<double, b2linalg::Vdense>& value,
387 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_value_d_icoor,
388 b2linalg::Index& dof_numbering,
389 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_value_d_dof,
390 b2linalg::Index& d_value_d_dof_dep_col) {
391 const int* edge_node = SHAPE::face_node[face_id - 1];
392 size_t dn[FACE_SHAPE::num_nodes];
393 {
394 size_t* ptr = &dn[0];
395 for (int k = 0; k != FACE_SHAPE::num_nodes; ++k) {
396 ptr = node_type::get_global_dof_numbering_s(ptr, nodes[edge_node[k]]);
397 }
398 }
399 if (!value.is_null()) {
400 b2linalg::Vector<double> shape(FACE_SHAPE::num_nodes);
401 FACE_SHAPE::eval_h(&internal_coor[0], shape);
402 value.resize(1);
403 double& tmp = value[0];
404 tmp = 0;
405 if (!dof.is_null()) {
406 for (int k = 0; k != FACE_SHAPE::num_nodes; ++k) {
407 tmp += shape[k] * dof(0, dn[k]);
408 }
409 }
410 }
411 if (!d_value_d_icoor.is_null()) {
412 b2linalg::Matrix<double> d_shape(2, FACE_SHAPE::num_nodes);
413 FACE_SHAPE::eval_h_derived(&internal_coor[0], d_shape);
414 d_value_d_icoor.resize(1, 2);
415 for (int i = 0; i != 2; ++i) {
416 double& tmp = d_value_d_icoor(0, i);
417 tmp = 0;
418 if (!dof.is_null()) {
419 for (int k = 0; k != FACE_SHAPE::num_nodes; ++k) {
420 tmp += d_shape(i, k) * dof(0, dn[k]);
421 }
422 }
423 }
424 }
425 if (!d_value_d_dof.is_null()) {
426 dof_numbering.resize(FACE_SHAPE::num_nodes);
427 std::copy(dn, dn + FACE_SHAPE::num_nodes, dof_numbering.begin());
428 b2linalg::Vector<double> shape(FACE_SHAPE::num_nodes);
429 FACE_SHAPE::eval_h(&internal_coor[0], shape);
430 d_value_d_dof.resize(1, FACE_SHAPE::num_nodes);
431 d_value_d_dof.set_zero();
432 for (int k = 0; k != FACE_SHAPE::num_nodes; ++k) { d_value_d_dof(0, k) = shape[k]; }
433 }
434 }
435
436 void face_field_value(
437 const int face_id, const std::string& field_name,
438 const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
439 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof, const double time,
440 b2linalg::Vector<double, b2linalg::Vdense>& value,
441 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_value_d_icoor,
442 b2linalg::Index& dof_numbering,
443 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_value_d_dof,
444 b2linalg::Index& d_value_d_dof_dep_col) override {
445 switch (SHAPE::face_type_index[face_id - 1]) {
446 case 0:
447 face_field_value_for_face<typename SHAPE::face_shape_type_0>(
448 face_id, field_name, internal_coor, dof, time, value, d_value_d_icoor,
449 dof_numbering, d_value_d_dof, d_value_d_dof_dep_col);
450 break;
451 case 1:
452 face_field_value_for_face<typename SHAPE::face_shape_type_1>(
453 face_id, field_name, internal_coor, dof, time, value, d_value_d_icoor,
454 dof_numbering, d_value_d_dof, d_value_d_dof_dep_col);
455 break;
456 }
457 }
458
459 void field_volume_integration(
460 const VDC& dof, const VDC& dofdot, const VDC& dofdotdot, const Field<double>& field,
461 b2linalg::Index& dof_numbering, VD& discretised_field,
462 MP& d_discretised_field_d_dof) override {
463 // First expand and initialize vector of concentrated
464 // forces.
465 if (!discretised_field.is_null()) {
466 discretised_field.resize(N_t::nb_nodes);
467 discretised_field.set_zero();
468 }
469
470 // Set the dof numbering scheme.
471 get_dof_numbering(dof_numbering);
472
473 // Get the nodes coordinates.
474 double node_coor[N_t::nb_nodes][3];
475 for (int k = 0; k != N_t::nb_nodes; ++k) { node_type::get_coor_s(node_coor[k], nodes[k]); }
476
477 // Loop over all volume integration points.
478 for (int ng = 0; ng != N_t::nb_gauss; ++ng) {
479 const typename N_t::OnePoint& N = N_t::array[ng];
480
481 // Compute the Jacobian of the natural coordinate
482 // J[j][i] = d x[j] / d r[i]
483 double J[nb_dimensions][nb_dimensions];
484 for (int j = 0; j != nb_dimensions; ++j) {
485 for (int i = 0; i != nb_dimensions; ++i) {
486 double v = 0;
487 for (int n = 0; n != N_t::nb_nodes; ++n) { v += N.dN[n][i] * node_coor[n][j]; }
488 J[j][i] = v;
489 }
490 }
491
492 // Compute the determinant and the inverse of J.
493 double inv_J[nb_dimensions][nb_dimensions];
494 double volume = invert_x_x(J, inv_J) * N.weight;
495 if (nb_dimensions == 2) {
496 double thickness;
497 if (AXISYMMETRIC) {
498 thickness = 0;
499 for (int n = 0; n != N_t::nb_nodes; ++n) {
500 thickness += N.N[n] * node_coor[n][1];
501 }
502 thickness = std::abs(thickness) * 2 * M_PIl;
503 } else {
504 thickness = HeatMaterialType<nb_dimensions>::get_thickness(
505 property, this,
506 b2linalg::Vector<double, b2linalg::Vdense_constref>(N_t::nb_nodes, N.N));
507 }
508 volume *= thickness;
509 }
510 if (volume <= 0) {
511 Exception() << "Negative Jacobian, element " << this->get_object_name() << THROW;
512 }
513
514 // The element-local coordinates.
515 const double lcoor[3] = {
516 N.IntCoor[0], N.IntCoor[1], nb_dimensions > 2 ? N.IntCoor[2] : 0};
517
518 // Interpolate coordinates of integration point.
519 double coor[nb_dimensions];
520 for (int i = 0; i != nb_dimensions; ++i) {
521 double v = 0;
522 for (int n = 0; n != N_t::nb_nodes; ++n) { v += N.N[n] * node_coor[n][i]; }
523 coor[i] = v;
524 }
525
526 // Obtain the value from the Field object.
527 b2linalg::Vector<double, b2linalg::Vdense> value;
528 field.get_value(
529 b2linalg::Vector<double, b2linalg::Vdense_constref>(3, lcoor),
530 b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_dimensions, coor),
531 b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_dimensions, coor),
532 b2linalg::Matrix<double, b2linalg::Mrectangle_constref>(
533 nb_dimensions, nb_dimensions, &J[0][0]),
534 b2linalg::Matrix<double, b2linalg::Mrectangle_constref>(
535 nb_dimensions, nb_dimensions, &J[0][0]),
536 value);
537 if (value.size() != 1) { UnimplementedError() << THROW; }
538 if (field.get_multiply_with_density()) { UnimplementedError() << THROW; }
539
540 if (!discretised_field.is_null()) {
541 // Distribute value to nodes (discretize).
542 for (int k = 0; k != N_t::nb_nodes; ++k) {
543 discretised_field[k] += volume * value[0] * N.N[k];
544 }
545 }
546 }
547 }
548
549 static type_t type;
550
551private:
552 // The nodes list.
553 Node* nodes[N_t::nb_nodes];
554
555 // The property
556 heat_material_type* property;
557
558 static N_t internal_n_t;
559};
560
561template <typename SHAPE, typename INTEGRATION, int NB_GAUSS, bool AXISYMMETRIC>
562typename ElementHeatConduction<SHAPE, INTEGRATION, NB_GAUSS, AXISYMMETRIC>::type_t
563 ElementHeatConduction<SHAPE, INTEGRATION, NB_GAUSS, AXISYMMETRIC>::type(
564 std::string(SHAPE::name) + ".HEAT.CONDUCTION"
565 + (SHAPE::num_dimensions == 2 ? (AXISYMMETRIC ? ".AXISYMMETRIC" : ".2D") : ""),
566 "", StringList(), element::module, &TypedElement<double>::type);
567
568template <typename SHAPE, typename INTEGRATION, int NB_GAUSS, bool AXISYMMETRIC>
569typename ElementHeatConduction<SHAPE, INTEGRATION, NB_GAUSS, AXISYMMETRIC>::N_t
570 ElementHeatConduction<SHAPE, INTEGRATION, NB_GAUSS, AXISYMMETRIC>::internal_n_t;
571
572class ElementHeatRadConvBase : public TypedElement<double> {
573public:
574 virtual bool has_radiation_property() const = 0;
575
576 virtual void get_radiation_property_at_internal_coor(
577 Model& model, const double coordinates[3], const double internal_coor[2],
578 const double temperature, double& heat_radiation_unit, double& diffuse_reflectivity,
579 b2linalg::Vector<double, b2linalg::Vdense>& dof_interpolation) = 0;
580};
581
582template <typename SHAPE, typename INTEGRATION, int NB_GAUSS, bool AXISYMMETRIC = false>
583class ElementHeatRadConv : public ElementHeatRadConvBase {
584public:
585 typedef GenericInterpolation<SHAPE, INTEGRATION, NB_GAUSS> N_t;
586 typedef node::HNode<
587 coordof::Coor<coordof::Trans<coordof::X, coordof::Y, coordof::Z>>,
588 coordof::Dof<coordof::Temperature>>
589 node_type;
590 static const int nb_dimensions = N_t::nb_dimensions;
591 typedef typename HeatMaterialType<nb_dimensions + 1>::heat_material_type heat_material_type;
592
593 void set_nodes(std::pair<int, Node* const*> nodes_) override {
594 if (nodes_.first != N_t::nb_nodes) {
595 Exception() << "This element has " << N_t::nb_nodes << " nodes." << THROW;
596 }
597 for (int i = 0; i != N_t::nb_nodes; ++i) { nodes[i] = nodes_.second[i]; }
598 }
599
600 std::pair<int, Node* const*> get_nodes() const override {
601 return std::pair<int, Node* const*>(N_t::nb_nodes, nodes);
602 }
603
604 void set_property(ElementProperty* property_) override {
605 property = dynamic_cast<heat_material_type*>(property_);
606 if (property == 0) {
607 TypeError() << "Bad property type " << typeid(*property_)
608 << " for heat element (forgot MID element option?)." << THROW;
609 }
610 }
611
612 const ElementProperty* get_property() const override { return property; }
613
614 void get_dof_numbering(b2linalg::Index& dof_numbering) override {
615 dof_numbering.resize(N_t::nb_nodes);
616 size_t* ptr = &dof_numbering[0];
617 for (int k = 0; k != N_t::nb_nodes; ++k) {
618 ptr = node_type::get_global_dof_numbering_s(ptr, nodes[k]);
619 }
620 }
621
622 const std::vector<VariableInfo> get_value_info() const override {
623 return std::vector<VariableInfo>(1, Element::nonlinear);
624 }
625
626 void get_value(
627 Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
628 const double time, const double delta_time,
629 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof,
630 GradientContainer* gradient_container, SolverHints* solver_hints,
631 b2linalg::Index& dof_numbering, b2linalg::Vector<double, b2linalg::Vdense>& value,
632 const std::vector<bool>& d_value_d_dof_flags,
633 std::vector<b2linalg::Matrix<double, b2linalg::Mpacked>>& d_value_d_dof,
634 b2linalg::Vector<double, b2linalg::Vdense>& d_value_d_time) override;
635
636 bool has_radiation_property() const override;
637
638 void get_radiation_property_at_internal_coor(
639 Model& model, const double coordinates[3], const double internal_coor[nb_dimensions],
640 const double temperature, double& heat_radiation_unit, double& diffuse_reflectivity,
641 b2linalg::Vector<double, b2linalg::Vdense>& dof_interpolation) override;
642
643 using type_t = ObjectTypeComplete<ElementHeatRadConv, typename TypedElement<double>::type_t>;
644
645 static type_t type;
646
647private:
648 // The nodes list.
649 Node* nodes[N_t::nb_nodes];
650
651 // The property
652 heat_material_type* property;
653
654 static N_t internal_n_t;
655};
656
657template <typename SHAPE, typename INTEGRATION, int NB_GAUSS, bool AXISYMMETRIC>
658typename ElementHeatRadConv<SHAPE, INTEGRATION, NB_GAUSS, AXISYMMETRIC>::type_t
659 ElementHeatRadConv<SHAPE, INTEGRATION, NB_GAUSS, AXISYMMETRIC>::type(
660 std::string(SHAPE::name) + ".HEAT.RADCONV"
661 + (SHAPE::num_dimensions == 1 ? (AXISYMMETRIC ? ".AXISYMMETRIC" : ".2D") : ""),
662 "", StringList(), element::module, &TypedElement<double>::type);
663
664template <typename SHAPE, typename INTEGRATION, int NB_GAUSS, bool AXISYMMETRIC>
665typename ElementHeatRadConv<SHAPE, INTEGRATION, NB_GAUSS, AXISYMMETRIC>::N_t
666 ElementHeatRadConv<SHAPE, INTEGRATION, NB_GAUSS, AXISYMMETRIC>::internal_n_t;
667
668} // namespace b2000
669
670template <typename SHAPE, typename INTEGRATION, int NB_GAUSS, bool AXISYMMETRIC>
671void b2000::ElementHeatConduction<SHAPE, INTEGRATION, NB_GAUSS, AXISYMMETRIC>::get_value(
672 Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
673 const double time, const double delta_time,
674 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof,
675 GradientContainer* gradient_container, SolverHints* solver_hints,
676 b2linalg::Index& dof_numbering, b2linalg::Vector<double, b2linalg::Vdense>& value,
677 const std::vector<bool>& d_value_d_dof_flags,
678 std::vector<b2linalg::Matrix<double, b2linalg::Mpacked>>& d_value_d_dof,
679 b2linalg::Vector<double, b2linalg::Vdense>& d_value_d_time) {
680 get_dof_numbering(dof_numbering);
681 if (!value.is_null()) {
682 value.resize(dof_numbering.size());
683 value.set_zero();
684 }
685 for (size_t i = 0; i != d_value_d_dof_flags.size(); ++i) {
686 if (d_value_d_dof_flags[i]) {
687 d_value_d_dof[i].resize(dof_numbering.size(), true);
688 d_value_d_dof[i].set_zero();
689 }
690 }
691 if (!d_value_d_time.is_null()) {
692 d_value_d_time.resize(dof_numbering.size());
693 d_value_d_time.set_zero();
694 }
695 const bool compute_d_value_d_dof = (d_value_d_dof_flags.size() > 0 && d_value_d_dof_flags[0]);
696 const bool compute_d_value_d_dofdot =
697 (d_value_d_dof_flags.size() > 1 && d_value_d_dof_flags[1]);
698
699 // Get the nodes coordinate
700 double NodeCoor[N_t::nb_nodes][3];
701 for (int k = 0; k != N_t::nb_nodes; ++k) { node_type::get_coor_s(NodeCoor[k], nodes[k]); }
702
703#ifdef LUMPED_CAPACITY
704 double volume = 0;
705#endif
706 // Iterate through the integration point
707 for (int ng = 0; ng != N_t::nb_gauss; ++ng) {
708 const typename N_t::OnePoint& N = N_t::array[ng];
709
710 // Compute the Jacobian of the natural coordinate
711 // J[j][i] = d x[j] / d r[i]
712 double J[nb_dimensions][nb_dimensions];
713 for (int j = 0; j != nb_dimensions; ++j) {
714 for (int i = 0; i != nb_dimensions; ++i) {
715 double v = 0;
716 for (int n = 0; n != N_t::nb_nodes; ++n) { v += N.dN[n][i] * NodeCoor[n][j]; }
717 J[j][i] = v;
718 }
719 }
720
721 // Compute the determinant and the inverse of J
722 double inv_J[nb_dimensions][nb_dimensions];
723 const double weight = invert_x_x(J, inv_J) * N.weight;
724#ifdef LUMPED_CAPACITY
725 volume += weight;
726#endif
727 if (weight <= 0) {
728 Exception() << "The initial configuration of the element " << this->get_object_name()
729 << " is invalid (negative Jacobian determinant)." << THROW;
730 }
731
732 // Compute the gradient of the interpolation
733 // B[j][i] = d N[j] / d x[i]
734 double B[N_t::nb_nodes][nb_dimensions];
735 blas::gemm(
736 'N', 'N', nb_dimensions, N_t::nb_nodes, nb_dimensions, 1.0, inv_J[0], nb_dimensions,
737 N.dN[0], nb_dimensions, 0.0, B[0], nb_dimensions);
738
739 double coor[3] = {};
740 for (int i = 0; i != nb_dimensions; ++i) {
741 double v = 0;
742 for (int n = 0; n != N_t::nb_nodes; ++n) { v += N.N[n] * NodeCoor[n][i]; }
743 coor[i] = v;
744 }
745 if (gradient_container) {
746 const GradientContainer::InternalElementPosition ip = {
747 N.IntCoor[0], N.IntCoor[1], nb_dimensions > 2 ? N.IntCoor[2] : 0, 0};
748 gradient_container->set_current_position(ip, weight);
749
750 static const GradientContainer::FieldDescription coor_descr = {
751 "COOR_IP",
752 nb_dimensions > 2 ? "x.y.z" : "x.y",
753 "Physical coordinate of the point "
754 "in the reference configuration",
755 nb_dimensions,
756 nb_dimensions,
757 1,
758 nb_dimensions,
759 false,
760 typeid(double)};
761 if (gradient_container->compute_field_value(coor_descr.name)) {
762 gradient_container->set_field_value(coor_descr, coor);
763 }
764 }
765
766 double heat_conduction[nb_dimensions];
767 double d_heat_conduction_d_temperature[nb_dimensions];
768 double d_heat_conduction_d_grad_temperature[(nb_dimensions * (nb_dimensions + 1)) / 2];
769 double heat_capacity;
770 double d_heat_capacity_d_temperature;
771 double d_heat_capacity_d_temperature_dot;
772 if (dof.size2() == 0) {
773 property->get_conduction_and_capacity_heat(
774 &model, this, coor,
775 b2linalg::Vector<double, b2linalg::Vdense_constref>(N_t::nb_nodes, N.N),
776 N.IntCoor, 0, J, 0, 0, 0, time, linear, heat_conduction,
777 (compute_d_value_d_dof ? d_heat_conduction_d_temperature : 0),
778 (compute_d_value_d_dof ? d_heat_conduction_d_grad_temperature : 0), heat_capacity,
779 d_heat_capacity_d_temperature, d_heat_capacity_d_temperature_dot,
780 gradient_container, solver_hints);
781 } else {
782 double temperature = 0;
783 double grad_temperature[nb_dimensions];
784 std::fill_n(grad_temperature, nb_dimensions, 0.0);
785 for (int n = 0; n != N_t::nb_nodes; ++n) {
786 const double v = dof(dof_numbering[n], 0);
787 temperature += v * N.N[n];
788 for (int i = 0; i != nb_dimensions; ++i) { grad_temperature[i] += v * B[n][i]; }
789 }
790 if (dof.size2() == 1) {
791 property->get_conduction_and_capacity_heat(
792 &model, this, coor,
793 b2linalg::Vector<double, b2linalg::Vdense_constref>(N_t::nb_nodes, N.N),
794 N.IntCoor, 0, J, temperature, 0, grad_temperature, time, linear,
795 heat_conduction,
796 (compute_d_value_d_dof ? d_heat_conduction_d_temperature : 0),
797 (compute_d_value_d_dof ? d_heat_conduction_d_grad_temperature : 0),
798 heat_capacity, d_heat_capacity_d_temperature,
799 d_heat_capacity_d_temperature_dot, gradient_container, solver_hints);
800 } else {
801 double temperature_dot = 0;
802 for (int n = 0; n != N_t::nb_nodes; ++n) {
803 temperature_dot += dof(dof_numbering[n], 1) * N.N[n];
804 }
805 property->get_conduction_and_capacity_heat(
806 &model, this, coor,
807 b2linalg::Vector<double, b2linalg::Vdense_constref>(N_t::nb_nodes, N.N),
808 N.IntCoor, 0, J, temperature, temperature_dot, grad_temperature, time, linear,
809 heat_conduction,
810 (compute_d_value_d_dof ? d_heat_conduction_d_temperature : 0),
811 (compute_d_value_d_dof ? d_heat_conduction_d_grad_temperature : 0),
812 heat_capacity, d_heat_capacity_d_temperature,
813 d_heat_capacity_d_temperature_dot, gradient_container, solver_hints);
814 }
815 }
816
817 if (nb_dimensions == 2) {
818 double thickness;
819 if (AXISYMMETRIC) {
820 thickness = 0;
821 for (int n = 0; n != N_t::nb_nodes; ++n) { thickness += N.N[n] * NodeCoor[n][1]; }
822 thickness = std::abs(thickness) * 2 * M_PIl;
823 } else {
824 thickness = HeatMaterialType<nb_dimensions>::get_thickness(
825 property, this,
826 b2linalg::Vector<double, b2linalg::Vdense_constref>(N_t::nb_nodes, N.N));
827 }
828 heat_conduction[0] *= thickness;
829 heat_conduction[1] *= thickness;
830 if (compute_d_value_d_dof) {
831 d_heat_conduction_d_temperature[0] *= thickness;
832 d_heat_conduction_d_temperature[1] *= thickness;
833 d_heat_conduction_d_grad_temperature[0] *= thickness;
834 d_heat_conduction_d_grad_temperature[1] *= thickness;
835 d_heat_conduction_d_grad_temperature[2] *= thickness;
836 }
837 heat_capacity *= thickness;
838 d_heat_capacity_d_temperature *= thickness;
839 d_heat_capacity_d_temperature_dot *= thickness;
840 }
841
842 if (!value.is_null()) {
843 for (int n = 0; n != N_t::nb_nodes; ++n) {
844 double tmp = 0;
845 for (int i = 0; i != nb_dimensions; ++i) { tmp += B[n][i] * heat_conduction[i]; }
846 value[n] += weight
847 * (tmp
848#ifndef LUMPED_CAPACITY
849 + N.N[n] * heat_capacity
850#endif
851 );
852 }
853 }
854
855 if (compute_d_value_d_dof) {
856 double tmp[N_t::nb_nodes];
857 for (int i = 0; i != N_t::nb_nodes; ++i) {
858 tmp[i] = B[i][0] * d_heat_conduction_d_temperature[0]
859 + B[i][1] * d_heat_conduction_d_temperature[1];
860 if (nb_dimensions == 3) { tmp[i] += B[i][2] * d_heat_conduction_d_temperature[2]; }
861 }
862 double* ptr = &d_value_d_dof[0](0, 0);
863 for (int i = 0; i != N_t::nb_nodes; ++i) {
864 if (nb_dimensions == 3) {
865 const double a0 = d_heat_conduction_d_grad_temperature[0] * B[i][0]
866 + d_heat_conduction_d_grad_temperature[1] * B[i][1]
867 + d_heat_conduction_d_grad_temperature[2] * B[i][2];
868 const double a1 = d_heat_conduction_d_grad_temperature[1] * B[i][0]
869 + d_heat_conduction_d_grad_temperature[3] * B[i][1]
870 + d_heat_conduction_d_grad_temperature[4] * B[i][2];
871 const double a2 = d_heat_conduction_d_grad_temperature[2] * B[i][0]
872 + d_heat_conduction_d_grad_temperature[4] * B[i][1]
873 + d_heat_conduction_d_grad_temperature[5] * B[i][2];
874 for (int j = 0; j <= i; ++j, ++ptr) {
875 *ptr += weight
876 * (B[j][0] * a0 + B[j][1] * a1 + B[j][2] * a2
877 + 0.5 * (tmp[i] * N.N[j] + tmp[j] * N.N[i])
878#ifndef LUMPED_CAPACITY
879 + d_heat_capacity_d_temperature * N.N[i] * N.N[j]
880#endif
881 );
882 }
883 } else {
884 const double a0 = d_heat_conduction_d_grad_temperature[0] * B[i][0]
885 + d_heat_conduction_d_grad_temperature[1] * B[i][1];
886 const double a1 = d_heat_conduction_d_grad_temperature[1] * B[i][0]
887 + d_heat_conduction_d_grad_temperature[2] * B[i][1];
888 for (int j = 0; j <= i; ++j, ++ptr) {
889 *ptr += weight
890 * (B[j][0] * a0 + B[j][1] * a1
891 + 0.5 * (tmp[i] * N.N[j] + tmp[j] * N.N[i])
892#ifndef LUMPED_CAPACITY
893 + d_heat_capacity_d_temperature * N.N[i] * N.N[j]
894#endif
895 );
896 }
897 }
898 }
899 }
900
901#ifndef LUMPED_CAPACITY
902 if (compute_d_value_d_dofdot) {
903 double* ptr = &d_value_d_dof[1](0, 0);
904 for (int i = 0; i != N_t::nb_nodes; ++i) {
905 for (int j = 0; j <= i; ++j, ++ptr) {
906 *ptr += weight * (N.N[i] * N.N[j] * d_heat_capacity_d_temperature_dot);
907 }
908 }
909 }
910#endif
911 }
912
913#ifdef LUMPED_CAPACITY
914 volume /= N_t::nb_nodes;
915 double* ptr = &d_value_d_dof[1](0, 0);
916 for (int n = 0; n != N_t::nb_nodes; ++n) {
917 double temperature = 0;
918 double temperature_dot = 0;
919 if (dof.size2() != 0) {
920 temperature = dof(dof_numbering[n], 0);
921 if (dof.size2() > 1) { temperature_dot = dof(dof_numbering[n], 1); }
922 }
923 double heat_capacity;
924 double d_heat_capacity_d_temperature;
925 double d_heat_capacity_d_temperature_dot;
926 double N[N_t::nb_nodes];
927 std::fill_n(N, N_t::nb_nodes, 0);
928 N[n] = 1;
929 double IntCoor[3] = {};
930 double J[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
931 property->get_conduction_and_capacity_heat(
932 &model, this, Vector<double, Vdense_constref>(N_t::nb_nodes, N), IntCoor, 0, J,
933 temperature, temperature_dot, 0, time, linear, 0, 0, 0, heat_capacity,
934 d_heat_capacity_d_temperature, d_heat_capacity_d_temperature_dot, 0);
935 if (!value.is_null()) { value[n] += volume * heat_capacity; }
936 if (compute_d_value_d_dofdot) {
937 *ptr += volume * d_heat_capacity_d_temperature_dot;
938 ptr += n + 2;
939 }
940 }
941#endif
942}
943
944template <typename SHAPE, typename INTEGRATION, int NB_GAUSS, bool AXISYMMETRIC>
945void b2000::ElementHeatRadConv<SHAPE, INTEGRATION, NB_GAUSS, AXISYMMETRIC>::get_value(
946 Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
947 const double time, const double delta_time,
948 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof,
949 GradientContainer* gradient_container, SolverHints* solver_hints,
950 b2linalg::Index& dof_numbering, b2linalg::Vector<double, b2linalg::Vdense>& value,
951 const std::vector<bool>& d_value_d_dof_flags,
952 std::vector<b2linalg::Matrix<double, b2linalg::Mpacked>>& d_value_d_dof,
953 b2linalg::Vector<double, b2linalg::Vdense>& d_value_d_time) {
954 get_dof_numbering(dof_numbering);
955 if (!value.is_null()) {
956 value.resize(dof_numbering.size());
957 value.set_zero();
958 }
959 for (size_t i = 0; i != d_value_d_dof_flags.size(); ++i) {
960 if (d_value_d_dof_flags[i]) {
961 d_value_d_dof[i].resize(dof_numbering.size(), true);
962 d_value_d_dof[i].set_zero();
963 }
964 }
965 if (!d_value_d_time.is_null()) {
966 d_value_d_time.resize(dof_numbering.size());
967 d_value_d_time.set_zero();
968 }
969 const bool compute_d_value_d_dof = (d_value_d_dof_flags.size() > 0 && d_value_d_dof_flags[0]);
970
971 // Get the nodes coordinate
972 double NodeCoor[N_t::nb_nodes][3];
973 for (int k = 0; k != N_t::nb_nodes; ++k) { node_type::get_coor_s(NodeCoor[k], nodes[k]); }
974
975 // Iterate through the integration point
976 for (int ng = 0; ng != N_t::nb_gauss; ++ng) {
977 const typename N_t::OnePoint& N = N_t::array[ng];
978
979 double weight;
980 {
981 // Compute the Jacobian of the natural coordinate
982 // J[j][i] = d x[j] / d r[i]
983 double J[nb_dimensions + 1][nb_dimensions];
984 for (int j = 0; j != nb_dimensions + 1; ++j) {
985 for (int i = 0; i != nb_dimensions; ++i) {
986 double v = 0;
987 for (int n = 0; n != N_t::nb_nodes; ++n) { v += N.dN[n][i] * NodeCoor[n][j]; }
988 J[j][i] = v;
989 }
990 }
991 if (nb_dimensions == 2) {
992 const double outer[3] = {
993 J[1][0] * J[2][1] - J[2][0] * J[1][1], J[2][0] * J[0][1] - J[0][0] * J[2][1],
994 J[0][0] * J[1][1] - J[1][0] * J[0][1]};
995
996 weight = std::sqrt(outer[0] * outer[0] + outer[1] * outer[1] + outer[2] * outer[2])
997 * N.weight;
998 } else {
999 weight = std::sqrt(J[0][0] * J[0][0] + J[1][0] * J[1][0]) * N.weight;
1000 }
1001 }
1002
1003 double coor[3] = {};
1004 for (int i = 0; i != nb_dimensions + 1; ++i) {
1005 double v = 0;
1006 for (int n = 0; n != N_t::nb_nodes; ++n) { v += N.N[n] * NodeCoor[n][i]; }
1007 coor[i] = v;
1008 }
1009 if (gradient_container) {
1010 const GradientContainer::InternalElementPosition ip = {
1011 N.IntCoor[0], nb_dimensions > 1 ? N.IntCoor[1] : 0, 0, 0};
1012 gradient_container->set_current_position(ip, weight);
1013
1014 static const GradientContainer::FieldDescription coor_descr = {
1015 "COOR_IP",
1016 nb_dimensions > 1 ? "x.y.z" : "x.y",
1017 "Physical coordinate of the point "
1018 "in the reference configuration",
1019 nb_dimensions + 1,
1020 nb_dimensions + 1,
1021 1,
1022 nb_dimensions + 1,
1023 false,
1024 typeid(double)};
1025 if (gradient_container->compute_field_value(coor_descr.name)) {
1026 gradient_container->set_field_value(coor_descr, coor);
1027 }
1028 }
1029
1030 double heat;
1031 double d_heat_d_temperature;
1032 double temperature = 0;
1033 if (dof.size2() == 0) {
1034 property->get_radiation_and_convection_heat(
1035 &model, this, coor,
1036 b2linalg::Vector<double, b2linalg::Vdense_constref>(N_t::nb_nodes, N.N),
1037 N.IntCoor, 0, time, linear, heat, d_heat_d_temperature, gradient_container,
1038 solver_hints);
1039 } else {
1040 for (int n = 0; n != N_t::nb_nodes; ++n) {
1041 temperature += dof(dof_numbering[n], 0) * N.N[n];
1042 }
1043 property->get_radiation_and_convection_heat(
1044 &model, this, coor,
1045 b2linalg::Vector<double, b2linalg::Vdense_constref>(N_t::nb_nodes, N.N),
1046 N.IntCoor, temperature, time, linear, heat, d_heat_d_temperature,
1047 gradient_container, solver_hints);
1048 }
1049
1050 if (nb_dimensions == 1) {
1051 double thickness;
1052 if (AXISYMMETRIC) {
1053 thickness = 0;
1054 for (int n = 0; n != N_t::nb_nodes; ++n) { thickness += N.N[n] * NodeCoor[n][1]; }
1055 thickness = std::abs(thickness) * 2 * M_PIl;
1056 } else {
1057 thickness = HeatMaterialType<nb_dimensions + 1>::get_thickness(
1058 property, this,
1059 b2linalg::Vector<double, b2linalg::Vdense_constref>(N_t::nb_nodes, N.N));
1060 }
1061 heat *= thickness;
1062 d_heat_d_temperature *= thickness;
1063 }
1064
1065 if (!value.is_null()) {
1066 for (int n = 0; n != N_t::nb_nodes; ++n) { value[n] += weight * N.N[n] * heat; }
1067 }
1068
1069 if (compute_d_value_d_dof) {
1070 double* ptr = &d_value_d_dof[0](0, 0);
1071 for (int i = 0; i != N_t::nb_nodes; ++i) {
1072 for (int j = 0; j <= i; ++j, ++ptr) {
1073 *ptr += weight * d_heat_d_temperature * N.N[i] * N.N[j];
1074 }
1075 }
1076 }
1077 }
1078}
1079
1080template <typename SHAPE, typename INTEGRATION, int NB_GAUSS, bool AXISYMMETRIC>
1081bool b2000::ElementHeatRadConv<SHAPE, INTEGRATION, NB_GAUSS, AXISYMMETRIC>::has_radiation_property()
1082 const {
1083 double emisivity, receptivity;
1084 return property->get_radiation_properties(
1085 0, this, 0, b2linalg::Vector<double, b2linalg::Vdense_constref>::null, 0, 0, emisivity,
1086 receptivity);
1087}
1088
1089template <typename SHAPE, typename INTEGRATION, int NB_GAUSS, bool AXISYMMETRIC>
1090void b2000::ElementHeatRadConv<SHAPE, INTEGRATION, NB_GAUSS, AXISYMMETRIC>::
1091 get_radiation_property_at_internal_coor(
1092 Model& model, const double coordinates[3], const double internal_coor[nb_dimensions],
1093 const double temperature, double& heat_radiation_unit, double& diffuse_reflectivity,
1094 b2linalg::Vector<double, b2linalg::Vdense>& dof_interpolation) {
1095 dof_interpolation.resize(N_t::nb_nodes);
1096 SHAPE::eval_h(internal_coor, dof_interpolation);
1097 property->get_radiation_properties(
1098 &model, this, coordinates, dof_interpolation, internal_coor, temperature,
1099 heat_radiation_unit, diffuse_reflectivity);
1100}
1101
1102#endif // B2_ELEMENT_HEAT_H_
THROW
#define THROW
Definition b2exception.H:198
b2tensor_calculus.H
b2000::Element::nonlinear
@ nonlinear
Definition b2element.H:619
b2000::Element::linear
@ linear
Definition b2element.H:615
b2000::Element::get_object_name
const std::string & get_object_name() const override
Definition b2element.H:220
b2000::TypedElement< double >::get_value
virtual void get_value(Model &model, const bool linear, const EquilibriumSolution equilibrium_solution, const double time, const double delta_time, const b2linalg::Matrix< double, b2linalg::Mrectangle_constref > &dof, GradientContainer *gradient_container, SolverHints *solver_hints, b2linalg::Index &dof_numbering, b2linalg::Vector< double, b2linalg::Vdense > &value, const std::vector< bool > &d_value_d_dof_flags, std::vector< b2linalg::Matrix< double, b2linalg::Mrectangle > > &d_value_d_dof, b2linalg::Vector< double, b2linalg::Vdense > &d_value_d_time)
Compute the internal forces and their derivatives of the element.
Definition b2element.H:1619
b2000::TypedElement< double >::field_volume_integration
virtual void field_volume_integration(const b2linalg::Vector< double, b2linalg::Vdense_constref > &dof, const b2linalg::Vector< double, b2linalg::Vdense_constref > &dofdot, const b2linalg::Vector< double, b2linalg::Vdense_constref > &dofdotdot, const Field< double > &f, b2linalg::Index &dof_numbering, b2linalg::Vector< double, b2linalg::Vdense > &discretised_field, b2linalg::Matrix< double, b2linalg::Mpacked > &d_discretised_field_d_dof)
Definition b2element.H:1256
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32
b2000::invert_x_x
T invert_x_x(const T a[2][2], T b[2][2])
Definition b2tensor_calculus.H:742
b2000::TypeError
GenericException< TypeError_name > TypeError
Definition b2exception.H:325
b2000::UnimplementedError
GenericException< UnimplementedError_name > UnimplementedError
Definition b2exception.H:314