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