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-2022 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
37namespace 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 {
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
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
288 zeta, lower_transition_geom_object_pt, s_transition_lo);
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
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
427 zeta, lower_sub_geom_object_pt, s_lo);
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
522 zeta, lower_sub_geom_object_pt, s_lo);
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
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
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
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
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
e
Definition: cfortran.h:571
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...
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.
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1271
void solve(DoubleVector &rhs)
Complete LU solve (replaces matrix by its LU decomposition and overwrites RHS with solution)....
Definition: matrices.cc:50
Class that is used to create FaceElement from bulk elements and to provide these FaceElement with a g...
A general Finite Element class.
Definition: elements.h:1313
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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:2218
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
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:381
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
A general mesh class.
Definition: mesh.h:67
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:243
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
Definition: mesh.h:833
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
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:407
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
virtual void setup_boundary_element_info()
Interface for function that is used to setup the boundary information (Empty virtual function – imple...
Definition: mesh.h:275
void remove_boundary_nodes()
Clear all pointers to boundary nodes.
Definition: mesh.cc:204
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
Definition: mesh.h:460
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:493
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:2394
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
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...
Vector< Spine * > Spine_pt
A Spine mesh contains a Vector of pointers to spines.
Definition: spines.h:616
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:669
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
double & fraction()
Set reference to fraction along spine.
Definition: spines.h:378
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:391
Spine *& spine_pt()
Access function to spine.
Definition: spines.h:372
unsigned & node_update_fct_id()
Access function to ID of node update function (within specific mesh)
Definition: spines.h:384
Spines are used for algebraic node update operations in free-surface fluid problems: They form the ba...
Definition: spines.h:64
Data *& spine_height_pt()
Access function to Data object that stores the spine height.
Definition: spines.h:156
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:225
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:273
double & height()
Access function to spine height.
Definition: spines.h:150
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...