eighth_sphere_mesh.template.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-2024 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#ifndef OOMPH_EIGHTH_SPHERE_MESH_TEMPLATE_CC
27#define OOMPH_EIGHTH_SPHERE_MESH_TEMPLATE_CC
28
29#include "eighth_sphere_mesh.template.h"
30
31
32namespace oomph
33{
34 //======================================================================
35 /// Constructor for the eighth of a sphere mesh: Pass timestepper;
36 /// defaults to static default timestepper.
37 //======================================================================
38 template<class ELEMENT>
39 EighthSphereMesh<ELEMENT>::EighthSphereMesh(const double& radius,
40 TimeStepper* time_stepper_pt)
41 : Radius(radius)
42 {
43 // Mesh can only be built with 3D Qelements.
44 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
45
46 // Set the number of boundaries
47 this->set_nboundary(4);
48
49 // Provide storage for the four elements
50 this->Element_pt.resize(4);
51
52 // Set the domain pointer: Pass the radius of the sphere
53 Domain_pt = new EighthSphereDomain(Radius);
54
55
56 Vector<double> s(3), s_fraction(3);
57 Vector<double> r(3);
58
59 // Create first element
60 //--------------------
61 this->Element_pt[0] = new ELEMENT;
62
63 // Give it nodes
64
65 // How many 1D nodes are there?
66 unsigned nnode1d = this->finite_element_pt(0)->nnode_1d();
67
68 // Loop over nodes
69 for (unsigned i = 0; i < nnode1d; i++)
70 {
71 for (unsigned j = 0; j < nnode1d; j++)
72 {
73 for (unsigned k = 0; k < nnode1d; k++)
74 {
75 unsigned jnod = k * nnode1d * nnode1d + j * nnode1d + i;
76
77 // If we're on a boundary, make a boundary node
78 if ((i == 0) || (j == 0) || (k == 0))
79 {
80 // Allocate the node according to the element's wishes
81 this->Node_pt.push_back(
82 this->finite_element_pt(0)->construct_boundary_node(
83 jnod, time_stepper_pt));
84 }
85 // Otherwise it's a normal node
86 else
87 {
88 // Allocate the node according to the element's wishes
89 this->Node_pt.push_back(this->finite_element_pt(0)->construct_node(
90 jnod, time_stepper_pt));
91 }
92
93 // Work out the node's coordinates in the finite element's local
94 // coordinate system
95 this->finite_element_pt(0)->local_fraction_of_node(jnod, s_fraction);
96
97
98 // Get the position of the node from macro element mapping
99 s[0] = -1.0 + 2.0 * s_fraction[0];
100 s[1] = -1.0 + 2.0 * s_fraction[1];
101 s[2] = -1.0 + 2.0 * s_fraction[2];
102 Domain_pt->macro_element_pt(0)->macro_map(s, r);
103
104 // Assign coordinates
105 this->finite_element_pt(0)->node_pt(jnod)->x(0) = r[0];
106 this->finite_element_pt(0)->node_pt(jnod)->x(1) = r[1];
107 this->finite_element_pt(0)->node_pt(jnod)->x(2) = r[2];
108
109 // Add the node to the appropriate boundary
110 if (i == 0)
111 add_boundary_node(0, this->finite_element_pt(0)->node_pt(jnod));
112 if (j == 0)
113 add_boundary_node(1, this->finite_element_pt(0)->node_pt(jnod));
114 if (k == 0)
115 add_boundary_node(2, this->finite_element_pt(0)->node_pt(jnod));
116 }
117 }
118 }
119
120
121 // Create a second element
122 //------------------------
123 this->Element_pt[1] = new ELEMENT;
124 ;
125
126 // Give it nodes
127
128 // Loop over nodes
129 for (unsigned i = 0; i < nnode1d; i++)
130 {
131 for (unsigned j = 0; j < nnode1d; j++)
132 {
133 for (unsigned k = 0; k < nnode1d; k++)
134 {
135 unsigned jnod = k * nnode1d * nnode1d + j * nnode1d + i;
136
137 // If node has not been yet created, create it
138 if (i > 0)
139 {
140 // If the node is on a boundary, make a boundary node
141 if ((i == nnode1d - 1) || (j == 0) || (k == 0))
142 {
143 this->Node_pt.push_back(
144 this->finite_element_pt(1)->construct_boundary_node(
145 jnod, time_stepper_pt));
146 }
147 // Otherwise make a normal node
148 else
149 {
150 this->Node_pt.push_back(
151 this->finite_element_pt(1)->construct_node(jnod,
152 time_stepper_pt));
153 }
154
155 // Work out the node's coordinates in the finite element's local
156 // coordinate system
157 this->finite_element_pt(1)->local_fraction_of_node(jnod,
158 s_fraction);
159
160 // Get the position of the node from macro element mapping
161 s[0] = -1.0 + 2.0 * s_fraction[0];
162 s[1] = -1.0 + 2.0 * s_fraction[1];
163 s[2] = -1.0 + 2.0 * s_fraction[2];
164 Domain_pt->macro_element_pt(1)->macro_map(s, r);
165
166 // Assign coordinate
167 this->finite_element_pt(1)->node_pt(jnod)->x(0) = r[0];
168 this->finite_element_pt(1)->node_pt(jnod)->x(1) = r[1];
169 this->finite_element_pt(1)->node_pt(jnod)->x(2) = r[2];
170
171 // Add the node to the approriate boundary
172 if (j == 0)
173 add_boundary_node(1, this->finite_element_pt(1)->node_pt(jnod));
174 if (k == 0)
175 add_boundary_node(2, this->finite_element_pt(1)->node_pt(jnod));
176 if (i == nnode1d - 1)
177 add_boundary_node(3, this->finite_element_pt(1)->node_pt(jnod));
178 }
179
180 // ...else use the node already created
181 else
182 {
183 this->finite_element_pt(1)->node_pt(jnod) =
184 this->finite_element_pt(0)->node_pt(jnod + nnode1d - 1);
185 }
186 }
187 }
188 }
189
190
191 // Create a third element
192 //------------------------
193 this->Element_pt[2] = new ELEMENT;
194
195 // Give it nodes
196
197 // Loop over nodes
198 for (unsigned i = 0; i < nnode1d; i++)
199 {
200 for (unsigned j = 0; j < nnode1d; j++)
201 {
202 for (unsigned k = 0; k < nnode1d; k++)
203 {
204 unsigned jnod = k * nnode1d * nnode1d + j * nnode1d + i;
205
206 // If the node has not been yet created, create it
207 if ((i < nnode1d - 1) && (j > 0))
208 {
209 // If it's on a boundary, make a boundary node
210 if ((i == 0) || (j == nnode1d - 1) || (k == 0))
211 {
212 // Allocate the node according to the element's wishes
213 this->Node_pt.push_back(
214 this->finite_element_pt(2)->construct_boundary_node(
215 jnod, time_stepper_pt));
216 }
217 // Otherwise allocate a normal node
218 else
219 {
220 // Allocate the node according to the element's wishes
221 this->Node_pt.push_back(
222 this->finite_element_pt(2)->construct_node(jnod,
223 time_stepper_pt));
224 }
225
226 // Work out the node's coordinates in the finite element's local
227 // coordinate system
228 this->finite_element_pt(2)->local_fraction_of_node(jnod,
229 s_fraction);
230
231 // Get the position of the node from macro element mapping
232 s[0] = -1.0 + 2.0 * s_fraction[0];
233 s[1] = -1.0 + 2.0 * s_fraction[1];
234 s[2] = -1.0 + 2.0 * s_fraction[2];
235 Domain_pt->macro_element_pt(2)->macro_map(s, r);
236
237 // Assign coordinates
238 this->finite_element_pt(2)->node_pt(jnod)->x(0) = r[0];
239 this->finite_element_pt(2)->node_pt(jnod)->x(1) = r[1];
240 this->finite_element_pt(2)->node_pt(jnod)->x(2) = r[2];
241
242 // Add the node to the appropriate boundary
243 if (i == 0)
244 add_boundary_node(0, this->finite_element_pt(2)->node_pt(jnod));
245 if (k == 0)
246 add_boundary_node(2, this->finite_element_pt(2)->node_pt(jnod));
247 if (j == nnode1d - 1)
248 add_boundary_node(3, this->finite_element_pt(2)->node_pt(jnod));
249 }
250
251 // ...else use the nodes already created
252 else
253 {
254 // If the node belongs to the element 0
255 if (j == 0)
256 this->finite_element_pt(2)->node_pt(jnod) =
257 this->finite_element_pt(0)->node_pt(jnod +
258 nnode1d * (nnode1d - 1));
259
260 // ...else it belongs to the element 1
261 else if (i == nnode1d - 1)
262 this->finite_element_pt(2)->node_pt(jnod) =
263 this->finite_element_pt(1)->node_pt(nnode1d * nnode1d * k + j +
264 i * nnode1d);
265 }
266 }
267 }
268 }
269
270
271 // Create the fourth element
272 //-------------------------
273 this->Element_pt[3] = new ELEMENT;
274
275 // Give it nodes
276
277 // Loop over nodes
278 for (unsigned i = 0; i < nnode1d; i++)
279 {
280 for (unsigned j = 0; j < nnode1d; j++)
281 {
282 for (unsigned k = 0; k < nnode1d; k++)
283 {
284 unsigned jnod = k * nnode1d * nnode1d + j * nnode1d + i;
285
286 // If the node has not been yet created, create it
287 if ((k > 0) && (i < nnode1d - 1) && (j < nnode1d - 1))
288 {
289 // If it's on a boundary, allocate a boundary node
290 if ((i == 0) || (j == 0) || (k == nnode1d - 1))
291 {
292 // Allocate the node according to the element's wishes
293 this->Node_pt.push_back(
294 this->finite_element_pt(3)->construct_boundary_node(
295 jnod, time_stepper_pt));
296 }
297 // Otherwise allocate a normal node
298 else
299 {
300 // Allocate the node according to the element's wishes
301 this->Node_pt.push_back(
302 this->finite_element_pt(3)->construct_node(jnod,
303 time_stepper_pt));
304 }
305
306 // Work out the node's coordinates in the finite element's local
307 // coordinate system
308 this->finite_element_pt(3)->local_fraction_of_node(jnod,
309 s_fraction);
310
311 // Get the position of the node from macro element mapping
312 s[0] = -1.0 + 2.0 * s_fraction[0];
313 s[1] = -1.0 + 2.0 * s_fraction[1];
314 s[2] = -1.0 + 2.0 * s_fraction[2];
315 Domain_pt->macro_element_pt(3)->macro_map(s, r);
316
317 // Assign coordinates
318 this->finite_element_pt(3)->node_pt(jnod)->x(0) = r[0];
319 this->finite_element_pt(3)->node_pt(jnod)->x(1) = r[1];
320 this->finite_element_pt(3)->node_pt(jnod)->x(2) = r[2];
321
322 // Add the node to the appropriate boundary
323 if (i == 0)
324 add_boundary_node(0, this->finite_element_pt(3)->node_pt(jnod));
325 if (j == 0)
326 add_boundary_node(1, this->finite_element_pt(3)->node_pt(jnod));
327 if (k == nnode1d - 1)
328 add_boundary_node(3, this->finite_element_pt(3)->node_pt(jnod));
329 }
330
331 // ...otherwise the node was already created: use it.
332 else
333 {
334 // if k=0 then the node belongs to element 0
335 if (k == 0)
336 {
337 this->finite_element_pt(3)->node_pt(jnod) =
338 this->finite_element_pt(0)->node_pt(jnod + nnode1d * nnode1d *
339 (nnode1d - 1));
340 }
341 else
342 {
343 // else if i==nnode1d-1 the node already exists in element 1
344 if (i == nnode1d - 1)
345 {
346 this->finite_element_pt(3)->node_pt(jnod) =
347 this->finite_element_pt(1)->node_pt(
348 k + i * nnode1d * nnode1d + j * nnode1d);
349 }
350 else
351 // else, the node exists in element 2
352 {
353 this->finite_element_pt(3)->node_pt(jnod) =
354 this->finite_element_pt(2)->node_pt(i + k * nnode1d +
355 j * nnode1d * nnode1d);
356 }
357 }
358 }
359 }
360 }
361 }
362
363 // Setup boundary element lookup schemes
364 setup_boundary_element_info();
365 }
366
367} // namespace oomph
368#endif
oomph::EighthSphereDomain
Eighth sphere as domain. Domain is parametrised by four macro elements.
Definition eighth_sphere_domain.h:42
oomph::EighthSphereMesh::EighthSphereMesh
EighthSphereMesh(const double &radius, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass radius and timestepper; defaults to static default timestepper.
Definition eighth_sphere_mesh.template.cc:39
oomph::EighthSphereMesh::Domain_pt
Domain * Domain_pt
Pointer to the domain.
Definition eighth_sphere_mesh.template.h:71
oomph::EighthSphereMesh::Radius
double Radius
Radius of the sphere.
Definition eighth_sphere_mesh.template.h:74
oomph::TetgenMesh
Unstructured tet mesh based on output from Tetgen: http://wias-berlin.de/software/tetgen/.
Definition tetgen_mesh.template.h:52
eighth_sphere_mesh.template.h
oomph
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition annular_domain.h:35