adaptive_bubble_in_channel.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 
31 //Generic routines
32 #include "generic.h"
33 
34 
35 // The equations
36 #include "navier_stokes.h"
37 #include "solid.h"
38 #include "constitutive.h"
39 #include "fluid_interface.h"
40 
41 // The mesh
42 #include "meshes/triangle_mesh.h"
43 
44 using namespace std;
45 using namespace oomph;
46 
47 
48 namespace oomph
49 {
50 
51 
52 //==============================================================
53 /// Overload TaylorHood element to modify output
54 //==============================================================
56  public virtual PseudoSolidNodeUpdateElement<TTaylorHoodElement<2>,
57  TPVDElement<2,3> >
58  {
59 
60  private:
61 
62  /// Storage for elemental error estimate -- used for post-processing
63  double Error;
64 
65  public:
66 
67  /// Constructor initialise error
69  {
70  Error=0.0;
71  }
72 
73  /// Set error value for post-processing
74  void set_error(const double& error){Error=error;}
75 
76  /// Return variable identifier
77  std::string variable_identifier()
78  {
79  std::string txt="VARIABLES=";
80  txt+="\"x\",";
81  txt+="\"y\",";
82  txt+="\"u\",";
83  txt+="\"v\",";
84  txt+="\"p\",";
85  txt+="\"du/dt\",";
86  txt+="\"dv/dt\",";
87  txt+="\"u_m\",";
88  txt+="\"v_m\",";
89  txt+="\"x_h1\",";
90  txt+="\"y_h1\",";
91  txt+="\"x_h2\",";
92  txt+="\"y_h2\",";
93  txt+="\"u_h1\",";
94  txt+="\"v_h1\",";
95  txt+="\"u_h2\",";
96  txt+="\"v_h2\",";
97  txt+="\"du/dx\",";
98  txt+="\"du/dy\",";
99  txt+="\"dv/dx\",";
100  txt+="\"dv/dy\",";
101  txt+="\"error\",";
102  txt+="\"size\",";
103  txt+="\n";
104  return txt;
105  }
106 
107 
108  /// Overload output function
109  void output(std::ostream &outfile,
110  const unsigned &nplot)
111  {
112 
113  // Assign dimension
114  unsigned el_dim=2;
115 
116  // Vector of local coordinates
117  Vector<double> s(el_dim);
118 
119  // Acceleration
120  Vector<double> dudt(el_dim);
121 
122  // Mesh elocity
123  Vector<double> mesh_veloc(el_dim,0.0);
124 
125  // Tecplot header info
126  outfile << tecplot_zone_string(nplot);
127 
128  // Find out how many nodes there are
129  unsigned n_node = nnode();
130 
131  //Set up memory for the shape functions
132  Shape psif(n_node);
133  DShape dpsifdx(n_node,el_dim);
134 
135  // Loop over plot points
136  unsigned num_plot_points=nplot_points(nplot);
137  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
138  {
139 
140  // Get local coordinates of plot point
141  get_s_plot(iplot,nplot,s);
142 
143  //Call the derivatives of the shape and test functions
144  dshape_eulerian(s,psif,dpsifdx);
145 
146  //Allocate storage
147  Vector<double> mesh_veloc(el_dim);
148  Vector<double> dudt(el_dim);
149  Vector<double> dudt_ALE(el_dim);
150  DenseMatrix<double> interpolated_dudx(el_dim,el_dim);
151 
152  //Initialise everything to zero
153  for(unsigned i=0;i<el_dim;i++)
154  {
155  mesh_veloc[i]=0.0;
156  dudt[i]=0.0;
157  dudt_ALE[i]=0.0;
158  for(unsigned j=0;j<el_dim;j++)
159  {
160  interpolated_dudx(i,j) = 0.0;
161  }
162  }
163 
164  //Calculate velocities and derivatives
165 
166  //Loop over directions
167  for(unsigned i=0;i<el_dim;i++)
168  {
169  //Get the index at which velocity i is stored
170  unsigned u_nodal_index = u_index_nst(i);
171  // Loop over nodes
172  for(unsigned l=0;l<n_node;l++)
173  {
174  dudt[i]+=du_dt_nst(l,u_nodal_index)*psif[l];
175  mesh_veloc[i]+=dnodal_position_dt(l,i)*psif[l];
176 
177  //Loop over derivative directions for velocity gradients
178  for(unsigned j=0;j<el_dim;j++)
179  {
180  interpolated_dudx(i,j) += nodal_value(l,u_nodal_index)*
181  dpsifdx(l,j);
182  }
183  }
184  }
185 
186 
187  // Get dudt in ALE form (incl mesh veloc)
188  for(unsigned i=0;i<el_dim;i++)
189  {
190  dudt_ALE[i]=dudt[i];
191  for (unsigned k=0;k<el_dim;k++)
192  {
193  dudt_ALE[i]-=mesh_veloc[k]*interpolated_dudx(i,k);
194  }
195  }
196 
197 
198  // Coordinates
199  for(unsigned i=0;i<el_dim;i++)
200  {
201  outfile << interpolated_x(s,i) << " ";
202  }
203 
204  // Velocities
205  for(unsigned i=0;i<el_dim;i++)
206  {
207  outfile << interpolated_u_nst(s,i) << " ";
208  }
209 
210  // Pressure
211  outfile << interpolated_p_nst(s) << " ";
212 
213  // Accelerations
214  for(unsigned i=0;i<el_dim;i++)
215  {
216  outfile << dudt_ALE[i] << " ";
217  }
218 
219  // Mesh velocity
220  for(unsigned i=0;i<el_dim;i++)
221  {
222  outfile << mesh_veloc[i] << " ";
223  }
224 
225  // History values of coordinates
226  unsigned n_prev=node_pt(0)->position_time_stepper_pt()->ntstorage();
227  for (unsigned t=1;t<n_prev;t++)
228  {
229  for(unsigned i=0;i<el_dim;i++)
230  {
231  outfile << interpolated_x(t,s,i) << " ";
232  }
233  }
234 
235  // History values of velocities
236  n_prev=node_pt(0)->time_stepper_pt()->ntstorage();
237  for (unsigned t=1;t<n_prev;t++)
238  {
239  for(unsigned i=0;i<el_dim;i++)
240  {
241  outfile << interpolated_u_nst(t,s,i) << " ";
242  }
243  }
244 
245  // Velocity gradients
246  for(unsigned i=0;i<el_dim;i++)
247  {
248  for(unsigned j=0;j<el_dim;j++)
249  {
250  outfile << interpolated_dudx(i,j) << " ";
251  }
252  }
253 
254  outfile << Error << " "
255  << size() << std::endl;
256  }
257 
258  // Write tecplot footer (e.g. FE connectivity lists)
259  write_tecplot_zone_footer(outfile,nplot);
260  }
261 
262 
263  /// Get square of L2 norm of velocity
265  {
266 
267  // Assign dimension
268  unsigned el_dim=2;
269  // Initalise
270  double sum=0.0;
271 
272  //Find out how many nodes there are
273  unsigned n_node = nnode();
274 
275  //Find the indices at which the local velocities are stored
276  unsigned u_nodal_index[el_dim];
277  for(unsigned i=0;i<el_dim;i++) {u_nodal_index[i] = u_index_nst(i);}
278 
279  //Set up memory for the velocity shape fcts
280  Shape psif(n_node);
281  DShape dpsidx(n_node,el_dim);
282 
283  //Number of integration points
284  unsigned n_intpt = integral_pt()->nweight();
285 
286  //Set the Vector to hold local coordinates
287  Vector<double> s(el_dim);
288 
289  //Loop over the integration points
290  for(unsigned ipt=0;ipt<n_intpt;ipt++)
291  {
292  //Assign values of s
293  for(unsigned i=0;i<el_dim;i++) s[i] = integral_pt()->knot(ipt,i);
294 
295  //Get the integral weight
296  double w = integral_pt()->weight(ipt);
297 
298  // Call the derivatives of the veloc shape functions
299  // (Derivs not needed but they are free)
300  double J = this->dshape_eulerian_at_knot(ipt,psif,dpsidx);
301 
302  //Premultiply the weights and the Jacobian
303  double W = w*J;
304 
305  //Calculate velocities
306  Vector<double> interpolated_u(el_dim,0.0);
307 
308  // Loop over nodes
309  for(unsigned l=0;l<n_node;l++)
310  {
311  //Loop over directions
312  for(unsigned i=0;i<el_dim;i++)
313  {
314  //Get the nodal value
315  double u_value = raw_nodal_value(l,u_nodal_index[i]);
316  interpolated_u[i] += u_value*psif[l];
317  }
318  }
319 
320  //Assemble square of L2 norm
321  for(unsigned i=0;i<el_dim;i++)
322  {
323  sum+=interpolated_u[i]*interpolated_u[i]*W;
324  }
325  }
326 
327  return sum;
328 
329  }
330 
331  };
332 
333 
334 //=======================================================================
335 /// Face geometry for element is the same as that for the underlying
336 /// wrapped element
337 //=======================================================================
338  template<>
339  class FaceGeometry<MyTaylorHoodElement>
340  : public virtual SolidTElement<1,3>
341  {
342  public:
343  FaceGeometry() : SolidTElement<1,3>() {}
344  };
345 
346 
347 //=======================================================================
348 /// Face geometry of Face geometry for element is the same
349 /// as that for the underlying
350 /// wrapped element
351 //=======================================================================
352  template<>
353  class FaceGeometry<FaceGeometry<MyTaylorHoodElement> >
354  : public virtual SolidPointElement
355  {
356  public:
357  FaceGeometry() : SolidPointElement() {}
358  };
359 
360 
361 } //End of namespace extension
362 
363 
364 
365 /////////////////////////////////////////////////////////////////
366 /////////////////////////////////////////////////////////////////
367 /////////////////////////////////////////////////////////////////
368 
369 
370 //==start_of_namespace==============================
371 /// Namespace for Problem Parameter
372 //==================================================
373  namespace Problem_Parameter
374  {
375  /// Doc info object
376  DocInfo Doc_info;
377 
378  /// Reynolds number
379  double Re=0.0;
380 
381  /// Capillary number
382  double Ca = 10.0;
383 
384  /// Pseudo-solid Poisson ratio
385  double Nu=0.3;
386 
387  /// Initial radius of bubble
388  double Radius = 0.25;
389 
390  /// \short Volume of the bubble (negative because it's outside the
391  /// fluid!)
392  double Volume = -MathematicalConstants::Pi*Radius*Radius;
393 
394  /// \short Scaling factor for inflow velocity (allows it to be switched off
395  /// to do hydrostatics)
397 
398  /// \short Length of the channel
399  double Length = 3.0;
400 
401  /// Constitutive law used to determine the mesh deformation
402  ConstitutiveLaw *Constitutive_law_pt=0;
403 
404  /// Trace file
405  ofstream Trace_file;
406 
407  /// \short File to document the norm of the solution (for validation
408  /// purposes -- triangle doesn't give fully reproducible results so
409  /// mesh generation/adaptation may generate slightly different numbers
410  /// of elements on different machines!)
411  ofstream Norm_file;
412 
413  } // end_of_namespace
414 
415 
416 
417 
418 ///////////////////////////////////////////////////////////
419 ///////////////////////////////////////////////////////////
420 ///////////////////////////////////////////////////////////
421 
422 
423 
424 //==start_of_problem_class============================================
425 /// Problem class to simulate inviscid bubble propagating along 2D channel
426 //====================================================================
427 template<class ELEMENT>
428 class BubbleInChannelProblem : public Problem
429 {
430 
431 public:
432 
433  /// Constructor
435 
436  /// Destructor
438  {
439  // Fluid timestepper
440  delete this->time_stepper_pt(0);
441 
442  // Kill data associated with outer boundary
443  unsigned n=Outer_boundary_polyline_pt->npolyline();
444  for (unsigned j=0;j<n;j++)
445  {
446  delete Outer_boundary_polyline_pt->polyline_pt(j);
447  }
448  delete Outer_boundary_polyline_pt;
449 
450  //Kill data associated with bubbles
451  unsigned n_bubble = Bubble_polygon_pt.size();
452  for(unsigned ibubble=0;ibubble<n_bubble;ibubble++)
453  {
454  unsigned n=Bubble_polygon_pt[ibubble]->npolyline();
455  for (unsigned j=0;j<n;j++)
456  {
457  delete Bubble_polygon_pt[ibubble]->polyline_pt(j);
458  }
459  delete Bubble_polygon_pt[ibubble];
460  }
461 
462  // Flush element of free surface elements
463  delete_free_surface_elements();
464  delete Free_surface_mesh_pt;
465  delete_volume_constraint_elements();
466  delete Volume_constraint_mesh_pt;
467 
468  // Delete error estimator
469  delete Fluid_mesh_pt->spatial_error_estimator_pt();
470 
471  // Delete fluid mesh
472  delete Fluid_mesh_pt;
473 
474  // Delete the global pressure bubble data
475  delete Bubble_pressure_data_pt;
476 
477  // Kill const eqn
479 
480  }
481 
482 
483  /// Actions before adapt: Wipe the mesh of free surface elements
485  {
486  // Kill the elements and wipe surface mesh
487  delete_free_surface_elements();
488  delete_volume_constraint_elements();
489 
490  // Rebuild the Problem's global mesh from its various sub-meshes
491  this->rebuild_global_mesh();
492 
493  }// end of actions_before_adapt
494 
495 
496  /// Actions after adapt: Rebuild the mesh of free surface elements
498  {
499  // Create the elements that impose the displacement constraint
500  create_free_surface_elements();
501  create_volume_constraint_elements();
502 
503  // Rebuild the Problem's global mesh from its various sub-meshes
504  this->rebuild_global_mesh();
505 
506  // Setup the problem again -- remember that fluid mesh has been
507  // completely rebuilt and its element's don't have any
508  // pointers to Re etc. yet
509  complete_problem_setup();
510 
511  }// end of actions_after_adapt
512 
513 
514  /// Update the after solve (empty)
516 
517  /// Update the problem specs before solve
519  {
520  //Reset the Lagrangian coordinates of the nodes to be the current
521  //Eulerian coordinates -- this makes the current configuration
522  //stress free
523  Fluid_mesh_pt->set_lagrangian_nodal_coordinates();
524  }
525 
526  /// \short Set boundary conditions and complete the build of all elements
527  void complete_problem_setup();
528 
529  /// Doc the solution
530  void doc_solution(const std::string& comment="");
531 
532  /// Compute the error estimates and assign to elements for plotting
533  void compute_error_estimate(double& max_err,
534  double& min_err);
535 
536 private:
537 
538 
539  /// \short Create free surface elements
540  void create_free_surface_elements();
541 
542  /// \short Delete free surface elements
544  {
545  // How many surface elements are in the surface mesh
546  unsigned n_element = Free_surface_mesh_pt->nelement();
547 
548  // Loop over the surface elements
549  for(unsigned e=0;e<n_element;e++)
550  {
551  // Kill surface element
552  delete Free_surface_mesh_pt->element_pt(e);
553  }
554 
555  // Wipe the mesh
556  Free_surface_mesh_pt->flush_element_and_node_storage();
557 
558  } // end of delete_free_surface_elements
559 
560 
561 /// Create elements that impose volume constraint on the bubble
562  void create_volume_constraint_elements();
563 
564  /// \short Delete volume constraint elements
566  {
567  // How many surface elements are in the surface mesh
568  unsigned n_element = Volume_constraint_mesh_pt->nelement();
569 
570  // Loop over the surface elements (but don't kill the volume constraint
571  // element (element 0))
572  unsigned first_el_to_be_killed=1;
573  for(unsigned e=first_el_to_be_killed;e<n_element;e++)
574  {
575  delete Volume_constraint_mesh_pt->element_pt(e);
576  }
577 
578  // Wipe the mesh
579  Volume_constraint_mesh_pt->flush_element_and_node_storage();
580 
581  } // end of delete_volume_constraint_elements
582 
583  /// Pointers to mesh of free surface elements
585 
586  /// Pointer to mesh containing elements that impose volume constraint
588 
589  /// Pointer to Fluid_mesh
590  RefineableSolidTriangleMesh<ELEMENT>* Fluid_mesh_pt;
591 
592  /// Vector storing pointer to the bubble polygons
593  Vector<TriangleMeshPolygon*> Bubble_polygon_pt;
594 
595  /// Triangle mesh polygon for outer boundary
596  TriangleMeshPolygon* Outer_boundary_polyline_pt;
597 
598  /// Pointer to a global bubble pressure datum
600 
601  /// Pointer to element that imposes volume constraint for bubble
602  VolumeConstraintElement* Vol_constraint_el_pt;
603 
604  /// Enumeration of mesh boundaries
605  enum
606  {
607  Inflow_boundary_id=0,
608  Upper_wall_boundary_id=1,
609  Outflow_boundary_id=2,
610  Bottom_wall_boundary_id=3,
611  First_bubble_boundary_id=4,
612  Second_bubble_boundary_id=5
613  };
614 
615 
616 }; // end_of_problem_class
617 
618 
619 //==start_constructor=====================================================
620 /// Constructor
621 //========================================================================
622 template<class ELEMENT>
624 {
625  // Output directory
626  Problem_Parameter::Doc_info.set_directory("RESLT");
627 
628  // Allocate the timestepper -- this constructs the Problem's
629  // time object with a sufficient amount of storage to store the
630  // previous timsteps.
631  this->add_time_stepper_pt(new BDF<2>);
632 
633 
634 
635  // Build volume constraint element: Pass pointer to double that
636  // specifies target volume, data that contains the "traded" pressure
637  // and the index of the traded pressure value within this Data item
638 
639 #ifdef GLOBAL_DATA
640 
641  // Create bubble pressure as global Data
642  Bubble_pressure_data_pt = new Data(1);
643  unsigned index_of_traded_pressure=0;
644  this->add_global_data(Bubble_pressure_data_pt);
645 
646 
647  Vol_constraint_el_pt=
648  new VolumeConstraintElement(&Problem_Parameter::Volume,
649  Bubble_pressure_data_pt,
650  index_of_traded_pressure);
651 
652  //Provide a reasonable initial guess for bubble pressure (hydrostatics):
653  Bubble_pressure_data_pt->set_value(index_of_traded_pressure,
656 #else
657 
658  // Build element and create pressure internally
659  Vol_constraint_el_pt=
660  new VolumeConstraintElement(&Problem_Parameter::Volume);
661 
662  // Which value stores the pressure?
663  unsigned index=Vol_constraint_el_pt->index_of_traded_pressure();
664 
665  // Pressure data
666  Bubble_pressure_data_pt=Vol_constraint_el_pt->p_traded_data_pt();
667 
668  // Assign initial value
669  Vol_constraint_el_pt->p_traded_data_pt()->
671 
672 #endif
673 
674 
675  // Build the boundary segments for outer boundary, consisting of
676  //--------------------------------------------------------------
677  // four separate polylines
678  //------------------------
679  Vector<TriangleMeshCurveSection*> boundary_polyline_pt(4);
680 
681  // Each polyline only has two vertices -- provide storage for their
682  // coordinates
683  Vector<Vector<double> > vertex_coord(2);
684  for(unsigned i=0;i<2;i++)
685  {
686  vertex_coord[i].resize(2);
687  }
688 
689  // First polyline: Inflow
690  vertex_coord[0][0]=0.0;
691  vertex_coord[0][1]=0.0;
692  vertex_coord[1][0]=0.0;
693  vertex_coord[1][1]=1.0;
694 
695  // Build the 1st boundary polyline
696  boundary_polyline_pt[0] = new TriangleMeshPolyLine(vertex_coord,
697  Inflow_boundary_id);
698 
699  // Second boundary polyline: Upper wall
700  vertex_coord[0][0]=0.0;
701  vertex_coord[0][1]=1.0;
702  vertex_coord[1][0]=Problem_Parameter::Length;
703  vertex_coord[1][1]=1.0;
704 
705  // Build the 2nd boundary polyline
706  boundary_polyline_pt[1] = new TriangleMeshPolyLine(vertex_coord,
707  Upper_wall_boundary_id);
708 
709  // Third boundary polyline: Outflow
710  vertex_coord[0][0]=Problem_Parameter::Length;
711  vertex_coord[0][1]=1.0;
712  vertex_coord[1][0]=Problem_Parameter::Length;
713  vertex_coord[1][1]=0.0;
714 
715  // Build the 3rd boundary polyline
716  boundary_polyline_pt[2] = new TriangleMeshPolyLine(vertex_coord,
717  Outflow_boundary_id);
718 
719  // Fourth boundary polyline: Bottom wall
720  vertex_coord[0][0]=Problem_Parameter::Length;
721  vertex_coord[0][1]=0.0;
722  vertex_coord[1][0]=0.0;
723  vertex_coord[1][1]=0.0;
724 
725  // Build the 4th boundary polyline
726  boundary_polyline_pt[3] = new TriangleMeshPolyLine(vertex_coord,
727  Bottom_wall_boundary_id);
728 
729  // Create the triangle mesh polygon for outer boundary
730  Outer_boundary_polyline_pt = new TriangleMeshPolygon(boundary_polyline_pt);
731 
732 
733  // Now define initial shape of bubble(s) with polygon
734  //---------------------------------------------------
735 
736  // We have one bubble
737  Bubble_polygon_pt.resize(1);
738 
739  // Place it smack in the middle of the channel
740  double x_center = 0.5*Problem_Parameter::Length;
741  double y_center = 0.5;
742  Ellipse * bubble_pt = new Ellipse(Problem_Parameter::Radius,
744 
745  // Intrinsic coordinate along GeomObject defining the bubble
746  Vector<double> zeta(1);
747 
748  // Position vector to GeomObject defining the bubble
749  Vector<double> coord(2);
750 
751  // Number of points defining bubble
752  unsigned npoints = 16;
753  double unit_zeta = MathematicalConstants::Pi/double(npoints-1);
754 
755  // This bubble is bounded by two distinct boundaries, each
756  // represented by its own polyline
757  Vector<TriangleMeshCurveSection*> bubble_polyline_pt(2);
758 
759  // Vertex coordinates
760  Vector<Vector<double> > bubble_vertex(npoints);
761 
762  // Create points on boundary
763  for(unsigned ipoint=0; ipoint<npoints;ipoint++)
764  {
765  // Resize the vector
766  bubble_vertex[ipoint].resize(2);
767 
768  // Get the coordinates
769  zeta[0]=unit_zeta*double(ipoint);
770  bubble_pt->position(zeta,coord);
771 
772  // Shift
773  bubble_vertex[ipoint][0]=coord[0]+x_center;
774  bubble_vertex[ipoint][1]=coord[1]+y_center;
775  }
776 
777  // Build the 1st bubble polyline
778  bubble_polyline_pt[0] = new TriangleMeshPolyLine(bubble_vertex,
779  First_bubble_boundary_id);
780 
781  // Second boundary of bubble
782  for(unsigned ipoint=0; ipoint<npoints;ipoint++)
783  {
784  // Resize the vector
785  bubble_vertex[ipoint].resize(2);
786 
787  // Get the coordinates
788  zeta[0]=(unit_zeta*double(ipoint))+MathematicalConstants::Pi;
789  bubble_pt->position(zeta,coord);
790 
791  // Shift
792  bubble_vertex[ipoint][0]=coord[0]+x_center;
793  bubble_vertex[ipoint][1]=coord[1]+y_center;
794  }
795 
796  // Build the 2nd bubble polyline
797  bubble_polyline_pt[1] = new TriangleMeshPolyLine(bubble_vertex,
798  Second_bubble_boundary_id);
799 
800 
801  // Define coordinates of a point inside the bubble
802  Vector<double> bubble_center(2);
803  bubble_center[0]=x_center;
804  bubble_center[1]=y_center;
805 
806 
807  // Create closed polygon from two polylines
808  Bubble_polygon_pt[0] = new TriangleMeshPolygon(
809  bubble_polyline_pt,
810  bubble_center);
811 
812  // Now build the mesh, based on the boundaries specified by
813  //---------------------------------------------------------
814  // polygons just created
815  //----------------------
816 
817  // Convert to "closed curve" objects
818  TriangleMeshClosedCurve* outer_closed_curve_pt=Outer_boundary_polyline_pt;
819  unsigned nb=Bubble_polygon_pt.size();
820  Vector<TriangleMeshClosedCurve*> bubble_closed_curve_pt(nb);
821  for (unsigned i=0;i<nb;i++)
822  {
823  bubble_closed_curve_pt[i]=Bubble_polygon_pt[i];
824  }
825 
826  // Target area for initial mesh
827  double uniform_element_area=0.2;
828 
829  // Use the TriangleMeshParameters object for gathering all
830  // the necessary arguments for the TriangleMesh object
831  TriangleMeshParameters triangle_mesh_parameters(
832  outer_closed_curve_pt);
833 
834  // Define the holes on the boundary
835  triangle_mesh_parameters.internal_closed_curve_pt() =
836  bubble_closed_curve_pt;
837 
838  // Define the maximum element areas
839  triangle_mesh_parameters.element_area() =
840  uniform_element_area;
841 
842  // Create the mesh
843  Fluid_mesh_pt =
844  new RefineableSolidTriangleMesh<ELEMENT>(
845  triangle_mesh_parameters, this->time_stepper_pt());
846 
847  // Set error estimator for bulk mesh
848  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
849  Fluid_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
850 
851  // Set targets for spatial adaptivity
852  Fluid_mesh_pt->max_permitted_error()=0.005;
853  Fluid_mesh_pt->min_permitted_error()=0.001;
854  Fluid_mesh_pt->max_element_size()=0.2;
855  Fluid_mesh_pt->min_element_size()=0.001;
856 
857  // Use coarser mesh during validation
858  if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
859  {
860  Fluid_mesh_pt->min_element_size()=0.01;
861  }
862 
863  // Output boundary and mesh initial mesh for information
864  this->Fluid_mesh_pt->output_boundaries("boundaries.dat");
865  this->Fluid_mesh_pt->output("mesh.dat");
866 
867  // Set boundary condition and complete the build of all elements
868  complete_problem_setup();
869 
870  // Construct the mesh of free surface elements
871  Free_surface_mesh_pt=new Mesh;
872  create_free_surface_elements();
873 
874  // Construct the mesh of elements that impose the volume constraint
875  Volume_constraint_mesh_pt = new Mesh;
876  create_volume_constraint_elements();
877 
878  // Combine meshes
879  //---------------
880 
881  // Add volume constraint sub mesh
882  this->add_sub_mesh(this->Volume_constraint_mesh_pt);
883 
884  // Add Fluid_mesh_pt sub meshes
885  this->add_sub_mesh(Fluid_mesh_pt);
886 
887  // Add Free_surface sub meshes
888  this->add_sub_mesh(this->Free_surface_mesh_pt);
889 
890  // Build global mesh
891  this->build_global_mesh();
892 
893  // Setup equation numbering scheme
894  cout <<"Number of equations: " << this->assign_eqn_numbers() << std::endl;
895 
896 } // end_of_constructor
897 
898 
899 //============start_of_create_free_surface_elements======================
900 /// Create elements that impose the kinematic and dynamic bcs
901 /// for the pseudo-solid fluid mesh
902 //=======================================================================
903 template<class ELEMENT>
905 {
906 
907  // Volume constraint element stores the Data item that stores
908  // the bubble pressure that is adjusted/traded to allow for
909  // volume conservation. Which value is the pressure stored in?
910  unsigned p_traded_index=Vol_constraint_el_pt->index_of_traded_pressure();
911 
912  //Loop over the free surface boundaries
913  unsigned nb=Fluid_mesh_pt->nboundary();
914  for(unsigned b=First_bubble_boundary_id;b<nb;b++)
915  {
916  // How many bulk fluid elements are adjacent to boundary b?
917  unsigned n_element = Fluid_mesh_pt->nboundary_element(b);
918 
919  // Loop over the bulk fluid elements adjacent to boundary b?
920  for(unsigned e=0;e<n_element;e++)
921  {
922  // Get pointer to the bulk fluid element that is
923  // adjacent to boundary b
924  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
925  Fluid_mesh_pt->boundary_element_pt(b,e));
926 
927  //Find the index of the face of element e along boundary b
928  int face_index = Fluid_mesh_pt->face_index_at_boundary(b,e);
929 
930  // Create new element
931  ElasticLineFluidInterfaceElement<ELEMENT>* el_pt =
932  new ElasticLineFluidInterfaceElement<ELEMENT>(
933  bulk_elem_pt,face_index);
934 
935  // Add it to the mesh
936  Free_surface_mesh_pt->add_element_pt(el_pt);
937 
938  //Add the appropriate boundary number
939  el_pt->set_boundary_number_in_bulk_mesh(b);
940 
941  //Specify the capillary number
942  el_pt->ca_pt() = &Problem_Parameter::Ca;
943 
944  // Specify the bubble pressure (pointer to Data object and
945  // index of value within that Data object that corresponds
946  // to the traded pressure
947  el_pt->set_external_pressure_data(
948  Vol_constraint_el_pt->p_traded_data_pt(),p_traded_index);
949  }
950  }
951 }
952 // end of create_free_surface_elements
953 
954 
955 
956 
957 
958 //============start_of_create_volume_constraint_elements=================
959 /// Create elements that impose volume constraint on the bubble
960 //=======================================================================
961 template<class ELEMENT>
963 {
964 
965  // Add volume constraint element to the mesh
966  Volume_constraint_mesh_pt->add_element_pt(Vol_constraint_el_pt);
967 
968  //Loop over the free surface boundaries
969  unsigned nb=Fluid_mesh_pt->nboundary();
970  for(unsigned b=First_bubble_boundary_id;b<nb;b++)
971  {
972  // How many bulk fluid elements are adjacent to boundary b?
973  unsigned n_element = Fluid_mesh_pt->nboundary_element(b);
974 
975  // Loop over the bulk fluid elements adjacent to boundary b?
976  for(unsigned e=0;e<n_element;e++)
977  {
978  // Get pointer to the bulk fluid element that is
979  // adjacent to boundary b
980  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
981  Fluid_mesh_pt->boundary_element_pt(b,e));
982 
983  //Find the index of the face of element e along boundary b
984  int face_index = Fluid_mesh_pt->face_index_at_boundary(b,e);
985 
986  // Create new element
987  ElasticLineVolumeConstraintBoundingElement<ELEMENT>* el_pt =
988  new ElasticLineVolumeConstraintBoundingElement<ELEMENT>(
989  bulk_elem_pt,face_index);
990 
991  //Set the "master" volume constraint element
992  el_pt->set_volume_constraint_element(Vol_constraint_el_pt);
993 
994  // Add it to the mesh
995  Volume_constraint_mesh_pt->add_element_pt(el_pt);
996  }
997  }
998 }
999 // end of create_volume_constraint_elements
1000 
1001 
1002 
1003 //==start_of_complete_problem_setup=======================================
1004 /// Set boundary conditions and complete the build of all elements
1005 //========================================================================
1006 template<class ELEMENT>
1008 {
1009  // Map to record if a given boundary is on a bubble or not
1010  map<unsigned,bool> is_on_bubble_bound;
1011 
1012  // Loop over the bubbles
1013  unsigned nbubble=Bubble_polygon_pt.size();
1014  for(unsigned ibubble=0;ibubble<nbubble;ibubble++)
1015  {
1016  // Get the vector all boundary IDs associated with the polylines that
1017  // make up the closed polygon
1018  Vector<unsigned> bubble_bound_id=this->Bubble_polygon_pt[ibubble]->
1019  polygon_boundary_id();
1020 
1021  // Get the number of boundary
1022  unsigned nbound=bubble_bound_id.size();
1023 
1024  // Fill in the map
1025  for(unsigned ibound=0;ibound<nbound;ibound++)
1026  {
1027  // This boundary...
1028  unsigned bound_id=bubble_bound_id[ibound];
1029 
1030  // ...is on the bubble
1031  is_on_bubble_bound[bound_id]=true;
1032  }
1033  } // points on bubble boundary located
1034 
1035  // Re-set the boundary conditions for fluid problem: All nodes are
1036  // free by default -- just pin the ones that have Dirichlet conditions
1037  // here.
1038  unsigned nbound=Fluid_mesh_pt->nboundary();
1039  for(unsigned ibound=0;ibound<nbound;ibound++)
1040  {
1041  unsigned num_nod=Fluid_mesh_pt->nboundary_node(ibound);
1042  for (unsigned inod=0;inod<num_nod;inod++)
1043  {
1044  // Get node
1045  Node* nod_pt=Fluid_mesh_pt->boundary_node_pt(ibound,inod);
1046 
1047  //Pin both velocities on inflow (0) and side boundaries (1 and 3)
1048  if((ibound==0) || (ibound==1) || (ibound==3))
1049  {
1050  nod_pt->pin(0);
1051  nod_pt->pin(1);
1052  }
1053 
1054  //If it's the outflow pin only the vertical velocity
1055  if(ibound==2) {nod_pt->pin(1);}
1056 
1057  // Pin pseudo-solid positions apart from bubble boundary which
1058  // we allow to move
1059  SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
1060  if(is_on_bubble_bound[ibound])
1061  {
1062  solid_node_pt->unpin_position(0);
1063  solid_node_pt->unpin_position(1);
1064  }
1065  else
1066  {
1067  solid_node_pt->pin_position(0);
1068  solid_node_pt->pin_position(1);
1069  }
1070  }
1071  } // end loop over boundaries
1072 
1073  // Complete the build of all elements so they are fully functional
1074  // Remember that adaptation for triangle meshes involves a complete
1075  // regneration of the mesh (rather than splitting as in tree-based
1076  // meshes where such parameters can be passed down from the father
1077  // element!)
1078  unsigned n_element = Fluid_mesh_pt->nelement();
1079  for(unsigned e=0;e<n_element;e++)
1080  {
1081  // Upcast from GeneralisedElement to the present element
1082  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
1083 
1084  // Set the Reynolds number
1085  el_pt->re_pt() = &Problem_Parameter::Re;
1086 
1087  // Set the Womersley number (same as Re since St=1)
1088  el_pt->re_st_pt() = &Problem_Parameter::Re;
1089 
1090  // Set the constitutive law for pseudo-elastic mesh deformation
1091  el_pt->constitutive_law_pt()=Problem_Parameter::Constitutive_law_pt;
1092  }
1093 
1094  // Re-apply boundary values on Dirichlet boundary conditions
1095  // (Boundary conditions are ignored when the solution is transferred
1096  // from the old to the new mesh by projection; this leads to a slight
1097  // change in the boundary values (which are, of course, never changed,
1098  // unlike the actual unknowns for which the projected values only
1099  // serve as an initial guess)
1100 
1101  // Set velocity and history values of velocity on walls
1102  nbound=this->Fluid_mesh_pt->nboundary();
1103  for(unsigned ibound=0;ibound<nbound;++ibound)
1104  {
1105  if ((ibound==Upper_wall_boundary_id)||
1106  (ibound==Bottom_wall_boundary_id)||
1107  (ibound==Outflow_boundary_id)||
1108  (ibound==Inflow_boundary_id))
1109  {
1110  // Loop over nodes on this boundary
1111  unsigned num_nod=this->Fluid_mesh_pt->nboundary_node(ibound);
1112  for (unsigned inod=0;inod<num_nod;inod++)
1113  {
1114  // Get node
1115  Node* nod_pt=this->Fluid_mesh_pt->boundary_node_pt(ibound,inod);
1116 
1117  // Get number of previous (history) values
1118  unsigned n_prev=nod_pt->time_stepper_pt()->nprev_values();
1119 
1120  // Velocity is and was zero at all previous times
1121  for (unsigned t=0;t<=n_prev;t++)
1122  {
1123  if (ibound!=Inflow_boundary_id)
1124  {
1125  // Parallel outflow
1126  if (ibound!=Outflow_boundary_id)
1127  {
1128  nod_pt->set_value(t,0,0.0);
1129  }
1130  nod_pt->set_value(t,1,0.0);
1131  }
1132 
1133  // Nodes have always been there...
1134  nod_pt->x(t,0)=nod_pt->x(0,0);
1135  nod_pt->x(t,1)=nod_pt->x(0,1);
1136  }
1137  }
1138  }
1139  }
1140 
1141  // Re-assign prescribed inflow velocity at inlet
1142  unsigned num_nod=this->Fluid_mesh_pt->nboundary_node(Inflow_boundary_id);
1143  for (unsigned inod=0;inod<num_nod;inod++)
1144  {
1145  // Get node
1146  Node* nod_pt=this->Fluid_mesh_pt->boundary_node_pt(Inflow_boundary_id,
1147  inod);
1148  //Now set the boundary velocity
1149  double y = nod_pt->x(1);
1150  nod_pt->set_value(0,Problem_Parameter::Inflow_veloc_magnitude*y*(1-y));
1151  }
1152 
1153 } // end of complete_problem_setup
1154 
1155 
1156 //==start_of_doc_solution=================================================
1157 /// Doc the solution
1158 //========================================================================
1159 template<class ELEMENT>
1160 void BubbleInChannelProblem<ELEMENT>::doc_solution(const std::string& comment)
1161 {
1162  oomph_info << "Docing step: " << Problem_Parameter::Doc_info.number()
1163  << std::endl;
1164 
1165  ofstream some_file;
1166  char filename[100];
1167  sprintf(filename,"%s/soln%i.dat",
1168  Problem_Parameter::Doc_info.directory().c_str(),
1169  Problem_Parameter::Doc_info.number());
1170 
1171  // Number of plot points
1172  unsigned npts;
1173  npts=5;
1174 
1175  // Compute errors and assign to each element for plotting
1176  double max_err;
1177  double min_err;
1178  compute_error_estimate(max_err,min_err);
1179 
1180  // Assemble square of L2 norm
1181  double square_of_l2_norm=0.0;
1182  unsigned nel=Fluid_mesh_pt->nelement();
1183  for (unsigned e=0;e<nel;e++)
1184  {
1185  square_of_l2_norm+=
1186  dynamic_cast<ELEMENT*>(this->Fluid_mesh_pt->element_pt(e))->
1187  square_of_l2_norm();
1188  }
1189  Problem_Parameter::Norm_file << sqrt(square_of_l2_norm) << std::endl;
1190 
1191 
1192  some_file.open(filename);
1193  some_file << dynamic_cast<ELEMENT*>(this->Fluid_mesh_pt->element_pt(0))
1194  ->variable_identifier();
1195  this->Fluid_mesh_pt->output(some_file,npts);
1196  some_file << "TEXT X = 25, Y = 78, CS=FRAME T = \"Global Step "
1197  << Problem_Parameter::Doc_info.number() << " "
1198  << comment << "\"\n";
1199  some_file.close();
1200 
1201 
1202 
1203  // Output boundaries
1204  sprintf(filename,"%s/boundaries%i.dat",
1205  Problem_Parameter::Doc_info.directory().c_str(),
1206  Problem_Parameter::Doc_info.number());
1207  some_file.open(filename);
1208  this->Fluid_mesh_pt->output_boundaries(some_file);
1209  some_file.close();
1210 
1211  // Get max/min area
1212  double max_area;
1213  double min_area;
1214  Fluid_mesh_pt->max_and_min_element_size(max_area, min_area);
1215 
1216  // Get total volume enclosed by face elements (ignore first one)
1217  double vol=0.0;
1219 
1220  // Write trace file
1222  << this->time_pt()->time() << " "
1223  << Fluid_mesh_pt->nelement() << " "
1224  << max_area << " "
1225  << min_area << " "
1226  << max_err << " "
1227  << min_err << " "
1228  << sqrt(square_of_l2_norm) << " "
1229  << vol << " "
1230  << std::endl;
1231 
1232  // Increment the doc_info number
1233  Problem_Parameter::Doc_info.number()++;
1234 
1235 } //end_of_doc_solution
1236 
1237 //========================================================================
1238 /// Compute error estimates and assign to elements for plotting
1239 //========================================================================
1240 template<class ELEMENT>
1242  double& min_err)
1243 {
1244  // Get error estimator
1245  ErrorEstimator* err_est_pt=Fluid_mesh_pt->spatial_error_estimator_pt();
1246 
1247  // Get/output error estimates
1248  unsigned nel=Fluid_mesh_pt->nelement();
1249  Vector<double> elemental_error(nel);
1250 
1251  // We need a dynamic cast, get_element_errors from the Fluid_mesh_pt
1252  // Dynamic cast is used because get_element_errors require a Mesh* ans
1253  // not a SolidMesh*
1254  Mesh* fluid_mesh_pt=dynamic_cast<Mesh*>(Fluid_mesh_pt);
1255  err_est_pt->get_element_errors(fluid_mesh_pt,
1256  elemental_error);
1257 
1258  // Set errors for post-processing and find extrema
1259  max_err=0.0;
1260  min_err=DBL_MAX;
1261  for (unsigned e=0;e<nel;e++)
1262  {
1263  dynamic_cast<MyTaylorHoodElement*>(Fluid_mesh_pt->element_pt(e))->
1264  set_error(elemental_error[e]);
1265 
1266  max_err=std::max(max_err,elemental_error[e]);
1267  min_err=std::min(min_err,elemental_error[e]);
1268  }
1269 
1270 }
1271 
1272 
1273 //==========start_of_main=====================================
1274 ///Driver code for moving bubble problem
1275 //============================================================
1276 int main(int argc, char **argv)
1277 {
1278 
1279  // Store command line arguments
1280  CommandLineArgs::setup(argc,argv);
1281 
1282  // Define possible command line arguments and parse the ones that
1283  // were actually specified
1284 
1285  // Validation?
1286  CommandLineArgs::specify_command_line_flag("--validation");
1287 
1288  // Parse command line
1289  CommandLineArgs::parse_and_assign();
1290 
1291  // Doc what has actually been specified on the command line
1292  CommandLineArgs::doc_specified_flags();
1293 
1294  // Create generalised Hookean constitutive equations
1296  new GeneralisedHookean(&Problem_Parameter::Nu);
1297 
1298  // Open trace file
1299  Problem_Parameter::Trace_file.open("RESLT/trace.dat");
1300 
1301  // Increase precision of output
1302  Problem_Parameter::Trace_file.precision(20);
1303 
1304  // Open norm file
1305  Problem_Parameter::Norm_file.open("RESLT/norm.dat");
1306 
1307 
1308  // Create problem in initial configuration
1310  problem;
1311 
1312 
1313  // Output the problem's state with the bubble in its
1314  // initial polygonal representation
1315  problem.doc_solution();
1316 
1317  // Before starting the time-integration we want to "inflate" it to form
1318  // a proper circular bubble. We do this by setting the inflow to zero
1319  // and doing a steady solve (with one adaptation)
1321 
1322  problem.steady_newton_solve(1);
1323 
1324  // If all went well, this should show us a nice circular bubble
1325  // in a stationary fluid
1326  problem.doc_solution();
1327 
1328  // Initialise timestepper
1329  double dt=0.025;
1330  problem.initialise_dt(dt);
1331 
1332  // Perform impulsive start from current state
1333  problem.assign_initial_values_impulsive();
1334 
1335 
1336  // Now switch on the inflow and re-assign the boundary conditions
1337  // (Call to complete_problem_setup() is a bit expensive given that we
1338  // we only want to set the inflow velocity but who cares -- it's just
1339  // a one off.
1341  problem.complete_problem_setup();
1342 
1343 
1344  // Solve problem on fixed mesh
1345  unsigned nstep=6;
1346  if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
1347  {
1348  nstep=2;
1349  oomph_info << "Remeshing after every second step during validation\n";
1350  }
1351  for (unsigned i=0;i<nstep;i++)
1352  {
1353  // Solve the problem
1354  problem.unsteady_newton_solve(dt);
1355  problem.doc_solution();
1356  } // done solution on fixed mesh
1357 
1358 
1359  // Now do a proper loop, doing nstep timesteps before adapting/remeshing
1360  // and repeating the lot ncycle times
1361  unsigned ncycle=1000;
1362  if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
1363  {
1364  ncycle=1;
1365  oomph_info << "Only doing one cycle during validation\n";
1366  }
1367 
1368  // Do the cycles
1369  for(unsigned j=0;j<ncycle;j++)
1370  {
1371  // Allow up to one level of refinement for next solve
1372  unsigned max_adapt=1;
1373 
1374  //Solve problem a few times
1375  for (unsigned i=0;i<nstep;i++)
1376  {
1377  // Solve the problem
1378  problem.unsteady_newton_solve(dt,max_adapt,false);
1379 
1380 
1381  // Build the label for doc and output solution
1382  std::stringstream label;
1383  label << "Adaptation " <<j << " Step "<< i;
1384  problem.doc_solution(label.str());
1385 
1386  // No more refinement for the next nstep steps
1387  max_adapt=0;
1388  }
1389 
1390  }
1391 
1392 
1393 } //End of main
double Ca
Capillary number.
Problem class to simulate inviscid bubble propagating along 2D channel.
double Error
Storage for elemental error estimate – used for post-processing.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of free surface elements.
void set_error(const double &error)
Set error value for post-processing.
MyTaylorHoodElement()
Constructor initialise error.
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of free surface elements.
double Inflow_veloc_magnitude
Scaling factor for inflow velocity (allows it to be switched off to do hydrostatics) ...
void compute_error_estimate(double &max_err, double &min_err)
Compute the error estimates and assign to elements for plotting.
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
TriangleMeshPolygon * Outer_boundary_polyline_pt
Triangle mesh polygon for outer boundary.
Overload TaylorHood element to modify output.
void doc_solution(const std::string &comment="")
Doc the solution.
Mesh * Free_surface_mesh_pt
Pointers to mesh of free surface elements.
void create_volume_constraint_elements()
Create elements that impose volume constraint on the bubble.
double square_of_l2_norm()
Get square of L2 norm of velocity.
ofstream Trace_file
Trace file.
double Radius
Initial radius of bubble.
void delete_free_surface_elements()
Delete free surface elements.
void create_free_surface_elements()
Create free surface elements.
Vector< TriangleMeshPolygon * > Bubble_polygon_pt
Vector storing pointer to the bubble polygons.
double Volume
Volume of the bubble (negative because it's outside the fluid!)
std::string variable_identifier()
Return variable identifier.
RefineableSolidTriangleMesh< ELEMENT > * Fluid_mesh_pt
Pointer to Fluid_mesh.
void delete_volume_constraint_elements()
Delete volume constraint elements.
double Nu
Pseudo-solid Poisson ratio.
Mesh * Volume_constraint_mesh_pt
Pointer to mesh containing elements that impose volume constraint.
void actions_before_newton_solve()
Update the problem specs before solve.
ofstream Norm_file
File to document the norm of the solution (for validation purposes – triangle doesn't give fully repr...
Data * Bubble_pressure_data_pt
Pointer to a global bubble pressure datum.
void output(std::ostream &outfile, const unsigned &nplot)
Overload output function.
void actions_after_newton_solve()
Update the after solve (empty)
double Re
Reynolds number.
int main(int argc, char **argv)
Driver code for moving bubble problem.
void complete_problem_setup()
Set boundary conditions and complete the build of all elements.
VolumeConstraintElement * Vol_constraint_el_pt
Pointer to element that imposes volume constraint for bubble.
double Length
Length of the channel.
DocInfo Doc_info
Doc info object.