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-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 #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 
32 namespace 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
eighth_sphere_mesh.template.h
oomph
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: annular_domain.h:35