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