macro_element.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2025 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#include "macro_element.h"
27#include "domain.h"
28
29namespace oomph
30{
31 //=================================================================
32 /// Get global position r(S) at discrete time level t.
33 /// t=0: Present time; t>0: previous timestep.
34 //=================================================================
35 void QMacroElement<2>::macro_map(const double& t,
36 const Vector<double>& s,
37 Vector<double>& r)
38 {
39 using namespace QuadTreeNames;
40
41 Vector<double> bound_N(2);
42 Vector<double> bound_S(2);
43 Vector<double> bound_W(2);
44 Vector<double> bound_E(2);
45
46 Vector<double> diff_N(2);
47 Vector<double> diff_S(2);
48 Vector<double> diff_W(2);
49 Vector<double> diff_E(2);
50
51 Vector<double> f_rect(2);
52
53 Vector<double> corner_SE(2);
54 Vector<double> corner_SW(2);
55 Vector<double> corner_NE(2);
56 Vector<double> corner_NW(2);
57
58 Vector<double> edge_N(2);
59 Vector<double> edge_S(2);
60
61 Vector<double> zeta(1);
62
63 // Determine Vectors to corners
64 zeta[0] = 1.0;
65 Domain_pt->macro_element_boundary(
66 t, Macro_element_number, QuadTreeNames::S, zeta, corner_SE);
67 zeta[0] = -1.0;
68 Domain_pt->macro_element_boundary(
69 t, Macro_element_number, QuadTreeNames::S, zeta, corner_SW);
70 zeta[0] = 1.0;
71 Domain_pt->macro_element_boundary(
72 t, Macro_element_number, QuadTreeNames::N, zeta, corner_NE);
73 zeta[0] = -1.0;
74 Domain_pt->macro_element_boundary(
75 t, Macro_element_number, QuadTreeNames::N, zeta, corner_NW);
76
77
78 // Get the position on the N/S/W/E edges
79 zeta[0] = s[0];
80 Domain_pt->macro_element_boundary(
81 t, Macro_element_number, QuadTreeNames::N, zeta, bound_N);
82 zeta[0] = s[0];
83 Domain_pt->macro_element_boundary(
84 t, Macro_element_number, QuadTreeNames::S, zeta, bound_S);
85 zeta[0] = s[1];
86 Domain_pt->macro_element_boundary(
87 t, Macro_element_number, QuadTreeNames::W, zeta, bound_W);
88 zeta[0] = s[1];
89 Domain_pt->macro_element_boundary(
90 t, Macro_element_number, QuadTreeNames::E, zeta, bound_E);
91
92
93 for (int i = 0; i < 2; i++)
94 {
95 // Position on the straight S/W edges of the rectangle formed
96 // by the corner points
97 edge_S[i] =
98 corner_SW[i] + (corner_SE[i] - corner_SW[i]) * 0.5 * (s[0] + 1.0);
99 edge_N[i] =
100 corner_NW[i] + (corner_NE[i] - corner_NW[i]) * 0.5 * (s[0] + 1.0);
101
102 // Position inside rectangle
103 f_rect[i] = edge_S[i] + (edge_N[i] - edge_S[i]) * 0.5 * (s[1] + 1.0);
104
105 // Get difference between curved edge and point in rectangle
106 diff_N[i] = bound_N[i] - f_rect[i];
107 diff_S[i] = bound_S[i] - f_rect[i];
108 diff_E[i] = bound_E[i] - f_rect[i];
109 diff_W[i] = bound_W[i] - f_rect[i];
110
111 // Map it...
112 r[i] = f_rect[i] + diff_S[i] * (1.0 - 0.5 * (s[1] + 1.0)) +
113 diff_N[i] * 0.5 * (s[1] + 1.0) +
114 diff_W[i] * (1.0 - 0.5 * (s[0] + 1.0)) +
115 diff_E[i] * 0.5 * (s[0] + 1.0);
116 }
117 }
118
119 //=================================================================
120 /// Get global position r(S) at discrete time level t.
121 /// t=0: Present time; t>0: previous timestep.
122 //=================================================================
123 void QMacroElement<2>::macro_map(const unsigned& t,
124 const Vector<double>& S,
125 Vector<double>& r)
126 {
127 using namespace QuadTreeNames;
128
129 Vector<double> bound_N(2);
130 Vector<double> bound_S(2);
131 Vector<double> bound_W(2);
132 Vector<double> bound_E(2);
133
134 Vector<double> diff_N(2);
135 Vector<double> diff_S(2);
136 Vector<double> diff_W(2);
137 Vector<double> diff_E(2);
138
139 Vector<double> f_rect(2);
140
141 Vector<double> corner_SE(2);
142 Vector<double> corner_SW(2);
143 Vector<double> corner_NE(2);
144 Vector<double> corner_NW(2);
145
146 Vector<double> edge_N(2);
147 Vector<double> edge_S(2);
148
149 Vector<double> zeta(1);
150
151 // Determine Vectors to corners
152 zeta[0] = 1.0;
153 Domain_pt->macro_element_boundary(
154 t, Macro_element_number, QuadTreeNames::S, zeta, corner_SE);
155 zeta[0] = -1.0;
156 Domain_pt->macro_element_boundary(
157 t, Macro_element_number, QuadTreeNames::S, zeta, corner_SW);
158 zeta[0] = 1.0;
159 Domain_pt->macro_element_boundary(
160 t, Macro_element_number, QuadTreeNames::N, zeta, corner_NE);
161 zeta[0] = -1.0;
162 Domain_pt->macro_element_boundary(
163 t, Macro_element_number, QuadTreeNames::N, zeta, corner_NW);
164
165
166 // Get the position on the N/S/W/E edges
167 zeta[0] = S[0];
168 Domain_pt->macro_element_boundary(
169 t, Macro_element_number, QuadTreeNames::N, zeta, bound_N);
170 zeta[0] = S[0];
171 Domain_pt->macro_element_boundary(
172 t, Macro_element_number, QuadTreeNames::S, zeta, bound_S);
173 zeta[0] = S[1];
174 Domain_pt->macro_element_boundary(
175 t, Macro_element_number, QuadTreeNames::W, zeta, bound_W);
176 zeta[0] = S[1];
177 Domain_pt->macro_element_boundary(
178 t, Macro_element_number, QuadTreeNames::E, zeta, bound_E);
179
180
181 for (int i = 0; i < 2; i++)
182 {
183 // Position on the straight S/W edges of the rectangle formed
184 // by the corner points
185 edge_S[i] =
186 corner_SW[i] + (corner_SE[i] - corner_SW[i]) * 0.5 * (S[0] + 1.0);
187 edge_N[i] =
188 corner_NW[i] + (corner_NE[i] - corner_NW[i]) * 0.5 * (S[0] + 1.0);
189
190 // Position inside rectangle
191 f_rect[i] = edge_S[i] + (edge_N[i] - edge_S[i]) * 0.5 * (S[1] + 1.0);
192
193 // Get difference between curved edge and point in rectangle
194 diff_N[i] = bound_N[i] - f_rect[i];
195 diff_S[i] = bound_S[i] - f_rect[i];
196 diff_E[i] = bound_E[i] - f_rect[i];
197 diff_W[i] = bound_W[i] - f_rect[i];
198
199 // Map it...
200 r[i] = f_rect[i] + diff_S[i] * (1.0 - 0.5 * (S[1] + 1.0)) +
201 diff_N[i] * 0.5 * (S[1] + 1.0) +
202 diff_W[i] * (1.0 - 0.5 * (S[0] + 1.0)) +
203 diff_E[i] * 0.5 * (S[0] + 1.0);
204 }
205 }
206
207
208 //=================================================================
209 /// Output all macro element boundaries as tecplot zones
210 //=================================================================
211 void QMacroElement<2>::output_macro_element_boundaries(std::ostream& outfile,
212 const unsigned& nplot)
213 {
214 using namespace QuadTreeNames;
215
216 Vector<double> s(1);
217 Vector<double> f(2);
218 // Dump at present time
219 unsigned t = 0;
220 for (unsigned idirect = N; idirect <= W; idirect++)
221 {
222 outfile << "ZONE I=" << nplot << std::endl;
223 for (unsigned j = 0; j < nplot; j++)
224 {
225 s[0] = -1.0 + 2.0 * double(j) / double(nplot - 1);
226 Domain_pt->macro_element_boundary(
227 t, Macro_element_number, idirect, s, f);
228 outfile << f[0] << " " << f[1] << std::endl;
229 }
230 }
231 }
232
233
234 //=============================================================================
235 /// Assembles the jacobian of the mapping from the macro coordinates to
236 /// the global coordinates
237 //=============================================================================
238 void QMacroElement<2>::assemble_macro_to_eulerian_jacobian(
239 const unsigned& t, const Vector<double>& S, DenseMatrix<double>& jacobian)
240 {
241 using namespace QuadTreeNames;
242
243 Vector<double> bound_N(2);
244 Vector<double> bound_S(2);
245 Vector<double> bound_W(2);
246 Vector<double> bound_E(2);
247
248 Vector<double> dbound_N(2);
249 Vector<double> dbound_S(2);
250 Vector<double> dbound_E(2);
251 Vector<double> dbound_W(2);
252
253 Vector<double> corner_SE(2);
254 Vector<double> corner_SW(2);
255 Vector<double> corner_NE(2);
256 Vector<double> corner_NW(2);
257
258 Vector<double> zeta(1);
259
260
261 // Determine Vectors to corners
262 zeta[0] = 1.0;
263 Domain_pt->macro_element_boundary(
264 t, Macro_element_number, QuadTreeNames::S, zeta, corner_SE);
265 zeta[0] = -1.0;
266 Domain_pt->macro_element_boundary(
267 t, Macro_element_number, QuadTreeNames::S, zeta, corner_SW);
268 zeta[0] = 1.0;
269 Domain_pt->macro_element_boundary(
270 t, Macro_element_number, QuadTreeNames::N, zeta, corner_NE);
271 zeta[0] = -1.0;
272 Domain_pt->macro_element_boundary(
273 t, Macro_element_number, QuadTreeNames::N, zeta, corner_NW);
274
275
276 // Get the position and first derivativeson the N/S/W/E edges
277 zeta[0] = S[0];
278 Domain_pt->macro_element_boundary(
279 t, Macro_element_number, QuadTreeNames::N, zeta, bound_N);
280 Domain_pt->dmacro_element_boundary(
281 t, Macro_element_number, QuadTreeNames::N, zeta, dbound_N);
282 zeta[0] = S[0];
283 Domain_pt->macro_element_boundary(
284 t, Macro_element_number, QuadTreeNames::S, zeta, bound_S);
285 Domain_pt->dmacro_element_boundary(
286 t, Macro_element_number, QuadTreeNames::S, zeta, dbound_S);
287 zeta[0] = S[1];
288 Domain_pt->macro_element_boundary(
289 t, Macro_element_number, QuadTreeNames::W, zeta, bound_W);
290 Domain_pt->dmacro_element_boundary(
291 t, Macro_element_number, QuadTreeNames::W, zeta, dbound_W);
292 zeta[0] = S[1];
293 Domain_pt->macro_element_boundary(
294 t, Macro_element_number, QuadTreeNames::E, zeta, bound_E);
295 Domain_pt->dmacro_element_boundary(
296 t, Macro_element_number, QuadTreeNames::E, zeta, dbound_E);
297
298
299 // dr0/dm0
300 jacobian(0, 0) =
301 0.25 * (corner_SW[0] - corner_SE[0] + corner_NW[0] - corner_NE[0] -
302 corner_NE[0] * S[1] + corner_NW[0] * S[1] + corner_SE[0] * S[1] -
303 corner_SW[0] * S[1]) +
304 0.5 * (dbound_S[0] + dbound_N[0] - bound_W[0] + bound_E[0] -
305 dbound_S[0] * S[1] + dbound_N[0] * S[1]);
306 // dr1/dm0
307 jacobian(0, 1) =
308 0.25 * (corner_SW[1] - corner_SE[1] + corner_NW[1] - corner_NE[1] -
309 corner_NE[1] * S[1] + corner_NW[1] * S[1] + corner_SE[1] * S[1] -
310 corner_SW[1] * S[1]) +
311 0.5 * (dbound_S[1] + dbound_N[1] - bound_W[1] + bound_E[1] -
312 dbound_S[1] * S[1] + dbound_N[1] * S[1]);
313 // dr0/dm1
314 jacobian(1, 0) =
315 0.25 * (corner_SW[0] + corner_SE[0] - corner_NW[0] - corner_NE[0] +
316 corner_SE[0] * S[0] - corner_SW[0] * S[0] - corner_NE[0] * S[0] +
317 corner_NW[0] * S[0]) +
318 0.5 * (-bound_S[0] + bound_N[0] + dbound_W[0] + dbound_E[0] -
319 dbound_W[0] * S[0] + dbound_E[0] * S[0]);
320 // dr1/dm1
321 jacobian(1, 1) =
322 0.25 * (corner_SW[1] + corner_SE[1] - corner_NW[1] - corner_NE[1] +
323 corner_SE[1] * S[0] - corner_SW[1] * S[0] - corner_NE[1] * S[0] +
324 corner_NW[1] * S[0]) +
325 0.5 * (-bound_S[1] + bound_N[1] + dbound_W[1] + dbound_E[1] -
326 dbound_W[1] * S[0] + dbound_E[1] * S[0]);
327 }
328
329
330 //=============================================================================
331 /// Assembles the second derivative jacobian of the mapping from the
332 /// macro coordinates to global coordinates x
333 //=============================================================================
334 void QMacroElement<2>::assemble_macro_to_eulerian_jacobian2(
335 const unsigned& t, const Vector<double>& S, DenseMatrix<double>& jacobian2)
336 {
337 using namespace QuadTreeNames;
338
339 Vector<double> bound_N(2);
340 Vector<double> bound_S(2);
341 Vector<double> bound_W(2);
342 Vector<double> bound_E(2);
343
344 Vector<double> dbound_N(2);
345 Vector<double> dbound_S(2);
346 Vector<double> dbound_E(2);
347 Vector<double> dbound_W(2);
348
349 Vector<double> d2bound_N(2);
350 Vector<double> d2bound_S(2);
351 Vector<double> d2bound_E(2);
352 Vector<double> d2bound_W(2);
353
354 Vector<double> corner_SE(2);
355 Vector<double> corner_SW(2);
356 Vector<double> corner_NE(2);
357 Vector<double> corner_NW(2);
358
359 Vector<double> zeta(1);
360
361
362 // Determine Vectors to corners
363 zeta[0] = 1.0;
364 Domain_pt->macro_element_boundary(
365 t, Macro_element_number, QuadTreeNames::S, zeta, corner_SE);
366 zeta[0] = -1.0;
367 Domain_pt->macro_element_boundary(
368 t, Macro_element_number, QuadTreeNames::S, zeta, corner_SW);
369 zeta[0] = 1.0;
370 Domain_pt->macro_element_boundary(
371 t, Macro_element_number, QuadTreeNames::N, zeta, corner_NE);
372 zeta[0] = -1.0;
373 Domain_pt->macro_element_boundary(
374 t, Macro_element_number, QuadTreeNames::N, zeta, corner_NW);
375
376
377 // Get the position and first derivativeson the N/S/W/E edges
378 zeta[0] = S[0];
379 Domain_pt->macro_element_boundary(
380 t, Macro_element_number, QuadTreeNames::N, zeta, bound_N);
381 Domain_pt->dmacro_element_boundary(
382 t, Macro_element_number, QuadTreeNames::N, zeta, dbound_N);
383 Domain_pt->d2macro_element_boundary(
384 t, Macro_element_number, QuadTreeNames::N, zeta, d2bound_N);
385 zeta[0] = S[0];
386 Domain_pt->macro_element_boundary(
387 t, Macro_element_number, QuadTreeNames::S, zeta, bound_S);
388 Domain_pt->dmacro_element_boundary(
389 t, Macro_element_number, QuadTreeNames::S, zeta, dbound_S);
390 Domain_pt->d2macro_element_boundary(
391 t, Macro_element_number, QuadTreeNames::S, zeta, d2bound_S);
392 zeta[0] = S[1];
393 Domain_pt->macro_element_boundary(
394 t, Macro_element_number, QuadTreeNames::W, zeta, bound_W);
395 Domain_pt->dmacro_element_boundary(
396 t, Macro_element_number, QuadTreeNames::W, zeta, dbound_W);
397 Domain_pt->d2macro_element_boundary(
398 t, Macro_element_number, QuadTreeNames::W, zeta, d2bound_W);
399 zeta[0] = S[1];
400 Domain_pt->macro_element_boundary(
401 t, Macro_element_number, QuadTreeNames::E, zeta, bound_E);
402 Domain_pt->dmacro_element_boundary(
403 t, Macro_element_number, QuadTreeNames::E, zeta, dbound_E);
404 Domain_pt->d2macro_element_boundary(
405 t, Macro_element_number, QuadTreeNames::E, zeta, d2bound_E);
406
407
408 // d2x0/dm0^2
409 jacobian2(0, 0) = 0.5 * (d2bound_S[0] + d2bound_N[0] - d2bound_S[0] * S[1] +
410 d2bound_N[0] * S[1]);
411 // d2x0/dm1^2
412 jacobian2(1, 0) = 0.5 * (d2bound_W[0] + d2bound_E[0] - d2bound_W[0] * S[0] +
413 d2bound_E[0] * S[0]);
414 // d2x0/dm0dm1
415 jacobian2(2, 0) =
416 0.25 * (-corner_NE[0] + corner_NW[0] + corner_SE[0] - corner_SW[0]) +
417 0.5 * (-dbound_W[0] + dbound_E[0] - dbound_S[0] + dbound_N[0]);
418 // d2x1/dm0^2
419 jacobian2(0, 1) = 0.5 * (d2bound_S[1] + d2bound_N[1] - d2bound_S[1] * S[1] +
420 d2bound_N[1] * S[1]);
421 // d2x1/dm1^2
422 jacobian2(1, 1) = 0.5 * (d2bound_W[1] + d2bound_E[1] - d2bound_W[1] * S[0] +
423 d2bound_E[1] * S[0]);
424 // d2x1/dm0dm1
425 jacobian2(2, 1) =
426 0.25 * (-corner_NE[1] + corner_NW[1] + corner_SE[1] - corner_SW[1]) +
427 0.5 * (-dbound_W[1] + dbound_E[1] - dbound_S[1] + dbound_N[1]);
428 }
429
430
431 //////////////////////////////////////////////////////////////////////
432 //////////////////////////////////////////////////////////////////////
433 //////////////////////////////////////////////////////////////////////
434
435
436 //=================================================================
437 /// Get global position r(S) at discrete time level t.
438 /// t=0: Present time; t>0: previous timestep.
439 //=================================================================
440 void QMacroElement<3>::macro_map(const unsigned& t,
441 const Vector<double>& S,
442 Vector<double>& r)
443 {
444 // get the eight corners
445 Vector<double> corner_LDB(3);
446 Vector<double> corner_RDB(3);
447 Vector<double> corner_LUB(3);
448 Vector<double> corner_RUB(3);
449 Vector<double> corner_LDF(3);
450 Vector<double> corner_RDF(3);
451 Vector<double> corner_LUF(3);
452 Vector<double> corner_RUF(3);
453
454 Vector<double> zeta(2);
455 zeta[0] = -1.0;
456 zeta[1] = -1.0;
457 Domain_pt->macro_element_boundary(
458 t, Macro_element_number, OcTreeNames::B, zeta, corner_LDB);
459 Domain_pt->macro_element_boundary(
460 t, Macro_element_number, OcTreeNames::U, zeta, corner_LUB);
461 Domain_pt->macro_element_boundary(
462 t, Macro_element_number, OcTreeNames::F, zeta, corner_LDF);
463 Domain_pt->macro_element_boundary(
464 t, Macro_element_number, OcTreeNames::R, zeta, corner_RDB);
465 zeta[0] = 1.0;
466 zeta[1] = 1.0;
467 Domain_pt->macro_element_boundary(
468 t, Macro_element_number, OcTreeNames::B, zeta, corner_RUB);
469 Domain_pt->macro_element_boundary(
470 t, Macro_element_number, OcTreeNames::D, zeta, corner_RDF);
471 Domain_pt->macro_element_boundary(
472 t, Macro_element_number, OcTreeNames::L, zeta, corner_LUF);
473 Domain_pt->macro_element_boundary(
474 t, Macro_element_number, OcTreeNames::R, zeta, corner_RUF);
475
476
477 // get the position of the 4 corners of the center slice
478 Vector<double> corner_LD(3);
479 Vector<double> corner_RD(3);
480 Vector<double> corner_LU(3);
481 Vector<double> corner_RU(3);
482
483 zeta[0] = -1.0;
484 zeta[1] = S[2];
485 Domain_pt->macro_element_boundary(
486 t, Macro_element_number, OcTreeNames::D, zeta, corner_LD);
487 Domain_pt->macro_element_boundary(
488 t, Macro_element_number, OcTreeNames::U, zeta, corner_LU);
489 zeta[0] = 1.0;
490 Domain_pt->macro_element_boundary(
491 t, Macro_element_number, OcTreeNames::D, zeta, corner_RD);
492 Domain_pt->macro_element_boundary(
493 t, Macro_element_number, OcTreeNames::U, zeta, corner_RU);
494
495 // get position on the B,F faces;
496 Vector<double> face_B(3);
497 Vector<double> face_F(3);
498
499 zeta[0] = S[0];
500 zeta[1] = S[1];
501 Domain_pt->macro_element_boundary(
502 t, Macro_element_number, OcTreeNames::B, zeta, face_B);
503 Domain_pt->macro_element_boundary(
504 t, Macro_element_number, OcTreeNames::F, zeta, face_F);
505
506
507 // get position on the edges of the middle slice
508 Vector<double> edge_mid_L(3);
509 Vector<double> edge_mid_R(3);
510 Vector<double> edge_mid_D(3);
511 Vector<double> edge_mid_U(3);
512 zeta[0] = S[0];
513 zeta[1] = S[2];
514 Domain_pt->macro_element_boundary(
515 t, Macro_element_number, OcTreeNames::U, zeta, edge_mid_U);
516
517 Domain_pt->macro_element_boundary(
518 t, Macro_element_number, OcTreeNames::D, zeta, edge_mid_D);
519 zeta[0] = S[1];
520 zeta[1] = S[2];
521 Domain_pt->macro_element_boundary(
522 t, Macro_element_number, OcTreeNames::L, zeta, edge_mid_L);
523
524 Domain_pt->macro_element_boundary(
525 t, Macro_element_number, OcTreeNames::R, zeta, edge_mid_R);
526
527 // get position on the edges of the back slice
528 Vector<double> edge_back_L(3);
529 Vector<double> edge_back_R(3);
530 Vector<double> edge_back_D(3);
531 Vector<double> edge_back_U(3);
532 zeta[0] = S[0];
533 zeta[1] = -1.0;
534 Domain_pt->macro_element_boundary(
535 t, Macro_element_number, OcTreeNames::U, zeta, edge_back_U);
536
537 Domain_pt->macro_element_boundary(
538 t, Macro_element_number, OcTreeNames::D, zeta, edge_back_D);
539 zeta[0] = S[1];
540 zeta[1] = -1.0;
541 Domain_pt->macro_element_boundary(
542 t, Macro_element_number, OcTreeNames::L, zeta, edge_back_L);
543
544 Domain_pt->macro_element_boundary(
545 t, Macro_element_number, OcTreeNames::R, zeta, edge_back_R);
546
547 // get position on the edges of the front slice
548 Vector<double> edge_front_L(3);
549 Vector<double> edge_front_R(3);
550 Vector<double> edge_front_D(3);
551 Vector<double> edge_front_U(3);
552 zeta[0] = S[0];
553 zeta[1] = 1.0;
554 Domain_pt->macro_element_boundary(
555 t, Macro_element_number, OcTreeNames::U, zeta, edge_front_U);
556
557 Domain_pt->macro_element_boundary(
558 t, Macro_element_number, OcTreeNames::D, zeta, edge_front_D);
559 zeta[0] = S[1];
560 zeta[1] = 1.0;
561 Domain_pt->macro_element_boundary(
562 t, Macro_element_number, OcTreeNames::L, zeta, edge_front_L);
563
564 Domain_pt->macro_element_boundary(
565 t, Macro_element_number, OcTreeNames::R, zeta, edge_front_R);
566
567
568 for (unsigned i = 0; i < 3; i++)
569 {
570 // Position on the middle slice
571 // ============================
572 double slice_mid;
573
574 // the points on the up and down edges of the middle "rectangular slice"
575 double ptUp, ptDo;
576 ptUp = corner_LU[i] + (corner_RU[i] - corner_LU[i]) * 0.5 * (S[0] + 1.0);
577 ptDo = corner_LD[i] + (corner_RD[i] - corner_LD[i]) * 0.5 * (S[0] + 1.0);
578 // position in the rectangular middle slice
579 slice_mid = ptDo + 0.5 * (1.0 + S[1]) * (ptUp - ptDo);
580
581 // get the differences to the edges of the middle slice
582 double diff_L, diff_R, diff_D, diff_U;
583 diff_L = edge_mid_L[i] - slice_mid;
584 diff_R = edge_mid_R[i] - slice_mid;
585 diff_D = edge_mid_D[i] - slice_mid;
586 diff_U = edge_mid_U[i] - slice_mid;
587
588 // Map it to get the position in the middle slice
589 slice_mid +=
590 diff_L * (1.0 - 0.5 * (S[0] + 1.0)) + diff_R * 0.5 * (S[0] + 1.0) +
591 diff_D * (1.0 - 0.5 * (S[1] + 1.0)) + diff_U * 0.5 * (S[1] + 1.0);
592
593
594 // Position on the back slice
595 //===========================
596 double slice_back;
597
598 // the points on the up and down edges of the back "rectangular slice"
599 ptUp =
600 corner_LUB[i] + (corner_RUB[i] - corner_LUB[i]) * 0.5 * (S[0] + 1.0);
601 ptDo =
602 corner_LDB[i] + (corner_RDB[i] - corner_LDB[i]) * 0.5 * (S[0] + 1.0);
603 // position in the rectangular back slice
604 slice_back = ptDo + 0.5 * (1.0 + S[1]) * (ptUp - ptDo);
605
606 // get the differences to the edges of the middle slice
607 diff_L = edge_back_L[i] - slice_back;
608 diff_R = edge_back_R[i] - slice_back;
609 diff_D = edge_back_D[i] - slice_back;
610 diff_U = edge_back_U[i] - slice_back;
611
612 // Map it to get the position in the back slice
613 slice_back +=
614 diff_L * (1.0 - 0.5 * (S[0] + 1.0)) + diff_R * 0.5 * (S[0] + 1.0) +
615 diff_D * (1.0 - 0.5 * (S[1] + 1.0)) + diff_U * 0.5 * (S[1] + 1.0);
616
617 // Position on the front slice
618 //============================
619 double slice_front;
620
621 // the points on the up and down edges of the back "rectangular slice"
622 ptUp =
623 corner_LUF[i] + (corner_RUF[i] - corner_LUF[i]) * 0.5 * (S[0] + 1.0);
624 ptDo =
625 corner_LDF[i] + (corner_RDF[i] - corner_LDF[i]) * 0.5 * (S[0] + 1.0);
626 // position in the rectangular back slice
627 slice_front = ptDo + 0.5 * (1.0 + S[1]) * (ptUp - ptDo);
628
629 // get the differences to the edges of the middle slice
630 diff_L = edge_front_L[i] - slice_front;
631 diff_R = edge_front_R[i] - slice_front;
632 diff_D = edge_front_D[i] - slice_front;
633 diff_U = edge_front_U[i] - slice_front;
634
635 // Map it to get the position in the back slice
636 slice_front +=
637 diff_L * (1.0 - 0.5 * (S[0] + 1.0)) + diff_R * 0.5 * (S[0] + 1.0) +
638 diff_D * (1.0 - 0.5 * (S[1] + 1.0)) + diff_U * 0.5 * (S[1] + 1.0);
639
640 // Get difference between the back and front slices and the actual
641 // boundary
642 // ========================================================================
643
644 double diff_back = face_B[i] - slice_back;
645 double diff_front = face_F[i] - slice_front;
646
647 // final map
648 //==========
649
650 r[i] = slice_mid + 0.5 * (1 + S[2]) * diff_front +
651 0.5 * (1 - S[2]) * diff_back;
652 }
653 }
654
655
656 //=================================================================
657 /// Output all macro element boundaries as tecplot zones
658 //=================================================================
659 void QMacroElement<3>::output_macro_element_boundaries(std::ostream& outfile,
660 const unsigned& nplot)
661 {
662 using namespace OcTreeNames;
663
664 Vector<double> s(2);
665 Vector<double> f(3);
666 // Dump at present time
667 unsigned t = 0;
668 for (unsigned idirect = L; idirect <= F; idirect++)
669 {
670 outfile << "ZONE I=" << nplot << ", J=" << nplot << std::endl;
671 for (unsigned i = 0; i < nplot; i++)
672 {
673 s[1] = -1.0 + 2.0 * double(i) / double(nplot - 1);
674 for (unsigned j = 0; j < nplot; j++)
675 {
676 s[0] = -1.0 + 2.0 * double(j) / double(nplot - 1);
677 Domain_pt->macro_element_boundary(
678 t, Macro_element_number, idirect, s, f);
679 outfile << f[0] << " " << f[1] << " " << f[2] << std::endl;
680 }
681 }
682 }
683 }
684
685} // namespace oomph
s
static char t char * s
Definition cfortran.h:568
i
cstr elem_len * i
Definition cfortran.h:603
t
char t
Definition cfortran.h:568
oomph::MacroElement::assemble_macro_to_eulerian_jacobian2
virtual void assemble_macro_to_eulerian_jacobian2(const unsigned &t, const Vector< double > &s, DenseMatrix< double > &jacobian2)
Assembles the second derivative jacobian of the mapping from the macro coordinates to the global coor...
Definition macro_element.h:186
oomph::MacroElement::output_macro_element_boundaries
virtual void output_macro_element_boundaries(std::ostream &outfile, const unsigned &nplot)=0
Output all macro element boundaries as tecplot zones.
oomph::MacroElement::macro_map
void macro_map(const Vector< double > &s, Vector< double > &r)
The mapping from local to global coordinates at the current time : r(s)
Definition macro_element.h:126
oomph::MacroElement::assemble_macro_to_eulerian_jacobian
virtual void assemble_macro_to_eulerian_jacobian(const unsigned &t, const Vector< double > &s, DenseMatrix< double > &jacobian)
the jacobian of the mapping from the macro coordinates to the global coordinates
Definition macro_element.h:170
oomph::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Definition Tadvection_diffusion_reaction_elements.h:66
oomph::TAdvectionDiffusionReactionElement::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement()
Constructor: Call constructors for TElement and AdvectionDiffusionReaction equations.
Definition Tadvection_diffusion_reaction_elements.h:70
macro_element.h
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition advection_diffusion_elements.cc:30