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