b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2element_field_transfer.H
1//------------------------------------------------------------------------
2// b2element_field_transfer.H --
3//
4//
5// written by Mathias Doreille
6// Thomas Blome <thomas.blome@dlr.de>
7//
8// Copyright (c) 2009 SMR Engineering & Development SA
9// 2502 Bienne, Switzerland
10//
11// Copyright (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
12// Linder Höhe, 51147 Köln
13//
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 B2ELEMENT_FIELD_TRANSFER_H_
22#define B2ELEMENT_FIELD_TRANSFER_H_
23
24#include <cmath>
25#include <string>
26
27#include "model/b2element.H"
28#include "model_imp/b2hnode.H"
29#include "utils/b2constants.H"
30#include "utils/b2tensor_calculus.H"
31
32namespace b2000 {
33
34namespace {
35const std::string geometry("GEOMETRY");
36const std::string temperature("TEMPERATURE");
37const std::string heat_flux("HEATFLUX");
38const std::string displacement("DISPLACEMENT");
39const std::string traction("TRACTION");
40const std::string all("ALL");
41} // namespace
42
43class CFDMaterial : virtual public ElementProperty {
44public:
45 double get_internal_pressure() { return internal_pressure; }
46
47protected:
48 void set_internal_pressure(double p) { internal_pressure = p; }
49
50private:
51 double internal_pressure;
52};
53
54// The Fields are "GEOMETRY", "DISPLACEMENT", "TRACTION"
55
56class ElementCFDPressureTraction2DFieldL2 : public TypedElement<double> {
57public:
58 std::pair<int, Node* const*> get_nodes() const override {
59 return std::pair<int, Node* const*>(2, nodes);
60 }
61
62 void set_nodes(std::pair<int, Node* const*> nodes_) override {
63 if (nodes_.first != 2) { Exception() << "This element must have 2 nodes." << THROW; }
64 for (int i = 0; i != 2; ++i) { nodes[i] = nodes_.second[i]; }
65 }
66
67 int edge_field_order(const int edge_id, const std::string& field_name) override { return 1; }
68
69 bool edge_field_linear_on_dof(const int edge_id, const std::string& field_name) override {
70 if (field_name == traction) { return false; }
71 return true;
72 }
73
74 int edge_field_polynomial_sub_edge(
75 const int edge_id, const std::string& field_name,
76 b2linalg::Vector<double, b2linalg::Vdense>& sub_nodes) override {
77 sub_nodes.clear();
78 if (field_name == all || field_name == traction) {
79 sub_nodes.reserve(3);
80 sub_nodes.push_back(0);
81 sub_nodes.push_back(0.5);
82 sub_nodes.push_back(1);
83 } else {
84 sub_nodes.reserve(2);
85 sub_nodes.push_back(0);
86 sub_nodes.push_back(1);
87 }
88 return 1;
89 }
90
91 void set_property(ElementProperty* property_) override {
92 property = dynamic_cast<CFDMaterial*>(property_);
93 if (property == nullptr) {
94 TypeError() << "Bad property type " << typeid(*property_)
95 << " for ElementCFD2DFieldL2 element."
96 << "(forgot MID or PID element attribute?)" << THROW;
97 }
98 }
99
100 const ElementProperty* get_property() const override { return property; }
101
102protected:
103 Node* nodes[2];
104
105 // The property
106 CFDMaterial* property;
107};
108
109class ElementCFDPressure2DFieldL2 : public ElementCFDPressureTraction2DFieldL2 {
110public:
111 typedef node::HNode<
112 coordof::Coor<coordof::Trans<coordof::X, coordof::Y>>,
113 coordof::Dof<coordof::DTrans<coordof::DX, coordof::DY>, coordof::Pressure>>
114 node_type;
115
116 void edge_geom(
117 const int edge_id, const double internal_coor, b2linalg::Vector<double>& geom,
118 b2linalg::Vector<double>& d_geom_d_icoor) override {
119 double node_coor[2][2];
120 for (int k = 0; k != 2; ++k) { node_type::get_coor_s(node_coor[k], nodes[k]); }
121 if (!geom.is_null()) {
122 geom.resize(2);
123 for (size_t i = 0; i != 2; ++i) {
124 geom[i] = node_coor[0][i] * (1 - internal_coor) + node_coor[1][i] * internal_coor;
125 }
126 }
127 if (!d_geom_d_icoor.is_null()) {
128 d_geom_d_icoor.resize(2);
129 for (size_t i = 0; i != 2; ++i) {
130 d_geom_d_icoor[i] = -node_coor[0][i] + node_coor[1][i];
131 }
132 }
133 }
134
135 void edge_field_value(
136 const int edge_id, const std::string& field_name, const double internal_coor,
137 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof, const double time,
138 b2linalg::Vector<double, b2linalg::Vdense>& value,
139 b2linalg::Vector<double, b2linalg::Vdense>& d_value_d_icoor,
140 b2linalg::Index& dof_numbering,
141 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_value_d_dof,
142 b2linalg::Index& d_value_d_dof_dep) override {
143 double node_coor[2][2];
144 for (int k = 0; k != 2; ++k) { node_type::get_coor_s(node_coor[k], nodes[k]); }
145 size_t dn[6];
146 {
147 size_t* ptr = dn;
148 for (int k = 0; k != 2; ++k) {
149 ptr = node_type::get_global_dof_numbering_s(ptr, nodes[k]);
150 }
151 }
152 if (field_name == displacement) {
153 if (!value.is_null()) {
154 value.resize(2);
155 for (int i = 0; i != 2; ++i) {
156 value[i] =
157 dof(dn[i], 0) * (1 - internal_coor) + dof(dn[i + 3], 0) * internal_coor;
158 }
159 }
160 if (!d_value_d_icoor.is_null()) {
161 d_value_d_icoor.resize(2);
162 for (int i = 0; i != 2; ++i) {
163 d_value_d_icoor[i] = -dof(dn[i], 0) + dof(dn[i + 3], 0);
164 }
165 }
166 if (!d_value_d_dof.is_null()) {
167 dof_numbering.resize(4);
168 dof_numbering[0] = dn[0];
169 dof_numbering[1] = dn[1];
170 dof_numbering[2] = dn[3];
171 dof_numbering[3] = dn[4];
172 d_value_d_dof.resize(2, 4);
173 d_value_d_dof.set_zero();
174 d_value_d_dof(0, 0) = 1 - internal_coor;
175 d_value_d_dof(0, 2) = internal_coor;
176 d_value_d_dof(1, 1) = 1 - internal_coor;
177 d_value_d_dof(1, 3) = internal_coor;
178 }
179 }
180 if (field_name == traction) {
181 const int pos = internal_coor < 0.5 ? 0 : 3;
182 double n[2] = {node_coor[0][1] - node_coor[1][1], node_coor[1][0] - node_coor[0][0]};
183 const double inv_L = 1 / std::sqrt(n[0] * n[0] + n[1] * n[1]);
184 double s = -property->get_internal_pressure();
185 if (dof.size2()) {
186 n[0] += dof(dn[1], 0) - dof(dn[4], 0);
187 n[1] += dof(dn[3], 0) - dof(dn[0], 0);
188 s += dof(dn[2 + pos], 0);
189 }
190 n[0] *= inv_L;
191 n[1] *= inv_L;
192
193 if (!value.is_null()) {
194 value.resize(2);
195 value[0] = n[0] * s;
196 value[1] = n[1] * s;
197 }
198 if (!d_value_d_icoor.is_null()) {
199 d_value_d_icoor.resize(2);
200 d_value_d_icoor[0] = 0;
201 d_value_d_icoor[1] = 0;
202 }
203 if (!d_value_d_dof.is_null()) {
204 dof_numbering.resize(5);
205 std::copy(dn, dn + 5, dof_numbering.begin());
206 dof_numbering[2] = dn[2 + pos];
207 d_value_d_dof.resize(2, 5);
208 d_value_d_dof(0, 0) = 0;
209 d_value_d_dof(0, 1) = inv_L * s;
210 d_value_d_dof(0, 2) = n[0];
211 d_value_d_dof(0, 3) = 0;
212 d_value_d_dof(0, 4) = -inv_L * s;
213 d_value_d_dof(1, 0) = -inv_L * s;
214 d_value_d_dof(1, 1) = 0;
215 d_value_d_dof(1, 2) = n[1];
216 d_value_d_dof(1, 3) = inv_L * s;
217 d_value_d_dof(1, 4) = 0;
218 }
219 }
220 }
221 using type_t = ObjectTypeComplete<ElementCFDPressure2DFieldL2, TypedElement<double>::type_t>;
222
223 static type_t type;
224};
225
226class ElementCFDTraction2DFieldL2 : public ElementCFDPressureTraction2DFieldL2 {
227public:
228 typedef node::HNode<
229 coordof::Coor<coordof::Trans<coordof::X, coordof::Y>>,
230 coordof::Dof<
231 coordof::DTrans<coordof::DX, coordof::DY>, coordof::Pressure,
232 coordof::DTrans<coordof::TX, coordof::TY>>>
233 node_type;
234
235 void edge_geom(
236 const int edge_id, const double internal_coor, b2linalg::Vector<double>& geom,
237 b2linalg::Vector<double>& d_geom_d_icoor) override {
238 double node_coor[2][2];
239 for (int k = 0; k != 2; ++k) { node_type::get_coor_s(node_coor[k], nodes[k]); }
240 if (!geom.is_null()) {
241 geom.resize(2);
242 for (size_t i = 0; i != 2; ++i) {
243 geom[i] = node_coor[0][i] * (1 - internal_coor) + node_coor[1][i] * internal_coor;
244 }
245 }
246 if (!d_geom_d_icoor.is_null()) {
247 d_geom_d_icoor.resize(2);
248 for (size_t i = 0; i != 2; ++i) {
249 d_geom_d_icoor[i] = -node_coor[0][i] + node_coor[1][i];
250 }
251 }
252 }
253
254 void edge_field_value(
255 const int edge_id, const std::string& field_name, const double internal_coor,
256 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof, const double time,
257 b2linalg::Vector<double, b2linalg::Vdense>& value,
258 b2linalg::Vector<double, b2linalg::Vdense>& d_value_d_icoor,
259 b2linalg::Index& dof_numbering,
260 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_value_d_dof,
261 b2linalg::Index& d_value_d_dof_dep) override {
262 double node_coor[2][2];
263 for (int k = 0; k != 2; ++k) { node_type::get_coor_s(node_coor[k], nodes[k]); }
264 size_t dn[10];
265 {
266 size_t* ptr = dn;
267 for (int k = 0; k != 2; ++k) {
268 ptr = node_type::get_global_dof_numbering_s(ptr, nodes[k]);
269 }
270 }
271 if (field_name == displacement) {
272 if (!value.is_null()) {
273 value.resize(2);
274 for (int i = 0; i != 2; ++i) {
275 value[i] =
276 dof(dn[i], 0) * (1 - internal_coor) + dof(dn[i + 5], 0) * internal_coor;
277 }
278 }
279 if (!d_value_d_icoor.is_null()) {
280 d_value_d_icoor.resize(2);
281 for (int i = 0; i != 2; ++i) {
282 d_value_d_icoor[i] = -dof(dn[i], 0) + dof(dn[i + 5], 0);
283 }
284 }
285 if (!d_value_d_dof.is_null()) {
286 dof_numbering.resize(4);
287 dof_numbering[0] = dn[0];
288 dof_numbering[1] = dn[1];
289 dof_numbering[2] = dn[5];
290 dof_numbering[3] = dn[6];
291 d_value_d_dof.resize(2, 4);
292 d_value_d_dof.set_zero();
293 d_value_d_dof(0, 0) = 1 - internal_coor;
294 d_value_d_dof(0, 2) = internal_coor;
295 d_value_d_dof(1, 1) = 1 - internal_coor;
296 d_value_d_dof(1, 3) = internal_coor;
297 }
298 }
299 if (field_name == traction) {
300 const int pos = internal_coor < 0.5 ? 0 : 5;
301 const double ip = -property->get_internal_pressure();
302 if (!value.is_null()) {
303 double n[2] = {
304 node_coor[1][0] - node_coor[0][0], node_coor[1][1] - node_coor[0][1]};
305 const double inv_L = 1 / std::sqrt(n[0] * n[0] + n[1] * n[1]);
306 if (dof.size2()) {
307 n[0] += dof(dn[5], 0) - dof(dn[0], 0);
308 n[1] += dof(dn[6], 0) - dof(dn[1], 0);
309 }
310 n[0] *= inv_L;
311 n[1] *= inv_L;
312 const double s = std::sqrt(n[0] * n[0] + n[1] * n[1]) * inv_L;
313 value.resize(2);
314 if (dof.size2() == 0) {
315 value[0] = -n[1] * ip;
316 value[1] = n[0] * ip;
317 } else {
318 const double ip_ep = ip + dof(dn[2 + pos], 0);
319 value[0] = dof(dn[3 + pos], 0) * s - n[1] * ip_ep;
320 value[1] = dof(dn[4 + pos], 0) * s + n[0] * ip_ep;
321 }
322 }
323 if (!d_value_d_icoor.is_null()) {
324 d_value_d_icoor.resize(2);
325 d_value_d_icoor[0] = 0;
326 d_value_d_icoor[1] = 0;
327 }
328 if (!d_value_d_dof.is_null()) {
329 dof_numbering.resize(7);
330 std::copy(dn, dn + 7, dof_numbering.begin());
331 for (int i = 2; i != 5; ++i) { dof_numbering[i] = dn[i + pos]; }
332
333 const double tmp0 = node_coor[1][1] - node_coor[0][1];
334 const double tmp1 = node_coor[1][0] - node_coor[0][0];
335 const double tmp2 = tmp1 + dof(dn[5], 0) - dof(dn[0], 0);
336 const double tmp4 = 1 / sqrt(tmp0 * tmp0 + tmp1 * tmp1);
337 const double tmp5 = tmp0 + dof(dn[6], 0) - dof(dn[1], 0);
338 const double tmp6 = std::sqrt(tmp5 * tmp5 + tmp2 * tmp2);
339 const double tmp7 = 1 / tmp6;
340 const double tmp8 = ip + dof(dn[2 + pos], 0);
341 const double tmp9 = -tmp8 * tmp4;
342 const double tmp10 = tmp8 * tmp4;
343 const double tmp11 = tmp4 * tmp6;
344 const double t_0 = dof(dn[3 + pos], 0);
345 const double t_1 = dof(dn[4 + pos], 0);
346
347 d_value_d_dof.resize(2, 7);
348 d_value_d_dof(0, 0) = -t_0 * tmp2 * tmp4 * tmp7;
349 d_value_d_dof(0, 1) = tmp10 - t_0 * tmp5 * tmp4 * tmp7;
350 d_value_d_dof(0, 2) = -tmp5 * tmp4;
351 d_value_d_dof(0, 3) = tmp11;
352 d_value_d_dof(0, 4) = 0;
353 d_value_d_dof(0, 5) = t_0 * tmp2 * tmp4 * tmp7;
354 d_value_d_dof(0, 6) = t_0 * tmp5 * tmp4 * tmp7 + tmp9;
355 d_value_d_dof(1, 0) = tmp9 - t_1 * tmp2 * tmp4 * tmp7;
356 d_value_d_dof(1, 1) = -t_1 * tmp5 * tmp4 * tmp7;
357 d_value_d_dof(1, 2) = tmp2 * tmp4;
358 d_value_d_dof(1, 3) = 0;
359 d_value_d_dof(1, 4) = tmp11;
360 d_value_d_dof(1, 5) = t_1 * tmp2 * tmp4 * tmp7 + tmp10;
361 d_value_d_dof(1, 6) = t_1 * tmp5 * tmp4 * tmp7;
362 }
363 }
364 }
365 using type_t = ObjectTypeComplete<ElementCFDTraction2DFieldL2, TypedElement<double>::type_t>;
366
367 static type_t type;
368};
369
373template <typename SHAPE, bool EULER>
374class ElementCFD3DField : public TypedElement<double> {
375public:
376 using node_type = node::HNode<
377 coordof::Coor<coordof::Trans<coordof::X, coordof::Y, coordof::Z>>,
378 coordof::Dof<
379 coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>, coordof::Pressure,
380 coordof::DTrans<coordof::TX, coordof::TY, coordof::TZ>>>;
381
382 using node_type_euler = node::HNode<
383 coordof::Coor<coordof::Trans<coordof::X, coordof::Y, coordof::Z>>,
384 coordof::Dof<coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>, coordof::Pressure>>;
385
386 using type_t = ObjectTypeComplete<ElementCFD3DField, TypedElement<double>::type_t>;
387
388 std::pair<int, Node* const*> get_nodes() const override {
389 return std::pair<int, Node* const*>(2, nodes);
390 }
391
392 void set_nodes(std::pair<int, Node* const*> nodes_) override {
393 if (nodes_.first != SHAPE::num_nodes) {
394 Exception() << "This element must have " << SHAPE::num_nodes << " nodes." << THROW;
395 }
396 for (int i = 0; i != SHAPE::num_nodes; ++i) { nodes[i] = nodes_.second[i]; }
397 }
398
399 int face_field_order(const int face_id, const std::string& field_name) override {
400 return SHAPE::order;
401 }
402
403 bool face_field_linear_on_dof(const int edge_id, const std::string& field_name) override {
404 if (field_name == traction) { return false; }
405 return true;
406 }
407
408 int face_field_polynomial_sub_edge(
409 const int face_id, const std::string& field_name,
410 b2linalg::Matrix<double, b2linalg::Mrectangle>& sub_nodes,
411 std::vector<Element::Triangle>& sub_faces) {
412 sub_faces.clear();
413 if (SHAPE::name == "Q4") {
414 if (field_name == all || field_name == traction) {
415 sub_nodes.resize(2, 9);
416 size_t k = 0;
417 for (int j = -1; j <= 1; ++j) {
418 for (int i = -1; i <= 1; ++i, ++k) {
419 sub_nodes(0, k) = i;
420 sub_nodes(1, k) = j;
421 }
422 }
423 sub_faces.reserve(8);
424 sub_faces.push_back(Element::Triangle(0, 1, 4));
425 sub_faces.push_back(Element::Triangle(0, 4, 3));
426 sub_faces.push_back(Element::Triangle(1, 2, 5));
427 sub_faces.push_back(Element::Triangle(1, 5, 4));
428 sub_faces.push_back(Element::Triangle(3, 4, 7));
429 sub_faces.push_back(Element::Triangle(3, 7, 6));
430 sub_faces.push_back(Element::Triangle(4, 5, 8));
431 sub_faces.push_back(Element::Triangle(4, 8, 7));
432 } else {
433 sub_nodes.resize(2, 4);
434 size_t k = 0;
435 for (int j = -1; j <= 1; j += 2) {
436 for (int i = -1; i <= 1; i += 2, ++k) {
437 sub_nodes(0, k) = i;
438 sub_nodes(1, k) = j;
439 }
440 }
441 sub_faces.reserve(2);
442 sub_faces.push_back(Element::Triangle(0, 1, 3));
443 sub_faces.push_back(Element::Triangle(0, 3, 2));
444 }
445 } else if (SHAPE::name == "T3") {
446 if (field_name == all || field_name == traction) {
447 sub_nodes.resize(2, 7);
448 sub_nodes.set_zero();
449 sub_nodes(0, 1) = 0.5;
450 sub_nodes(0, 2) = 1;
451 sub_nodes(0, 3) = 0.5;
452 sub_nodes(1, 3) = 0.5;
453 sub_nodes(1, 4) = 1;
454 sub_nodes(1, 5) = 0.5;
455 sub_nodes(0, 6) = 1.0 / 3;
456 sub_nodes(1, 6) = 1.0 / 3;
457 sub_faces.reserve(6);
458 sub_faces.push_back(Element::Triangle(0, 1, 6));
459 sub_faces.push_back(Element::Triangle(1, 2, 6));
460 sub_faces.push_back(Element::Triangle(2, 3, 6));
461 sub_faces.push_back(Element::Triangle(3, 4, 6));
462 sub_faces.push_back(Element::Triangle(4, 5, 6));
463 sub_faces.push_back(Element::Triangle(5, 0, 6));
464 } else {
465 sub_nodes.resize(2, 3);
466 sub_nodes.set_zero();
467 sub_nodes(0, 1) = 1;
468 sub_nodes(1, 2) = 1;
469 sub_faces.reserve(1);
470 sub_faces.push_back(Element::Triangle(0, 1, 2));
471 }
472 }
473 return SHAPE::order;
474 }
475
476 void face_geom(
477 const int face_id,
478 const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
479 b2linalg::Vector<double>& geom,
480 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_geom_d_icoor) override {
481 double node_coor[SHAPE::num_nodes][3];
482 for (int k = 0; k != SHAPE::num_nodes; ++k) {
483 if (EULER) {
484 node_type_euler::get_coor_s(node_coor[k], nodes[k]);
485 } else {
486 node_type::get_coor_s(node_coor[k], nodes[k]);
487 }
488 }
489
490 if (!geom.is_null()) {
491 b2linalg::Vector<double> shape(SHAPE::num_nodes);
492 SHAPE::eval_h(&internal_coor[0], shape);
493 geom.resize(3);
494 for (size_t i = 0; i != 3; ++i) {
495 double tmp = 0;
496 for (int k = 0; k != SHAPE::num_nodes; ++k) { tmp += shape[k] * node_coor[k][i]; }
497 geom[i] = tmp;
498 }
499 }
500
501 if (!d_geom_d_icoor.is_null()) {
502 b2linalg::Matrix<double> d_shape(2, SHAPE::num_nodes);
503 SHAPE::eval_h_derived(&internal_coor[0], d_shape);
504 d_geom_d_icoor.resize(3, 2);
505 for (size_t i = 0; i != 3; ++i) {
506 for (size_t j = 0; j != 2; ++j) {
507 double tmp = 0;
508 for (int k = 0; k != SHAPE::num_nodes; ++k) {
509 tmp += d_shape(j, k) * node_coor[k][i];
510 }
511 d_geom_d_icoor(i, j) = tmp;
512 }
513 }
514 }
515 }
516
517 void edge_field_value(
518 const int edge_id, const std::string& field_name, const double internal_coor,
519 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof, const double time,
520 b2linalg::Vector<double, b2linalg::Vdense>& value,
521 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_value_d_icoor,
522 b2linalg::Index& dof_numbering,
523 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_value_d_dof,
524 b2linalg::Index& d_value_d_dof_dep) {
525 size_t dn[SHAPE::num_nodes * (EULER ? 4 : 7)];
526 double node_coor[SHAPE::num_nodes][3];
527 {
528 size_t* ptr = dn;
529 for (int k = 0; k != SHAPE::num_nodes; ++k) {
530 if (EULER) {
531 ptr = node_type_euler::get_global_dof_numbering_s(ptr, nodes[k]);
532 node_type_euler::get_coor_s(node_coor[k], nodes[k]);
533 } else {
534 ptr = node_type::get_global_dof_numbering_s(ptr, nodes[k]);
535 node_type::get_coor_s(node_coor[k], nodes[k]);
536 }
537 }
538 }
539
540 if (field_name == displacement) {
541 if (!value.is_null()) {
542 b2linalg::Vector<double> shape(SHAPE::num_nodes);
543 SHAPE::eval_h(internal_coor, shape);
544 value.resize(3);
545 for (size_t i = 0; i != 3; ++i) {
546 double tmp = 0;
547 for (int k = 0; k != SHAPE::num_nodes; ++k) {
548 tmp += shape[k] * dof(dn[k * (EULER ? 4 : 7) + i], 0);
549 }
550 value[i] = tmp;
551 }
552 }
553 if (!d_value_d_icoor.is_null()) {
554 b2linalg::Matrix<double> d_shape(2, SHAPE::num_nodes);
555 SHAPE::eval_h_derived(internal_coor, d_shape);
556 d_value_d_icoor.resize(3, 2);
557 for (size_t i = 0; i != 3; ++i) {
558 for (size_t j = 0; j != 2; ++j) {
559 double tmp = 0;
560 for (int k = 0; k != SHAPE::num_nodes; ++k) {
561 tmp += d_shape(j, k) * dof(dn[k * (EULER ? 4 : 7) + i], 0);
562 }
563 d_value_d_icoor(i, j) = tmp;
564 }
565 }
566 }
567 if (!d_value_d_dof.is_null()) {
568 dof_numbering.resize(SHAPE::num_nodes * 3);
569 b2linalg::Vector<double> shape(SHAPE::num_nodes);
570 SHAPE::eval_h(internal_coor, shape);
571 d_value_d_dof.resize(3, dof_numbering.size());
572 d_value_d_dof.set_zero();
573 size_t i_ptr = 0;
574 for (int k = 0; k != SHAPE::num_nodes; ++k) {
575 for (size_t i = 0; i != 3; ++i, ++i_ptr) {
576 dof_numbering[i_ptr] = dn[k * (EULER ? 4 : 7) + i];
577 d_value_d_dof(i, i_ptr + i) = shape[k];
578 }
579 }
580 }
581 }
582
583 if (field_name == traction) {
584 // search the node affected to the internal dof
585 int pos_node;
586 Exception() << THROW; // the pos_node computation is false
587 /*
588 if (SHAPE::name == "T3") {
589 if (internal_coor[0] - internal_coor[1] < 0) {
590 if (2 * internal_coor[0] + internal_coor[1] < 1)
591 pos_node = 0;
592 else
593 pos_node = 1;
594 } else {
595 if (internal_coor[0] + 2 * internal_coor[1] < 1)
596 pos_node = 0;
597 else
598 pos_node = 2;
599 }
600 } else if (SHAPE::name == "Q4")
601 pos_node = internal_coor[1] < 0 ? (internal_coor[0] < 0 ? 0 : 1) : (internal_coor[0]
602 < 0 ? 3 : 2);
603 */
604
605 // compute the normal to the deformed surface and the
606 // scale of the surface
607 double n[3];
608 double scale;
609 double d_n_d_dof[SHAPE::num_nodes * 3][3];
610 double d_scale_d_dof[SHAPE::num_nodes * 3];
611 {
612 b2linalg::Matrix<double> d_shape(2, SHAPE::num_nodes);
613 SHAPE::eval_h_derived(internal_coor, d_shape);
614 double d_geom_d_icoor_tmp[2][3];
615 for (size_t i = 0; i != 3; ++i) {
616 for (size_t j = 0; j != 2; ++j) {
617 double tmp = 0;
618 for (int k = 0; k != SHAPE::num_nodes; ++k) {
619 tmp += d_shape(j, k) * node_coor[k][i];
620 }
621 d_geom_d_icoor_tmp[j][i] = tmp;
622 }
623 }
624 double inv_undef_area =
625 1 / norm_outer_product_3(d_geom_d_icoor_tmp[0], d_geom_d_icoor_tmp[1]);
626 if (dof.size2()) {
627 for (size_t i = 0; i != 3; ++i) {
628 for (size_t j = 0; j != 2; ++j) {
629 double tmp = 0;
630 for (int k = 0; k != SHAPE::num_nodes; ++k) {
631 tmp += d_shape(j, k) * dof(dn[k * (EULER ? 4 : 7) + i], 0);
632 }
633 d_geom_d_icoor_tmp[j][i] += tmp;
634 }
635 }
636 outer_product_3(d_geom_d_icoor_tmp[0], d_geom_d_icoor_tmp[1], n);
637 scale = std::sqrt(dot_3(n, n)) * inv_undef_area;
638 } else {
639 scale = 1;
640 }
641 for (int i = 0; i != 3; ++i) { n[i] *= inv_undef_area; }
642 if (!d_value_d_dof.is_null()) {
643 for (int k = 0; k != SHAPE::num_nodes; ++k) {
644 for (int j = 0; j != 3; ++j) {
645 double tmp_u[3] = {};
646 tmp_u[j] = d_shape(0, k);
647 double tmp_v[3] = {};
648 tmp_v[j] = d_shape(1, k);
649 double* res = d_n_d_dof + k * 3 + j;
650 outer_product_3(d_geom_d_icoor_tmp[0], tmp_v, res);
651 outer_product_add_3(tmp_u, d_geom_d_icoor_tmp[1], res);
652 for (int i = 0; i != 3; ++i) { res[i] *= inv_undef_area; }
653 }
654 }
655 double s = 2 / scale;
656 for (int k = 0; k != SHAPE::num_nodes * 3; ++k) {
657 double tmp = 0;
658 for (int i = 0; i != 3; ++i) { tmp += n[i] * d_n_d_dof[k][i]; }
659 d_scale_d_dof[k] = s * tmp;
660 }
661 }
662 }
663
664 double p = property->get_internal_pressure();
665 if (dof.size2()) { p -= dof(dn[(EULER ? 4 : 7) * pos_node], 0); }
666
667 if (!value.is_null()) {
668 value.resize(3);
669 if (dof.size2()) {
670 const size_t s_trac = pos_node * (EULER ? 4 : 7) + 4;
671 for (int i = 0; i != 3; ++i) {
672 value[i] = scale * dof(dn[s_trac + i], 0) + n[i] * p;
673 }
674 } else {
675 for (int i = 0; i != 3; ++i) { value[i] = n[i] * p; }
676 }
677 }
678
679 if (!d_value_d_icoor.is_null()) {
680 d_value_d_icoor.resize(3, 2);
681 d_value_d_icoor.set_zero();
682 }
683
684 if (!d_value_d_dof.is_null()) {
685 dof_numbering.resize(SHAPE::num_nodes * 3 + (EULER ? 1 : 4));
686 {
687 size_t i = 0;
688 for (int k = 0; k != SHAPE::num_nodes; ++k, i += 3) {
689 std::copy(
690 dn + k * (EULER ? 4 : 7), dn + k * (EULER ? 4 : 7) + 3,
691 &dof_numbering[i]);
692 }
693 std::copy(
694 dn + pos_node * (EULER ? 4 : 7) + 3,
695 dn + (pos_node + 1) * (EULER ? 4 : 7), &dof_numbering[i]);
696 }
697
698 const size_t s_trac = pos_node * (EULER ? 4 : 7) + 4;
699 d_value_d_dof.resize(3, dof_numbering.size());
700 for (int k = 0; k != SHAPE::num_nodes; ++k) {
701 for (int j = 0; j != 3; ++j) {
702 for (int i = 0; i != 3; ++i) {
703 d_value_d_dof(i, k * 3 + j) =
704 d_scale_d_dof[k * 3 + j] * dof(dn[s_trac + i], 0)
705 + d_n_d_dof[k * 3 + j][i] * p;
706 }
707 }
708 }
709 for (int i = 0; i != 3; ++i) {
710 d_value_d_dof(i, SHAPE::num_nodes * 3) = -n[i];
711 for (int j = 0; j != 3; ++j) {
712 d_value_d_dof(i, SHAPE::num_nodes * 3 + 1 + j) = scale;
713 }
714 }
715 }
716 }
717 }
718
719 void set_property(ElementProperty* property_) override {
720 property = dynamic_cast<CFDMaterial*>(property_);
721 if (property == nullptr) {
722 TypeError() << "Bad property type " << typeid(*property_)
723 << " for ElementCFD2DFieldL2 element." << THROW;
724 }
725 }
726
727 const ElementProperty* get_property() const override { return property; }
728
729 static type_t type;
730
731protected:
732 Node* nodes[SHAPE::num_nodes];
733
734 // The property
735 CFDMaterial* property;
736};
737
742template <bool FLUX, bool AXISYMMETRIC>
743class ElementCFDHeat2DFieldL2Constant : public TypedElement<double> {
744public:
745 using node_type =
746 node::HNode<coordof::Coor<coordof::Trans<coordof::X, coordof::Y>>, coordof::Dof<>>;
747
748 using type_t =
749 ObjectTypeComplete<ElementCFDHeat2DFieldL2Constant, TypedElement<double>::type_t>;
750
751 std::pair<int, Node* const*> get_nodes() const override {
752 return std::pair<int, Node* const*>(2, nodes);
753 }
754
755 void set_nodes(std::pair<int, Node* const*> nodes_) override {
756 if (nodes_.first != 2) { Exception() << "This element must have 2 nodes." << THROW; }
757 for (int i = 0; i != 2; ++i) { nodes[i] = nodes_.second[i]; }
758 }
759
760 int get_number_of_dof() const override { return FLUX ? 2 : 1; }
761
762 size_t set_global_dof_numbering(size_t index) override {
763 dof_index = index;
764 return index + (FLUX ? 2 : 1);
765 }
766
767 std::pair<size_t, size_t> get_global_dof_numbering() const override {
768 return std::pair<size_t, size_t>(dof_index, dof_index + (FLUX ? 2 : 1));
769 }
770
771 int edge_field_order(const int edge_id, const std::string& field_name) override {
772 return AXISYMMETRIC && FLUX ? 2 : 0;
773 }
774
775 bool edge_field_linear_on_dof(const int edge_id, const std::string& field_name) override {
776 return true;
777 }
778
779 int edge_field_polynomial_sub_edge(
780 const int edge_id, const std::string& field_name,
781 b2linalg::Vector<double, b2linalg::Vdense>& sub_nodes) override {
782 sub_nodes.clear();
783 sub_nodes.reserve(2);
784 sub_nodes.push_back(0);
785 sub_nodes.push_back(1);
786 return 0;
787 }
788
789 void edge_geom(
790 const int edge_id, const double internal_coor, b2linalg::Vector<double>& geom,
791 b2linalg::Vector<double>& d_geom_d_icoor) override {
792 double node_coor[2][2];
793 for (int k = 0; k != 2; ++k) { node_type::get_coor_s(node_coor[k], nodes[k]); }
794 if (!geom.is_null()) {
795 geom.resize(2);
796 for (size_t i = 0; i != 2; ++i) {
797 geom[i] = node_coor[0][i] * (1 - internal_coor) + node_coor[1][i] * internal_coor;
798 }
799 }
800 if (!d_geom_d_icoor.is_null()) {
801 d_geom_d_icoor.resize(2);
802 for (size_t i = 0; i != 2; ++i) {
803 d_geom_d_icoor[i] = -node_coor[0][i] + node_coor[1][i];
804 }
805 }
806 }
807
808 void edge_field_value(
809 const int edge_id, const std::string& field_name, const double internal_coor,
810 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof, const double time,
811 b2linalg::Vector<double, b2linalg::Vdense>& value,
812 b2linalg::Vector<double, b2linalg::Vdense>& d_value_d_icoor,
813 b2linalg::Index& dof_numbering,
814 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_value_d_dof,
815 b2linalg::Index& d_value_d_dof_dep) override {
816 if (!FLUX || field_name == "Temperature") {
817 if (!value.is_null()) {
818 value.resize(1);
819 if (dof.size2() > 0) {
820 value[0] = dof(dof_index, 0);
821 } else {
822 value[0] = 0;
823 }
824 }
825 if (!d_value_d_icoor.is_null()) {
826 d_value_d_icoor.resize(1);
827 d_value_d_icoor[0] = 0;
828 }
829 if (!d_value_d_dof.is_null()) {
830 dof_numbering.resize(1);
831 dof_numbering[0] = dof_index;
832 d_value_d_dof.resize(1, 1);
833 d_value_d_dof(0, 0) = 1;
834 }
835 } else if (FLUX && field_name == "HeatFlux") {
836 double tickness, lenght;
837 if (AXISYMMETRIC) {
838 double node_coor[2][2];
839 for (int k = 0; k != 2; ++k) { node_type::get_coor_s(node_coor[k], nodes[k]); }
840 tickness =
841 2 * M_PIl
842 * (internal_coor * node_coor[0][1] + (1 - internal_coor) * node_coor[1][1]);
843 lenght = node_coor[0][1] - node_coor[1][1];
844 } else {
845 tickness = 1;
846 }
847 if (!value.is_null()) {
848 value.resize(1);
849 if (dof.size2() > 0) {
850 value[0] = dof(dof_index + 1, 0) * tickness;
851 } else {
852 value[0] = 0;
853 }
854 }
855 if (!d_value_d_icoor.is_null()) {
856 d_value_d_icoor.resize(1);
857 if (AXISYMMETRIC && dof.size2() > 0) {
858 d_value_d_icoor[0] = dof(dof_index + 1, 0) * 2 * M_PIl * lenght;
859 } else {
860 d_value_d_icoor[0] = 0;
861 }
862 }
863 if (!d_value_d_dof.is_null()) {
864 dof_numbering.resize(1);
865 dof_numbering[0] = dof_index + 1;
866 d_value_d_dof.resize(1, 1);
867 d_value_d_dof(0, 0) = tickness;
868 }
869 }
870 }
871
872 static type_t type;
873
874private:
875 Node* nodes[2];
876 size_t dof_index;
877};
878
879} // namespace b2000
880
881#endif /* B2ELEMENT_FIELD_TRANSFER_H_ */
b2000::details::sqrt
csda< T > sqrt(const csda< T > &a)
Square root of a csda number.
Definition b2csda.H:396
THROW
#define THROW
Definition b2exception.H:198
b2tensor_calculus.H
b2000::ElementCFD3DField
Definition b2element_field_transfer.H:374
b2000::ElementCFD3DField::get_nodes
std::pair< int, Node *const * > get_nodes() const override
Definition b2element_field_transfer.H:388
b2000::ElementCFD3DField::get_property
const ElementProperty * get_property() const override
Definition b2element_field_transfer.H:727
b2000::ElementCFD3DField::set_nodes
void set_nodes(std::pair< int, Node *const * > nodes_) override
Definition b2element_field_transfer.H:392
b2000::ElementCFD3DField::face_geom
void face_geom(const int face_id, const b2linalg::Vector< double, b2linalg::Vdense_constref > &internal_coor, b2linalg::Vector< double > &geom, b2linalg::Matrix< double, b2linalg::Mrectangle > &d_geom_d_icoor) override
Definition b2element_field_transfer.H:476
b2000::ElementCFD3DField::face_field_order
int face_field_order(const int face_id, const std::string &field_name) override
Definition b2element_field_transfer.H:399
b2000::ElementCFD3DField::set_property
void set_property(ElementProperty *property_) override
Definition b2element_field_transfer.H:719
b2000::ElementCFD3DField::face_field_linear_on_dof
bool face_field_linear_on_dof(const int edge_id, const std::string &field_name) override
Definition b2element_field_transfer.H:403
b2000::ElementCFDHeat2DFieldL2Constant
Definition b2element_field_transfer.H:743
b2000::ElementCFDHeat2DFieldL2Constant::edge_field_order
int edge_field_order(const int edge_id, const std::string &field_name) override
Definition b2element_field_transfer.H:771
b2000::ElementCFDHeat2DFieldL2Constant::get_number_of_dof
int get_number_of_dof() const override
Definition b2element_field_transfer.H:760
b2000::ElementCFDHeat2DFieldL2Constant::edge_field_value
void edge_field_value(const int edge_id, const std::string &field_name, const double internal_coor, const b2linalg::Matrix< double, b2linalg::Mrectangle_constref > &dof, const double time, b2linalg::Vector< double, b2linalg::Vdense > &value, b2linalg::Vector< double, b2linalg::Vdense > &d_value_d_icoor, b2linalg::Index &dof_numbering, b2linalg::Matrix< double, b2linalg::Mrectangle > &d_value_d_dof, b2linalg::Index &d_value_d_dof_dep) override
Definition b2element_field_transfer.H:808
b2000::ElementCFDHeat2DFieldL2Constant::edge_field_linear_on_dof
bool edge_field_linear_on_dof(const int edge_id, const std::string &field_name) override
Definition b2element_field_transfer.H:775
b2000::ElementCFDHeat2DFieldL2Constant::set_nodes
void set_nodes(std::pair< int, Node *const * > nodes_) override
Definition b2element_field_transfer.H:755
b2000::ElementCFDHeat2DFieldL2Constant::set_global_dof_numbering
size_t set_global_dof_numbering(size_t index) override
Definition b2element_field_transfer.H:762
b2000::ElementCFDHeat2DFieldL2Constant::get_global_dof_numbering
std::pair< size_t, size_t > get_global_dof_numbering() const override
Definition b2element_field_transfer.H:767
b2000::ElementCFDHeat2DFieldL2Constant::get_nodes
std::pair< int, Node *const * > get_nodes() const override
Definition b2element_field_transfer.H:751
b2000::ElementCFDHeat2DFieldL2Constant::edge_geom
void edge_geom(const int edge_id, const double internal_coor, b2linalg::Vector< double > &geom, b2linalg::Vector< double > &d_geom_d_icoor) override
Definition b2element_field_transfer.H:789
b2000::ElementCFDHeat2DFieldL2Constant::edge_field_polynomial_sub_edge
int edge_field_polynomial_sub_edge(const int edge_id, const std::string &field_name, b2linalg::Vector< double, b2linalg::Vdense > &sub_nodes) override
Definition b2element_field_transfer.H:779
b2000::ElementProperty
Definition b2element.H:71
b2000::Exception
Definition b2exception.H:131
b2000::Node
Definition b2node.H:53
b2000::ObjectTypeComplete
Definition b2object.H:415
b2000::TypedElement
Definition b2element.H:1054
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32
b2000::norm_outer_product_3
T norm_outer_product_3(const T a[3], const T b[3])
Definition b2tensor_calculus.H:407
b2000::outer_product_3
void outer_product_3(const T1 a[3], const T2 b[3], T3 c[3])
Definition b2tensor_calculus.H:389
b2000::TypeError
GenericException< TypeError_name > TypeError
Definition b2exception.H:325
b2000::outer_product_add_3
void outer_product_add_3(const T1 a[3], const T2 b[3], T3 c[3])
Definition b2tensor_calculus.H:398
b2000::dot_3
T dot_3(const T a[3], const T b[3])
Definition b2tensor_calculus.H:328