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-2023 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 
29 namespace 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::Domain::macro_element_boundary
virtual void macro_element_boundary(const unsigned &t, const unsigned &i_macro, const unsigned &i_direct, const Vector< double > &s, Vector< double > &f)=0
Vector representation of the i_macro-th macro element boundary i_direct (e.g. N/S/W/E in 2D) at time ...
oomph::Domain::dmacro_element_boundary
virtual void dmacro_element_boundary(const unsigned &t, const unsigned &i_macro, const unsigned &i_direct, const Vector< double > &s, Vector< double > &f)
Vector representation of the i_macro-th macro element boundary derivatives i_direct (e....
Definition: domain.h:212
oomph::Domain::d2macro_element_boundary
virtual void d2macro_element_boundary(const unsigned &t, const unsigned &i_macro, const unsigned &i_direct, const Vector< double > &s, Vector< double > &f)
Vector representation of the i_macro-th macro element boundary second derivatives i_direct (e....
Definition: domain.h:258
oomph::MacroElement::Domain_pt
Domain * Domain_pt
Pointer to domain.
Definition: macro_element.h:236
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::Macro_element_number
unsigned Macro_element_number
What is the number of the current macro element within its domain.
Definition: macro_element.h:239
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
macro_element.h
oomph
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: advection_diffusion_elements.cc:30