two_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_TWO_LAYER_SPINE_MESH_TEMPLATE_CC
31 #define OOMPH_TWO_LAYER_SPINE_MESH_TEMPLATE_CC
32 
35 
36 
37 namespace oomph
38 {
39 
40 
41 //===========================================================================
42 /// Constuctor for spine 2D mesh: Pass number of elements in x-direction,
43 /// number of elements in y-direction in bottom and top layer, respectively,
44 /// axial length and height of top and bottom layers, and pointer to
45 /// timestepper (defaults to Static timestepper).
46 ///
47 /// The mesh contains two layers of elements (of type ELEMENT;
48 /// e.g SpineElement<QCrouzeixRaviartElement<2>)
49 /// and an interfacial layer of corresponding Spine interface elements
50 /// of type INTERFACE_ELEMENT, e.g.
51 /// SpineLineFluidInterfaceElement<ELEMENT> for 2D planar
52 /// problems.
53 //===========================================================================
54 template<class ELEMENT>
56  const unsigned &nx, const unsigned &ny1, const unsigned &ny2,
57  const double &lx, const double &h1, const double &h2,
58  TimeStepper* time_stepper_pt) :
59  RectangularQuadMesh<ELEMENT >(nx,ny1+ny2,0.0,lx,0.0,h1+h2,false,false,
60  time_stepper_pt)
61 {
62  // Mesh can only be built with 2D Qelements.
63  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
64 
65  //Mesh can only be built with spine elements
66  MeshChecker::assert_geometric_element<SpineFiniteElement,ELEMENT>(2);
67 
68  // We've called the "generic" constructor for the RectangularQuadMesh
69  // which doesn't do much...
70  // Now set up the parameters that characterise the mesh geometry
71  // then call the build function
72 
73  // Number of elements in bottom and top layers
74  Ny1 = ny1;
75  Ny2 = ny2;
76 
77  // Set height of upper and lower layers
78  H1 = h1;
79  H2 = h2;
80 
81  // Now build the mesh:
82  build_two_layer_mesh(time_stepper_pt);
83 
84 }
85 
86 
87 
88 
89 //===========================================================================
90 /// Constuctor for spine 2D mesh: Pass number of elements in x-direction,
91 /// number of elements in y-direction in bottom and top layer, respectively,
92 /// axial length and height of top and bottom layers, a boolean
93 /// flag to make the mesh periodic in the x-direction, and pointer to
94 /// timestepper (defaults to Static timestepper).
95 ///
96 /// The mesh contains two layers of elements (of type ELEMENT;
97 /// e.g SpineElement<QCrouzeixRaviartElement<2>)
98 /// and an interfacial layer of corresponding Spine interface elements
99 /// of type INTERFACE_ELEMENT, e.g.
100 /// SpineLineFluidInterfaceElement<ELEMENT> for 2D planar
101 /// problems.
102 //===========================================================================
103 template<class ELEMENT>
105  const unsigned &nx, const unsigned &ny1, const unsigned &ny2,
106  const double &lx, const double &h1, const double &h2,
107  const bool& periodic_in_x, TimeStepper* time_stepper_pt) :
108  RectangularQuadMesh<ELEMENT >(nx,ny1+ny2,0.0,lx,0.0,h1+h2,periodic_in_x,
109  false,time_stepper_pt)
110 {
111  // Mesh can only be built with 2D Qelements.
112  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
113 
114  //Mesh can only be built with spine elements
115  MeshChecker::assert_geometric_element<SpineFiniteElement,ELEMENT>(2);
116 
117  // We've called the "generic" constructor for the RectangularQuadMesh
118  // which doesn't do much...
119  // Now set up the parameters that characterise the mesh geometry
120  // then call the constructor
121 
122  // Number of elements in bottom and top layers
123  Ny1 = ny1;
124  Ny2 = ny2;
125 
126  // Set height of upper and lower layers
127  H1 = h1;
128  H2 = h2;
129 
130  // Now build the mesh:
131  build_two_layer_mesh(time_stepper_pt);
132 
133 }
134 
135 
136 
137 
138 //===========================================================================
139 /// Constuctor for spine 2D mesh: Pass number of elements in x-direction,
140 /// number of elements in y-direction in bottom and top layer, respectively,
141 /// axial length and height of top and bottom layers, a boolean
142 /// flag to make the mesh periodic in the x-direction, a boolean flag to
143 /// specify whether or not to call the "build_two_layer_mesh" function,
144 /// and pointer to timestepper (defaults to Static timestepper).
145 ///
146 /// The mesh contains two layers of elements (of type ELEMENT;
147 /// e.g SpineElement<QCrouzeixRaviartElement<2>)
148 /// and an interfacial layer of corresponding Spine interface elements
149 /// of type INTERFACE_ELEMENT, e.g.
150 /// SpineLineFluidInterfaceElement<ELEMENT> for 2D planar
151 /// problems.
152 //===========================================================================
153 template<class ELEMENT>
155  const unsigned &nx, const unsigned &ny1, const unsigned &ny2,
156  const double &lx, const double &h1, const double &h2,
157  const bool& periodic_in_x, const bool& build_mesh,
158  TimeStepper* time_stepper_pt) :
159  RectangularQuadMesh<ELEMENT >(nx,ny1+ny2,0.0,lx,0.0,h1+h2,periodic_in_x,
160  false,time_stepper_pt)
161 {
162  // Mesh can only be built with 2D Qelements.
163  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
164 
165  //Mesh can only be built with spine elements
166  MeshChecker::assert_geometric_element<SpineFiniteElement,ELEMENT>(2);
167 
168  // We've called the "generic" constructor for the RectangularQuadMesh
169  // which doesn't do much...
170  // Now set up the parameters that characterise the mesh geometry
171  // then call the constructor
172 
173  // Number of elements in bottom and top layers
174  Ny1 = ny1;
175  Ny2 = ny2;
176 
177  // Set height of upper and lower layers
178  H1 = h1;
179  H2 = h2;
180 
181  // Only build the mesh here if build_mesh=true
182  // This is useful when calling this constructor from a derived class
183  // (such as Axisym2x6TwoLayerSpineMesh) where the mesh building
184  // needs to be called from *its* constructor and this constructor is
185  // only used to pass arguments to the RectangularQuadMesh constructor.
186  if(build_mesh) { build_two_layer_mesh(time_stepper_pt); }
187 
188 }
189 
190 
191 //==================================================================
192 /// \short The spacing function for the x co-ordinate, which is the
193 /// same as the default function.
194 //==================================================================
195 template<class ELEMENT>
197 x_spacing_function(unsigned xelement, unsigned xnode,
198  unsigned yelement, unsigned ynode)
199 {
200  //Calculate the values of equal increments in nodal values
201  double xstep = (this->Xmax-this->Xmin)/((this->Np-1)*this->Nx);
202  //Return the appropriate value
203  return (this->Xmin + xstep*((this->Np-1)*xelement + xnode));
204 }
205 
206 //==================================================================
207 /// \short The spacing function for the y co-ordinates, which splits
208 /// the region into two regions (1 and 2), according to the
209 /// heights H1 and H2, with Ny1 and Ny2 elements respectively.
210 template<class ELEMENT>
212 y_spacing_function(unsigned xelement, unsigned xnode,
213  unsigned yelement, unsigned ynode)
214 {
215  //Set up some spacing parameters
216  //The lower region a starts at Ymin
218  //The upper region a ends at Ymax
220  //Number of nodes per element
221  unsigned n_p = RectangularQuadMesh<ELEMENT >::Np;
222  //The lower region starts at Ymin
223  double y1init = Ymin;
224  //The upper region starts at H1 - Ymin
225  double y2init = H1 - Ymin;
226  //Calculate the space between each node in each region,
227  //Assumming uniform spacing
228  //Lower region has a length (H1-Ymin)
229  double y1step = (H1-Ymin)/((n_p-1)*Ny1);
230  //Upper region has a length (Ymax-H1)
231  double y2step = (Ymax-H1)/((n_p-1)*Ny2);
232 
233  //Now return the actual node position, it's different in the two
234  //regions, of course
235  if(yelement < Ny1)
236  {
237  return (y1init + y1step*((n_p-1)*yelement + ynode));
238  }
239  else
240  {
241  return (y2init + y2step*((n_p-1)*(yelement-Ny1) + ynode));
242  }
243 }
244 
245 //===========================================================================
246 /// Helper function that actually builds the two-layer spine mesh
247 /// based on the parameters set in the various constructors
248 //===========================================================================
249 template<class ELEMENT>
251  TimeStepper* time_stepper_pt)
252 {
253  // Build the underlying quad mesh:
255 
256  // Reset the number of boundaries
257  set_nboundary(7);
258 
259  //Set up the pointers to elements in the upper and lower fluid
260  //Calculate numbers of nodes in upper and lower regions
261  unsigned long n_lower = this->Nx*Ny1;
262  unsigned long n_upper = this->Nx*Ny2;
263  //Loop over lower elements and push back
264  Lower_layer_element_pt.reserve(n_lower);
265  for(unsigned e=0;e<n_lower;e++)
266  {
267  Lower_layer_element_pt.push_back(this->finite_element_pt(e));
268  }
269  //Loop over upper elements and push back
270  Upper_layer_element_pt.reserve(n_upper);
271  for(unsigned e=n_lower;e<(n_lower+n_upper);e++)
272  {
273  Upper_layer_element_pt.push_back(this->finite_element_pt(e));
274  }
275 
276  //Set the elements adjacent to the interface
277  Interface_lower_boundary_element_pt.resize(this->Nx);
278  Interface_upper_boundary_element_pt.resize(this->Nx);
279  {
280  unsigned count_lower=this->Nx*(this->Ny1-1);
281  unsigned count_upper= count_lower + this->Nx;
282  for(unsigned e=0;e<this->Nx;e++)
283  {
284  Interface_lower_boundary_element_pt[e] =
285  this->finite_element_pt(count_lower); ++count_lower;
286  Interface_upper_boundary_element_pt[e] =
287  this->finite_element_pt(count_upper); ++count_upper;
288  }
289  }
290 
291  //Relabel the boundary nodes
292  //Storage for the boundary coordinates that will be transferred directly
293  Vector<double> b_coord;
294  {
295  //Store the interface level
296  const double y_interface = H1 - this->Ymin;
297 
298  //Nodes on boundary 3 are now on boundaries 4 and 5
299  unsigned n_boundary_node = this->nboundary_node(3);
300  //Loop over the nodes remove them from boundary 3
301  //and add them to boundary 4 or 5 depending on their y coordinate
302  for(unsigned n=0;n<n_boundary_node;n++)
303  {
304  //Cache pointer to the node
305  Node* const nod_pt = this->boundary_node_pt(3,n);
306  //Get the boundary coordinates if set
307  if(boundary_coordinate_exists(3))
308  {
309  b_coord.resize(nod_pt->ncoordinates_on_boundary(3));
310  nod_pt->get_coordinates_on_boundary(3,b_coord);
311  //Indicate that the boundary coordinates are to be set on the
312  //new boundaries
313  this->Boundary_coordinate_exists[4]=true;
314  this->Boundary_coordinate_exists[5]=true;
315  }
316 
317  //Now remove the node from the boundary
318  nod_pt->remove_from_boundary(3);
319 
320  //Find the height of the node
321  double y = nod_pt->x(1);
322  //If it's above or equal to the interface, it's on boundary 4
323  if(y >= y_interface)
324  {
325  this->add_boundary_node(4,nod_pt);
326  //Add the boundary coordinate if it has been set up
327  if(this->Boundary_coordinate_exists[4])
328  {
329  nod_pt->set_coordinates_on_boundary(4,b_coord);
330  }
331  }
332  //otherwise it's on boundary 5
333  if(y <= y_interface)
334  {
335  this->add_boundary_node(5,nod_pt);
336  //Add the boundary coordinate if it has been set up
337  if(this->Boundary_coordinate_exists[5])
338  {
339  nod_pt->set_coordinates_on_boundary(5,b_coord);
340  }
341  }
342  }
343 
344  //Now clear the boundary node information from boundary 3
345  this->Boundary_node_pt[3].clear();
346 
347  //Relabel the nodes on boundary 2 to be on boundary 3
348  n_boundary_node = this->nboundary_node(2);
349  //Loop over the nodes remove them from boundary 2
350  //and add them to boundary 3
351  for(unsigned n=0;n<n_boundary_node;n++)
352  {
353  //Cache pointer to the node
354  Node* const nod_pt = this->boundary_node_pt(2,n);
355  //Get the boundary coordinates if set
356  if(this->boundary_coordinate_exists(2))
357  {
358  b_coord.resize(nod_pt->ncoordinates_on_boundary(2));
359  nod_pt->get_coordinates_on_boundary(2,b_coord);
360  this->Boundary_coordinate_exists[3]=true;
361  }
362 
363  //Now remove the node from the boundary 2
364  nod_pt->remove_from_boundary(2);
365  //and add to boundary 3
366  this->add_boundary_node(3,nod_pt);
367  if(this->Boundary_coordinate_exists[3])
368  {
369  nod_pt->set_coordinates_on_boundary(3,b_coord);
370  }
371  }
372 
373  //Clear the information from boundary 2
374  this->Boundary_node_pt[2].clear();
375 
376  std::list<Node*> nodes_to_be_removed;
377 
378  //Nodes on boundary 1 are now on boundaries 1 and 2
379  n_boundary_node = this->nboundary_node(1);
380  //Loop over the nodes remove some of them from boundary 1
381  for(unsigned n=0;n<n_boundary_node;n++)
382  {
383  //Cache pointer to the node
384  Node* const nod_pt = this->boundary_node_pt(1,n);
385 
386  //Find the height of the node
387  double y = nod_pt->x(1);
388  //If it's above or equal to the interface it's on boundary 2
389  if(y >= y_interface)
390  {
391  //Get the boundary coordinates if set
392  if(this->boundary_coordinate_exists(1))
393  {
394  b_coord.resize(nod_pt->ncoordinates_on_boundary(1));
395  nod_pt->get_coordinates_on_boundary(1,b_coord);
396  this->Boundary_coordinate_exists[2]=true;
397  }
398 
399  //Now remove the node from the boundary 1 if above interace
400  if(y > y_interface)
401  {
402  nodes_to_be_removed.push_back(nod_pt);
403  }
404  //Always add to boundary 2
405  this->add_boundary_node(2,nod_pt);
406  //Add the boundary coordinate if it has been set up
407  if(this->Boundary_coordinate_exists[2])
408  {
409  nod_pt->set_coordinates_on_boundary(2,b_coord);
410  }
411  }
412  }
413 
414  //Loop over the nodes that are to be removed from boundary 1 and remove
415  //them
416  for(std::list<Node*>::iterator it=nodes_to_be_removed.begin();
417  it!=nodes_to_be_removed.end();++it)
418  {
419  this->remove_boundary_node(1,*it);
420  }
421  nodes_to_be_removed.clear();
422  }
423  //Allocate memory for the spines and fractions along spines
424  //---------------------------------------------------------
425 
426  //Read out number of linear points in the element
427  unsigned n_p = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
428 
429 
430  //Add the nodes on the interface to the boundary 6
431  //Storage for boundary coordinates (x-coordinate)
432  b_coord.resize(1);
433  this->Boundary_coordinate_exists[6];
434  //Starting index of the nodes
435  unsigned n_start=0;
436  for(unsigned e=0;e<this->Nx;e++)
437  {
438  //If we are past the
439  if(e>0) {n_start=1;}
440  //Get the layer of elements just above the interface
441  FiniteElement* el_pt = this->finite_element_pt(this->Nx*this->Ny1+e);
442  //The first n_p nodes lie on the boundary
443  for(unsigned n=n_start;n<n_p;n++)
444  {
445  Node* nod_pt = el_pt->node_pt(n);
446  this->convert_to_boundary_node(nod_pt);
447  this->add_boundary_node(6,nod_pt);
448  b_coord[0] = nod_pt->x(0);
449  nod_pt->set_coordinates_on_boundary(6,b_coord);
450  }
451  }
452 
453  //Allocate store for the spines:
454  if (this->Xperiodic)
455  {
456  Spine_pt.reserve((n_p-1)*this->Nx);
457  }
458  else
459  {
460  Spine_pt.reserve((n_p-1)*this->Nx+1);
461  }
462 
463 
464  //FIRST SPINE
465  //-----------
466 
467  //Element 0
468  //Node 0
469  //Assign the new spine of height of the lower layer
470  Spine* new_spine_pt=new Spine(H1);
471  Spine_pt.push_back(new_spine_pt);
472 
473  // Get pointer to node
474  SpineNode* nod_pt=element_node_pt(0,0);
475 
476  //Set the pointer to the spine
477  nod_pt->spine_pt() = new_spine_pt;
478  //Set the fraction
479  nod_pt->fraction() = 0.0;
480  // Pointer to the mesh that implements the update fct
481  nod_pt->spine_mesh_pt() = this;
482  // ID of update function within this mesh: 0 = lower; 1 = upper
483  nod_pt->node_update_fct_id() = 0;
484 
485  //Loop vertically along the spine
486  //Loop over the elements in fluid 1
487  for(unsigned long i=0;i<Ny1;i++)
488  {
489  //Loop over the vertical nodes, apart from the first
490  for(unsigned l1=1;l1<n_p;l1++)
491  {
492  // Get pointer to node
493  SpineNode* nod_pt=element_node_pt(i*this->Nx,l1*n_p);
494  //Set the pointer to the spine
495  nod_pt->spine_pt() = new_spine_pt;
496  //Set the fraction
497  nod_pt->fraction() = (nod_pt->x(1)-this->Ymin)/(H1);
498  // Pointer to the mesh that implements the update fct
499  nod_pt->spine_mesh_pt() = this;
500  // ID of update function within this mesh: 0 = lower; 1 = upper
501  nod_pt->node_update_fct_id() = 0;
502  }
503  }
504 
505  //Loop over the elements in fluid 2
506  for(unsigned long i=0;i<Ny2;i++)
507  {
508  //Loop over vertical nodes, apart from the first
509  for(unsigned l1=1;l1<n_p;l1++)
510  {
511  // Get pointer to node
512  SpineNode* nod_pt=element_node_pt((Ny1+i)*this->Nx,l1*n_p);
513 
514  //Set the pointer to the spine
515  nod_pt->spine_pt() = new_spine_pt;
516  //Set the fraction
517  nod_pt->fraction() =(nod_pt->x(1)-(this->Ymin+H1))/(H2);
518  // Pointer to the mesh that implements the update fct
519  nod_pt->spine_mesh_pt() = this;
520  // ID of update function within this mesh: 0 = lower; 1 = upper
521  nod_pt->node_update_fct_id() = 1;
522  }
523  }
524 
525 
526  //LOOP OVER OTHER SPINES
527  //----------------------
528 
529  //Now loop over the elements horizontally
530  for(unsigned long j=0;j<this->Nx;j++)
531  {
532  //Loop over the nodes in the elements horizontally, ignoring
533  //the first column
534 
535  // Last spine needs special treatment in x-periodic meshes:
536  unsigned n_pmax=n_p;
537  if ((this->Xperiodic)&&(j==this->Nx-1)) n_pmax=n_p-1;
538 
539  for(unsigned l2=1;l2<n_pmax;l2++)
540  {
541  //Assign the new spine with length the height of the lower layer
542  new_spine_pt=new Spine(H1);
543  Spine_pt.push_back(new_spine_pt);
544 
545  // Get pointer to node
546  SpineNode* nod_pt=element_node_pt(j,l2);
547 
548  //Set the pointer to the spine
549  nod_pt->spine_pt() = new_spine_pt;
550  //Set the fraction
551  nod_pt->fraction() = 0.0;
552  // Pointer to the mesh that implements the update fct
553  nod_pt->spine_mesh_pt() = this;
554  // ID of update function within this mesh: 0 = lower; 1 = upper
555  nod_pt->node_update_fct_id() = 0;
556 
557  //Loop vertically along the spine
558  //Loop over the elements in fluid 1
559  for(unsigned long i=0;i<Ny1;i++)
560  {
561  //Loop over the vertical nodes, apart from the first
562  for(unsigned l1=1;l1<n_p;l1++)
563  {
564  // Get pointer to node
565  SpineNode* nod_pt=element_node_pt(i*this->Nx+j,l1*n_p+l2);
566  //Set the pointer to the spine
567  nod_pt->spine_pt() = new_spine_pt;
568  //Set the fraction
569  nod_pt->fraction() = (nod_pt->x(1)-this->Ymin)/H1;
570  // Pointer to the mesh that implements the update fct
571  nod_pt->spine_mesh_pt() = this;
572  // ID of update function within this mesh: 0 = lower; 1 = upper
573  nod_pt->node_update_fct_id() = 0;
574  }
575  }
576 
577  //Loop over the elements in fluid 2
578  for(unsigned long i=0;i<Ny2;i++)
579  {
580  //Loop over vertical nodes, apart from the first
581  for(unsigned l1=1;l1<n_p;l1++)
582  {
583  // Get pointer to node
584  SpineNode* nod_pt=element_node_pt((Ny1+i)*this->Nx+j,l1*n_p+l2);
585 
586  //Set the pointer to the spine
587  nod_pt->spine_pt() = new_spine_pt;
588  //Set the fraction
589  nod_pt->fraction() = (nod_pt->x(1)-(this->Ymin+H1))/H2;
590  // Pointer to the mesh that implements the update fct
591  nod_pt->spine_mesh_pt() = this;
592  // ID of update function within this mesh: 0 = lower; 1 = upper
593  nod_pt->node_update_fct_id() = 1;
594  }
595  }
596  }
597  }
598 
599 
600  // Last spine needs special treatment for periodic meshes
601  // because it's the same as the first one...
602  if (this->Xperiodic)
603  {
604  // Last spine is the same as first one...
605  Spine* final_spine_pt=Spine_pt[0];
606 
607  // Get pointer to node
608  SpineNode* nod_pt=element_node_pt((this->Nx-1),(n_p-1));
609 
610  //Set the pointer to the spine
611  nod_pt->spine_pt() = final_spine_pt;
612  //Set the fraction to be the same as for the nodes on the first row
613  nod_pt->fraction() = element_node_pt(0,0)->fraction();
614  // Pointer to the mesh that implements the update fct
615  nod_pt->spine_mesh_pt() = element_node_pt(0,0)->spine_mesh_pt();
616  // ID of update function within this mesh: 0 = lower; 1 = upper
617  nod_pt->node_update_fct_id() = element_node_pt(0,0)->node_update_fct_id();
618 
619  //Now loop vertically along the spine
620  for(unsigned i=0;i<(Ny1+Ny2);i++)
621  {
622  //Loop over the vertical nodes, apart from the first
623  for(unsigned l1=1;l1<n_p;l1++)
624  {
625  // Get pointer to node
626  SpineNode* nod_pt =
627  element_node_pt(i*this->Nx+(this->Nx-1),l1*n_p+(n_p-1));
628 
629  //Set the pointer to the spine
630  nod_pt->spine_pt() = final_spine_pt;
631  //Set the fraction to be the same as in first row
632  nod_pt->fraction() = element_node_pt(i*this->Nx,l1*n_p)->fraction();
633  // ID of update function within this mesh: 0 = lower; 1 = upper
634  nod_pt->node_update_fct_id() =
635  element_node_pt(i*this->Nx,l1*n_p)->node_update_fct_id();
636  // Pointer to the mesh that implements the update fct
637  nod_pt->spine_mesh_pt() = element_node_pt(i*this->Nx,l1*n_p)
638  ->spine_mesh_pt();
639  }
640  }
641  }
642 
643 
644 
645  //Assign the 1D Line elements
646  //---------------------------
647 
648  //Get the present number of elements
649  /*unsigned long element_count = Element_pt.size();
650 
651  //Loop over the horizontal elements
652  for(unsigned i=0;i<this->Nx;i++)
653  {
654  //Construct a new 1D line element on the face on which the local
655  //coordinate 1 is fixed at its max. value (1) -- Face 2
656  FiniteElement *interface_element_element_pt =
657  new INTERFACE_ELEMENT(finite_element_pt(this->Nx*(Ny1-1)+i),2);
658 
659  //Push it back onto the stack
660  Element_pt.push_back(interface_element_element_pt);
661 
662  //Push it back onto the stack of interface elements
663  Interface_element_pt.push_back(interface_element_element_pt);
664 
665  element_count++;
666  }*/
667 
668 //Setup the boundary information
669 this->setup_boundary_element_info();
670 
671 }
672 
673 
674 //======================================================================
675 /// Reorder the elements, so we loop over them vertically first
676 /// (advantageous in "wide" domains if a frontal solver is used).
677 //======================================================================
678 /*template <class ELEMENT, >
679 void TwoLayerSpineMesh<ELEMENT,INTERFACE_ELEMENT>::element_reorder()
680 {
681 
682  unsigned Nx = this->Nx;
683  //Find out how many elements there are
684  unsigned long Nelement = nelement();
685  //Find out how many fluid elements there are
686  unsigned long Nfluid = Nx*(Ny1+Ny2);
687  //Create a dummy array of elements
688  Vector<FiniteElement *> dummy;
689 
690  //Loop over the elements in horizontal order
691  for(unsigned long j=0;j<Nx;j++)
692  {
693  //Loop over the elements in lower layer vertically
694  for(unsigned long i=0;i<Ny1;i++)
695  {
696  //Push back onto the new stack
697  dummy.push_back(finite_element_pt(Nx*i + j));
698  }
699 
700  //Push back the line element onto the stack
701  dummy.push_back(finite_element_pt(Nfluid+j));
702 
703  //Loop over the elements in upper layer vertically
704  for(unsigned long i=Ny1;i<(Ny1+Ny2);i++)
705  {
706  //Push back onto the new stack
707  dummy.push_back(finite_element_pt(Nx*i + j));
708  }
709  }
710 
711  //Now copy the reordered elements into the element_pt
712  for(unsigned long e=0;e<Nelement;e++)
713  {
714  Element_pt[e] = dummy[e];
715  }
716 
717  }*/
718 
719 }
720 #endif
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Set the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interf...
Definition: nodes.cc:2301
virtual void build_two_layer_mesh(TimeStepper *time_stepper_pt)
Helper function to actually build the two-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
A general Finite Element class.
Definition: elements.h:1271
double y_spacing_function(unsigned xelement, unsigned xnode, unsigned yelement, unsigned ynode)
The spacing function for the y co-ordinates with three regions in each fluid.
unsigned & node_update_fct_id()
Access function to ID of node update function (within specific mesh)
Definition: spines.h:353
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
e
Definition: cfortran.h:575
virtual void get_coordinates_on_boundary(const unsigned &b, const unsigned &k, Vector< double > &boundary_zeta)
Return the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interface provides run-time error checking.
Definition: nodes.cc:2287
virtual unsigned ncoordinates_on_boundary(const unsigned &b)
Get the number of boundary coordinates on mesh boundary b. Broken virtual interface provides run-time...
Definition: nodes.cc:2271
virtual void remove_from_boundary(const unsigned &b)
Broken interface for removing the node from the mesh boundary b Here to provide error reporting...
Definition: nodes.cc:2258
TwoLayerSpineMesh(const unsigned &nx, const unsigned &ny1, const unsigned &ny2, const double &lx, const double &h1, const double &h2, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in x-direction, number of elements in y-direction in bottom and ...
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
unsigned Ny2
Number of elements in upper layer.
unsigned Ny1
Number of elements in lower layer.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double H1
Height of the lower layer.
double H2
Height of the upper layer.
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
double x_spacing_function(unsigned xelement, unsigned xnode, unsigned yelement, unsigned ynode)
The spacing function for the x co-ordinates with two regions.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219