extruded_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 // Necessary oomph-lib headers
27 #include "extruded_macro_element.h"
28 #include "extruded_domain.h"
29 
30 namespace oomph
31 {
32  //=================================================================
33  /// Get global position r(s) at the continuous time, t
34  //=================================================================
35  void QExtrudedMacroElement<3>::macro_map(const unsigned& t,
36  const Vector<double>& s,
37  Vector<double>& r)
38  {
39  // Make sure that t=0 otherwise this doesn't make sense
40  if (t != 0)
41  {
42  // Create an output stream
43  std::ostringstream error_message_stream;
44 
45  // Create an error message
46  error_message_stream << "This output function outputs a space-time\n"
47  << "representation of the solution. As such, it\n"
48  << "does not make sense to output the solution\n"
49  << "at a previous time level!" << std::endl;
50 
51  // Throw an error
52  throw OomphLibError(error_message_stream.str(),
53  OOMPH_CURRENT_FUNCTION,
54  OOMPH_EXCEPTION_LOCATION);
55  }
56 
57  // Storage for the eight corners
58  Vector<double> corner_LDB(3, 0.0);
59  Vector<double> corner_RDB(3, 0.0);
60  Vector<double> corner_LUB(3, 0.0);
61  Vector<double> corner_RUB(3, 0.0);
62  Vector<double> corner_LDF(3, 0.0);
63  Vector<double> corner_RDF(3, 0.0);
64  Vector<double> corner_LUF(3, 0.0);
65  Vector<double> corner_RUF(3, 0.0);
66 
67  // Lagrangian coordinates of a point on a 2D surface in 3D space
68  Vector<double> zeta(2, 0.0);
69 
70  // First coordinate of the bottom-left coordinates of a surface
71  zeta[0] = -1.0;
72 
73  // Second coordinate of the bottom-left coordinates of a surface
74  zeta[1] = -1.0;
75 
76  // Calculate the space-time coordinates of the LDB corner
78  t, Macro_element_number, OcTreeNames::B, zeta, corner_LDB);
79 
80  // Calculate the space-time coordinates of the LUB corner
82  t, Macro_element_number, OcTreeNames::U, zeta, corner_LUB);
83 
84  // Calculate the space-time coordinates of the LDF corner
86  t, Macro_element_number, OcTreeNames::F, zeta, corner_LDF);
87 
88  // Calculate the space-time coordinates of the RDB corner
90  t, Macro_element_number, OcTreeNames::R, zeta, corner_RDB);
91 
92  // First coordinate of the the top-right coordinates of a surface
93  zeta[0] = 1.0;
94 
95  // Second coordinate of the the top-right coordinates of a surface
96  zeta[1] = 1.0;
97 
98  // Calculate the space-time coordinates of the RUB corner
100  t, Macro_element_number, OcTreeNames::B, zeta, corner_RUB);
101 
102  // Calculate the space-time coordinates of the RDF corner
104  t, Macro_element_number, OcTreeNames::D, zeta, corner_RDF);
105 
106  // Calculate the space-time coordinates of the LUF corner
108  t, Macro_element_number, OcTreeNames::L, zeta, corner_LUF);
109 
110  // Calculate the space-time coordinates of the RUF corner
112  t, Macro_element_number, OcTreeNames::R, zeta, corner_RUF);
113 
114  // Get the position of the 4 corners of the center slice
115  Vector<double> corner_LD(3, 0.0);
116  Vector<double> corner_RD(3, 0.0);
117  Vector<double> corner_LU(3, 0.0);
118  Vector<double> corner_RU(3, 0.0);
119 
120  // Set the coordinates of the point on a surface
121  zeta[0] = -1.0;
122  zeta[1] = s[2];
123 
125  t, Macro_element_number, OcTreeNames::D, zeta, corner_LD);
127  t, Macro_element_number, OcTreeNames::U, zeta, corner_LU);
128  zeta[0] = 1.0;
130  t, Macro_element_number, OcTreeNames::D, zeta, corner_RD);
132  t, Macro_element_number, OcTreeNames::U, zeta, corner_RU);
133 
134  // Get position on the B,F faces;
135  Vector<double> face_B(3);
136  Vector<double> face_F(3);
137 
138  zeta[0] = s[0];
139  zeta[1] = s[1];
141  t, Macro_element_number, OcTreeNames::B, zeta, face_B);
143  t, Macro_element_number, OcTreeNames::F, zeta, face_F);
144 
145 
146  // Get position on the edges of the middle slice
147  Vector<double> edge_mid_L(3);
148  Vector<double> edge_mid_R(3);
149  Vector<double> edge_mid_D(3);
150  Vector<double> edge_mid_U(3);
151  zeta[0] = s[0];
152  zeta[1] = s[2];
154  t, Macro_element_number, OcTreeNames::U, zeta, edge_mid_U);
155 
157  t, Macro_element_number, OcTreeNames::D, zeta, edge_mid_D);
158  zeta[0] = s[1];
159  zeta[1] = s[2];
161  t, Macro_element_number, OcTreeNames::L, zeta, edge_mid_L);
162 
164  t, Macro_element_number, OcTreeNames::R, zeta, edge_mid_R);
165 
166  // Get position on the edges of the back slice
167  Vector<double> edge_back_L(3);
168  Vector<double> edge_back_R(3);
169  Vector<double> edge_back_D(3);
170  Vector<double> edge_back_U(3);
171  zeta[0] = s[0];
172  zeta[1] = -1.0;
174  t, Macro_element_number, OcTreeNames::U, zeta, edge_back_U);
175 
177  t, Macro_element_number, OcTreeNames::D, zeta, edge_back_D);
178  zeta[0] = s[1];
179  zeta[1] = -1.0;
181  t, Macro_element_number, OcTreeNames::L, zeta, edge_back_L);
182 
184  t, Macro_element_number, OcTreeNames::R, zeta, edge_back_R);
185 
186  // Get position on the edges of the front slice
187  Vector<double> edge_front_L(3);
188  Vector<double> edge_front_R(3);
189  Vector<double> edge_front_D(3);
190  Vector<double> edge_front_U(3);
191  zeta[0] = s[0];
192  zeta[1] = 1.0;
194  t, Macro_element_number, OcTreeNames::U, zeta, edge_front_U);
195 
197  t, Macro_element_number, OcTreeNames::D, zeta, edge_front_D);
198  zeta[0] = s[1];
199  zeta[1] = 1.0;
201  t, Macro_element_number, OcTreeNames::L, zeta, edge_front_L);
202 
204  t, Macro_element_number, OcTreeNames::R, zeta, edge_front_R);
205 
206  // The number of dimensions (=space+time)
207  unsigned n_dim = 3;
208 
209  // Loop over the coordinate directions
210  for (unsigned i = 0; i < n_dim; i++)
211  {
212  //-----------------------------
213  // Position on the middle slice
214  //-----------------------------
215  double slice_mid;
216 
217  // The points on the up and down edges of the middle "rectangular slice"
218  double point_up = 0.0;
219  double point_down = 0.0;
220 
221  // Calculate the point on the upper edge
222  point_up =
223  corner_LU[i] + (corner_RU[i] - corner_LU[i]) * 0.5 * (s[0] + 1.0);
224 
225  // Calculate the point on the lower edge
226  point_down =
227  corner_LD[i] + (corner_RD[i] - corner_LD[i]) * 0.5 * (s[0] + 1.0);
228 
229  // Position in the rectangular middle slice
230  slice_mid = point_down + 0.5 * (1.0 + s[1]) * (point_up - point_down);
231 
232  // Get the differences to the edges of the middle slice
233  double diff_L, diff_R, diff_D, diff_U;
234  diff_L = edge_mid_L[i] - slice_mid;
235  diff_R = edge_mid_R[i] - slice_mid;
236  diff_D = edge_mid_D[i] - slice_mid;
237  diff_U = edge_mid_U[i] - slice_mid;
238 
239  // Map it to get the position in the middle slice
240  slice_mid +=
241  (diff_L * (1.0 - 0.5 * (s[0] + 1.0)) + diff_R * 0.5 * (s[0] + 1.0) +
242  diff_D * (1.0 - 0.5 * (s[1] + 1.0)) + diff_U * 0.5 * (s[1] + 1.0));
243 
244  //---------------------------
245  // Position on the back slice
246  //---------------------------
247  double slice_back;
248 
249  // Calculate the point on the upper edge of the back "rectangular slice"
250  point_up =
251  corner_LUB[i] + (corner_RUB[i] - corner_LUB[i]) * 0.5 * (s[0] + 1.0);
252 
253  // Calculate the point on the lower edge of the back "rectangular slice"
254  point_down =
255  corner_LDB[i] + (corner_RDB[i] - corner_LDB[i]) * 0.5 * (s[0] + 1.0);
256 
257  // Position in the rectangular back slice
258  slice_back = point_down + 0.5 * (1.0 + s[1]) * (point_up - point_down);
259 
260  // Get the differences to the edges of the middle slice
261  diff_L = edge_back_L[i] - slice_back;
262  diff_R = edge_back_R[i] - slice_back;
263  diff_D = edge_back_D[i] - slice_back;
264  diff_U = edge_back_U[i] - slice_back;
265 
266  // Map it to get the position in the back slice
267  slice_back +=
268  (diff_L * (1.0 - 0.5 * (s[0] + 1.0)) + diff_R * 0.5 * (s[0] + 1.0) +
269  diff_D * (1.0 - 0.5 * (s[1] + 1.0)) + diff_U * 0.5 * (s[1] + 1.0));
270 
271  //----------------------------
272  // Position on the front slice
273  //----------------------------
274  double slice_front;
275 
276  // Calculate the point on the upper edge of the back "rectangular slice"
277  point_up =
278  corner_LUF[i] + (corner_RUF[i] - corner_LUF[i]) * 0.5 * (s[0] + 1.0);
279 
280  // Calculate the point on the lower edge of the back "rectangular slice"
281  point_down =
282  corner_LDF[i] + (corner_RDF[i] - corner_LDF[i]) * 0.5 * (s[0] + 1.0);
283 
284  // Position in the rectangular back slice
285  slice_front = point_down + 0.5 * (1.0 + s[1]) * (point_up - point_down);
286 
287  // Get the differences to the edges of the middle slice
288  diff_L = edge_front_L[i] - slice_front;
289  diff_R = edge_front_R[i] - slice_front;
290  diff_D = edge_front_D[i] - slice_front;
291  diff_U = edge_front_U[i] - slice_front;
292 
293  // Map it to get the position in the back slice
294  slice_front +=
295  (diff_L * (1.0 - 0.5 * (s[0] + 1.0)) + diff_R * 0.5 * (s[0] + 1.0) +
296  diff_D * (1.0 - 0.5 * (s[1] + 1.0)) + diff_U * 0.5 * (s[1] + 1.0));
297 
298  //-------------------------------------------------------------------------
299  // Get difference between the back and front slices and the actual
300  // boundary
301  //-------------------------------------------------------------------------
302  double diff_back = face_B[i] - slice_back;
303  double diff_front = face_F[i] - slice_front;
304 
305  //----------
306  // Final map
307  //----------
308  r[i] = slice_mid + 0.5 * (1.0 + s[2]) * diff_front +
309  0.5 * (1.0 - s[2]) * diff_back;
310  }
311  } // End of macro_map
312 
313  //=================================================================
314  /// Output all macro element boundaries as tecplot zones
315  //=================================================================
317  std::ostream& outfile, const unsigned& nplot)
318  {
319  // Use the OcTree enumeration for corners, edges and faces
320  using namespace OcTreeNames;
321 
322  // Storage for the local coordinates (of a point on a face)
323  Vector<double> s(2);
324 
325  // Storage for the global coordinates
326  Vector<double> x(3);
327 
328  // Loop over the faces
329  for (unsigned idirect = L; idirect <= F; idirect++)
330  {
331  // Output the header
332  outfile << "ZONE I=" << nplot << ", J=" << nplot << std::endl;
333 
334  // Loop over the plot points in the second surface direction
335  for (unsigned i = 0; i < nplot; i++)
336  {
337  // Calculate the second local coordinate associated with this plot point
338  s[1] = -1.0 + 2.0 * double(i) / double(nplot - 1);
339 
340  // Loop over the plot points in the first surface direction
341  for (unsigned j = 0; j < nplot; j++)
342  {
343  // Calculate the first local coordinate associated with this plot
344  // point
345  s[0] = -1.0 + 2.0 * double(j) / double(nplot - 1);
346 
347  // To make the extrusion machinery consistent with the Domain
348  // machinery a time level argument has to be provided. For our
349  // purposes we set this to zero to ensure that the appropriate output
350  // is provided.
351  unsigned t = 0;
352 
353  // Get the global coordinates associated with these local coordinates
355  t, Macro_element_number, idirect, s, x);
356 
357  // Output the global (space-time) coordinates
358  outfile << x[0] << " " << x[1] << " " << x[2] << std::endl;
359  }
360  } // for (unsigned i=0;i<nplot;i++)
361  } // for (unsigned idirect=L;idirect<=F;idirect++)
362  } // End of output_macro_element_boundaries
363 } // End of namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
void macro_element_boundary(const unsigned &time, const unsigned &i_macro, const unsigned &i_direct, const Vector< double > &s, Vector< double > &x)
Vector representation of the i_macro-th macro element boundary i_direct (e.g. N/S/W/E in 2D spatial =...
ExtrudedDomain * Extruded_domain_pt
Pointer to the extruded domain.
virtual void output_macro_element_boundaries(std::ostream &outfile, const unsigned &nplot)=0
Output all macro element boundaries as tecplot zones.
void macro_map(const Vector< double > &s, Vector< double > &r)
The mapping from local to global coordinates at the current time : r(s)
unsigned Macro_element_number
What is the number of the current macro element within its domain.
An OomphLibError object which should be thrown when an run-time error is encountered....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...