interface_elements.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 // Non-inline functions for fluid free surface elements
27 
28 // OOMPH-LIB headers
29 #include "interface_elements.h"
30 #include "../generic/integral.h"
31 
32 
33 namespace oomph
34 {
35  /// ///////////////////////////////////////////////////////////////
36  /// ///////////////////////////////////////////////////////////////
37  /// ///////////////////////////////////////////////////////////////
38 
39 
40  //=========================================================================
41  /// Set a pointer to the desired contact angle. Optional boolean
42  /// (defaults to true)
43  /// chooses strong imposition via hijacking (true) or weak imposition
44  /// via addition to momentum equation (false). The default strong imposition
45  /// is appropriate for static contact problems.
46  //=========================================================================
48  const bool& strong)
49  {
50  // Set the pointer to the contact angle
51  Contact_angle_pt = angle_pt;
52 
53  // If we are hijacking the kinematic condition (the default)
54  // to do the strong (pointwise form of the contact-angle condition)
55  if (strong)
56  {
57  // Remember what we're doing
59 
60  // Hijack the bulk element residuals
61  dynamic_cast<FluidInterfaceElement*>(bulk_element_pt())
62  ->hijack_kinematic_conditions(Bulk_node_number);
63  }
64  // Otherwise, we'll impose it weakly via the momentum equations.
65  // This will require that the appropriate velocity node is unpinned,
66  // which is why this is a bad choice for static contact problems in which
67  // there is a no-slip condition on the wall. In that case, the momentum
68  // equation is never assembled and so the contact angle condition is not
69  // applied unless we use the strong version above.
70  else
71  {
73  }
74  }
75 
76 
77  /// ///////////////////////////////////////////////////////////////
78  /// ///////////////////////////////////////////////////////////////
79  /// ///////////////////////////////////////////////////////////////
80 
81 
82  //=========================================================================
83  /// Add contribution to element's residual vector and Jacobian
84  //=========================================================================
87  Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
88  {
89  // Let's get the info from the parent
90  FiniteElement* parent_pt = bulk_element_pt();
91 
92  // Find the dimension of the problem
93  unsigned spatial_dim = this->nodal_dimension();
94 
95  // Outer unit normal to the wall
96  Vector<double> wall_normal(spatial_dim);
97 
98  // Outer unit normal to the free surface
99  Vector<double> unit_normal(spatial_dim);
100 
101  // Storage for the coordinate
102  Vector<double> x(spatial_dim);
103 
104  // Find the dimension of the parent
105  unsigned n_dim = parent_pt->dim();
106 
107  // Dummy local coordinate, of size zero
108  Vector<double> s_local(0);
109 
110  // Get the x coordinate
111  this->interpolated_x(s_local, x);
112 
113  // Get the unit normal to the wall
114  wall_unit_normal(x, wall_normal);
115 
116  // Find the local coordinates in the parent
117  Vector<double> s_parent(n_dim);
118  this->get_local_coordinate_in_bulk(s_local, s_parent);
119 
120  // Just get the outer unit normal
121  dynamic_cast<FaceElement*>(parent_pt)->outer_unit_normal(s_parent,
122  unit_normal);
123 
124  // Find the dot product of the two vectors
125  double dot = 0.0;
126  for (unsigned i = 0; i < spatial_dim; i++)
127  {
128  dot += unit_normal[i] * wall_normal[i];
129  }
130 
131  // Get the value of sigma from the parent
132  double sigma_local =
133  dynamic_cast<FluidInterfaceElement*>(parent_pt)->sigma(s_parent);
134 
135  // Are we doing the weak form replacement
136  if (Contact_angle_flag == 2)
137  {
138  // Get the wall tangent vector
139  Vector<double> wall_tangent(spatial_dim);
140  wall_tangent[0] = -wall_normal[1];
141  wall_tangent[1] = wall_normal[0];
142 
143  // Get the capillary number
144  double ca_local = ca();
145 
146  // Just add the appropriate contribution to the momentum equations
147  for (unsigned i = 0; i < 2; i++)
148  {
149  int local_eqn = nodal_local_eqn(0, this->U_index_interface_boundary[i]);
150  if (local_eqn >= 0)
151  {
152  residuals[local_eqn] +=
153  (sigma_local / ca_local) * (sin(contact_angle()) * wall_normal[i] +
154  cos(contact_angle()) * wall_tangent[i]);
155  }
156  }
157  }
158  // Otherwise [strong imposition (by hijacking) of contact angle or
159  // "no constraint at all"], add the appropriate contribution to
160  // the momentum equation
161  else
162  {
163  // Need to find the current outer normal from the surface
164  // which does not necessarily correspond to an imposed angle.
165  // It is whatever it is...
166  Vector<double> m(spatial_dim);
167  this->outer_unit_normal(s_local, m);
168 
169  // Get the capillary number
170  double ca_local = ca();
171 
172  // Just add the appropriate contribution to the momentum equations
173  // This will, of course, not be added if the equation is pinned
174  //(no slip)
175  for (unsigned i = 0; i < 2; i++)
176  {
177  int local_eqn = nodal_local_eqn(0, this->U_index_interface_boundary[i]);
178  if (local_eqn >= 0)
179  {
180  residuals[local_eqn] += (sigma_local / ca_local) * m[i];
181  }
182  }
183  }
184 
185  // If we are imposing the contact angle strongly (by hijacking)
186  // overwrite the kinematic equation
187  if (Contact_angle_flag == 1)
188  {
189  // Read out the kinematic equation number
190  int local_eqn = kinematic_local_eqn(0);
191 
192  // Note that because we have outer unit normals for the free surface
193  // and the wall, the cosine of the contact angle is equal to
194  // MINUS the dot product computed above
195  if (local_eqn >= 0)
196  {
197  residuals[local_eqn] = cos(contact_angle()) + dot;
198  }
199  // NOTE: The jacobian entries will be computed automatically
200  // by finite differences.
201  }
202 
203  // Dummy arguments
204  Shape psif(1);
205  DShape dpsifds(1, 1);
206  Vector<double> interpolated_n(1);
207  double W = 0.0;
208 
209  // Now add the additional contributions
211  residuals, jacobian, flag, psif, dpsifds, interpolated_n, W);
212  }
213 
214 
215  /// ///////////////////////////////////////////////////////////////
216  /// ///////////////////////////////////////////////////////////////
217  /// ///////////////////////////////////////////////////////////////
218 
219 
220  //=========================================================================
221  /// Add contribution to element's residual vector and Jacobian
222  //=========================================================================
225  Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
226  {
227  // Let's get the info from the parent
228  FiniteElement* parent_pt = bulk_element_pt();
229 
230  // Find the dimension of the problem
231  unsigned spatial_dim = this->nodal_dimension();
232 
233  // Outer unit normal to the wall
234  Vector<double> wall_normal(spatial_dim);
235 
236  // Outer unit normal to the free surface
237  Vector<double> unit_normal(spatial_dim);
238 
239  // Find the dimension of the parent
240  unsigned n_dim = parent_pt->dim();
241 
242  // Find the local coordinates in the parent
243  Vector<double> s_parent(n_dim);
244 
245  // Storage for the shape functions
246  unsigned n_node = this->nnode();
247  Shape psi(n_node);
248  DShape dpsids(n_node, 1);
249  Vector<double> s_local(1);
250 
251  // Loop over intergration points
252  unsigned n_intpt = this->integral_pt()->nweight();
253  for (unsigned ipt = 0; ipt < n_intpt; ++ipt)
254  {
255  // Get the local coordinate of the integration point
256  s_local[0] = this->integral_pt()->knot(ipt, 0);
257  get_local_coordinate_in_bulk(s_local, s_parent);
258 
259  // Get the local shape functions
260  this->dshape_local(s_local, psi, dpsids);
261 
262  // Zero the position
263  Vector<double> x(spatial_dim, 0.0);
264 
265  // Now construct the position and the tangent
266  Vector<double> interpolated_t1(spatial_dim, 0.0);
267  for (unsigned n = 0; n < n_node; n++)
268  {
269  const double psi_local = psi(n);
270  const double dpsi_local = dpsids(n, 0);
271  for (unsigned i = 0; i < spatial_dim; i++)
272  {
273  double pos = this->nodal_position(n, i);
274  interpolated_t1[i] += pos * dpsi_local;
275  x[i] += pos * psi_local;
276  }
277  }
278 
279  // Now we can calculate the Jacobian term
280  double t_length = 0.0;
281  for (unsigned i = 0; i < spatial_dim; ++i)
282  {
283  t_length += interpolated_t1[i] * interpolated_t1[i];
284  }
285  double W = std::sqrt(t_length) * this->integral_pt()->weight(ipt);
286 
287  // Imposition of contact angle in weak form
288  if (Contact_angle_flag == 2)
289  {
290  // Get the outer unit normal of the entire interface
291  dynamic_cast<FaceElement*>(parent_pt)->outer_unit_normal(s_parent,
292  unit_normal);
293 
294  // Calculate the wall normal
295  wall_unit_normal(x, wall_normal);
296 
297  // Find the dot product of the two
298  double dot = 0.0;
299  for (unsigned i = 0; i < spatial_dim; i++)
300  {
301  dot += unit_normal[i] * wall_normal[i];
302  }
303 
304  // Find the projection of the outer normal of the surface into the plane
305  Vector<double> binorm(spatial_dim);
306  for (unsigned i = 0; i < spatial_dim; i++)
307  {
308  binorm[i] = unit_normal[i] - dot * wall_normal[i];
309  }
310 
311  // Get the value of sigma from the parent
312  const double sigma_local =
313  dynamic_cast<FluidInterfaceElement*>(parent_pt)->sigma(s_parent);
314 
315  // Get the capillary number
316  const double ca_local = ca();
317 
318  // Get the contact angle
319  const double theta = contact_angle();
320 
321  // Add the contributions to the momentum equation
322 
323  // Loop over the shape functions
324  for (unsigned l = 0; l < n_node; l++)
325  {
326  // Loop over the velocity components
327  for (unsigned i = 0; i < 3; i++)
328  {
329  // Get the equation number for the momentum equation
330  int local_eqn =
331  this->nodal_local_eqn(l, this->U_index_interface_boundary[i]);
332 
333  // If it's not a boundary condition
334  if (local_eqn >= 0)
335  {
336  // Add the surface-tension contribution to the momentum equation
337  residuals[local_eqn] +=
338  (sigma_local / ca_local) *
339  (sin(theta) * wall_normal[i] + cos(theta) * binorm[i]) *
340  psi(l) * W;
341  }
342  }
343  }
344  }
345  // Otherwise [strong imposition (by hijacking) of contact angle or
346  // "no constraint at all"], add the appropriate contribution to
347  // the momentum equation
348  else
349  {
350  // Storage for the outer vector
351  Vector<double> m(3);
352 
353  // Get the outer unit normal of the line
354  this->outer_unit_normal(s_local, m);
355 
356  // Get the value of sigma from the parent
357  const double sigma_local =
358  dynamic_cast<FluidInterfaceElement*>(parent_pt)->sigma(s_parent);
359 
360  // Get the capillary number
361  const double ca_local = ca();
362 
363  // Loop over the shape functions
364  for (unsigned l = 0; l < n_node; l++)
365  {
366  // Loop over the velocity components
367  for (unsigned i = 0; i < 3; i++)
368  {
369  // Get the equation number for the momentum equation
370  int local_eqn =
371  this->nodal_local_eqn(l, this->U_index_interface_boundary[i]);
372 
373  // If it's not a boundary condition
374  if (local_eqn >= 0)
375  {
376  // Add the surface-tension contribution to the momentum equation
377  residuals[local_eqn] +=
378  m[i] * (sigma_local / ca_local) * psi(l) * W;
379  }
380  }
381  }
382  } // End of the line integral terms
383 
384 
385  // If we are imposing the contact angle strongly (by hijacking)
386  // overwrite the kinematic equation
387  if (Contact_angle_flag == 1)
388  {
389  // Get the outer unit normal of the whole interface
390  dynamic_cast<FaceElement*>(parent_pt)->outer_unit_normal(s_parent,
391  unit_normal);
392 
393  // Calculate the wall normal
394  wall_unit_normal(x, wall_normal);
395 
396  // Find the dot product
397  double dot = 0.0;
398  for (unsigned i = 0; i < spatial_dim; i++)
399  {
400  dot += unit_normal[i] * wall_normal[i];
401  }
402 
403  // Loop over the test functions
404  for (unsigned l = 0; l < n_node; l++)
405  {
406  // Read out the kinematic equation number
407  int local_eqn = kinematic_local_eqn(l);
408 
409  // Note that because we have outer unit normals for the free surface
410  // and the wall, the cosine of the contact angle is equal to
411  // MINUS the dot product
412  if (local_eqn >= 0)
413  {
414  residuals[local_eqn] += (cos(contact_angle()) + dot) * psi(l) * W;
415  }
416  // NOTE: The jacobian entries will be computed automatically
417  // by finite differences.
418  }
419  } // End of strong form of contact angle condition
420 
421  // Add any additional residual contributions
423  residuals, jacobian, flag, psi, dpsids, unit_normal, W);
424  }
425  }
426 
427 
428  /// //////////////////////////////////////////////////////////////
429  /// //////////////////////////////////////////////////////////////
430  /// //////////////////////////////////////////////////////////////
431 
432 
433  //============================================================
434  /// Default value for physical constant (static)
435  //============================================================
437 
438 
439  //================================================================
440  /// Calculate the i-th velocity component at local coordinate s
441  //================================================================
442  double FluidInterfaceElement::interpolated_u(const Vector<double>& s,
443  const unsigned& i)
444  {
445  // Find number of nodes
446  unsigned n_node = FiniteElement::nnode();
447 
448  // Storage for the local shape function
449  Shape psi(n_node);
450 
451  // Get values of shape function at local coordinate s
452  this->shape(s, psi);
453 
454  // Initialise value of u
455  double interpolated_u = 0.0;
456 
457  // Loop over the local nodes and sum
458  for (unsigned l = 0; l < n_node; l++)
459  {
460  interpolated_u += u(l, i) * psi(l);
461  }
462 
463  return (interpolated_u);
464  }
465 
466  //===========================================================================
467  /// Calculate the contribution to the residuals from the interface
468  /// implemented generically with geometric information to be
469  /// added from the specific elements
470  //========================================================================
472  Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
473  {
474  // Find out how many nodes there are
475  unsigned n_node = this->nnode();
476 
477  // Set up memeory for the shape functions
478  Shape psif(n_node);
479 
480  // Find out the number of surface coordinates
481  const unsigned el_dim = this->dim();
482  // We have el_dim local surface coordinates
483  DShape dpsifds(n_node, el_dim);
484 
485  // Find the dimension of the problem
486  // Note that this will return 2 for the axisymmetric case
487  // which is correct in the sense that no
488  // terms will be added in the azimuthal direction
489  const unsigned n_dim = this->nodal_dimension();
490 
491  // Storage for the surface gradient
492  // and divergence terms
493  // These will actually be identical
494  // apart from in the axisymmetric case
495  DShape dpsifdS(n_node, n_dim);
496  DShape dpsifdS_div(n_node, n_dim);
497 
498  // Set the value of n_intpt
499  unsigned n_intpt = this->integral_pt()->nweight();
500 
501  // Get the value of the Capillary number
502  double Ca = ca();
503 
504  // Get the value of the Strouhal numer
505  double St = st();
506 
507  // Get the value of the external pressure
508  double p_ext = pext();
509 
510  // Integers used to hold the local equation numbers and local unknowns
511  int local_eqn = 0, local_unknown = 0;
512 
513  // Storage for the local coordinate
514  Vector<double> s(el_dim);
515 
516  // Loop over the integration points
517  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
518  {
519  // Get the value of the local coordiantes at the integration point
520  for (unsigned i = 0; i < el_dim; i++)
521  {
522  s[i] = this->integral_pt()->knot(ipt, i);
523  }
524 
525  // Get the integral weight
526  double W = this->integral_pt()->weight(ipt);
527 
528  // Call the derivatives of the shape function
529  this->dshape_local_at_knot(ipt, psif, dpsifds);
530 
531  // Define and zero the tangent Vectors and local velocities
532  Vector<double> interpolated_x(n_dim, 0.0);
533  Vector<double> interpolated_u(n_dim, 0.0);
534  Vector<double> interpolated_dx_dt(n_dim, 0.0);
535  ;
536  DenseMatrix<double> interpolated_t(el_dim, n_dim, 0.0);
537 
538  // Loop over the shape functions
539  for (unsigned l = 0; l < n_node; l++)
540  {
541  const double psi_ = psif(l);
542  // Loop over directional components
543  for (unsigned i = 0; i < n_dim; i++)
544  {
545  // Coordinate
546  interpolated_x[i] += this->nodal_position(l, i) * psi_;
547 
548  // Calculate velocity of mesh
549  interpolated_dx_dt[i] += this->dnodal_position_dt(l, i) * psi_;
550 
551  // Calculate the tangent vectors
552  for (unsigned j = 0; j < el_dim; j++)
553  {
554  interpolated_t(j, i) += this->nodal_position(l, i) * dpsifds(l, j);
555  }
556 
557  // Calculate velocity and tangent vector
558  interpolated_u[i] += u(l, i) * psi_;
559  }
560  }
561 
562 
563  // Calculate the surface gradient and divergence
564  double J = this->compute_surface_derivatives(
565  psif, dpsifds, interpolated_t, interpolated_x, dpsifdS, dpsifdS_div);
566  // Get the normal vector
567  //(ALH: This could be more efficient because
568  // there is some recomputation of shape functions here)
569  Vector<double> interpolated_n(n_dim);
570  this->outer_unit_normal(s, interpolated_n);
571 
572  // Now also get the (possible variable) surface tension
573  double Sigma = this->sigma(s);
574 
575  // Loop over the shape functions
576  for (unsigned l = 0; l < n_node; l++)
577  {
578  // Loop over the velocity components
579  for (unsigned i = 0; i < n_dim; i++)
580  {
581  // Get the equation number for the momentum equation
582  local_eqn = this->nodal_local_eqn(l, this->U_index_interface[i]);
583 
584  // If it's not a boundary condition
585  if (local_eqn >= 0)
586  {
587  // Add the surface-tension contribution to the momentum equation
588  residuals[local_eqn] -= (Sigma / Ca) * dpsifdS_div(l, i) * J * W;
589 
590  // If the element is a free surface, add in the external pressure
591  if (Pext_data_pt != 0)
592  {
593  // External pressure term
594  residuals[local_eqn] -=
595  p_ext * interpolated_n[i] * psif(l) * J * W;
596 
597  // Add in the Jacobian term for the external pressure
598  // The correct area is included in the length of the normal
599  // vector
600  if (flag)
601  {
602  local_unknown = pext_local_eqn();
603  if (local_unknown >= 0)
604  {
605  jacobian(local_eqn, local_unknown) -=
606  interpolated_n[i] * psif(l) * J * W;
607  }
608  }
609  } // End of pressure contribution
610  }
611  } // End of contribution to momentum equation
612 
613 
614  // Kinematic BC
615  local_eqn = kinematic_local_eqn(l);
616  if (local_eqn >= 0)
617  {
618  // Assemble the kinematic condition
619  // The correct area is included in the normal vector
620  for (unsigned k = 0; k < n_dim; k++)
621  {
622  residuals[local_eqn] +=
623  (interpolated_u[k] - St * interpolated_dx_dt[k]) *
624  interpolated_n[k] * psif(l) * J * W;
625  }
626 
627  // Add in the jacobian
628  if (flag)
629  {
630  // Loop over shape functions
631  for (unsigned l2 = 0; l2 < n_node; l2++)
632  {
633  // Loop over the components
634  for (unsigned i2 = 0; i2 < n_dim; i2++)
635  {
636  local_unknown =
637  this->nodal_local_eqn(l2, this->U_index_interface[i2]);
638  // If it's a non-zero dof add
639  if (local_unknown >= 0)
640  {
641  jacobian(local_eqn, local_unknown) +=
642  psif(l2) * interpolated_n[i2] * psif(l) * J * W;
643  }
644  }
645  }
646  } // End of Jacobian contribution
647  }
648  } // End of loop over shape functions
649 
650 
651  // Add additional contribution required from the implementation
652  // of the node update (e.g. Lagrange multpliers etc)
654  jacobian,
655  flag,
656  psif,
657  dpsifds,
658  dpsifdS,
659  dpsifdS_div,
660  s,
661  interpolated_x,
662  interpolated_n,
663  W,
664  J);
665 
666  } // End of loop over integration points
667  }
668 
669  //========================================================
670  /// Overload the output functions generically
671  //=======================================================
672  void FluidInterfaceElement::output(std::ostream& outfile,
673  const unsigned& n_plot)
674  {
675  const unsigned el_dim = this->dim();
676  const unsigned n_dim = this->nodal_dimension();
677  const unsigned n_velocity = this->U_index_interface.size();
678  // Set output Vector
679  Vector<double> s(el_dim);
680 
681  // Tecplot header info
682  outfile << tecplot_zone_string(n_plot);
683 
684  // Loop over plot points
685  unsigned num_plot_points = nplot_points(n_plot);
686  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
687  {
688  // Get local coordinates of pliot point
689  get_s_plot(iplot, n_plot, s);
690 
691  // Output the x,y,u,v
692  for (unsigned i = 0; i < n_dim; i++)
693  outfile << this->interpolated_x(s, i) << " ";
694  for (unsigned i = 0; i < n_velocity; i++)
695  outfile << interpolated_u(s, i) << " ";
696 
697  // Output a dummy pressure
698  outfile << 0.0 << "\n";
699  }
700 
701  write_tecplot_zone_footer(outfile, n_plot);
702 
703  outfile << "\n";
704  }
705 
706 
707  //===========================================================================
708  /// Overload the output function
709  //===========================================================================
710  void FluidInterfaceElement::output(FILE* file_pt, const unsigned& n_plot)
711  {
712  const unsigned el_dim = this->dim();
713  const unsigned n_dim = this->nodal_dimension();
714  const unsigned n_velocity = this->U_index_interface.size();
715  // Set output Vector
716  Vector<double> s(el_dim);
717 
718  // Tecplot header info
719  fprintf(file_pt, "%s", tecplot_zone_string(n_plot).c_str());
720 
721  // Loop over plot points
722  unsigned num_plot_points = nplot_points(n_plot);
723  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
724  {
725  // Get local coordinates of plot point
726  get_s_plot(iplot, n_plot, s);
727 
728  // Coordinates
729  for (unsigned i = 0; i < n_dim; i++)
730  {
731  fprintf(file_pt, "%g ", interpolated_x(s, i));
732  }
733 
734  // Velocities
735  for (unsigned i = 0; i < n_velocity; i++)
736  {
737  fprintf(file_pt, "%g ", interpolated_u(s, i));
738  }
739 
740  // Dummy Pressure
741  fprintf(file_pt, "%g \n", 0.0);
742  }
743  fprintf(file_pt, "\n");
744 
745  // Write tecplot footer (e.g. FE connectivity lists)
746  write_tecplot_zone_footer(file_pt, n_plot);
747  }
748 
749 
750  /// //////////////////////////////////////////////////////////////////////
751  /// //////////////////////////////////////////////////////////////////////
752  /// //////////////////////////////////////////////////////////////////////
753 
754  //====================================================================
755  /// Specialise the surface derivatives for the line interface case
756  //===================================================================
758  const Shape& psi,
759  const DShape& dpsids,
760  const DenseMatrix<double>& interpolated_t,
761  const Vector<double>& interpolated_x,
762  DShape& dpsidS,
763  DShape& dpsidS_div)
764  {
765  const unsigned n_shape = psi.nindex1();
766  const unsigned n_dim = 2;
767 
768  // Calculate the only entry of the surface
769  // metric tensor (square length of the tangent vector)
770  double a11 = interpolated_t(0, 0) * interpolated_t(0, 0) +
771  interpolated_t(0, 1) * interpolated_t(0, 1);
772 
773  // Now set the derivative terms of the shape functions
774  for (unsigned l = 0; l < n_shape; l++)
775  {
776  for (unsigned i = 0; i < n_dim; i++)
777  {
778  dpsidS(l, i) = dpsids(l, 0) * interpolated_t(0, i) / a11;
779  }
780  }
781 
782  // The surface divergence is the same as the surface
783  // gradient operator
784  dpsidS_div = dpsidS;
785 
786  // Return the jacobian
787  return sqrt(a11);
788  }
789 
790 
791  /// /////////////////////////////////////////////////////////////////////
792  /// /////////////////////////////////////////////////////////////////////
793  /// /////////////////////////////////////////////////////////////////////
794 
795  //====================================================================
796  /// Specialise the surface derivatives for the axisymmetric interface case
797  //===================================================================
799  const Shape& psi,
800  const DShape& dpsids,
801  const DenseMatrix<double>& interpolated_t,
802  const Vector<double>& interpolated_x,
803  DShape& dpsidS,
804  DShape& dpsidS_div)
805  {
806  // Initially the same as the 2D case
807  const unsigned n_shape = psi.nindex1();
808  const unsigned n_dim = 2;
809 
810  // Calculate the only entry of the surface
811  // metric tensor (square length of the tangent vector)
812  double a11 = interpolated_t(0, 0) * interpolated_t(0, 0) +
813  interpolated_t(0, 1) * interpolated_t(0, 1);
814 
815  // Now set the derivative terms of the shape functions
816  for (unsigned l = 0; l < n_shape; l++)
817  {
818  for (unsigned i = 0; i < n_dim; i++)
819  {
820  dpsidS(l, i) = dpsids(l, 0) * interpolated_t(0, i) / a11;
821  // Set the standard components of the divergence
822  dpsidS_div(l, i) = dpsidS(l, i);
823  }
824  }
825 
826  const double r = interpolated_x[0];
827 
828  // The surface divergence is different because we need
829  // to take account of variation of the base vectors
830  for (unsigned l = 0; l < n_shape; l++)
831  {
832  dpsidS_div(l, 0) += psi(l) / r;
833  }
834 
835  // Return the jacobian, needs to be multiplied by r
836  return r * sqrt(a11);
837  }
838 
839 
840  /// /////////////////////////////////////////////////////////////////////
841  /// /////////////////////////////////////////////////////////////////////
842  /// /////////////////////////////////////////////////////////////////////
843 
844  //====================================================================
845  /// Specialise the surface derivatives for 2D surface case
846  //===================================================================
848  const Shape& psi,
849  const DShape& dpsids,
850  const DenseMatrix<double>& interpolated_t,
851  const Vector<double>& interpolated_x,
852  DShape& dpsidS,
853  DShape& dpsidS_div)
854  {
855  const unsigned n_shape = psi.nindex1();
856  const unsigned n_dim = 3;
857 
858  // Calculate the local metric tensor
859  // The dot product of the two tangent vectors
860  double amet[2][2];
861  for (unsigned al = 0; al < 2; al++)
862  {
863  for (unsigned be = 0; be < 2; be++)
864  {
865  // Initialise to zero
866  amet[al][be] = 0.0;
867  // Add the dot product contributions
868  for (unsigned i = 0; i < n_dim; i++)
869  {
870  amet[al][be] += interpolated_t(al, i) * interpolated_t(be, i);
871  }
872  }
873  }
874 
875  // Work out the determinant
876  double det_a = amet[0][0] * amet[1][1] - amet[0][1] * amet[1][0];
877  // Also work out the contravariant metric
878  double aup[2][2];
879  aup[0][0] = amet[1][1] / det_a;
880  aup[0][1] = -amet[0][1] / det_a;
881  aup[1][0] = -amet[1][0] / det_a;
882  aup[1][1] = amet[0][0] / det_a;
883 
884 
885  // Now construct the surface gradient terms
886  for (unsigned l = 0; l < n_shape; l++)
887  {
888  // Do some pre-computation
889  const double dpsi_temp[2] = {
890  aup[0][0] * dpsids(l, 0) + aup[0][1] * dpsids(l, 1),
891  aup[1][0] * dpsids(l, 0) + aup[1][1] * dpsids(l, 1)};
892 
893  for (unsigned i = 0; i < n_dim; i++)
894  {
895  dpsidS(l, i) = dpsi_temp[0] * interpolated_t(0, i) +
896  dpsi_temp[1] * interpolated_t(1, i);
897  }
898  }
899 
900  // The divergence operator is the same
901  dpsidS_div = dpsidS;
902 
903  // Return the jacobian
904  return sqrt(det_a);
905  }
906 
907 } // namespace oomph
double compute_surface_derivatives(const Shape &psi, const DShape &dpsids, const DenseMatrix< double > &interpolated_t, const Vector< double > &interpolated_x, DShape &surface_gradient, DShape &surface_divergence)
Fill in the specific surface derivative calculations.
Vector< unsigned > U_index_interface_boundary
Index at which the i-th velocity component is stored in the element's nodes.
virtual int kinematic_local_eqn(const unsigned &n)=0
Function that is used to determine the local equation number of the kinematic equation associated wit...
void wall_unit_normal(const Vector< double > &x, Vector< double > &normal)
Function that returns the unit normal of the bounding wall directed out of the fluid.
double * Contact_angle_pt
Pointer to the desired value of the contact angle (if any)
virtual void add_additional_residual_contributions_interface_boundary(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag, const Shape &psif, const DShape &dpsifds, const Vector< double > &interpolated_n, const double &W)
Empty helper function to calculate the additional contributions arising from the node update strategy...
double & contact_angle()
Return value of the contact angle.
void set_contact_angle(double *const &angle_pt, const bool &strong=true)
Set a pointer to the desired contact angle. Optional boolean (defaults to true) chooses strong imposi...
double ca()
Return the value of the capillary number.
unsigned Contact_angle_flag
Flag used to determine whether the contact angle is to be used (0 if not), and whether it will be app...
Base class establishing common interfaces and functions for all Navier-Stokes-like fluid interface el...
virtual double compute_surface_derivatives(const Shape &psi, const DShape &dpsids, const DenseMatrix< double > &interpolated_t, const Vector< double > &interpolated_x, DShape &dpsidS, DShape &dpsidS_div)=0
Compute the surface gradient and surface divergence operators given the shape functions,...
virtual 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)
Helper function to calculate the additional contributions to the resisuals and Jacobian that arise fr...
const double & st() const
The value of the Strouhal number.
double u(const unsigned &j, const unsigned &i)
Return the i-th velocity component at local node j.
virtual int kinematic_local_eqn(const unsigned &n)=0
Access function that returns the local equation number for the (scalar) kinematic equation associated...
Vector< unsigned > U_index_interface
Nodal index at which the i-th velocity component is stored.
double pext() const
Return the value of the external pressure.
virtual double sigma(const Vector< double > &s_local)
Virtual function that specifies the non-dimensional surface tension as a function of local position w...
double interpolated_u(const Vector< double > &s, const unsigned &i)
Calculate the i-th velocity component at the local coordinate s.
Data * Pext_data_pt
Pointer to the Data item that stores the external pressure.
void output(std::ostream &outfile)
Overload the output function.
const double & ca() const
The value of the Capillary number.
int pext_local_eqn()
Access function for the local equation number that corresponds to the external pressure.
static double Default_Physical_Constant_Value
Default value for physical constants.
virtual void fill_in_generic_residual_contribution_interface(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Helper function to calculate the residuals and (if flag==1) the Jacobian of the equations....
double compute_surface_derivatives(const Shape &psi, const DShape &dpsids, const DenseMatrix< double > &interpolated_t, const Vector< double > &interpolated_x, DShape &surface_gradient, DShape &surface_divergence)
Fill in the specific surface derivative calculations.
void fill_in_generic_residual_contribution_interface_boundary(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Overload the helper function to calculate the residuals and (if flag==true) the Jacobian – this funct...
void fill_in_generic_residual_contribution_interface_boundary(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Overload the helper function to calculate the residuals and (if flag==1) the Jacobian – this function...
double compute_surface_derivatives(const Shape &psi, const DShape &dpsids, const DenseMatrix< double > &interpolated_t, const Vector< double > &interpolated_x, DShape &surface_gradient, DShape &surface_divergence)
Fill in the specific surface derivative calculations.