axisym_fluid_traction_elements.h
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: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 May 2016) $
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 //Header file for elements that are used to integrate fluid tractions
31 //This includes the guts (i.e. equations) because we want to inline them
32 //for faster operation, although it slows down the compilation!
33 #ifndef OOMPH_AXISYM_FLUID_TRACTION_ELEMENTS_HEADER
34 #define OOMPH_AXISYM_FLUID_TRACTION_ELEMENTS_HEADER
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38 #include <oomph-lib-config.h>
39 #endif
40 
41 
42 //OOMPH-LIB headers
43 #include "../generic/shape.h"
44 #include "../generic/elements.h"
45 #include "../generic/element_with_external_element.h"
46 
47 
48 namespace oomph
49 {
50 
51 
52 //=======================================================================
53 /// Namespace containing the zero traction function for axisymmetric
54 /// Navier Stokes traction elements
55 //=======================================================================
56 namespace AxisymmetricNavierStokesTractionElementHelper
57  {
58 
59  //=======================================================================
60  /// Default load function (zero traction)
61  //=======================================================================
62  extern void Zero_traction_fct(const double& time,
63  const Vector<double> &x,
64  const Vector<double>& N,
65  Vector<double>& load);
66 
67  }
68 
69 
70 //======================================================================
71 /// A class for elements that allow the imposition of an applied traction
72 /// in the axisym Navier Stokes eqns.
73 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
74 /// class and and thus, we can be generic enough without the need to have
75 /// a separate equations class.
76 //======================================================================
77 template <class ELEMENT>
79 public virtual FaceGeometry<ELEMENT>,
80  public virtual FaceElement
81  {
82  protected:
83 
84  /// Index at which the i-th velocity component is stored
86 
87  /// \short Pointer to an imposed traction function. Arguments:
88  /// Eulerian coordinate; outer unit normal;
89  /// applied traction. (Not all of the input arguments will be
90  /// required for all specific load functions but the list should
91  /// cover all cases)
92  void (*Traction_fct_pt)(const double& time,
93  const Vector<double> &x,
94  const Vector<double> &n,
95  Vector<double> &result);
96 
97 
98  /// \short Get the traction vector: Pass number of integration point (dummy),
99  /// Eulerian coordinate and normal vector and return the load vector
100  /// (not all of the input arguments will be
101  /// required for all specific load functions but the list should
102  /// cover all cases). This function is virtual so it can be
103  /// overloaded for FSI.
104  virtual void get_traction(const double& time,
105  const unsigned& intpt,
106  const Vector<double>& x,
107  const Vector<double>& n,
108  Vector<double>& traction) const
109  {
110  Traction_fct_pt(time,x,n,traction);
111  }
112 
113 
114  /// \short Helper function that actually calculates the residuals
115  // This small level of indirection is required to avoid calling
116  // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
117  // which causes all kinds of pain if overloading later on
119  Vector<double> &residuals);
120 
121 
122  public:
123 
124  /// \short Constructor, which takes a "bulk" element and the
125  /// value of the index and its limit
127  (FiniteElement* const &element_pt, const int &face_index) :
129  {
130  //Attach the geometrical information to the element. N.B. This function
131  //also assigns nbulk_value from the required_nvalue of the bulk element
132  element_pt->build_face_element(face_index,this);
133 
134  //Find the dimension of the problem
135  unsigned n_dim = element_pt->nodal_dimension();
136 
137  //Find the index at which the velocity unknowns are stored
138  ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
139  this->U_index_axisymmetric_nst_traction.resize(n_dim+1);
140  for(unsigned i=0;i<n_dim+1;i++)
141  {
143  cast_element_pt->u_index_axi_nst(i);
144  }
145 
146  // Zero traction
150  }
151 
152 
153  /// Reference to the traction function pointer
154  void (* &traction_fct_pt())(const double& time,
155  const Vector<double>& x,
156  const Vector<double>& n,
158  {return Traction_fct_pt;}
159 
160 
161  /// Return the residuals
163  {
165  }
166 
167 
168 
169  /// Fill in contribution from Jacobian
171  DenseMatrix<double> &jacobian)
172  {
173  //Call the residuals
175  }
176 
177  /// Specify the value of nodal zeta from the face geometry
178  /// \short The "global" intrinsic coordinate of the element when
179  /// viewed as part of a geometric object should be given by
180  /// the FaceElement representation, by default (needed to break
181  /// indeterminacy if bulk element is SolidElement)
182  double zeta_nodal(const unsigned &n, const unsigned &k,
183  const unsigned &i) const
184  {return FaceElement::zeta_nodal(n,k,i);}
185 
186  /// \short Output function
187  void output(std::ostream &outfile)
188  {
189  unsigned nplot=5;
190  output(outfile,nplot);
191  }
192 
193  /// \short Number of scalars/fields output by this element. Reimplements
194  /// broken virtual function in base class.
195  unsigned nscalar_paraview() const
196  {
197  //Number of dimensions
198  unsigned n_dim = this->nodal_dimension();
199 
200  return 2*(n_dim+1);
201  }
202 
203  /// \short Write values of the k-th scalar field at the plot points. Needs
204  /// to be implemented for each new specific element type.
205  void scalar_value_paraview(std::ofstream& file_out,
206  const unsigned& k,
207  const unsigned& nplot) const
208  {
209  //Number of dimensions
210  unsigned n_dim = this->nodal_dimension();
211 
212  //Find out how many nodes there are
213  const unsigned n_node = nnode();
214 
215  // Get continuous time from timestepper of first node
216  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
217 
218  //Set up memory for the shape functions
219  Shape psi(n_node);
220 
221  // Local and global coordinates
222  Vector<double> s(n_dim-1);
224 
225  // Loop over plot points
226  unsigned num_plot_points=this->nplot_points_paraview(nplot);
227  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
228  {
229  // Get local coordinates of plot point
230  this->get_s_plot(iplot,nplot,s);
231 
232  // Outer unit normal
233  Vector<double> unit_normal(n_dim);
234  outer_unit_normal(s,unit_normal);
235 
236  //Find the shape functions
237  shape(s,psi);
238 
239  //Initialise to zero
240  for(unsigned i=0;i<n_dim;i++)
241  {
242  interpolated_x[i] = 0.0;
243  }
244 
245  //Calculate stuff
246  for(unsigned l=0;l<n_node;l++)
247  {
248  //Loop over directions
249  for(unsigned i=0;i<n_dim;i++)
250  {
251  interpolated_x[i] += this->nodal_position(l,i)*psi[l];
252  }
253  }
254 
255  //Get the imposed traction
257 
258  //Dummy integration point
259  unsigned ipt=0;
260  get_traction(time,ipt,interpolated_x,unit_normal,traction);
261 
262  // Traction components
263  if (k<n_dim+1)
264  {
265  file_out << traction[k] << std::endl;
266  }
267  // Advection Diffusion
268  else if (k<2*n_dim+1 && k>=n_dim+1)
269  {
270  file_out << unit_normal[k] << std::endl;
271  }
272  // Never get here
273  else
274  {
275  std::stringstream error_stream;
276  error_stream
277  << "Axisymmetric Fluid Traction Navier-Stokes Elements only store "
278  << 2*(n_dim+1) << " fields " << std::endl;
279  throw OomphLibError(
280  error_stream.str(),
281  OOMPH_CURRENT_FUNCTION,
282  OOMPH_EXCEPTION_LOCATION);
283  }
284  }
285  }
286 
287  /// \short Name of the i-th scalar field. Default implementation
288  /// returns V1 for the first one, V2 for the second etc. Can (should!) be
289  /// overloaded with more meaningful names in specific elements.
290  std::string scalar_name_paraview(const unsigned& i) const
291  {
292  //Number of dimensions
293  unsigned n_dim = this->nodal_dimension();
294 
295  // Traction components
296  if (i<n_dim+1)
297  {
298  return "Traction component "+StringConversion::to_string(i);
299  }
300  // Normals
301  else if (i<2*n_dim+1 && i>=n_dim+1)
302  {
303  return "Normal "+StringConversion::to_string(i%(n_dim+1));
304  }
305  // Never get here
306  else
307  {
308  std::stringstream error_stream;
309  error_stream
310  << "Axisymmetric Fluid Traction Navier-Stokes Elements only store "
311  << 2*(n_dim+1) << " fields " << std::endl;
312  throw OomphLibError(
313  error_stream.str(),
314  OOMPH_CURRENT_FUNCTION,
315  OOMPH_EXCEPTION_LOCATION);
316  }
317  }
318 
319  /// \short Output function
320  void output(std::ostream &outfile, const unsigned &n_plot)
321  {
322  //Number of dimensions
323  unsigned n_dim = this->nodal_dimension();
324 
325  //Find out how many nodes there are
326  const unsigned n_node = nnode();
327 
328  // Get continuous time from timestepper of first node
329  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
330 
331  //Set up memory for the shape functions
332  Shape psi(n_node);
333 
334  // Local and global coordinates
335  Vector<double> s(n_dim-1);
337 
338  // Tecplot header info
339  outfile << this->tecplot_zone_string(n_plot);
340 
341  // Loop over plot points
342  unsigned num_plot_points=this->nplot_points(n_plot);
343  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
344  {
345  // Get local coordinates of plot point
346  this->get_s_plot(iplot,n_plot,s);
347 
348  // Outer unit normal
349  Vector<double> unit_normal(n_dim);
350  outer_unit_normal(s,unit_normal);
351 
352  //Find the shape functions
353  shape(s,psi);
354 
355  //Initialise to zero
356  for(unsigned i=0;i<n_dim;i++)
357  {
358  interpolated_x[i] = 0.0;
359  }
360 
361  //Calculate stuff
362  for(unsigned l=0;l<n_node;l++)
363  {
364  //Loop over directions
365  for(unsigned i=0;i<n_dim;i++)
366  {
367  interpolated_x[i] += this->nodal_position(l,i)*psi[l];
368  }
369  }
370 
371  //Get the imposed traction
373 
374  //Dummy integration point
375  unsigned ipt=0;
376  get_traction(time,ipt,interpolated_x,unit_normal,traction);
377 
378  //Output the x,y,..
379  for(unsigned i=0;i<n_dim;i++)
380  {
381  outfile << interpolated_x[i] << " ";
382  }
383 
384  //Output the traction components
385  for(unsigned i=0;i<n_dim+1;i++)
386  {
387  outfile << traction[i] << " ";
388  }
389 
390  // Output normal
391  for(unsigned i=0;i<n_dim;i++)
392  {
393  outfile << unit_normal[i] << " ";
394  }
395  outfile << std::endl;
396 
397  }
398  }
399 
400  /// \short C_style output function
401  void output(FILE* file_pt)
402  {FiniteElement::output(file_pt);}
403 
404  /// \short C-style output function
405  void output(FILE* file_pt, const unsigned &n_plot)
406  {FiniteElement::output(file_pt,n_plot);}
407 
408 
409  /// \short Compute traction vector at specified local coordinate
410  /// Should only be used for post-processing; ignores dependence
411  /// on integration point!
412  void traction(const double &time,
413  const Vector<double>& s,
415 
416  };
417 
418 ///////////////////////////////////////////////////////////////////////
419 ///////////////////////////////////////////////////////////////////////
420 ///////////////////////////////////////////////////////////////////////
421 
422 //=====================================================================
423 /// Compute traction vector at specified local coordinate
424 /// Should only be used for post-processing; ignores dependence
425 /// on integration point!
426 //=====================================================================
427 template<class ELEMENT>
429  traction(const double& time, const Vector<double>& s, Vector<double>& traction)
430  {
431  unsigned n_dim = this->nodal_dimension();
432 
433  // Position vector
434  Vector<double> x(n_dim);
435  interpolated_x(s,x);
436 
437  // Outer unit normal (only in r and z direction!)
438  Vector<double> unit_normal(n_dim);
439  outer_unit_normal(s,unit_normal);
440 
441  // Dummy
442  unsigned ipt=0;
443 
444  // Traction vector
445  get_traction(time,ipt,x,unit_normal,traction);
446 
447  }
448 
449 
450 //=====================================================================
451 /// Return the residuals for the
452 /// AxisymmetricNavierStokesTractionElement equations
453 //=====================================================================
454 template<class ELEMENT>
457  Vector<double> &residuals)
458  {
459 
460  //Find out how many nodes there are
461  unsigned n_node = nnode();
462 
463  // Get continuous time from timestepper of first node
464  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
465 
466  #ifdef PARANOID
467  //Find out how many positional dofs there are
468  unsigned n_position_type = this->nnodal_position_type();
469  if(n_position_type != 1)
470  {
471  throw OomphLibError(
472  "AxisymmetricNavierStokes is not yet implemented for more than one position type",
473  OOMPH_CURRENT_FUNCTION,
474  OOMPH_EXCEPTION_LOCATION);
475  }
476 #endif
477 
478  //Find out the dimension of the node
479  unsigned n_dim = this->nodal_dimension();
480 
481  //Cache the nodal indices at which the velocity components are stored
482  unsigned u_nodal_index[n_dim+1];
483  for(unsigned i=0;i<n_dim+1;i++)
484  {
485  u_nodal_index[i] = this->U_index_axisymmetric_nst_traction[i];
486  }
487 
488  //Integer to hold the local equation number
489  int local_eqn=0;
490 
491  //Set up memory for the shape functions
492  //Note that in this case, the number of lagrangian coordinates is always
493  //equal to the dimension of the nodes
494  Shape psi(n_node);
495  DShape dpsids(n_node,n_dim-1);
496 
497  //Set the value of n_intpt
498  unsigned n_intpt = integral_pt()->nweight();
499 
500  //Loop over the integration points
501  for(unsigned ipt=0;ipt<n_intpt;ipt++)
502  {
503  //Get the integral weight
504  double w = integral_pt()->weight(ipt);
505 
506  //Only need to call the local derivatives
507  dshape_local_at_knot(ipt,psi,dpsids);
508 
509  //Calculate the Eulerian and Lagrangian coordinates
510  Vector<double> interpolated_x(n_dim,0.0);
511 
512  //Also calculate the surface Vectors (derivatives wrt local coordinates)
513  DenseMatrix<double> interpolated_A(n_dim-1,n_dim,0.0);
514 
515  //Calculate positions and derivatives
516  for(unsigned l=0;l<n_node;l++)
517  {
518  //Loop over directions
519  for(unsigned i=0;i<n_dim;i++)
520  {
521  //Calculate the Eulerian coords
522  const double x_local = nodal_position(l,i);
523  interpolated_x[i] += x_local*psi(l);
524 
525  //Loop over LOCAL derivative directions, to calculate the tangent(s)
526  for(unsigned j=0;j<n_dim-1;j++)
527  {
528  interpolated_A(j,i) += x_local*dpsids(l,j);
529  }
530  }
531  }
532 
533  //Now find the local metric tensor from the tangent Vectors
534  DenseMatrix<double> A(n_dim-1);
535  for(unsigned i=0;i<n_dim-1;i++)
536  {
537  for(unsigned j=0;j<n_dim-1;j++)
538  {
539  //Initialise surface metric tensor to zero
540  A(i,j) = 0.0;
541 
542  //Take the dot product
543  for(unsigned k=0;k<n_dim;k++)
544  {
545  A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
546  }
547  }
548  }
549 
550  //Get the outer unit normal
551  Vector<double> interpolated_normal(n_dim);
552  outer_unit_normal(ipt,interpolated_normal);
553 
554  //Find the determinant of the metric tensor
555  double Adet =0.0;
556  switch(n_dim)
557  {
558  case 2:
559  Adet = A(0,0);
560  break;
561  case 3:
562  Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
563  break;
564  default:
565  throw
567  "Wrong dimension in AxisymmetricNavierStokesTractionElement",
568  "AxisymmetricNavierStokesTractionElement::fill_in_contribution_to_residuals()",
569  OOMPH_EXCEPTION_LOCATION);
570  }
571 
572  //Premultiply the weights and the square-root of the determinant of
573  //the metric tensor
574  double W = w*sqrt(Adet);
575 
576  //Now calculate the load
577  Vector<double> traction(n_dim+1);
578  get_traction(time,
579  ipt,
580  interpolated_x,
581  interpolated_normal,
582  traction);
583 
584  //Loop over the test functions, nodes of the element
585  for(unsigned l=0;l<n_node;l++)
586  {
587  //Loop over the velocity components
588  for(unsigned i=0;i<n_dim+1;i++)
589  {
590  // Equation number
591  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
592  /*IF it's not a boundary condition*/
593  if(local_eqn >= 0)
594  {
595  //Add the loading terms to the residuals
596  residuals[local_eqn] -= traction[i]*psi(l)*interpolated_x[0]*W;
597  }
598  }
599  } //End of loop over shape functions
600  } //End of loop over integration points
601 
602  }
603 
604 
605 
606 
607 /////////////////////////////////////////////////////////////////////////
608 /////////////////////////////////////////////////////////////////////////
609 /////////////////////////////////////////////////////////////////////////
610 
611 
612 //=======================================================================
613 /// Namespace containing the default Strouhal number of axisymmetric
614 /// linearised FSI.
615 //=======================================================================
616 namespace LinearisedFSIAxisymmetricNStNoSlipBCHelper
617  {
618 
619  /// Default for fluid Strouhal number
620  extern double Default_strouhal_number;
621 
622  }
623 
624 
625 
626 //======================================================================
627 /// \short A class for elements that allow the imposition of the linearised
628 /// FSI no slip condition from an adjacent linearly elastic axisymmetric
629 /// solid. The element geometry is obtained from the FaceGeometry<ELEMENT>
630 /// policy class.
631 //======================================================================
632  template <class FLUID_BULK_ELEMENT, class SOLID_BULK_ELEMENT>
634  public virtual FaceGeometry<FLUID_BULK_ELEMENT>,
635  public virtual FaceElement,
636  public virtual ElementWithExternalElement
637  {
638 
639  public:
640 
641 
642  /// \short Constructor, takes the pointer to the "bulk" element and the
643  /// face index identifying the face to which the element is attached.
644  /// The optional identifier can be used
645  /// to distinguish the additional nodal values created by
646  /// this element from thos created by other FaceElements.
648  FiniteElement* const &bulk_el_pt,
649  const int& face_index,
650  const unsigned &id=0);
651 
652  /// Broken copy constructor
655  {
657  "LinearisedFSIAxisymmetricNStNoSlipBCElementElement");
658  }
659 
660  /// Broken assignment operator
661 //Commented out broken assignment operator because this can lead to a conflict warning
662 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
663 //realise that two separate implementations of the broken function are the same and so,
664 //quite rightly, it shouts.
665  /*void operator=(const LinearisedFSIAxisymmetricNStNoSlipBCElementElement&)
666  {
667  BrokenCopy::broken_assign(
668  "LinearisedFSIAxisymmetricNStNoSlipBCElementElement");
669  }*/
670 
671 
672  /// \short Access function for the pointer to the fluid Strouhal number
673  /// (if not set, St defaults to 1)
674  double* &st_pt() {return St_pt;}
675 
676  /// Access function for the fluid Strouhal number
677  double st() const
678  {
679  return *St_pt;
680  }
681 
682  /// Add the element's contribution to its residual vector
684  {
685  //Call the generic residuals function with flag set to 0
686  //using a dummy matrix argument
689  }
690 
691 
692  /// \short Add the element's contribution to its residual vector and its
693  /// Jacobian matrix
695  DenseMatrix<double> &jacobian)
696  {
697  //Call the generic routine with the flag set to 1
699  (residuals,jacobian,1);
700 
701  //Derivatives w.r.t. external data
703  }
704 
705  /// Output function
706  void output(std::ostream &outfile)
707  {
708  //Dummy
709  unsigned nplot=0;
710  output(outfile,nplot);
711  }
712 
713  /// Output function: Output at Gauss points; n_plot is ignored.
714  void output(std::ostream &outfile, const unsigned &n_plot)
715  {
716  outfile << "ZONE\n";
717 
718  //Get the value of Nintpt
719  const unsigned n_intpt = integral_pt()->nweight();
720 
721  //Set the Vector to hold local coordinates
722  Vector<double> s_int(Dim-1);
723 
724  //Loop over the integration points
725  for(unsigned ipt=0;ipt<n_intpt;ipt++)
726  {
727 
728  //Assign values of s
729  for(unsigned i=0;i<(Dim-1);i++)
730  {
731  s_int[i] = integral_pt()->knot(ipt,i);
732  }
733 
734  // Get boundary coordinate
735  Vector<double> zeta(1);
736  interpolated_zeta(s_int,zeta);
737 
738  // Get velocity from adjacent solid
739  SOLID_BULK_ELEMENT* ext_el_pt=dynamic_cast<SOLID_BULK_ELEMENT*>(
740  external_element_pt(0,ipt));
742  Vector<double> dudt(3);
743  ext_el_pt->interpolated_du_dt_axisymmetric_linear_elasticity(s_ext,dudt);
744 
745  // Output
746  outfile << ext_el_pt->interpolated_x(s_ext,0) << " "
747  << ext_el_pt->interpolated_x(s_ext,1) << " "
748  << dudt[0] << " "
749  << dudt[1] << " "
750  << dudt[2] << " "
751  << zeta[0] << std::endl;
752  }
753  }
754 
755 
756  /// C-style output function
757  void output(FILE* file_pt)
759 
760  /// C-style output function
761  void output(FILE* file_pt, const unsigned &n_plot)
763 
764 
765 
766  protected:
767 
768  /// \short Function to compute the shape and test functions and to return
769  /// the Jacobian of mapping between local and global (Eulerian)
770  /// coordinates
771  inline double shape_and_test(const Vector<double> &s, Shape &psi, Shape &test)
772  const
773  {
774  //Find number of nodes
775  unsigned n_node = nnode();
776 
777  //Get the shape functions
778  shape(s,psi);
779 
780  //Set the test functions to be the same as the shape functions
781  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
782 
783  //Return the value of the jacobian
784  return J_eulerian(s);
785  }
786 
787 
788  /// \short Function to compute the shape and test functions and to return
789  /// the Jacobian of mapping between local and global (Eulerian)
790  /// coordinates
791  inline double shape_and_test_at_knot(const unsigned &ipt,
792  Shape &psi, Shape &test)
793  const
794  {
795  //Find number of nodes
796  unsigned n_node = nnode();
797 
798  //Get the shape functions
799  shape_at_knot(ipt,psi);
800 
801  //Set the test functions to be the same as the shape functions
802  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
803 
804  //Return the value of the jacobian
805  return J_eulerian_at_knot(ipt);
806  }
807 
808 
809 
810  private:
811 
812 
813  /// \short Add the element's contribution to its residual vector.
814  /// flag=1(or 0): do (or don't) compute the contribution to the
815  /// Jacobian as well.
817  Vector<double> &residuals, DenseMatrix<double> &jacobian,
818  const unsigned& flag);
819 
820  ///The spatial dimension of the problem
821  unsigned Dim;
822 
823  ///The index at which the unknowns are stored at the nodes
825 
826  /// Lagrange Id
827  unsigned Id;
828 
829  /// Pointer to fluid Strouhal number
830  double* St_pt;
831 
832  };
833 
834 //////////////////////////////////////////////////////////////////////
835 //////////////////////////////////////////////////////////////////////
836 //////////////////////////////////////////////////////////////////////
837 
838 
839 
840 //===========================================================================
841 /// Constructor, takes the pointer to the "bulk" element, and the
842 /// face index that identifies the face of the bulk element to which
843 /// this face element is to be attached.
844 /// The optional identifier can be used
845 /// to distinguish the additional nodal values created by
846 /// this element from thos created by other FaceElements.
847 //===========================================================================
848  template <class FLUID_BULK_ELEMENT, class SOLID_BULK_ELEMENT>
850  SOLID_BULK_ELEMENT>::
852  FiniteElement* const &bulk_el_pt,
853  const int &face_index,
854  const unsigned &id) :
855  FaceGeometry<FLUID_BULK_ELEMENT>(), FaceElement()
856  {
857  // Set source element storage: one interaction with an external element
858  // that provides the velocity of the adjacent linear elasticity
859  // element
860  this->set_ninteraction(1);
861 
862  // Store the ID of the FaceElement -- this is used to distinguish
863  // it from any others
864  Id=id;
865 
866  // Initialise pointer to fluid Strouhal number. Defaults to 1
868 
869  // Let the bulk element build the FaceElement, i.e. setup the pointers
870  // to its nodes (by referring to the appropriate nodes in the bulk
871  // element), etc.
872  bulk_el_pt->build_face_element(face_index,this);
873 
874  // Extract the dimension of the problem from the dimension of
875  // the first node
876  Dim = this->node_pt(0)->ndim();
877 
878  //Read the index from the (cast) bulk element.
879  U_index_fsi_no_slip_axisym.resize(3);
880  for (unsigned i=0;i<3;i++)
881  {
883  dynamic_cast<FLUID_BULK_ELEMENT*>(bulk_el_pt)->u_index_axi_nst(i);
884  }
885 
886  // We need Dim+1 additional values for each FaceElement node
887  // to store the Lagrange multipliers.
888  Vector<unsigned> n_additional_values(nnode(), Dim+1);
889 
890  // Now add storage for Lagrange multipliers and set the map containing
891  // the position of the first entry of this face element's
892  // additional values.
893  add_additional_values(n_additional_values,id);
894 
895  }
896 
897 
898 //===========================================================================
899 /// \short Helper function to compute the element's residual vector and
900 /// the Jacobian matrix.
901 //===========================================================================
902 template <class FLUID_BULK_ELEMENT, class SOLID_BULK_ELEMENT>
904  SOLID_BULK_ELEMENT>:: fill_in_generic_residual_contribution_fsi_no_slip_axisym(
905  Vector<double> &residuals, DenseMatrix<double> &jacobian,
906  const unsigned& flag)
907  {
908  //Find out how many nodes there are
909  const unsigned n_node = nnode();
910 
911  //Set up memory for the shape and test functions
912  Shape psif(n_node), testf(n_node);
913 
914  //Set the value of Nintpt
915  const unsigned n_intpt = integral_pt()->nweight();
916 
917  //Set the Vector to hold local coordinates
918  Vector<double> s(Dim-1);
919 
920  // Cache the Strouhal number
921  const double local_st=st();
922 
923  //Integers to hold the local equation and unknown numbers
924  int local_eqn=0;
925 
926  //Loop over the integration points
927  //--------------------------------
928  for(unsigned ipt=0;ipt<n_intpt;ipt++)
929  {
930 
931  //Assign values of s
932  for(unsigned i=0;i<(Dim-1);i++) {s[i] = integral_pt()->knot(ipt,i);}
933 
934  //Get the integral weight
935  double w = integral_pt()->weight(ipt);
936 
937  //Find the shape and test functions and return the Jacobian
938  //of the mapping
939  double J = shape_and_test(s,psif,testf);
940 
941  //Premultiply the weights and the Jacobian
942  double W = w*J;
943 
944  //Calculate the Lagrange multiplier and the fluid veloc
945  Vector<double> lambda(Dim+1,0.0);
946  Vector<double> fluid_veloc(Dim+1,0.0);
947 
948  // Loop over nodes
949  for(unsigned j=0;j<n_node;j++)
950  {
951  Node* nod_pt=node_pt(j);
952 
953  // Cast to a boundary node
954  BoundaryNodeBase *bnod_pt=dynamic_cast<BoundaryNodeBase*>(node_pt(j));
955 
956  // Get the index of the first nodal value associated with
957  // this FaceElement
958  unsigned first_index=
960 
961  //Assemble
962  for(unsigned i=0;i<Dim+1;i++)
963  {
964  lambda[i]+=nod_pt->value(first_index+i)*psif(j);
965  fluid_veloc[i]+=nod_pt->value(U_index_fsi_no_slip_axisym[i])*psif(j);
966  }
967  }
968 
969  // Get velocity from adjacent solid
970  SOLID_BULK_ELEMENT* ext_el_pt=dynamic_cast<SOLID_BULK_ELEMENT*>(
971  external_element_pt(0,ipt));
972  Vector<double> s_ext(external_element_local_coord(0,ipt));
973  Vector<double> dudt(3);
974  ext_el_pt->interpolated_du_dt_axisymmetric_linear_elasticity(s_ext,dudt);
975 
976 
977  //Now add to the appropriate equations
978 
979  //Loop over the test functions
980  for(unsigned l=0;l<n_node;l++)
981  {
982 
983  // Loop over directions
984  for (unsigned i=0;i<Dim+1;i++)
985  {
986 
987  // Add contribution to bulk Navier Stokes equations where
988  //-------------------------------------------------------
989  // the Lagrange multiplier acts as a traction
990  //-------------------------------------------
991  local_eqn = nodal_local_eqn(l,U_index_fsi_no_slip_axisym[i]);
992 
993  /*IF it's not a boundary condition*/
994  if(local_eqn >= 0)
995  {
996  //Add the Lagrange multiplier "traction" to the bulk
997  residuals[local_eqn] -= lambda[i]*testf[l]*W;
998 
999 
1000  //Jacobian entries
1001  if(flag)
1002  {
1003  //Loop over the lagrange multiplier unknowns
1004  for(unsigned l2=0;l2<n_node;l2++)
1005  {
1006  // Cast to a boundary node
1007  BoundaryNodeBase *bnod_pt =
1008  dynamic_cast<BoundaryNodeBase*>(node_pt(l2));
1009 
1010  // Local unknown
1011  int local_unknown=nodal_local_eqn
1012  (l2,bnod_pt->
1013  index_of_first_value_assigned_by_face_element(Id)+i);
1014 
1015  // If it's not pinned
1016  if (local_unknown>=0)
1017  {
1018  jacobian(local_eqn,local_unknown) -=
1019  psif[l2]*testf[l]*W;
1020  }
1021  }
1022  }
1023  }
1024 
1025  // Now do the Lagrange multiplier equations
1026  //-----------------------------------------
1027  // Cast to a boundary node
1028  BoundaryNodeBase *bnod_pt =
1029  dynamic_cast<BoundaryNodeBase*>(node_pt(l));
1030 
1031  // Local eqn number:
1032  int local_eqn=nodal_local_eqn
1034 
1035  // If it's not pinned
1036  if (local_eqn>=0)
1037  {
1038  residuals[local_eqn]+=(local_st*dudt[i]-fluid_veloc[i])*testf(l)*W;
1039 
1040  //Jacobian entries
1041  if(flag)
1042  {
1043  // Loop over the velocity unknowns [derivs w.r.t. to
1044  // wall velocity taken care of by fd-ing
1045  for(unsigned l2=0;l2<n_node;l2++)
1046  {
1047  int local_unknown =
1048  nodal_local_eqn(l2,U_index_fsi_no_slip_axisym[i]);
1049 
1050  /*IF it's not a boundary condition*/
1051  if(local_unknown >= 0)
1052  {
1053  jacobian(local_eqn,local_unknown) -=
1054  psif[l2]*testf[l]*W;
1055  }
1056  }
1057  }
1058 
1059  }
1060 
1061  }
1062  }
1063  }
1064  }
1065 
1066 }
1067 
1068 
1069 
1070 #endif
void Zero_traction_fct(const double &time, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Default load function (zero traction)
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void fill_in_contribution_to_residuals_axisymmetric_nst_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
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 st() const
Access function for the fluid Strouhal number.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation...
Definition: nodes.cc:2314
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
void scalar_value_paraview(std::ofstream &file_out, const unsigned &k, const unsigned &nplot) const
Write values of the k-th scalar field at the plot points. Needs to be implemented for each new specif...
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5121
AxisymmetricNavierStokesTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
void fill_in_generic_residual_contribution_fsi_no_slip_axisym(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add the element's contribution to its residual vector. flag=1(or 0): do (or don't) compute the contri...
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class...
cstr elem_len * i
Definition: cfortran.h:607
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
void output(std::ostream &outfile)
Output function.
LinearisedFSIAxisymmetricNStNoSlipBCElementElement(const LinearisedFSIAxisymmetricNStNoSlipBCElementElement &dummy)
Broken copy constructor.
void fill_in_jacobian_from_external_interaction_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the jacobian from all external interaction degrees of freedom (geometr...
void output(FILE *file_pt)
C_style output 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
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation.
LinearisedFSIAxisymmetricNStNoSlipBCElementElement(FiniteElement *const &bulk_el_pt, const int &face_index, const unsigned &id=0)
Constructor, takes the pointer to the "bulk" element and the face index identifying the face to which...
A general Finite Element class.
Definition: elements.h:1271
unsigned index_of_first_value_assigned_by_face_element(const unsigned &face_id=0) const
Return the index of the first value associated with the i-th face element value. If no argument is sp...
Definition: nodes.h:1924
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void set_ninteraction(const unsigned &n_interaction)
Set the number of interactions in the element This function is usually called in the specific element...
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2358
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
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 unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2700
void(*&)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction) traction_fct_pt()
Reference to the traction function pointer.
void add_additional_values(const Vector< unsigned > &nadditional_values, const unsigned &id)
Helper function adding additional values for the unknowns associated with the FaceElement. This function also sets the map containing the position of the first entry of this face element's additional values.The inputs are the number of additional values and the face element's ID. Note the same number of additonal values are allocated at ALL nodes.
Definition: elements.h:4181
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:66
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and its Jacobian matrix.
void traction(const double &time, const Vector< double > &s, Vector< double > &traction)
Compute traction vector at specified local coordinate Should only be used for post-processing; ignore...
Vector< unsigned > U_index_axisymmetric_nst_traction
Index at which the i-th velocity component is stored.
A class that contains the information required by Nodes that are located on Mesh boundaries. A BoundaryNode of a particular type is obtained by combining a given Node with this class. By differentiating between Nodes and BoundaryNodes we avoid a lot of un-necessary storage in the bulk Nodes.
Definition: nodes.h:1854
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
FaceElement()
Constructor: Initialise all appropriate member data.
Definition: elements.h:4201
double Default_strouhal_number
Default for fluid Strouhal number.
double shape_and_test_at_knot(const unsigned &ipt, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return.
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:540
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.
double shape_and_test(const Vector< double > &s, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3158
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:5033
void output(std::ostream &outfile)
Output with default number of plot points.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary, when viewed as part of a compound geometric object is specified using the boundary coordinate defined by the mesh. Note: Boundary coordinates will have been set up when creating the underlying mesh, and their values will have been stored at the nodes.
Definition: elements.h:4251
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
Vector< unsigned > U_index_fsi_no_slip_axisym
The index at which the unknowns are stored at the nodes.
A class for elements that allow the imposition of the linearised FSI no slip condition from an adjace...
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
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point...
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
Definition: elements.cc:4480
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 void get_traction(const double &time, const unsigned &intpt, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction) const
Get the traction vector: Pass number of integration point (dummy), Eulerian coordinate and normal vec...
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
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Pointer to an imposed traction function. Arguments: Eulerian coordinate; outer unit normal; applied t...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: Output at Gauss points; n_plot is ignored.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
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...
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:228
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types)...