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