young_laplace_contact_angle_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//====================================================================
31 #include "young_laplace_elements.h"
33 
34 
35 namespace oomph
36 {
37 
38 //===========================================================================
39 /// Constructor, takes the pointer to the "bulk" element and the
40 /// index of the face to which the element is attached.
41 //===========================================================================
42  template<class ELEMENT>
45  const int &face_index) :
46  FaceGeometry<ELEMENT>(), FaceElement()
47  {
48 
49  // Let the bulk element build the FaceElement, i.e. setup the pointers
50  // to its nodes (by referring to the appropriate nodes in the bulk
51  // element), etc.
52  bulk_el_pt->build_face_element(face_index,this);
53 
54  // Initialise the prescribed contact angle pointer to zero
56 
57 #ifdef PARANOID
58  // Extract the dimension of the problem from the dimension of
59  // the first node
60  unsigned dim = node_pt(0)->ndim();
61  if (dim!=2)
62  {
63  throw OomphLibError(
64  "YoungLaplaceContactAngleElement only work with 2D nodes",
65  OOMPH_CURRENT_FUNCTION,
66  OOMPH_EXCEPTION_LOCATION);
67  }
68 #endif
69 
70  }
71 
72 
73 
74 
75 
76 //================================================================
77 /// Compute the element's contribution to the residual vector
78 //================================================================
79  template<class ELEMENT>
82  {
83  //Find out how many nodes there are
84  unsigned n_node = nnode();
85 
86  //Set up memory for the shape functions
87  Shape psi(n_node);
88 
89  //Number of integration points
90  unsigned n_intpt = integral_pt()->nweight();
91 
92  // Set the Vector to hold local coordinates
93  // Note: We need the coordinate itself below for the evaluation
94  // of the contact line vectors even though we use the *at_knot
95  // version for the various shape-function-related functions
96  Vector<double> s(1);
97 
98  //Integers to hold the local equation and unknown numbers
99  int local_eqn=0;
100 
101  //Loop over the integration points
102  for(unsigned ipt=0;ipt<n_intpt;ipt++)
103  {
104  //Assign value of s
105  s[0] = integral_pt()->knot(ipt,0);
106 
107  //Get the integral weight
108  double w = integral_pt()->weight(ipt);
109 
110  //Find the shape functions
111  shape_at_knot(ipt,psi);
112 
113  //Get the prescribed cos_gamma
114  double cos_gamma=prescribed_cos_gamma();
115 
116  // Get the various contact line vectors
117  Vector<double> tangent(3);
118  Vector<double> normal(3);
119  Vector<double> spine(3);
120  double norm_of_drds;
121  contact_line_vectors(s,tangent,normal,spine,norm_of_drds);
122 
123  // Get beta factor:
124 
125  // Cross product of spine and tangent to contact line is
126  // the wall normal
127  Vector<double> wall_normal(3);
128  ELEMENT::cross_product(spine,tangent,wall_normal);
129 
130  // Normalise
131  double norm=ELEMENT::two_norm(wall_normal);
132  for (unsigned i=0;i<3;i++) wall_normal[i]/=norm;
133 
134  // Take cross product with tangent to get the normal to the
135  // contact line parallel to wall
136  Vector<double> normal_to_contact_line_parallel_to_wall(3);
137  ELEMENT::cross_product(tangent,wall_normal,
138  normal_to_contact_line_parallel_to_wall);
139 
140  double beta=ELEMENT::scalar_product(spine,
141  normal_to_contact_line_parallel_to_wall);
142 
143 
144  //Now add to the appropriate equations
145 
146  //Loop over the test functions
147  for(unsigned l=0;l<n_node;l++)
148  {
149  local_eqn = u_local_eqn(l);
150  /*IF it's not a boundary condition*/
151  if(local_eqn >= 0)
152  {
153  //Add to residual:
154  residuals[local_eqn] -= beta*cos_gamma*psi[l]*norm_of_drds*w;
155  }
156  }
157  }
158 
159 }
160 
161 
162 
163 
164 //========================================================================
165 /// Get the actual contact angle
166 //========================================================================
167  template<class ELEMENT>
169  const Vector<double>& s)
170  {
171  // Get pointer to bulk element
172  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
173 
174  double cos_gamma=0.0;
175 
176  // Spine case
177  if ( bulk_elem_pt->use_spines() )
178  {
179  // Get various contact line vectors
180  Vector<double> tangent(3);
181  Vector<double> normal(3);
182  Vector<double> spine(3);
183  double norm_of_drds;
184  contact_line_vectors(s,tangent,normal,spine,norm_of_drds);
185 
186  // Get the wall normal: Both the tangent to the contact line
187  // and the spine vector are tangential to the wall:
188  Vector<double> wall_normal(3);
189  Vector<double> axe_ez(3,0.0);
190  axe_ez[2]=1.0;
191  ELEMENT::cross_product(axe_ez,tangent,wall_normal);
192 
193  // Normalise
194  double norm=0.0;
195  for (unsigned i=0;i<3;i++) norm+=wall_normal[i]*wall_normal[i];
196  for (unsigned i=0;i<3;i++)
197  {
198  wall_normal[i]/=sqrt(norm);
199  }
200 
201  // Get the auxiliary unit vector that's normal to
202  // the contact line and tangent to the wall
203  Vector<double> aux(3);
204  ELEMENT::cross_product(tangent,wall_normal,aux);
205 
206  // Cos of contact angle is dot product with wall normal
207  cos_gamma=ELEMENT::scalar_product(aux,normal);
208 
209  }
210 
211  // Cartesian case
212  else
213  {
214 
215  // Get local coordinates in bulk element by copy construction
216  Vector<double> s_bulk(local_coordinate_in_bulk(s));
217 
218  // Number of nodes in bulk element
219  unsigned nnode_bulk=bulk_elem_pt->nnode();
220 
221 
222 #ifdef PARANOID
223  // Dimension of (= number of local coordinates in) bulk element
224  unsigned dim_bulk=bulk_elem_pt->dim();
225 
226  if (dim_bulk!=2)
227  {
228  throw OomphLibError(
229  "YoungLaplaceContactAngleElements only work with 2D bulk elements",
230  OOMPH_CURRENT_FUNCTION,
231  OOMPH_EXCEPTION_LOCATION);
232  }
233 #endif
234 
235  //Set up memory for the shape functions
236  Shape psi(nnode_bulk);
237  DShape dpsidzeta(nnode_bulk,2);
238 
239  //Call the derivatives of the shape and test functions
240  //in the bulk -- must do this via s because this point
241  //is not an integration point the bulk element!
242  (void)bulk_elem_pt->dshape_eulerian(s_bulk,psi,dpsidzeta);
243 
244  // Get the gradient at s
245  Vector<double> gradient_u(2,0.0);
246 
247  //Calculate function value and derivatives:
248  //-----------------------------------------
249  // Loop over nodes
250  for(unsigned l=0;l<nnode_bulk;l++)
251  {
252  // Loop over directions
253  for(unsigned j=0;j<2;j++)
254  { gradient_u[j]+= bulk_elem_pt->u(l)*dpsidzeta(l,j); }
255  }
256 
257  // Get the outer unit normal to boundary
258  Vector<double> outer_normal(2,0.0);
259  outer_unit_normal(s,outer_normal);
260 
261  // Compute the cosinus of the angle
262  double gradient_norm_2=
263  ELEMENT::two_norm(gradient_u)*
264  ELEMENT::two_norm(gradient_u);
265  cos_gamma=ELEMENT::scalar_product(gradient_u,outer_normal)/
266  sqrt(1+gradient_norm_2);
267 
268  }
269 
270  return cos_gamma;
271 
272 }
273 
274 
275 
276 //========================================================================
277 /// Get unit tangent and normal to contact line and the spine itself (this
278 /// allows the wall normal to be worked out by a cross product). Final
279 /// argument gives the norm of dR/ds where R is the vector to the
280 /// contact line and s the local coordinate in the element.
281 //========================================================================
282 template<class ELEMENT>
284  const Vector<double>& s,
285  Vector<double>& tangent,
286  Vector<double>& normal,
287  Vector<double>& spine,
288  double& norm_of_drds)
289 {
290 
291  // Get pointer to bulk element
292  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
293 
294  // Dimension of (= number of local coordinates in) bulk element
295  unsigned dim_bulk=bulk_elem_pt->dim();
296 
297 #ifdef PARANOID
298  if (dim_bulk!=2)
299  {
300  throw OomphLibError(
301  "YoungLaplaceContactAngleElements only work with 2D bulk elements",
302  OOMPH_CURRENT_FUNCTION,
303  OOMPH_EXCEPTION_LOCATION);
304  }
305 #endif
306 
307  // Which local coordinate is const in bulk element?
308  unsigned s_fixed_index_in_bulk;
309 
310  DenseMatrix<double> ds_bulk_ds_face(2,1);
311  this->get_ds_bulk_ds_face(s,ds_bulk_ds_face,
312  s_fixed_index_in_bulk);
313 
314  // Get local coordinates in bulk element by copy construction
315  Vector<double> s_bulk(local_coordinate_in_bulk(s));
316 
317  // Number of nodes in bulk element
318  unsigned nnode_bulk=bulk_elem_pt->nnode();
319 
320  // Get the shape functions and their derivatives w.r.t.
321  // to the local coordinates in the bulk element
322  Shape psi_bulk(nnode_bulk);
323  DShape dpsi_bulk(nnode_bulk,dim_bulk);
324  bulk_elem_pt->dshape_local(s_bulk,psi_bulk,dpsi_bulk);
325 
326  // Displacement along spine
327  double interpolated_u=0.0;
328 
329  // Derivative of u w.r.t. tangential and pseudo-normal bulk coordinate
330  double interpolated_du_ds_tangent=0.0;
331  double interpolated_du_ds_pseudo_normal=0.0;
332 
333  // Global intrinsic coordinates of current point for evaluation of
334  // spine and spine base
335  Vector<double> interpolated_zeta(dim_bulk,0.0);
336 
337 
338  // Derivatives of zeta (in bulk) w.r.t. to tangent and pseudo-normal
339  // local coordinate
340  Vector<double> interpolated_dzeta_ds_tangent(dim_bulk);
341  Vector<double> interpolated_dzeta_ds_pseudo_normal(dim_bulk);
342 
343  // Loop over nodes
344  for(unsigned l=0;l<nnode_bulk;l++)
345  {
346  interpolated_u+=bulk_elem_pt->u(l)*psi_bulk(l);
347 
348  // Chain rule for tangent derivative
349  double aux=0.0;
350  for (unsigned i=0;i<dim_bulk;i++)
351  {
352  aux+=dpsi_bulk(l,i)*ds_bulk_ds_face(i,0);
353  }
354  interpolated_du_ds_tangent+=bulk_elem_pt->u(l)*aux;
355 
356  // Straight derivative w.r.t. one of the bulk coordinates
357  // that's fixed along the face
358  interpolated_du_ds_pseudo_normal +=
359  bulk_elem_pt->u(l)*dpsi_bulk(l,s_fixed_index_in_bulk);
360 
361  // Loop over directions
362  for(unsigned j=0;j<dim_bulk;j++)
363  {
364  interpolated_zeta[j]+=bulk_elem_pt->nodal_position(l,j)*psi_bulk(l);
365 
366  // Chain rule for tangent derivative
367  double aux=0.0;
368  for (unsigned i=0;i<dim_bulk;i++)
369  {
370  aux+=dpsi_bulk(l,i)*ds_bulk_ds_face(i,0);
371  }
372  interpolated_dzeta_ds_tangent[j] +=
373  bulk_elem_pt->nodal_position(l,j)*aux;
374 
375  // Straight derivative w.r.t. one of the bulk coordinates
376  // that's fixed along the face
377  interpolated_dzeta_ds_pseudo_normal[j] +=
378  bulk_elem_pt->nodal_position(l,j)*dpsi_bulk(l,s_fixed_index_in_bulk);
379 
380  }
381  }
382 
383  // Auxiliary vector (tangent to non-fixed bulk coordinate but
384  // not necessarily normal to contact line)
385  Vector<double> aux_vector(3);
386  double tang_norm=0.0;
387  double aux_norm=0.0;
388 
389  if ( bulk_elem_pt->use_spines() )
390  {
391  // Get the spine and spine base vector at this point from the bulk element
392  Vector<double> spine_base(3,0.0);
393  Vector< Vector<double> > dspine;
394  ELEMENT::allocate_vector_of_vectors(2,3,dspine);
395  Vector< Vector<double> > dspine_base;
396  ELEMENT::allocate_vector_of_vectors(2,3,dspine_base);
397  bulk_elem_pt->get_spine(interpolated_zeta,spine,dspine);
398  bulk_elem_pt->get_spine_base(interpolated_zeta,spine_base,dspine_base);
399 
400  // Derivative of spine and spine base w.r.t. local coordinate in
401  // FaceElement:
402  Vector<double> dspine_ds_tangent(3,0.0);
403  Vector<double> dspine_base_ds_tangent(3,0.0);
404  Vector<double> dspine_ds_pseudo_normal(3,0.0);
405  Vector<double> dspine_base_ds_pseudo_normal(3,0.0);
406  for (unsigned i=0;i<3;i++)
407  {
408 
409  dspine_ds_tangent[i]+=
410  dspine[0][i]*interpolated_dzeta_ds_tangent[0]+
411  dspine[1][i]*interpolated_dzeta_ds_tangent[1];
412 
413  dspine_base_ds_tangent[i]+=
414  dspine_base[0][i]*interpolated_dzeta_ds_tangent[0]+
415  dspine_base[1][i]*interpolated_dzeta_ds_tangent[1];
416 
417  dspine_ds_pseudo_normal[i]+=
418  dspine[0][i]*interpolated_dzeta_ds_pseudo_normal[0]+
419  dspine[1][i]*interpolated_dzeta_ds_pseudo_normal[1];
420 
421  dspine_base_ds_pseudo_normal[i]+=
422  dspine_base[0][i]*interpolated_dzeta_ds_pseudo_normal[0]+
423  dspine_base[1][i]*interpolated_dzeta_ds_pseudo_normal[1];
424  }
425 
426  // Auxiliary vector (tangent to non-fixed bulk coordinate but
427  // not necessarily normal to contact line)
428  for (unsigned i=0;i<3;i++)
429  {
430  tangent[i]=
431  dspine_base_ds_tangent[i]+interpolated_du_ds_tangent*spine[i]+
432  interpolated_u*dspine_ds_tangent[i];
433  tang_norm+=tangent[i]*tangent[i];
434 
435 
436  aux_vector[i]=
437  dspine_base_ds_pseudo_normal[i]+
438  interpolated_du_ds_pseudo_normal*spine[i]+
439  interpolated_u*dspine_ds_pseudo_normal[i];
440  aux_norm+=aux_vector[i]*aux_vector[i];
441  }
442 
443  }
444 
445  // Cartesian case
446  else
447  {
448  for(unsigned i=0;i<2;i++)
449  {
450  tangent[i]= interpolated_dzeta_ds_tangent[i];
451  tang_norm+=tangent[i]*tangent[i];
452  aux_vector[i]= interpolated_dzeta_ds_pseudo_normal[i];
453  aux_norm+=aux_vector[i]*aux_vector[i];
454 
455  }
456 
457  tangent[2]=interpolated_du_ds_tangent;
458  tang_norm+=tangent[2]*tangent[2];
459  aux_vector[2]=interpolated_du_ds_pseudo_normal;
460  aux_norm+=aux_vector[2]*aux_vector[2];
461 
462  }
463 
464  // Norm of the rate of change
465  norm_of_drds=sqrt(tang_norm);
466 
467  // Normalise
468  double tang_norm_fact=1.0/sqrt(tang_norm);
469  double aux_norm_fact=1.0/sqrt(aux_norm);
470  for (unsigned i=0;i<3;i++)
471  {
472  tangent[i]*=tang_norm_fact;
473  aux_vector[i]*=aux_norm_fact;
474  }
475 
476 
477  // Normal to meniscus is the cross product between the
478  // two contact line vectors:
479  Vector<double> meniscus_normal(3);
480  ELEMENT::cross_product(tangent,aux_vector,meniscus_normal);
481 
482  // Calculate the norm
483  double men_norm_fact=0.0;
484  for (unsigned i=0;i<3;i++)
485  {
486  men_norm_fact+=meniscus_normal[i]*meniscus_normal[i];
487  }
488 
489  // Normalise and adjust direction
490  double sign=-double(normal_sign());
491  for (unsigned i=0;i<3;i++)
492  {
493  meniscus_normal[i]*=sign/sqrt(men_norm_fact);
494  }
495 
496  // The actual (bi) normal to the contact line is given
497  // by another cross product
498  ELEMENT::cross_product(meniscus_normal,tangent,normal);
499 
500 }
501 
502 
503 
504 
505 
506 //============================================================
507 // Build the required elements
508 //============================================================
512 
516 
517 }
cstr elem_len * i
Definition: cfortran.h:607
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
A general Finite Element class.
Definition: elements.h:1271
A class for elements that allow the imposition of an contact angle bcs for Young Laplace elements...
static char t char * s
Definition: cfortran.h:572
double * Prescribed_cos_gamma_pt
Pointer to prescribed cos gamma.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
void contact_line_vectors(const Vector< double > &s, Vector< double > &tangent, Vector< double > &normal)
Get unit tangent and normal vectors to contact line.
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
double actual_cos_contact_angle(const Vector< double > &s)
Compute cosinus of actual contact angle.