helmholtz_bc_elements.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1282 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-01-16 08:27:53 +0000 (Mon, 16 Jan 2017) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 // Header file for elements that are used to apply Sommerfeld
31 // boundary conditions to the Helmholtz equations
32 #ifndef OOMPH_HELMHOLTZ_BC_ELEMENTS_HEADER
33 #define OOMPH_HELMHOLTZ_BC_ELEMENTS_HEADER
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37 #include <oomph-lib-config.h>
38 #endif
39 
40 #include "math.h"
41 #include <complex>
42 
43 // Get the Bessel functions
44 #include "oomph_crbond_bessel.h"
45 
46 
47 namespace oomph
48 {
49 ///////////////////////////////////////////////////////////////////////
50 ///////////////////////////////////////////////////////////////////////
51 // Collection of the Bessel functions used in the Helmholtz problem
52 ///////////////////////////////////////////////////////////////////////
53 ///////////////////////////////////////////////////////////////////////
54 
55 //====================================================================
56 /// Namespace to provide Hankel function of the first kind
57 /// and various orders -- needed for Helmholtz computations.
58 //====================================================================
59 namespace Hankel_functions_for_helmholtz_problem
60 {
61 
62 //====================================================================
63 /// Compute Hankel function of the first kind of orders 0...n and
64 /// its derivates at coordinate x. The function returns the vector
65 /// then its derivative.
66 //====================================================================
67  void Hankel_first(const unsigned& n, const double& x,
68  Vector<std::complex<double> >& h,
69  Vector<std::complex<double> >& hp)
70  {
71  int n_actual = 0;
72  Vector<double> jn(n+1),yn(n+1),jnp(n+1), ynp(n+1);
73  CRBond_Bessel::bessjyna(int(n),x,n_actual,&jn[0],&yn[0],
74  &jnp[0],&ynp[0]);
75 #ifdef PARANOID
76  if (n_actual!=int(n))
77  {
78  std::ostringstream error_stream;
79  error_stream << "CRBond_Bessel::bessjyna() only computed "
80  << n_actual << " rather than " << n
81  << " Bessel functions.\n";
82  throw OomphLibError(error_stream.str(),
83  OOMPH_CURRENT_FUNCTION,
84  OOMPH_EXCEPTION_LOCATION);
85  }
86 #endif
87  for (unsigned i=0;i<n;i++)
88  {
89  h[i] =std::complex<double>(jn[i], yn[i]);
90  hp[i]=std::complex<double>(jnp[i],ynp[i]);
91  }
92  } // End of Hankel_first
93 
94 //====================================================================
95 /// Compute Hankel function of the first kind of orders 0...n and
96 /// its derivates at coordinate x. The function returns the vector
97 /// then its derivative (complex version). This functionality is only
98 /// required in the computation of the solution for the complex-
99 /// shifted Laplacian preconditioner.
100 //====================================================================
101  void CHankel_first(const unsigned& n,
102  const std::complex<double>& x,
103  Vector<std::complex<double> >& h,
104  Vector<std::complex<double> >& hp)
105  {
106  // Set the highest order actually calculated
107  int n_actual=0;
108 
109  // Create a vector for the Bessel function of the 1st kind
110  Vector<std::complex<double> > jn(n+1);
111 
112  // Create a vector for the Bessel function of the 2nd kind
113  Vector<std::complex<double> > yn(n+1);
114 
115  // Create a vector for the Bessel function (1st kind) derivative
116  Vector<std::complex<double> > jnp(n+1);
117 
118  // Create a vector for the Bessel function (2nd kind) derivative
119  Vector<std::complex<double> > ynp(n+1);
120 
121  // Call the (complex) Bessel function to calculate the solution
122  CRBond_Bessel::cbessjyna(int(n),x,n_actual,&jn[0],
123  &yn[0],&jnp[0],&ynp[0]);
124 
125 #ifdef PARANOID
126  // Tell the user if they tried to calculate higher order terms
127  if (n_actual!=int(n))
128  {
129  // Create an output stream
130  std::ostringstream error_message_stream;
131 
132  // Create the error message
133  error_message_stream << "CRBond_Bessel::cbessjyna() only computed "
134  << n_actual << " rather than " << n
135  << " Bessel functions.\n";
136 
137  // Throw the error message
138  throw OomphLibError(error_message_stream.str(),
139  OOMPH_CURRENT_FUNCTION,
140  OOMPH_EXCEPTION_LOCATION);
141  }
142 #endif
143 
144  // Loop over the number of terms requested (only the first entry
145  // has *actually* been calculated)
146  for (unsigned i=0;i<n;i++)
147  {
148  // Set the entries in the Hankel function vector
149  h[i]=std::complex<double>(jn[i].real()-yn[i].imag(),
150  jn[i].imag()+yn[i].real());
151 
152  // Set the entries in the Hankel function derivative vector
153  hp[i]=std::complex<double>(jnp[i].real()-ynp[i].imag(),
154  jnp[i].imag()+ynp[i].real());
155  }
156  } // End of Hankel_first
157 } // End of namespace
158 
159 
160 
161 /////////////////////////////////////////////////////////////////////
162 /////////////////////////////////////////////////////////////////////
163 /////////////////////////////////////////////////////////////////////
164 
165 
166 //======================================================================
167 /// \short A class for elements that allow the approximation of the
168 /// Sommerfeld radiation BC.
169 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
170 /// policy class.
171 //======================================================================
172 template <class ELEMENT>
173  class HelmholtzBCElementBase : public virtual FaceGeometry<ELEMENT>,
174  public virtual FaceElement
175  {
176  public:
177 
178  /// \short Constructor, takes the pointer to the "bulk" element and the
179  /// index of the face to which the element is attached.
180  HelmholtzBCElementBase(FiniteElement* const &bulk_el_pt,
181  const int& face_index);
182 
183  ///\short Broken empty constructor
185  {
186  throw OomphLibError(
187  "Don't call empty constructor for HelmholtzBCElementBase",
188  OOMPH_CURRENT_FUNCTION,
189  OOMPH_EXCEPTION_LOCATION);
190  }
191 
192  /// Broken copy constructor
194  {
195  BrokenCopy::broken_copy("HelmholtzBCElementBase");
196  }
197 
198  /// Broken assignment operator
199 //Commented out broken assignment operator because this can lead to a conflict warning
200 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
201 //realise that two separate implementations of the broken function are the same and so,
202 //quite rightly, it shouts.
203  /*void operator=(const HelmholtzBCElementBase&)
204  {
205  BrokenCopy::broken_assign("HelmholtzBCElementBase");
206  }*/
207 
208 
209  /// \short Specify the value of nodal zeta from the face geometry
210  /// The "global" intrinsic coordinate of the element when
211  /// viewed as part of a geometric object should be given by
212  /// the FaceElement representation, by default (needed to break
213  /// indeterminacy if bulk element is SolidElement)
214  double zeta_nodal(const unsigned &n, const unsigned &k,
215  const unsigned &i) const
216  {return FaceElement::zeta_nodal(n,k,i);}
217 
218 
219  /// Output function -- forward to broken version in FiniteElement
220  /// until somebody decides what exactly they want to plot here...
221  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
222 
223  /// \short Output function -- forward to broken version in FiniteElement
224  /// until somebody decides what exactly they want to plot here...
225  void output(std::ostream &outfile, const unsigned &n_plot)
226  {FiniteElement::output(outfile,n_plot);}
227 
228  /// C-style output function -- forward to broken version in FiniteElement
229  /// until somebody decides what exactly they want to plot here...
230  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
231 
232  /// \short C-style output function -- forward to broken version in
233  /// FiniteElement until somebody decides what exactly they want to plot
234  /// here...
235  void output(FILE* file_pt, const unsigned &n_plot)
236  {FiniteElement::output(file_pt,n_plot);}
237 
238  /// \short Return the index at which the real/imag unknown value
239  /// is stored.
240  virtual inline std::complex<unsigned> u_index_helmholtz() const
241  {
242  return std::complex<unsigned>(U_index_helmholtz.real(),
243  U_index_helmholtz.imag());
244  }
245 
246  /// \short Compute the element's contribution to the time-averaged
247  /// radiated power over the artificial boundary
249  {
250  // Dummy output file
251  std::ofstream outfile;
252  return global_power_contribution(outfile);
253  }
254 
255  /// \short Compute the element's contribution to the time-averaged
256  /// radiated power over the artificial boundary. Also output the
257  /// power density as a fct of the polar angle in the specified
258  ///output file if it's open.
259  double global_power_contribution(std::ofstream& outfile)
260  {
261  // pointer to the corresponding bulk element
262  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
263 
264  // Number of nodes in bulk element
265  unsigned nnode_bulk=bulk_elem_pt->nnode();
266  const unsigned n_node_local = nnode();
267 
268  //get the dim of the bulk and local nodes
269  const unsigned bulk_dim= bulk_elem_pt->dim();
270  const unsigned local_dim= this->dim();
271 
272  //Set up memory for the shape and test functions
273  Shape psi(n_node_local);
274 
275  //Set up memory for the shape functions
276  Shape psi_bulk(nnode_bulk);
277  DShape dpsi_bulk_dx(nnode_bulk,bulk_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  double power=0.0;
288 
289  // Output?
290  if (outfile.is_open())
291  {
292  outfile << "ZONE\n";
293  }
294 
295  //Loop over the integration points
296  //--------------------------------
297  for(unsigned ipt=0;ipt<n_intpt;ipt++)
298  {
299  //Assign values of s
300  for(unsigned i=0;i<local_dim;i++)
301  {
302  s[i] = integral_pt()->knot(ipt,i);
303  }
304  //get the outer_unit_ext vector
305  this->outer_unit_normal(s,unit_normal);
306 
307  //Get the integral weight
308  double w = integral_pt()->weight(ipt);
309 
310  // Get jacobian of mapping
311  double J=J_eulerian(s);
312 
313  //Premultiply the weights and the Jacobian
314  double W = w*J;
315 
316  // Get local coordinates in bulk element by copy construction
318 
319  //Call the derivatives of the shape functions
320  //in the bulk -- must do this via s because this point
321  //is not an integration point the bulk element!
322  (void)bulk_elem_pt->dshape_eulerian(s_bulk,psi_bulk,dpsi_bulk_dx);
323  this->shape(s,psi);
324 
325  // Derivs of Eulerian coordinates w.r.t. local coordinates
326  std::complex<double> dphi_dn(0.0,0.0);
327  Vector<std::complex <double> > interpolated_dphidx(bulk_dim);
328  std::complex<double> interpolated_phi(0.0,0.0);
329  Vector<double> x(bulk_dim);
330 
331  //Calculate function value and derivatives:
332  //-----------------------------------------
333  // Loop over nodes
334  for(unsigned l=0;l<nnode_bulk;l++)
335  {
336  //Get the nodal value of the helmholtz unknown
337  const std::complex<double> phi_value(
338  bulk_elem_pt->nodal_value(l,bulk_elem_pt->u_index_helmholtz().real()),
339  bulk_elem_pt->nodal_value(l,bulk_elem_pt->u_index_helmholtz().imag()));
340 
341  //Loop over directions
342  for(unsigned i=0;i<bulk_dim;i++)
343  {
344  interpolated_dphidx[i] += phi_value*dpsi_bulk_dx(l,i);
345  }
346  } // End of loop over the bulk_nodes
347 
348  for(unsigned l=0;l<n_node_local;l++)
349  {
350  //Get the nodal value of the helmholtz unknown
351  const std::complex<double> phi_value(
352  nodal_value(l,u_index_helmholtz().real()),
353  nodal_value(l,u_index_helmholtz().imag()));
354 
355  interpolated_phi += phi_value*psi(l);
356  }
357 
358  //define dphi_dn
359  for(unsigned i=0;i<bulk_dim;i++)
360  {
361  dphi_dn += interpolated_dphidx[i]*unit_normal[i];
362  }
363 
364  // Power density
365  double integrand=0.5*
366  (interpolated_phi.real()*dphi_dn.imag()-
367  interpolated_phi.imag()*dphi_dn.real());
368 
369  // Output?
370  if (outfile.is_open())
371  {
372  interpolated_x(s,x);
373  double phi=atan2(x[1],x[0]);
374  outfile << x[0] << " "
375  << x[1] << " "
376  << phi << " "
377  << integrand << "\n";
378  }
379 
380  // ...add to integral
381  power+=integrand*W;
382  }
383 
384  return power;
385  }
386 
387 
388 
389  /// \short Compute element's contribution to Fourier components of the
390  /// solution -- length of vector indicates number of terms to be computed.
392  Vector<std::complex<double> >& a_coeff_pos,
393  Vector<std::complex<double> >& a_coeff_neg)
394  {
395 
396 #ifdef PARANOID
397  if (a_coeff_pos.size()!=a_coeff_neg.size())
398  {
399  std::ostringstream error_stream;
400  error_stream << "a_coeff_pos and a_coeff_neg must have "
401  << "the same size. \n";
402  throw OomphLibError(
403  error_stream.str(),
404  OOMPH_CURRENT_FUNCTION,
405  OOMPH_EXCEPTION_LOCATION);
406  }
407 #endif
408 
409  // define the imaginary number
410  const std::complex<double> I(0.0,1.0);
411 
412  //Find out how many nodes there are
413  const unsigned n_node = this->nnode();
414 
415  //Set up memory for the shape functions
416  Shape psi(n_node);
417  DShape dpsi(n_node,1);
418 
419  //Set the value of n_intpt
420  const unsigned n_intpt=this->integral_pt()->nweight();
421 
422  //Set the Vector to hold local coordinates
423  Vector<double> s(this->Dim-1);
424 
425  // Initialise
426  unsigned n=a_coeff_pos.size();
427  for (unsigned i=0;i<n;i++)
428  {
429  a_coeff_pos[i]=std::complex<double>(0.0,0.0);
430  a_coeff_neg[i]=std::complex<double>(0.0,0.0);
431  }
432 
433  //Loop over the integration points
434  //--------------------------------
435  for(unsigned ipt=0;ipt<n_intpt;ipt++)
436  {
437  //Assign values of s
438  for(unsigned i=0;i<(this->Dim-1);i++)
439  {
440  s[i]=this->integral_pt()->knot(ipt,i);
441  }
442 
443  //Get the integral weight
444  double w=this->integral_pt()->weight(ipt);
445 
446  // Get the shape functions
447  this->dshape_local(s,psi,dpsi);
448 
449  // Eulerian coordinates at Gauss point
450  Vector<double> interpolated_x(this->Dim,0.0);
451 
452  // Derivs of Eulerian coordinates w.r.t. local coordinates
453  Vector<double> interpolated_dxds(this->Dim);
454  std::complex<double> interpolated_u(0.0,0.0);
455 
456  // Assemble x and its derivs
457  for(unsigned l=0;l<n_node;l++)
458  {
459  //Loop over directions
460  for(unsigned i=0;i<this->Dim;i++)
461  {
462  interpolated_x[i]+=this->nodal_position(l,i)*psi[l];
463  interpolated_dxds[i]+=this->nodal_position(l,i)*dpsi(l,0);
464  }
465 
466  //Get the nodal value of the helmholtz unknown
467  const std::complex<double> u_value(
468  this->nodal_value(l,this->U_index_helmholtz.real()),
469  this->nodal_value(l,this->U_index_helmholtz.imag()));
470 
471  //Add to the interpolated value
472  interpolated_u += u_value*psi(l);
473  } // End of loop over the nodes
474 
475  // calculate the integral
476  //-----------------------
477 
478  // Get polar angle
479  double phi=atan2(interpolated_x[1],interpolated_x[0]);
480 
481  //define dphi_ds=(-yx'+y'x)/(x^2+y^2)
482  double denom =(interpolated_x[0]*interpolated_x[0])+
483  (interpolated_x[1]*interpolated_x[1]);
484  double nom =-interpolated_dxds[1]*interpolated_x[0]+
485  interpolated_dxds[0]*interpolated_x[1];
486  double dphi_ds=std::fabs(nom/denom);
487 
488  // Positive coefficients
489  for (unsigned i=0;i<n;i++)
490  {
491  a_coeff_pos[i]+=interpolated_u*exp(-I*phi*double(i))*dphi_ds*w;
492  }
493  // Negative coefficients
494  for (unsigned i=1;i<n;i++)
495  {
496  a_coeff_neg[i]+=interpolated_u*exp(I*phi*double(i))*dphi_ds*w;
497  }
498 
499  }//End of loop over integration points
500 
501  }
502 
503 
504 
505 
506  protected:
507 
508  /// \short Function to compute the shape and test functions and to return
509  /// the Jacobian of mapping between local and global (Eulerian)
510  /// coordinates
511  inline double test_only(const Vector<double> &s, Shape &test) const
512  {
513  //Get the shape functions
514  shape(s,test);
515 
516  //Return the value of the jacobian
517  return J_eulerian(s);
518  }
519 
520  /// \short Function to compute the shape, test functions and derivates
521  /// and to return
522  /// the Jacobian of mapping between local and global (Eulerian)
523  /// coordinates
524  inline double d_shape_and_test_local(const Vector<double> &s, Shape &psi,
525  Shape &test,
526  DShape &dpsi_ds,DShape &dtest_ds)
527  const
528  {
529  //Find number of nodes
530  unsigned n_node = nnode();
531 
532  //Get the shape functions
533  dshape_local(s,psi,dpsi_ds);
534 
535  //Set the test functions to be the same as the shape functions
536  for(unsigned i=0;i<n_node;i++)
537  {
538  for(unsigned j=0;j<(Dim-1);j++)
539  {
540  test[i] = psi[i];
541  dtest_ds(i,j)= dpsi_ds(i,j);
542  }
543  }
544  //Return the value of the jacobian
545  return J_eulerian(s);
546  }
547 
548  /// \short The index at which the real and imag part of the unknown is stored
549  /// at the nodes
550  std::complex<unsigned> U_index_helmholtz;
551 
552  ///The spatial dimension of the problem
553  unsigned Dim;
554 
555 
556  };
557 
558 //////////////////////////////////////////////////////////////////////
559 //////////////////////////////////////////////////////////////////////
560 //////////////////////////////////////////////////////////////////////
561 
562 
563 ///=================================================================
564 /// Mesh for DtN boundary condition elements -- provides
565 /// functionality to apply Sommerfeld radiation condtion
566 /// via DtN BC
567 ///=================================================================
568 template<class ELEMENT>
569  class HelmholtzDtNMesh : public virtual Mesh
570 {
571  public:
572 
573  /// Constructor: Specify radius of outer boundary and number of
574  /// Fourier terms used in the computation of the gamma integral
576  const unsigned& nfourier_terms) :
577  Outer_radius(outer_radius), Nfourier_terms(nfourier_terms)
578  {}
579 
580  /// \short Compute and store the gamma integral at all integration
581  /// points of the constituent elements.
582  void setup_gamma();
583 
584  /// \short Compute Fourier components of the solution -- length of
585  /// vector indicates number of terms to be computed.
586  void compute_fourier_components(Vector<std::complex<double> >& a_coeff_pos,
587  Vector<std::complex<double> >& a_coeff_neg);
588 
589  /// \short Gamma integral evaluated at Gauss points
590  /// for specified element
592  {
593  return Gamma_at_gauss_point[el_pt];
594  }
595 
596  /// \short Derivative of Gamma integral w.r.t global unknown, evaluated
597  /// at Gauss points for specified element
600  {
601  return D_Gamma_at_gauss_point[el_pt];
602  }
603 
604  /// \short The outer radius
605  double &outer_radius()
606  {
607  return Outer_radius ;
608  }
609 
610  /// \short Number of Fourier terms used in the computation of the
611  /// gamma integral
612  unsigned& nfourier_terms()
613  {
614  return Nfourier_terms;
615  }
616 
617  private:
618 
619  /// Outer radius
620  double Outer_radius;
621 
622 /// Nbr of Fourier terms used in the Gamma computation
623  unsigned Nfourier_terms;
624 
625 
626  /// \short Container to store the gamma integral for given Gauss point
627  /// and element
628  std::map<FiniteElement*,Vector<std::complex<double> > > Gamma_at_gauss_point;
629 
630 
631  /// \short Container to store the derivate of Gamma integral w.r.t
632  /// global unknown evaluated at Gauss points for specified element
633  std::map<FiniteElement*,Vector<std::map<unsigned,std::complex<double> > > >
635 
636 };
637 
638 
639 
640 /////////////////////////////////////////////////////////////////////
641 /////////////////////////////////////////////////////////////////////
642 /////////////////////////////////////////////////////////////////////
643 
644 
645 //=============================================================
646 /// Absorbing BC element for approximation imposition of
647 /// Sommerfeld radiation condition
648 //==============================================================
649 template<class ELEMENT>
651 {
652 
653  public:
654 
655 
656  /// \short Construct element from specification of bulk element and
657  /// face index.
659  const int& face_index) :
660  HelmholtzBCElementBase<ELEMENT>(bulk_el_pt,face_index)
661  {
662  // Initialise pointers
663  Outer_radius_pt=0;
664 
665  // Initialise order of absorbing boundary condition
666  ABC_order_pt=0;
667  }
668 
669  /// Pointer to order of absorbing boundary condition
670  unsigned*& abc_order_pt()
671  {
672  return ABC_order_pt;
673  }
674 
675 
676  /// Add the element's contribution to its residual vector
678  {
679  //Call the generic residuals function with flag set to 0
680  //using a dummy matrix argument
683  }
684 
685  /// \short Add the element's contribution to its residual vector and its
686  /// Jacobian matrix
688  DenseMatrix<double> &jacobian)
689  {
690  //Call the generic routine with the flag set to 1
692  }
693 
694 
695  /// Get pointer to radius of outer boundary (must be a cirle)
696  double*& outer_radius_pt()
697  {
698  return Outer_radius_pt;
699  }
700 
701  private:
702 
703  /// \short Compute the element's residual vector and the (
704  /// Jacobian matrix.
705  /// Overloaded version, using the abc approximation
707  Vector<double> &residuals, DenseMatrix<double> &jacobian,
708  const unsigned& flag)
709  {
710 
711 #ifdef PARANOID
712 
713  if (ABC_order_pt==0)
714  {
715  throw OomphLibError(
716  "Order of ABC hasn't been set!",
717  OOMPH_CURRENT_FUNCTION,
718  OOMPH_EXCEPTION_LOCATION);
719  }
720 
721  if (this->Outer_radius_pt==0)
722  {
723  throw OomphLibError(
724  "Pointer to outer radius hasn't been set!",
725  OOMPH_CURRENT_FUNCTION,
726  OOMPH_EXCEPTION_LOCATION);
727  }
728 
729 #endif
730 
731  //Find out how many nodes there are
732  const unsigned n_node = this->nnode();
733 
734  //Set up memory for the shape and test functions
735  Shape psi(n_node), test(n_node);
736  DShape dpsi_ds (n_node,this->Dim-1),
737  dtest_ds(n_node,this->Dim-1),
738  dtest_dS(n_node,this->Dim-1),
739  dpsi_dS (n_node,this->Dim-1);
740 
741  //Set the value of Nintpt
742  const unsigned n_intpt = this->integral_pt()->nweight();
743 
744  //Set the Vector to hold local coordinates
745  Vector<double> s(this->Dim-1);
746 
747  //Integers to hold the local equation and unknown numbers
748  int local_eqn_real=0,local_unknown_real=0;
749  int local_eqn_imag=0,local_unknown_imag=0;
750 
751  // Define the problem parameters
752  double R=*Outer_radius_pt;
753  double k=sqrt(dynamic_cast<ELEMENT*>(this->bulk_element_pt())->k_squared());
754 
755  //Loop over the integration points
756  //--------------------------------
757  for(unsigned ipt=0;ipt<n_intpt;ipt++)
758  {
759  //Assign values of s
760  for(unsigned i=0;i<(this->Dim-1);i++)
761  {
762  s[i] = this->integral_pt()->knot(ipt,i);
763  }
764 
765  //Get the integral weight
766  double w = this->integral_pt()->weight(ipt);
767 
768  //Find the shape test functions and derivates; return the Jacobian
769  //of the mapping between local and global (Eulerian)
770  // coordinates
771  double J = this->d_shape_and_test_local(s,psi,test,dpsi_ds,dtest_ds);
772 
773  //Premultiply the weights and the Jacobian
774  double W = w*J;
775  // get the inverse of jacibian
776  double inv_J=1/J;
777 
778  //Need to find position to feed into flux function,
779  //initialise to zero
780  std::complex<double> interpolated_u(0.0,0.0);
781  std::complex<double> du_dS(0.0,0.0);
782 
783  //Calculate velocities and derivatives:loop over the nodes
784  for(unsigned l=0;l<n_node;l++)
785  {
786  // Loop over real and imag part
787  //Get the nodal value of the helmholtz unknown
788  const std::complex<double> u_value(
789  this->nodal_value(l,this->U_index_helmholtz.real()),
790  this->nodal_value(l,this->U_index_helmholtz.imag()));
791 
792  interpolated_u += u_value*psi[l];
793 
794  du_dS += u_value*dpsi_ds(l,0)*inv_J;
795 
796  // Get the value of dtest_dS
797  dtest_dS(l,0)=dtest_ds(l,0)*inv_J;
798  // Get the value of dpsif_dS
799  dpsi_dS(l,0)=dpsi_ds(l,0)*inv_J;
800  }
801 
802  // use ABC first order approximation
803  if (*ABC_order_pt==1)
804  {
805  //Now add to the appropriate equations:use second order approximation
806  //Loop over the test functions
807  for(unsigned l=0;l<n_node;l++)
808  {
809  local_eqn_real = this->nodal_local_eqn
810  (l,this->U_index_helmholtz.real());
811  local_eqn_imag = this->nodal_local_eqn
812  (l,this->U_index_helmholtz.imag());
813 
814  // first, calculate the real part contrubution
815  //-----------------------
816  //IF it's not a boundary condition
817  if(local_eqn_real >= 0)
818  {
819  //Add the second order terms:compute manually real and imag part
820 
821  residuals[local_eqn_real] +=
822  (k*interpolated_u.imag()+(0.5/R)
823  *interpolated_u.real())*test[l]*W;
824 
825  // Calculate the jacobian
826  //-----------------------
827  if(flag)
828  {
829  //Loop over the shape functions again
830  for(unsigned l2=0;l2<n_node;l2++)
831  {
832  local_unknown_real = this->nodal_local_eqn(
833  l2,this->U_index_helmholtz.real());
834  local_unknown_imag = this->nodal_local_eqn(
835  l2,this->U_index_helmholtz.imag());
836  //If at a non-zero degree of freedom add in the entry
837  if(local_unknown_real >= 0)
838  {
839  jacobian(local_eqn_real,local_unknown_real)
840  +=(0.5/R)*psi[l2]*test[l]*W;
841  }
842 
843  //If at a non-zero degree of freedom add in the entry
844  if(local_unknown_imag >= 0)
845  {
846  jacobian(local_eqn_real,local_unknown_imag)
847  += k *psi[l2]*test[l]*W;
848 
849  }
850  }
851  }
852  }// end of local_eqn_real
853 
854  // second, calculate the imag part contrubution
855  //-----------------------
856  //IF it's not a boundary condition
857  if(local_eqn_imag >= 0)
858  {
859  //Add the second order terms contibution to the residual
860  residuals[local_eqn_imag] +=
861  (-k*interpolated_u.real()+(0.5/R)
862  *interpolated_u.imag())*test[l]*W;
863 
864 
865  // Calculate the jacobian
866  //-----------------------
867  if(flag)
868  {
869  //Loop over the shape functions again
870  for(unsigned l2=0;l2<n_node;l2++)
871  {
872  local_unknown_real = this->nodal_local_eqn(
873  l2,this->U_index_helmholtz.real());
874  local_unknown_imag = this->nodal_local_eqn(
875  l2,this->U_index_helmholtz.imag());
876  //If at a non-zero degree of freedom add in the entry
877  if(local_unknown_real >= 0)
878  {
879  jacobian(local_eqn_imag,local_unknown_real)
880  += (-k) *psi[l2]*test[l]*W;
881  }
882 
883  //If at a non-zero degree of freedom add in the entry
884  if(local_unknown_imag >= 0)
885  {
886  jacobian(local_eqn_imag,local_unknown_imag)
887  +=(0.5/R)*psi[l2]*test[l]*W;
888  }
889  }
890  }
891  }
892  }// End of loop over the nodes
893  }
894 
895  //:use second order approximation
896  if(*ABC_order_pt==2)
897  {
898  //Now add to the appropriate equations:use second order approximation
899  //Loop over the test functions
900  for(unsigned l=0;l<n_node;l++)
901  {
902  local_eqn_real = this->nodal_local_eqn
903  (l,this->U_index_helmholtz.real());
904  local_eqn_imag = this->nodal_local_eqn
905  (l,this->U_index_helmholtz.imag());
906 
907  // first, calculate the real part contrubution
908  //-----------------------
909  //IF it's not a boundary condition
910  if(local_eqn_real >= 0)
911  {
912  //Add the second order terms:compute manually real and imag part
913 
914  residuals[local_eqn_real] +=
915  (k*interpolated_u.imag()+(0.5/R)
916  *interpolated_u.real())*test[l]*W
917  + ((0.125/(k*R*R))*interpolated_u.imag())*test[l]*W;
918 
919  residuals[local_eqn_real] +=
920  (-0.5/k)*du_dS.imag()*dtest_dS(l,0)*W;
921 
922  // Calculate the jacobian
923  //-----------------------
924  if(flag)
925  {
926  //Loop over the shape functions again
927  for(unsigned l2=0;l2<n_node;l2++)
928  {
929  local_unknown_real = this->nodal_local_eqn(
930  l2,this->U_index_helmholtz.real());
931  local_unknown_imag = this->nodal_local_eqn(
932  l2,this->U_index_helmholtz.imag());
933  //If at a non-zero degree of freedom add in the entry
934  if(local_unknown_real >= 0)
935  {
936  jacobian(local_eqn_real,local_unknown_real)
937  +=(0.5/R)*psi[l2]*test[l]*W;
938  }
939 
940  //If at a non-zero degree of freedom add in the entry
941  if(local_unknown_imag >= 0)
942  {
943  jacobian(local_eqn_real,local_unknown_imag)
944  += k *psi[l2]*test[l]*W
945  + (0.125/(k*R*R))*psi[l2]*test[l]*W;
946 
947  jacobian(local_eqn_real,local_unknown_imag)
948  +=(-0.5/k)*dpsi_dS(l2,0)*dtest_dS(l,0)*W;
949  }
950  }
951  }
952  }// end of local_eqn_real
953 
954  // second, calculate the imag part contrubution
955  //-----------------------
956  //IF it's not a boundary condition
957  if(local_eqn_imag >= 0)
958  {
959  //Add the second order terms contibution to the residual
960  residuals[local_eqn_imag] +=
961  (-k*interpolated_u.real()+(0.5/R)
962  *interpolated_u.imag())*test[l]*W
963  + ((-0.125/(k*R*R))*interpolated_u.real())*test[l]*W;
964 
965  residuals[local_eqn_imag] +=
966  (0.5/k)*du_dS.real()*dtest_dS(l,0)*W;
967 
968  // Calculate the jacobian
969  //-----------------------
970  if(flag)
971  {
972  //Loop over the shape functions again
973  for(unsigned l2=0;l2<n_node;l2++)
974  {
975  local_unknown_real = this->nodal_local_eqn(
976  l2,this->U_index_helmholtz.real());
977  local_unknown_imag = this->nodal_local_eqn(
978  l2,this->U_index_helmholtz.imag());
979  //If at a non-zero degree of freedom add in the entry
980  if(local_unknown_real >= 0)
981  {
982  jacobian(local_eqn_imag,local_unknown_real)
983  += (-k) *psi[l2]*test[l]*W
984  -(0.125/(k*R*R))*psi[l2]*test[l]*W;
985 
986  jacobian(local_eqn_imag,local_unknown_real)
987  +=(0.5/k)*dpsi_dS(l2,0)*dtest_dS(l,0)*W;
988  }
989 
990  //If at a non-zero degree of freedom add in the entry
991  if(local_unknown_imag >= 0)
992  {
993  jacobian(local_eqn_imag,local_unknown_imag)
994  +=(0.5/R)*psi[l2]*test[l]*W;
995  }
996  }
997  }
998  }
999  }// End of loop over the nodes
1000  }
1001 
1002  if(*ABC_order_pt==3)
1003  {
1004  //Now add to the appropriate equations:use second order approximation
1005  //Loop over the test functions
1006  for(unsigned l=0;l<n_node;l++)
1007  {
1008  local_eqn_real = this->nodal_local_eqn
1009  (l,this->U_index_helmholtz.real());
1010  local_eqn_imag = this->nodal_local_eqn
1011  (l,this->U_index_helmholtz.imag());
1012 
1013  // first, calculate the real part contrubution
1014  //-----------------------
1015  //IF it's not a boundary condition
1016  if(local_eqn_real >= 0)
1017  {
1018  //Add the second order terms:compute manually real and imag part
1019  residuals[local_eqn_real] +=
1020  ((k*(1+0.125/(k*k*R*R)))*interpolated_u.imag()
1021  +(0.5/R-0.125/(k*k*R*R*R))*interpolated_u.real())
1022  *test[l]*W;
1023 
1024  residuals[local_eqn_real] +=
1025  ((-0.5/k)*du_dS.imag()+(0.5/(k*k*R))*du_dS.real())
1026  *dtest_dS(l,0)*W;
1027 
1028  // Calculate the jacobian
1029  //-----------------------
1030  if(flag)
1031  {
1032  //Loop over the shape functions again
1033  for(unsigned l2=0;l2<n_node;l2++)
1034  {
1035  local_unknown_real = this->nodal_local_eqn(
1036  l2,this->U_index_helmholtz.real());
1037  local_unknown_imag = this->nodal_local_eqn(
1038  l2,this->U_index_helmholtz.imag());
1039  //If at a non-zero degree of freedom add in the entry
1040  if(local_unknown_real >= 0)
1041  {
1042  jacobian(local_eqn_real,local_unknown_real)
1043  +=(0.5/R-0.125/(k*k*R*R*R))*psi[l2]*test[l]*W;
1044 
1045  jacobian(local_eqn_real,local_unknown_real)
1046  +=(0.5/(k*k*R))*dpsi_dS(l2,0)*dtest_dS(l,0)*W;
1047  }
1048 
1049  //If at a non-zero degree of freedom add in the entry
1050  if(local_unknown_imag >= 0)
1051  {
1052  jacobian(local_eqn_real,local_unknown_imag)
1053  +=(k*(1+0.125/(k*k*R*R)))*psi[l2]*test[l]*W;
1054 
1055  jacobian(local_eqn_real,local_unknown_imag)
1056  +=(-0.5/k)*dpsi_dS(l2,0)*dtest_dS(l,0)*W;
1057  }
1058  }
1059  }
1060  }// end of local_eqn_real
1061 
1062  // second, calculate the imag part contrubution
1063  //-----------------------
1064  //IF it's not a boundary condition
1065  if(local_eqn_imag >= 0)
1066  {
1067  //Add the second order terms contibution to the residual
1068  residuals[local_eqn_imag] +=
1069  ((-k*(1+0.125/(k*k*R*R)))*interpolated_u.real()
1070  +(0.5/R-0.125/(k*k*R*R*R))*interpolated_u.imag())
1071  *test[l]*W;
1072 
1073  residuals[local_eqn_imag] +=
1074  ((0.5/k)*du_dS.real()+(0.5/(k*k*R))*du_dS.imag())
1075  *dtest_dS(l,0)*W;
1076 
1077  // Calculate the jacobian
1078  //-----------------------
1079  if(flag)
1080  {
1081  //Loop over the shape functions again
1082  for(unsigned l2=0;l2<n_node;l2++)
1083  {
1084  local_unknown_real = this->nodal_local_eqn(
1085  l2,this->U_index_helmholtz.real());
1086  local_unknown_imag = this->nodal_local_eqn(
1087  l2,this->U_index_helmholtz.imag());
1088  //If at a non-zero degree of freedom add in the entry
1089  if(local_unknown_real >= 0)
1090  {
1091  jacobian(local_eqn_imag,local_unknown_real)
1092  +=(-k*(1+0.125/(k*k*R*R)))*psi[l2]*test[l]*W;
1093 
1094  jacobian(local_eqn_imag,local_unknown_real)
1095  +=(0.5/k)*dpsi_dS(l2,0)*dtest_dS(l,0)*W;
1096  }
1097  //If at a non-zero degree of freedom add in the entry
1098  if(local_unknown_imag >= 0)
1099  {
1100  jacobian(local_eqn_imag,local_unknown_imag)
1101  +=(0.5/R-0.125/(k*k*R*R*R))*psi[l2]*test[l]*W;
1102 
1103  jacobian(local_eqn_imag,local_unknown_imag)
1104  +=(0.5/(k*k*R))*dpsi_dS(l2,0)*dtest_dS(l,0)*W;
1105  }
1106  }
1107  }
1108  }
1109  }// End of loop over the nodes
1110  }
1111  } //End of loop over int_pt
1112 
1113  } // End of fill_in_generic_residual_contribution_helmholtz_flux
1114 
1115  private:
1116 
1117  /// Pointer to radius of outer boundary (must be a circle!)
1119 
1120  /// Pointer to order of absorbing boundary condition
1121  unsigned* ABC_order_pt;
1122 
1123 };
1124 
1125 //=============================================================
1126 /// FaceElement used to apply Sommerfeld radiation conditon
1127 /// via Dirichlet to Neumann map.
1128 //==============================================================
1129 template<class ELEMENT>
1131 {
1132 
1133  public:
1134 
1135  /// \short Construct element from specification of bulk element and
1136  /// face index.
1138  const int& face_index) :
1139  HelmholtzBCElementBase<ELEMENT>(bulk_el_pt,face_index)
1140  {}
1141 
1142  /// Add the element's contribution to its residual vector
1144  {
1145  //Call the generic residuals function with flag set to 0
1146  //using a dummy matrix argument
1148  (residuals,GeneralisedElement::Dummy_matrix,0);
1149  }
1150 
1151  /// \short Add the element's contribution to its residual vector and its
1152  /// Jacobian matrix
1154  DenseMatrix<double> &jacobian)
1155  {
1156  //Call the generic routine with the flag set to 1
1158  (residuals,jacobian,1);
1159  }
1160 
1161  /// \short Compute the contribution of the element
1162  /// to the Gamma integral and its derivates w.r.t
1163  /// to global unknows; the function takes the wavenumber and the polar
1164  /// angle phi as input
1165  void compute_gamma_contribution(const double& phi,const int& n,
1166  std::complex<double>& gamma_con,
1167  std::map<unsigned,std::complex<double> >&
1168  d_gamma_con);
1169 
1170 
1171  /// \short Access function to mesh of all DtN boundary condition elements
1172  /// (needed to get access to gamma values)
1174  {
1175  return Outer_boundary_mesh_pt;
1176  }
1177 
1178  /// \short Set mesh of all DtN boundary condition elements
1181  {
1182  Outer_boundary_mesh_pt=mesh_pt;
1183  }
1184 
1185 
1186  /// \short Complete the setup of additional dependencies arising
1187  /// through the far-away interaction with other nodes in
1188  /// Outer_boundary_mesh_pt.
1190  {
1191  // Create a set of all nodes
1192  std::set<Node*> node_set;
1193  unsigned nel=Outer_boundary_mesh_pt->nelement();
1194  for (unsigned e=0;e<nel;e++)
1195  {
1196  FiniteElement* el_pt=Outer_boundary_mesh_pt->finite_element_pt(e);
1197  unsigned nnod=el_pt->nnode();
1198  for (unsigned j=0;j<nnod;j++)
1199  {
1200  Node* nod_pt=el_pt->node_pt(j);
1201 
1202  // Don't add copied nodes
1203  if (!(nod_pt->is_a_copy()))
1204  {
1205  node_set.insert(nod_pt);
1206  }
1207  }
1208  }
1209  // Now erase the current element's own nodes
1210  unsigned nnod=this->nnode();
1211  for (unsigned j=0;j<nnod;j++)
1212  {
1213  Node* nod_pt=this->node_pt(j);
1214  node_set.erase(nod_pt);
1215 
1216  // If the element's node is a copy then its "master" will
1217  // already have been added in the set above -- remove the
1218  // master to avoid double counting eqn numbers
1219  if (nod_pt->is_a_copy())
1220  {
1221  node_set.erase(nod_pt->copied_node_pt());
1222  }
1223  }
1224 
1225  // Now declare these nodes to be the element's external Data
1226  for (std::set<Node*>::iterator it=node_set.begin();
1227  it!=node_set.end();it++)
1228  {
1229  this->add_external_data(*it);
1230  }
1231  }
1232 
1233 
1234  private:
1235 
1236  /// \short Compute the element's residual vector
1237  /// Jacobian matrix.
1238  /// Overloaded version, using the gamma computed in the mesh
1240  (Vector<double> &residuals, DenseMatrix<double> &jacobian,
1241  const unsigned& flag)
1242  {
1243  //Find out how many nodes there are
1244  const unsigned n_node = this->nnode();
1245 
1246  //Set up memory for the shape and test functions
1247  Shape test(n_node);
1248 
1249  //Set the value of Nintpt
1250  const unsigned n_intpt = this->integral_pt()->nweight();
1251 
1252  //Set the Vector to hold local coordinates
1253  Vector<double> s(this->Dim-1);
1254 
1255  //Integers to hold the local equation and unknown numbers
1256  int local_eqn_real=0,local_unknown_real=0,global_eqn_real=0,
1257  local_eqn_imag=0,local_unknown_imag=0,global_eqn_imag=0;
1258  int external_global_eqn_real=0, external_unknown_real=0,
1259  external_global_eqn_imag=0, external_unknown_imag=0;
1260 
1261 
1262  // Get the gamma value for the current integration point
1263  // from the mesh
1265  gamma(Outer_boundary_mesh_pt->gamma_at_gauss_point(this));
1266 
1268  d_gamma(Outer_boundary_mesh_pt->d_gamma_at_gauss_point(this));
1269 
1270  //Loop over the integration points
1271  //--------------------------------
1272  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1273  {
1274  //Assign values of s
1275  for(unsigned i=0;i<(this->Dim-1);i++)
1276  {
1277  s[i] = this->integral_pt()->knot(ipt,i);
1278  }
1279 
1280  //Get the integral weight
1281  double w = this->integral_pt()->weight(ipt);
1282 
1283  //Find the shape test functions and derivates; return the Jacobian
1284  //of the mapping between local and global (Eulerian)
1285  // coordinates
1286  double J = this->test_only(s,test);
1287 
1288  //Premultiply the weights and the Jacobian
1289  double W = w*J;
1290 
1291  //Now add to the appropriate equations
1292  //Loop over the test functions:loop over the nodes
1293  for(unsigned l=0;l<n_node;l++)
1294  {
1295  local_eqn_real = this->nodal_local_eqn
1296  (l,this->U_index_helmholtz.real());
1297  local_eqn_imag = this->nodal_local_eqn
1298  (l,this->U_index_helmholtz.imag());
1299 
1300  //IF it's not a boundary condition
1301  if(local_eqn_real >= 0)
1302  {
1303  //Add the gamma contribution in this int_point to the res
1304  residuals[local_eqn_real] -=gamma[ipt].real()*test[l]*W;
1305 
1306  // Calculate the jacobian
1307  //-----------------------
1308  if(flag)
1309  {
1310  //Loop over the shape functions again
1311  for(unsigned l2=0;l2<n_node;l2++)
1312  {
1313  // Add the contribution of the local data
1314  local_unknown_real = this->nodal_local_eqn(
1315  l2,this->U_index_helmholtz.real());
1316 
1317  local_unknown_imag = this->nodal_local_eqn(
1318  l2,this->U_index_helmholtz.imag());
1319 
1320  //If at a non-zero degree of freedom add in the entry
1321  if(local_unknown_real >= 0)
1322  {
1323  global_eqn_real=this->eqn_number(local_unknown_real);
1324 
1325  // Add the first order terms contribution
1326  jacobian(local_eqn_real,local_unknown_real)
1327  -=d_gamma[ipt][global_eqn_real].real()*test[l]*W;
1328  }
1329  if(local_unknown_imag >= 0)
1330  {
1331  global_eqn_imag=this->eqn_number(local_unknown_imag);
1332 
1333  // Add the first order terms contribution
1334  jacobian(local_eqn_real,local_unknown_imag)
1335  -=d_gamma[ipt][global_eqn_imag].real()*test[l]*W;
1336  }
1337  } // End of loop over nodes l2
1338 
1339  // Add the contribution of the external data
1340  unsigned n_ext_data=this->nexternal_data();
1341  //Loop over the shape functions again
1342  for(unsigned l2=0;l2<n_ext_data;l2++)
1343  {
1344  // Add the contribution of the local data
1345  external_unknown_real = this->external_local_eqn(
1346  l2,this->U_index_helmholtz.real());
1347 
1348  external_unknown_imag = this->external_local_eqn(
1349  l2,this->U_index_helmholtz.imag());
1350 
1351  //If at a non-zero degree of freedom add in the entry
1352  if(external_unknown_real >= 0)
1353  {
1354  external_global_eqn_real=this->eqn_number(external_unknown_real);
1355 
1356  // Add the first order terms contribution
1357  jacobian(local_eqn_real,external_unknown_real)
1358  -=d_gamma[ipt][external_global_eqn_real].real()*test[l]*W;
1359  }
1360  if(external_unknown_imag >= 0)
1361  {
1362  external_global_eqn_imag=this->eqn_number(external_unknown_imag);
1363 
1364  // Add the first order terms contribution
1365  jacobian(local_eqn_real,external_unknown_imag)
1366  -=d_gamma[ipt][external_global_eqn_imag].real()*test[l]*W;
1367  }
1368  } // End of loop over external data
1369  }// End of flag
1370  }// end of local_eqn_real
1371 
1372  if(local_eqn_imag >= 0)
1373  {
1374  //Add the gamma contribution in this int_point to the res
1375  residuals[local_eqn_imag] -=gamma[ipt].imag()*test[l]*W;
1376 
1377  // Calculate the jacobian
1378  //-----------------------
1379  if(flag)
1380  {
1381  //Loop over the shape functions again
1382  for(unsigned l2=0;l2<n_node;l2++)
1383  {
1384  // Add the contribution of the local data
1385  local_unknown_real = this->nodal_local_eqn(
1386  l2,this->U_index_helmholtz.real());
1387 
1388  local_unknown_imag = this->nodal_local_eqn(
1389  l2,this->U_index_helmholtz.imag());
1390 
1391  //If at a non-zero degree of freedom add in the entry
1392  if(local_unknown_real >= 0)
1393  {
1394  global_eqn_real=this->eqn_number(local_unknown_real);
1395 
1396  // Add the first order terms contribution
1397  jacobian(local_eqn_imag,local_unknown_real)
1398  -=d_gamma[ipt][global_eqn_real].imag()*test[l]*W;
1399  }
1400  if(local_unknown_imag >= 0)
1401  {
1402  global_eqn_imag=this->eqn_number(local_unknown_imag);
1403 
1404  // Add the first order terms contribution
1405  jacobian(local_eqn_imag,local_unknown_imag)
1406  -=d_gamma[ipt][global_eqn_imag].imag()*test[l]*W;
1407  }
1408  } // End of loop over nodes l2
1409 
1410  // Add the contribution of the external data
1411  unsigned n_ext_data=this->nexternal_data();
1412  //Loop over the shape functions again
1413  for(unsigned l2=0;l2<n_ext_data;l2++)
1414  {
1415  // Add the contribution of the local data
1416  external_unknown_real = this->external_local_eqn(
1417  l2,this->U_index_helmholtz.real());
1418 
1419  external_unknown_imag = this->external_local_eqn(
1420  l2,this->U_index_helmholtz.imag());
1421 
1422  //If at a non-zero degree of freedom add in the entry
1423  if(external_unknown_real >= 0)
1424  {
1425  external_global_eqn_real=this->eqn_number(external_unknown_real);
1426 
1427  // Add the first order terms contribution
1428  jacobian(local_eqn_imag,external_unknown_real)
1429  -=d_gamma[ipt][external_global_eqn_real].imag()*test[l]*W;
1430  }
1431  if(external_unknown_imag >= 0)
1432  {
1433  external_global_eqn_imag=this->eqn_number(external_unknown_imag);
1434 
1435  // Add the first order terms contribution
1436  jacobian(local_eqn_imag,external_unknown_imag)
1437  -=d_gamma[ipt][external_global_eqn_imag].imag()*test[l]*W;
1438  }
1439  } // End of loop over external data
1440  }// End of flag
1441  } // end of local_eqn_imag
1442  }// end of llop over yhe node
1443  } //End of loop over int_pt
1444  } // End of fill_in_generic_residual_contribution_helmholtz_flux
1445 
1446 
1447  /// \short Pointer to mesh of all DtN boundary condition elements
1448  /// (needed to get access to gamma values)
1450 
1451 };
1452 
1453 
1454 ////////////////////////////////////////////////////////////////
1455 ////////////////////////////////////////////////////////////////
1456 ////////////////////////////////////////////////////////////////
1457 
1458 
1459 
1460 //===========start_compute_gamma_contribution==================
1461 /// \short compute the contribution of the element
1462 /// to the Gamma integral and its derivates w.r.t
1463 /// to global unknows; the function takes wavenumber n
1464 /// and polar angle phi as input.
1465 //==============================================================
1466 template<class ELEMENT>
1469  const double& phi,
1470  const int& n,
1471  std::complex<double>& gamma_con,
1472  std::map<unsigned,std::complex<double> >& d_gamma_con)
1473  {
1474  // define the imaginary number
1475  const std::complex<double> I(0.0,1.0);
1476 
1477  //Find out how many nodes there are
1478  const unsigned n_node = this->nnode();
1479 
1480  //Set up memory for the shape functions
1481  Shape psi(n_node);
1482  DShape dpsi(n_node,1);
1483 
1484  // initialise the variable
1485  int local_unknown_real=0, local_unknown_imag=0;
1486  int global_unknown_real=0,global_unknown_imag=0;
1487 
1488  //Set the value of n_intpt
1489  const unsigned n_intpt=this->integral_pt()->nweight();
1490 
1491  //Set the Vector to hold local coordinates
1492  Vector<double> s(this->Dim-1);
1493 
1494  // Initialise
1495  gamma_con=std::complex<double>(0.0,0.0);
1496  d_gamma_con.clear();
1497 
1498  //Loop over the integration points
1499  //--------------------------------
1500  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1501  {
1502  //Assign values of s
1503  for(unsigned i=0;i<(this->Dim-1);i++)
1504  {
1505  s[i]=this->integral_pt()->knot(ipt,i);
1506  }
1507 
1508  //Get the integral weight
1509  double w=this->integral_pt()->weight(ipt);
1510 
1511  // Get the shape functions
1512  this->dshape_local(s,psi,dpsi);
1513 
1514  // Eulerian coordinates at Gauss point
1515  Vector<double> interpolated_x(this->Dim,0.0);
1516 
1517  // Derivs of Eulerian coordinates w.r.t. local coordinates
1518  Vector<double> interpolated_dxds(this->Dim);
1519  std::complex<double> interpolated_u(0.0,0.0);
1520 
1521  // Assemble x and its derivs
1522  for(unsigned l=0;l<n_node;l++)
1523  {
1524  //Loop over directions
1525  for(unsigned i=0;i<this->Dim;i++)
1526  {
1527  interpolated_x[i]+=this->nodal_position(l,i)*psi[l];
1528  interpolated_dxds[i]+=this->nodal_position(l,i)*dpsi(l,0);
1529  }
1530 
1531  //Get the nodal value of the helmholtz unknown
1532  std::complex<double> u_value(
1533  this->nodal_value(l,this->U_index_helmholtz.real()),
1534  this->nodal_value(l,this->U_index_helmholtz.imag()));
1535 
1536  interpolated_u += u_value*psi(l);
1537  } // End of loop over the nodes
1538 
1539  // calculate the integral
1540  //-----------------------
1541  // define the variable phi_p
1542  double phi_p=atan2(interpolated_x[1],interpolated_x[0]);
1543 
1544  //define dphi_ds=(-yx'+y'x)/(x^2+y^2)
1545  double denom =(interpolated_x[0]*interpolated_x[0])+
1546  (interpolated_x[1]*interpolated_x[1]);
1547  double nom =-interpolated_dxds[1]*interpolated_x[0]+
1548  interpolated_dxds[0]*interpolated_x[1];
1549  double dphi_ds=std::fabs(nom/denom);
1550 
1551  // compute the element contribution to gamma
1552  // ALH: The awkward construction with pow and the static_cast is to
1553  // avoid a floating point error on my machine when running unoptimised
1554  // (no idea why!)
1555  gamma_con+=(dphi_ds)*w*pow(exp(I*(phi-phi_p)),static_cast<double>(n))
1556  *interpolated_u;
1557 
1558  // compute the contribution to each node to the map
1559  for(unsigned l=0;l<n_node;l++)
1560  {
1561  // Add the contribution of the real local data
1562  local_unknown_real = this->nodal_local_eqn(
1563  l,this->U_index_helmholtz.real());
1564  if (local_unknown_real >= 0)
1565  {
1566  global_unknown_real=this->eqn_number(local_unknown_real);
1567  d_gamma_con[global_unknown_real]+=
1568  (dphi_ds)*w*exp(I*(phi-phi_p)*double(n))*psi(l);
1569  }
1570 
1571  // Add the contribution of the imag local data
1572  local_unknown_imag = this->nodal_local_eqn(
1573  l,this->U_index_helmholtz.imag());
1574  if (local_unknown_imag >= 0)
1575  {
1576  global_unknown_imag=this->eqn_number(local_unknown_imag);
1577  // ALH: The awkward construction with pow and the static_cast is to
1578  // avoid a floating point error on my machine when running unoptimised
1579  // (no idea why!)
1580  d_gamma_con[global_unknown_imag]+=
1581  I* (dphi_ds)*w*pow(exp(I*(phi-phi_p)),static_cast<double>(n))*psi(l);
1582  }
1583  }// end of loop over the node
1584  }//End of loop over integration points
1585 
1586  }
1587 
1588 
1589 /////////////////////////////////////////////////////////////////////
1590 /////////////////////////////////////////////////////////////////////
1591 /////////////////////////////////////////////////////////////////////
1592 
1593 
1594 //===========================================================================
1595 /// \short Namespace for checking radius of nodes on (assumed to be circular)
1596 /// DtN boundary
1597 //===========================================================================
1598 namespace ToleranceForHelmholtzOuterBoundary
1599 {
1600  /// \short Relative tolerance to within radius of points on DtN boundary
1601  /// are allowed to deviate from specified value
1602  extern double Tol;
1603 
1604 }
1605 
1606 
1607 //===========================================================================
1608 /// \short Namespace for checking radius of nodes on (assumed to be circular)
1609 /// DtN boundary
1610 //===========================================================================
1611 namespace ToleranceForHelmholtzOuterBoundary
1612 {
1613  /// \short Relative tolerance to within radius of points on DtN boundary
1614  /// are allowed to deviate from specified value
1615  double Tol=1.0e-3;
1616 
1617 }
1618 
1619 /////////////////////////////////////////////////////////////////////
1620 /////////////////////////////////////////////////////////////////////
1621 /////////////////////////////////////////////////////////////////////
1622 
1623 
1624 ///================================================================
1625  /// \short Compute Fourier components of the solution -- length of
1626  /// vector indicates number of terms to be computed.
1627 ///================================================================
1628 template<class ELEMENT>
1630  Vector<std::complex<double> >& a_coeff_pos,
1631  Vector<std::complex<double> >& a_coeff_neg)
1632 {
1633 
1634 #ifdef PARANOID
1635  if (a_coeff_pos.size()!=a_coeff_neg.size())
1636  {
1637  std::ostringstream error_stream;
1638  error_stream << "a_coeff_pos and a_coeff_neg must have "
1639  << "the same size. \n";
1640  throw OomphLibError(
1641  error_stream.str(),
1642  OOMPH_CURRENT_FUNCTION,
1643  OOMPH_EXCEPTION_LOCATION);
1644  }
1645 #endif
1646 
1647  // Initialise
1648  unsigned n=a_coeff_pos.size();
1649  Vector<std::complex<double> > el_a_coeff_pos(n);
1650  Vector<std::complex<double> > el_a_coeff_neg(n);
1651  for (unsigned i=0;i<n;i++)
1652  {
1653  a_coeff_pos[i]=std::complex<double>(0.0,0.0);
1654  a_coeff_neg[i]=std::complex<double>(0.0,0.0);
1655  }
1656 
1657  //Loop over elements e
1658  unsigned nel=this->nelement();
1659  for (unsigned e=0;e<nel;e++)
1660  {
1661  // Get a pointer to element
1663  dynamic_cast<HelmholtzBCElementBase<ELEMENT>*>(this->element_pt(e));
1664 
1665  // Compute contribution
1666  el_pt->compute_contribution_to_fourier_components(el_a_coeff_pos,
1667  el_a_coeff_neg);
1668 
1669  // Add to coefficients
1670  for (unsigned i=0;i<n;i++)
1671  {
1672  a_coeff_pos[i]+=el_a_coeff_pos[i];
1673  a_coeff_neg[i]+=el_a_coeff_neg[i];
1674  }
1675  }
1676 
1677 }
1678 
1679 
1680 
1681 
1682 ///================================================================
1683 /// Compute and store the gamma integral and derivates
1684 // /w.r.t global unknows at all integration points
1685 /// of the mesh's constituent elements
1686 //================================================================
1687 template<class ELEMENT>
1689 {
1690 
1691 #ifdef PARANOID
1692  {
1693  // Loop over elements e
1694  unsigned nel=this->nelement();
1695  for (unsigned e=0;e<nel;e++)
1696  {
1697  FiniteElement* fe_pt=finite_element_pt(e);
1698  unsigned nnod=fe_pt->nnode();
1699  for (unsigned j=0;j<nnod;j++)
1700  {
1701  Node* nod_pt=fe_pt->node_pt(j);
1702 
1703  // Extract nodal coordinates from node:
1704  Vector<double> x(2);
1705  x[0]=nod_pt->x(0);
1706  x[1]=nod_pt->x(1);
1707 
1708  // Evaluate the radial distance
1709  double r=sqrt(x[0]*x[0]+x[1]*x[1]);
1710 
1711  // Check
1712  if (Outer_radius==0.0)
1713  {
1714  throw OomphLibError(
1715  "Outer radius for DtN BC must not be zero!",
1716  OOMPH_CURRENT_FUNCTION,
1717  OOMPH_EXCEPTION_LOCATION);
1718  }
1719 
1720  if(std::fabs((r-this->Outer_radius)/Outer_radius)
1722  {
1723  std::ostringstream error_stream;
1724  error_stream << "Node at " << x[0] << " " << x[1]
1725  << " has radius " << r << " which does not "
1726  << " agree with \nspecified outer radius "
1727  << this->Outer_radius << " within relative tolerance "
1729  << ".\nYou can adjust the tolerance via\n"
1730  << "ToleranceForHelmholtzOuterBoundary::Tol\n"
1731  << "or recompile without PARANOID.\n";
1732  throw OomphLibError(
1733  error_stream.str(),
1734  OOMPH_CURRENT_FUNCTION,
1735  OOMPH_EXCEPTION_LOCATION);
1736  }
1737  }
1738  }
1739  }
1740 #endif
1741 
1742 
1743  // Get Helmholtz parameter from bulk element
1746  (this->element_pt(0));
1747  double k=sqrt(dynamic_cast<ELEMENT*>(el_pt->bulk_element_pt())->k_squared());
1748 
1749  // Precompute factors in sum
1750  Vector<std::complex<double> > h_a(Nfourier_terms), hp_a(Nfourier_terms),
1751  q(Nfourier_terms);
1753  Hankel_first(Nfourier_terms,Outer_radius*k,h_a,hp_a);
1754  for (unsigned i=0;i<Nfourier_terms;i++)
1755  {
1756  q[i]=(k/(2.0*MathematicalConstants::Pi))*(hp_a[i]/h_a[i]);
1757  }
1758 
1759  //first loop over elements e
1760  unsigned nel=this->nelement();
1761  for (unsigned e=0;e<nel;e++)
1762  {
1763  // Get a pointer to element
1766  (this->element_pt(e));
1767 
1768  //Set the value of n_intpt
1769  const unsigned n_intpt =el_pt->integral_pt()->nweight();
1770 
1771  // initialise gamma integral and its derivatives
1772  Vector<std::complex<double> > gamma_vector(
1773  n_intpt,std::complex<double>(0.0,0.0));
1775  d_gamma_vector(n_intpt);
1776 
1777  //Loop over the integration points
1778  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1779  {
1780  //Allocate and initialise coordinate
1781  unsigned ndim_local=el_pt->dim();
1782  Vector<double> x(ndim_local+1,0.0);
1783 
1784  //Set the Vector to hold local coordinates
1785  Vector<double> s(ndim_local,0.0);
1786  for(unsigned i=0;i<ndim_local;i++)
1787  {
1788  s[i]=el_pt->integral_pt()->knot(ipt,i);
1789  }
1790 
1791  //Get the coordinates of the integration point
1792  el_pt->interpolated_x(s,x);
1793 
1794  // Polar angle
1795  double phi=atan2(x[1],x[0]);
1796 
1797  // Elemental contribution to gamma integral and its derivative
1798  std::complex<double> gamma_con_p(0.0,0.0),gamma_con_n(0.0,0.0);
1799  std::map<unsigned,std::complex<double> > d_gamma_con_p,d_gamma_con_n;
1800 
1801  // loop over the Fourier terms
1802  for (unsigned nn=0;nn<Nfourier_terms;nn++)
1803  {
1804  //Second loop over the element
1805  //to evaluate the complete integral
1806  for (unsigned ee=0;ee<nel;ee++)
1807  {
1810  (this->element_pt(ee));
1811 
1812  // contribution of the positive term in the sum
1814  phi,int(nn),gamma_con_p,d_gamma_con_p) ;
1815 
1816  // contribution of the negative term in the sum
1818  phi,-int(nn),gamma_con_n,d_gamma_con_n) ;
1819 
1820  unsigned n_node=eel_pt->nnode();
1821  if (nn==0)
1822  {
1823  gamma_vector[ipt]+=q[nn]*gamma_con_p;
1824  for(unsigned l=0;l<n_node;l++)
1825  {
1826  // Add the contribution of the real local data
1827  int local_unknown_p_real=eel_pt->nodal_local_eqn(
1828  l,eel_pt->u_index_helmholtz().real());
1829  if (local_unknown_p_real >= 0)
1830  {
1831  int global_unknown_p_real=eel_pt->eqn_number(
1832  local_unknown_p_real);
1833  d_gamma_vector[ipt][global_unknown_p_real]+=
1834  q[nn]*d_gamma_con_p[global_unknown_p_real];
1835  }
1836 
1837  // Add the contribution of the imag local data
1838  int local_unknown_p_imag=
1839  eel_pt->nodal_local_eqn(
1840  l,eel_pt->u_index_helmholtz().imag());
1841 
1842  if (local_unknown_p_imag >= 0)
1843  {
1844  int global_unknown_p_imag=eel_pt->eqn_number(
1845  local_unknown_p_imag);
1846 
1847  d_gamma_vector[ipt][global_unknown_p_imag]+=
1848  q[nn]*d_gamma_con_p[global_unknown_p_imag];
1849  }
1850  }// end of loop over the node
1851  }//End of if
1852  else
1853  {
1854  gamma_vector[ipt]+=q[nn]*(gamma_con_p+gamma_con_n);
1855  for(unsigned l=0;l<n_node;l++)
1856  {
1857  // Add the contribution of the real local data
1858  int local_unknown_real=eel_pt->nodal_local_eqn(
1859  l,eel_pt->u_index_helmholtz().real());
1860  if (local_unknown_real >= 0)
1861  {
1862  int global_unknown_real=eel_pt->eqn_number(local_unknown_real);
1863  d_gamma_vector[ipt][global_unknown_real]+=
1864  q[nn]*(d_gamma_con_p[global_unknown_real]+
1865  d_gamma_con_n[global_unknown_real]);
1866  }
1867  // Add the contribution of the imag local data
1868  int local_unknown_imag=eel_pt->nodal_local_eqn(
1869  l,eel_pt->u_index_helmholtz().imag());
1870  if (local_unknown_imag >= 0)
1871  {
1872  int global_unknown_imag=eel_pt->eqn_number(local_unknown_imag);
1873  d_gamma_vector[ipt][global_unknown_imag]+=
1874  q[nn]*(d_gamma_con_p[global_unknown_imag]+
1875  d_gamma_con_n[global_unknown_imag]);
1876  }
1877  }// end of loop over the node
1878  }//End of else
1879  }// End of second loop over the elements
1880  }// End of loop over Fourier terms
1881  }// end of loop over integration point
1882 
1883  // Store it in map
1884  Gamma_at_gauss_point[el_pt]=gamma_vector;
1885  D_Gamma_at_gauss_point[el_pt]=d_gamma_vector;
1886 
1887  }// end of first loop over element
1888 }
1889 
1890 //===========================================================================
1891 /// Constructor, takes the pointer to the "bulk" element and the face index
1892 //===========================================================================
1893 template<class ELEMENT>
1896  const int &face_index) :
1897  FaceGeometry<ELEMENT>(), FaceElement()
1898  {
1899 #ifdef PARANOID
1900  {
1901  //Check that the element is not a refineable 3d element
1902  ELEMENT* elem_pt = new ELEMENT;
1903  //If it's three-d
1904  if(elem_pt->dim()==3)
1905  {
1906  //Is it refineable
1907  if(dynamic_cast<RefineableElement*>(elem_pt))
1908  {
1909  //Issue a warning
1911  "This flux element will not"
1912  "work correctly if nodes are hanging\n",
1913  "HelmholtzBCElementBase::Constructor",
1914  OOMPH_EXCEPTION_LOCATION);
1915  }
1916  }
1917  }
1918 #endif
1919 
1920  // Let the bulk element build the FaceElement, i.e. setup the pointers
1921  // to its nodes (by referring to the appropriate nodes in the bulk
1922  // element), etc.
1923  bulk_el_pt->build_face_element(face_index,this);
1924 
1925  // Extract the dimension of the problem from the dimension of
1926  // the first node
1927  Dim = this->node_pt(0)->ndim();
1928 
1929  //Set up U_index_helmholtz. Initialise to zero, which probably won't change
1930  //in most cases, oh well, the price we pay for generality
1931  U_index_helmholtz = std::complex<unsigned>(0,1);
1932 
1933  //Cast to the appropriate HelmholtzEquation so that we can
1934  //find the index at which the variable is stored
1935  //We assume that the dimension of the full problem is the same
1936  //as the dimension of the node, if this is not the case you will have
1937  //to write custom elements, sorry
1938  switch(Dim)
1939  {
1940  //One dimensional problem
1941  case 1:
1942  {
1943  HelmholtzEquations<1>* eqn_pt =
1944  dynamic_cast<HelmholtzEquations<1>*>(bulk_el_pt);
1945  //If the cast has failed die
1946  if(eqn_pt==0)
1947  {
1948  std::string error_string =
1949  "Bulk element must inherit from HelmholtzEquations.";
1950  error_string +=
1951  "Nodes are one dimensional, but cannot cast the bulk element to\n";
1952  error_string += "HelmholtzEquations<1>\n.";
1953  error_string +=
1954  "If you desire this functionality, you must implement it yourself\n";
1955 
1956  throw OomphLibError(error_string,
1957  OOMPH_CURRENT_FUNCTION,
1958  OOMPH_EXCEPTION_LOCATION);
1959  }
1960  //Otherwise read out the value
1961  else
1962  {
1963  //Read the index from the (cast) bulk element
1965  }
1966  }
1967  break;
1968 
1969  //Two dimensional problem
1970  case 2:
1971  {
1972  HelmholtzEquations<2>* eqn_pt =
1973  dynamic_cast<HelmholtzEquations<2>*>(bulk_el_pt);
1974  //If the cast has failed die
1975  if(eqn_pt==0)
1976  {
1977  std::string error_string =
1978  "Bulk element must inherit from HelmholtzEquations.";
1979  error_string +=
1980  "Nodes are two dimensional, but cannot cast the bulk element to\n";
1981  error_string += "HelmholtzEquations<2>\n.";
1982  error_string +=
1983  "If you desire this functionality, you must implement it yourself\n";
1984 
1985  throw OomphLibError(error_string,
1986  OOMPH_CURRENT_FUNCTION,
1987  OOMPH_EXCEPTION_LOCATION);
1988  }
1989  else
1990  {
1991  //Read the index from the (cast) bulk element
1993  }
1994  }
1995 
1996  break;
1997 
1998  //Three dimensional problem
1999  case 3:
2000  {
2001  HelmholtzEquations<3>* eqn_pt =
2002  dynamic_cast<HelmholtzEquations<3>*>(bulk_el_pt);
2003  //If the cast has failed die
2004  if(eqn_pt==0)
2005  {
2006  std::string error_string =
2007  "Bulk element must inherit from HelmholtzEquations.";
2008  error_string +=
2009  "Nodes are three dimensional, but cannot cast the bulk element to\n";
2010  error_string += "HelmholtzEquations<3>\n.";
2011  error_string +=
2012  "If you desire this functionality, you must implement it yourself\n";
2013 
2014  throw OomphLibError(error_string,
2015  OOMPH_CURRENT_FUNCTION,
2016  OOMPH_EXCEPTION_LOCATION);
2017  }
2018  else
2019  {
2020  //Read the index from the (cast) bulk element
2022  }
2023  }
2024  break;
2025 
2026  //Any other case is an error
2027  default:
2028  std::ostringstream error_stream;
2029  error_stream << "Dimension of node is " << Dim
2030  << ". It should be 1,2, or 3!" << std::endl;
2031 
2032  throw OomphLibError(error_stream.str(),
2033  OOMPH_CURRENT_FUNCTION,
2034  OOMPH_EXCEPTION_LOCATION);
2035  break;
2036  }
2037  }
2038 
2039 
2040 }
2041 
2042 #endif
2043 
HelmholtzBCElementBase()
Broken empty constructor.
double global_power_contribution()
Compute the element's contribution to the time-averaged radiated power over the artificial boundary...
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
unsigned * ABC_order_pt
Pointer to order of absorbing boundary condition.
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
double & outer_radius()
The outer radius.
HelmholtzDtNMesh(const double &outer_radius, const unsigned &nfourier_terms)
void Hankel_first(const unsigned &n, const double &x, Vector< std::complex< double > > &h, Vector< std::complex< double > > &hp)
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 complete_setup_of_dependencies()
Complete the setup of additional dependencies arising through the far-away interaction with other nod...
const double Pi
50 digits from maple
virtual Node * copied_node_pt() const
Return pointer to copied node (null if the current node is not a copy – always the case here; it's ov...
Definition: nodes.h:1042
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
virtual std::complex< unsigned > u_index_helmholtz() const
Broken assignment operator.
A general Finite Element class.
Definition: elements.h:1271
HelmholtzDtNMesh< ELEMENT > * Outer_boundary_mesh_pt
Pointer to mesh of all DtN boundary condition elements (needed to get access to gamma values) ...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void compute_contribution_to_fourier_components(Vector< std::complex< double > > &a_coeff_pos, Vector< std::complex< double > > &a_coeff_neg)
Compute element's contribution to Fourier components of the solution – length of vector indicates num...
void CHankel_first(const unsigned &n, const std::complex< double > &x, Vector< std::complex< double > > &h, Vector< std::complex< double > > &hp)
void output(std::ostream &outfile)
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
e
Definition: cfortran.h:575
void fill_in_generic_residual_contribution_helmholtz_DtN_bc(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute the element's residual vector Jacobian matrix. Overloaded version, using the gamma computed i...
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
double global_power_contribution(std::ofstream &outfile)
Compute the element's contribution to the time-averaged radiated power over the artificial boundary...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Definition: elements.cc:293
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:66
Vector< std::map< unsigned, std::complex< double > > > & d_gamma_at_gauss_point(FiniteElement *el_pt)
Derivative of Gamma integral w.r.t global unknown, evaluated at Gauss points for specified element...
std::complex< unsigned > U_index_helmholtz
The index at which the real and imag part of the unknown is stored at the nodes.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
HelmholtzDtNBoundaryElement(FiniteElement *const &bulk_el_pt, const int &face_index)
Construct element from specification of bulk element and face index.
virtual bool is_a_copy() const
Return a boolean to indicate whether the Data objact contains any copied values. A base Data object c...
Definition: nodes.h:255
A class for elements that allow the approximation of the Sommerfeld radiation BC. The element geometr...
double d_shape_and_test_local(const Vector< double > &s, Shape &psi, Shape &test, DShape &dpsi_ds, DShape &dtest_ds) const
Function to compute the shape, test functions and derivates and to return the Jacobian of mapping bet...
const std::complex< double > I(0.0, 1.0)
The imaginary unit.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
static char t char * s
Definition: cfortran.h:572
double test_only(const Vector< double > &s, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
std::map< FiniteElement *, Vector< std::complex< double > > > Gamma_at_gauss_point
Container to store the gamma integral for given Gauss point and element.
double * Outer_radius_pt
Pointer to radius of outer boundary (must be a circle!)
void setup_gamma()
Compute and store the gamma integral at all integration points of the constituent elements...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
double Outer_radius
Outer radius.
void compute_fourier_components(Vector< std::complex< double > > &a_coeff_pos, Vector< std::complex< double > > &a_coeff_neg)
Compute Fourier components of the solution – length of vector indicates number of terms to be compute...
void fill_in_generic_residual_contribution_helmholtz_abc(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute the element's residual vector and the ( Jacobian matrix. Overloaded version, using the abc approximation.
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
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Broken assignment operator.
virtual std::complex< unsigned > u_index_helmholtz() const
Return the index at which the real/imag unknown value is stored.
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
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.
unsigned & nfourier_terms()
Number of Fourier terms used in the computation of the gamma integral.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary, when viewed as part of a compound geometric object is specified using the boundary coordinate defined by the mesh. Note: Boundary coordinates will have been set up when creating the underlying mesh, and their values will have been stored at the nodes.
Definition: elements.h:4251
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2215
HelmholtzBCElementBase(const HelmholtzBCElementBase &dummy)
Broken copy constructor.
HelmholtzAbsorbingBCElement(FiniteElement *const &bulk_el_pt, const int &face_index)
Construct element from specification of bulk element and face index.
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
unsigned *& abc_order_pt()
Pointer to order of absorbing boundary condition.
std::map< FiniteElement *, Vector< std::map< unsigned, std::complex< double > > > > D_Gamma_at_gauss_point
Container to store the derivate of Gamma integral w.r.t global unknown evaluated at Gauss points for ...
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
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 Nfourier_terms
Nbr of Fourier terms used in the Gamma computation.
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
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
unsigned nexternal_data() const
Return the number of external data objects.
Definition: elements.h:831
double Tol
Relative tolerance to within radius of points on DtN boundary are allowed to deviate from specified v...
Vector< std::complex< double > > & gamma_at_gauss_point(FiniteElement *el_pt)
Gamma integral evaluated at Gauss points for specified element.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and its Jacobian matrix.
void compute_gamma_contribution(const double &phi, const int &n, std::complex< double > &gamma_con, std::map< unsigned, std::complex< double > > &d_gamma_con)
Compute the contribution of the element to the Gamma integral and its derivates w.r.t to global unknows; the function takes the wavenumber and the polar angle phi as input.
unsigned Dim
The spatial dimension of the problem.
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...
HelmholtzDtNMesh< ELEMENT > * outer_boundary_mesh_pt() const
Access function to mesh of all DtN boundary condition elements (needed to get access to gamma values)...
A general mesh class.
Definition: mesh.h:74
int external_local_eqn(const unsigned &i, const unsigned &j)
Return the local equation number corresponding to the j-th value stored at the i-th external data...
Definition: elements.h:313
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
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
Definition: elements.h:1914
double *& outer_radius_pt()
Get pointer to radius of outer boundary (must be a cirle)
void set_outer_boundary_mesh_pt(HelmholtzDtNMesh< ELEMENT > *mesh_pt)
Set mesh of all DtN boundary condition elements.