bretherton_spine_mesh.h
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-2025 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_HEADER
27#define OOMPH_BRETHERTON_SPINE_MESH_HEADER
28
29// The mesh
32
33namespace oomph
34{
35 //============================================================================
36 /// Mesh for 2D Bretherton problem -- based on single layer mesh. Templated
37 /// by spine-ified Navier-Stokes element type (e.g.
38 /// SpineElement<QCrouzeixRaviartElement<2> > and the corresponding
39 /// interface element (e.g.
40 /// SpineLineFluidInterfaceElement<SpineElement<QCrouzeixRaviartElement<2> > >
41 ///
42 //============================================================================
43 template<class ELEMENT, class INTERFACE_ELEMENT>
45 {
46 public:
47 /// Constructor. Arguments:
48 /// - nx1: Number of elements along wall in deposited film region
49 /// - nx2: Number of elements along wall in horizontal transition region
50 /// - nx3: Number of elements along wall in channel region
51 /// - nhalf: Number of elements in vertical transition region (there are
52 /// twice as many elements across the channel)
53 /// - nh: Number of elements across the deposited film
54 /// - h: Thickness of deposited film
55 /// - zeta0: Start coordinate on wall
56 /// - zeta1: Wall coordinate of start of transition region
57 /// - zeta2: Wall coordinate of end of liquid filled region (inflow)
58 /// - lower_wall_pt: Pointer to geometric object that represents the lower
59 /// wall
60 /// - upper_wall_pt: Pointer to geometric object that represents the upper
61 /// wall
62 /// - time_stepper_pt: Pointer to timestepper; defaults to Steady
64 const unsigned& nx1,
65 const unsigned& nx2,
66 const unsigned& nx3,
67 const unsigned& nh,
68 const unsigned& nhalf,
69 const double& h,
72 const double& zeta_start,
73 const double& zeta_transition_start,
74 const double& zeta_transition_end,
75 const double& zeta_end,
76 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper);
77
78 /// Access functions for pointers to interface elements
79 FiniteElement*& interface_element_pt(const unsigned long& i)
80 {
81 return Interface_element_pt[i];
82 }
83
84 /// Number of elements on interface
85 unsigned long ninterface_element() const
86 {
88 }
89
90 /// Access functions for pointers to elements in bulk
91 FiniteElement*& bulk_element_pt(const unsigned long& i)
92 {
93 return Bulk_element_pt[i];
94 }
95
96 /// Number of elements in bulk
97 unsigned long nbulk() const
98 {
99 return Bulk_element_pt.size();
100 }
101
102 /// Number of free-surface spines (i.e. excluding the dummy spines
103 /// in the channel region)
105 {
106 unsigned np = this->finite_element_pt(0)->nnode_1d();
107 return 2 * (Nx1 + Nx2 + Nhalf) * (np - 1);
108 }
109
110 /// Recalculate the spine lengths after repositioning
115
116 /// Reposition the spines in response to changes in geometry
118 const double& zeta_lo_transition_end,
119 const double& zeta_up_transition_start,
120 const double& zeta_up_transition_end);
121
122 /// Pin all spines so the mesh can be used for computation
123 /// without free surfaces
125 {
126 unsigned n_spine = this->nspine();
127 for (unsigned i = 0; i < n_spine; i++)
128 {
129 this->spine_pt(i)->spine_height_pt()->pin(0);
130 }
131 }
132
133 /// General node update function implements pure virtual function
134 /// defined in SpineMesh base class and performs specific update
135 /// actions, depending on the node update fct id stored for each node.
137 {
138 unsigned id = spine_node_pt->node_update_fct_id();
139 switch (id)
140 {
141 case 0:
143 break;
144
145 case 1:
147 break;
148
149 case 2:
151 break;
152
153 case 3:
155 break;
156
157 case 4:
159 break;
160
161 case 5:
163 break;
164
165 case 6:
167 break;
168
169 default:
170 std::ostringstream error_message;
171 error_message << "Incorrect spine update id " << id << std::endl;
172
173 throw OomphLibError(error_message.str(),
176 }
177 }
178
179
180 /// Pointer to control element (just under the symmetry line
181 /// near the bubble tip, so the bubble tip is located at
182 /// s=[1.0,1.0] in this element.
184 {
185 return Control_element_pt;
186 }
187
188
189 /// Read the value of the spine centre's vertical fraction
191 {
193 }
194
195 /// Set the pointer to the spine centre's vertial fraction
200
201 protected:
202 /// Update function for the deposited film region in the
203 /// lower part of the domain: Vertical spines
205 {
206 // Get fraction along the spine
207 double w = spine_node_pt->fraction();
208 // Get spine height
209 double h = spine_node_pt->h();
210
211 // Get wall coordinate
213 s_lo[0] = spine_node_pt->spine_pt()->geom_parameter(0);
214
215 // Get position vector to wall
217 spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_lo, r_wall_lo);
218
219 // Update the nodal position
220 spine_node_pt->x(0) = r_wall_lo[0];
221 spine_node_pt->x(1) = r_wall_lo[1] + w * h;
222 }
223
224
225 /// Update function for the horizontal transitition region in the
226 /// lower part of the domain: Spine points from wall to origin
228 {
229 // Get fraction along the spine
230 double w = spine_node_pt->fraction();
231 // Get spine height
232 double h = spine_node_pt->h();
233
234 // Get wall coordinate
236 s_lo[0] = spine_node_pt->spine_pt()->geom_parameter(0);
237
238 // Get position vector to wall
240 spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_lo, r_wall_lo);
241
242 // Get the spine centre
244 s_transition_lo[0] = spine_node_pt->spine_pt()->geom_parameter(1);
245 s_transition_up[0] = spine_node_pt->spine_pt()->geom_parameter(2);
247 spine_node_pt->spine_pt()->geom_object_pt(1)->position(s_transition_lo,
249 spine_node_pt->spine_pt()->geom_object_pt(2)->position(s_transition_up,
251
253 // Horizontal position is always halfway between the walls
254 spine_centre[0] = 0.5 * (r_transition_lo[0] + r_transition_up[0]);
255 // Vertical spine centre is given by fraction between the walls
256 spine_centre[1] =
257 r_transition_lo[1] +
259
260 // Get vector twoards spine origin
261 Vector<double> N(2);
262 N[0] = spine_centre[0] - r_wall_lo[0];
263 N[1] = spine_centre[1] - r_wall_lo[1];
264 double inv_length = 1.0 / sqrt(N[0] * N[0] + N[1] * N[1]);
265 // Update the nodal position
266 spine_node_pt->x(0) = r_wall_lo[0] + w * h * N[0] * inv_length;
267 spine_node_pt->x(1) = r_wall_lo[1] + w * h * N[1] * inv_length;
268 }
269
270
271 /// Update function for the vertical transitition region in the
272 /// lower part of the domain: Spine points to origin
274 {
275 // Get fraction along the spine
276 double w = spine_node_pt->fraction();
277 // Get spine height
278 double h = spine_node_pt->h();
279
280 // Get local coordinates
281 Vector<double> s_lo(1), s_up(1);
282 s_lo[0] = spine_node_pt->spine_pt()->geom_parameter(0);
283 s_up[0] = spine_node_pt->spine_pt()->geom_parameter(1);
284
285 // Get position vector to wall
286 Vector<double> r_lo(2), r_up(2);
287 spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_lo, r_lo);
288 spine_node_pt->spine_pt()->geom_object_pt(1)->position(s_up, r_up);
289
290 // Get fraction along vertical line across
291 double vertical_fraction = spine_node_pt->spine_pt()->geom_parameter(2);
292
293 // Origin of spine
295 S0[0] = r_lo[0] + vertical_fraction * (r_up[0] - r_lo[0]);
296 S0[1] = r_lo[1] + vertical_fraction * (r_up[1] - r_lo[1]);
297
298 // Get the spine centre
300 s_transition_lo[0] = spine_node_pt->spine_pt()->geom_parameter(3);
301 s_transition_up[0] = spine_node_pt->spine_pt()->geom_parameter(4);
303 spine_node_pt->spine_pt()->geom_object_pt(2)->position(s_transition_lo,
305 spine_node_pt->spine_pt()->geom_object_pt(3)->position(s_transition_up,
307
309 // Horizontal position is always halfway between the walls
310 spine_centre[0] = 0.5 * (r_transition_lo[0] + r_transition_up[0]);
311 // Vertical spine centre is given by fraction between the walls
312 spine_centre[1] =
313 r_transition_lo[1] +
315
316 // Get vector towards origin
317 Vector<double> N(2);
318 N[0] = spine_centre[0] - S0[0];
319 N[1] = spine_centre[1] - S0[1];
320
321 double inv_length = 1.0 / sqrt(N[0] * N[0] + N[1] * N[1]);
322 // Update the nodal position
323 spine_node_pt->x(0) = S0[0] + w * h * N[0] * inv_length;
324 spine_node_pt->x(1) = S0[1] + w * h * N[1] * inv_length;
325 }
326
327
328 /// Update function for the vertical transitition region in the
329 /// upper part of the domain: Spine points to origin
331 {
332 // Get fraction along the spine
333 double w = spine_node_pt->fraction();
334 // Get spine height
335 double h = spine_node_pt->h();
336
337 // Get local coordinates
338 Vector<double> s_lo(1), s_up(1);
339 s_lo[0] = spine_node_pt->spine_pt()->geom_parameter(0);
340 s_up[0] = spine_node_pt->spine_pt()->geom_parameter(1);
341
342 // Get position vector to wall
343 Vector<double> r_lo(2), r_up(2);
344 spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_lo, r_lo);
345 spine_node_pt->spine_pt()->geom_object_pt(1)->position(s_up, r_up);
346
347 // Get fraction along vertical line across
348 double vertical_fraction = spine_node_pt->spine_pt()->geom_parameter(2);
349
350 // Origin of spine
352 S0[0] = r_lo[0] + vertical_fraction * (r_up[0] - r_lo[0]);
353 S0[1] = r_lo[1] + vertical_fraction * (r_up[1] - r_lo[1]);
354
355 // Get the spine centre
357 s_transition_lo[0] = spine_node_pt->spine_pt()->geom_parameter(3);
358 s_transition_up[0] = spine_node_pt->spine_pt()->geom_parameter(4);
360 spine_node_pt->spine_pt()->geom_object_pt(2)->position(s_transition_lo,
362 spine_node_pt->spine_pt()->geom_object_pt(3)->position(s_transition_up,
364
366 // Horizontal position is always halfway between the walls
367 spine_centre[0] = 0.5 * (r_transition_lo[0] + r_transition_up[0]);
368 // Vertical spine centre is given by fraction between the walls
369 spine_centre[1] =
370 r_transition_lo[1] +
372
373 // Get vector towards origin
374 Vector<double> N(2);
375 N[0] = spine_centre[0] - S0[0];
376 N[1] = spine_centre[1] - S0[1];
377
378 double inv_length = 1.0 / sqrt(N[0] * N[0] + N[1] * N[1]);
379 // Update the nodal position
380 spine_node_pt->x(0) = S0[0] + w * h * N[0] * inv_length;
381 spine_node_pt->x(1) = S0[1] + w * h * N[1] * inv_length;
382 }
383
384
385 /// Update function for the horizontal transitition region in the
386 /// upper part of the domain: Spine points towards origin
388 {
389 // Get fraction along the spine
390 double w = spine_node_pt->fraction();
391 // Get spine height
392 double h = spine_node_pt->h();
393
394 // Get wall coordinate
396 s_up[0] = spine_node_pt->spine_pt()->geom_parameter(0);
397
398 // Get position vector to wall
400 spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_up, r_wall_up);
401
402 // Get the spine centre
404 s_transition_lo[0] = spine_node_pt->spine_pt()->geom_parameter(1);
405 s_transition_up[0] = spine_node_pt->spine_pt()->geom_parameter(2);
407 spine_node_pt->spine_pt()->geom_object_pt(1)->position(s_transition_lo,
409 spine_node_pt->spine_pt()->geom_object_pt(2)->position(s_transition_up,
411
413 // Horizontal position is always halfway between the walls
414 spine_centre[0] = 0.5 * (r_transition_lo[0] + r_transition_up[0]);
415 // Vertical spine centre is given by fraction between the walls
416 spine_centre[1] =
417 r_transition_lo[1] +
419
420 // Get vector towards origin
421 Vector<double> N(2);
422 N[0] = spine_centre[0] - r_wall_up[0];
423 N[1] = spine_centre[1] - r_wall_up[1];
424 double inv_length = 1.0 / sqrt(N[0] * N[0] + N[1] * N[1]);
425
426 // Update the nodal position
427 spine_node_pt->x(0) = r_wall_up[0] + w * h * N[0] * inv_length;
428 spine_node_pt->x(1) = r_wall_up[1] + w * h * N[1] * inv_length;
429 }
430
431
432 /// Update function for the deposited film region in the
433 /// upper part of the domain: Vertical spines
435 {
436 // Get fraction along the spine
437 double w = spine_node_pt->fraction();
438 // Get spine height
439 double h = spine_node_pt->h();
440
441 // Get wall coordinate
443 s_up[0] = spine_node_pt->spine_pt()->geom_parameter(0);
444
445 // Get position vector to wall
447 spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_up, r_wall_up);
448
449 // Update the nodal position
450 spine_node_pt->x(0) = r_wall_up[0];
451 spine_node_pt->x(1) = r_wall_up[1] - w * h;
452 }
453
454
455 /// Update function for the nodes in the channel region ahead
456 /// of the finger tip: Nodes are evenly distributed along vertical
457 /// lines between the top and bottom walls
459 {
460 // Get fraction along the spine
461 double w = spine_node_pt->fraction();
462
463 // Get upper and lower local coordinates
464 Vector<double> s_lo(1), s_up(1);
465 s_lo[0] = spine_node_pt->spine_pt()->geom_parameter(0);
466 s_up[0] = spine_node_pt->spine_pt()->geom_parameter(1);
467
468 // Get position vector to lower wall
469 Vector<double> r_lo(2), r_up(2);
470 spine_node_pt->spine_pt()->geom_object_pt(0)->position(s_lo, r_lo);
471 spine_node_pt->spine_pt()->geom_object_pt(1)->position(s_up, r_up);
472
473 // Update the nodal position
474 spine_node_pt->x(0) = r_lo[0] + w * (r_up[0] - r_lo[0]);
475 spine_node_pt->x(1) = r_lo[1] + w * (r_up[1] - r_lo[1]);
476 }
477
478
479 /// Initial reordering elements of the elements -- before
480 /// the channel mesh is added. Vertical stacks of elements, each topped by
481 /// their interface element -- this is (currently) identical to the
482 /// version in the SingleLayerSpineMesh but it's important
483 /// that the element reordering is maintained in exactly this form
484 /// for the rest of the mesh generation process to work properly.
485 /// Therefore we keep a copy of the function in here.
487
488 /// Number of elements along wall in deposited film region
489 unsigned Nx1;
490
491 /// Number of elements along wall in horizontal transition region
492 unsigned Nx2;
493
494 /// Number of elements along wall in channel region
495 unsigned Nx3;
496
497 /// Number of elements in vertical transition region (there are
498 /// twice as many elements across the channel)
499 unsigned Nhalf;
500
501 /// Number of elements across the deposited film
502 unsigned Nh;
503
504 /// Thickness of deposited film
505 double H;
506
507 /// Pointer to geometric object that represents the upper wall
509
510 /// Pointer to geometric object that represents the lower wall
512
513 /// Start coordinate on wall
515
516 /// Wall coordinate of end of liquid filled region (inflow)
517 double Zeta_end;
518
519 /// Wall coordinate of start of the transition region
521
522 /// Wall coordinate of end of transition region
524
525 /// Pointer to vertical fraction of the spine centre
527
528 /// Default spine fraction
530
531 /// Pointer to control element (just under the symmetry line
532 /// near the bubble tip; the bubble tip is located at s=[1.0,1.0]
533 /// in this element
535
536 /// Vector of pointers to element in the fluid layer
538
539 /// Vector of pointers to interface elements
541 };
542
543} // namespace oomph
544
546#endif
cstr elem_len * i
Definition cfortran.h:603
Mesh for 2D Bretherton problem – based on single layer mesh. Templated by spine-ified Navier-Stokes e...
FiniteElement *& interface_element_pt(const unsigned long &i)
Access functions for pointers to interface elements.
void spine_node_update_film_lower(SpineNode *spine_node_pt)
Update function for the deposited film region in the lower part of the domain: Vertical spines.
GeomObject * Lower_wall_pt
Pointer to geometric object that represents the lower wall.
void spine_node_update_horizontal_transition_lower(SpineNode *spine_node_pt)
Update function for the horizontal transitition region in the lower part of the domain: Spine points ...
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.
double Zeta_start
Start coordinate on wall.
double Zeta_transition_end
Wall coordinate of end of transition region.
unsigned Nhalf
Number of elements in vertical transition region (there are twice as many elements across the channel...
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 Nx2
Number of elements along wall in horizontal transition region.
unsigned Nx3
Number of elements along wall in channel region.
FiniteElement *& bulk_element_pt(const unsigned long &i)
Access functions for pointers to elements in bulk.
double Default_spine_centre_fraction
Default spine fraction.
void set_spine_centre_fraction_pt(double *const &fraction_pt)
Set the pointer to the spine centre's vertial fraction.
void spine_node_update_film_upper(SpineNode *spine_node_pt)
Update function for the deposited film region in the upper part of the domain: Vertical spines.
double H
Thickness of deposited film.
ELEMENT * control_element_pt()
Pointer to control element (just under the symmetry line near the bubble tip, so the bubble tip is lo...
unsigned long nbulk() const
Number of elements in bulk.
unsigned Nh
Number of elements across the deposited film.
void pin_all_spines()
Pin all spines so the mesh can be used for computation without free surfaces.
void initial_element_reorder()
Initial reordering elements of the elements – before the channel mesh is added. Vertical stacks of el...
void spine_node_update_channel(SpineNode *spine_node_pt)
Update function for the nodes in the channel region ahead of the finger tip: Nodes are evenly distrib...
double spine_centre_fraction() const
Read the value of the spine centre's vertical fraction.
void spine_node_update_vertical_transition_upper(SpineNode *spine_node_pt)
Update function for the vertical transitition region in the upper part of the domain: Spine points to...
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...
void spine_node_update(SpineNode *spine_node_pt)
General node update function implements pure virtual function defined in SpineMesh base class and per...
void spine_node_update_vertical_transition_lower(SpineNode *spine_node_pt)
Update function for the vertical transitition region in the lower part of the domain: Spine points to...
void spine_node_update_horizontal_transition_upper(SpineNode *spine_node_pt)
Update function for the horizontal transitition region in the upper part of the domain: Spine points ...
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.
unsigned nfree_surface_spines()
Number of free-surface spines (i.e. excluding the dummy spines in the channel region)
unsigned long ninterface_element() const
Number of elements on interface.
double * Spine_centre_fraction_pt
Pointer to vertical fraction of the spine centre.
unsigned Nx1
Number of elements along wall in deposited film region.
double Zeta_transition_start
Wall coordinate of start of the transition region.
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition nodes.h:385
A general Finite Element class.
Definition elements.h:1317
void position(const Vector< double > &zeta, Vector< double > &r) const
Return the parametrised position of the FiniteElement in its incarnation as a GeomObject,...
Definition elements.h:2680
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition elements.cc:4320
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:2222
A geometric object is an object that provides a parametrised description of its shape via the functio...
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors.
Definition mesh.h:75
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition mesh.h:477
An OomphLibError object which should be thrown when an run-time error is encountered....
Single-layer spine mesh class derived from standard 2D mesh. The mesh contains a layer of spinified f...
Spine *& spine_pt(const unsigned long &i)
Return the i-th spine in the mesh.
Definition spines.h:623
unsigned long nspine() const
Return the number of spines in the mesh.
Definition spines.h:635
Class for nodes that live on spines. The assumption is that each Node lies at a fixed fraction on a s...
Definition spines.h:328
Data *& spine_height_pt()
Access function to Data object that stores the spine height.
Definition spines.h:156
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).