bretherton_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_BRETHERTON_SPINE_MESH_TEMPLATE_CC
27 #define OOMPH_BRETHERTON_SPINE_MESH_TEMPLATE_CC
28 
30 #include "../generic/mesh_as_geometric_object.h"
31 #include "../generic/face_element_as_geometric_object.h"
32 
35 
36 
37 namespace oomph
38 {
39  //======================================================================
40  /// Constructor. Arguments:
41  /// - nx1: Number of elements along wall in deposited film region
42  /// - nx2: Number of elements along wall in horizontal transition region
43  /// - nx3: Number of elements along wall in channel region
44  /// - nhalf: Number of elements in vertical transition region (there are
45  /// twice as many elements across the channel)
46  /// - nh: Number of elements across the deposited film
47  /// - h: Thickness of deposited film
48  /// - zeta0: Start coordinate on wall
49  /// - zeta1: Wall coordinate of start of transition region
50  /// - zeta2: Wall coordinate of end of liquid filled region (inflow)
51  /// - lower_wall_pt: Pointer to geometric object that represents the lower
52  /// wall
53  /// - upper_wall_pt: Pointer to geometric object that represents the upper
54  /// wall
55  /// - time_stepper_pt: Pointer to timestepper; defaults to Static
56  //======================================================================
57  template<class ELEMENT, class INTERFACE_ELEMENT>
59  const unsigned& nx1,
60  const unsigned& nx2,
61  const unsigned& nx3,
62  const unsigned& nh,
63  const unsigned& nhalf,
64  const double& h,
65  GeomObject* lower_wall_pt,
66  GeomObject* upper_wall_pt,
67  const double& zeta_start,
68  const double& zeta_transition_start,
69  const double& zeta_transition_end,
70  const double& zeta_end,
71  TimeStepper* time_stepper_pt)
72  : SingleLayerSpineMesh<ELEMENT>(
73  2 * (nx1 + nx2 + nhalf), nh, 1.0, h, time_stepper_pt),
74  Nx1(nx1),
75  Nx2(nx2),
76  Nx3(nx3),
77  Nhalf(nhalf),
78  Nh(nh),
79  H(h),
80  Upper_wall_pt(upper_wall_pt),
81  Lower_wall_pt(lower_wall_pt),
82  Zeta_start(zeta_start),
83  Zeta_end(zeta_end),
84  Zeta_transition_start(zeta_transition_start),
85  Zeta_transition_end(zeta_transition_end),
86  Spine_centre_fraction_pt(&Default_spine_centre_fraction),
87  Default_spine_centre_fraction(0.5)
88  {
89  // Add the created elements to the bulk domain
90  unsigned n_bulk = this->nelement();
91  Bulk_element_pt.resize(n_bulk);
92  for (unsigned e = 0; e < n_bulk; e++)
93  {
94  Bulk_element_pt[e] = this->finite_element_pt(e);
95  }
96 
97  // Mesh can only be built with 2D Qelements.
98  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
99 
100  // Mesh can only be built with spine elements
101  MeshChecker::assert_geometric_element<SpineFiniteElement, ELEMENT>(2);
102 
103  // Initialise start of transition region to zero
104  // Zeta_transition_start = 0.0;
105 
106  // Length of deposited film region
107  double llayer = Zeta_transition_start - Zeta_start;
108 
109  // Length of transition region
111 
112  // Work out radius of circular cap from lower and upper wall
113  Vector<double> r_wall_lo(2), r_wall_up(2);
114  Vector<double> zeta(1), s_lo(1), s_up(1);
115  GeomObject *lower_sub_geom_object_pt = 0, *upper_sub_geom_object_pt = 0;
116 
117  GeomObject* lower_transition_geom_object_pt = 0;
118  GeomObject* upper_transition_geom_object_pt = 0;
119  Vector<double> s_transition_lo(1), s_transition_up(1);
120  Vector<double> spine_centre(2);
121 
122  // Find position of lower and upper walls at start of transition region
123  zeta[0] = Zeta_transition_start;
124  Lower_wall_pt->position(zeta, r_wall_lo);
125  Upper_wall_pt->position(zeta, r_wall_up);
126  // Radius is the difference between the film thickness and the distance
127  // to the lower wall
128  double radius = -r_wall_lo[1] - H;
129 
130  // Check to non-symmetric mesh
131  if (std::fabs(r_wall_lo[1] + r_wall_up[1]) > 1.0e-4)
132  {
133  oomph_info << "\n\n=================================================== "
134  << std::endl;
135  oomph_info << "Warning: " << std::endl;
136  oomph_info << "-------- " << std::endl;
137  oomph_info << " " << std::endl;
138  oomph_info << "Upper and lower walls are not symmetric at zeta=0"
139  << std::endl;
140  oomph_info << "Your initial mesh will look very odd." << std::endl;
141  oomph_info << "y-coordinates of walls at zeta=0 are: " << r_wall_lo[1]
142  << " and " << r_wall_up[1] << std::endl
143  << std::endl;
144  oomph_info << "===================================================\n\n "
145  << std::endl;
146  }
147 
148 
149  // Add the interface elements
150  // Loop over the horizontal elements
151  {
152  unsigned n_x = 2 * (nx1 + nx2 + nhalf);
153  unsigned n_y = nh;
154  for (unsigned i = 0; i < n_x; i++)
155  {
156  // Construct a new 1D line element on the face on which the local
157  // coordinate 1 is fixed at its max. value (1) --- This is face 2
158  FiniteElement* interface_element_element_pt = new INTERFACE_ELEMENT(
159  this->finite_element_pt(n_x * (n_y - 1) + i), 2);
160 
161  // Push it back onto the stack
162  this->Element_pt.push_back(interface_element_element_pt);
163 
164  // Push it back onto the stack of interface elements
165  this->Interface_element_pt.push_back(interface_element_element_pt);
166  }
167  }
168 
169  // Reorder elements: Vertical stacks of elements, each topped by
170  // their interface element -- this is (currently) identical to the
171  // version in the SingleLayerSpineMesh but it's important
172  // that element reordering is maintained in exactly this form
173  // so to be on the safe side, we move the function in here.
175 
176  // Store pointer to control element
177  Control_element_pt = dynamic_cast<ELEMENT*>(
178  this->element_pt((nx1 + nx2 + nhalf) * (nh + 1) - 2));
179 
180  // Temporary storage for boundary lookup scheme
181  Vector<std::set<Node*>> set_boundary_node_pt(6);
182 
183  // Boundary 1 -> 3; 2 -> 4; 3 -> 5
184  for (unsigned ibound = 1; ibound <= 3; ibound++)
185  {
186  unsigned numnod = this->nboundary_node(ibound);
187  for (unsigned j = 0; j < numnod; j++)
188  {
189  set_boundary_node_pt[ibound + 2].insert(
190  this->boundary_node_pt(ibound, j));
191  }
192  }
193 
194  // Get number of nodes per element
195  unsigned nnod = this->finite_element_pt(0)->nnode();
196 
197  // Get number of nodes along element edge
198  unsigned np = this->finite_element_pt(0)->nnode_1d();
199 
200  // Initialise number of elements in previous regions:
201  unsigned n_prev_elements = 0;
202 
203  // We've now built the straight single-layer mesh and need to change the
204  // the update functions for all nodes so that the domain
205  // deforms into the Bretherton shape
206 
207  // Loop over elements in lower deposited film region
208  // -------------------------------------------------
209  {
210  // Increments in wall coordinate
211  double dzeta_el = llayer / double(nx1);
212  double dzeta_node = llayer / double(nx1 * (np - 1));
213 
214  // Loop over elements horizontally
215  for (unsigned i = 0; i < nx1; i++)
216  {
217  // Start of wall coordinate
218  double zeta_lo = Zeta_start + double(i) * dzeta_el;
219 
220  // Loop over elements vertically
221  for (unsigned j = 0; j < nh; j++)
222  {
223  // Work out element number in overall mesh
224  unsigned e = n_prev_elements + i * (nh + 1) + j;
225 
226  // Get pointer to element
227  FiniteElement* el_pt = this->finite_element_pt(e);
228 
229  // Loop over its nodes
230  for (unsigned i0 = 0; i0 < np; i0++)
231  {
232  for (unsigned i1 = 0; i1 < np; i1++)
233  {
234  // Node number:
235  unsigned n = i0 * np + i1;
236 
237  // Get spine node
238  SpineNode* nod_pt = dynamic_cast<SpineNode*>(el_pt->node_pt(n));
239 
240  // Set update fct id
241  nod_pt->node_update_fct_id() = 0;
242 
243  // Provide spine with additional storage for wall coordinate
244  // and wall geom object:
245  if (i0 == 0)
246  {
247  // Get the Lagrangian coordinate in the Lower Wall
248  zeta[0] = zeta_lo + double(i1) * dzeta_node;
249  // Get the geometric object and local coordinate
250  Lower_wall_pt->locate_zeta(
251  zeta, lower_sub_geom_object_pt, s_lo);
252 
253  // The local coordinate is a geometric parameter
254  // This needs to be set (rather than added) because the
255  // same spine may be visited more than once
256  Vector<double> parameters(1, s_lo[0]);
257  nod_pt->spine_pt()->set_geom_parameter(parameters);
258 
259  // Adjust spine height
260  nod_pt->spine_pt()->height() = H;
261 
262  // The sub geom object is one (and only) geom object
263  // for spine:
264  Vector<GeomObject*> geom_object_pt(1);
265  geom_object_pt[0] = lower_sub_geom_object_pt;
266 
267  // Pass geom object(s) to spine
268  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
269 
270  // Push the node back onto boundaries
271  if (j == 0) set_boundary_node_pt[0].insert(nod_pt);
272  }
273  }
274  }
275  }
276  }
277 
278  // Increment number of previous elements
279  n_prev_elements += nx1 * (nh + 1);
280  }
281 
282  {
283  // Calculate the centre for the spine nodes in the transition region
284  zeta[0] = Zeta_transition_start;
285  // Get the geometric objects on the walls at the start of the transition
286  // region
287  Lower_wall_pt->locate_zeta(
288  zeta, lower_transition_geom_object_pt, s_transition_lo);
289  Upper_wall_pt->locate_zeta(
290  zeta, upper_transition_geom_object_pt, s_transition_up);
291 
292  // Find the Eulerian coordinates of the walls at the transition region
293  lower_transition_geom_object_pt->position(s_transition_lo, r_wall_lo);
294  upper_transition_geom_object_pt->position(s_transition_up, r_wall_up);
295 
296  // Take the average of these positions to define the origin of the spines
297  // in the transition region Horizontal position is always halfway
298  spine_centre[0] = 0.5 * (r_wall_lo[0] + r_wall_up[0]);
299 
300  // Vertical Position is given by a specified fraction
301  // between the upper and lower walls
302  spine_centre[1] =
303  r_wall_lo[1] + spine_centre_fraction() * (r_wall_up[1] - r_wall_lo[1]);
304  }
305 
306 
307  // Loop over elements in lower horizontal transition region
308  // --------------------------------------------------------
309  {
310  // Increments in wall coordinate
311  double dzeta_el = d / double(nx2);
312  double dzeta_node = d / double(nx2 * (np - 1));
313 
314  // Loop over elements horizontally
315  for (unsigned i = 0; i < nx2; i++)
316  {
317  // Start of wall coordinate
318  double zeta_lo = Zeta_transition_start + double(i) * dzeta_el;
319 
320  // Loop over elements vertically
321  for (unsigned j = 0; j < nh; j++)
322  {
323  // Work out element number in overall mesh
324  unsigned e = n_prev_elements + i * (nh + 1) + j;
325 
326  // Get pointer to element
327  FiniteElement* el_pt = this->finite_element_pt(e);
328 
329  // Loop over its nodes
330  for (unsigned i0 = 0; i0 < np; i0++)
331  {
332  for (unsigned i1 = 0; i1 < np; i1++)
333  {
334  // Node number:
335  unsigned n = i0 * np + i1;
336 
337  // Get spine node
338  SpineNode* nod_pt = dynamic_cast<SpineNode*>(el_pt->node_pt(n));
339 
340  // Set update id
341  nod_pt->node_update_fct_id() = 1;
342 
343  // Provide spine with additional storage for wall coordinate
344  if (i0 == 0)
345  {
346  // Get the Lagrangian coordinate in the Lower Wall
347  zeta[0] = zeta_lo + double(i1) * dzeta_node;
348  // Get the sub geometric object and local coordinate
349  Lower_wall_pt->locate_zeta(
350  zeta, lower_sub_geom_object_pt, s_lo);
351 
352  // Pass geometric parameter to the spine
353  Vector<double> parameters(3);
354  parameters[0] = s_lo[0];
355  parameters[1] = s_transition_lo[0];
356  parameters[2] = s_transition_up[0];
357  nod_pt->spine_pt()->set_geom_parameter(parameters);
358 
359  // Get position vector to wall
360  lower_sub_geom_object_pt->position(s_lo, r_wall_lo);
361 
362  // Get normal vector towards origin
363  Vector<double> N(2);
364  N[0] = spine_centre[0] - r_wall_lo[0];
365  N[1] = spine_centre[1] - r_wall_lo[1];
366  double length = sqrt(N[0] * N[0] + N[1] * N[1]);
367  nod_pt->spine_pt()->height() = length - radius;
368 
369  // Lower sub geom object is one (and only) geom object
370  // for spine:
371  Vector<GeomObject*> geom_object_pt(3);
372  geom_object_pt[0] = lower_sub_geom_object_pt;
373  geom_object_pt[1] = lower_transition_geom_object_pt;
374  geom_object_pt[2] = upper_transition_geom_object_pt;
375 
376  // Pass geom object(s) to spine
377  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
378 
379  // Push the node back onto boundaries
380  if (j == 0) set_boundary_node_pt[0].insert(nod_pt);
381  }
382  }
383  }
384  }
385  }
386 
387  // Increment number of previous elements
388  n_prev_elements += nx2 * (nh + 1);
389  }
390 
391  // Loop over elements in lower vertical transition region
392  // --------------------------------------------------------
393  {
394  for (unsigned i = 0; i < nhalf; i++)
395  {
396  // Loop over elements vertically
397  for (unsigned j = 0; j < nh; j++)
398  {
399  // Work out element number in overall mesh
400  unsigned e = n_prev_elements + i * (nh + 1) + j;
401 
402  // Get pointer to element
403  FiniteElement* el_pt = this->finite_element_pt(e);
404 
405  // Loop over its nodes
406  for (unsigned i0 = 0; i0 < np; i0++)
407  {
408  for (unsigned i1 = 0; i1 < np; i1++)
409  {
410  // Node number:
411  unsigned n = i0 * np + i1;
412 
413  // Get spine node
414  SpineNode* nod_pt = dynamic_cast<SpineNode*>(el_pt->node_pt(n));
415 
416  // Set update id
417  nod_pt->node_update_fct_id() = 2;
418 
419  // Provide spine with additional storage for fraction along
420  // vertical line
421  if (i0 == 0)
422  {
423  // Get position vectors to wall
424  zeta[0] = Zeta_transition_end;
425  // Get the sub geometric objects
426  Lower_wall_pt->locate_zeta(
427  zeta, lower_sub_geom_object_pt, s_lo);
428  Upper_wall_pt->locate_zeta(
429  zeta, upper_sub_geom_object_pt, s_up);
430 
431  lower_sub_geom_object_pt->position(s_lo, r_wall_lo);
432  upper_sub_geom_object_pt->position(s_up, r_wall_up);
433 
434  // Set vertical fraction
435  double vertical_fraction =
436  (double(i) + double(i1) / double(np - 1)) /
437  double(2.0 * nhalf);
438 
439  // Add the geometric parameters in order
440  Vector<double> parameters(5);
441  parameters[0] = s_lo[0];
442  parameters[1] = s_up[0];
443  parameters[2] = vertical_fraction;
444  parameters[3] = s_transition_lo[0];
445  parameters[4] = s_transition_up[0];
446  nod_pt->spine_pt()->set_geom_parameter(parameters);
447 
448  // Origin of spine
449  Vector<double> S0(2);
450  S0[0] = r_wall_lo[0] +
451  vertical_fraction * (r_wall_up[0] - r_wall_lo[0]);
452  S0[1] = r_wall_lo[1] +
453  vertical_fraction * (r_wall_up[1] - r_wall_lo[1]);
454 
455  // Get normal vector towards origin
456  Vector<double> N(2);
457  N[0] = spine_centre[0] - S0[0];
458  N[1] = spine_centre[1] - S0[1];
459 
460  double length = sqrt(N[0] * N[0] + N[1] * N[1]);
461  nod_pt->spine_pt()->height() = length - radius;
462 
463  // Lower and Upper wall sub geom objects affect spine:
464  Vector<GeomObject*> geom_object_pt(4);
465  geom_object_pt[0] = lower_sub_geom_object_pt;
466  geom_object_pt[1] = upper_sub_geom_object_pt;
467  geom_object_pt[2] = lower_transition_geom_object_pt;
468  geom_object_pt[3] = upper_transition_geom_object_pt;
469 
470  // Pass geom object(s) to spine
471  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
472 
473  // Push the node back onto boundaries
474  if (j == 0) set_boundary_node_pt[1].insert(nod_pt);
475  }
476  }
477  }
478  }
479  }
480 
481  // Increment number of previous elements
482  n_prev_elements += nhalf * (nh + 1);
483  }
484 
485 
486  // Loop over elements in upper vertical transition region
487  // ------------------------------------------------------
488  {
489  for (unsigned i = 0; i < nhalf; i++)
490  {
491  // Loop over elements vertically
492  for (unsigned j = 0; j < nh; j++)
493  {
494  // Work out element number in overall mesh
495  unsigned e = n_prev_elements + i * (nh + 1) + j;
496 
497  // Get pointer to element
498  FiniteElement* el_pt = this->finite_element_pt(e);
499 
500  // Loop over its nodes
501  for (unsigned i0 = 0; i0 < np; i0++)
502  {
503  for (unsigned i1 = 0; i1 < np; i1++)
504  {
505  // Node number:
506  unsigned n = i0 * np + i1;
507 
508  // Get spine node
509  SpineNode* nod_pt = dynamic_cast<SpineNode*>(el_pt->node_pt(n));
510 
511  // Set update id
512  nod_pt->node_update_fct_id() = 3;
513 
514  // Provide spine with additional storage for fraction along
515  // vertical line
516  if (i0 == 0)
517  {
518  // Get position vectors to wall
519  zeta[0] = Zeta_transition_end;
520  // Get the sub geometric objects
521  Lower_wall_pt->locate_zeta(
522  zeta, lower_sub_geom_object_pt, s_lo);
523  Upper_wall_pt->locate_zeta(
524  zeta, upper_sub_geom_object_pt, s_up);
525 
526  lower_sub_geom_object_pt->position(s_lo, r_wall_lo);
527  upper_sub_geom_object_pt->position(s_up, r_wall_up);
528 
529  // Set vertical fraction
530  double vertical_fraction =
531  0.5 + (double(i) + double(i1) / double(np - 1)) /
532  double(2.0 * nhalf);
533 
534  // Add the geometric parameters in order
535  Vector<double> parameters(5);
536  parameters[0] = s_lo[0];
537  parameters[1] = s_up[0];
538  parameters[2] = vertical_fraction;
539  parameters[3] = s_transition_lo[0];
540  parameters[4] = s_transition_up[0];
541  nod_pt->spine_pt()->set_geom_parameter(parameters);
542 
543  // Origin of spine
544  Vector<double> S0(2);
545  S0[0] = r_wall_lo[0] +
546  vertical_fraction * (r_wall_up[0] - r_wall_lo[0]);
547  S0[1] = r_wall_lo[1] +
548  vertical_fraction * (r_wall_up[1] - r_wall_lo[1]);
549 
550  // Get normal vector towards origin
551  Vector<double> N(2);
552  N[0] = spine_centre[0] - S0[0];
553  N[1] = spine_centre[1] - S0[1];
554 
555  double length = sqrt(N[0] * N[0] + N[1] * N[1]);
556  nod_pt->spine_pt()->height() = length - radius;
557 
558  // Upper and Lower wall geom objects affect spine
559  Vector<GeomObject*> geom_object_pt(4);
560  geom_object_pt[0] = lower_sub_geom_object_pt;
561  geom_object_pt[1] = upper_sub_geom_object_pt;
562  geom_object_pt[2] = lower_transition_geom_object_pt;
563  geom_object_pt[3] = upper_transition_geom_object_pt;
564 
565  // Pass geom object(s) to spine
566  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
567 
568  // Push the node back onto boundaries
569  if (j == 0) set_boundary_node_pt[1].insert(nod_pt);
570  }
571  }
572  }
573  }
574  }
575  // Increment number of previous elements
576  n_prev_elements += nhalf * (nh + 1);
577  }
578 
579 
580  // Loop over elements in upper horizontal transition region
581  // --------------------------------------------------------
582  {
583  // Increments in wall coordinate
584  double dzeta_el = d / double(nx2);
585  double dzeta_node = d / double(nx2 * (np - 1));
586 
587  // Loop over elements horizontally
588  for (unsigned i = 0; i < nx2; i++)
589  {
590  // Start of wall coordinate
591  double zeta_lo = Zeta_transition_end - double(i) * dzeta_el;
592 
593  // Loop over elements vertically
594  for (unsigned j = 0; j < nh; j++)
595  {
596  // Work out element number in overall mesh
597  unsigned e = n_prev_elements + i * (nh + 1) + j;
598 
599  // Get pointer to element
600  FiniteElement* el_pt = this->finite_element_pt(e);
601 
602  // Loop over its nodes
603  for (unsigned i0 = 0; i0 < np; i0++)
604  {
605  for (unsigned i1 = 0; i1 < np; i1++)
606  {
607  // Node number:
608  unsigned n = i0 * np + i1;
609 
610  // Get spine node
611  SpineNode* nod_pt = dynamic_cast<SpineNode*>(el_pt->node_pt(n));
612 
613  // Set update id
614  nod_pt->node_update_fct_id() = 4;
615 
616  // Provide spine with additional storage for wall coordinate
617  if (i0 == 0)
618  {
619  // Get the Lagrangian coordinate in the Upper wall
620  zeta[0] = zeta_lo - double(i1) * dzeta_node;
621  // Get the sub geometric object and local coordinate
622  Upper_wall_pt->locate_zeta(
623  zeta, upper_sub_geom_object_pt, s_up);
624 
625  // Pass geometric parameter to spine
626  Vector<double> parameters(3);
627  parameters[0] = s_up[0];
628  parameters[1] = s_transition_lo[0];
629  parameters[2] = s_transition_up[0];
630  nod_pt->spine_pt()->set_geom_parameter(parameters);
631 
632  // Get position vector to wall
633  upper_sub_geom_object_pt->position(s_up, r_wall_up);
634 
635  // Get normal vector towards origin
636  Vector<double> N(2);
637  N[0] = spine_centre[0] - r_wall_up[0];
638  N[1] = spine_centre[1] - r_wall_up[1];
639  double length = sqrt(N[0] * N[0] + N[1] * N[1]);
640  nod_pt->spine_pt()->height() = length - radius;
641 
642  // upper wall sub geom object is one (and only) geom object
643  // for spine:
644  Vector<GeomObject*> geom_object_pt(3);
645  geom_object_pt[0] = upper_sub_geom_object_pt;
646  geom_object_pt[1] = lower_transition_geom_object_pt;
647  geom_object_pt[2] = upper_transition_geom_object_pt;
648 
649  // Pass geom object(s) to spine
650  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
651 
652  // Push the node back onto boundaries
653  if (j == 0) set_boundary_node_pt[2].insert(nod_pt);
654  }
655  }
656  }
657  }
658  }
659 
660  // Increment number of previous elements
661  n_prev_elements += nx2 * (nh + 1);
662  }
663 
664 
665  // Loop over elements in upper deposited film region
666  // -------------------------------------------------
667  {
668  // Increments in wall coordinate
669  double dzeta_el = llayer / double(nx1);
670  double dzeta_node = llayer / double(nx1 * (np - 1));
671 
672  // Loop over elements horizontally
673  for (unsigned i = 0; i < nx1; i++)
674  {
675  // Start of wall coordinate
676  double zeta_lo = Zeta_transition_start - double(i) * dzeta_el;
677 
678  // Loop over elements vertically
679  for (unsigned j = 0; j < nh; j++)
680  {
681  // Work out element number in overall mesh
682  unsigned e = n_prev_elements + i * (nh + 1) + j;
683 
684  // Get pointer to element
685  FiniteElement* el_pt = this->finite_element_pt(e);
686 
687  // Loop over its nodes
688  for (unsigned i0 = 0; i0 < np; i0++)
689  {
690  for (unsigned i1 = 0; i1 < np; i1++)
691  {
692  // Node number:
693  unsigned n = i0 * np + i1;
694 
695  // Get spine node
696  SpineNode* nod_pt = dynamic_cast<SpineNode*>(el_pt->node_pt(n));
697 
698  // Set update id
699  nod_pt->node_update_fct_id() = 5;
700 
701  // Provide spine with additional storage for wall coordinate
702  if (i0 == 0)
703  {
704  // Get the Lagrangian coordinate in the Upper wall
705  zeta[0] = zeta_lo - double(i1) * dzeta_node;
706  // Get the geometric object and local coordinate
707  Upper_wall_pt->locate_zeta(
708  zeta, upper_sub_geom_object_pt, s_up);
709 
710  // The local coordinate is a geometric parameter
711  Vector<double> parameters(1, s_up[0]);
712  nod_pt->spine_pt()->set_geom_parameter(parameters);
713 
714  // Adjust spine height
715  nod_pt->spine_pt()->height() = H;
716 
717  // upper sub geom object is one (and only) geom object
718  // for spine:
719  Vector<GeomObject*> geom_object_pt(1);
720  geom_object_pt[0] = upper_sub_geom_object_pt;
721 
722  // Pass geom object(s) to spine
723  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
724 
725  // Push the node back onto boundaries
726  if (j == 0) set_boundary_node_pt[2].insert(nod_pt);
727  }
728  }
729  }
730  }
731  }
732  }
733 
734  // Wipe the boundary lookup schemes
735  this->remove_boundary_nodes();
736  this->set_nboundary(6);
737 
738  // Copy from sets to vectors
739  for (unsigned ibound = 0; ibound < 6; ibound++)
740  {
741  typedef std::set<Node*>::iterator IT;
742  for (IT it = set_boundary_node_pt[ibound].begin();
743  it != set_boundary_node_pt[ibound].end();
744  it++)
745  {
746  this->add_boundary_node(ibound, *it);
747  }
748  }
749 
750 
751  // Loop over the free surface boundary (4) and set a boundary coordinate
752  {
753  // Boundary coordinate
754  Vector<double> zeta(1);
755  unsigned n_boundary_node = this->nboundary_node(4);
756  for (unsigned n = 0; n < n_boundary_node; ++n)
757  {
758  zeta[0] = this->boundary_node_pt(4, n)->x(0);
759  this->boundary_node_pt(4, n)->set_coordinates_on_boundary(4, zeta);
760  }
761  }
762 
763 
764  // Now add the rectangular mesh in the channel ahead of the finger tip
765  //--------------------------------------------------------------------
766 
767  // Build a temporary version of a SimpleRectangularQuadMesh as
768  // a unit square
771  Nx3, 2 * nhalf, 1.0, 1.0, time_stepper_pt);
772 
773  // Wipe the boundary information in the auxilliary mesh
774  aux_mesh_pt->remove_boundary_nodes();
775 
776  // Copy all nodes in the auxiliary mesh into a set (from where they
777  // can easily be removed)
778  nnod = aux_mesh_pt->nnode();
779  std::set<Node*> aux_node_pt;
780  for (unsigned i = 0; i < nnod; i++)
781  {
782  aux_node_pt.insert(aux_mesh_pt->node_pt(i));
783  }
784 
785  // Loop over elements in first column and set pointers
786  // to the nodes in their first column to the ones that already exist
787 
788  // Boundary node number for first node in element
789  unsigned first_bound_node = 0;
790 
791  // Loop over elements
792  for (unsigned e = 0; e < 2 * nhalf; e++)
793  {
794  FiniteElement* el_pt = aux_mesh_pt->finite_element_pt(e * Nx3);
795  // Loop over first column of nodes
796  for (unsigned i = 0; i < np; i++)
797  {
798  // Node number in element
799  unsigned n = i * np;
800  // Remove node from set and kill it
801  if ((e < 1) || (i > 0))
802  {
803  aux_node_pt.erase(el_pt->node_pt(n));
804  delete el_pt->node_pt(n);
805  }
806  // Set pointer to existing node
807  el_pt->node_pt(n) = this->boundary_node_pt(1, first_bound_node + i);
808  }
809 
810  // Increment first node number
811  first_bound_node += np - 1;
812  }
813 
814  // Now add all the remaining nodes in the auxiliary mesh
815  // to the actual one
816  typedef std::set<Node*>::iterator IT;
817  for (IT it = aux_node_pt.begin(); it != aux_node_pt.end(); it++)
818  {
819  this->Node_pt.push_back(*it);
820  }
821 
822  // Store number of elements before the new bit gets added:
823  unsigned nelement_orig = this->Element_pt.size();
824 
825  // Add the elements to the actual mesh
826  unsigned nelem = aux_mesh_pt->nelement();
827  for (unsigned e = 0; e < nelem; e++)
828  {
829  this->Element_pt.push_back(aux_mesh_pt->element_pt(e));
830  // Don't forget to add them to the bulk
831  this->Bulk_element_pt.push_back(aux_mesh_pt->finite_element_pt(e));
832  }
833 
834  // Now wipe the boundary node storage scheme for the
835  // nodes that used to be at the end of the domain:
836  this->remove_boundary_nodes(1);
837 
838 
839  // FIRST SPINE
840  //-----------
841 
842  // Element 0
843  // Node 0
844  // Assign the new spine with unit height (pinned dummy -- never used)
845  Spine* new_spine_pt = new Spine(1.0);
846  this->Spine_pt.push_back(new_spine_pt);
847  new_spine_pt->spine_height_pt()->pin(0);
848 
849  // Get pointer to node
850  SpineNode* nod_pt = this->element_node_pt(nelement_orig + 0, 0);
851  // Set the pointer to node
852  nod_pt->spine_pt() = new_spine_pt;
853  // Set the fraction
854  nod_pt->fraction() = 0.0;
855  // Pointer to the mesh that implements the update fct
856  nod_pt->spine_mesh_pt() = this;
857  // ID for the update function
858  nod_pt->node_update_fct_id() = 6;
859 
860  // Need to get the local coordinates for the upper and lower wall
861  zeta[0] = Zeta_transition_end;
862  // Get the sub geometric objects
863  Lower_wall_pt->locate_zeta(zeta, lower_sub_geom_object_pt, s_lo);
864  Upper_wall_pt->locate_zeta(zeta, upper_sub_geom_object_pt, s_up);
865 
866  lower_sub_geom_object_pt->position(s_lo, r_wall_lo);
867  upper_sub_geom_object_pt->position(s_up, r_wall_up);
868 
869  // Pass additional data to spine
870  Vector<double> parameters(2);
871  parameters[0] = s_lo[0];
872  parameters[1] = s_up[0];
873  new_spine_pt->set_geom_parameter(parameters);
874 
875  // Lower and upper wall sub geom objects affect update
876  // for spine:
877  Vector<GeomObject*> geom_object_pt(2);
878  geom_object_pt[0] = lower_sub_geom_object_pt;
879  geom_object_pt[1] = upper_sub_geom_object_pt;
880 
881  // Pass geom object(s) to spine
882  new_spine_pt->set_geom_object_pt(geom_object_pt);
883 
884  // Loop vertically along the spine
885  // Loop over the elements
886  for (unsigned long i = 0; i < 2 * nhalf; i++)
887  {
888  // Loop over the vertical nodes, apart from the first
889  for (unsigned l1 = 1; l1 < np; l1++)
890  {
891  // Get pointer to node
892  SpineNode* nod_pt =
893  this->element_node_pt(nelement_orig + i * Nx3, l1 * np);
894  // Set the pointer to the spine
895  nod_pt->spine_pt() = new_spine_pt;
896  // Set the fraction
897  nod_pt->fraction() =
898  (double(i) + double(l1) / double(np - 1)) / double(2 * nhalf);
899  // Pointer to the mesh that implements the update fct
900  nod_pt->spine_mesh_pt() = this;
901  // ID for the update function
902  nod_pt->node_update_fct_id() = 6;
903  }
904  }
905 
906 
907  // LOOP OVER OTHER SPINES
908  //----------------------
909 
910  // Now loop over the elements horizontally
911  for (unsigned long j = 0; j < Nx3; j++)
912  {
913  // Loop over the nodes in the elements horizontally, ignoring
914  // the first column
915  for (unsigned l2 = 1; l2 < np; l2++)
916  {
917  // Assign the new spine with unit height (pinned dummy; never used)
918  new_spine_pt = new Spine(1.0);
919  this->Spine_pt.push_back(new_spine_pt);
920  new_spine_pt->spine_height_pt()->pin(0);
921 
922  // Get the node
923  SpineNode* nod_pt = this->element_node_pt(nelement_orig + j, l2);
924  // Set the pointer to the spine
925  nod_pt->spine_pt() = new_spine_pt;
926  // Set the fraction
927  nod_pt->fraction() = 0.0;
928  // Pointer to the mesh that implements the update fct
929  nod_pt->spine_mesh_pt() = this;
930  // ID for the update function
931  nod_pt->node_update_fct_id() = 6;
932 
933  // Add to boundary lookup scheme
934  this->add_boundary_node(0, nod_pt);
935  if ((j == Nx3 - 1) && (l2 == np - 1))
936  {
937  this->add_boundary_node(1, nod_pt);
938  }
939 
940  // Increment in wall coordinate
941  double dzeta_el = (Zeta_end - Zeta_transition_end) / double(Nx3);
942  double dzeta_nod = dzeta_el / double(np - 1);
943 
944  // Get wall coordinate
945  zeta[0] = Zeta_transition_end + j * dzeta_el + l2 * dzeta_nod;
946 
947  // Get the sub geometric objects
948  Lower_wall_pt->locate_zeta(zeta, lower_sub_geom_object_pt, s_lo);
949  Upper_wall_pt->locate_zeta(zeta, upper_sub_geom_object_pt, s_up);
950 
951  lower_sub_geom_object_pt->position(s_lo, r_wall_lo);
952  upper_sub_geom_object_pt->position(s_up, r_wall_up);
953 
954  // Add geometric parameters to spine
955  Vector<double> parameters(2);
956  parameters[0] = s_lo[0];
957  parameters[1] = s_up[0];
958  new_spine_pt->set_geom_parameter(parameters);
959 
960  // Lower and upper sub geom objects affect update
961  // for spine:
962  Vector<GeomObject*> geom_object_pt(2);
963  geom_object_pt[0] = lower_sub_geom_object_pt;
964  geom_object_pt[1] = upper_sub_geom_object_pt;
965 
966  // Pass geom object(s) to spine
967  new_spine_pt->set_geom_object_pt(geom_object_pt);
968 
969 
970  // Loop vertically along the spine
971  // Loop over the elements
972  for (unsigned long i = 0; i < 2 * nhalf; i++)
973  {
974  // Loop over the vertical nodes, apart from the first
975  for (unsigned l1 = 1; l1 < np; l1++)
976  {
977  // Get node
978  SpineNode* nod_pt =
979  this->element_node_pt(nelement_orig + i * Nx3 + j, l1 * np + l2);
980  // Set the pointer to the spine
981  nod_pt->spine_pt() = new_spine_pt;
982  // Set the fraction
983  nod_pt->fraction() =
984  (double(i) + double(l1) / double(np - 1)) / double(2 * nhalf);
985  // Pointer to the mesh that implements the update fct
986  nod_pt->spine_mesh_pt() = this;
987  // ID for the update function
988  nod_pt->node_update_fct_id() = 6;
989 
990  // Add to boundary lookup scheme
991  if ((j == Nx3 - 1) && (l2 == np - 1))
992  {
993  this->add_boundary_node(1, nod_pt);
994  }
995 
996  // Add to boundary lookup scheme
997  if ((i == 2 * nhalf - 1) && (l1 == np - 1))
998  {
999  this->add_boundary_node(2, nod_pt);
1000  }
1001  }
1002  }
1003  }
1004  }
1005 
1006  // (Re-)setup lookup scheme that establishes which elements are located
1007  // on the mesh boundaries
1008  this->setup_boundary_element_info();
1009 
1010  // Flush the storage for elements and nodes in the auxiliary mesh
1011  // so it can be safely deleted
1012  aux_mesh_pt->flush_element_and_node_storage();
1013 
1014  // Kill the auxiliary mesh
1015  delete aux_mesh_pt;
1016  }
1017 
1018 
1019  //======================================================================
1020  /// Reorder elements: Vertical stacks of elements, each topped by
1021  /// their interface element -- this is (currently) identical to the
1022  /// version in the SingleLayerSpineMesh but it's important
1023  /// that element reordering is maintained in exactly this form
1024  /// so to be on the safe side, we move the function in here.
1025  //======================================================================
1026  template<class ELEMENT, class INTERFACE_ELEMENT>
1027  void BrethertonSpineMesh<ELEMENT,
1028  INTERFACE_ELEMENT>::initial_element_reorder()
1029  {
1030  unsigned Nx = this->Nx;
1031  unsigned Ny = this->Ny;
1032  // Find out how many elements there are
1033  unsigned long Nelement = this->nelement();
1034  // Find out how many fluid elements there are
1035  unsigned long Nfluid = Nx * Ny;
1036  // Create a dummy array of elements
1037  Vector<FiniteElement*> dummy;
1038 
1039  // Loop over the elements in horizontal order
1040  for (unsigned long j = 0; j < Nx; j++)
1041  {
1042  // Loop over the elements in lower layer vertically
1043  for (unsigned long i = 0; i < Ny; i++)
1044  {
1045  // Push back onto the new stack
1046  dummy.push_back(this->finite_element_pt(Nx * i + j));
1047  }
1048 
1049  // Push back the line element onto the stack
1050  dummy.push_back(this->finite_element_pt(Nfluid + j));
1051  }
1052 
1053  // Now copy the reordered elements into the element_pt
1054  for (unsigned long e = 0; e < Nelement; e++)
1055  {
1056  this->Element_pt[e] = dummy[e];
1057  }
1058  }
1059 
1060  //=======================================================================
1061  /// Calculate the distance from the spine base to the free surface, i.e.
1062  /// the spine height.
1063  //=======================================================================
1064  template<class ELEMENT, class INTERFACE_ELEMENT>
1066  find_distance_to_free_surface(GeomObject* const& fs_geom_object_pt,
1067  Vector<double>& initial_zeta,
1068  const Vector<double>& spine_base,
1069  const Vector<double>& spine_end)
1070  {
1071  // Now we need to find the intersection points
1072  double lambda = 0.5;
1073 
1074  Vector<double> dx(2);
1075  Vector<double> new_free_x(2);
1076  DenseDoubleMatrix jacobian(2, 2, 0.0);
1077  Vector<double> spine_x(2);
1078  Vector<double> free_x(2);
1079 
1080  double maxres = 100.0;
1081 
1082  unsigned count = 0;
1083  // Let's Newton method it
1084  do
1085  {
1086  count++;
1087  // Find the spine's location
1088  for (unsigned i = 0; i < 2; ++i)
1089  {
1090  spine_x[i] = spine_base[i] + lambda * (spine_end[i] - spine_base[i]);
1091  }
1092 
1093  // Find the free_surface location
1094  fs_geom_object_pt->position(initial_zeta, free_x);
1095 
1096  for (unsigned i = 0; i < 2; ++i)
1097  {
1098  dx[i] = spine_x[i] - free_x[i];
1099  }
1100 
1101  // Calculate the entries of the jacobian matrix
1102  jacobian(0, 0) = (spine_end[0] - spine_base[0]);
1103  jacobian(1, 0) = (spine_end[1] - spine_base[1]);
1104 
1105  // Calculate the others by finite differences
1106  double FD_Jstep = 1.0e-6;
1107  double old_zeta = initial_zeta[0];
1108  initial_zeta[0] += FD_Jstep;
1109  fs_geom_object_pt->position(initial_zeta, new_free_x);
1110 
1111  for (unsigned i = 0; i < 2; ++i)
1112  {
1113  jacobian(i, 1) = (free_x[i] - new_free_x[i]) / FD_Jstep;
1114  }
1115 
1116  // Now let's solve the damn thing
1117  jacobian.solve(dx);
1118 
1119  lambda -= dx[0];
1120  initial_zeta[0] = old_zeta - dx[1];
1121 
1122  // How are we doing
1123 
1124  for (unsigned i = 0; i < 2; ++i)
1125  {
1126  spine_x[i] = spine_base[i] + lambda * (spine_end[i] - spine_base[i]);
1127  }
1128  fs_geom_object_pt->position(initial_zeta, free_x);
1129 
1130  for (unsigned i = 0; i < 2; ++i)
1131  {
1132  dx[i] = spine_x[i] - free_x[i];
1133  }
1134  maxres = std::fabs(dx[0]) > std::fabs(dx[1]) ? std::fabs(dx[0]) :
1135  std::fabs(dx[1]);
1136 
1137  if (count > 100)
1138  {
1139  throw OomphLibError("Failed to converge after 100 steps",
1140  OOMPH_CURRENT_FUNCTION,
1141  OOMPH_EXCEPTION_LOCATION);
1142  }
1143 
1144  } while (maxres > 1.0e-8);
1145 
1146 
1147  oomph_info << "DONE " << initial_zeta[0] << std::endl;
1148  double spine_length = sqrt(pow((spine_base[0] - spine_end[0]), 2.0) +
1149  pow((spine_base[1] - spine_end[1]), 2.0));
1150 
1151  return (lambda * spine_length); // Divided by spine length
1152  }
1153 
1154  //=======================================================================
1155  /// Reposition the spines that emenate from the lower wall
1156  //=======================================================================
1157  template<class ELEMENT, class INTERFACE_ELEMENT>
1159  const double& zeta_lo_transition_start,
1160  const double& zeta_lo_transition_end,
1161  const double& zeta_up_transition_start,
1162  const double& zeta_up_transition_end)
1163  {
1164  // Firstly create a geometric object of the free surface
1165  Mesh* fs_mesh_pt = new Mesh;
1166  this->template build_face_mesh<ELEMENT, FaceElementAsGeomObject>(
1167  4, fs_mesh_pt);
1168 
1169  // Loop over these new face elements and set the boundary number
1170  // of the bulk mesh
1171  unsigned n_face_element = fs_mesh_pt->nelement();
1172  // Loop over the elements
1173  for (unsigned e = 0; e < n_face_element; e++)
1174  {
1175  // Cast the element pointer to the correct thing!
1176  dynamic_cast<FaceElementAsGeomObject<ELEMENT>*>(fs_mesh_pt->element_pt(e))
1177  ->set_boundary_number_in_bulk_mesh(4);
1178  }
1179 
1180  // Now make a single geometric object that represents the collection of
1181  // geometric objects that form the boundary of the bulk mesh. Two
1182  // Eulerian coordinates, one intrinsic coordinate.
1183  MeshAsGeomObject* fs_geom_object_pt = new MeshAsGeomObject(fs_mesh_pt);
1184 
1185 
1186  // Length of deposited film region
1187  double llayer_lower = zeta_lo_transition_start - Zeta_start;
1188  double llayer_upper = zeta_up_transition_start - Zeta_start;
1189 
1190  // Length of transition region
1191  double d_lower = zeta_lo_transition_end - zeta_lo_transition_start;
1192  double d_upper = zeta_up_transition_end - zeta_up_transition_start;
1193 
1194  // Work out radius of circular cap from lower and upper wall
1195  Vector<double> r_wall_lo(2), r_wall_up(2);
1196  Vector<double> zeta(1), s_lo(1), s_up(1);
1197  GeomObject *lower_sub_geom_object_pt = 0, *upper_sub_geom_object_pt = 0;
1198 
1199  GeomObject* lower_transition_geom_object_pt = 0;
1200  GeomObject* upper_transition_geom_object_pt = 0;
1201  Vector<double> s_transition_lo(1), s_transition_up(1);
1202  Vector<double> spine_centre(2);
1203 
1204  // Get number of nodes along element edge
1205  unsigned np = this->finite_element_pt(0)->nnode_1d();
1206 
1207  // Calculate the centre for the spine nodes in the transition region
1208  {
1209  // Get the geometric objects on the walls at the start of the transition
1210  // region
1211  // Lower wall
1212  zeta[0] = zeta_lo_transition_start;
1213  Lower_wall_pt->locate_zeta(
1214  zeta, lower_transition_geom_object_pt, s_transition_lo);
1215  // Upper wall
1216  zeta[0] = zeta_up_transition_start;
1217  Upper_wall_pt->locate_zeta(
1218  zeta, upper_transition_geom_object_pt, s_transition_up);
1219 
1220  // Find the Eulerian coordinates of the walls at the transition region
1221  lower_transition_geom_object_pt->position(s_transition_lo, r_wall_lo);
1222  upper_transition_geom_object_pt->position(s_transition_up, r_wall_up);
1223 
1224  // Take the average of these positions to define the origin of the spines
1225  // in the transition region Horizontal position is always halfway
1226  spine_centre[0] = 0.5 * (r_wall_lo[0] + r_wall_up[0]);
1227 
1228  // Vertical Position is given by a specified fraction
1229  // between the upper and lower walls
1230  spine_centre[1] =
1231  r_wall_lo[1] + spine_centre_fraction() * (r_wall_up[1] - r_wall_lo[1]);
1232  }
1233 
1234 
1235  // Initialise number of elements in previous regions:
1236  unsigned n_prev_elements = 0;
1237 
1238  // Storage for the end of the spines
1239  Vector<double> spine_end(2);
1240  Vector<double> fs_zeta(1, 0.0);
1241 
1242  // Loop over elements in lower deposited film region
1243  // -------------------------------------------------
1244  {
1245  oomph_info << "LOWER FILM " << std::endl;
1246  // Increments in wall coordinate
1247  double dzeta_el = llayer_lower / double(Nx1);
1248  double dzeta_node = llayer_lower / double(Nx1 * (np - 1));
1249 
1250  // Loop over elements horizontally
1251  for (unsigned i = 0; i < Nx1; i++)
1252  {
1253  // Start of wall coordinate
1254  double zeta_lo = Zeta_start + double(i) * dzeta_el;
1255 
1256  // Work out element number in overall mesh
1257  unsigned e = n_prev_elements + i * (Nh + 1);
1258 
1259  // Get pointer to lower element
1260  FiniteElement* el_pt = this->finite_element_pt(e);
1261 
1262  // Loop over its nodes "horizontally"
1263  for (unsigned i1 = 0; i1 < (np - 1); i1++)
1264  {
1265  // Get spine node
1266  SpineNode* nod_pt = dynamic_cast<SpineNode*>(el_pt->node_pt(i1));
1267 
1268  // Get the Lagrangian coordinate in the Lower Wall
1269  zeta[0] = zeta_lo + double(i1) * dzeta_node;
1270  // Reset the boundary coordinate
1271  nod_pt->set_coordinates_on_boundary(0, zeta);
1272  // Get the geometric object and local coordinate
1273  Lower_wall_pt->locate_zeta(zeta, lower_sub_geom_object_pt, s_lo);
1274 
1275  // The local coordinate is a geometric parameter
1276  // This needs to be set (rather than added) because the
1277  // same spine may be visited more than once
1278  Vector<double> parameters(1, s_lo[0]);
1279  nod_pt->spine_pt()->set_geom_parameter(parameters);
1280 
1281  // The sub geom object is one (and only) geom object
1282  // for spine:
1283  Vector<GeomObject*> geom_object_pt(1);
1284  geom_object_pt[0] = lower_sub_geom_object_pt;
1285 
1286  // Pass geom object(s) to spine
1287  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
1288 
1289  // Get the wall position at the bottom of the spine
1290  lower_sub_geom_object_pt->position(s_lo, r_wall_lo);
1291  // The end of the spine is vertically above the base
1292  spine_end[0] = r_wall_lo[0];
1293  spine_end[1] = spine_centre[1];
1294  nod_pt->spine_pt()->height() = find_distance_to_free_surface(
1295  fs_geom_object_pt, fs_zeta, r_wall_lo, spine_end);
1296  }
1297  }
1298 
1299  // Increment number of previous elements
1300  n_prev_elements += Nx1 * (Nh + 1);
1301  }
1302 
1303 
1304  // Loop over elements in lower horizontal transition region
1305  // --------------------------------------------------------
1306  {
1307  oomph_info << "LOWER HORIZONTAL TRANSITION " << std::endl;
1308  // Increments in wall coordinate
1309  double dzeta_el = d_lower / double(Nx2);
1310  double dzeta_node = d_lower / double(Nx2 * (np - 1));
1311 
1312  // Loop over elements horizontally
1313  for (unsigned i = 0; i < Nx2; i++)
1314  {
1315  // Start of wall coordinate
1316  double zeta_lo = zeta_lo_transition_start + double(i) * dzeta_el;
1317 
1318  // Work out element number in overall mesh
1319  unsigned e = n_prev_elements + i * (Nh + 1);
1320 
1321  // Get pointer to element
1322  FiniteElement* el_pt = this->finite_element_pt(e);
1323 
1324  // Loop over its nodes
1325  for (unsigned i1 = 0; i1 < (np - 1); i1++)
1326  {
1327  // Get spine node
1328  SpineNode* nod_pt = dynamic_cast<SpineNode*>(el_pt->node_pt(i1));
1329 
1330  // Get the Lagrangian coordinate in the Lower Wall
1331  zeta[0] = zeta_lo + double(i1) * dzeta_node;
1332  // Reset the boundary coordinate
1333  nod_pt->set_coordinates_on_boundary(0, zeta);
1334  // Get the sub geometric object and local coordinate
1335  Lower_wall_pt->locate_zeta(zeta, lower_sub_geom_object_pt, s_lo);
1336 
1337  // Pass geometric parameter to the spine
1338  Vector<double> parameters(3);
1339  parameters[0] = s_lo[0];
1340  parameters[1] = s_transition_lo[0];
1341  parameters[2] = s_transition_up[0];
1342  nod_pt->spine_pt()->set_geom_parameter(parameters);
1343 
1344  // Lower sub geom object is one (and only) geom object
1345  // for spine:
1346  Vector<GeomObject*> geom_object_pt(3);
1347  geom_object_pt[0] = lower_sub_geom_object_pt;
1348  geom_object_pt[1] = lower_transition_geom_object_pt;
1349  geom_object_pt[2] = upper_transition_geom_object_pt;
1350 
1351  // Pass geom object(s) to spine
1352  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
1353 
1354  // Get position vector to wall
1355  lower_sub_geom_object_pt->position(s_lo, r_wall_lo);
1356  // The end of the spine is the spine centre,so the height is easy(ish)
1357  nod_pt->spine_pt()->height() = find_distance_to_free_surface(
1358  fs_geom_object_pt, fs_zeta, r_wall_lo, spine_centre);
1359  }
1360  }
1361 
1362  // Increment number of previous elements
1363  n_prev_elements += Nx2 * (Nh + 1);
1364  }
1365 
1366  // Loop over elements in lower vertical transition region
1367  // --------------------------------------------------------
1368  {
1369  oomph_info << "LOWER VERTICAL TRANSITION " << std::endl;
1370  for (unsigned i = 0; i < Nhalf; i++)
1371  {
1372  // Work out element number in overall mesh
1373  unsigned e = n_prev_elements + i * (Nh + 1);
1374 
1375  // Get pointer to element
1376  FiniteElement* el_pt = this->finite_element_pt(e);
1377 
1378  // Loop over its nodes
1379  for (unsigned i1 = 0; i1 < (np - 1); i1++)
1380  {
1381  // Get spine node
1382  // Note that I have to loop over the second row of nodes
1383  // because the first row are updated in region 6 and so
1384  // you get the wrong spines from them (doh!)
1385  SpineNode* nod_pt = dynamic_cast<SpineNode*>(el_pt->node_pt(np + i1));
1386 
1387  // Get position vectors to wall
1388  zeta[0] = zeta_lo_transition_end;
1389  // Get the sub geometric objects
1390  Lower_wall_pt->locate_zeta(zeta, lower_sub_geom_object_pt, s_lo);
1391  zeta[0] = zeta_up_transition_end;
1392  Upper_wall_pt->locate_zeta(zeta, upper_sub_geom_object_pt, s_up);
1393 
1394  lower_sub_geom_object_pt->position(s_lo, r_wall_lo);
1395  upper_sub_geom_object_pt->position(s_up, r_wall_up);
1396 
1397  // Set vertical fraction
1398  double vertical_fraction =
1399  (double(i) + double(i1) / double(np - 1)) / double(2.0 * Nhalf);
1400 
1401  // Add the geometric parameters in order
1402  Vector<double> parameters(5);
1403  parameters[0] = s_lo[0];
1404  parameters[1] = s_up[0];
1405  parameters[2] = vertical_fraction;
1406  parameters[3] = s_transition_lo[0];
1407  parameters[4] = s_transition_up[0];
1408  nod_pt->spine_pt()->set_geom_parameter(parameters);
1409 
1410  // Origin of spine
1411  Vector<double> S0(2);
1412  S0[0] =
1413  r_wall_lo[0] + vertical_fraction * (r_wall_up[0] - r_wall_lo[0]);
1414  S0[1] =
1415  r_wall_lo[1] + vertical_fraction * (r_wall_up[1] - r_wall_lo[1]);
1416 
1417  // Lower and Upper wall sub geom objects affect spine:
1418  Vector<GeomObject*> geom_object_pt(4);
1419  geom_object_pt[0] = lower_sub_geom_object_pt;
1420  geom_object_pt[1] = upper_sub_geom_object_pt;
1421  geom_object_pt[2] = lower_transition_geom_object_pt;
1422  geom_object_pt[3] = upper_transition_geom_object_pt;
1423 
1424  // Pass geom object(s) to spine
1425  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
1426 
1427  // Calculate the spine height
1428  nod_pt->spine_pt()->height() = find_distance_to_free_surface(
1429  fs_geom_object_pt, fs_zeta, S0, spine_centre);
1430  }
1431  }
1432 
1433  // Increment number of previous elements
1434  n_prev_elements += Nhalf * (Nh + 1);
1435  }
1436 
1437 
1438  // Loop over elements in upper vertical transition region
1439  // --------------------------------------------------------
1440  {
1441  oomph_info << "UPPER VERTICAL TRANSITION" << std::endl;
1442  for (unsigned i = 0; i < Nhalf; i++)
1443  {
1444  // Work out element number in overall mesh
1445  unsigned e = n_prev_elements + i * (Nh + 1);
1446 
1447  // Get pointer to element
1448  FiniteElement* el_pt = this->finite_element_pt(e);
1449 
1450  // Loop over its nodes
1451  for (unsigned i1 = 0; i1 < (np - 1); i1++)
1452  {
1453  // Get spine node
1454  // Note that I have to loop over the second row of nodes
1455  // because the first row are updated in region 6 and so
1456  // you get the wrong spines from them (doh!)
1457  SpineNode* nod_pt = dynamic_cast<SpineNode*>(el_pt->node_pt(np + i1));
1458 
1459  // Get position vectors to wall
1460  zeta[0] = zeta_lo_transition_end;
1461  // Get the sub geometric objects
1462  Lower_wall_pt->locate_zeta(zeta, lower_sub_geom_object_pt, s_lo);
1463  zeta[0] = zeta_up_transition_end;
1464  Upper_wall_pt->locate_zeta(zeta, upper_sub_geom_object_pt, s_up);
1465 
1466  lower_sub_geom_object_pt->position(s_lo, r_wall_lo);
1467  upper_sub_geom_object_pt->position(s_up, r_wall_up);
1468 
1469 
1470  // Set vertical fraction
1471  double vertical_fraction =
1472  0.5 +
1473  (double(i) + double(i1) / double(np - 1)) / double(2.0 * Nhalf);
1474 
1475  // Add the geometric parameters in order
1476  Vector<double> parameters(5);
1477  parameters[0] = s_lo[0];
1478  parameters[1] = s_up[0];
1479  parameters[2] = vertical_fraction;
1480  parameters[3] = s_transition_lo[0];
1481  parameters[4] = s_transition_up[0];
1482  nod_pt->spine_pt()->set_geom_parameter(parameters);
1483 
1484  // Origin of spine
1485  Vector<double> S0(2);
1486  S0[0] =
1487  r_wall_lo[0] + vertical_fraction * (r_wall_up[0] - r_wall_lo[0]);
1488  S0[1] =
1489  r_wall_lo[1] + vertical_fraction * (r_wall_up[1] - r_wall_lo[1]);
1490 
1491  // Lower and Upper wall sub geom objects affect spine:
1492  Vector<GeomObject*> geom_object_pt(4);
1493  geom_object_pt[0] = lower_sub_geom_object_pt;
1494  geom_object_pt[1] = upper_sub_geom_object_pt;
1495  geom_object_pt[2] = lower_transition_geom_object_pt;
1496  geom_object_pt[3] = upper_transition_geom_object_pt;
1497 
1498  // Pass geom object(s) to spine
1499  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
1500 
1501  // Calculate the spine height
1502  nod_pt->spine_pt()->height() = find_distance_to_free_surface(
1503  fs_geom_object_pt, fs_zeta, S0, spine_centre);
1504  }
1505  }
1506 
1507  // Increment number of previous elements
1508  n_prev_elements += Nhalf * (Nh + 1);
1509  }
1510 
1511 
1512  // Loop over elements in upper horizontal transition region
1513  // --------------------------------------------------------
1514  {
1515  oomph_info << "UPPER HORIZONTAL TRANSITION " << std::endl;
1516  // Increments in wall coordinate
1517  double dzeta_el = d_upper / double(Nx2);
1518  double dzeta_node = d_upper / double(Nx2 * (np - 1));
1519 
1520  // Loop over elements horizontally
1521  for (unsigned i = 0; i < Nx2; i++)
1522  {
1523  // Start of wall coordinate
1524  double zeta_lo = zeta_up_transition_end - double(i) * dzeta_el;
1525 
1526  // Work out element number in overall mesh
1527  unsigned e = n_prev_elements + i * (Nh + 1);
1528 
1529  // Get pointer to element
1530  FiniteElement* el_pt = this->finite_element_pt(e);
1531 
1532  // Loop over its nodes
1533  for (unsigned i1 = 0; i1 < (np - 1); i1++)
1534  {
1535  // Get spine node (same comment)
1536  SpineNode* nod_pt = dynamic_cast<SpineNode*>(el_pt->node_pt(np + i1));
1537 
1538  // Get the Lagrangian coordinate in the Lower Wall
1539  zeta[0] = zeta_lo - double(i1) * dzeta_node;
1540  // Reset the boundary coordinate
1541  el_pt->node_pt(i1)->set_coordinates_on_boundary(2, zeta);
1542  // Get the sub geometric object and local coordinate
1543  Upper_wall_pt->locate_zeta(zeta, upper_sub_geom_object_pt, s_up);
1544 
1545  // Pass geometric parameter to the spine
1546  Vector<double> parameters(3);
1547  parameters[0] = s_up[0];
1548  parameters[1] = s_transition_lo[0];
1549  parameters[2] = s_transition_up[0];
1550  nod_pt->spine_pt()->set_geom_parameter(parameters);
1551 
1552  // Lower sub geom object is one (and only) geom object
1553  // for spine:
1554  Vector<GeomObject*> geom_object_pt(3);
1555  geom_object_pt[0] = upper_sub_geom_object_pt;
1556  geom_object_pt[1] = lower_transition_geom_object_pt;
1557  geom_object_pt[2] = upper_transition_geom_object_pt;
1558 
1559  // Pass geom object(s) to spine
1560  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
1561 
1562 
1563  // Get position vector to wall
1564  upper_sub_geom_object_pt->position(s_up, r_wall_up);
1565  // Find spine height
1566  nod_pt->spine_pt()->height() = find_distance_to_free_surface(
1567  fs_geom_object_pt, fs_zeta, r_wall_up, spine_centre);
1568  }
1569  }
1570 
1571  // Increment number of previous elements
1572  n_prev_elements += Nx2 * (Nh + 1);
1573  }
1574 
1575 
1576  // Loop over elements in upper deposited film region
1577  // -------------------------------------------------
1578  {
1579  oomph_info << "UPPER THIN FILM" << std::endl;
1580  // Increments in wall coordinate
1581  double dzeta_el = llayer_upper / double(Nx1);
1582  double dzeta_node = llayer_upper / double(Nx1 * (np - 1));
1583 
1584  // Loop over elements horizontally
1585  for (unsigned i = 0; i < Nx1; i++)
1586  {
1587  // Start of wall coordinate
1588  double zeta_lo = zeta_up_transition_start - double(i) * dzeta_el;
1589 
1590  // Work out element number in overall mesh
1591  unsigned e = n_prev_elements + i * (Nh + 1);
1592 
1593  // Get pointer to element
1594  FiniteElement* el_pt = this->finite_element_pt(e);
1595 
1596  // Loop over its nodes "horizontally"
1597  for (unsigned i1 = 0; i1 < (np - 1); i1++)
1598  {
1599  // Get spine node
1600  SpineNode* nod_pt = dynamic_cast<SpineNode*>(el_pt->node_pt(i1));
1601 
1602  // Get the Lagrangian coordinate in the Upper wall
1603  zeta[0] = zeta_lo - double(i1) * dzeta_node;
1604  // Reset coordinate on boundary
1605  nod_pt->set_coordinates_on_boundary(2, zeta);
1606  // Get the geometric object and local coordinate
1607  Upper_wall_pt->locate_zeta(zeta, upper_sub_geom_object_pt, s_up);
1608 
1609  // The local coordinate is a geometric parameter
1610  Vector<double> parameters(1, s_up[0]);
1611  nod_pt->spine_pt()->set_geom_parameter(parameters);
1612 
1613  // upper sub geom object is one (and only) geom object
1614  // for spine:
1615  Vector<GeomObject*> geom_object_pt(1);
1616  geom_object_pt[0] = upper_sub_geom_object_pt;
1617 
1618  // Pass geom object(s) to spine
1619  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
1620 
1621  // Get the wall position at the bottom of the spine
1622  upper_sub_geom_object_pt->position(s_up, r_wall_up);
1623  spine_end[0] = r_wall_up[0];
1624  spine_end[1] = spine_centre[1];
1625  // Find the new spine height
1626  nod_pt->spine_pt()->height() = find_distance_to_free_surface(
1627  fs_geom_object_pt, fs_zeta, r_wall_up, spine_end);
1628  }
1629  }
1630 
1631 
1632  // Increment number of previous elements
1633  n_prev_elements += Nx1 * (Nh + 1);
1634  }
1635 
1636 
1637  // Additional mesh
1638  {
1639  unsigned e = n_prev_elements;
1640 
1641  // Get pointer to node
1642  SpineNode* nod_pt = this->element_node_pt(e, 0);
1643 
1644  // Need to get the local coordinates for the upper and lower wall
1645  zeta[0] = zeta_lo_transition_end;
1646  // Get the sub geometric objects
1647  Lower_wall_pt->locate_zeta(zeta, lower_sub_geom_object_pt, s_lo);
1648  zeta[0] = zeta_up_transition_end;
1649  Upper_wall_pt->locate_zeta(zeta, upper_sub_geom_object_pt, s_up);
1650 
1651  lower_sub_geom_object_pt->position(s_lo, r_wall_lo);
1652  upper_sub_geom_object_pt->position(s_up, r_wall_up);
1653 
1654  // Pass additional data to spine
1655  Vector<double> parameters(2);
1656  parameters[0] = s_lo[0];
1657  parameters[1] = s_up[0];
1658  nod_pt->spine_pt()->set_geom_parameter(parameters);
1659 
1660  // Lower and upper wall sub geom objects affect update
1661  // for spine:
1662  Vector<GeomObject*> geom_object_pt(2);
1663  geom_object_pt[0] = lower_sub_geom_object_pt;
1664  geom_object_pt[1] = upper_sub_geom_object_pt;
1665 
1666  // Pass geom object(s) to spine
1667  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
1668 
1669  // LOOP OVER OTHER SPINES
1670  //----------------------
1671 
1672  // Now loop over the elements horizontally
1673  for (unsigned long j = 0; j < Nx3; j++)
1674  {
1675  unsigned e = n_prev_elements + j;
1676 
1677  // Loop over the nodes in the elements horizontally, ignoring
1678  // the first column
1679  for (unsigned l2 = 0; l2 < np; l2++)
1680  {
1681  // Get pointer to node at the base of the spine
1682  SpineNode* nod_pt = this->element_node_pt(e, l2);
1683 
1684  // Increment in wall coordinate
1685  double dzeta_el_lower =
1686  (Zeta_end - zeta_lo_transition_end) / double(Nx3);
1687  double dzeta_nod_lower = dzeta_el_lower / double(np - 1);
1688 
1689  double dzeta_el_upper =
1690  (Zeta_end - zeta_up_transition_end) / double(Nx3);
1691  double dzeta_nod_upper = dzeta_el_upper / double(np - 1);
1692 
1693  // Get wall coordinate
1694  zeta[0] =
1695  zeta_lo_transition_end + j * dzeta_el_lower + l2 * dzeta_nod_lower;
1696  // Reset the boundary coordinate
1697  nod_pt->set_coordinates_on_boundary(0, zeta);
1698 
1699  // Get the sub geometric objects
1700  Lower_wall_pt->locate_zeta(zeta, lower_sub_geom_object_pt, s_lo);
1701 
1702  zeta[0] =
1703  zeta_up_transition_end + j * dzeta_el_upper + l2 * dzeta_nod_upper;
1704  // Reset the upper boundary coordinate
1705  this->element_node_pt(e + Nx3 * (2 * Nhalf - 1), np * (np - 1) + l2)
1706  ->set_coordinates_on_boundary(2, zeta);
1707  Upper_wall_pt->locate_zeta(zeta, upper_sub_geom_object_pt, s_up);
1708 
1709  lower_sub_geom_object_pt->position(s_lo, r_wall_lo);
1710  upper_sub_geom_object_pt->position(s_up, r_wall_up);
1711 
1712  // Add geometric parameters to spine
1713  Vector<double> parameters(2);
1714  parameters[0] = s_lo[0];
1715  parameters[1] = s_up[0];
1716  nod_pt->spine_pt()->set_geom_parameter(parameters);
1717 
1718  // Lower and upper sub geom objects affect update
1719  // for spine:
1720  Vector<GeomObject*> geom_object_pt(2);
1721  geom_object_pt[0] = lower_sub_geom_object_pt;
1722  geom_object_pt[1] = upper_sub_geom_object_pt;
1723 
1724  // Pass geom object(s) to spine
1725  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
1726  }
1727  }
1728  }
1729 
1730 
1731  // Clean up all the memory
1732  // Can delete the Geometric object
1733  delete fs_geom_object_pt;
1734  // Need to be careful with the FaceMesh, because we can't delete the nodes
1735  // Loop
1736  for (unsigned e = n_face_element; e > 0; e--)
1737  {
1738  delete fs_mesh_pt->element_pt(e - 1);
1739  fs_mesh_pt->element_pt(e - 1) = 0;
1740  }
1741  // Now clear all element and node storage (should maybe fine-grain this)
1742  fs_mesh_pt->flush_element_and_node_storage();
1743  // Finally delete the mesh
1744  delete fs_mesh_pt;
1745  }
1746 
1747 } // namespace oomph
1748 
1749 #endif
Mesh for 2D Bretherton problem – based on single layer mesh. Templated by spine-ified Navier-Stokes e...
GeomObject * Lower_wall_pt
Pointer to geometric object that represents the lower wall.
void reposition_spines(const double &zeta_lo_transition_start, const double &zeta_lo_transition_end, const double &zeta_up_transition_start, const double &zeta_up_transition_end)
Reposition the spines in response to changes in geometry.
BrethertonSpineMesh(const unsigned &nx1, const unsigned &nx2, const unsigned &nx3, const unsigned &nh, const unsigned &nhalf, const double &h, GeomObject *lower_wall_pt, GeomObject *upper_wall_pt, const double &zeta_start, const double &zeta_transition_start, const double &zeta_transition_end, const double &zeta_end, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor. Arguments:
double Zeta_start
Start coordinate on wall.
double Zeta_transition_end
Wall coordinate of end of transition region.
GeomObject * Upper_wall_pt
Pointer to geometric object that represents the upper wall.
Vector< FiniteElement * > Bulk_element_pt
Vector of pointers to element in the fluid layer.
unsigned Nx3
Number of elements along wall in channel region.
double H
Thickness of deposited film.
void initial_element_reorder()
Initial reordering elements of the elements – before the channel mesh is added. Vertical stacks of el...
double spine_centre_fraction() const
Read the value of the spine centre's vertical fraction.
Vector< FiniteElement * > Interface_element_pt
Vector of pointers to interface elements.
ELEMENT * Control_element_pt
Pointer to control element (just under the symmetry line near the bubble tip; the bubble tip is locat...
double Zeta_end
Wall coordinate of end of liquid filled region (inflow)
double find_distance_to_free_surface(GeomObject *const &fs_geom_object_pt, Vector< double > &initial_zeta, const Vector< double > &spine_base, const Vector< double > &spine_end)
Recalculate the spine lengths after repositioning.
double Zeta_transition_start
Wall coordinate of start of the transition region.
Simple rectangular 2D Quad mesh class. Nx : number of elements in the x direction.
Single-layer spine mesh class derived from standard 2D mesh. The mesh contains a layer of spinified f...
////////////////////////////////////////////////////////////////////// //////////////////////////////...