axisym_linear_elasticity_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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header file for elements that are used to apply surface loads to
31 //the equations of axisymmetric linear elasticity
32 
33 #ifndef OOMPH_AXISYMMETRIC_LINEAR_ELASTICITY_TRACTION_ELEMENTS_HEADER
34 #define OOMPH_AXISYMMETRIC_LINEAR_ELASTICITY_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/Qelements.h"
44 #include "src/generic/Qelements.h"
46 
47 
48 namespace oomph
49 {
50 
51 
52 //=======================================================================
53 /// Namespace containing the zero traction function for
54 /// axisymmetric linear elasticity traction elements
55 //=======================================================================
56 namespace AxisymmetricLinearElasticityTractionElementHelper
57  {
58 
59  //=======================================================================
60  /// Default load function (zero traction)
61  //=======================================================================
62  void Zero_traction_fct(const double& time,
63  const Vector<double> &x,
64  const Vector<double>& N,
65  Vector<double>& load)
66  {
67  unsigned n_dim=load.size();
68  for (unsigned i=0;i<n_dim;i++) {load[i]=0.0;}
69  }
70 
71  }
72 
73 
74 //======================================================================
75 /// A class for elements that allow the imposition of an applied traction
76 /// in the equations of axisymmetric linear elasticity.
77 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
78 /// class and and thus, we can be generic enough without the need to have
79 /// a separate equations class.
80 //======================================================================
81 template <class ELEMENT>
83 public virtual FaceGeometry<ELEMENT>,
84  public virtual FaceElement
85  {
86  protected:
87 
88  /// Index at which the i-th displacement component is stored
91 
92  /// \short Pointer to an imposed traction function. Arguments:
93  /// Eulerian coordinate; outer unit normal;
94  /// applied traction. (Not all of the input arguments will be
95  /// required for all specific load functions but the list should
96  /// cover all cases)
97  void (*Traction_fct_pt)(const double& time,
98  const Vector<double> &x,
99  const Vector<double> &n,
100  Vector<double> &result);
101 
102 
103  /// \short Get the traction vector: Pass number of integration point (dummy),
104  /// Eulerian coordinate and normal vector and return the load vector
105  /// (not all of the input arguments will be
106  /// required for all specific load functions but the list should
107  /// cover all cases). This function is virtual so it can be
108  /// overloaded for FSI.
109  virtual void get_traction(const double& time,
110  const unsigned& intpt,
111  const Vector<double>& x,
112  const Vector<double>& n,
114  {
115  Traction_fct_pt(time,x,n,traction);
116  }
117 
118 
119  /// \short Helper function that actually calculates the residuals
120  // This small level of indirection is required to avoid calling
121  // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
122  // which causes all kinds of pain if overloading later on
124  Vector<double> &residuals);
125 
126 
127  public:
128 
129  /// \short Constructor, which takes a "bulk" element and the
130  /// value of the index and its limit
132  (FiniteElement* const &element_pt, const int &face_index) :
134  {
135  //Attach the geometrical information to the element. N.B. This function
136  //also assigns nbulk_value from the required_nvalue of the bulk element
137  element_pt->build_face_element(face_index,this);
138 
139  //Find the dimension of the problem
140  unsigned n_dim = element_pt->nodal_dimension();
141 
142  //Find the index at which the displacement unknowns are stored
143  ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
145  for(unsigned i=0;i<n_dim+1;i++)
146  {
147  this->
149  cast_element_pt->
150  u_index_axisymmetric_linear_elasticity(i);
151  }
152 
153  // Zero traction
157  }
158 
159 
160  /// Reference to the traction function pointer
161  void (* &traction_fct_pt())(const double& time,
162  const Vector<double>& x,
163  const Vector<double>& n,
165  {return Traction_fct_pt;}
166 
167 
168  /// Return the residuals
170  {
172  (residuals);
173  }
174 
175 
176 
177  /// Fill in contribution from Jacobian
179  DenseMatrix<double> &jacobian)
180  {
181  //Call the residuals
183  (residuals);
184  }
185 
186  /// Specify the value of nodal zeta from the face geometry
187  /// \short The "global" intrinsic coordinate of the element when
188  /// viewed as part of a geometric object should be given by
189  /// the FaceElement representation, by default (needed to break
190  /// indeterminacy if bulk element is SolidElement)
191  double zeta_nodal(const unsigned &n, const unsigned &k,
192  const unsigned &i) const
193  {return FaceElement::zeta_nodal(n,k,i);}
194 
195  /// \short Output function
196  void output(std::ostream &outfile)
197  {
198  unsigned nplot=5;
199  output(outfile,nplot);
200  }
201 
202  /// \short Output function
203  void output(std::ostream &outfile, const unsigned &n_plot)
204  {
205  //Number of dimensions
206  unsigned n_dim = this->nodal_dimension();
207 
208  //Find out how many nodes there are
209  const unsigned n_node = nnode();
210 
211  // Get continuous time from timestepper of first node
212  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
213 
214  //Set up memory for the shape functions
215  Shape psi(n_node);
216 
217  // Local and global coordinates
218  Vector<double> s(n_dim-1);
220 
221  // Tecplot header info
222  outfile << this->tecplot_zone_string(n_plot);
223 
224  // Loop over plot points
225  unsigned num_plot_points=this->nplot_points(n_plot);
226  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
227  {
228  // Get local coordinates of plot point
229  this->get_s_plot(iplot,n_plot,s);
230 
231  // Outer unit normal
232  Vector<double> unit_normal(n_dim);
233  outer_unit_normal(s,unit_normal);
234 
235  //Find the shape functions
236  shape(s,psi);
237 
238  //Initialise to zero
239  for(unsigned i=0;i<n_dim;i++)
240  {
241  interpolated_x[i] = 0.0;
242  }
243 
244  //Calculate stuff
245  for(unsigned l=0;l<n_node;l++)
246  {
247  //Loop over directions
248  for(unsigned i=0;i<n_dim;i++)
249  {
250  interpolated_x[i] += nodal_position(l,i)*psi[l];
251  }
252  }
253 
254  //Get the imposed flux
256 
257  //Dummy integration point
258  unsigned ipt=0;
259 
260  get_traction(time,ipt,interpolated_x,unit_normal,traction);
261 
262  //Output the x,y,..
263  for(unsigned i=0;i<n_dim;i++)
264  {
265  outfile << interpolated_x[i] << " ";
266  }
267 
268  //Output the traction components
269  for(unsigned i=0;i<n_dim+1;i++)
270  {
271  outfile << traction[i] << " ";
272  }
273 
274  // Output normal
275  for(unsigned i=0;i<n_dim;i++)
276  {
277  outfile << unit_normal[i] << " ";
278  }
279  outfile << std::endl;
280 
281  }
282  }
283 
284  /// \short C_style output function
285  void output(FILE* file_pt)
286  {FiniteElement::output(file_pt);}
287 
288  /// \short C-style output function
289  void output(FILE* file_pt, const unsigned &n_plot)
290  {FiniteElement::output(file_pt,n_plot);}
291 
292 
293  /// \short Compute traction vector at specified local coordinate
294  /// Should only be used for post-processing; ignores dependence
295  /// on integration point!
296  void traction(const double &time,
297  const Vector<double>& s,
299 
300  };
301 
302 ///////////////////////////////////////////////////////////////////////
303 ///////////////////////////////////////////////////////////////////////
304 ///////////////////////////////////////////////////////////////////////
305 
306 //=====================================================================
307 /// Compute traction vector at specified local coordinate
308 /// Should only be used for post-processing; ignores dependence
309 /// on integration point!
310 //=====================================================================
311 template<class ELEMENT>
313  traction(const double& time, const Vector<double>& s, Vector<double>& traction)
314  {
315  unsigned n_dim = this->nodal_dimension();
316 
317  // Position vector
318  Vector<double> x(n_dim);
319  interpolated_x(s,x);
320 
321  // Outer unit normal (only in r and z direction!)
322  Vector<double> unit_normal(n_dim);
323  outer_unit_normal(s,unit_normal);
324 
325  // Dummy
326  unsigned ipt=0;
327 
328  // Traction vector
329  get_traction(time,ipt,x,unit_normal,traction);
330 
331  }
332 
333 
334 //=====================================================================
335 /// Return the residuals for the
336 /// AxisymmetricLinearElasticityTractionElement equations
337 //=====================================================================
338 template<class ELEMENT>
341  Vector<double> &residuals)
342  {
343 
344  //Find out how many nodes there are
345  unsigned n_node = nnode();
346 
347  // Get continuous time from timestepper of first node
348  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
349 
350  #ifdef PARANOID
351  //Find out how many positional dofs there are
352  unsigned n_position_type = this->nnodal_position_type();
353  if(n_position_type != 1)
354  {
355  throw OomphLibError(
356  "AxisymmetricLinearElasticity is not yet implemented for more than one position type",
357  OOMPH_CURRENT_FUNCTION,
358  OOMPH_EXCEPTION_LOCATION);
359  }
360 #endif
361 
362  //Find out the dimension of the node
363  unsigned n_dim = this->nodal_dimension();
364 
365  //Cache the nodal indices at which the displacement components are stored
366  unsigned u_nodal_index[n_dim+1];
367  for(unsigned i=0;i<n_dim+1;i++)
368  {
369  u_nodal_index[i]=this->U_index_axisymmetric_linear_elasticity_traction[i];
370  }
371 
372  //Integer to hold the local equation number
373  int local_eqn=0;
374 
375  //Set up memory for the shape functions
376  //Note that in this case, the number of lagrangian coordinates is always
377  //equal to the dimension of the nodes
378  Shape psi(n_node);
379  DShape dpsids(n_node,n_dim-1);
380 
381  //Set the value of n_intpt
382  unsigned n_intpt = integral_pt()->nweight();
383 
384  //Loop over the integration points
385  for(unsigned ipt=0;ipt<n_intpt;ipt++)
386  {
387  //Get the integral weight
388  double w = integral_pt()->weight(ipt);
389 
390  //Only need to call the local derivatives
391  dshape_local_at_knot(ipt,psi,dpsids);
392 
393  //Calculate the Eulerian and Lagrangian coordinates
394  Vector<double> interpolated_x(n_dim,0.0);
395 
396  //Also calculate the surface Vectors (derivatives wrt local coordinates)
397  DenseMatrix<double> interpolated_A(n_dim-1,n_dim,0.0);
398 
399  //Calculate displacements and derivatives
400  for(unsigned l=0;l<n_node;l++)
401  {
402  //Loop over directions
403  for(unsigned i=0;i<n_dim;i++)
404  {
405  //Calculate the Eulerian coords
406  const double x_local = nodal_position(l,i);
407  interpolated_x[i] += x_local*psi(l);
408 
409  //Loop over LOCAL derivative directions, to calculate the tangent(s)
410  for(unsigned j=0;j<n_dim-1;j++)
411  {
412  interpolated_A(j,i) += x_local*dpsids(l,j);
413  }
414  }
415  }
416 
417  //Now find the local metric tensor from the tangent Vectors
418  DenseMatrix<double> A(n_dim-1);
419  for(unsigned i=0;i<n_dim-1;i++)
420  {
421  for(unsigned j=0;j<n_dim-1;j++)
422  {
423  //Initialise surface metric tensor to zero
424  A(i,j) = 0.0;
425 
426  //Take the dot product
427  for(unsigned k=0;k<n_dim;k++)
428  {
429  A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
430  }
431  }
432  }
433 
434  //Get the outer unit normal
435  Vector<double> interpolated_normal(n_dim);
436  outer_unit_normal(ipt,interpolated_normal);
437 
438  //Find the determinant of the metric tensor
439  double Adet =0.0;
440  switch(n_dim)
441  {
442  case 2:
443  Adet = A(0,0);
444  break;
445  case 3:
446  Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
447  break;
448  default:
449  throw
451  "Wrong dimension in AxisymmetricLinearElasticityTractionElement",
452  "AxisymmetricLinearElasticityTractionElement::fill_in_contribution_to_residuals()",
453  OOMPH_EXCEPTION_LOCATION);
454  }
455 
456  //Premultiply the weights and the square-root of the determinant of
457  //the metric tensor
458  double W = w*sqrt(Adet);
459 
460  //Now calculate the load
461  Vector<double> traction(n_dim+1);
462  get_traction(time,
463  ipt,
464  interpolated_x,
465  interpolated_normal,
466  traction);
467 
468  //Loop over the test functions, nodes of the element
469  for(unsigned l=0;l<n_node;l++)
470  {
471  //Loop over the displacement components
472  for(unsigned i=0;i<n_dim+1;i++)
473  {
474  // Equation number
475  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
476  /*IF it's not a boundary condition*/
477  if(local_eqn >= 0)
478  {
479  //Add the loading terms to the residuals
480  residuals[local_eqn] -= traction[i]*psi(l)*interpolated_x[0]*W;
481  }
482  }
483  } //End of loop over shape functions
484  } //End of loop over integration points
485 
486  }
487 
488 
489 ///////////////////////////////////////////////////////////////////////
490 ///////////////////////////////////////////////////////////////////////
491 ///////////////////////////////////////////////////////////////////////
492 
493 
494 
495 //======================================================================
496 /// A class for elements that allow the imposition of an applied traction
497 /// in the equations of axisymmetric linear elasticity from an adjacent
498 /// axisymmetric Navier Stokes element in a linearised FSI problem.
499 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
500 /// class and and thus, we can be generic enough without the need to have
501 /// a separate equations class.
502 //======================================================================
503  template <class ELASTICITY_BULK_ELEMENT, class NAVIER_STOKES_BULK_ELEMENT>
505  public virtual FaceGeometry<ELASTICITY_BULK_ELEMENT>,
506  public virtual FaceElement,
507  public virtual ElementWithExternalElement
508  {
509 
510  protected:
511 
512  /// \short Pointer to the ratio, \f$ Q \f$ , of the stress used to
513  /// non-dimensionalise the fluid stresses to the stress used to
514  /// non-dimensionalise the solid stresses.
515  double *Q_pt;
516 
517  /// \short Static default value for the ratio of stress scales
518  /// used in the fluid and solid equations (default is 1.0)
519  static double Default_Q_Value;
520 
521  /// Index at which the i-th displacement component is stored
523 
524  /// \short Helper function that actually calculates the residuals
525  // This small level of indirection is required to avoid calling
526  // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
527  // which causes all kinds of pain if overloading later on
529  Vector<double> &residuals);
530 
531  public:
532 
533  /// \short Constructor, which takes a "bulk" element and the
534  /// value of the index and its limit
536  FiniteElement* const &element_pt,
537  const int &face_index) :
538  FaceGeometry<ELASTICITY_BULK_ELEMENT>(), FaceElement(),
540  {
541  // Set source element storage: one interaction with an external
542  // element -- the Navier Stokes bulk element that provides the traction
543  this->set_ninteraction(1);
544 
545  //Attach the geometrical information to the element. N.B. This function
546  //also assigns nbulk_value from the required_nvalue of the bulk element
547  element_pt->build_face_element(face_index,this);
548 
549  //Find the dimension of the problem
550  unsigned n_dim = element_pt->nodal_dimension();
551 
552  //Find the index at which the displacement unknowns are stored
553  ELASTICITY_BULK_ELEMENT* cast_element_pt =
554  dynamic_cast<ELASTICITY_BULK_ELEMENT*>(element_pt);
555  this->U_index_axisym_fsi_traction.resize(n_dim+1);
556  for(unsigned i=0;i<n_dim+1;i++)
557  {
559  cast_element_pt->u_index_axisymmetric_linear_elasticity(i);
560  }
561  }
562 
563 
564  /// Return the residuals
566  {
568  }
569 
570 
571 
572  /// Fill in contribution from Jacobian
574  DenseMatrix<double> &jacobian)
575  {
576  //Call the residuals
578 
579  //Derivatives w.r.t. external data
581  }
582 
583  /// \short Return the ratio of the stress scales used to non-dimensionalise
584  /// the fluid and elasticity equations. E.g.
585  /// \f$ Q = (\omega a)^2 \rho/E \f$, i.e. the
586  /// ratio between the inertial fluid stress and the solid's elastic
587  /// modulus E.
588  const double &q() const {return *Q_pt;}
589 
590  /// \short Return a pointer the ratio of stress scales used to
591  /// non-dimensionalise the fluid and solid equations.
592  double* &q_pt() {return Q_pt;}
593 
594 
595  /// \short Output function
596  void output(std::ostream &outfile)
597  {
598  /// Dummy
599  unsigned nplot=0;
600  output(outfile,nplot);
601  }
602 
603  /// \short Output function: Plot traction etc at Gauss points
604  /// nplot is ignored.
605  void output(std::ostream &outfile, const unsigned &n_plot)
606  {
607  // Dimension
608  unsigned n_dim = this->nodal_dimension();
609 
610  // Get FSI parameter
611  const double q_value=q();
612 
613  // Storage for traction (includes swirl!)
614  Vector<double> traction(n_dim+1);
615 
616  outfile << "ZONE\n";
617 
618  //Set the value of n_intpt
619  unsigned n_intpt = integral_pt()->nweight();
620 
621  //Loop over the integration points
622  for(unsigned ipt=0;ipt<n_intpt;ipt++)
623  {
624  Vector<double> s_int(n_dim-1);
625  s_int[0]=integral_pt()->knot(ipt,0);
626 
627  //Get the outer unit normal
628  Vector<double> interpolated_normal(n_dim);
629  outer_unit_normal(ipt,interpolated_normal);
630 
631  // Boundary coordinate
632  Vector<double> zeta(1);
633  interpolated_zeta(s_int,zeta);
634 
635  // Get bulk element for traction
636  NAVIER_STOKES_BULK_ELEMENT* ext_el_pt=
637  dynamic_cast<NAVIER_STOKES_BULK_ELEMENT*>
638  (this->external_element_pt(0,ipt));
639  Vector<double> s_ext(this->external_element_local_coord(0,ipt));
640 
641  // Get traction from bulk element (on fluid scale)
642  ext_el_pt->traction(s_ext,
643  interpolated_normal,
644  traction);
645 
646  outfile << ext_el_pt->interpolated_x(s_ext,0) << " "
647  << ext_el_pt->interpolated_x(s_ext,1) << " "
648  << q_value*traction[0] << " "
649  << q_value*traction[1] << " "
650  << q_value*traction[2] << " "
651  << interpolated_normal[0] << " "
652  << interpolated_normal[1] << " "
653  << zeta[0] << std::endl;
654  }
655  }
656 
657  /// \short C_style output function
658  void output(FILE* file_pt)
660 
661  /// \short C-style output function
662  void output(FILE* file_pt, const unsigned &n_plot)
664 
665  };
666 
667 ///////////////////////////////////////////////////////////////////////
668 ///////////////////////////////////////////////////////////////////////
669 ///////////////////////////////////////////////////////////////////////
670 
671 
672 
673 //=================================================================
674 /// Static default value for the ratio of stress scales
675 /// used in the fluid and solid equations (default is 1.0)
676 //=================================================================
677 template <class ELASTICITY_BULK_ELEMENT, class NAVIER_STOKES_BULK_ELEMENT>
678 double FSIAxisymmetricLinearElasticityTractionElement<
679  ELASTICITY_BULK_ELEMENT, NAVIER_STOKES_BULK_ELEMENT>::Default_Q_Value=1.0;
680 
681 
682 //=====================================================================
683 /// Return the residuals
684 //=====================================================================
685 template <class ELASTICITY_BULK_ELEMENT, class NAVIER_STOKES_BULK_ELEMENT>
686  void FSIAxisymmetricLinearElasticityTractionElement<
687  ELASTICITY_BULK_ELEMENT, NAVIER_STOKES_BULK_ELEMENT>::
688  fill_in_contribution_to_residuals_axisym_fsi_traction(
689  Vector<double> &residuals)
690  {
691 
692  //Find out how many nodes there are
693  unsigned n_node = nnode();
694 
695 #ifdef PARANOID
696  //Find out how many positional dofs there are
697  unsigned n_position_type = this->nnodal_position_type();
698  if(n_position_type != 1)
699  {
700  throw OomphLibError(
701  "LinearElasticity is not yet implemented for more than one position type.",
702  OOMPH_CURRENT_FUNCTION,
703  OOMPH_EXCEPTION_LOCATION);
704  }
705 #endif
706 
707  //Find out the dimension of the node
708  unsigned n_dim = this->nodal_dimension();
709 
710  //Cache the nodal indices at which the displacement components are stored
711  unsigned u_nodal_index[n_dim+1];
712  for(unsigned i=0;i<n_dim+1;i++)
713  {
714  u_nodal_index[i] = this->U_index_axisym_fsi_traction[i];
715  }
716 
717  //Integer to hold the local equation number
718  int local_eqn=0;
719 
720  //Set up memory for the shape functions
721  Shape psi(n_node);
722  DShape dpsids(n_node,n_dim-1);
723 
724  // Get FSI parameter
725  const double q_value=q();
726 
727  // Storage for traction
728  Vector<double> traction(3);
729 
730  //Set the value of n_intpt
731  unsigned n_intpt = integral_pt()->nweight();
732 
733  //Loop over the integration points
734  for(unsigned ipt=0;ipt<n_intpt;ipt++)
735  {
736  //Get the integral weight
737  double w = integral_pt()->weight(ipt);
738 
739  //Only need to call the local derivatives
740  dshape_local_at_knot(ipt,psi,dpsids);
741 
742  //Calculate the coordinates
743  Vector<double> interpolated_x(n_dim,0.0);
744 
745  //Also calculate the surface tangent vectors
746  DenseMatrix<double> interpolated_A(n_dim-1,n_dim,0.0);
747 
748  //Calculate displacements and derivatives
749  for(unsigned l=0;l<n_node;l++)
750  {
751  //Loop over directions
752  for(unsigned i=0;i<n_dim;i++)
753  {
754  //Calculate the Eulerian coords
755  const double x_local = nodal_position(l,i);
756  interpolated_x[i] += x_local*psi(l);
757 
758  //Loop over LOCAL derivative directions, to calculate the tangent(s)
759  for(unsigned j=0;j<n_dim-1;j++)
760  {
761  interpolated_A(j,i) += x_local*dpsids(l,j);
762  }
763  }
764  }
765 
766  //Now find the local metric tensor from the tangent Vectors
767  DenseMatrix<double> A(n_dim-1);
768  for(unsigned i=0;i<n_dim-1;i++)
769  {
770  for(unsigned j=0;j<n_dim-1;j++)
771  {
772  //Initialise surface metric tensor to zero
773  A(i,j) = 0.0;
774 
775  //Take the dot product
776  for(unsigned k=0;k<n_dim;k++)
777  {
778  A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
779  }
780  }
781  }
782 
783  //Get the outer unit normal
784  Vector<double> interpolated_normal(n_dim);
785  outer_unit_normal(ipt,interpolated_normal);
786 
787  //Find the determinant of the metric tensor
788  double Adet =0.0;
789  switch(n_dim)
790  {
791  case 2:
792  Adet = A(0,0);
793  break;
794  case 3:
795  Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
796  break;
797  default:
798  throw
800  "Wrong dimension in TimeHarmonicLinElastLoadedByPressureElement",
801  "TimeHarmonicLinElastLoadedByPressureElement::fill_in_contribution_to_residuals()",
802  OOMPH_EXCEPTION_LOCATION);
803  }
804 
805  //Premultiply the weights and the square-root of the determinant of
806  //the metric tensor
807  double W = w*sqrt(Adet);
808 
809  // Get bulk element for traction
810  NAVIER_STOKES_BULK_ELEMENT* ext_el_pt=
811  dynamic_cast<NAVIER_STOKES_BULK_ELEMENT*>
812  (this->external_element_pt(0,ipt));
813  Vector<double> s_ext(this->external_element_local_coord(0,ipt));
814 
815  // Get traction from bulk element
816  ext_el_pt->traction(s_ext,
817  interpolated_normal,
818  traction);
819 
820  //Loop over the test functions, nodes of the element
821  for(unsigned l=0;l<n_node;l++)
822  {
823  //Loop over the displacement components
824  for(unsigned i=0;i<n_dim+1;i++)
825  {
826  // Equation number
827  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
828  /*IF it's not a boundary condition*/
829  if(local_eqn >= 0)
830  {
831  //Add the loading terms (multiplied by fsi scaling factor)
832  //to the residuals
833  residuals[local_eqn] -=
834  q_value*traction[i]*psi(l)*interpolated_x[0]*W;
835  }
836  }
837  } //End of loop over shape functions
838 
839  } //End of loop over integration points
840  }
841 
842 
843 
844 ///////////////////////////////////////////////////////////////////////
845 ///////////////////////////////////////////////////////////////////////
846 ///////////////////////////////////////////////////////////////////////
847 
848 
849 
850 
851 }
852 
853 #endif
Vector< unsigned > U_index_axisymmetric_linear_elasticity_traction
Index at which the i-th displacement component is stored.
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.
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
FSIAxisymmetricLinearElasticityTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
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...
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
virtual void get_traction(const double &time, const unsigned &intpt, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Get the traction vector: Pass number of integration point (dummy), Eulerian coordinate and normal vec...
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 fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
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 fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
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
A general Finite Element class.
Definition: elements.h:1271
void set_ninteraction(const unsigned &n_interaction)
Set the number of interactions in the element This function is usually called in the specific element...
void fill_in_contribution_to_residuals_axisymmetric_linear_elasticity_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
Vector< unsigned > U_index_axisym_fsi_traction
Index at which the i-th displacement component is stored.
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2358
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
void(*&)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction) traction_fct_pt()
Reference to the traction function pointer.
double *& q_pt()
Return a pointer the ratio of stress scales used to non-dimensionalise the fluid and solid equations...
const double & q() const
Return the ratio of the stress scales used to non-dimensionalise the fluid and elasticity equations...
FaceElement()
Constructor: Initialise all appropriate member data.
Definition: elements.h:4201
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 ...
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:540
AxisymmetricLinearElasticityTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
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.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
double * Q_pt
Pointer to the ratio, , of the stress used to non-dimensionalise the fluid stresses to the stress us...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
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
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...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
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
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...
void Zero_traction_fct(const double &time, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Default load function (zero traction)
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: Plot traction etc at Gauss points nplot is ignored.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
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 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
void fill_in_contribution_to_residuals_axisym_fsi_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
static double Default_Q_Value
Static default value for the ratio of stress scales used in the fluid and solid equations (default is...
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...