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