pml_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: 1259 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-11-08 23:52:03 +0000 (Tue, 08 Nov 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_GEN_HH_TIME_HARMONIC_LIN_ELAST_HEADER
31 #define OOMPH_GEN_HH_TIME_HARMONIC_LIN_ELAST_HEADER
32 
33 
34 //Oomph-lib includes
35 #include "generic.h"
36 #include "pml_helmholtz.h"
38 
39 namespace oomph
40 {
41 
42 
43 //======================================================================
44 /// A class for elements that allow the imposition of an applied traction
45 /// in the equations of time-harmonic linear elasticity from a
46 /// PMLHelmholtz potential (interpreted as a displacement potential
47 /// for the fluid in a quasi-steady, linearised FSI problem.)
48 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
49 /// class and thus, we can be generic enough without the need to have
50 /// a separate equations class.
51 //======================================================================
52  template <class ELASTICITY_BULK_ELEMENT, class HELMHOLTZ_BULK_ELEMENT>
54  public virtual FaceGeometry<ELASTICITY_BULK_ELEMENT>,
55  public virtual FaceElement,
56  public virtual ElementWithExternalElement
57  {
58 
59  protected:
60 
61  /// \short Pointer to the ratio, \f$ Q \f$ , of the stress used to
62  /// non-dimensionalise the fluid stresses to the stress used to
63  /// non-dimensionalise the solid stresses.
64  double *Q_pt;
65 
66  /// \short Static default value for the ratio of stress scales
67  /// used in the fluid and solid equations (default is 1.0)
68  static double Default_Q_Value;
69 
70  /// Index at which the i-th displacement component is stored
73 
74  /// \short Helper function that actually calculates the residuals
75  // This small level of indirection is required to avoid calling
76  // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
77  // which causes all kinds of pain if overloading later on
79  Vector<double> &residuals);
80 
81  public:
82 
83  /// \short Constructor, which takes a "bulk" element and the
84  /// value of the index and its limit
86  FiniteElement* const &element_pt,
87  const int &face_index) :
88  FaceGeometry<ELASTICITY_BULK_ELEMENT>(), FaceElement(),
90  {
91 
92  //Attach the geometrical information to the element. N.B. This function
93  //also assigns nbulk_value from the required_nvalue of the bulk element
94  element_pt->build_face_element(face_index,this);
95 
96 #ifdef PARANOID
97  {
98  //Check that the element is not a refineable 3d element
99  ELASTICITY_BULK_ELEMENT* elem_pt =
100  dynamic_cast<ELASTICITY_BULK_ELEMENT*>(element_pt);
101  //If it's three-d
102  if(elem_pt->dim()==3)
103  {
104  //Is it refineable
105  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
106  if(ref_el_pt!=0)
107  {
108  if (this->has_hanging_nodes())
109  {
110  throw OomphLibError(
111  "This flux element will not work correctly if nodes are hanging\n",
112  OOMPH_CURRENT_FUNCTION,
113  OOMPH_EXCEPTION_LOCATION);
114  }
115  }
116  }
117  }
118 #endif
119 
120  // Set source element storage: one interaction with an external
121  // element -- the PMLHelmholtz bulk element that provides
122  // 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);
133  for(unsigned i=0;i<n_dim;i++)
134  {
136  cast_element_pt->u_index_time_harmonic_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);
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  for (unsigned i=0;i<n_dim-1;i++)
203  {
204  s_int[i]=integral_pt()->knot(ipt,i);
205  }
206 
207  //Get the outer unit normal
208  Vector<double> interpolated_normal(n_dim);
209  outer_unit_normal(ipt,interpolated_normal);
210 
211  // Boundary coordinate
212  Vector<double> zeta(1);
213  interpolated_zeta(s_int,zeta);
214 
215  // Get external element for potential
216  HELMHOLTZ_BULK_ELEMENT* ext_el_pt=
217  dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(external_element_pt(0,ipt));
219  std::complex<double> u_helmholtz=
220  ext_el_pt->interpolated_u_pml_helmholtz(s_ext);
221 
222  // Traction: Pressure is proportional to POSITIVE potential
223  ext_el_pt->interpolated_u_pml_helmholtz(s_ext);
224  traction[0]=-q_value*interpolated_normal[0]*u_helmholtz;
225  traction[1]=-q_value*interpolated_normal[1]*u_helmholtz;
226 
227  outfile << ext_el_pt->interpolated_x(s_ext,0) << " "
228  << ext_el_pt->interpolated_x(s_ext,1) << " "
229  << traction[0].real() << " "
230  << traction[1].real() << " "
231  << traction[0].imag() << " "
232  << traction[1].imag() << " "
233  << interpolated_normal[0] << " "
234  << interpolated_normal[1] << " "
235  << u_helmholtz.real() << " "
236  << u_helmholtz.imag() << " "
237  << interpolated_x(s_int,0) << " "
238  << interpolated_x(s_int,1) << " "
239  << sqrt(pow(ext_el_pt->interpolated_x(s_ext,0)-
240  interpolated_x(s_int,0),2)+
241  pow(ext_el_pt->interpolated_x(s_ext,1)-
242  interpolated_x(s_int,1),2)) << " "
243  << zeta[0] << std::endl;
244  }
245  }
246 
247 
248  /// \short C_style output function
249  void output(FILE* file_pt)
251 
252  /// \short C-style output function
253  void output(FILE* file_pt, const unsigned &n_plot)
255 
256  /// \short Compute the global_flux_contribution through the template
257  /// elasticity elements : we compute u.n
258  std::complex<double> global_flux_contribution_from_solid()
259  {
260  // Dummy output file
261  std::ofstream outfile;
262  return global_flux_contribution_from_solid(outfile);
263  }
264 
265  /// \short Compute the element's contribution to the integral of
266  /// the flux over the artificial boundary. Also output the
267  /// flux in the specified output file if it's open.
268  std::complex<double> global_flux_contribution_from_solid
269  (std::ofstream& outfile)
270  {
271  // pointer to the corresponding elasticity bulk element
272  ELASTICITY_BULK_ELEMENT* bulk_elem_pt =
273  dynamic_cast<ELASTICITY_BULK_ELEMENT*>(this->bulk_element_pt());
274 
275  //get the dim of the bulk and local nodes
276  const unsigned bulk_dim= bulk_elem_pt->dim();
277  const unsigned local_dim=this->dim();
278 
279  //Set up memory for the outer unit normal
280  Vector<double> unit_normal(bulk_dim);
281 
282  //Set the value of n_intpt
283  const unsigned n_intpt = integral_pt()->nweight();
284 
285  //Set the Vector to hold local coordinates
286  Vector<double> s(local_dim);
287  std::complex<double> flux(0.0,0.0);
288  Vector<std::complex<double> > u(bulk_dim);
289 
290  // Output?
291  if (outfile.is_open())
292  {
293  outfile << "ZONE\n";
294  }
295 
296  //Loop over the integration points
297  //--------------------------------
298  for(unsigned ipt=0;ipt<n_intpt;ipt++)
299  {
300  //Assign values of s
301  for(unsigned i=0;i<local_dim;i++)
302  {
303  s[i] = integral_pt()->knot(ipt,i);
304  }
305  //get the outer_unit_ext vector
306  this->outer_unit_normal(s,unit_normal);
307 
308  //Get the integral weight
309  double w = integral_pt()->weight(ipt);
310 
311  // Get jacobian of mapping
312  double J=J_eulerian(s);
313 
314  //Premultiply the weights and the Jacobian
315  double W = w*J;
316 
317  // Get local coordinates in bulk element by copy construction
319 
320  // Get the displacement from the bulk element
321  bulk_elem_pt->interpolated_u_time_harmonic_linear_elasticity(s_bulk,u);
322 
323  // Compute normal displacement u.n
324  std::complex<double> u_dot_n(0.0,0.0);
325  for(unsigned i=0;i<bulk_dim;i++)
326  {
327  u_dot_n += u[i]*unit_normal[i];
328  }
329 
330  // Output?
331  if (outfile.is_open())
332  {
333  Vector<double> x(bulk_dim);
334  interpolated_x(s,x);
335  outfile << x[0] << " "
336  << x[1] << " "
337  << u_dot_n.real() << " "
338  << u_dot_n.imag() << "\n";
339  }
340 
341  // ...add to integral
342  flux+=u_dot_n*W;
343  }
344 
345  return flux;
346  }
347 
348 
349  /// \short Compute the global_flux_contribution through the helmholtz
350  /// elements : we compute dphidn.n
352  {
353  // Dummy output file
354  std::ofstream outfile;
356  }
357 
358  /// \short Compute the element's contribution to the integral of
359  /// the flux over the artificial boundary. Also output the
360  /// flux in the specified output file if it's open.
361  std::complex<double> global_flux_contribution_from_helmholtz
362  (std::ofstream& outfile)
363  {
364  // Get the dim of this element
365  const unsigned local_dim=this->dim();
366 
367  //Set the value of n_intpt
368  const unsigned n_intpt = integral_pt()->nweight();
369 
370  //Set the Vector to hold local coordinates
371  Vector<double> s(local_dim);
372  std::complex<double> flux(0.0,0.0);
373 
374  // Output?
375  if (outfile.is_open())
376  {
377  outfile << "ZONE\n";
378  }
379 
380  //Loop over the integration points
381  //--------------------------------
382  for(unsigned ipt=0;ipt<n_intpt;ipt++)
383  {
384  //Assign values of s
385  for(unsigned i=0;i<local_dim;i++)
386  {
387  s[i] = integral_pt()->knot(ipt,i);
388  }
389 
390  // Get local coordinates in bulk element by copy construction
392 
393  // Get external element for potential
394  HELMHOLTZ_BULK_ELEMENT* ext_el_pt=
395  dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(external_element_pt(0,ipt));
396 
397  // number of nodes in ext element
398  unsigned nnode_ext=ext_el_pt->nnode();
399 
400  // get the dim of the external nodes
401  const unsigned ext_dim=ext_el_pt->dim();
402 
403  // Set up memory for the external shape function
404  Shape psi_ext(nnode_ext);
405  DShape dpsi_ext_dx(nnode_ext,ext_dim);
406 
407  //Call the derivatives of the shape functions
408  //in the external element -- must do this via s because this point
409  //is not an integration point the bulk element!
410  (void) ext_el_pt->dshape_eulerian(s_bulk,psi_ext,dpsi_ext_dx);
411 
412  //Set up memory for the outer unit normal
413  Vector<double> unit_normal(ext_dim);
414 
415  //get the outer_unit_ext vector
416  this->outer_unit_normal(s,unit_normal);
417 
418  //Get the integral weight
419  double w = integral_pt()->weight(ipt);
420 
421  // Get jacobian of mapping
422  double J=J_eulerian(s);
423 
424  //Premultiply the weights and the Jacobian
425  double W = w*J;
426 
427  // Global gradient
429  interpolated_dphidx(ext_dim,std::complex<double>(0.0,0.0));
430  for(unsigned l=0;l<nnode_ext;l++)
431  {
432  //Get the nodal value of the helmholtz unknown
433  const std::complex<double> phi_value(
434  ext_el_pt->nodal_value(l,ext_el_pt->u_index_helmholtz().real()),
435  ext_el_pt->nodal_value(l,ext_el_pt->u_index_helmholtz().imag()));
436 
437  //Loop over directions
438  for(unsigned i=0;i<ext_dim;i++)
439  {
440  interpolated_dphidx[i] += phi_value*dpsi_ext_dx(l,i);
441  }
442  }
443 
444  //define dphi_dn
445  std::complex<double> dphi_dn(0.0,0.0);
446  for(unsigned i=0;i<ext_dim;i++)
447  {
448  dphi_dn += interpolated_dphidx[i]*unit_normal[i];
449  }
450 
451 
452 #ifdef PARANOID
453  double max_permitted_discrepancy=1.0e-10;
454  double diff=sqrt(pow(ext_el_pt->interpolated_x(s_bulk,0)-
455  interpolated_x(s,0),2)+
456  pow(ext_el_pt->interpolated_x(s_bulk,1)-
457  interpolated_x(s,1),2));
458  if (diff>max_permitted_discrepancy)
459  {
460  std::ostringstream error_stream;
461  error_stream
462  << "Position computed by external el and face element differ by "
463  << diff << " which is more than the acceptable tolerance "
464  << max_permitted_discrepancy << std::endl;
465  throw OomphLibError(
466  error_stream.str(),
467  OOMPH_CURRENT_FUNCTION,
468  OOMPH_EXCEPTION_LOCATION);
469  }
470 #endif
471 
472  // Output?
473  if (outfile.is_open())
474  {
475  Vector<double> x(ext_dim);
476  interpolated_x(s,x);
477  outfile << x[0] << " "
478  << x[1] << " "
479  << dphi_dn.real() << " "
480  << dphi_dn.imag() << " "
481  << "\n";
482  }
483 
484  // ...add to integral
485  flux+=dphi_dn*W;
486  }
487 
488  return flux;
489  }
490 
491 
492 
493  };
494 
495 
496 ///////////////////////////////////////////////////////////////////////
497 ///////////////////////////////////////////////////////////////////////
498 ///////////////////////////////////////////////////////////////////////
499 
500 
501 
502 //=================================================================
503 /// Static default value for the ragoogletio of stress scales
504 /// used in the fluid and solid equations (default is 1.0)
505 //=================================================================
506 template <class ELASTICITY_BULK_ELEMENT, class HELMHOLTZ_BULK_ELEMENT>
507 double TimeHarmonicLinElastLoadedByPMLHelmholtzPressureBCElement<
508  ELASTICITY_BULK_ELEMENT, HELMHOLTZ_BULK_ELEMENT>::Default_Q_Value=1.0;
509 
510 
511 //=====================================================================
512 /// Return the residuals
513 //=====================================================================
514  template <class ELASTICITY_BULK_ELEMENT, class HELMHOLTZ_BULK_ELEMENT>
515  void TimeHarmonicLinElastLoadedByPMLHelmholtzPressureBCElement<
516  ELASTICITY_BULK_ELEMENT, HELMHOLTZ_BULK_ELEMENT>::
517  fill_in_contribution_to_residuals_helmholtz_traction(
518  Vector<double> &residuals)
519  {
520 
521  //Find out how many nodes there are
522  unsigned n_node = nnode();
523 
524 #ifdef PARANOID
525  //Find out how many positional dofs there are
526  unsigned n_position_type = this->nnodal_position_type();
527  if(n_position_type != 1)
528  {
529  throw OomphLibError(
530  "LinearElasticity is not yet implemented for more than one position type.",
531  OOMPH_CURRENT_FUNCTION,
532  OOMPH_EXCEPTION_LOCATION);
533  }
534 #endif
535 
536  //Find out the dimension of the node
537  const unsigned n_dim = this->nodal_dimension();
538 
539  //Cache the nodal indices at which the displacement components are stored
540  std::complex<unsigned> u_nodal_index[n_dim];
541  for(unsigned i=0;i<n_dim;i++)
542  {
543  u_nodal_index[i] =
544  this->U_index_time_harmonic_linear_elasticity_helmholtz_traction[i];
545  }
546 
547  //Integer to hold the local equation number
548  int local_eqn=0;
549 
550  //Set up memory for the shape functions
551  Shape psi(n_node);
552  DShape dpsids(n_node,n_dim-1);
553 
554  // Get FSI parameter
555  const double q_value=q();
556 
557  // Storage for traction
558  Vector<std::complex<double> > traction(2);
559 
560  //Set the value of n_intpt
561  unsigned n_intpt = integral_pt()->nweight();
562 
563  //Loop over the integration points
564  for(unsigned ipt=0;ipt<n_intpt;ipt++)
565  {
566  //Get the integral weight
567  double w = integral_pt()->weight(ipt);
568 
569  //Only need to call the local derivatives
570  dshape_local_at_knot(ipt,psi,dpsids);
571 
572  //Calculate the coordinates
573  Vector<double> interpolated_x(n_dim,0.0);
574 
575  //Also calculate the surface tangent vectors
576  DenseMatrix<double> interpolated_A(n_dim-1,n_dim,0.0);
577 
578  //Calculate displacements and derivatives
579  for(unsigned l=0;l<n_node;l++)
580  {
581  //Loop over directions
582  for(unsigned i=0;i<n_dim;i++)
583  {
584  //Calculate the Eulerian coords
585  const double x_local = nodal_position(l,i);
586  interpolated_x[i] += x_local*psi(l);
587 
588  //Loop over LOCAL derivative directions, to calculate the tangent(s)
589  for(unsigned j=0;j<n_dim-1;j++)
590  {
591  interpolated_A(j,i) += x_local*dpsids(l,j);
592  }
593  }
594  }
595 
596  //Now find the local metric tensor from the tangent Vectors
597  DenseMatrix<double> A(n_dim-1);
598  for(unsigned i=0;i<n_dim-1;i++)
599  {
600  for(unsigned j=0;j<n_dim-1;j++)
601  {
602  //Initialise surface metric tensor to zero
603  A(i,j) = 0.0;
604 
605  //Take the dot product
606  for(unsigned k=0;k<n_dim;k++)
607  {
608  A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
609  }
610  }
611  }
612 
613  //Get the outer unit normal
614  Vector<double> interpolated_normal(n_dim);
615  outer_unit_normal(ipt,interpolated_normal);
616 
617  //Find the determinant of the metric tensor
618  double Adet =0.0;
619  switch(n_dim)
620  {
621  case 2:
622  Adet = A(0,0);
623  break;
624  case 3:
625  Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
626  break;
627  default:
628  throw
630  "Wrong dimension in TimeHarmonicLinElastLoadedByPressureElement",
631  OOMPH_CURRENT_FUNCTION,
632  OOMPH_EXCEPTION_LOCATION);
633  }
634 
635  //Premultiply the weights and the square-root of the determinant of
636  //the metric tensor
637  double W = w*sqrt(Adet);
638 
639  // Get bulk element for potential
640  HELMHOLTZ_BULK_ELEMENT* ext_el_pt=
641  dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(external_element_pt(0,ipt));
642  Vector<double> s_ext(external_element_local_coord(0,ipt));
643 
644  // Traction: Pressure is proportional to POSITIVE potential
645  std::complex<double> u_helmholtz=ext_el_pt->
646  interpolated_u_pml_helmholtz(s_ext);
647  traction[0]=-q_value*interpolated_normal[0]*u_helmholtz;
648  traction[1]=-q_value*interpolated_normal[1]*u_helmholtz;
649 
650  //Loop over the test functions, nodes of the element
651  for(unsigned l=0;l<n_node;l++)
652  {
653  //Loop over the displacement components
654  for(unsigned i=0;i<n_dim;i++)
655  {
656 
657  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i].real());
658  /*IF it's not a boundary condition*/
659  if(local_eqn >= 0)
660  {
661  //Add the loading terms to the residuals
662  residuals[local_eqn] -= traction[i].real()*psi(l)*W;
663  }
664 
665  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i].imag());
666  /*IF it's not a boundary condition*/
667  if(local_eqn >= 0)
668  {
669  //Add the loading terms to the residuals
670  residuals[local_eqn] -= traction[i].imag()*psi(l)*W;
671  }
672 
673  } //End of if not boundary condition
674  } //End of loop over shape functions
675  } //End of loop over integration points
676 
677  }
678 
679 
680 
681 
682 ///////////////////////////////////////////////////////////////////////
683 ///////////////////////////////////////////////////////////////////////
684 ///////////////////////////////////////////////////////////////////////
685 
686 
687 
688 //======================================================================
689 /// \short A class for elements that allow the imposition of an
690 /// prescribed flux (determined from the normal displacements of an
691 /// adjacent linearly elastic solid. Normal derivative for displacement
692 /// potential is given by normal displacement of adjacent solid multiplies
693 /// by FSI parameter (q = k^2 B/E).
694 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
695 /// policy class.
696 //======================================================================
697  template <class HELMHOLTZ_BULK_ELEMENT, class ELASTICITY_BULK_ELEMENT>
699  public virtual FaceGeometry<HELMHOLTZ_BULK_ELEMENT>,
700  public virtual FaceElement,
701  public virtual ElementWithExternalElement
702  {
703 
704  public:
705 
706 
707  /// \short Constructor, takes the pointer to the "bulk" element and the
708  /// face index identifying the face to which the element is attached.
710  FiniteElement* const &bulk_el_pt,
711  const int& face_index);
712 
713  /// Broken copy constructor
716  {
718  ("PMLHelmholtzFluxFromNormalDisplacementBCElement");
719  }
720 
721  /// Broken assignment operator
722 //Commented out broken assignment operator because this can lead to a conflict warning
723 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
724 //realise that two separate implementations of the broken function are the same and so,
725 //quite rightly, it shouts.
726  /*void operator=(const
727  PMLHelmholtzFluxFromNormalDisplacementBCElement&)
728  {
729  BrokenCopy::broken_assign
730  ("PMLHelmholtzFluxFromNormalDisplacementBCElement");
731  }*/
732 
733 
734  /// Add the element's contribution to its residual vector
736  {
737  //Call the generic residuals function with flag set to 0
738  //using a dummy matrix argument
741  }
742 
743 
744  /// \short Add the element's contribution to its residual vector and its
745  /// Jacobian matrix
747  DenseMatrix<double> &jacobian)
748  {
749  //Call the generic routine with the flag set to 1
751  (residuals,jacobian,1);
752 
753  //Derivatives w.r.t. external data
755  }
756 
757  /// Output function
758  void output(std::ostream &outfile)
759  {
760  //Dummy
761  unsigned nplot=0;
762  output(outfile,nplot);
763  }
764 
765  /// Output function: flux etc at Gauss points; n_plot is ignored.
766  void output(std::ostream &outfile, const unsigned &n_plot)
767  {
768  outfile << "ZONE\n";
769 
770  //Get the value of Nintpt
771  const unsigned n_intpt = integral_pt()->nweight();
772 
773  //Set the Vector to hold local coordinates
774  Vector<double> s_int(Dim-1);
775 
776  //Loop over the integration points
777  for(unsigned ipt=0;ipt<n_intpt;ipt++)
778  {
779 
780  //Assign values of s
781  for(unsigned i=0;i<(Dim-1);i++)
782  {
783  s_int[i] = integral_pt()->knot(ipt,i);
784  }
785 
786  //Get the unit normal
787  Vector<double> interpolated_normal(Dim);
788  outer_unit_normal(ipt,interpolated_normal);
789 
790  Vector<double> zeta(1);
791  interpolated_zeta(s_int,zeta);
792 
793  // Get displacements
794  ELASTICITY_BULK_ELEMENT* ext_el_pt=
795  dynamic_cast<ELASTICITY_BULK_ELEMENT*>(
796  external_element_pt(0,ipt));
798  Vector<std::complex<double> > displ(2);
799  ext_el_pt->interpolated_u_time_harmonic_linear_elasticity(s_ext,displ);
800 
801  // Convert into flux BC: This takes the dot product of the
802  // actual displacement with the flux element's own outer
803  // unit normal so the plus sign is OK.
804  std::complex<double> flux=(displ[0]*interpolated_normal[0]+
805  displ[1]*interpolated_normal[1]);
806 
807  // Output
808  outfile << ext_el_pt->interpolated_x(s_ext,0) << " "
809  << ext_el_pt->interpolated_x(s_ext,1) << " "
810  << flux.real()*interpolated_normal[0] << " "
811  << flux.real()*interpolated_normal[1] << " "
812  << flux.imag()*interpolated_normal[0] << " "
813  << flux.imag()*interpolated_normal[1] << " "
814  << interpolated_normal[0] << " "
815  << interpolated_normal[1] << " "
816  << flux.real() << " "
817  << flux.imag() << " "
818  << interpolated_x(s_int,0) << " "
819  << interpolated_x(s_int,1) << " "
820  << sqrt(pow(ext_el_pt->interpolated_x(s_ext,0)-
821  interpolated_x(s_int,0),2)+
822  pow(ext_el_pt->interpolated_x(s_ext,1)-
823  interpolated_x(s_int,1),2)) << " "
824  << zeta[0] << std::endl;
825  }
826  }
827 
828 
829  /// C-style output function
830  void output(FILE* file_pt)
832 
833  /// C-style output function
834  void output(FILE* file_pt, const unsigned &n_plot)
836 
837 
838 
839  protected:
840 
841  /// \short Function to compute the shape and test functions and to return
842  /// the Jacobian of mapping between local and global (Eulerian)
843  /// coordinates
844  inline double shape_and_test(const Vector<double> &s, Shape &psi,
845  Shape &test)
846  const
847  {
848  //Find number of nodes
849  unsigned n_node = nnode();
850 
851  //Get the shape functions
852  shape(s,psi);
853 
854  //Set the test functions to be the same as the shape functions
855  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
856 
857  //Return the value of the jacobian
858  return J_eulerian(s);
859  }
860 
861 
862  /// \short Function to compute the shape and test functions and to return
863  /// the Jacobian of mapping between local and global (Eulerian)
864  /// coordinates
865  inline double shape_and_test_at_knot(const unsigned &ipt,
866  Shape &psi, Shape &test)
867  const
868  {
869  //Find number of nodes
870  unsigned n_node = nnode();
871 
872  //Get the shape functions
873  shape_at_knot(ipt,psi);
874 
875  //Set the test functions to be the same as the shape functions
876  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
877 
878  //Return the value of the jacobian
879  return J_eulerian_at_knot(ipt);
880  }
881 
882 
883 
884  private:
885 
886 
887  /// \short Add the element's contribution to its residual vector.
888  /// flag=1(or 0): do (or don't) compute the contribution to the
889  /// Jacobian as well.
891  Vector<double> &residuals, DenseMatrix<double> &jacobian,
892  const unsigned& flag);
893 
894  ///The spatial dimension of the problem
895  unsigned Dim;
896 
897  ///The index at which the unknown is stored at the nodes
898  std::complex<unsigned> U_index_helmholtz_from_displacement;
899 
900 
901  };
902 
903 //////////////////////////////////////////////////////////////////////
904 //////////////////////////////////////////////////////////////////////
905 //////////////////////////////////////////////////////////////////////
906 
907 
908 
909 //===========================================================================
910 /// Constructor, takes the pointer to the "bulk" element, and the
911 /// face index that identifies the face of the bulk element to which
912 /// this face element is to be attached.
913 //===========================================================================
914  template <class HELMHOLTZ_BULK_ELEMENT, class ELASTICITY_BULK_ELEMENT>
916  <HELMHOLTZ_BULK_ELEMENT, ELASTICITY_BULK_ELEMENT>::
918  FiniteElement* const &bulk_el_pt,
919  const int &face_index) :
920  FaceGeometry<HELMHOLTZ_BULK_ELEMENT>(), FaceElement()
921  {
922 
923  // Let the bulk element build the FaceElement, i.e. setup the pointers
924  // to its nodes (by referring to the appropriate nodes in the bulk
925  // element), etc.
926  bulk_el_pt->build_face_element(face_index,this);
927 
928 #ifdef PARANOID
929  {
930  //Check that the element is not a refineable 3d element
931  HELMHOLTZ_BULK_ELEMENT* elem_pt =
932  dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(bulk_el_pt);
933  //If it's three-d
934  if(elem_pt->dim()==3)
935  {
936  //Is it refineable
937  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
938  if(ref_el_pt!=0)
939  {
940  if (this->has_hanging_nodes())
941  {
942  throw OomphLibError(
943  "This flux element will not work correctly if nodes are hanging\n",
944  OOMPH_CURRENT_FUNCTION,
945  OOMPH_EXCEPTION_LOCATION);
946  }
947  }
948  }
949  }
950 #endif
951 
952  // Set source element storage: one interaction with an external element
953  // that provides the displacement of the adjacent linear elasticity
954  // element
955  this->set_ninteraction(1);
956 
957  // Extract the dimension of the problem from the dimension of
958  // the first node
959  Dim = this->node_pt(0)->ndim();
960 
961  //Set up U_index_helmholtz_displacement. Initialise to zero, which
962  // probably won't change in most cases, oh well, the price we
963  // pay for generality
964  U_index_helmholtz_from_displacement = std::complex<unsigned>(0,0);
965 
966  //Cast to the appropriate PMLHelmholtzEquation so that we can
967  //find the index at which the variable is stored
968  //We assume that the dimension of the full problem is the same
969  //as the dimension of the node, if this is not the case you will have
970  //to write custom elements, sorry
971  switch(Dim)
972  {
973  //One dimensional problem
974  case 1:
975  {
976  PMLHelmholtzEquations<1>* eqn_pt =
977  dynamic_cast<PMLHelmholtzEquations<1>*>(bulk_el_pt);
978  //If the cast has failed die
979  if(eqn_pt==0)
980  {
981  std::string error_string =
982  "Bulk element must inherit from PMLHelmholtzEquations.";
983  error_string +=
984  "Nodes are one dimensional, but cannot cast the bulk element to\n";
985  error_string += "PMLHelmholtzEquations<1>\n.";
986  error_string +=
987  "If you desire this functionality, you must implement it yourself\n";
988 
989  throw OomphLibError(
990  error_string,
991  OOMPH_CURRENT_FUNCTION,
992  OOMPH_EXCEPTION_LOCATION);
993  }
994  //Otherwise read out the value
995  else
996  {
997  //Read the index from the (cast) bulk element
999  }
1000  }
1001  break;
1002 
1003  //Two dimensional problem
1004  case 2:
1005  {
1006  PMLHelmholtzEquations<2>* eqn_pt =
1007  dynamic_cast<PMLHelmholtzEquations<2>*>(bulk_el_pt);
1008  //If the cast has failed die
1009  if(eqn_pt==0)
1010  {
1011  std::string error_string =
1012  "Bulk element must inherit from PMLHelmholtzEquations.";
1013  error_string +=
1014  "Nodes are two dimensional, but cannot cast the bulk element to\n";
1015  error_string += "PMLHelmholtzEquations<2>\n.";
1016  error_string +=
1017  "If you desire this functionality, you must implement it yourself\n";
1018 
1019  throw OomphLibError(
1020  error_string,
1021  OOMPH_CURRENT_FUNCTION,
1022  OOMPH_EXCEPTION_LOCATION);
1023  }
1024  else
1025  {
1026  //Read the index from the (cast) bulk element.
1028  }
1029  }
1030  break;
1031 
1032  //Three dimensional problem
1033  case 3:
1034  {
1035  PMLHelmholtzEquations<3>* eqn_pt =
1036  dynamic_cast<PMLHelmholtzEquations<3>*>(bulk_el_pt);
1037  //If the cast has failed die
1038  if(eqn_pt==0)
1039  {
1040  std::string error_string =
1041  "Bulk element must inherit from PMLHelmholtzEquations.";
1042  error_string +=
1043  "Nodes are three dimensional, but cannot cast the bulk element to\n";
1044  error_string += "PMLHelmholtzEquations<3>\n.";
1045  error_string +=
1046  "If you desire this functionality, you must implement it yourself\n";
1047 
1048  throw OomphLibError(
1049  error_string,
1050  OOMPH_CURRENT_FUNCTION,
1051  OOMPH_EXCEPTION_LOCATION);
1052 
1053  }
1054  else
1055  {
1056  //Read the index from the (cast) bulk element.
1058  }
1059  }
1060  break;
1061 
1062  //Any other case is an error
1063  default:
1064  std::ostringstream error_stream;
1065  error_stream
1066  << "Dimension of node is " << Dim
1067  << ". It should be 1,2, or 3!" << std::endl;
1068 
1069  throw OomphLibError(
1070  error_stream.str(),
1071  OOMPH_CURRENT_FUNCTION,
1072  OOMPH_EXCEPTION_LOCATION);
1073  break;
1074  }
1075  }
1076 
1077 
1078 //===========================================================================
1079 /// \short Helper function to compute the element's residual vector and
1080 /// the Jacobian matrix.
1081 //===========================================================================
1082  template <class HELMHOLTZ_BULK_ELEMENT, class ELASTICITY_BULK_ELEMENT>
1084  <HELMHOLTZ_BULK_ELEMENT,ELASTICITY_BULK_ELEMENT>::
1085  fill_in_generic_residual_contribution_helmholtz_flux_from_displacement(
1086  Vector<double> &residuals, DenseMatrix<double> &jacobian,
1087  const unsigned& flag)
1088  {
1089  //Find out how many nodes there are
1090  const unsigned n_node = nnode();
1091 
1092  //Set up memory for the shape and test functions
1093  Shape psif(n_node), testf(n_node);
1094 
1095  //Set the value of Nintpt
1096  const unsigned n_intpt = integral_pt()->nweight();
1097 
1098  //Set the Vector to hold local coordinates
1099  Vector<double> s(Dim-1);
1100 
1101  //Integers to hold the local equation and unknown numbers
1102  int local_eqn=0;
1103 
1104  // Locally cache the index at which the variable is stored
1105  const std::complex<unsigned> u_index_helmholtz =
1106  U_index_helmholtz_from_displacement;
1107 
1108  //Loop over the integration points
1109  //--------------------------------
1110  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1111  {
1112 
1113  //Assign values of s
1114  for(unsigned i=0;i<(Dim-1);i++) {s[i] = integral_pt()->knot(ipt,i);}
1115 
1116  //Get the integral weight
1117  double w = integral_pt()->weight(ipt);
1118 
1119  //Find the shape and test functions and return the Jacobian
1120  //of the mapping
1121  double J = shape_and_test(s,psif,testf);
1122 
1123  //Premultiply the weights and the Jacobian
1124  double W = w*J;
1125 
1126  //Need to find position to feed into flux function, initialise to zero
1127  Vector<double> interpolated_x(Dim,0.0);
1128 
1129  //Calculate velocities and derivatives
1130  for(unsigned l=0;l<n_node;l++)
1131  {
1132  //Loop over velocity components
1133  for(unsigned i=0;i<Dim;i++)
1134  {
1135  interpolated_x[i] += nodal_position(l,i)*psif[l];
1136  }
1137  }
1138 
1139  //Get the outer unit normal
1140  Vector<double> interpolated_normal(2);
1141  outer_unit_normal(ipt,interpolated_normal);
1142 
1143  // Get displacements
1144  ELASTICITY_BULK_ELEMENT* ext_el_pt=dynamic_cast<ELASTICITY_BULK_ELEMENT*>(
1145  external_element_pt(0,ipt));
1146  Vector<double> s_ext(external_element_local_coord(0,ipt));
1147  Vector<std::complex<double> > displ(2);
1148  ext_el_pt->interpolated_u_time_harmonic_linear_elasticity(s_ext,displ);
1149 
1150  // Convert into flux BC: This takes the dot product of the
1151  // actual displacement with the flux element's own outer
1152  // unit normal so the plus sign is OK.
1153  std::complex<double> flux=(displ[0]*interpolated_normal[0]+
1154  displ[1]*interpolated_normal[1]);
1155 
1156  //Now add to the appropriate equations
1157 
1158  //Loop over the test functions
1159  for(unsigned l=0;l<n_node;l++)
1160  {
1161  local_eqn = nodal_local_eqn(l,u_index_helmholtz.real());
1162  /*IF it's not a boundary condition*/
1163  if(local_eqn >= 0)
1164  {
1165  //Add the prescribed flux terms
1166  residuals[local_eqn] -= flux.real()*testf[l]*W;
1167 
1168  // Imposed traction doesn't depend upon the solution,
1169  // --> the Jacobian is always zero, so no Jacobian
1170  // terms are required
1171  }
1172 
1173  local_eqn = nodal_local_eqn(l,u_index_helmholtz.imag());
1174  /*IF it's not a boundary condition*/
1175  if(local_eqn >= 0)
1176  {
1177  //Add the prescribed flux terms
1178  residuals[local_eqn] -= flux.imag()*testf[l]*W;
1179 
1180  // Imposed traction doesn't depend upon the solution,
1181  // --> the Jacobian is always zero, so no Jacobian
1182  // terms are required
1183  }
1184  }
1185  }
1186  }
1187 }
1188 
1189 #endif
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
std::complex< unsigned > U_index_helmholtz_from_displacement
The index at which the unknown is stored at the nodes.
static double Default_Q_Value
Static default value for the ratio of stress scales used in the fluid and solid equations (default is...
PMLHelmholtzFluxFromNormalDisplacementBCElement(const PMLHelmholtzFluxFromNormalDisplacementBCElement &dummy)
Broken copy constructor.
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
PMLHelmholtzFluxFromNormalDisplacementBCElement(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...
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...
A general Finite Element class.
Definition: elements.h:1271
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 fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and its Jacobian matrix.
A class for elements that allow the imposition of an prescribed flux (determined from the normal disp...
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
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: flux etc at Gauss points; n_plot is ignored.
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:66
std::complex< double > global_flux_contribution_from_solid()
Compute the global_flux_contribution through the template elasticity elements : we compute u...
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 ...
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.
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
TimeHarmonicLinElastLoadedByPMLHelmholtzPressureBCElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: Plot traction etc at Gauss points nplot is ignored.
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
Vector< std::complex< unsigned > > U_index_time_harmonic_linear_elasticity_helmholtz_traction
Index at which the i-th displacement component is stored.
void output(std::ostream &outfile)
Output with default number of plot points.
virtual std::complex< unsigned > u_index_helmholtz() const
Broken assignment operator.
std::complex< double > global_flux_contribution_from_helmholtz()
Compute the global_flux_contribution through the helmholtz elements : we compute dphidn.n.
double *& q_pt()
Return a pointer the ratio of stress scales used to non-dimensionalise the fluid and solid equations...
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.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
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.
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4470
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
Vector< double > local_coordinate_in_bulk(const Vector< double > &s) const
Return vector of local coordinates in bulk element, given the local coordinates in this FaceElement...
Definition: elements.cc:6032
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
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.
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
double * Q_pt
Pointer to the ratio, , of the stress used to non-dimensionalise the fluid stresses to the stress us...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Broken assignment operator.
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_residuals_helmholtz_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.