single_layer_spine_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_SINGLE_LAYER_SPINE_MESH_TEMPLATE_CC
31 #define OOMPH_SINGLE_LAYER_SPINE_MESH_TEMPLATE_CC
32 
35 
36 
37 namespace oomph
38 {
39 
40 
41 
42 //===========================================================================
43 /// Constructor for spine 2D mesh: Pass number of elements in x-direction,
44 /// number of elements in y-direction, axial length and height of layer,
45 /// and pointer to timestepper (defaults to Static timestepper).
46 ///
47 /// The mesh must be called with spinified elements and it
48 /// constructs the spines and contains the information on how to update
49 /// the nodal positions within the mesh as a function of the spine lengths.
50 /// Equations that determine the spine heights (even if they are pinned)
51 /// must be specified externally or else there will be problems.
52 
53 //===========================================================================
54 template<class ELEMENT>
56  const unsigned &nx, const unsigned &ny,
57  const double &lx, const double &h, TimeStepper* time_stepper_pt) :
58  RectangularQuadMesh<ELEMENT >(nx,ny,0.0,lx,0.0,h,false,false,
59  time_stepper_pt)
60 {
61  // Mesh can only be built with 2D Qelements.
62  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
63 
64  //Mesh can only be built with spine elements
65  MeshChecker::assert_geometric_element<SpineFiniteElement,ELEMENT>(2);
66 
67  // We've called the "generic" constructor for the RectangularQuadMesh
68  // which doesn't do much...
69 
70  // Now build the mesh:
71  build_single_layer_mesh(time_stepper_pt);
72 }
73 
74 
75 
76 //===========================================================================
77 /// Constuctor for spine 2D mesh: Pass number of elements in x-direction,
78 /// number of elements in y-direction, axial length and height of layer,
79 /// a boolean flag to make the mesh periodic in the x-direction,
80 /// and pointer to timestepper (defaults to Static timestepper).
81 ///
82 /// The mesh must be called with spinified elements and it
83 /// constructs the spines and contains the information on how to update
84 /// the nodal positions within the mesh as a function of the spine lengths.
85 /// Equations that determine the spine heights (even if they are pinned)
86 /// must be specified externally or else there will be problems.
87 //===========================================================================
88 template<class ELEMENT>
90  const unsigned &nx, const unsigned &ny,
91  const double &lx, const double &h, const bool& periodic_in_x,
92  TimeStepper* time_stepper_pt) :
93  RectangularQuadMesh<ELEMENT >(nx,ny,0.0,lx,0.0,h,periodic_in_x,false,
94  time_stepper_pt)
95 {
96  // Mesh can only be built with 2D Qelements.
97  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
98 
99  //Mesh can only be built with spine elements
100  MeshChecker::assert_geometric_element<SpineFiniteElement,ELEMENT>(2);
101 
102 
103  // We've called the "generic" constructor for the RectangularQuadMesh
104  // which doesn't do much...
105 
106  // Now build the mesh:
107  build_single_layer_mesh(time_stepper_pt);
108 }
109 
110 
111 
112 
113 
114 //===========================================================================
115 /// Helper function that actually builds the single-layer spine mesh
116 /// based on the parameters set in the various constructors
117 //===========================================================================
118 template<class ELEMENT>
120  TimeStepper* time_stepper_pt)
121 {
122 
123  // Mesh can only be built with 2D Qelements.
124  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
125 
126  // Build the underlying quad mesh:
128 
129  //Read out the number of elements in the x-direction
130  unsigned n_x = this->Nx;
131  unsigned n_y = this->Ny;
132 
133  //Allocate memory for the spines and fractions along spines
134  //---------------------------------------------------------
135 
136  //Read out number of linear points in the element
137  unsigned n_p = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
138 
139  //Allocate store for the spines:
140  if (this->Xperiodic)
141  {
142  Spine_pt.reserve((n_p-1)*n_x);
143  }
144  else
145  {
146  Spine_pt.reserve((n_p-1)*n_x+1);
147  }
148 
149 
150  //FIRST SPINE
151  //-----------
152 
153  //Element 0
154  //Node 0
155  //Assign the new spine with unit length
156  Spine* new_spine_pt=new Spine(1.0);
157  Spine_pt.push_back(new_spine_pt);
158 
159 
160  // Get pointer to node
161  SpineNode* nod_pt=element_node_pt(0,0);
162  //Set the pointer to the spine
163  nod_pt->spine_pt() = new_spine_pt;
164  //Set the fraction
165  nod_pt->fraction() = 0.0;
166  // Pointer to the mesh that implements the update fct
167  nod_pt->spine_mesh_pt() = this;
168 
169  //Loop vertically along the spine
170  //Loop over the elements
171  for(unsigned long i=0;i<n_y;i++)
172  {
173  //Loop over the vertical nodes, apart from the first
174  for(unsigned l1=1;l1<n_p;l1++)
175  {
176  // Get pointer to node
177  SpineNode* nod_pt=element_node_pt(i*n_x,l1*n_p);
178  //Set the pointer to the spine
179  nod_pt->spine_pt() = new_spine_pt;
180  //Set the fraction
181  nod_pt->fraction()=(double(i)+double(l1)/double(n_p-1))/double(n_y);
182  // Pointer to the mesh that implements the update fct
183  nod_pt->spine_mesh_pt() = this;
184  }
185  }
186 
187 
188  //LOOP OVER OTHER SPINES
189  //----------------------
190 
191  //Now loop over the elements horizontally
192  for(unsigned long j=0;j<n_x;j++)
193  {
194  //Loop over the nodes in the elements horizontally, ignoring
195  //the first column
196 
197  // Last spine needs special treatment in x-periodic meshes:
198  unsigned n_pmax=n_p;
199  if ((this->Xperiodic)&&(j==n_x-1)) n_pmax=n_p-1;
200 
201  for(unsigned l2=1;l2<n_pmax;l2++)
202  {
203  //Assign the new spine with unit height
204  new_spine_pt=new Spine(1.0);
205  Spine_pt.push_back(new_spine_pt);
206 
207  // Get the node
208  SpineNode* nod_pt=element_node_pt(j,l2);
209  //Set the pointer to spine
210  nod_pt->spine_pt() = new_spine_pt;
211  //Set the fraction
212  nod_pt->fraction() = 0.0;
213  // Pointer to the mesh that implements the update fct
214  nod_pt->spine_mesh_pt() = this;
215 
216  //Loop vertically along the spine
217  //Loop over the elements
218  for(unsigned long i=0;i<n_y;i++)
219  {
220  //Loop over the vertical nodes, apart from the first
221  for(unsigned l1=1;l1<n_p;l1++)
222  {
223  // Get the node
224  SpineNode* nod_pt=element_node_pt(i*n_x+j,l1*n_p+l2);
225  //Set the pointer to the spine
226  nod_pt->spine_pt() = new_spine_pt;
227  //Set the fraction
228  nod_pt->fraction()=(double(i)+double(l1)/double(n_p-1))/double(n_y);
229  // Pointer to the mesh that implements the update fct
230  nod_pt->spine_mesh_pt() = this;
231  }
232  }
233  }
234  }
235 
236 
237  // Last spine needs special treatment for periodic meshes
238  // because it's the same as the first one...
239  if (this->Xperiodic)
240  {
241 
242  // Last spine is the same as first one...
243  Spine* final_spine_pt=Spine_pt[0];
244 
245  // Get the node
246  SpineNode* nod_pt=element_node_pt((n_x-1),(n_p-1));
247 
248  //Set the pointer for the first node
249  nod_pt->spine_pt() = final_spine_pt;
250  //Set the fraction to be the same as for the nodes on the first row
251  nod_pt->fraction() = element_node_pt(0,0)->fraction();
252  // Pointer to the mesh that implements the update fct
253  nod_pt->spine_mesh_pt() = element_node_pt(0,0)->spine_mesh_pt();
254 
255  //Now loop vertically along the spine
256  for(unsigned i=0;i<n_y;i++)
257  {
258  //Loop over the vertical nodes, apart from the first
259  for(unsigned l1=1;l1<n_p;l1++)
260  {
261  // Get the node
262  SpineNode* nod_pt=element_node_pt(i*n_x+(n_x-1),l1*n_p+(n_p-1));
263 
264  //Set the pointer to the spine
265  nod_pt->spine_pt() = final_spine_pt;
266  //Set the fraction to be the same as in first row
267  nod_pt->fraction() = element_node_pt(i*n_x,l1*n_p)->fraction();
268  // Pointer to the mesh that implements the update fct
269  nod_pt->spine_mesh_pt() =
270  element_node_pt(i*n_x,l1*n_p)->spine_mesh_pt();
271  }
272  }
273 
274  }
275 }
276 
277 }
278 #endif
SingleLayerSpineMesh(const unsigned &nx, const unsigned &ny, const double &lx, const double &h, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in x-direction, number of elements in y-direction, axial length, height of layer, and pointer to timestepper (defaults to Steady timestepper)
virtual void build_single_layer_mesh(TimeStepper *time_stepper_pt)
Helper function to actually build the single-layer spine mesh (called from various constructors) ...
cstr elem_len * i
Definition: cfortran.h:607
double & fraction()
Set reference to fraction along spine.
Definition: spines.h:350
void build_mesh(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh construction function: contains all the hard work.
Spine *& spine_pt()
Access function to spine.
Definition: spines.h:347
SpineMesh *& spine_mesh_pt()
Access function to Pointer to SpineMesh that this node is a part of and which implements the node upd...
Definition: spines.h:357
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219