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-2022 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
28#include "extruded_domain.h"
29
30namespace oomph
31{
32 //=================================================================
33 /// Get global position r(s) at the continuous time, t
34 //=================================================================
36 const Vector<double>& s,
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....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...