channel_spine_mesh.template.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #ifndef OOMPH_CHANNEL_SPINE_MESH_TEMPLATE_CC
31 #define OOMPH_CHANNEL_SPINE_MESH_TEMPLATE_CC
32 
35 
36 
37 namespace oomph
38 {
39 
40 
41 
42 //===========================================================================
43 /// Constructor for spine 2D mesh: Pass number of elements in x-direction,
44 /// number of elements in y-direction, axial length and height of layer,
45 /// and pointer to timestepper (defaults to Static timestepper).
46 ///
47 /// The mesh contains a layer of spinified fluid elements (of type ELEMENT;
48 /// e.g SpineElement<QCrouzeixRaviartElement<2>)
49 /// and a surface layer of corresponding Spine interface elements
50 /// of type INTERFACE_ELEMENT, e.g.
51 /// SpineLineFluidInterfaceElement<ELEMENT> for 2D planar
52 /// problems.
53 //===========================================================================
54 template<class ELEMENT>
56  const unsigned &nx0,
57  const unsigned &nx1,
58  const unsigned &nx2,
59  const unsigned &ny,
60  const double &lx0,
61  const double &lx1,
62  const double &lx2,
63  const double &h,
64  GeomObject* wall_pt,
65  TimeStepper* time_stepper_pt) :
66  RectangularQuadMesh<ELEMENT >(nx0+nx1+nx2,ny,0.0,lx0+lx1+lx2,
67  0.0,h,false,false,time_stepper_pt),
68  Nx0(nx0), Nx1(nx1), Nx2(nx2),
69  Lx0(lx0), Lx1(lx1), Lx2(lx2),
70  Wall_pt(wall_pt)
71 {
72 
73  // Mesh can only be built with 2D Qelements.
74  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
75 
76  //Mesh can only be built with spine elements
77  MeshChecker::assert_geometric_element<SpineFiniteElement,ELEMENT>(2);
78 
79  // We've called the "generic" constructor for the RectangularQuadMesh
80  // which doesn't do much...
81 
82  // Build the straight line object
84 
85  // Now build the mesh:
86  build_channel_spine_mesh(time_stepper_pt);
87 
88 }
89 
90 
91 
92 //===========================================================================
93 /// Constuctor for spine 2D mesh: Pass number of elements in x-direction,
94 /// number of elements in y-direction, axial length and height of layer,
95 /// a boolean flag to make the mesh periodic in the x-direction,
96 /// and pointer to timestepper (defaults to Static timestepper).
97 ///
98 /// The mesh contains a layer of elements (of type ELEMENT)
99 /// and a surface layer of corresponding Spine interface elements (of type
100 /// SpineLineFluidInterfaceElement<ELEMENT>).
101 //===========================================================================
102 template<class ELEMENT>
104  const unsigned &nx0,
105  const unsigned &nx1,
106  const unsigned &nx2,
107  const unsigned &ny,
108  const double &lx0,
109  const double &lx1,
110  const double &lx2,
111  const double &h,
112  GeomObject* wall_pt,
113  const bool& periodic_in_x,
114  TimeStepper* time_stepper_pt) :
115  RectangularQuadMesh<ELEMENT >(nx0+nx1+nx2,ny,0.0,lx0+lx1+lx2,0.0,h,
116  periodic_in_x,false,time_stepper_pt),
117  Nx0(nx0), Nx1(nx1), Nx2(nx2),
118  Lx0(lx0), Lx1(lx1), Lx2(lx2),
119  Wall_pt(wall_pt)
120 {
121  // Mesh can only be built with 2D Qelements.
122  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
123 
124  //Mesh can only be built with spine elements
125  MeshChecker::assert_geometric_element<SpineFiniteElement,ELEMENT>(2);
126 
127  // We've called the "generic" constructor for the RectangularQuadMesh
128  // which doesn't do much...
129 
130  // Build the straight line object
132 
133  // Now build the mesh:
134  build_channel_spine_mesh(time_stepper_pt);
135 
136 }
137 
138 
139 
140 
141 
142 //===========================================================================
143 /// Helper function that actually builds the channel-spine mesh
144 /// based on the parameters set in the various constructors
145 //===========================================================================
146 template<class ELEMENT>
148  TimeStepper* time_stepper_pt)
149 {
150  // Build the underlying quad mesh:
152 
153  // Read out the number of elements in the x-direction and y-direction
154  // and in each of the left, centre and right regions
155  unsigned n_x = this->Nx;
156  unsigned n_y = this->Ny;
157  unsigned n_x0 = this->Nx0;
158  unsigned n_x1 = this->Nx1;
159  unsigned n_x2 = this->Nx2;
160 
161  // Set up the pointers to elements in the left region
162  unsigned nleft=n_x0*n_y;;
163  Left_element_pt.reserve(nleft);
164  unsigned ncentre=n_x1*n_y;;
165  Centre_element_pt.reserve(ncentre);
166  unsigned nright=n_x2*n_y;;
167  Right_element_pt.reserve(nright);
168  for(unsigned irow=0;irow<n_y;irow++)
169  {
170  for (unsigned e=0;e<n_x0;e++)
171  {
172  Left_element_pt.push_back(
173  this->finite_element_pt(irow*n_x+e));
174  }
175  for (unsigned e=0;e<n_x1;e++)
176  {
177  Centre_element_pt.push_back(
178  this->finite_element_pt(irow*n_x+(n_x0+e)));
179  }
180  for (unsigned e=0;e<n_x2;e++)
181  {
182  Right_element_pt.push_back(
183  this->finite_element_pt(irow*n_x+(n_x0+n_x1+e)));
184  }
185  }
186 
187 #ifdef PARANOID
188  // Check that we have the correct number of elements
189  if (nelement()!=nleft+ncentre+nright)
190  {
191  throw OomphLibError("Incorrect number of element pointers!",
192  OOMPH_CURRENT_FUNCTION,
193  OOMPH_EXCEPTION_LOCATION);
194  }
195 #endif
196 
197  // Allocate memory for the spines and fractions along spines
198  //---------------------------------------------------------
199 
200  // Read out number of linear points in the element
201  unsigned n_p = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
202 
203  unsigned nspine;
204  // Allocate store for the spines:
205  if (this->Xperiodic)
206  {
207  nspine = (n_p-1)*n_x;
208  Spine_pt.reserve(nspine);
209  // Number of spines in each region
210  // NOTE that boundary spines are in both regions
211  Nleft_spine = (n_p-1)*n_x0+1;
212  Ncentre_spine = (n_p-1)*n_x1+1;
213  Nright_spine = (n_p-1)*n_x2;
214  }
215  else
216  {
217  nspine = (n_p-1)*n_x+1;
218  Spine_pt.reserve(nspine);
219  // Number of spines in each region
220  // NOTE that boundary spines are in both regions
221  Nleft_spine = (n_p-1)*n_x0+1;
222  Ncentre_spine = (n_p-1)*n_x1+1;
223  Nright_spine = (n_p-1)*n_x2+1;
224  }
225 
226  // end Allocating memory
227 
228 
229  // set up the vectors of geometric data & objects for building spines
230  Vector<double> r_wall(2), zeta(1), s_wall(1);
231  GeomObject* geometric_object_pt=0;
232 
233  // LEFT REGION
234  // ===========
235 
236  // SPINES IN LEFT REGION
237  // ---------------------
238 
239  // Set up zeta increments
240  double zeta_lo=0.0;
241  double dzeta=Lx0/n_x0;
242 
243  // Initialise number of elements in previous regions:
244  unsigned n_prev_elements=0;
245 
246 
247  //FIRST SPINE
248  //-----------
249 
250  //Element 0
251  //Node 0
252  //Assign the new spine with unit length
253  Spine* new_spine_pt=new Spine(1.0);
254  new_spine_pt->spine_height_pt()->pin(0);
255  Spine_pt.push_back(new_spine_pt);
256 
257  // Get pointer to node
258  SpineNode* nod_pt=element_node_pt(0,0);
259  //Set the pointer to the spine
260  nod_pt->spine_pt() = new_spine_pt;
261  //Set the fraction
262  nod_pt->fraction() = 0.0;
263  // Pointer to the mesh that implements the update fct
264  nod_pt->spine_mesh_pt() = this;
265  // Set update fct id
266  nod_pt->node_update_fct_id()=0;
267 
268  // Provide spine with additional storage for wall coordinate
269  // and wall geom object:
270 
271  {
272  // Get the Lagrangian coordinate in the Lower Wall
273  zeta[0] = 0.0;
274  //Get the geometric object and local coordinate
275  Straight_wall_pt->locate_zeta(zeta,geometric_object_pt,s_wall);
276 
277  //The local coordinate is a geometric parameter
278  //This needs to be set (rather than added) because the
279  //same spine may be visited more than once
280  Vector<double> parameters(1,s_wall[0]);
281  nod_pt->spine_pt()->set_geom_parameter(parameters);
282 
283  // Get position of wall
284  Straight_wall_pt->position(s_wall,r_wall);
285 
286  // Adjust spine height
287  nod_pt->spine_pt()->height()=r_wall[1];
288 
289  // The sub geom object is one (and only) geom object
290  // for spine:
291  Vector<GeomObject*> geom_object_pt(1);
292  geom_object_pt[0] = geometric_object_pt;
293 
294  // Pass geom object(s) to spine
295  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
296  }
297 
298  //Loop vertically along the spine
299  //Loop over the elements
300  for(unsigned long i=0;i<n_y;i++)
301  {
302  //Loop over the vertical nodes, apart from the first
303  for(unsigned l1=1;l1<n_p;l1++)
304  {
305  // Get pointer to node
306  SpineNode* nod_pt=element_node_pt(i*n_x,l1*n_p);
307  //Set the pointer to the spine
308  nod_pt->spine_pt() = new_spine_pt;
309  //Set the fraction
310  nod_pt->fraction()=(double(i)+double(l1)/double(n_p-1))/double(n_y);
311  // Pointer to the mesh that implements the update fct
312  nod_pt->spine_mesh_pt() = this;
313  // Set update fct id
314  nod_pt->node_update_fct_id()=0;
315  }
316  } // end loop over elements
317 
318 
319  //LOOP OVER OTHER SPINES
320  //----------------------
321 
322  //Now loop over the elements horizontally
323  for(unsigned long j=0;j<n_x0;j++)
324  {
325  //Loop over the nodes in the elements horizontally, ignoring
326  //the first column
327  unsigned n_pmax=n_p;
328  for(unsigned l2=1;l2<n_pmax;l2++)
329  {
330  //Assign the new spine with unit height
331  new_spine_pt=new Spine(1.0);
332  new_spine_pt->spine_height_pt()->pin(0);
333  Spine_pt.push_back(new_spine_pt);
334 
335  // Get the node
336  SpineNode* nod_pt=element_node_pt(j,l2);
337  //Set the pointer to spine
338  nod_pt->spine_pt() = new_spine_pt;
339  //Set the fraction
340  nod_pt->fraction() = 0.0;
341  // Pointer to the mesh that implements the update fct
342  nod_pt->spine_mesh_pt() = this;
343  // Set update fct id
344  nod_pt->node_update_fct_id()=0;
345 
346  {
347  // Provide spine with additional storage for wall coordinate
348  // and wall geom object:
349 
350  // Get the Lagrangian coordinate in the Lower Wall
351  zeta[0] = zeta_lo + double(j)*dzeta + double(l2)*dzeta/(n_p-1.0);
352  //Get the geometric object and local coordinate
353  Straight_wall_pt->locate_zeta(zeta,geometric_object_pt,s_wall);
354 
355  //The local coordinate is a geometric parameter
356  //This needs to be set (rather than added) because the
357  //same spine may be visited more than once
358  Vector<double> parameters(1,s_wall[0]);
359  nod_pt->spine_pt()->set_geom_parameter(parameters);
360 
361  // Get position of wall
362  Straight_wall_pt->position(s_wall,r_wall);
363 
364  // Adjust spine height
365  nod_pt->spine_pt()->height()=r_wall[1];
366 
367  // The sub geom object is one (and only) geom object
368  // for spine:
369  Vector<GeomObject*> geom_object_pt(1);
370  geom_object_pt[0] = geometric_object_pt;
371 
372  // Pass geom object(s) to spine
373  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
374  }
375 
376  //Loop vertically along the spine
377  //Loop over the elements
378  for(unsigned long i=0;i<n_y;i++)
379  {
380  //Loop over the vertical nodes, apart from the first
381  for(unsigned l1=1;l1<n_p;l1++)
382  {
383  // Get the node
384  SpineNode* nod_pt=element_node_pt(i*n_x+j,l1*n_p+l2);
385  //Set the pointer to the spine
386  nod_pt->spine_pt() = new_spine_pt;
387  //Set the fraction
388  nod_pt->fraction()=(double(i)+double(l1)/double(n_p-1))/double(n_y);
389  // Pointer to the mesh that implements the update fct
390  nod_pt->spine_mesh_pt() = this;
391  // Set update fct id
392  nod_pt->node_update_fct_id()=0;
393  }
394  }
395  }
396  }
397 
398  // Increment number of previous elements
399  n_prev_elements+=n_x0*n_y;
400 
401 
402  // CENTRE REGION
403  // ===========
404 
405  zeta_lo=Lx0;
406  dzeta=Lx1/n_x1;
407 
408  // SPINES IN LEFT REGION
409  // ---------------------
410 
411  //LOOP OVER OTHER SPINES
412  //----------------------
413 
414  //Now loop over the elements horizontally
415  for(unsigned long j=n_x0;j<n_x0+n_x1;j++)
416  {
417  //Loop over the nodes in the elements horizontally, ignoring
418  //the first column
419  unsigned n_pmax=n_p;
420  for(unsigned l2=1;l2<n_pmax;l2++)
421  {
422  //Assign the new spine with unit height
423  new_spine_pt=new Spine(1.0);
424  new_spine_pt->spine_height_pt()->pin(0);
425  Spine_pt.push_back(new_spine_pt);
426 
427  // Get the node
428  SpineNode* nod_pt=element_node_pt(j,l2);
429  //Set the pointer to spine
430  nod_pt->spine_pt() = new_spine_pt;
431  //Set the fraction
432  nod_pt->fraction() = 0.0;
433  // Pointer to the mesh that implements the update fct
434  nod_pt->spine_mesh_pt() = this;
435  // Set update fct id
436  nod_pt->node_update_fct_id()=1;
437 
438  {
439  // Provide spine with additional storage for wall coordinate
440  // and wall geom object:
441 
442  // Get the Lagrangian coordinate in the Lower Wall
443  zeta[0] = zeta_lo + double(j-n_x0)*dzeta + double(l2)*dzeta/(n_p-1.0);
444  //Get the geometric object and local coordinate
445  Wall_pt->locate_zeta(zeta,geometric_object_pt,s_wall);
446 
447  //The local coordinate is a geometric parameter
448  //This needs to be set (rather than added) because the
449  //same spine may be visited more than once
450  Vector<double> parameters(1,s_wall[0]);
451  nod_pt->spine_pt()->set_geom_parameter(parameters);
452 
453  // Get position of wall
454  Wall_pt->position(s_wall,r_wall);
455 
456  // Adjust spine height
457  nod_pt->spine_pt()->height()=r_wall[1];
458 
459  // The sub geom object is one (and only) geom object
460  // for spine:
461  Vector<GeomObject*> geom_object_pt(1);
462  geom_object_pt[0] = geometric_object_pt;
463 
464  // Pass geom object(s) to spine
465  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
466  }
467 
468  //Loop vertically along the spine
469  //Loop over the elements
470  for(unsigned long i=0;i<n_y;i++)
471  {
472  //Loop over the vertical nodes, apart from the first
473  for(unsigned l1=1;l1<n_p;l1++)
474  {
475  // Get the node
476  SpineNode* nod_pt=element_node_pt(i*n_x+j,l1*n_p+l2);
477  //Set the pointer to the spine
478  nod_pt->spine_pt() = new_spine_pt;
479  //Set the fraction
480  nod_pt->fraction()=(double(i)+double(l1)/double(n_p-1))/double(n_y);
481  // Pointer to the mesh that implements the update fct
482  nod_pt->spine_mesh_pt() = this;
483  // Set update fct id
484  nod_pt->node_update_fct_id()=1;
485  }
486  }
487  }
488  }
489 
490  // Increment number of previous elements
491  n_prev_elements+=n_x1*n_y;
492 
493 
494  // RIGHT REGION
495  // ===========
496 
497  // SPINES IN RIGHT REGION
498  // ----------------------
499 
500  // Set up zeta increments
501  zeta_lo=Lx0+Lx1;
502  dzeta=Lx2/n_x2;
503 
504  //LOOP OVER OTHER SPINES
505  //----------------------
506 
507  //Now loop over the elements horizontally
508  for(unsigned long j=n_x0+n_x1;j<n_x0+n_x1+n_x2;j++)
509  {
510  //Need to copy last spine in previous element over to first spine
511  //in next elements, for all elements other than the first
512  if(j>0)
513  {
514  // For first spine, re-assign the geometric objects, since
515  // we treat it as part of the right region.
516  if (j==n_x0+n_x1)
517  {
518  SpineNode* nod_pt=element_node_pt(j,0);
519  // Set update fct id
520  nod_pt->node_update_fct_id()=2;
521  {
522  // Provide spine with additional storage for wall coordinate
523  // and wall geom object:
524 
525  // Get the Lagrangian coordinate in the Lower Wall
526  zeta[0] = zeta_lo + double(j-n_x0-n_x1)*dzeta;
527  //Get the geometric object and local coordinate
528  Straight_wall_pt->locate_zeta(zeta,geometric_object_pt,s_wall);
529 
530  //The local coordinate is a geometric parameter
531  //This needs to be set (rather than added) because the
532  //same spine may be visited more than once
533  Vector<double> parameters(1,s_wall[0]);
534  nod_pt->spine_pt()->set_geom_parameter(parameters);
535 
536  // Get position of wall
537  Straight_wall_pt->position(s_wall,r_wall);
538 
539  // Adjust spine height
540  nod_pt->spine_pt()->height()=r_wall[1];
541 
542  // The sub geom object is one (and only) geom object
543  // for spine:
544  Vector<GeomObject*> geom_object_pt(1);
545  geom_object_pt[0] = geometric_object_pt;
546 
547  // Pass geom object(s) to spine
548  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
549  }
550  }
551  }
552  //Loop over the nodes in the elements horizontally, ignoring
553  //the first column
554 
555  // Last spine needs special treatment in x-periodic meshes:
556  unsigned n_pmax=n_p;
557  if ((this->Xperiodic)&&(j==n_x-1)) n_pmax=n_p-1;
558 
559  for(unsigned l2=1;l2<n_pmax;l2++)
560  {
561  //Assign the new spine with unit height
562  new_spine_pt=new Spine(1.0);
563  new_spine_pt->spine_height_pt()->pin(0);
564  Spine_pt.push_back(new_spine_pt);
565 
566  // Get the node
567  SpineNode* nod_pt=element_node_pt(j,l2);
568  //Set the pointer to spine
569  nod_pt->spine_pt() = new_spine_pt;
570  //Set the fraction
571  nod_pt->fraction() = 0.0;
572  // Pointer to the mesh that implements the update fct
573  nod_pt->spine_mesh_pt() = this;
574  // Set update fct id
575  nod_pt->node_update_fct_id()=2;
576 
577  {
578  // Provide spine with additional storage for wall coordinate
579  // and wall geom object:
580 
581  // Get the Lagrangian coordinate in the Lower Wall
582  zeta[0] = zeta_lo + double(j-n_x0-n_x1)*dzeta + double(l2)*dzeta/(n_p-1.0);
583  //Get the geometric object and local coordinate
584  Straight_wall_pt->locate_zeta(zeta,geometric_object_pt,s_wall);
585 
586  //The local coordinate is a geometric parameter
587  //This needs to be set (rather than added) because the
588  //same spine may be visited more than once
589  Vector<double> parameters(1,s_wall[0]);
590  nod_pt->spine_pt()->set_geom_parameter(parameters);
591 
592  // Get position of wall
593  Straight_wall_pt->position(s_wall,r_wall);
594 
595  // Adjust spine height
596  nod_pt->spine_pt()->height()=r_wall[1];
597 
598  // The sub geom object is one (and only) geom object
599  // for spine:
600  Vector<GeomObject*> geom_object_pt(1);
601  geom_object_pt[0] = geometric_object_pt;
602 
603  // Pass geom object(s) to spine
604  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
605  }
606 
607  //Loop vertically along the spine
608  //Loop over the elements
609  for(unsigned long i=0;i<n_y;i++)
610  {
611  //Loop over the vertical nodes, apart from the first
612  for(unsigned l1=1;l1<n_p;l1++)
613  {
614  // Get the node
615  SpineNode* nod_pt=element_node_pt(i*n_x+j,l1*n_p+l2);
616  //Set the pointer to the spine
617  nod_pt->spine_pt() = new_spine_pt;
618  //Set the fraction
619  nod_pt->fraction()=(double(i)+double(l1)/double(n_p-1))/double(n_y);
620  // Pointer to the mesh that implements the update fct
621  nod_pt->spine_mesh_pt() = this;
622  // Set update fct id
623  nod_pt->node_update_fct_id()=2;
624  }
625  }
626  }
627  }
628 
629  // Increment number of previous elements
630  n_prev_elements+=n_x2*n_y;
631 
632  // Last spine needs special treatment for periodic meshes
633  // because it's the same as the first one...
634  if (this->Xperiodic)
635  {
636  // Last spine is the same as first one...
637  Spine* final_spine_pt=Spine_pt[0];
638 
639  // Get the node
640  SpineNode* nod_pt=element_node_pt((n_x-1),(n_p-1));
641 
642  //Set the pointer for the first node
643  nod_pt->spine_pt() = final_spine_pt;
644  //Set the fraction to be the same as for the nodes on the first row
645  nod_pt->fraction() = element_node_pt(0,0)->fraction();
646  // Pointer to the mesh that implements the update fct
647  nod_pt->spine_mesh_pt() = element_node_pt(0,0)->spine_mesh_pt();
648 
649  //Now loop vertically along the spine
650  for(unsigned i=0;i<n_y;i++)
651  {
652  //Loop over the vertical nodes, apart from the first
653  for(unsigned l1=1;l1<n_p;l1++)
654  {
655  // Get the node
656  SpineNode* nod_pt=element_node_pt(i*n_x+(n_x-1),l1*n_p+(n_p-1));
657 
658  //Set the pointer to the spine
659  nod_pt->spine_pt() = final_spine_pt;
660  //Set the fraction to be the same as in first row
661  nod_pt->fraction() = element_node_pt(i*n_x,l1*n_p)->fraction();
662  // Pointer to the mesh that implements the update fct
663  nod_pt->spine_mesh_pt() =
664  element_node_pt(i*n_x,l1*n_p)->spine_mesh_pt();
665  }
666  }
667 
668  }
669 
670 
671 
672 } // end of build_channel_spine_mesh
673 
674 
675 //======================================================================
676 /// Reorder the elements, so we loop over them vertically first
677 /// (advantageous in "wide" domains if a frontal solver is used).
678 //======================================================================
679 template <class ELEMENT>
681 {
682 
683  unsigned n_x = this->Nx;
684  unsigned n_y = this->Ny;
685  //Find out how many elements there are
686  unsigned long Nelement = nelement();
687  //Find out how many fluid elements there are
688  unsigned long Nfluid = n_x*n_y;
689  //Create a dummy array of elements
691 
692  //Loop over the elements in horizontal order
693  for(unsigned long j=0;j<n_x;j++)
694  {
695  //Loop over the elements in lower layer vertically
696  for(unsigned long i=0;i<n_y;i++)
697  {
698  //Push back onto the new stack
699  dummy.push_back(finite_element_pt(n_x*i + j));
700  }
701 
702  //Push back the line element onto the stack
703  dummy.push_back(finite_element_pt(Nfluid+j));
704  }
705 
706  //Now copy the reordered elements into the element_pt
707  for(unsigned long e=0;e<Nelement;e++)
708  {
709  Element_pt[e] = dummy[e];
710  }
711 
712 } // end of element_reorder
713 
714 
715 }
716 
717 #endif
void set_geom_parameter(const Vector< double > &geom_parameter)
Set vector of geometric parameters that are involved in the node update operations for this Spine...
Definition: spines.h:257
void element_reorder()
Reorder the elements so we loop over them vertically first (advantageous in "wide" domains if a front...
cstr elem_len * i
Definition: cfortran.h:607
double & height()
Access function to spine height.
Definition: spines.h:158
double & fraction()
Set reference to fraction along spine.
Definition: spines.h:350
ChannelSpineMesh(const unsigned &nx0, const unsigned &nx1, const unsigned &nx2, const unsigned &ny, const double &lx0, const double &lx1, const double &lx2, const double &h, GeomObject *wall_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in x-direction in regions 0,1 and 2, number of elements in y-dir...
unsigned & node_update_fct_id()
Access function to ID of node update function (within specific mesh)
Definition: spines.h:353
e
Definition: cfortran.h:575
virtual void build_channel_spine_mesh(TimeStepper *time_stepper_pt)
Helper function to actually build the channel-spine mesh (called from various constructors) ...
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
void build_mesh(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh construction function: contains all the hard work.
GeomObject * Straight_wall_pt
GeomObject for the straight upper wall.
Data *& spine_height_pt()
Access function to Data object that stores the spine height.
Definition: spines.h:164
Spine *& spine_pt()
Access function to spine.
Definition: spines.h:347
void set_geom_object_pt(const Vector< GeomObject * > &geom_object_pt)
Set vector of (pointers to) geometric objects that is involved in the node update operations for this...
Definition: spines.h:218
SpineMesh *& spine_mesh_pt()
Access function to Pointer to SpineMesh that this node is a part of and which implements the node upd...
Definition: spines.h:357
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219