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