fourier_decomposed_helmholtz_time_harmonic_linear_elasticity_interaction.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 #ifndef OOMPH_FD_HH_TIME_HARMONIC_LIN_ELAST_HEADER
31 #define OOMPH_FD_HH_TIME_HARMONIC_LIN_ELAST_HEADER
32 
33 
34 
35 //Oomph-lib includes
36 #include "generic.h"
39 
40 namespace oomph
41 {
42 
43 
44 //======================================================================
45 /// A class for elements that allow the imposition of an applied traction
46 /// in the equations of time-harmonic linear elasticity from a Helmholtz
47 /// potential (interpreted as a displacement potential for the fluid in a
48 /// quasi-steady, linearised FSI problem.)
49 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
50 /// class and and thus, we can be generic enough without the need to have
51 /// a separate equations class.
52 //======================================================================
53  template <class ELASTICITY_BULK_ELEMENT, class HELMHOLTZ_BULK_ELEMENT>
55  : public virtual FaceGeometry<ELASTICITY_BULK_ELEMENT>,
56  public virtual FaceElement,
57  public virtual ElementWithExternalElement
58  {
59 
60  protected:
61 
62  /// \short Pointer to the ratio, \f$ Q \f$ , of the stress used to
63  /// non-dimensionalise the fluid stresses to the stress used to
64  /// non-dimensionalise the solid stresses.
65  double *Q_pt;
66 
67  /// \short Static default value for the ratio of stress scales
68  /// used in the fluid and solid equations (default is 1.0)
69  static double Default_Q_Value;
70 
71  /// Index at which the i-th displacement component is stored
74 
75  /// \short Helper function that actually calculates the residuals
76  // This small level of indirection is required to avoid calling
77  // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
78  // which causes all kinds of pain if overloading later on
80  Vector<double> &residuals);
81 
82  public:
83 
84  /// \short Constructor, which takes a "bulk" element and the
85  /// value of the index and its limit
87  FiniteElement* const &element_pt,
88  const int &face_index) :
89  FaceGeometry<ELASTICITY_BULK_ELEMENT>(), FaceElement(),
91  {
92 
93  //Attach the geometrical information to the element. N.B. This function
94  //also assigns nbulk_value from the required_nvalue of the bulk element
95  element_pt->build_face_element(face_index,this);
96 
97 #ifdef PARANOID
98  {
99  //Check that the element is not a refineable 3d element
100  ELASTICITY_BULK_ELEMENT* elem_pt =
101  dynamic_cast<ELASTICITY_BULK_ELEMENT*>(element_pt);
102  //If it's three-d
103  if(elem_pt->dim()==3)
104  {
105  //Is it refineable
106  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
107  if(ref_el_pt!=0)
108  {
109  if (this->has_hanging_nodes())
110  {
111  throw OomphLibError(
112  "This flux element will not work correctly if nodes are hanging\n",
113  OOMPH_CURRENT_FUNCTION,
114  OOMPH_EXCEPTION_LOCATION);
115  }
116  }
117  }
118  }
119 #endif
120 
121  // Set source element storage: one interaction with an external
122  // element -- the Helmholtz bulk element that provides the pressure
123  this->set_ninteraction(1);
124 
125  //Find the dimension of the problem
126  unsigned n_dim = element_pt->nodal_dimension();
127 
128  //Find the index at which the displacement unknowns are stored
129  ELASTICITY_BULK_ELEMENT* cast_element_pt =
130  dynamic_cast<ELASTICITY_BULK_ELEMENT*>(element_pt);
132  resize(n_dim+1);
133  for(unsigned i=0;i<n_dim+1;i++)
134  {
136  cast_element_pt->u_index_time_harmonic_fourier_decomposed_linear_elasticity(i);
137  }
138  }
139 
140 
141  /// Return the residuals
143  {
145  }
146 
147 
148 
149  /// Fill in contribution from Jacobian
151  DenseMatrix<double> &jacobian)
152  {
153  //Call the residuals
155 
156  //Derivatives w.r.t. external data
158  }
159 
160  /// \short Return the ratio of the stress scales used to non-dimensionalise
161  /// the fluid and elasticity equations. E.g.
162  /// \f$ Q = (\omega a)^2 \rho/E \f$, i.e. the
163  /// ratio between the inertial fluid stress and the solid's elastic
164  /// modulus E.
165  const double &q() const {return *Q_pt;}
166 
167  /// \short Return a pointer the ratio of stress scales used to
168  /// non-dimensionalise the fluid and solid equations.
169  double* &q_pt() {return Q_pt;}
170 
171 
172  /// \short Output function
173  void output(std::ostream &outfile)
174  {
175  /// Dummy
176  unsigned nplot=0;
177  output(outfile,nplot);
178  }
179 
180  /// \short Output function: Plot traction etc at Gauss points
181  /// nplot is ignored.
182  void output(std::ostream &outfile, const unsigned &n_plot)
183  {
184  // Dimension
185  unsigned n_dim = this->nodal_dimension();
186 
187  // Storage for traction
188  Vector<std::complex<double> > traction(n_dim+1);
189 
190  // Get FSI parameter
191  const double q_value=q();
192 
193  outfile << "ZONE\n";
194 
195  //Set the value of n_intpt
196  unsigned n_intpt = integral_pt()->nweight();
197 
198  //Loop over the integration points
199  for(unsigned ipt=0;ipt<n_intpt;ipt++)
200  {
201  Vector<double> s_int(n_dim-1);
202  s_int[0]=integral_pt()->knot(ipt,0);
203 
204  //Get the outer unit normal
205  Vector<double> interpolated_normal(n_dim);
206  outer_unit_normal(ipt,interpolated_normal);
207 
208  // Boundary coordinate
209  Vector<double> zeta(1);
210  interpolated_zeta(s_int,zeta);
211 
212 
213  // Get bulk element for potential
214  HELMHOLTZ_BULK_ELEMENT* ext_el_pt=
215  dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(external_element_pt(0,ipt));
217  std::complex<double> u_helmholtz=
218  ext_el_pt->interpolated_u_fourier_decomposed_helmholtz(s_ext);
219 
220  // Traction: Pressure is proportional to POSITIVE potential
221  ext_el_pt->interpolated_u_fourier_decomposed_helmholtz(s_ext);
222  for(unsigned i=0;i<n_dim;i++)
223  {
224  traction[i]=-q_value*interpolated_normal[i]*u_helmholtz;
225  }
226 
227  outfile << ext_el_pt->interpolated_x(s_ext,0) << " "
228  << ext_el_pt->interpolated_x(s_ext,1) << " "
229  << traction[0] << " "
230  << traction[1] << " "
231  << interpolated_normal[0] << " "
232  << interpolated_normal[1] << " "
233  << u_helmholtz << " "
234  << zeta[0] << std::endl;
235  }
236  }
237 
238  /// \short C_style output function
239  void output(FILE* file_pt)
241 
242  /// \short C-style output function
243  void output(FILE* file_pt, const unsigned &n_plot)
245 
246  };
247 
248 ///////////////////////////////////////////////////////////////////////
249 ///////////////////////////////////////////////////////////////////////
250 ///////////////////////////////////////////////////////////////////////
251 
252 
253 
254 //=================================================================
255 /// Static default value for the ratio of stress scales
256 /// used in the fluid and solid equations (default is 1.0)
257 //=================================================================
258 template <class ELASTICITY_BULK_ELEMENT, class HELMHOLTZ_BULK_ELEMENT>
259 double FourierDecomposedTimeHarmonicLinElastLoadedByHelmholtzPressureBCElement<
260  ELASTICITY_BULK_ELEMENT, HELMHOLTZ_BULK_ELEMENT>::Default_Q_Value=1.0;
261 
262 
263 //=====================================================================
264 /// Return the residuals
265 //=====================================================================
266  template <class ELASTICITY_BULK_ELEMENT, class HELMHOLTZ_BULK_ELEMENT>
267  void FourierDecomposedTimeHarmonicLinElastLoadedByHelmholtzPressureBCElement<
268  ELASTICITY_BULK_ELEMENT, HELMHOLTZ_BULK_ELEMENT>::
269  fill_in_contribution_to_residuals_helmholtz_traction(
270  Vector<double> &residuals)
271  {
272 
273  //Find out how many nodes there are
274  unsigned n_node = nnode();
275 
276 #ifdef PARANOID
277  //Find out how many positional dofs there are
278  unsigned n_position_type = this->nnodal_position_type();
279  if(n_position_type != 1)
280  {
281  throw OomphLibError(
282  "LinearElasticity is not yet implemented for more than one position type.",
283  OOMPH_CURRENT_FUNCTION,
284  OOMPH_EXCEPTION_LOCATION);
285  }
286 #endif
287 
288  //Find out the dimension of the node
289  const unsigned n_dim = this->nodal_dimension();
290 
291  //Cache the nodal indices at which the displacement components are stored
292  std::vector<std::complex<unsigned> > u_nodal_index(n_dim+1);
293  for(unsigned i=0;i<n_dim+1;i++)
294  {
295  u_nodal_index[i] =
296  this->U_index_fourier_decomposed_time_harmonic_linear_elasticity_helmholtz_traction[i];
297  }
298 
299  //Integer to hold the local equation number
300  int local_eqn=0;
301 
302  //Set up memory for the shape functions
303  Shape psi(n_node);
304  DShape dpsids(n_node,n_dim-1);
305 
306  // Get FSI parameter
307  const double q_value=q();
308 
309  // Storage for traction
310  Vector<std::complex<double> > traction(n_dim+1);
311 
312  //Set the value of n_intpt
313  unsigned n_intpt = integral_pt()->nweight();
314 
315  //Loop over the integration points
316  for(unsigned ipt=0;ipt<n_intpt;ipt++)
317  {
318  //Get the integral weight
319  double w = integral_pt()->weight(ipt);
320 
321  //Only need to call the local derivatives
322  dshape_local_at_knot(ipt,psi,dpsids);
323 
324  //Calculate the coordinates
325  Vector<double> interpolated_x(n_dim,0.0);
326 
327  //Also calculate the surface tangent vectors
328  DenseMatrix<double> interpolated_A(n_dim-1,n_dim,0.0);
329 
330  //Calculate displacements and derivatives
331  for(unsigned l=0;l<n_node;l++)
332  {
333  //Loop over directions
334  for(unsigned i=0;i<n_dim;i++)
335  {
336  //Calculate the Eulerian coords
337  const double x_local = nodal_position(l,i);
338  interpolated_x[i] += x_local*psi(l);
339 
340  //Loop over LOCAL derivative directions, to calculate the tangent(s)
341  for(unsigned j=0;j<n_dim-1;j++)
342  {
343  interpolated_A(j,i) += x_local*dpsids(l,j);
344  }
345  }
346  }
347 
348  //First component for the residuals
349  double r=interpolated_x[0];
350 
351  //Now find the local metric tensor from the tangent Vectors
352  DenseMatrix<double> A(n_dim-1);
353  for(unsigned i=0;i<n_dim-1;i++)
354  {
355  for(unsigned j=0;j<n_dim-1;j++)
356  {
357  //Initialise surface metric tensor to zero
358  A(i,j) = 0.0;
359 
360  //Take the dot product
361  for(unsigned k=0;k<n_dim;k++)
362  {
363  A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
364  }
365  }
366  }
367 
368  //Get the outer unit normal
369  Vector<double> interpolated_normal(n_dim);
370  outer_unit_normal(ipt,interpolated_normal);
371 
372  //Find the determinant of the metric tensor
373  double Adet =0.0;
374  switch(n_dim)
375  {
376  case 2:
377  Adet = A(0,0);
378  break;
379  case 3:
380  Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
381  break;
382  default:
383  throw
385  "Wrong dimension in TimeHarmonicLinElastLoadedByPressureElement",
386  "TimeHarmonicLinElastLoadedByPressureElement::fill_in_contribution_to_residuals()",
387  OOMPH_EXCEPTION_LOCATION);
388  }
389 
390  //Premultiply the weights and the square-root of the determinant of
391  //the metric tensor
392  double W = w*sqrt(Adet);
393 
394  // Get bulk element for potential
395  HELMHOLTZ_BULK_ELEMENT* ext_el_pt=
396  dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(external_element_pt(0,ipt));
397  Vector<double> s_ext(external_element_local_coord(0,ipt));
398 
399  // Traction: Pressure is proportional to POSITIVE potential
400  std::complex<double> u_helmholtz=ext_el_pt->interpolated_u_fourier_decomposed_helmholtz(s_ext);
401  for(unsigned i=0;i<n_dim;i++)
402  {
403  traction[i]=-q_value*interpolated_normal[i]*u_helmholtz;
404  }
405  //Loop over the test functions, nodes of the element
406  for(unsigned l=0;l<n_node;l++)
407  {
408  //Loop over the displacement components
409  for(unsigned i=0;i<n_dim+1;i++)
410  {
411 
412  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i].real());
413  /*IF it's not a boundary condition*/
414  if(local_eqn >= 0)
415  {
416  //Add the loading terms to the residuals
417  residuals[local_eqn] -= traction[i].real()*psi(l)*r*W;
418  }
419 
420  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i].imag());
421  /*IF it's not a boundary condition*/
422  if(local_eqn >= 0)
423  {
424  //Add the loading terms to the residuals
425  residuals[local_eqn] -= traction[i].imag()*psi(l)*r*W;
426  }
427 
428  } //End of if not boundary condition
429  } //End of loop over shape functions
430  } //End of loop over integration points
431 
432  }
433 
434 
435 
436 /////////////////////////////////////////////////////////////////////////
437 /////////////////////////////////////////////////////////////////////////
438 /////////////////////////////////////////////////////////////////////////
439 
440 
441 
442 //======================================================================
443 /// \short A class for elements that allow the imposition of an
444 /// prescribed flux (determined from the normal displacements of an
445 /// adjacent linearly elastic solid. Normal derivative for displacement
446 /// potential is given by normal displacement of adjacent solid multiplies
447 /// by FSI parameter (q = k^2 B/E).
448 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
449 /// policy class.
450 //======================================================================
451  template <class HELMHOLTZ_BULK_ELEMENT, class ELASTICITY_BULK_ELEMENT>
453  public virtual FaceGeometry<HELMHOLTZ_BULK_ELEMENT>,
454  public virtual FaceElement,
455  public virtual ElementWithExternalElement
456  {
457 
458  public:
459 
460 
461  /// \short Constructor, takes the pointer to the "bulk" element and the
462  /// face index identifying the face to which the element is attached.
464  FiniteElement* const &bulk_el_pt,
465  const int& face_index);
466 
467  /// Broken copy constructor
470  {
472  "FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement");
473  }
474 
475  /// Broken assignment operator
476 //Commented out broken assignment operator because this can lead to a conflict warning
477 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
478 //realise that two separate implementations of the broken function are the same and so,
479 //quite rightly, it shouts.
480  /*void operator=(
481  const FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement&)
482  {
483  BrokenCopy::broken_assign(
484  "FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement");
485  }*/
486 
487 
488  /// Add the element's contribution to its residual vector
490  {
491  //Call the generic residuals function with flag set to 0
492  //using a dummy matrix argument
495  }
496 
497 
498  /// \short Add the element's contribution to its residual vector and its
499  /// Jacobian matrix
501  DenseMatrix<double> &jacobian)
502  {
503  //Call the generic routine with the flag set to 1
505  (residuals,jacobian,1);
506 
507  //Derivatives w.r.t. external data
509  }
510 
511  /// Output function
512  void output(std::ostream &outfile)
513  {
514  //Dummy
515  unsigned nplot=0;
516  output(outfile,nplot);
517  }
518 
519  /// Output function: flux etc at Gauss points; n_plot is ignored.
520  void output(std::ostream &outfile, const unsigned &n_plot)
521  {
522  outfile << "ZONE\n";
523 
524  //Get the value of Nintpt
525  const unsigned n_intpt = integral_pt()->nweight();
526 
527  //Set the Vector to hold local coordinates
528  Vector<double> s_int(Dim-1);
529 
530  //Loop over the integration points
531  for(unsigned ipt=0;ipt<n_intpt;ipt++)
532  {
533 
534  //Assign values of s
535  for(unsigned i=0;i<(Dim-1);i++)
536  {
537  s_int[i] = integral_pt()->knot(ipt,i);
538  }
539 
540  //Get the unit normal
541  Vector<double> interpolated_normal(Dim);
542  outer_unit_normal(ipt,interpolated_normal);
543 
544  Vector<double> zeta(1);
545  interpolated_zeta(s_int,zeta);
546  // Get displacements
547  ELASTICITY_BULK_ELEMENT* ext_el_pt=dynamic_cast<ELASTICITY_BULK_ELEMENT*>(
548  external_element_pt(0,ipt));
550  Vector<std::complex<double> > displ(Dim+1);
551  ext_el_pt->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(s_ext,displ);
552 
553  // Convert into flux BC: This takes the dot product of the
554  // actual displacement with the flux element's own outer
555  // unit normal so the plus sign is OK.
556  std::complex<double> flux=(displ[0]*interpolated_normal[0]+
557  displ[1]*interpolated_normal[1]);
558 
559  // Output
560  outfile << ext_el_pt->interpolated_x(s_ext,0) << " "
561  << ext_el_pt->interpolated_x(s_ext,1) << " "
562  << flux*interpolated_normal[0] << " "
563  << flux*interpolated_normal[1] << " "
564  << interpolated_normal[0] << " "
565  << interpolated_normal[1] << " "
566  << flux << " "
567  << zeta[0] << std::endl;
568  }
569  }
570 
571 
572  /// C-style output function
573  void output(FILE* file_pt)
575 
576  /// C-style output function
577  void output(FILE* file_pt, const unsigned &n_plot)
579 
580 
581 
582  protected:
583 
584  /// \short Function to compute the shape and test functions and to return
585  /// the Jacobian of mapping between local and global (Eulerian)
586  /// coordinates
587  inline double shape_and_test(const Vector<double> &s, Shape &psi, Shape &test)
588  const
589  {
590  //Find number of nodes
591  unsigned n_node = nnode();
592 
593  //Get the shape functions
594  shape(s,psi);
595 
596  //Set the test functions to be the same as the shape functions
597  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
598 
599  //Return the value of the jacobian
600  return J_eulerian(s);
601  }
602 
603 
604  /// \short Function to compute the shape and test functions and to return
605  /// the Jacobian of mapping between local and global (Eulerian)
606  /// coordinates
607  inline double shape_and_test_at_knot(const unsigned &ipt,
608  Shape &psi, Shape &test)
609  const
610  {
611  //Find number of nodes
612  unsigned n_node = nnode();
613 
614  //Get the shape functions
615  shape_at_knot(ipt,psi);
616 
617  //Set the test functions to be the same as the shape functions
618  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
619 
620  //Return the value of the jacobian
621  return J_eulerian_at_knot(ipt);
622  }
623 
624 
625 
626  private:
627 
628 
629  /// \short Add the element's contribution to its residual vector.
630  /// flag=1(or 0): do (or don't) compute the contribution to the
631  /// Jacobian as well.
633  Vector<double> &residuals, DenseMatrix<double> &jacobian,
634  const unsigned& flag);
635 
636  ///The spatial dimension of the problem
637  unsigned Dim;
638 
639  ///The index at which the unknown is stored at the nodes
640  std::complex<unsigned> U_index_helmholtz_from_displacement;
641 
642 
643  };
644 
645 //////////////////////////////////////////////////////////////////////
646 //////////////////////////////////////////////////////////////////////
647 //////////////////////////////////////////////////////////////////////
648 
649 
650 
651 //===========================================================================
652 /// Constructor, takes the pointer to the "bulk" element, and the
653 /// face index that identifies the face of the bulk element to which
654 /// this face element is to be attached.
655 //===========================================================================
656  template <class HELMHOLTZ_BULK_ELEMENT, class ELASTICITY_BULK_ELEMENT>
658  <HELMHOLTZ_BULK_ELEMENT, ELASTICITY_BULK_ELEMENT>::
660  FiniteElement* const &bulk_el_pt,
661  const int &face_index) :
662  FaceGeometry<HELMHOLTZ_BULK_ELEMENT>(), FaceElement()
663  {
664 
665  // Let the bulk element build the FaceElement, i.e. setup the pointers
666  // to its nodes (by referring to the appropriate nodes in the bulk
667  // element), etc.
668  bulk_el_pt->build_face_element(face_index,this);
669 
670 #ifdef PARANOID
671  {
672  //Check that the element is not a refineable 3d element
673  HELMHOLTZ_BULK_ELEMENT* elem_pt =
674  dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(bulk_el_pt);
675  //If it's three-d
676  if(elem_pt->dim()==3)
677  {
678  //Is it refineable
679  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
680  if(ref_el_pt!=0)
681  {
682  if (this->has_hanging_nodes())
683  {
684  throw OomphLibError(
685  "This flux element will not work correctly if nodes are hanging\n",
686  OOMPH_CURRENT_FUNCTION,
687  OOMPH_EXCEPTION_LOCATION);
688  }
689  }
690  }
691  }
692 #endif
693 
694  // Set source element storage: one interaction with an external element
695  // that provides the displacement of the adjacent linear elasticity
696  // element
697  this->set_ninteraction(1);
698 
699 
700  // Extract the dimension of the problem from the dimension of
701  // the first node
702  Dim = bulk_el_pt->nodal_dimension();//this->node_pt(0)->ndim();// icidav
703 
704  //Set up U_index_helmholtz_displacement. Initialise to zero, which
705  // probably won't change in most cases, oh well, the price we
706  // pay for generality
707  U_index_helmholtz_from_displacement = std::complex<unsigned>(0,0);
708 
709  //Cast to the appropriate HelmholtzEquation so that we can
710  //find the index at which the variable is stored
711  //We assume that the dimension of the full problem is the same
712  //as the dimension of the node, if this is not the case you will have
713  //to write custom elements, sorry
714 
716  dynamic_cast<FourierDecomposedHelmholtzEquations*>(bulk_el_pt);
717  //If the cast has failed die
718  if(eqn_pt==0)
719  {
720  std::string error_string =
721  "Bulk element must inherit from HelmholtzEquations.";
722  error_string +=
723  "Nodes are three dimensional, but cannot cast the bulk element to\n";
724  error_string += "HelmholtzEquations<3>\n.";
725  error_string +=
726  "If you desire this functionality, you must implement it yourself\n";
727 
728  throw OomphLibError(
729  error_string,
730  OOMPH_CURRENT_FUNCTION,
731  OOMPH_EXCEPTION_LOCATION);
732 
733  }
734  else
735  {
736  //Read the index from the (cast) bulk element.
737  U_index_helmholtz_from_displacement = eqn_pt->u_index_fourier_decomposed_helmholtz();
738  }
739 
740 
741  }
742 
743 
744 //===========================================================================
745 /// \short Helper function to compute the element's residual vector and
746 /// the Jacobian matrix.
747 //===========================================================================
748  template <class HELMHOLTZ_BULK_ELEMENT, class ELASTICITY_BULK_ELEMENT>
750  <HELMHOLTZ_BULK_ELEMENT,ELASTICITY_BULK_ELEMENT>::
751  fill_in_generic_residual_contribution_helmholtz_flux_from_displacement(
752  Vector<double> &residuals, DenseMatrix<double> &jacobian,
753  const unsigned& flag)
754  {
755  //Find out how many nodes there are
756  const unsigned n_node = nnode();
757 
758  //Set up memory for the shape and test functions
759  Shape psif(n_node), testf(n_node);
760 
761  //Set the value of Nintpt
762  const unsigned n_intpt = integral_pt()->nweight();
763 
764  //Set the Vector to hold local coordinates
765  Vector<double> s(Dim-1);
766 
767  //Integers to hold the local equation and unknown numbers
768  int local_eqn=0;
769 
770  // Locally cache the index at which the variable is stored
771  const std::complex<unsigned> u_index_helmholtz =
772  U_index_helmholtz_from_displacement;
773 
774  //Loop over the integration points
775  //--------------------------------
776  for(unsigned ipt=0;ipt<n_intpt;ipt++)
777  {
778 
779  //Assign values of s
780  for(unsigned i=0;i<(Dim-1);i++) {s[i] = integral_pt()->knot(ipt,i);}
781 
782  //Get the integral weight
783  double w = integral_pt()->weight(ipt);
784 
785  //Find the shape and test functions and return the Jacobian
786  //of the mapping
787  double J = shape_and_test(s,psif,testf);
788 
789  //Premultiply the weights and the Jacobian
790  double W = w*J;
791 
792  //Need to find position to feed into flux function, initialise to zero
793  Vector<double> interpolated_x(Dim,0.0);
794 
795  //Calculate velocities and derivatives
796  for(unsigned l=0;l<n_node;l++)
797  {
798  //Loop over velocity components
799  for(unsigned i=0;i<Dim;i++)
800  {
801  interpolated_x[i] += nodal_position(l,i)*psif[l];
802  }
803  }
804 
805  //first component
806  double r = interpolated_x[0];
807  //Get the outer unit normal
808  Vector<double> interpolated_normal(Dim);
809  outer_unit_normal(ipt,interpolated_normal);
810 
811  // Get displacements
812  ELASTICITY_BULK_ELEMENT* ext_el_pt=dynamic_cast<ELASTICITY_BULK_ELEMENT*>(
813  external_element_pt(0,ipt));
814  Vector<double> s_ext(external_element_local_coord(0,ipt));
815  Vector<std::complex<double> > displ(Dim+1);
816  ext_el_pt->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(s_ext,displ);
817 
818  // Convert into flux BC: This takes the dot product of the
819  // actual displacement with the flux element's own outer
820  // unit normal so the plus sign is OK.
821  std::complex<double> flux=(displ[0]*interpolated_normal[0]+
822  displ[1]*interpolated_normal[1]);
823 
824  //Now add to the appropriate equations
825 
826  //Loop over the test functions
827  for(unsigned l=0;l<n_node;l++)
828  {
829  local_eqn = nodal_local_eqn(l,u_index_helmholtz.real());
830  /*IF it's not a boundary condition*/
831  if(local_eqn >= 0)
832  {
833  //Add the prescribed flux terms
834  residuals[local_eqn] -= flux.real()*testf[l]*r*W;
835 
836  // Imposed traction doesn't depend upon the solution,
837  // --> the Jacobian is always zero, so no Jacobian
838  // terms are required
839  }
840 
841  local_eqn = nodal_local_eqn(l,u_index_helmholtz.imag());
842  /*IF it's not a boundary condition*/
843  if(local_eqn >= 0)
844  {
845  //Add the prescribed flux terms
846  residuals[local_eqn] -= flux.imag()*testf[l]*r*W;
847 
848  // Imposed traction doesn't depend upon the solution,
849  // --> the Jacobian is always zero, so no Jacobian
850  // terms are required
851  }
852  }
853  }
854  }
855 }
856 
857 
858 #endif
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: Plot traction etc at Gauss points nplot is ignored.
std::complex< unsigned > U_index_helmholtz_from_displacement
The index at which the unknown is stored at the nodes.
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...
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
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
const double & q() const
Return the ratio of the stress scales used to non-dimensionalise the fluid and elasticity equations...
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_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...
FourierDecomposedTimeHarmonicLinElastLoadedByHelmholtzPressureBCElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
Vector< std::complex< unsigned > > U_index_fourier_decomposed_time_harmonic_linear_elasticity_helmholtz_traction
Index at which the i-th displacement component is stored.
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 ...
A general Finite Element class.
Definition: elements.h:1271
void fill_in_contribution_to_residuals_helmholtz_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...
A class for elements that allow the imposition of an prescribed flux (determined from the normal disp...
FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement(const FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement &dummy)
Broken copy constructor.
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 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
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.
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_generic_residual_contribution_helmholtz_flux_from_displacement(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...
double * Q_pt
Pointer to the ratio, , of the stress used to non-dimensionalise the fluid stresses to the stress us...
static char t char * s
Definition: cfortran.h:572
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.
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
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 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.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: flux etc at Gauss points; n_plot is ignored.
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...
double *& q_pt()
Return a pointer the ratio of stress scales used to non-dimensionalise the fluid and solid equations...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
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
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2344
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
FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement(FiniteElement *const &bulk_el_pt, const int &face_index)
Constructor, takes the pointer to the "bulk" element and the face index identifying the face to which...
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