rayleigh_instability_insoluble_surfactant.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-2023 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 // Driver for axisymmetric single-layer fluid problem. Plateau-Rayleigh
27 // instability (unstable if H>2*pi*R -> forming drops)
28 // in the presence of an insoluble surfactant. The problem is described in
29 // A 2-D model of Rayleigh instability in capillary tubes --- surfactant
30 // effects by D. Campana, J. Di Paolo & F. A. Saita, Int. J. Multiphase Flow,
31 // vol 30, pp 431--454, (2004).
32 // The initial version of this code was developed by Aman Rajvardhan.
33 
34 // Generic oomph-lib header
35 #include "generic.h"
36 
37 // Axisymmetric Navier-Stokes headers
38 #include "axisym_navier_stokes.h"
39 
40 // Interface headers
41 #include "fluid_interface.h"
42 
43 // The mesh, including horizontal spines
44 #include "meshes/horizontal_single_layer_spine_mesh.h"
45 
46 //Use the oomph and std namespaces
47 using namespace oomph;
48 using namespace std;
49 
50 //==start_of_namespace===================================================
51 /// Namespace for physical parameters
52 /// The parameter values are chosen to be those used in Figures 8, 9
53 /// in Campana et al.
54 //=======================================================================
56 {
57 
58  //Film thickness parameter
59  double Film_Thickness = 0.2;
60 
61  /// Reynolds number
62  double Re = 40.0;
63 
64  /// Womersley number
65  double ReSt = Re; // (St = 1)
66 
67  /// Product of Reynolds number and inverse of Froude number
68  double ReInvFr = 0.0; // (Fr = 0)
69 
70  /// Capillary number
71  double Ca = pow(Film_Thickness,3.0);
72 
73  /// External pressure
74  double P_ext = 0.0;
75 
76  /// Direction of gravity
77  Vector<double> G(3);
78 
79  /// Wavelength of the domain
80  double Alpha = 1.047;
81 
82  /// Free surface cosine deformation parameter
83  double Epsilon = 1.0e-3;
84 
85  /// Surface Elasticity number (weak case)
86  double Beta = 3.6e-3;
87 
88  /// Surface Peclet number
89  double Peclet_S = 4032.0;
90 
91  /// Sufrace Peclet number multiplied by Strouhal number
92  double Peclet_St_S = 1.0;
93 
94  /// Pvd file -- a wrapper for all the different
95  /// vtu output files plus information about continuous time
96  /// to facilitate animations in paraview
97  ofstream Pvd_file;
98 
99 } // End of namespace
100 
101 
102 
103 namespace oomph
104 {
105 
106 //==================================================================
107 /// Spine-based Marangoni surface tension elements that add
108 /// a linear dependence on concentration
109 /// of a surface chemical to the surface tension,
110 /// which decreases with increasing concentration.
111 /// The non-dimensionalisation is the same as Campana et al (2004)
112 /// but we may wish to revisit this.
113 //=================================================================
114 template<class ELEMENT>
117 {
118 private:
119  /// Pointer to an Elasticity number
120  double *Beta_pt;
121 
122  /// Pointer to Surface Peclet number
123  double *Peclet_S_pt;
124 
125  /// Pointer to the surface Peclect Strouhal number
127 
128  /// Index at which the surfactant concentration is stored at the
129  /// nodes
130  unsigned C_index;
131 
132  /// Default value of the physical constants
134 
135 protected:
136 
137  /// Get the surfactant concentration
138  double interpolated_C(const Vector<double> &s)
139  {
140  //Find number of nodes
141  unsigned n_node = this->nnode();
142 
143  //Get the nodal index at which the unknown is stored
144  const unsigned c_index = C_index;
145 
146  //Local shape function
147  Shape psi(n_node);
148 
149  //Find values of shape function
150  this->shape(s,psi);
151 
152  //Initialise value of C
153  double C = 0.0;
154 
155  //Loop over the local nodes and sum
156  for(unsigned l=0;l<n_node;l++)
157  {
158  C += this->nodal_value(l,c_index)*psi(l);
159  }
160 
161  return(C);
162  }
163 
164  /// The time derivative of the surface concentration
165  double dcdt_surface(const unsigned &l) const
166  {
167  // Get the data's timestepper
168  TimeStepper* time_stepper_pt= this->node_pt(l)->time_stepper_pt();
169 
170  //Initialise dudt
171  double dcdt=0.0;
172  //Loop over the timesteps, if there is a non Steady timestepper
173  if (time_stepper_pt->type()!="Steady")
174  {
175  //Find the index at which the variable is stored
176  const unsigned c_index = C_index;
177 
178  // Number of timsteps (past & present)
179  const unsigned n_time = time_stepper_pt->ntstorage();
180 
181  for(unsigned t=0;t<n_time;t++)
182  {
183  dcdt += time_stepper_pt->weight(1,t)*this->nodal_value(t,l,c_index);
184  }
185  }
186  return dcdt;
187  }
188 
189  /// The surface tension function is linear in the
190  /// concentration with constant of proportionality equal
191  /// to the elasticity number.
192  double sigma(const Vector<double> &s)
193  {
194  //Find the number of shape functions
195  const unsigned n_node = this->nnode();
196  //Now get the shape fuctions at the local coordinate
197  Shape psi(n_node);
198  this->shape(s,psi);
199 
200  //Now interpolate the surfactant concentration
201  double C=0.0;
202  for(unsigned l=0;l<n_node;l++)
203  {
204  C += this->nodal_value(l,C_index)*psi(l);
205  }
206 
207  //Get the Elasticity number
208  double Beta = this->beta();
209  //Return the variable surface tension
210  return (1.0 - Beta*(C-1.0));
211  } // End of sigma
212 
213  /// Fill in the contribution to the residuals
214  /// Calculate the contribution to the jacobian
215  void fill_in_contribution_to_jacobian(Vector<double> &residuals,DenseMatrix<double> &jacobian)
216  {
217  //Call the generic routine with the flag set to 1
218  this->fill_in_generic_residual_contribution_interface(residuals,jacobian,1);
219  {
220  //Use finite differences to handle concentration variations in the
221  //bulk equations
222  const unsigned n_node = this->nnode();
223  //Find the number of dofs in the element
224  const unsigned n_dof = this->ndof();
225  //Create newres vector
226  Vector<double> newres(n_dof);
227 
228  //Integer storage for local unknown
229  int local_unknown=0;
230 
231  //Use the default finite difference step
232  const double fd_step = this->Default_fd_jacobian_step;
233 
234  //Loop over the nodes again
235  for(unsigned n=0;n<n_node;n++)
236  {
237  //Get the number of values stored at the node
238  unsigned c_index = this->C_index;
239 
240  //Get the local equation number
241  local_unknown = this->nodal_local_eqn(n,c_index);
242  //If it's not pinned
243  if(local_unknown >= 0)
244  {
245  //Store a pointer to the nodal data value
246  double *value_pt = this->node_pt(n)->value_pt(c_index);
247 
248  //Save the old value of the Nodal data
249  double old_var = *value_pt;
250 
251  //Increment the value of the Nodal data
252  *value_pt += fd_step;
253 
254  //Calculate the new residuals
255  this->get_residuals(newres);
256 
257  //Do finite differences
258  for(unsigned m=0;m<n_dof;m++)
259  {
260  double sum = (newres[m] - residuals[m])/fd_step;
261  //Stick the entry into the Jacobian matrix
262  jacobian(m,local_unknown) = sum;
263  }
264 
265  //Reset the Nodal data
266  *value_pt = old_var;
267  }
268  }
269  }
270 
271  //Call the generic routine to handle the spine variables
272  SpineElement<FaceGeometry<ELEMENT> >::fill_in_jacobian_from_geometric_data(jacobian);
273  }
274 
275 
276  /// Overload the Helper function to calculate the residuals and
277  /// jacobian entries. This particular function ensures that the
278  /// additional entries are calculated inside the integration loop
279  void add_additional_residual_contributions_interface(Vector<double> &residuals, DenseMatrix<double> &jacobian,
280  const unsigned &flag,const Shape &psif, const DShape &dpsifds,
281  const DShape &dpsifdS, const DShape &dpsifdS_div,
282  const Vector<double> &s,
283  const Vector<double> &interpolated_x, const Vector<double> &interpolated_n,
284  const double &W,const double &J)
285  {
286  //Flag to control whether the Campana formulation (false)
287  // or our own (true) is used
288  bool Integrated_curvature = true;
289 
290  //Find out how many nodes there are
291  unsigned n_node = this->nnode();
292 
293  //Storage for the local equation numbers and unknowns
294  int local_eqn = 0, local_unknown = 0;
295 
296  //Surface advection-diffusion equation
297 
298  //Find the index at which the concentration is stored
299  unsigned c_index = this->C_index;
300  Vector<unsigned> u_index = this->U_index_interface;
301 
302  //Read out the surface peclect number
303  const double Pe_s = this->peclet_s();
304  //const double PeSt_s = this->peclet_strouhal_s();
305 
306  //Read out the radial position
307  const double r = interpolated_x[0];
308 
309  //Rescale the jacobian to be the "raw" version
310  const double J_raw = J/r;
311 
312  //Now calculate the concentration at this point
313  //Assuming the same shape functions are used (which they are)
314  double interpolated_C = 0.0;
315  double interpolated_dCds = 0.0;
316  double dCdt = 0.0;
317  //The tangent vector
318  const unsigned ndim = this->node_pt(0)->ndim();
319  Vector<double> interpolated_tangent(ndim,0.0);
320  Vector<double> interpolated_u(ndim,0.0);
321  Vector<double> mesh_velocity(ndim,0.0);
322  Vector<double> interpolated_duds(ndim,0.0);
323  if(ndim+1 != u_index.size())
324  {
325  throw OomphLibError("Dimension Incompatibility",
326  OOMPH_CURRENT_FUNCTION,
327  OOMPH_EXCEPTION_LOCATION);
328  }
329 
330  for(unsigned l=0;l<n_node;l++)
331  {
332  const double psi = psif(l);
333  const double dpsi = dpsifds(l,0);
334  interpolated_C += this->nodal_value(l,c_index)*psi;
335  interpolated_dCds += this->nodal_value(l,c_index)*dpsi;
336  dCdt += dcdt_surface(l)*psi;
337  for(unsigned i=0;i<ndim;i++)
338  {
339  interpolated_tangent[i] += this->nodal_position(l,i)*dpsi;
340  interpolated_u[i] += this->nodal_value(l,u_index[i])*psi;
341  interpolated_duds[i] += this->nodal_value(l,u_index[i])*dpsi;
342  }
343  //Mesh Velocity
344  for(unsigned j=0;j<ndim;j++)
345  {
346  mesh_velocity[j] += this->dnodal_position_dt(l,j)*psi;
347  }
348  }
349 
350 
351  double u_tangent = 0.0, t_length = 0.0;
352  for(unsigned i=0;i<ndim;i++)
353  {
354  u_tangent += interpolated_u[i]*interpolated_tangent[i];
355  t_length += interpolated_tangent[i]*interpolated_tangent[i];
356  }
357 
358  //Work out the second derivative of position
359  Vector<double> d2xds(2,0.0);
360  for(unsigned i=0;i<2;i++)
361  {
362  d2xds[i] = this->nodal_position(0,i) +
363  this->nodal_position(2,i) - 2.0*this->nodal_position(1,i);
364  }
365 
366  //Let's do the first component of the curvature
367  double k1 = 0.0;
368  for(unsigned i=0;i<2;i++)
369  {
370  k1 += (d2xds[i]/(J_raw*J_raw) - interpolated_tangent[i]*(
371  interpolated_tangent[0]*d2xds[0]
372  + interpolated_tangent[1]*d2xds[1])/(J_raw*J_raw*J_raw*J_raw))*interpolated_n[i];
373  }
374 
375  //Second component of the curvature
376  double k2 = - (interpolated_n[0] / r);
377 
378  //Now we add the flux term to the appropriate residuals
379  for(unsigned l=0;l<n_node;l++)
380  {
381  //Read out the apprporiate local equation
382  local_eqn = this->nodal_local_eqn(l,c_index);
383 
384  //If not a boundary condition
385  if(local_eqn >= 0)
386  {
387  //Time derivative
388  residuals[local_eqn] += dCdt*psif(l)*r*W*J_raw;
389 
390  double tmp = 0.0;
391  //Compute our version in which second derivatives are not required
392  if(Integrated_curvature)
393  {
394  for(unsigned i=0;i<2;i++)
395  {
396  tmp +=
397  (interpolated_u[i] - mesh_velocity[i])*interpolated_tangent[i];
398  }
399  //First Advection term
400  residuals[local_eqn] += tmp*interpolated_dCds*psif(l)*r*W/J_raw;
401  //Additional term from axisymmetric formulation
402  residuals[local_eqn] += interpolated_C*interpolated_u[0]*psif(l)*W*J_raw;
403  //Second Advection term
404  residuals[local_eqn] += interpolated_C*
405  (interpolated_duds[0]*interpolated_tangent[0]
406  + interpolated_duds[1]*interpolated_tangent[1])*r*W*psif(l)/J_raw;
407  }
408  //This is the Campana formulation
409  else
410  {
411  //ALE term
412  for(unsigned i=0;i<2;i++)
413  {
414  tmp += -mesh_velocity[i]*interpolated_tangent[i];
415  }
416  residuals[local_eqn] += tmp*interpolated_dCds*psif(l)*r*W/J_raw;
417  // Curvature (normal velocity) term
418  residuals[local_eqn] -=
419  interpolated_C*(k1 + k2)*psif(l)*r*W*J_raw*(
420  interpolated_u[0]*interpolated_n[0]
421  + interpolated_u[1]*interpolated_n[1]);
422  //Integrated by parts tangential advection term
423  residuals[local_eqn] -= interpolated_C*(
424  interpolated_tangent[0]*interpolated_u[0] +
425  interpolated_tangent[1]*interpolated_u[1])*dpsifds(l,0)*r*W/J_raw;
426  }
427 
428  //Diffusion term
429  residuals[local_eqn] += (1.0/Pe_s)*interpolated_dCds*dpsifds(l,0)*r*W/J_raw;
430 
431  //We also need to worry about the jacobian terms
432  if(flag)
433  {
434  //Loop over the nodes again
435  for(unsigned l2=0;l2<n_node;l2++)
436  {
437  //Get the time stepper
438  TimeStepper* time_stepper_pt=this->node_pt(l2)->time_stepper_pt();
439 
440  //Get the unknown c_index
441  local_unknown =this->nodal_local_eqn(l2,c_index);
442 
443  if(local_unknown >=0)
444  {
445  jacobian(local_eqn,local_unknown) +=
446  time_stepper_pt->weight(1,0)*psif(l2)*psif(l)*r*W*J_raw;
447 
448  if(Integrated_curvature)
449  {
450  jacobian(local_eqn,local_unknown) += ((interpolated_u[0] - mesh_velocity[0])*interpolated_tangent[0]
451  + (interpolated_u[1] - mesh_velocity[1])*interpolated_tangent[1])*dpsifds(l2,0)
452  *psif(l)*r*W/J_raw;
453 
454 
455  jacobian(local_eqn,local_unknown) += psif(l2)*interpolated_u[0]*psif(l)*W*J_raw;
456 
457  jacobian(local_eqn,local_unknown) += psif(l2)*(interpolated_tangent[0]*interpolated_duds[0]
458  + interpolated_tangent[1]*interpolated_duds[1])*psif(l)*r*W/J_raw;
459  }
460  else
461  {
462  jacobian(local_eqn,local_unknown) -=
463  (mesh_velocity[0]*interpolated_tangent[0] + mesh_velocity[1]*interpolated_tangent[1])*dpsifds(l2,0)*psif(l)*r*W/J_raw;
464 
465  jacobian(local_eqn,local_unknown) -= psif(l2)*(k1 + k2)*psif(l)*r*W*J_raw*(interpolated_u[0]*interpolated_n[0]
466  + interpolated_u[1]*interpolated_n[1]);
467 
468  jacobian(local_eqn,local_unknown) -= psif(l2)*(interpolated_tangent[0]*interpolated_u[0] + interpolated_tangent[1]*
469  interpolated_u[1])*dpsifds(l,0)*r*W/J_raw;
470  }
471 
472  jacobian(local_eqn,local_unknown) += (1.0/Pe_s)*dpsifds(l2,0)*dpsifds(l,0)*r*W/J_raw;
473 
474  }
475 
476 
477  //Loop over the velocity components
478  for(unsigned i2=0;i2<ndim;i2++)
479  {
480 
481  //Get the unknown
482  local_unknown = this->nodal_local_eqn(l2,u_index[i2]);
483 
484 
485  //If not a boundary condition
486  if(local_unknown >= 0)
487  {
488 
489  // jacobian(local_eqn,local_unknown) += time_stepper_pt->weight(1,0)*psif(l2)*PeSt_s*psif(l)*r*W*J_raw;
490  if(Integrated_curvature)
491  {
492  jacobian(local_eqn,local_unknown) +=
493  interpolated_dCds*psif(l2)*interpolated_tangent[i2]*psif(l)*r*W/J_raw;
494 
495  if(i2==0)
496  {
497  jacobian(local_eqn,local_unknown) += interpolated_C*psif(l2)*W*J_raw;
498  }
499  }
500  else
501  {
502  jacobian(local_eqn,local_unknown) -= interpolated_C*(k1 + k2)*psif(l)*r*W*J_raw*psif(l2)*interpolated_n[i2];
503 
504  jacobian(local_eqn,local_unknown) -= interpolated_C*interpolated_tangent[i2]*psif(l2)*dpsifds(l,0)*r*W/J_raw;
505  }
506 
507  jacobian(local_eqn,local_unknown) +=
508  interpolated_C*interpolated_tangent[i2]*dpsifds(l2,0)*psif(l)*r*W/J_raw;
509  }
510  }
511  }
512  }
513  }
514  } //End of loop over the nodes
515 
516  }
517 
518 
519  /// Add the element's contribution to its residuals vector,
520  /// jacobian matrix and mass matrix
521  void fill_in_contribution_to_jacobian_and_mass_matrix( Vector<double> &residuals, DenseMatrix<double> &jacobian,
522  DenseMatrix<double> &mass_matrix)
523  {
524  //Add the contribution to the jacobian
525  this->fill_in_contribution_to_jacobian(residuals,jacobian);
526  //No mass matrix terms, but should probably do kinematic bit here
527  }
528 
529 public:
530  /// Constructor that passes the bulk element and face index down
531  /// to the underlying
532  SpineAxisymmetricMarangoniSurfactantFluidInterfaceElement(FiniteElement* const &element_pt, const int &face_index) : SpineAxisymmetricFluidInterfaceElement<ELEMENT>
533  (element_pt,face_index)
534  {
535  //Initialise the values
536  Beta_pt = &Default_Physical_Constant_Value;
537  Peclet_S_pt = &Default_Physical_Constant_Value;
538  Peclet_Strouhal_S_pt = &Default_Physical_Constant_Value;
539 
540  //Add the additional surfactant terms to these surface elements
541 
542  //Read out the number of nodes on the face
543  //For some reason I need to specify the this pointer here(!)
544  unsigned n_node_face = this->nnode();
545  //Set the additional data values in the face
546  //There is one additional values at each node --- the lagrange multiplier
547  Vector<unsigned> additional_data_values(n_node_face);
548  for(unsigned i=0;i<n_node_face;i++) additional_data_values[i] = 1;
549  //Resize the data arrays accordingly
550  this->resize_nodes(additional_data_values);
551 
552  //The C_index is the new final value
553  //Minor HACK HERE
554  C_index = this->node_pt(0)->nvalue()-1;
555 }
556 
557  /// Return the Elasticity number
558  double beta() {return *Beta_pt;}
559 
560  /// Return the surface peclect number
561  double peclet_s() {return *Peclet_S_pt;}
562 
563  /// Return the surface peclect strouhal number
564  double peclet_strouhal_s() {return *Peclet_Strouhal_S_pt;}
565 
566  /// Access function for pointer to the Elasticity number
567  double* &beta_pt() {return Beta_pt;}
568 
569  /// Access function for pointer to the surface Peclet number
570  double* &peclet_s_pt() {return Peclet_S_pt;}
571 
572  /// Access function for pointer to the surface Peclet x Strouhal number
573  double* &peclet_strouhal_s_pt() {return Peclet_Strouhal_S_pt;}
574 
575  void output(std::ostream &outfile, const unsigned &n_plot)
576  {
577  outfile.precision(16);
578 
579  //Set output Vector
580  Vector<double> s(1);
581 
582  //Tecplot header info
583  outfile << "ZONE I=" << n_plot << std::endl;
584 
585  const unsigned n_node = this->nnode();
586  const unsigned dim = this->dim();
587 
588  Shape psi(n_node);
589  DShape dpsi(n_node,dim);
590 
591  //Loop over plot points
592  for(unsigned l=0;l<n_plot;l++)
593  {
594  s[0] = -1.0 + l*2.0/(n_plot-1);
595 
596  this->dshape_local(s,psi,dpsi);
597  Vector<double> interpolated_tangent(2,0.0);
598  for(unsigned l=0;l<n_node;l++)
599  {
600  const double dpsi_ = dpsi(l,0);
601  for(unsigned i=0;i<2;i++)
602  {
603  interpolated_tangent[i] += this->nodal_position(l,i)*dpsi_;
604  }
605  }
606 
607  //Tangent
608  double t_vel = (this->interpolated_u(s,0)*interpolated_tangent[0] + this->interpolated_u(s,1)*interpolated_tangent[1])/
609  sqrt(interpolated_tangent[0]*interpolated_tangent[0] + interpolated_tangent[1]*interpolated_tangent[1]);
610 
611 
612  //Output the x,y,u,v
613  for(unsigned i=0;i<2;i++) outfile << this->interpolated_x(s,i) << " ";
614  for(unsigned i=0;i<2;i++) outfile << this->interpolated_u(s,i) << " ";
615  //Output a dummy pressure
616  outfile << 0.0 << " ";
617  //Output the concentration
618  outfile << interpolated_C(s) << " ";
619  //Output the interfacial tension
620  outfile << sigma(s) << " " << t_vel << std::endl;
621  }
622  outfile << std::endl;
623  }
624 
625 
626  /// Compute the concentration intergated over the area
627  double integrate_c() const
628  {
629  //Find out how many nodes there are
630  unsigned n_node = this->nnode();
631 
632  //Set up memeory for the shape functions
633  Shape psif(n_node);
634  DShape dpsifds(n_node,1);
635 
636  //Set the value of n_intpt
637  unsigned n_intpt = this->integral_pt()->nweight();
638 
639  //Storage for the local coordinate
640  Vector<double> s(1);
641 
642  //Storage for the total mass
643  double mass = 0.0;
644 
645  //Loop over the integration points
646  for(unsigned ipt=0;ipt<n_intpt;ipt++)
647  {
648  //Get the local coordinate at the integration point
649  s[0] = this->integral_pt()->knot(ipt,0);
650 
651  //Get the integral weight
652  double W = this->integral_pt()->weight(ipt);
653 
654  //Call the derivatives of the shape function
655  this->dshape_local_at_knot(ipt,psif,dpsifds);
656 
657  //Define and zero the tangent Vectors and local velocities
658  Vector<double> interpolated_x(2,0.0);
659  Vector<double> interpolated_t(2,0.0);
660  double interpolated_c = 0.0;
661 
662 
663  //Loop over the shape functions to compute concentration and tangent
664  for(unsigned l=0;l<n_node;l++)
665  {
666  interpolated_c += this->nodal_value(l,C_index)*psif(l);
667  //Calculate the tangent vector
668  for(unsigned i=0;i<2;i++)
669  {
670  interpolated_x[i] += this->nodal_position(l,i)*psif(l);
671  interpolated_t[i] += this->nodal_position(l,i)*dpsifds(l,0);
672  }
673  }
674 
675  //The first positional coordinate is the radial coordinate
676  double r = interpolated_x[0];
677 
678  //Calculate the length of the tangent Vector
679  double tlength = interpolated_t[0]*interpolated_t[0] +
680  interpolated_t[1]*interpolated_t[1];
681 
682  //Set the Jacobian of the line element
683  double J = sqrt(tlength);
684 
685  mass += interpolated_c*r*W*J;
686  }
687  return mass;
688  }
689 
690 };
691 
692 
693 //Define the default physical value to be one
694 template<class ELEMENT>
696 
697 }
698 
699 
700 namespace oomph
701 {
702 
703 //======================================================================
704 /// Inherit from the standard Horizontal single-layer SpineMesh
705 /// and modify the spine_node_update() function so that it is appropriate
706 /// for an annular film, rather than a fluid fibre.
707 //======================================================================
708 template <class ELEMENT>
710  public HorizontalSingleLayerSpineMesh<ELEMENT>
711 {
712 
713 public:
714 
715  /// Constructor: Pass number of elements in x-direction, number of
716  /// elements in y-direction, radial extent, axial length , and pointer
717  /// to timestepper (defaults to Steady timestepper)
718  MyHorizontalSingleLayerSpineMesh(const unsigned &nx,
719  const unsigned &ny,
720  const double &lx,
721  const double &ly,
722  TimeStepper* time_stepper_pt=
723  &Mesh::Default_TimeStepper) :
724  HorizontalSingleLayerSpineMesh<ELEMENT>(nx,ny,lx,ly,time_stepper_pt) {}
725 
726  /// Node update function assumed spines rooted at the wall
727  /// fixed to be at r=1 and directed inwards to r=0.
728  virtual void spine_node_update(SpineNode* spine_node_pt)
729  {
730  //Get fraction along the spine
731  double W = spine_node_pt->fraction();
732 
733  //Get spine height
734  double H = spine_node_pt->h();
735 
736  //Set the value of y
737  spine_node_pt->x(0) = 1.0-(1.0-W)*H;
738  }
739 
740 };
741 
742 }
743 
744 
745 //==start_of_problem_class===============================================
746 /// Single axisymmetric fluid interface problem including the
747 /// transport of an insoluble surfactant.
748 //=======================================================================
749 template<class ELEMENT, class TIMESTEPPER>
750 class InterfaceProblem : public Problem
751 {
752 
753 public:
754 
755  /// Constructor: Pass the number of elements in radial and axial directions
756  /// and the length of the domain in the z direction)
757  InterfaceProblem(const unsigned &n_r, const unsigned &n_z,
758  const double &l_z);
759 
760  /// Destructor (empty)
762 
763  /// Spine heights/lengths are unknowns in the problem so their values get
764  /// corrected during each Newton step. However, changing their value does
765  /// not automatically change the nodal positions, so we need to update all
766  /// of them here.
768  {
769  Bulk_mesh_pt->node_update();
770  }
771 
772  /// Set initial conditions: Set all nodal velocities to zero and
773  /// initialise the previous velocities to correspond to an impulsive
774  /// start
776  {
777  // Determine number of nodes in mesh
778  const unsigned n_node = Bulk_mesh_pt->nnode();
779 
780  // Loop over all nodes in mesh
781  for(unsigned n=0;n<n_node;n++)
782  {
783  // Loop over the three velocity components
784  for(unsigned i=0;i<3;i++)
785  {
786  // Set velocity component i of node n to zero
787  Bulk_mesh_pt->node_pt(n)->set_value(i,0.0);
788  }
789  }
790 
791  // Initialise the previous velocity values for timestepping
792  // corresponding to an impulsive start
793  assign_initial_values_impulsive();
794 
795  } // End of set_initial_condition
796 
797 
798  /// The global temporal error norm, based on the movement of the nodes
799  /// in the radial direction only (because that's the only direction
800  /// in which they move!)
802  {
803  //Temp
804  double global_error = 0.0;
805 
806  //Find out how many nodes there are in the problem
807  const unsigned n_node = Bulk_mesh_pt->nnode();
808 
809  //Loop over the nodes and calculate the errors in the positions
810  for(unsigned n=0;n<n_node;n++)
811  {
812  //Set the dimensions to be restricted to the radial direction only
813  const unsigned n_dim = 1;
814  //Set the position error to zero
815  double node_position_error = 0.0;
816  //Loop over the dimensions
817  for(unsigned i=0;i<n_dim;i++)
818  {
819  //Get position error
820  double error =
821  Bulk_mesh_pt->node_pt(n)->position_time_stepper_pt()->
822  temporal_error_in_position(Bulk_mesh_pt->node_pt(n),i);
823 
824  //Add the square of the individual error to the position error
825  node_position_error += error*error;
826  }
827 
828  //Divide the position error by the number of dimensions
829  node_position_error /= n_dim;
830  //Now add to the global error
831  global_error += node_position_error;
832  }
833 
834  //Now the global error must be divided by the number of nodes
835  global_error /= n_node;
836 
837  //Return the square root of the errr
838  return sqrt(global_error);
839  }
840 
841 
842  /// Access function for the specific mesh
844 
845  /// Mesh for the free surface (interface) elements
847 
848  /// Doc the solution
849  void doc_solution(DocInfo &doc_info);
850 
851  /// Do unsteady run up to maximum time t_max with given timestep dt
852  void unsteady_run(const double &t_max, const double &dt);
853 
854  /// Compute the total mass of the insoluble surfactant
856  {
857  //Initialise to zero
858  double mass = 0.0;
859 
860  // Determine number of 1D interface elements in mesh
861  const unsigned n_interface_element = Interface_mesh_pt->nelement();
862 
863  // Loop over the interface elements
864  for(unsigned e=0;e<n_interface_element;e++)
865  {
866  // Upcast from GeneralisedElement to the present element
869  (Interface_mesh_pt->element_pt(e));
870  //Add contribution from each element
871  mass += el_pt->integrate_c();
872  }
873  return mass;
874  } // End of compute_total_mass
875 
876 
877 private:
878 
879  /// Deform the mesh/free surface to a prescribed function
880  void deform_free_surface(const double &epsilon)
881  {
882  // Determine number of spines in mesh
883  const unsigned n_spine = Bulk_mesh_pt->nspine();
884 
885  // Loop over spines in mesh
886  for(unsigned i=0;i<n_spine;i++)
887  {
888 
889  // Determine z coordinate of spine
890  double z_value = Bulk_mesh_pt->boundary_node_pt(1,i)->x(1);
891 
892  // Set spine height
893  Bulk_mesh_pt->spine_pt(i)->height() =
895  + epsilon*cos(Global_Physical_Variables::Alpha*z_value));
896 
897  } // End of loop over spines
898 
899  // Update nodes in bulk mesh
900  Bulk_mesh_pt->node_update();
901 
902  } // End of deform_free_surface
903 
904  /// Trace file
905  ofstream Trace_file;
906 
907 }; // End of problem class
908 
909 
910 
911 //==start_of_constructor==================================================
912 /// Constructor for single fluid interface problem
913 //========================================================================
914 template<class ELEMENT, class TIMESTEPPER>
916 InterfaceProblem(const unsigned &n_r,
917  const unsigned &n_z,
918  const double &l_z)
919 
920 {
921  // Allocate the timestepper (this constructs the time object as well)
922  add_time_stepper_pt(new TIMESTEPPER(true));
923 
924  // Build and assign mesh (the "false" boolean flag tells the mesh
925  // constructor that the domain is not periodic in r)
926  Bulk_mesh_pt = new MyHorizontalSingleLayerSpineMesh<ELEMENT>(n_r,n_z,1.0,l_z,time_stepper_pt());
927 
928  //Create "surface mesh" that will only contain the interface elements
929  Interface_mesh_pt = new Mesh;
930  {
931  // How many bulk elements are adjacent to boundary b?
932  // Boundary 1 is the outer boundary
933  // The boundaries are labelled
934  // 2
935  // 3 1
936  // 0
937 
938  unsigned n_element = Bulk_mesh_pt->nboundary_element(3);
939 
940  // Loop over the bulk elements adjacent to boundary b?
941  for(unsigned e=0;e<n_element;e++)
942  {
943  // Get pointer to the bulk element that is adjacent to boundary b
944  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
945  Bulk_mesh_pt->boundary_element_pt(3,e));
946 
947  // Find the index of the face of element e along boundary b
948  int face_index = Bulk_mesh_pt->face_index_at_boundary(3,e);
949 
950  // Build the corresponding free surface element
953 
954  //Add the prescribed-flux element to the surface mesh
955  Interface_mesh_pt->add_element_pt(interface_element_pt);
956 
957  } //end of loop over bulk elements adjacent to boundary b
958  }
959 
960 
961  // Add the two sub meshes to the problem
962  add_sub_mesh(Bulk_mesh_pt);
963  add_sub_mesh(Interface_mesh_pt);
964 
965  // Combine all submeshes into a single Mesh
966  build_global_mesh();
967 
968 
969  // --------------------------------------------
970  // Set the boundary conditions for this problem
971  // --------------------------------------------
972 
973  //Pin all azimuthal velocities
974  {
975  unsigned n_node = this->Bulk_mesh_pt->nnode();
976  for(unsigned n=0;n<n_node;++n)
977  {
978  this->Bulk_mesh_pt->node_pt(n)->pin(2);
979  }
980  }
981 
982  // All nodes are free by default -- just pin the ones that have
983  // Dirichlet conditions here
984  unsigned ibound = 3;
985  unsigned n_node = Bulk_mesh_pt->nboundary_node(ibound);
986  for(unsigned n=0;n<n_node;n++)
987  {
988  Bulk_mesh_pt->boundary_node_pt(ibound,n)->set_value(3,1.0);
989  }
990 
991  // Determine number of mesh boundaries
992  const unsigned n_boundary = Bulk_mesh_pt->nboundary();
993 
994  // Loop over mesh boundaries
995  for(unsigned b=0;b<n_boundary;b++)
996  {
997  // Determine number of nodes on boundary b
998  const unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
999 
1000  // Loop over nodes on boundary b
1001  for(unsigned n=0;n<n_node;n++)
1002  {
1003  // Pin azimuthal velocity on bounds
1004  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(2);
1005 
1006  // Pin velocity on the outer wall
1007  if(b==1)
1008  {
1009  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
1010  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1);
1011  }
1012  // Pin axial velocity on bottom and top walls (no penetration)
1013  if(b==0 || b==2)
1014  {
1015  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1);
1016  }
1017  } // End of loop over nodes on boundary b
1018  } // End of loop over mesh boundaries
1019 
1020  // ----------------------------------------------------------------
1021  // Complete the problem setup to make the elements fully functional
1022  // ----------------------------------------------------------------
1023 
1024  // Determine number of bulk elements in mesh
1025  const unsigned n_bulk = Bulk_mesh_pt->nelement();
1026 
1027  // Loop over the bulk elements
1028  for(unsigned e=0;e<n_bulk;e++)
1029  {
1030  // Upcast from GeneralisedElement to the present element
1031  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
1032 
1033  // Set the Reynolds number
1034  el_pt->re_pt() = &Global_Physical_Variables::Re;
1035 
1036  // Set the Womersley number
1037  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
1038 
1039  // Set the product of the Reynolds number and the inverse of the
1040  // Froude number
1041  el_pt->re_invfr_pt() = &Global_Physical_Variables::ReInvFr;
1042 
1043  // Set the direction of gravity
1044  el_pt->g_pt() = &Global_Physical_Variables::G;
1045 
1046  } // End of loop over bulk elements
1047 
1048  // Create a Data object whose single value stores the external pressure
1049  Data* external_pressure_data_pt = new Data(1);
1050 
1051  // Pin and set the external pressure to some arbitrary value
1052  double p_ext = Global_Physical_Variables::P_ext;
1053 
1054  external_pressure_data_pt->pin(0);
1055  external_pressure_data_pt->set_value(0,p_ext);
1056 
1057  // Determine number of 1D interface elements in mesh
1058  const unsigned n_interface_element = Interface_mesh_pt->nelement();
1059 
1060  // Loop over the interface elements
1061  for(unsigned e=0;e<n_interface_element;e++)
1062  {
1063  // Upcast from GeneralisedElement to the present element
1066  (Interface_mesh_pt->element_pt(e));
1067 
1068  // Set the Capillary number
1070 
1071  // Pass the Data item that contains the single external pressure value
1072  el_pt->set_external_pressure_data(external_pressure_data_pt);
1073 
1074  // Set the surface elasticity number
1076 
1077  // Set the surface peclect number
1079 
1080  // Set the surface peclect number multiplied by strouhal number
1082 
1083  } // End of loop over interface elements
1084 
1085  // Setup equation numbering scheme
1086  cout << "Number of equations: " << assign_eqn_numbers() << std::endl;
1087 
1088 } // End of constructor
1089 
1090 
1091 
1092 //==start_of_doc_solution=================================================
1093 /// Doc the solution
1094 //========================================================================
1095 template<class ELEMENT, class TIMESTEPPER>
1097 doc_solution(DocInfo &doc_info)
1098 {
1099 
1100  // Output the time
1101  double t= time_pt()->time();
1102  cout << "Time is now " << t << std::endl;
1103 
1104  // Document in trace file
1105  Trace_file << time_pt()->time() << " "
1106  << Bulk_mesh_pt->spine_pt(0)->height()
1107  << " " << this->compute_total_mass() << std::endl;
1108 
1109  ofstream some_file;
1110  char filename[100];
1111 
1112  // Set number of plot points (in each coordinate direction)
1113  const unsigned npts = 5;
1114 
1115  // Open solution output file
1116  //sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
1117  // doc_info.number());
1118  //some_file.open(filename);
1119 
1120  // Output solution to file
1121  //Bulk_mesh_pt->output(some_file,npts);
1122  //some_file.close();
1123  //Put interface in separate file
1124  sprintf(filename,"%s/int%i.dat",doc_info.directory().c_str(),
1125  doc_info.number());
1126  some_file.open(filename);
1127  Interface_mesh_pt->output(some_file,npts);
1128  some_file.close();
1129 
1130  // Write file as a tecplot text object...
1131 // some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
1132  // << time_pt()->time() << "\"";
1133  // ...and draw a horizontal line whose length is proportional
1134  // to the elapsed time
1135  //some_file << "GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
1136  //some_file << "1" << std::endl;
1137  //some_file << "2" << std::endl;
1138  //some_file << " 0 0" << std::endl;
1139  //some_file << time_pt()->time()*20.0 << " 0" << std::endl;
1140 
1141  // Close solution output file
1142 // some_file.close();
1143 
1144  // Output solution to file in paraview format
1145  //sprintf(filename,"%s/soln%i.vtu",doc_info.directory().c_str(),
1146  // doc_info.number());
1147  //some_file.open(filename);
1148  //Bulk_mesh_pt->output_paraview(some_file,npts);
1149  //some_file.close();
1150 
1151  // Write pvd information
1152  string file_name="soln"+StringConversion::to_string(doc_info.number())
1153  +".vtu";
1154  ParaviewHelper::write_pvd_information(Global_Physical_Variables::Pvd_file,
1155  file_name,t);
1156 
1157 } // End of doc_solution
1158 
1159 
1160 
1161 //==start_of_unsteady_run=================================================
1162 /// Perform run up to specified time t_max with given timestep dt
1163 //========================================================================
1164 template<class ELEMENT, class TIMESTEPPER>
1166 unsteady_run(const double &t_max, const double &dt)
1167 {
1168 
1169  // Set value of epsilon
1170  double epsilon = Global_Physical_Variables::Epsilon;
1171 
1172  // Deform the mesh/free surface
1173  deform_free_surface(epsilon);
1174 
1175  // Initialise DocInfo object
1176  DocInfo doc_info;
1177 
1178  // Set output directory
1179  doc_info.set_directory("RESLT");
1180 
1181  // Initialise counter for solutions
1182  doc_info.number()=0;
1183 
1184  // Open trace file
1185  char filename[100];
1186  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
1187  Trace_file.open(filename);
1188 
1189  // Initialise trace file
1190  Trace_file << "time" << ", "
1191  << "edge spine height" << ", "
1192  << "contact angle left" << ", "
1193  << "contact angle right" << ", " << std::endl;
1194 
1195  // Initialise timestep
1196  initialise_dt(dt);
1197 
1198  // Set initial conditions
1199  set_initial_condition();
1200 
1201  // Determine number of timesteps
1202  const unsigned n_timestep = unsigned(t_max/dt);
1203 
1204  // Open pvd file -- a wrapper for all the different
1205  // vtu output files plus information about continuous time
1206  // to facilitate animations in paraview
1207  sprintf(filename,"%s/soln.pvd",doc_info.directory().c_str());
1208  Global_Physical_Variables::Pvd_file.open(filename);
1209  ParaviewHelper::write_pvd_header(Global_Physical_Variables::Pvd_file);
1210 
1211  // Doc initial solution
1212  doc_solution(doc_info);
1213 
1214  // Increment counter for solutions
1215  doc_info.number()++;
1216 
1217  //double dt_desired = dt;
1218 
1219  // Timestepping loop
1220  for(unsigned t=1;t<=n_timestep;t++)
1221  {
1222  // Output current timestep to screen
1223  cout << "\nTimestep " << t << " of " << n_timestep << std::endl;
1224 
1225  // Take one fixed timestep
1226  unsteady_newton_solve(dt);
1227 
1228  //double dt_actual =
1229  // adaptive_unsteady_newton_solve(dt_desired,1.0e-6);
1230  //dt_desired = dt_actual;
1231 
1232 
1233  // Doc solution
1234  doc_solution(doc_info);
1235 
1236  // Increment counter for solutions
1237  doc_info.number()++;
1238  } // End of timestepping loop
1239 
1240  // write footer and close pvd file
1241  ParaviewHelper::write_pvd_footer(Global_Physical_Variables::Pvd_file);
1243 
1244 } // End of unsteady_run
1245 
1246 
1247 /// ////////////////////////////////////////////////////////////////////
1248 /// ////////////////////////////////////////////////////////////////////
1249 /// ////////////////////////////////////////////////////////////////////
1250 
1251 
1252 //==start_of_main=========================================================
1253 /// Driver code for single fluid axisymmetric horizontal interface problem
1254 //========================================================================
1255 int main(int argc, char* argv[])
1256 {
1257 
1258  // Store command line arguments
1259  CommandLineArgs::setup(argc,argv);
1260 
1261  /// Maximum time
1262  double t_max = 1000.0;
1263 
1264  /// Duration of timestep
1265  const double dt = 0.1;
1266 
1267  // If we are doing validation run, use smaller number of timesteps
1268  if(CommandLineArgs::Argc>1)
1269  {
1270  t_max = 0.5;
1271  }
1272 
1273  // Number of elements in radial (r) direction
1274  const unsigned n_r = 10;
1275 
1276  // Number of elements in axial (z) direction
1277  const unsigned n_z = 80;
1278 
1279  // Height of domain
1280  const double l_z = MathematicalConstants::Pi/Global_Physical_Variables::Alpha;
1281 
1282  // Set direction of gravity (vertically downwards)
1286 
1287  // Set up the spine test problem with AxisymmetricQCrouzeixRaviartElements,
1288  // using the BDF<2> timestepper
1290  problem(n_r,n_z,l_z);
1291 
1292  // Run the unsteady simulation
1293  problem.unsteady_run(t_max,dt);
1294 } // End of main
1295 
Single fluid interface problem including transport of an insoluble surfactant.
void set_initial_condition()
Set initial conditions: Set all nodal velocities to zero and initialise the previous velocities to co...
void unsteady_run(const unsigned &nstep)
Run an unsteady simulation with specified number of steps.
InterfaceProblem(const unsigned &n_r, const unsigned &n_y, const unsigned &n_theta, const double &r_min, const double &r_max, const double &l_y, const double &theta_max)
Constructor: Pass number of elements in x and y directions. Also lengths of the domain in x- and y-di...
void deform_free_surface(const double &epsilon)
Deform the mesh/free surface to a prescribed function.
void doc_solution(DocInfo &doc_info)
Doc the solution.
double global_temporal_error_norm()
The global temporal error norm, based on the movement of the nodes in the radial direction only (beca...
MyHorizontalSingleLayerSpineMesh< ELEMENT > * Bulk_mesh_pt
Access function for the specific mesh.
double compute_total_mass()
Compute the total mass of the insoluble surfactant.
void actions_before_newton_convergence_check()
Spine heights/lengths are unknowns in the problem so their values get corrected during each Newton st...
Mesh * Interface_mesh_pt
Mesh for the free surface (interface) elements.
double *& ca_pt()
Pointer to the Capillary number.
void set_external_pressure_data(Data *external_pressure_data_pt)
Set the Data that contains the single pressure value that specifies the "external pressure" for the i...
Inherit from the standard Horizontal single-layer SpineMesh and modify the spine_node_update() functi...
MyHorizontalSingleLayerSpineMesh(const unsigned &nx, const unsigned &ny, const double &lx, const double &ly, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in x-direction, number of elements in y-direction,...
virtual void spine_node_update(SpineNode *spine_node_pt)
Node update function assumed spines rooted at the wall fixed to be at r=1 and directed inwards to r=0...
Spine-based Marangoni surface tension elements that add a linear dependence on concentration of a sur...
void add_additional_residual_contributions_interface(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag, const Shape &psif, const DShape &dpsifds, const DShape &dpsifdS, const DShape &dpsifdS_div, const Vector< double > &s, const Vector< double > &interpolated_x, const Vector< double > &interpolated_n, const double &W, const double &J)
Overload the Helper function to calculate the residuals and jacobian entries. This particular functio...
static double Default_Physical_Constant_Value
Default value of the physical constants.
double sigma(const Vector< double > &s)
The surface tension function is linear in the concentration with constant of proportionality equal to...
double dcdt_surface(const unsigned &l) const
The time derivative of the surface concentration.
unsigned C_index
Index at which the surfactant concentration is stored at the nodes.
double interpolated_C(const Vector< double > &s)
Get the surfactant concentration.
double *& peclet_strouhal_s_pt()
Access function for pointer to the surface Peclet x Strouhal number.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in the contribution to the residuals Calculate the contribution to the jacobian.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Add the element's contribution to its residuals vector, jacobian matrix and mass matrix.
SpineAxisymmetricMarangoniSurfactantFluidInterfaceElement(FiniteElement *const &element_pt, const int &face_index)
Constructor that passes the bulk element and face index down to the underlying.
double *& beta_pt()
Access function for pointer to the Elasticity number.
double *& peclet_s_pt()
Access function for pointer to the surface Peclet number.
double integrate_c() const
Compute the concentration intergated over the area.
Namepspace for global parameters, chosen from Campana et al. as in the axisymmetric problem.
double ReSt
Womersley = Reynolds times Strouhal.
double Epsilon
Free surface cosine deformation parameter.
double Beta
Surface Elasticity number (weak case)
ofstream Pvd_file
Pvd file – a wrapper for all the different vtu output files plus information about continuous time to...
double Alpha
Wavelength of the domain.
double ReInvFr
Product of Reynolds and Froude number.
double Peclet_St_S
Sufrace Peclet number multiplied by Strouhal number.
Vector< double > G(3)
Direction of gravity.
int main(int argc, char *argv[])
//////////////////////////////////////////////////////////////////// ////////////////////////////////...