b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2element_aeroacoustic.H
1//------------------------------------------------------------------------
2// b2element_aeroacoustic.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) 2006-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_AEROACOUSTIC_H_
22#define B2_ELEMENT_AEROACOUSTIC_H_
23
24#include <complex>
25
26#include "elements/properties/b2aeroacoustic.H"
27#include "model/b2element.H"
28#include "model_imp/b2hnode.H"
29#include "utils/b2exception.H"
30#include "utils/b2linear_algebra.H"
31#include "utils/b2tensor_calculus.H"
32
33#define EULER_2 1
34#define NB_POINT_IN_CONSTRAINT 4
35
36namespace b2000 {
37
38template <int NB_NODES, int NB_GAUSS, int NB_GAUSS_FACE>
39class GenericVolumeInterpolation {
40public:
41 static const int nb_nodes;
42 static const int nb_gauss;
43
44 // Node coordinate
45 double NodeCoor[NB_NODES][3];
46
47 // Integration point coordinate
48 double IntCoor[NB_GAUSS][3];
49
50 // Integration point weight
51 double weight[NB_GAUSS];
52
53 // The shape functions at the Gauss points.
54 double N[NB_GAUSS][NB_NODES];
55
56 // The derivative of shape function at the Gauss points.
57 double dN[NB_GAUSS][NB_NODES][3];
58
59 static const int nb_gauss_face;
60
61 // The integration point coordinate of the face 1
62 double IntCoor_face1[NB_GAUSS_FACE][3];
63
64 // Integration point weight
65 double weight_face1[NB_GAUSS_FACE];
66
67 // The shape functions at the Gauss points.
68 double N_face1[NB_GAUSS_FACE][NB_NODES];
69
70 // The derivative of shape function at the Gauss points.
71 double dN_face1[NB_GAUSS_FACE][NB_NODES][3];
72
73 // The centre of face 1
74 double IntCoor_face1_c[3];
75
76 // The shape function and the derivative at the centre of face 1
77 double N_face1_c[NB_NODES];
78 double dN_face1_c[NB_NODES][3];
79};
80
81template <int NB_NODES, int NB_GAUSS, int NB_GAUSS_FACE>
82const int GenericVolumeInterpolation<NB_NODES, NB_GAUSS, NB_GAUSS_FACE>::nb_nodes = NB_NODES;
83
84template <int NB_NODES, int NB_GAUSS, int NB_GAUSS_FACE>
85const int GenericVolumeInterpolation<NB_NODES, NB_GAUSS, NB_GAUSS_FACE>::nb_gauss = NB_GAUSS;
86
87template <int NB_NODES, int NB_GAUSS, int NB_GAUSS_FACE>
88const int GenericVolumeInterpolation<NB_NODES, NB_GAUSS, NB_GAUSS_FACE>::nb_gauss_face =
89 NB_GAUSS_FACE;
90
91template <typename INTERPOLATION>
92class ElementAeroAcoustic : public b2000::TypedElement<std::complex<double>> {
93public:
94 using node_type = node::HNode<
95 coordof::Coor<coordof::Trans<coordof::X, coordof::Y, coordof::Z>>,
96 coordof::Dof<
97 coordof::Density, coordof::DTrans<coordof::VX, coordof::VY, coordof::VZ>,
98 coordof::Temperature>>;
99
100 using type_t = ObjectTypeComplete<
101 ElementAeroAcoustic, typename TypedElement<std::complex<double>>::type_t>;
102
103 int get_number_of_dof() const override { return INTERPOLATION::nb_nodes; }
104
105 void set_nodes(std::pair<int, Node* const*> nodes_) override {
106 if (nodes_.first != INTERPOLATION::nb_nodes) {
107 Exception() << "This element has " << INTERPOLATION::nb_nodes << " nodes." << THROW;
108 }
109 for (int i = 0; i != INTERPOLATION::nb_nodes; ++i) { nodes[i] = nodes_.second[i]; }
110 }
111
112 std::pair<int, Node* const*> get_nodes() const override {
113 return std::pair<int, Node* const*>(INTERPOLATION::nb_nodes, nodes);
114 }
115
116 void set_property(ElementProperty* property_) override {
117 property = dynamic_cast<AeroAcousticPropertiesInterface*>(property_);
118 if (property == nullptr) {
119 TypeError() << "Bad or missing element property type for aeroacoustic element "
120 << get_object_name() << "(forgot MID or PID element attribute?) " << THROW;
121 }
122 }
123
124 const ElementProperty* get_property() const override { return property; }
125
126 void get_dof_numbering(b2linalg::Index& dof_numbering) override {
127 dof_numbering.resize(5 * INTERPOLATION::nb_nodes);
128 size_t* ptr = &dof_numbering[0];
129 for (int k = 0; k != INTERPOLATION::nb_nodes; ++k) {
130 ptr = nodes[k]->get_global_dof_numbering(ptr);
131 }
132 }
133
134 const std::vector<VariableInfo> get_value_info() const override {
135 return std::vector<VariableInfo>(1, Element::nonlinear);
136 }
137
138 void get_static_linear_value(
139 Model& model, b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& value,
140 b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle>& d_value_d_dof);
141
142 void get_value(
143 Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
144 const double time, const double delta_time,
145 const b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>& dof,
146 GradientContainer* gradient_container, SolverHints* solver_hints,
147 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& value,
148 std::vector<b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle>>& d_value_d_dof,
149 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& d_value_d_time) {
150 get_static_linear_value(
151 model, value,
152 d_value_d_dof.size() == 0 || d_value_d_dof[0].size1() == 0
153 ? b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle>::null
154 : d_value_d_dof[0]);
155 }
156
157 static type_t type;
158
159protected:
160 // The nodes list.
161 Node* nodes[INTERPOLATION::nb_nodes];
162
163 // The property
164 AeroAcousticPropertiesInterface* property;
165
166 static const INTERPOLATION N;
167
168 const bool symmetric{false};
169};
170
171template <typename INTERPOLATION>
172const INTERPOLATION ElementAeroAcoustic<INTERPOLATION>::N;
173
174template <typename INTERPOLATION>
175class ElementAeroAcousticRadiation : public ElementAeroAcoustic<INTERPOLATION> {
176public:
177 using parent = ElementAeroAcoustic<INTERPOLATION>;
178
179 using type_t = ObjectTypeComplete<
180 ElementAeroAcousticRadiation, typename TypedElement<std::complex<double>>::type_t>;
181
182 std::pair<int, Element::VariableInfo> get_constraint_info() {
183 return std::pair<int, Element::VariableInfo>(5 * NB_POINT_IN_CONSTRAINT, Element::linear);
184 }
185
186 void get_constraint(
187 Model& model, const bool linear, double time,
188 const b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>& dof,
189 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& constraint,
190 b2linalg::Matrix<std::complex<double>, b2linalg::Mcompressed_col>& d_constraint_d_dof,
191 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& d_constraint_d_time);
192
193 void get_static_linear_value(
194 Model& model, b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& value,
195 b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle>& d_value_d_dof);
196
197 static type_t type;
198};
199
200template <typename INTERPOLATION>
201class ElementAeroAcousticContacte : public ElementAeroAcoustic<INTERPOLATION> {
202public:
203 using parent = ElementAeroAcoustic<INTERPOLATION>;
204
205 using type_t = ObjectTypeComplete<
206 ElementAeroAcousticContacte, typename TypedElement<std::complex<double>>::type_t>;
207
208 std::pair<int, Element::VariableInfo> get_constraint_info() const {
209#ifdef EULER_2
210 return std::pair<int, Element::VariableInfo>(3 * NB_POINT_IN_CONSTRAINT, Element::linear);
211#else
212 return std::pair<int, Element::VariableInfo>(2 * NB_POINT_IN_CONSTRAINT, Element::linear);
213#endif
214 }
215
216 void get_constraint(
217 Model& model, const bool linear, double time,
218 const b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>& dof,
219 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& constraint,
220 b2linalg::Matrix<std::complex<double>, b2linalg::Mcompressed_col>& d_constraint_d_dof,
221 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& d_constraint_d_time);
222
223 void get_static_linear_value(
224 Model& model, b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& value,
225 b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle>& d_value_d_dof) {
226 const size_t s = INTERPOLATION::nb_nodes * 5;
227 value.resize(s);
228 value.set_zero();
229 d_value_d_dof.resize(s, s);
230 d_value_d_dof.set_zero();
231 }
232
233 static type_t type;
234};
235
236template <typename INTERPOLATION>
237class ElementAeroAcousticSymmetry : public ElementAeroAcoustic<INTERPOLATION> {
238public:
239 using parent = ElementAeroAcoustic<INTERPOLATION>;
240
241 using type_t = ObjectTypeComplete<
242 ElementAeroAcousticSymmetry, typename TypedElement<std::complex<double>>::type_t>;
243
244 std::pair<int, Element::VariableInfo> get_constraint_info() const {
245#ifdef EULER_2
246 return std::pair<int, Element::VariableInfo>(5 * NB_POINT_IN_CONSTRAINT, Element::linear);
247#else
248 return std::pair<int, Element::VariableInfo>(4 * NB_POINT_IN_CONSTRAINT, Element::linear);
249#endif
250 }
251
252 void get_constraint(
253 Model& model, const bool linear, double time,
254 const b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>& dof,
255 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& constraint,
256 b2linalg::Matrix<std::complex<double>, b2linalg::Mcompressed_col>& d_constraint_d_dof,
257 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& d_constraint_d_time);
258
259 void get_static_linear_value(
260 Model& model, b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& value,
261 b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle>& d_value_d_dof) {
262 const size_t s = INTERPOLATION::nb_nodes * 5;
263 value.resize(s);
264 value.set_zero();
265 d_value_d_dof.resize(s, s);
266 d_value_d_dof.set_zero();
267 }
268
269 static type_t type;
270};
271
272} // namespace b2000
273
274template <typename INTERPOLATION>
275void b2000::ElementAeroAcoustic<INTERPOLATION>::get_static_linear_value(
276 Model& model, b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& value,
277 b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle>& d_value_d_dof) {
278 const int number_of_dof = INTERPOLATION::nb_nodes * 5;
279 value.resize(number_of_dof);
280 value.set_zero();
281 d_value_d_dof.resize(number_of_dof, number_of_dof);
282 d_value_d_dof.set_zero();
283
284 // Get the nodes coordinate
285 double NodeCoor[INTERPOLATION::nb_nodes][3];
286 for (int k = 0; k != INTERPOLATION::nb_nodes; ++k) {
287 node_type::get_coor_s(NodeCoor[k], nodes[k]);
288 }
289
290 // Get the constante
291#ifdef EULER_2
292 double hcr = 1.4; // Heat capacity ratio
293 std::complex<double> omega;
294 {
295 double c_K;
296 double c_R;
297 double c_omega;
298 double c_r;
299 property->get_constant(c_K, c_R, c_r, c_omega);
300 omega = std::complex<double>(0, c_omega);
301 }
302#else
303 double c_K;
304 double c_cv;
305 double c_R;
306 std::complex<double> omega;
307 {
308 double c_omega;
309 double c_r;
310 property->get_constant(c_K, c_R, c_r, c_omega);
311 omega = std::complex<double>(0, c_omega);
312 c_cv = c_R / (c_K - 1);
313 }
314#endif
315
316 // Iterate through the integration point
317 for (int ng = 0; ng != INTERPOLATION::nb_gauss; ++ng) {
318 // Compute the Jacobian of the natural coordinate
319 // J[j][i] = d x[j] / d r[i]
320 double J[3][3];
321 for (int j = 0; j != 3; ++j) {
322 for (int i = 0; i != 3; ++i) {
323 double v = 0;
324 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
325 v += N.dN[ng][n][i] * NodeCoor[n][j];
326 }
327 J[j][i] = v;
328 }
329 }
330
331 // Compute the determinant and the inverse of J
332 double inv_J[3][3];
333 const double weight = invert_3_3(J, inv_J) * N.weight[ng];
334
335 // Compute the gradient of the interpolation
336 // B[j][i] = d N[j] / d x[i]
337 double B[INTERPOLATION::nb_nodes][3];
338 blas::gemm(
339 'N', 'N', 3, INTERPOLATION::nb_nodes, 3, 1.0, inv_J[0], 3, N.dN[ng][0], 3, 0.0, B[0],
340 3);
341
342 // Get the cfd value at the integration point
343 double W0[5];
344 double d_W0[5][3]; // d_w0[j][i] = d W0[j] / d x[i]
345 property->get_cfd_value_and_gradient(
346 model, INTERPOLATION::nb_nodes, nodes, N.N[ng], B, W0, d_W0);
347
348#ifdef EULER_2
349 {
350 const double Bg[5][5] = {
351 {d_W0[1][0] + d_W0[2][1] + d_W0[3][2], d_W0[0][0], d_W0[0][1], d_W0[0][2], 0},
352 {(W0[1] * d_W0[1][0] + W0[2] * d_W0[1][1] + W0[3] * d_W0[1][2]) / W0[0],
353 d_W0[1][0], d_W0[1][1], d_W0[1][2], 0},
354 {(W0[1] * d_W0[2][0] + W0[2] * d_W0[2][1] + W0[3] * d_W0[2][2]) / W0[0],
355 d_W0[2][0], d_W0[2][1], d_W0[2][2], 0},
356 {(W0[1] * d_W0[3][0] + W0[2] * d_W0[3][1] + W0[3] * d_W0[3][2]) / W0[0],
357 d_W0[3][0], d_W0[3][1], d_W0[3][2], 0},
358 {0, d_W0[4][0], d_W0[4][1], d_W0[4][2],
359 hcr * (d_W0[1][0] + d_W0[2][1] + d_W0[3][2])}};
360
361 const double A2 = hcr * W0[4];
362 const double inv_dense = 1 / W0[0];
363 for (int nj = 0; nj != INTERPOLATION::nb_nodes; ++nj) {
364 const int nnj = nj * 5;
365 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
366 const int nni = ni * 5;
367 {
368 const double v = weight * N.N[ng][ni];
369 const double vx = v * B[nj][0];
370 const double vy = v * B[nj][1];
371 const double vz = v * B[nj][2];
372 const double A3 = W0[1] * vx + W0[2] * vy + W0[3] * vz;
373 d_value_d_dof(nni, nnj) += A3;
374 d_value_d_dof(nni + 1, nnj + 1) += A3;
375 d_value_d_dof(nni + 2, nnj + 2) += A3;
376 d_value_d_dof(nni + 3, nnj + 3) += A3;
377 d_value_d_dof(nni + 4, nnj + 4) += A3;
378 d_value_d_dof(nni + 4, nnj + 1) += A2 * vx;
379 d_value_d_dof(nni + 4, nnj + 2) += A2 * vy;
380 d_value_d_dof(nni + 4, nnj + 3) += A2 * vz;
381 d_value_d_dof(nni, nnj + 1) += W0[0] * vx;
382 d_value_d_dof(nni, nnj + 2) += W0[0] * vy;
383 d_value_d_dof(nni, nnj + 3) += W0[0] * vz;
384 d_value_d_dof(nni + 1, nnj + 4) += inv_dense * vx;
385 d_value_d_dof(nni + 2, nnj + 4) += inv_dense * vy;
386 d_value_d_dof(nni + 3, nnj + 4) += inv_dense * vz;
387 }
388 {
389 const double v = weight * N.N[ng][ni] * N.N[ng][nj];
390 for (int j = 0; j != 5; ++j) {
391 d_value_d_dof(nni + j, nnj + j) += omega * v;
392 for (int i = 0; i != 5; ++i) {
393 d_value_d_dof(nni + i, nnj + j) += v * Bg[i][j];
394 }
395 }
396 }
397 }
398 }
399 }
400#else
401 {
402 // d_value_d_dof += - Bt_x * L * B_x - Bt_y * L * B_y - Bt_z * L * B_z
403 // d_value_d_dof += - Ht dL/dx B_x - Ht dL/dy B_y - Ht dL/dz B_z
404 // d_value_d_dof += Ht A_x B_x + Ht A_y B_y * Ht A_z B_z
405 // d_value_d_dof += Ht Bg H ; d_value_d_dofdof = Ht H
406 // d_value_d_dof += i * omega * Ht H
407 const double Bg[5][5] = {
408 {d_W0[1][0] + d_W0[2][1] + d_W0[3][2], d_W0[0][0], d_W0[0][1], d_W0[0][2], 0},
409 {(c_R * d_W0[4][0] + W0[1] * d_W0[1][0] + W0[2] * d_W0[1][1] + W0[3] * d_W0[1][2])
410 / W0[0],
411 d_W0[1][0], d_W0[1][1], d_W0[1][2], c_R / W0[0] * d_W0[0][0]},
412 {(c_R * d_W0[4][1] + W0[1] * d_W0[2][0] + W0[2] * d_W0[2][1] + W0[3] * d_W0[2][2])
413 / W0[0],
414 d_W0[2][0], d_W0[2][1], d_W0[2][2], c_R / W0[0] * d_W0[0][1]},
415 {(c_R * d_W0[4][2] + W0[1] * d_W0[3][0] + W0[2] * d_W0[3][1] + W0[3] * d_W0[3][2])
416 / W0[0],
417 d_W0[3][0], d_W0[3][1], d_W0[3][2], c_R / W0[0] * d_W0[0][2]},
418 {c_R * W0[4] / (c_cv * W0[0]) * (d_W0[1][0] + d_W0[2][1] + d_W0[3][2])
419 - (W0[1] * d_W0[4][0] + W0[2] * d_W0[4][1] + W0[3] * d_W0[4][2]) / W0[0],
420 d_W0[4][0], d_W0[4][1], d_W0[4][2],
421 c_R / c_cv * (d_W0[1][0] + d_W0[2][1] + d_W0[3][2])}};
422
423 const double L = -weight * c_K / (c_cv * W0[0]);
424 const double dL = -weight * c_K / (c_cv * W0[0] * W0[0]);
425 const double dL_dx = dL * d_W0[0][0];
426 const double dL_dy = dL * d_W0[0][1];
427 const double dL_dz = dL * d_W0[0][2];
428 const double A1 = c_R * W0[4] / W0[0];
429 const double A2 = c_R * c_cv / W0[4];
430 for (int nj = 0; nj != INTERPOLATION::nb_nodes; ++nj) {
431 const int nnj = nj * 5;
432 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
433 const int nni = ni * 5;
434 d_value_d_dof(nni + 4, nnj + 4) +=
435 N.N[ng][ni] * (dL_dx * B[nj][0] + dL_dy * B[nj][1] + dL_dz * B[nj][2])
436 + L * (B[ni][0] * B[nj][0] + B[ni][1] * B[nj][1] + B[ni][2] * B[nj][2]);
437 {
438 const double v = weight * N.N[ng][ni];
439 const double vx = v * B[nj][0];
440 const double vy = v * B[nj][1];
441 const double vz = v * B[nj][2];
442 const double A3 = W0[1] * vx + W0[2] * vy + W0[3] * vz;
443 d_value_d_dof(nni, nnj) += A3;
444 d_value_d_dof(nni + 1, nnj + 1) += A3;
445 d_value_d_dof(nni + 2, nnj + 2) += A3;
446 d_value_d_dof(nni + 3, nnj + 3) += A3;
447 d_value_d_dof(nni + 4, nnj + 4) += A3;
448 d_value_d_dof(nni + 1, nnj) += A1 * vx;
449 d_value_d_dof(nni + 2, nnj) += A1 * vy;
450 d_value_d_dof(nni + 3, nnj) += A1 * vz;
451 d_value_d_dof(nni + 4, nnj + 1) += A2 * vx;
452 d_value_d_dof(nni + 4, nnj + 2) += A2 * vy;
453 d_value_d_dof(nni + 4, nnj + 3) += A2 * vz;
454 d_value_d_dof(nni, nnj + 1) += W0[0] * vx;
455 d_value_d_dof(nni, nnj + 2) += W0[0] * vy;
456 d_value_d_dof(nni, nnj + 3) += W0[0] * vz;
457 d_value_d_dof(nni + 1, nnj + 4) += c_R * vx;
458 d_value_d_dof(nni + 2, nnj + 4) += c_R * vy;
459 d_value_d_dof(nni + 3, nnj + 4) += c_R * vz;
460 }
461 {
462 const double v = weight * N.N[ng][ni] * N.N[ng][nj];
463 for (int j = 0; j != 5; ++j) {
464 d_value_d_dof(nni + j, nnj + j) += omega * v;
465 for (int i = 0; i != 5; ++i) {
466 d_value_d_dof(nni + i, nnj + j) += v * Bg[i][j];
467 }
468 }
469 }
470 }
471 }
472 }
473#endif
474 }
475}
476
477template <typename INTERPOLATION>
478void b2000::ElementAeroAcousticRadiation<INTERPOLATION>::get_static_linear_value(
479 Model& model, b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& value,
480 b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle>& d_value_d_dof) {
481#ifdef EULER_2
482 const size_t s = INTERPOLATION::nb_nodes * 5;
483 value.resize(s);
484 value.set_zero();
485 d_value_d_dof.resize(s, s);
486 d_value_d_dof.set_zero();
487#else
488 const int number_of_dof = INTERPOLATION::nb_nodes * 5;
489 value.resize(number_of_dof);
490 value.set_zero();
491 d_value_d_dof.resize(number_of_dof, number_of_dof);
492 d_value_d_dof.set_zero();
493
494 // Get the nodes coordinate
495 double NodeCoor[INTERPOLATION::nb_nodes][3];
496 for (int k = 0; k != INTERPOLATION::nb_nodes; ++k) {
497 parent::node_type::get_coor_s(NodeCoor[k], parent::nodes[k]);
498 }
499
500 // Get the constant
501 double c_K;
502 double c_cv;
503 {
504 double c_omega;
505 double c_R;
506 double c_r;
507 parent::property->get_constant(c_K, c_R, c_r, c_omega);
508 c_cv = c_R / (c_K - 1);
509 }
510
511 // Iterate through the face 1 integration point
512 for (int ng = 0; ng != INTERPOLATION::nb_gauss_face; ++ng) {
513 // Compute the Jacobian of the natural coordinate
514 // J[j][i] = d x[j] / d r[i]
515 double J[3][3];
516 for (int j = 0; j != 3; ++j) {
517 for (int i = 0; i != 3; ++i) {
518 double v = 0;
519 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
520 v += parent::N.dN_face1[ng][n][i] * NodeCoor[n][j];
521 }
522 J[j][i] = v;
523 }
524 }
525
526 // Compute the determinant and the inverse of J
527 double inv_J[3][3];
528 const double weight = invert_3_3(J, inv_J) * parent::N.weight_face1[ng];
529
530 // Compute the gradient of the interpolation
531 // B[j][i] = d N[j] / d x[i]
532 double B[INTERPOLATION::nb_nodes][3];
533 blas::gemm(
534 'N', 'N', 3, INTERPOLATION::nb_nodes, 3, 1.0, inv_J[0], 3, parent::N.dN_face1[ng][0],
535 3, 0.0, B[0], 3);
536
537 // Compute the normal to the face 5
538 double Nface[3] = {
539 J[1][1] * J[2][0] - J[2][1] * J[1][0], J[2][1] * J[0][0] - J[0][1] * J[2][0],
540 J[0][1] * J[1][0] - J[1][1] * J[0][0]};
541 normalise_3(Nface);
542
543 // Compute B_n
544 double Bnorm[INTERPOLATION::nb_nodes];
545 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
546 Bnorm[n] = B[n][0] * Nface[0] + B[n][1] * Nface[1] + B[n][2] * Nface[2];
547 }
548
549 // Get the cfd value at the integration point
550 double W0[5];
551 parent::property->get_cfd_value(
552 model, INTERPOLATION::nb_nodes, parent::nodes, parent::N.N_face1[ng], W0);
553
554 const double L = weight * c_K / (c_cv * W0[0]);
555
556 for (int nj = 0; nj != INTERPOLATION::nb_nodes; ++nj) {
557 const int nnj = nj + 5;
558 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
559 d_value_d_dof(ni * 5 + 4, nnj + 4) += parent::N.N_face1[ng][ni] * L * Bnorm[nj];
560 }
561 }
562 }
563#endif
564}
565
566template <typename INTERPOLATION>
567void b2000::ElementAeroAcousticRadiation<INTERPOLATION>::get_constraint(
568 Model& model, const bool linear, double time,
569 const b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>& dof,
570 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& constraint,
571 b2linalg::Matrix<std::complex<double>, b2linalg::Mcompressed_col>& trans_d_constraint_d_dof,
572 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& d_constraint_d_time) {
573#if NB_POINT_IN_CONSTRAINT == 1
574 // Get the nodes coordinate and the centre face coordinate
575 double NodeCoor[INTERPOLATION::nb_nodes][3];
576 double face_coor[3] = {};
577 for (int k = 0; k != INTERPOLATION::nb_nodes; ++k) {
578 parent::node_type::get_coor_s(NodeCoor[k], parent::nodes[k]);
579 face_coor[0] += parent::N.N_face1_c[k] * NodeCoor[k][0];
580 face_coor[1] += parent::N.N_face1_c[k] * NodeCoor[k][1];
581 face_coor[2] += parent::N.N_face1_c[k] * NodeCoor[k][2];
582 }
583
584 // Get the constant
585 double c_omega;
586 double vg;
587 double rayon;
588 double r[3];
589 {
590 double c_K;
591 double c_R;
592 double c_r;
593 parent::property->get_constant(c_K, c_R, c_r, c_omega);
594
595 // Get the cfd value at the integration point
596 double W0[5];
597 parent::property->get_cfd_value(
598 model, INTERPOLATION::nb_nodes, parent::nodes, parent::N.N_face1_c, W0);
599 double mean_sound_speed = std::sqrt(c_K * c_r * W0[4]);
600 double noise_source_pos[3];
601 parent::property->get_radiation_value(noise_source_pos);
602
603 r[0] = face_coor[0] - noise_source_pos[0];
604 r[1] = face_coor[1] - noise_source_pos[1];
605 r[2] = face_coor[2] - noise_source_pos[2];
606 rayon = std::sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]);
607 const double inv_r = 1 / rayon;
608 r[0] *= inv_r;
609 r[1] *= inv_r;
610 r[2] *= inv_r;
611 const double u_er = r[0] * W0[1] + r[0] * W0[2] + r[0] * W0[3];
612 vg = u_er
613 + std::sqrt(
614 mean_sound_speed * mean_sound_speed - W0[1] * W0[1] - W0[2] * W0[2]
615 - W0[3] * W0[3] + u_er * u_er);
616 }
617
618 // Compute the Jacobian of the natural coordinate
619 // J[j][i] = d x[j] / d r[i]
620 double J[3][3];
621 for (int j = 0; j != 3; ++j) {
622 for (int i = 0; i != 3; ++i) {
623 double v = 0;
624 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
625 v += parent::N.dN_face1_c[n][i] * NodeCoor[n][j];
626 }
627 J[j][i] = v;
628 }
629 }
630
631 // Compute the determinant and the inverse of J
632 double inv_J[3][3];
633 invert_3_3(J, inv_J);
634
635 // Compute the gradient of the interpolation
636 // B[j][i] = d N[j] / d x[i]
637 double B[INTERPOLATION::nb_nodes][3];
638 blas::gemm(
639 'N', 'N', 3, INTERPOLATION::nb_nodes, 3, 1.0, inv_J[0], 3, parent::N.dN_face1_c[0], 3,
640 0.0, B[0], 3);
641
642 // Compute B_r
643 double Bnorm[INTERPOLATION::nb_nodes];
644 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
645 Bnorm[n] = B[n][0] * r[0] + B[n][1] * r[1] + B[n][2] * r[2];
646 }
647
648 constraint.resize(5);
649 constraint.set_zero();
650
651 // Compute R_r = i * omega * H + vg * B_s + vg / rayon * H_s
652 trans_d_constraint_d_dof.resize(INTERPOLATION::nb_nodes * 5, 5, INTERPOLATION::nb_nodes * 5);
653 size_t* colind;
654 size_t* rowind;
655 std::complex<double>* value;
656 trans_d_constraint_d_dof.get_values(colind, rowind, value);
657 colind[0] = 0;
658 for (int i = 0; i != 5; ++i) { colind[i + 1] = INTERPOLATION::nb_nodes * i; }
659 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
660 const int nni = ni * 5;
661 const std::complex<double> v(
662 vg * (Bnorm[ni] + parent::N.N_face1_c[ni] / rayon),
663 c_omega * parent::N.N_face1_c[ni]);
664 for (int i = 0; i != 5; ++i) {
665 const size_t ptr = colind[i + 1]++;
666 rowind[ptr] = nni + i;
667 value[ptr] = v;
668 }
669 }
670#else
671 double noise_source_pos[3];
672 parent::property->get_radiation_value(noise_source_pos);
673 double c_omega;
674 double c_K;
675 double c_r;
676 {
677 double c_R;
678 parent::property->get_constant(c_K, c_R, c_r, c_omega);
679 }
680
681 // Get the nodes coordinate
682 double NodeCoor[INTERPOLATION::nb_nodes][3];
683 for (int k = 0; k != INTERPOLATION::nb_nodes; ++k) {
684 parent::node_type::get_coor_s(NodeCoor[k], parent::nodes[k]);
685 }
686
687 constraint.resize(5 * NB_POINT_IN_CONSTRAINT);
688 constraint.set_zero();
689
690 // Compute R_r = i * omega * H + vg * B_s + vg / rayon * H_s
691 trans_d_constraint_d_dof.resize(
692 INTERPOLATION::nb_nodes * 5, 5 * NB_POINT_IN_CONSTRAINT,
693 INTERPOLATION::nb_nodes * 5 * NB_POINT_IN_CONSTRAINT);
694 size_t* colind;
695 size_t* rowind;
696 std::complex<double>* value;
697 trans_d_constraint_d_dof.get_values(colind, rowind, value);
698 colind[0] = 0;
699 for (int i = 0; i != 5 * NB_POINT_IN_CONSTRAINT; ++i) {
700 colind[i + 1] = INTERPOLATION::nb_nodes * i;
701 }
702
703 for (int ng = 0; ng != INTERPOLATION::nb_gauss_face; ++ng) {
704 // Get the centre face coordinate
705 double face_coor[3] = {};
706 for (int k = 0; k != INTERPOLATION::nb_nodes; ++k) {
707 face_coor[0] += parent::N.N_face1[ng][k] * NodeCoor[k][0];
708 face_coor[1] += parent::N.N_face1[ng][k] * NodeCoor[k][1];
709 face_coor[2] += parent::N.N_face1[ng][k] * NodeCoor[k][2];
710 }
711
712 // Get the constant
713 double vg;
714 double rayon;
715 double r[3];
716 {
717 // Get the cfd value at the integration point
718 double W0[5];
719 parent::property->get_cfd_value(
720 model, INTERPOLATION::nb_nodes, parent::nodes, parent::N.N_face1[ng], W0);
721 double mean_sound_speed = std::sqrt(c_K * c_r * W0[4]);
722
723 r[0] = face_coor[0] - noise_source_pos[0];
724 r[1] = face_coor[1] - noise_source_pos[1];
725 r[2] = face_coor[2] - noise_source_pos[2];
726 rayon = std::sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]);
727 const double inv_r = 1 / rayon;
728 r[0] *= inv_r;
729 r[1] *= inv_r;
730 r[2] *= inv_r;
731 const double u_er = r[0] * W0[1] + r[0] * W0[2] + r[0] * W0[3];
732 vg = u_er
733 + std::sqrt(
734 mean_sound_speed * mean_sound_speed - W0[1] * W0[1] - W0[2] * W0[2]
735 - W0[3] * W0[3] + u_er * u_er);
736 }
737
738 // Compute the Jacobian of the natural coordinate
739 // J[j][i] = d x[j] / d r[i]
740 double J[3][3];
741 for (int j = 0; j != 3; ++j) {
742 for (int i = 0; i != 3; ++i) {
743 double v = 0;
744 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
745 v += parent::N.dN_face1[ng][n][i] * NodeCoor[n][j];
746 }
747 J[j][i] = v;
748 }
749 }
750
751 // Compute the determinant and the inverse of J
752 double inv_J[3][3];
753 invert_3_3(J, inv_J);
754
755 // Compute the gradient of the interpolation
756 // B[j][i] = d N[j] / d x[i]
757 double B[INTERPOLATION::nb_nodes][3];
758 blas::gemm(
759 'N', 'N', 3, INTERPOLATION::nb_nodes, 3, 1.0, inv_J[0], 3, parent::N.dN_face1[ng][0],
760 3, 0.0, B[0], 3);
761
762 // Compute B_r
763 double Bnorm[INTERPOLATION::nb_nodes];
764 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
765 Bnorm[n] = B[n][0] * r[0] + B[n][1] * r[1] + B[n][2] * r[2];
766 }
767
768 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
769 const int nni = ni * 5;
770 const std::complex<double> v(
771 vg * (Bnorm[ni] + parent::N.N_face1_c[ni] / rayon),
772 c_omega * parent::N.N_face1[ng][ni]);
773 for (int i = 0; i != 5; ++i) {
774 const size_t ptr = colind[5 * ng + i + 1]++;
775 rowind[ptr] = nni + i;
776 value[ptr] = v;
777 }
778 }
779 }
780#endif
781}
782
783template <typename INTERPOLATION>
784void b2000::ElementAeroAcousticContacte<INTERPOLATION>::get_constraint(
785 Model& model, const bool linear, double time,
786 const b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>& dof,
787 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& constraint,
788 b2linalg::Matrix<std::complex<double>, b2linalg::Mcompressed_col>& trans_d_constraint_d_dof,
789 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& d_constraint_d_time) {
790#if NB_POINT_IN_CONSTRAINT == 1
791 // Get the nodes coordinate
792 double NodeCoor[INTERPOLATION::nb_nodes][3];
793 for (int k = 0; k != INTERPOLATION::nb_nodes; ++k) {
794 parent::node_type::get_coor_s(NodeCoor[k], parent::nodes[k]);
795 }
796
797 // Compute the Jacobian of the natural coordinate
798 // J[j][i] = d x[j] / d r[i]
799 double J[3][3];
800 for (int j = 0; j != 3; ++j) {
801 for (int i = 0; i != 3; ++i) {
802 double v = 0;
803 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
804 v += parent::N.dN_face1_c[n][i] * NodeCoor[n][j];
805 }
806 J[j][i] = v;
807 }
808 }
809
810 // Compute the normal to the face 5
811 double Nface[3] = {
812 J[1][1] * J[2][0] - J[2][1] * J[1][0], J[2][1] * J[0][0] - J[0][1] * J[2][0],
813 J[0][1] * J[1][0] - J[1][1] * J[0][0]};
814 normalise_3(Nface);
815
816 // Compute the determinant and the inverse of J
817 double inv_J[3][3];
818 invert_3_3(J, inv_J);
819
820 // Compute the gradient of the interpolation
821 // B[j][i] = d N[j] / d x[i]
822 double B[INTERPOLATION::nb_nodes][3];
823 blas::gemm(
824 'N', 'N', 3, INTERPOLATION::nb_nodes, 3, 1.0, inv_J[0], 3, parent::N.dN_face1_c[0], 3,
825 0.0, B[0], 3);
826
827 // Compute B_n
828 double Bnorm[INTERPOLATION::nb_nodes];
829 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
830 Bnorm[n] = B[n][0] * Nface[0] + B[n][1] * Nface[1] + B[n][2] * Nface[2];
831 }
832
833#ifdef EULER_2
834 constraint.resize(3);
835 constraint.set_zero();
836 trans_d_constraint_d_dof.resize(INTERPOLATION::nb_nodes * 5, 3, INTERPOLATION::nb_nodes * 5);
837 size_t* colind;
838 size_t* rowind;
839 std::complex<double>* value;
840 trans_d_constraint_d_dof.get_values(colind, rowind, value);
841 colind[0] = 0;
842 size_t ptr = 0;
843 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
844 int nni = ni * 5;
845 rowind[ptr] = nni;
846 value[ptr++] = Bnorm[ni];
847 }
848 colind[1] = ptr;
849 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
850 int nni = ni * 5;
851 rowind[ptr] = ++nni;
852 value[ptr++] = Nface[0] * parent::N.N_face1_c[ni];
853
854 rowind[ptr] = ++nni;
855 value[ptr++] = Nface[1] * parent::N.N_face1_c[ni];
856
857 rowind[ptr] = ++nni;
858 value[ptr++] = Nface[2] * parent::N.N_face1_c[ni];
859 }
860 colind[2] = ptr;
861 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
862 int nni = ni * 5;
863 rowind[ptr] = nni + 4;
864 value[ptr++] = Bnorm[ni];
865 }
866 colind[3] = ptr;
867#else
868 // Compute R_c = (0, n_x, n_y, n_z) * H_s
869 constraint.resize(2);
870 constraint.set_zero();
871
872 trans_d_constraint_d_dof.resize(INTERPOLATION::nb_nodes * 5, 2, INTERPOLATION::nb_nodes * 4);
873 size_t* colind;
874 size_t* rowind;
875 std::complex<double>* value;
876 trans_d_constraint_d_dof.get_values(colind, rowind, value);
877 colind[0] = 0;
878 size_t ptr = 0;
879 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
880 int nni = ni * 5;
881 rowind[ptr] = nni;
882 value[ptr++] = Bnorm[ni];
883 }
884 colind[1] = ptr;
885 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
886 int nni = ni * 5;
887 rowind[ptr] = ++nni;
888 value[ptr++] = Nface[0] * parent::N.N_face1_c[ni];
889
890 rowind[ptr] = ++nni;
891 value[ptr++] = Nface[1] * parent::N.N_face1_c[ni];
892
893 rowind[ptr] = ++nni;
894 value[ptr++] = Nface[2] * parent::N.N_face1_c[ni];
895 }
896 colind[2] = ptr;
897#endif
898
899#else
900 // Get the nodes coordinate
901 double NodeCoor[INTERPOLATION::nb_nodes][3];
902 for (int k = 0; k != INTERPOLATION::nb_nodes; ++k) {
903 parent::node_type::get_coor_s(NodeCoor[k], parent::nodes[k]);
904 }
905
906 constraint.resize(3 * NB_POINT_IN_CONSTRAINT);
907 constraint.set_zero();
908
909 trans_d_constraint_d_dof.resize(
910 INTERPOLATION::nb_nodes * 5, 3 * NB_POINT_IN_CONSTRAINT,
911 INTERPOLATION::nb_nodes * 5 * NB_POINT_IN_CONSTRAINT);
912 size_t* colind;
913 size_t* rowind;
914 std::complex<double>* value;
915 trans_d_constraint_d_dof.get_values(colind, rowind, value);
916 colind[0] = 0;
917 size_t ptr = 0;
918
919 for (int ng = 0; ng != INTERPOLATION::nb_gauss_face; ++ng) {
920 // Compute the Jacobian of the natural coordinate
921 // J[j][i] = d x[j] / d r[i]
922 double J[3][3];
923 for (int j = 0; j != 3; ++j) {
924 for (int i = 0; i != 3; ++i) {
925 double v = 0;
926 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
927 v += parent::N.dN_face1[ng][n][i] * NodeCoor[n][j];
928 }
929 J[j][i] = v;
930 }
931 }
932
933 // Compute the normal to the face 5
934 double Nface[3] = {
935 J[1][1] * J[2][0] - J[2][1] * J[1][0], J[2][1] * J[0][0] - J[0][1] * J[2][0],
936 J[0][1] * J[1][0] - J[1][1] * J[0][0]};
937 normalise_3(Nface);
938
939 // Compute the determinant and the inverse of J
940 double inv_J[3][3];
941 invert_3_3(J, inv_J);
942
943 // Compute the gradient of the interpolation
944 // B[j][i] = d N[j] / d x[i]
945 double B[INTERPOLATION::nb_nodes][3];
946 blas::gemm(
947 'N', 'N', 3, INTERPOLATION::nb_nodes, 3, 1.0, inv_J[0], 3, parent::N.dN_face1[ng][0],
948 3, 0.0, B[0], 3);
949
950 // Compute B_n
951 double Bnorm[INTERPOLATION::nb_nodes];
952 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
953 Bnorm[n] = B[n][0] * Nface[0] + B[n][1] * Nface[1] + B[n][2] * Nface[2];
954 }
955
956 // Compute R_c = (0, n_x, n_y, n_z) * H_s
957 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
958 int nni = ni * 5;
959 rowind[ptr] = nni;
960 value[ptr++] = Bnorm[ni];
961 }
962 *(++colind) = ptr;
963 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
964 int nni = ni * 5;
965 rowind[ptr] = ++nni;
966 value[ptr++] = Nface[0] * parent::N.N_face1[ng][ni];
967
968 rowind[ptr] = ++nni;
969 value[ptr++] = Nface[1] * parent::N.N_face1[ng][ni];
970
971 rowind[ptr] = ++nni;
972 value[ptr++] = Nface[2] * parent::N.N_face1[ng][ni];
973 }
974 *(++colind) = ptr;
975 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni) {
976 int nni = ni * 5;
977 rowind[ptr] = nni + 4;
978 value[ptr++] = Bnorm[ni];
979 }
980 *(++colind) = ptr;
981 }
982#endif
983}
984
985template <typename INTERPOLATION>
986void b2000::ElementAeroAcousticSymmetry<INTERPOLATION>::get_constraint(
987 Model& model, const bool linear, double time,
988 const b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>& dof,
989 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& constraint,
990 b2linalg::Matrix<std::complex<double>, b2linalg::Mcompressed_col>& trans_d_constraint_d_dof,
991 b2linalg::Vector<std::complex<double>, b2linalg::Vdense>& d_constraint_d_time) {
992#if NB_POINT_IN_CONSTRAINT == 1
993 // Get the nodes coordinate
994 double NodeCoor[INTERPOLATION::nb_nodes][3];
995 for (int k = 0; k != INTERPOLATION::nb_nodes; ++k) {
996 parent::node_type::get_coor_s(NodeCoor[k], parent::nodes[k]);
997 }
998
999 // Compute the Jacobian of the natural coordinate
1000 // J[j][i] = d x[j] / d r[i]
1001 double J[3][3];
1002 for (int j = 0; j != 3; ++j) {
1003 for (int i = 0; i != 3; ++i) {
1004 double v = 0;
1005 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
1006 v += parent::N.dN_face1_c[n][i] * NodeCoor[n][j];
1007 }
1008 J[j][i] = v;
1009 }
1010 }
1011
1012 // Compute the normal to the face 5
1013 double Nface[3] = {
1014 J[1][1] * J[2][0] - J[2][1] * J[1][0], J[2][1] * J[0][0] - J[0][1] * J[2][0],
1015 J[0][1] * J[1][0] - J[1][1] * J[0][0]};
1016 normalise_3(Nface);
1017
1018 // Compute the determinant and the inverse of J
1019 double inv_J[3][3];
1020 invert_3_3(J, inv_J);
1021
1022 // Compute the gradient of the interpolation
1023 // B[j][i] = d N[j] / d x[i]
1024 double B[INTERPOLATION::nb_nodes][3];
1025 blas::gemm(
1026 'N', 'N', 3, INTERPOLATION::nb_nodes, 3, 1.0, inv_J[0], 3, parent::N.dN_face1_c[0], 3,
1027 0.0, B[0], 3);
1028
1029 // Compute B_n
1030 double Bnorm[INTERPOLATION::nb_nodes];
1031 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
1032 Bnorm[n] = B[n][0] * Nface[0] + B[n][1] * Nface[1] + B[n][2] * Nface[2];
1033 }
1034
1035#ifdef EULER_2
1036 constraint.resize(5);
1037 constraint.set_zero();
1038
1039 // Compute R_c = (0, n_x, n_y, n_z) * H_s
1040 trans_d_constraint_d_dof.resize(INTERPOLATION::nb_nodes * 5, 5, INTERPOLATION::nb_nodes * 5);
1041 size_t* colind;
1042 size_t* rowind;
1043 std::complex<double>* value;
1044 trans_d_constraint_d_dof.get_values(colind, rowind, value);
1045 size_t ptr = 0;
1046 for (int i = 0; i != 5; ++i) {
1047 colind[i] = ptr;
1048 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni, ++ptr) {
1049 rowind[ptr] = ni * 5 + i;
1050 value[ptr] = Bnorm[ni];
1051 }
1052 }
1053 colind[5] = ptr;
1054#else
1055 constraint.resize(4);
1056 constraint.set_zero();
1057
1058 // Compute R_c = (0, n_x, n_y, n_z) * H_s
1059 trans_d_constraint_d_dof.resize(INTERPOLATION::nb_nodes * 5, 4, INTERPOLATION::nb_nodes * 4);
1060 size_t* colind;
1061 size_t* rowind;
1062 std::complex<double>* value;
1063 trans_d_constraint_d_dof.get_values(colind, rowind, value);
1064 size_t ptr = 0;
1065 for (int i = 0; i != 4; ++i) {
1066 colind[i] = ptr;
1067 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni, ++ptr) {
1068 rowind[ptr] = ni * 5 + i;
1069 value[ptr] = Bnorm[ni];
1070 }
1071 }
1072 colind[4] = ptr;
1073#endif
1074
1075#else
1076 // Get the nodes coordinate
1077 double NodeCoor[INTERPOLATION::nb_nodes][3];
1078 for (int k = 0; k != INTERPOLATION::nb_nodes; ++k) {
1079 parent::node_type::get_coor_s(NodeCoor[k], parent::nodes[k]);
1080 }
1081
1082 constraint.resize(5 * NB_POINT_IN_CONSTRAINT);
1083 constraint.set_zero();
1084
1085 // Compute R_c = (0, n_x, n_y, n_z) * H_s
1086 trans_d_constraint_d_dof.resize(
1087 INTERPOLATION::nb_nodes * 5, 5 * NB_POINT_IN_CONSTRAINT,
1088 INTERPOLATION::nb_nodes * 5 * NB_POINT_IN_CONSTRAINT);
1089 size_t* colind;
1090 size_t* rowind;
1091 std::complex<double>* value;
1092 trans_d_constraint_d_dof.get_values(colind, rowind, value);
1093 size_t ptr = 0;
1094 colind[0] = 0;
1095
1096 // Iterate through the face 1 integration point
1097 for (int ng = 0; ng != INTERPOLATION::nb_gauss_face; ++ng) {
1098 // Compute the Jacobian of the natural coordinate
1099 // J[j][i] = d x[j] / d r[i]
1100 double J[3][3];
1101 for (int j = 0; j != 3; ++j) {
1102 for (int i = 0; i != 3; ++i) {
1103 double v = 0;
1104 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
1105 v += parent::N.dN_face1[ng][n][i] * NodeCoor[n][j];
1106 }
1107 J[j][i] = v;
1108 }
1109 }
1110
1111 // Compute the normal to the face 5
1112 double Nface[3] = {
1113 J[1][1] * J[2][0] - J[2][1] * J[1][0], J[2][1] * J[0][0] - J[0][1] * J[2][0],
1114 J[0][1] * J[1][0] - J[1][1] * J[0][0]};
1115 normalise_3(Nface);
1116
1117 // Compute the determinant and the inverse of J
1118 double inv_J[3][3];
1119 invert_3_3(J, inv_J);
1120
1121 // Compute the gradient of the interpolation
1122 // B[j][i] = d N[j] / d x[i]
1123 double B[INTERPOLATION::nb_nodes][3];
1124 blas::gemm(
1125 'N', 'N', 3, INTERPOLATION::nb_nodes, 3, 1.0, inv_J[0], 3, parent::N.dN_face1[ng][0],
1126 3, 0.0, B[0], 3);
1127
1128 // Compute B_n
1129 double Bnorm[INTERPOLATION::nb_nodes];
1130 for (int n = 0; n != INTERPOLATION::nb_nodes; ++n) {
1131 Bnorm[n] = B[n][0] * Nface[0] + B[n][1] * Nface[1] + B[n][2] * Nface[2];
1132 }
1133
1134 for (int i = 0; i != 5; ++i) {
1135 for (int ni = 0; ni != INTERPOLATION::nb_nodes; ++ni, ++ptr) {
1136 rowind[ptr] = ni * 5 + i;
1137 value[ptr] = Bnorm[ni];
1138 }
1139 *(++colind) = ptr;
1140 }
1141 }
1142#endif
1143}
1144
1145#endif // B2_ELEMENT_AEROACOUSTIC_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::Element::get_constraint_info
virtual std::pair< int, VariableInfo > get_constraint_info()
Definition b2element.H:688
b2000::TypedElement
Definition b2element.H:1054
b2000::TypedElement< std::complex< 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< std::complex< double >, b2linalg::Mrectangle_constref > &dof, GradientContainer *gradient_container, SolverHints *solver_hints, b2linalg::Index &dof_numbering, b2linalg::Vector< std::complex< double >, b2linalg::Vdense > &value, const std::vector< bool > &d_value_d_dof_flags, std::vector< b2linalg::Matrix< std::complex< double >, b2linalg::Mrectangle > > &d_value_d_dof, b2linalg::Vector< std::complex< double >, b2linalg::Vdense > &d_value_d_time)
Compute the internal forces and their derivatives of the element.
Definition b2element.H:1619
b2000::TypedElement< std::complex< double > >::get_constraint
virtual void get_constraint(Model &model, const bool linear, const double time, const b2linalg::Matrix< std::complex< double >, b2linalg::Mrectangle_constref > &dof, b2linalg::Index &dof_numbering, b2linalg::Vector< std::complex< double >, b2linalg::Vdense > &constraint, b2linalg::Matrix< std::complex< double >, b2linalg::Mcompressed_col > &trans_d_constraint_d_dof, b2linalg::Vector< std::complex< double >, b2linalg::Vdense > &d_constraint_d_time)
Definition b2element.H:1096
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32
b2000::normalise_3
T normalise_3(T a[3])
Definition b2tensor_calculus.H:418
b2000::TypeError
GenericException< TypeError_name > TypeError
Definition b2exception.H:325
b2000::invert_3_3
T invert_3_3(const T a[3][3], T b[3][3])
Definition b2tensor_calculus.H:699