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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 
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>
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
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];
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];
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];
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];
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
367 
368 }
369 
370 }
371 #endif
void setup_boundary_element_info()
Definition: brick_mesh.h:223
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:194
Domain * Domain_pt
Pointer to the domain.
EighthSphereMesh(const double &radius, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass radius and timestepper; defaults to static default timestepper. ...
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:246
cstr elem_len * i
Definition: cfortran.h:607
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
double Radius
Radius of the sphere.
MacroElement * macro_element_pt(const unsigned &i)
Access to i-th macro element.
Definition: domain.h:100
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2139
virtual void local_fraction_of_node(const unsigned &j, Vector< double > &s_fraction)
Get the local fraction of the node j in the element A dumb, but correct default implementation is pro...
Definition: elements.cc:3129
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
static char t char * s
Definition: cfortran.h:572
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
void macro_map(const Vector< double > &s, Vector< double > &r)
The mapping from local to global coordinates at the current time : r(s)
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:197
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:456
Eighth sphere as domain. Domain is parametrised by four macro elements.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219