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