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-2022 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
33namespace 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
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 =
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 =
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 //================================================================
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,
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
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
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.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
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:4403
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:6006
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:4528
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4735
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:6384
A general Finite Element class.
Definition: elements.h:1313
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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:3239
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:3161
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...
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:1432
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
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:3186
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:2317
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:1981
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
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:3174
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:2333
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...
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.
double & contact_angle()
Return value of the contact angle.
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 & ca() const
The value of the Capillary number.
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.
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....
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
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...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
unsigned nindex1() const
Return the range of index 1 of the shape function object.
Definition: shape.h:259
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 dot(const Vector< double > &a, const Vector< double > &b)
Probably not always best/fastest because not optimised for dimension but useful...
Definition: Vector.h:291
//////////////////////////////////////////////////////////////////// ////////////////////////////////...