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