fourier_decomposed_helmholtz_bc_elements.h

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
// Header file for elements that are used to apply Sommerfeld
// boundary conditions to the Fourier decomposed Helmholtz equations
#ifndef OOMPH_FOURIER_DECOMPOSED_HELMHOLTZ_BC_ELEMENTS_HEADER
#define OOMPH_FOURIER_DECOMPOSED_HELMHOLTZ_BC_ELEMENTS_HEADER


// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif

#include "math.h"
#include<complex>

// Get the Bessel functions
#include "oomph_crbond_bessel.h"

#include "fourier_decomposed_helmholtz_elements.h"

namespace oomph
{




/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////


//======================================================================
/// \short A class for elements that allow the approximation of the 
/// Sommerfeld radiation BC for Fourier decomposed Helmholtz equations.
/// The element geometry is obtained from the  FaceGeometry<ELEMENT> 
/// policy class.
//======================================================================
template <class ELEMENT>
 class FourierDecomposedHelmholtzBCElementBase : 
public virtual FaceGeometry<ELEMENT>, 
 public virtual FaceElement
 {
   public:
  
  /// \short Constructor, takes the pointer to the "bulk" element and the 
  /// index of the face to which the element is attached.
  FourierDecomposedHelmholtzBCElementBase(FiniteElement* const &bulk_el_pt, 
                         const int& face_index); 
  
  ///\short Broken empty constructor
  FourierDecomposedHelmholtzBCElementBase()
   {
    throw OomphLibError(
     "Don't call empty constructor for FourierDecomposedHelmholtzBCElementBase",
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }
  
  /// Broken copy constructor
  FourierDecomposedHelmholtzBCElementBase(const FourierDecomposedHelmholtzBCElementBase& dummy) 
   { 
    BrokenCopy::broken_copy("FourierDecomposedHelmholtzBCElementBase");
   } 
  
  /// Broken assignment operator
//Commented out broken assignment operator because this can lead to a conflict warning
//when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
//realise that two separate implementations of the broken function are the same and so,
//quite rightly, it shouts.
  /*void operator=(const FourierDecomposedHelmholtzBCElementBase&) 
   {
    BrokenCopy::broken_assign("FourierDecomposedHelmholtzBCElementBase");
    }*/
  

  /// \short Specify the value of nodal zeta from the face geometry
  /// The "global" intrinsic coordinate of the element when
  /// viewed as part of a geometric object should be given by
  /// the FaceElement representation, by default (needed to break
  /// indeterminacy if bulk element is SolidElement)
  double zeta_nodal(const unsigned &n, const unsigned &k,           
                    const unsigned &i) const 
  {return FaceElement::zeta_nodal(n,k,i);}     
  

  /// Output function -- forward to broken version in FiniteElement
  /// until somebody decides what exactly they want to plot here...
  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
  
  /// \short Output function -- forward to broken version in FiniteElement
  /// until somebody decides what exactly they want to plot here...
  void output(std::ostream &outfile, const unsigned &n_plot)
  {FiniteElement::output(outfile,n_plot);}
  
  /// C-style output function -- forward to broken version in FiniteElement
  /// until somebody decides what exactly they want to plot here...
  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
  
  /// \short C-style output function -- forward to broken version in 
  /// FiniteElement until somebody decides what exactly they want to plot 
  /// here...
  void output(FILE* file_pt, const unsigned &n_plot)
  {FiniteElement::output(file_pt,n_plot);}

  /// \short Return the index at which the real/imag unknown value
  /// is stored.
  virtual inline std::complex<unsigned> u_index_fourier_decomposed_helmholtz() 
   const 
   {
    return std::complex<unsigned>(U_index_fourier_decomposed_helmholtz.real(),
                                  U_index_fourier_decomposed_helmholtz.imag());
   }
  
  /// \short Compute the element's contribution to the time-averaged 
  /// radiated power over the artificial boundary
  double global_power_contribution()
  {   
   // Dummy output file
   std::ofstream outfile;
   return global_power_contribution(outfile);
  }
  
  /// \short Compute the element's contribution to the time-averaged 
  /// radiated power over the artificial boundary. Also output the
  /// power density as a fct of the zenith angle in the specified 
  ///output file if it's open.
  double global_power_contribution(std::ofstream& outfile)
  {   
   // pointer to the corresponding bulk element
   ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt()); 
   
   // Number of nodes in bulk element
   unsigned nnode_bulk=bulk_elem_pt->nnode();
   const unsigned n_node_local = this->nnode();

   //get the dim of the bulk and local nodes
   const unsigned bulk_dim= bulk_elem_pt->dim();    
   const unsigned local_dim= this->node_pt(0)->ndim();
   
   //Set up memory for the shape and test functions
   Shape psi(n_node_local);
   
   //Set up memory for the shape functions
   Shape psi_bulk(nnode_bulk);
   DShape dpsi_bulk_dx(nnode_bulk,bulk_dim);
   
   //Set up memory for the outer unit normal
   Vector< double > unit_normal(bulk_dim);    
   
   //Set the value of n_intpt
   const unsigned n_intpt = integral_pt()->nweight();
   
   //Set the Vector to hold local coordinates
   Vector<double> s(local_dim-1);
   double power=0.0;    
   
   // Output?
   if (outfile.is_open())
    {
     outfile << "ZONE\n";
    }
   
   //Loop over the integration points
   //--------------------------------
   for(unsigned ipt=0;ipt<n_intpt;ipt++)
    { 
     //Assign values of s
     for(unsigned i=0;i<(local_dim-1);i++)
      {
       s[i] = integral_pt()->knot(ipt,i);
      }
     //get the outer_unit_ext vector      
     this->outer_unit_normal(s,unit_normal); 
     
     //Get the integral weight
     double w = integral_pt()->weight(ipt);
     
     // Get jacobian of mapping
     double J=J_eulerian(s);
     
     //Premultiply the weights and the Jacobian
     double W = w*J;
     
     // Get local coordinates in bulk element by copy construction
     Vector<double> s_bulk(local_coordinate_in_bulk(s));
     
     //Call the derivatives of the shape  functions
     //in the bulk -- must do this via s because this point
     //is not an integration point the bulk element!
     (void)bulk_elem_pt->dshape_eulerian(s_bulk,psi_bulk,dpsi_bulk_dx);
     this->shape(s,psi);
     
     // Derivs of Eulerian coordinates w.r.t. local coordinates
     std::complex<double>  dphi_dn(0.0,0.0);
     Vector<std::complex <double> > interpolated_dphidx(bulk_dim);
     std::complex<double> interpolated_phi(0.0,0.0);
     Vector<double> x(bulk_dim);
     
     //Calculate function value and derivatives:
     //-----------------------------------------
     // Loop over nodes
     for(unsigned l=0;l<nnode_bulk;l++) 
      {
       //Get the nodal value of the helmholtz unknown
       const std::complex<double> phi_value(
        bulk_elem_pt->nodal_value(
         l,bulk_elem_pt->u_index_fourier_decomposed_helmholtz().real()),
        bulk_elem_pt->nodal_value(
         l,bulk_elem_pt->u_index_fourier_decomposed_helmholtz().imag())
        );
       
       //Loop over directions
       for(unsigned i=0;i<bulk_dim;i++)
        {
         interpolated_dphidx[i] += phi_value*dpsi_bulk_dx(l,i);
        }
      } // End of loop over the bulk_nodes
     
     for(unsigned l=0;l<n_node_local;l++) 
      {
       //Get the nodal value of the Helmholtz unknown
       const std::complex<double> phi_value(
        nodal_value(l,u_index_fourier_decomposed_helmholtz().real()),
        nodal_value(l,u_index_fourier_decomposed_helmholtz().imag()));
       
       interpolated_phi += phi_value*psi(l);
      }
     
     //define dphi_dn 
     for(unsigned i=0;i<bulk_dim;i++)
      {
       dphi_dn += interpolated_dphidx[i]*unit_normal[i];
      }
     
     // Power density
     double integrand=
      (interpolated_phi.real()*dphi_dn.imag()-
       interpolated_phi.imag()*dphi_dn.real());
     
     interpolated_x(s,x);
     double theta=atan2(x[0],x[1]);
     // Output?
     if (outfile.is_open())
      {
       outfile << x[0] << " "
               << x[1] << " "
               << theta << " "
               << integrand << "\n";
      }
     
     // ...add to integral
     power+=MathematicalConstants::Pi*x[0]*integrand*W;
    }  

   return  power;
  }
  
   protected:
  
  /// \short Function to compute the test functions and to return 
  /// the Jacobian of mapping between local and global (Eulerian)
  /// coordinates
  inline double shape_and_test(const Vector<double> &s, 
                               Shape& psi, Shape &test) const
  {
   //Get the shape functions
   shape(s,test);
   
   unsigned nnod=nnode();
   for (unsigned i=0;i<nnod;i++)
    {
     psi[i]=test[i];
    }
   
   //Return the value of the jacobian
   return J_eulerian(s);
  }
  
  /// \short Function to compute the shape, test functions and derivates 
  /// and to return
  /// the Jacobian of mapping between local and global (Eulerian)
  /// coordinates
  inline double d_shape_and_test_local(const Vector<double> &s, Shape &psi, 
                                       Shape &test,
                                       DShape &dpsi_ds,DShape &dtest_ds)
   const
  {
   //Find number of nodes
   unsigned n_node = nnode();
   
   //Get the shape functions
   dshape_local(s,psi,dpsi_ds);
   
   //Set the test functions to be the same as the shape functions
   for(unsigned i=0;i<n_node;i++) 
    {
     for(unsigned j=0;j<1;j++) 
      {      
       test[i] = psi[i];
       dtest_ds(i,j)= dpsi_ds(i,j);
      }
    }
   //Return the value of the jacobian
   return J_eulerian(s);
  }
  
  /// \short The index at which the real and imag part of the unknown is stored 
  /// at the nodes
  std::complex<unsigned> U_index_fourier_decomposed_helmholtz;
    
 }; 


//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////


///=================================================================
/// Mesh for DtN boundary condition elements -- provides
/// functionality to apply Sommerfeld radiation condtion
/// via DtN BC.
///=================================================================
template<class ELEMENT>
 class FourierDecomposedHelmholtzDtNMesh : public virtual Mesh
{
  public:
 
 /// Constructor: Specify radius of outer boundary and number of
 /// terms used in the computation of the gamma integral
  FourierDecomposedHelmholtzDtNMesh(const double& outer_radius, 
                                    const unsigned& n_terms) : 
 Outer_radius(outer_radius), N_terms(n_terms)
 {}
 
 /// \short Compute and store the gamma integral at all integration
 /// points of the constituent elements.
 void setup_gamma();
  
 /// \short Gamma integral evaluated at Gauss points 
 /// for specified element
 Vector<std::complex<double> >& gamma_at_gauss_point(FiniteElement* el_pt) 
  {                                          
   return Gamma_at_gauss_point[el_pt];
  }
 
 /// \short Derivative of Gamma integral w.r.t global unknown, evaluated 
 /// at Gauss points for specified element
 Vector<std::map<unsigned,std::complex<double> > >
  &d_gamma_at_gauss_point(FiniteElement* el_pt) 
  {                                          
   return D_Gamma_at_gauss_point[el_pt];
  }
 
 /// \short The outer radius  
 double &outer_radius() 
 {                                          
  return Outer_radius ;
 }

 /// \short Number of terms used in the computation of the
 /// gamma integral 
 unsigned& n_terms() 
 {                                          
  return N_terms;
 }
 
  private:
 
 /// Outer radius
 double Outer_radius;
 
 /// Number of terms used in the  Gamma computation
 unsigned N_terms;
 
 
 /// \short Container to store the gamma integral for given Gauss point
 /// and element
 std::map<FiniteElement*,Vector<std::complex<double> > > Gamma_at_gauss_point;
 
 
 /// \short Container to store the derivate of Gamma integral w.r.t 
 /// global unknown evaluated at Gauss points for specified element
 std::map<FiniteElement*,Vector<std::map<unsigned,std::complex<double> > > > 
  D_Gamma_at_gauss_point;
 
};

/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////

//=============================================================
/// FaceElement used to apply Sommerfeld radiation conditon
/// via Dirichlet to Neumann map. 
//==============================================================  
template<class ELEMENT>
class FourierDecomposedHelmholtzDtNBoundaryElement : 
public FourierDecomposedHelmholtzBCElementBase<ELEMENT>
{
 
  public:
 
  /// \short Construct element from specification of bulk element and
  /// face index.
  FourierDecomposedHelmholtzDtNBoundaryElement(
   FiniteElement* const &bulk_el_pt, 
   const int& face_index) : 
 FourierDecomposedHelmholtzBCElementBase<ELEMENT>(bulk_el_pt,face_index)
  {}
 
 /// Add the element's contribution to its residual vector
 inline void fill_in_contribution_to_residuals(Vector<double> &residuals)
 {
  //Call the generic residuals function with flag set to 0
  //using a dummy matrix argument
  fill_in_generic_residual_contribution_fourier_decomposed_helmholtz_DtN_bc
   (residuals,GeneralisedElement::Dummy_matrix,0);
 }
 

 /// \short Add the element's contribution to its residual vector and its
 /// Jacobian matrix
 inline void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                              DenseMatrix<double> &jacobian)
 {
  //Call the generic routine with the flag set to 1
  fill_in_generic_residual_contribution_fourier_decomposed_helmholtz_DtN_bc
   (residuals,jacobian,1);
 }
 
 /// \short Compute the contribution of the element 
 /// to the Gamma integral and its derivates w.r.t 
 /// to global unknows; the function takes the wavenumber (for gamma integral,
 /// not the one from the Fourier decomposition of the Helmholtz equations!)
 /// and the polar angle theta as input.
 void compute_gamma_contribution(const double& theta,const unsigned& n, 
                                 std::complex<double>& gamma_con,
                                 std::map<unsigned,std::complex<double> >& 
                                 d_gamma_con);
 
 
 /// \short Access function to mesh of all DtN boundary condition elements
 /// (needed to get access to gamma values)
 FourierDecomposedHelmholtzDtNMesh<ELEMENT>* outer_boundary_mesh_pt() const
  {
   return Outer_boundary_mesh_pt;
  }
 
 /// \short Set mesh of all DtN boundary condition elements
 void set_outer_boundary_mesh_pt
  (FourierDecomposedHelmholtzDtNMesh<ELEMENT>* mesh_pt)
 {
  Outer_boundary_mesh_pt=mesh_pt;
 }
 
 
 /// \short Complete the setup of additional dependencies arising
 /// through the far-away interaction with other nodes in 
 /// Outer_boundary_mesh_pt.
 void complete_setup_of_dependencies()
 {
  // Create a set of all nodes
  std::set<Node*> node_set;
  unsigned nel=Outer_boundary_mesh_pt->nelement();
  for (unsigned e=0;e<nel;e++)
   {
    FiniteElement* el_pt=Outer_boundary_mesh_pt->finite_element_pt(e);
    unsigned nnod=el_pt->nnode();
    for (unsigned j=0;j<nnod;j++)
     {
      Node* nod_pt=el_pt->node_pt(j);
      
      // Don't add copied nodes
      if (!(nod_pt->is_a_copy()))
       {
        node_set.insert(nod_pt);
       }
     }
   }
  // Now erase the current element's own nodes
  unsigned nnod=this->nnode();
  for (unsigned j=0;j<nnod;j++)
   {
    Node* nod_pt=this->node_pt(j);
    node_set.erase(nod_pt);
    
    // If the element's node is a copy then its "master" will 
    // already have been added in the set above -- remove the
    // master to avoid double counting eqn numbers
    if (nod_pt->is_a_copy())
     {
      node_set.erase(nod_pt->copied_node_pt());     
     }
   }
  
  // Now declare these nodes to be the element's external Data
  for (std::set<Node*>::iterator it=node_set.begin();
       it!=node_set.end();it++)
   {
    this->add_external_data(*it);
   }
 }


 
  private:
 
 /// \short Compute the element's residual vector  
 /// Jacobian matrix.
 /// Overloaded version, using the gamma computed in the mesh
 void fill_in_generic_residual_contribution_fourier_decomposed_helmholtz_DtN_bc
  (Vector<double> &residuals, DenseMatrix<double> &jacobian, 
   const unsigned& flag)
 {
  //Find out how many nodes there are
  const unsigned n_node = this->nnode();
  
  //Set up memory for the shape and test functions
  Shape test(n_node);
  Shape psi(n_node);
  
  //Set the value of Nintpt
  const unsigned n_intpt = this->integral_pt()->nweight();
  
  //Set the Vector to hold local coordinates
  Vector<double> s(1);
  
  //Integers to hold the local equation and unknown numbers
  int local_eqn_real=0,local_unknown_real=0,global_unk_real=0,
   local_eqn_imag=0,local_unknown_imag=0,global_unk_imag=0;
  int external_global_unk_real=0, external_unknown_real=0,
   external_global_unk_imag=0, external_unknown_imag=0;
  
  
  // Get the gamma value for the current integration point
  // from the mesh
  Vector<std::complex<double> > 
   gamma(Outer_boundary_mesh_pt->gamma_at_gauss_point(this));
  
  Vector<std::map<unsigned,std::complex<double> > >
   d_gamma(Outer_boundary_mesh_pt->d_gamma_at_gauss_point(this));
  
  //Loop over the integration points
  //--------------------------------
  for(unsigned ipt=0;ipt<n_intpt;ipt++)
   {
    //Assign values of s
    for(unsigned i=0;i<1;i++)
     {
      s[i] = this->integral_pt()->knot(ipt,i);
     }
    
    //Get the integral weight
    double w = this->integral_pt()->weight(ipt);
    
    //Find the shape test functions and derivates; return the Jacobian
    //of the mapping between local and global (Eulerian)
    // coordinates
    double J = this->shape_and_test(s,psi,test);
    
    //Premultiply the weights and the Jacobian
    double W = w*J;
    
    // Build up radius
    double r = 0.0;
    for (unsigned j=0;j<n_node;j++)
     {
      r+=this->node_pt(j)->x(0)*psi(j);
     }

    //Now add to the appropriate equations
    //Loop over the test functions:loop over the nodes
    for(unsigned l=0;l<n_node;l++)
     {
      local_eqn_real = this->nodal_local_eqn
       (l,this->U_index_fourier_decomposed_helmholtz.real());

      local_eqn_imag = this->nodal_local_eqn
       (l,this->U_index_fourier_decomposed_helmholtz.imag());
      
      //IF it's not a boundary condition
      if(local_eqn_real >= 0)
       {
        //Add the gamma contribution in this int_point to the res
        residuals[local_eqn_real] +=gamma[ipt].real()*test[l]*r*W; 
        
        // Calculate the jacobian
        //-----------------------
        if(flag)
         {
          //Loop over the shape functions again
          for(unsigned l2=0;l2<n_node;l2++)
           { 
            // Add the contribution of the local data
            local_unknown_real = this->nodal_local_eqn(
             l2,this->U_index_fourier_decomposed_helmholtz.real());
            
            local_unknown_imag = this->nodal_local_eqn(
              l2,this->U_index_fourier_decomposed_helmholtz.imag());
             
             //If at a non-zero degree of freedom add in the entry
             if(local_unknown_real >= 0)
              {
               // Add the contribution
               global_unk_real=this->eqn_number(local_unknown_real);
               jacobian(local_eqn_real,local_unknown_real)
                +=d_gamma[ipt][global_unk_real].real()*test[l]*r*W; 
              }
             if(local_unknown_imag >= 0)
              {
               // Add the contribution
               global_unk_imag=this->eqn_number(local_unknown_imag);
               jacobian(local_eqn_real,local_unknown_imag)
                +=d_gamma[ipt][global_unk_imag].real()*test[l]*r*W; 
              }
            } // End of loop over nodes l2  
           
           // Add the contribution of the external data
           unsigned n_ext_data=this->nexternal_data();
           //Loop over the shape functions again
           for(unsigned l2=0;l2<n_ext_data;l2++)
            { 
             // Add the contribution of the local data
             external_unknown_real = this->external_local_eqn(
              l2,this->U_index_fourier_decomposed_helmholtz.real());
             
             external_unknown_imag = this->external_local_eqn(
              l2,this->U_index_fourier_decomposed_helmholtz.imag());
             
             //If at a non-zero degree of freedom add in the entry
             if(external_unknown_real >= 0)
              {
               // Add the contribution
               external_global_unk_real=this->eqn_number(external_unknown_real);
               jacobian(local_eqn_real,external_unknown_real)
                +=d_gamma[ipt][external_global_unk_real].real()*test[l]*r*W;
              }
             if(external_unknown_imag >= 0)
              {
               // Add the contribution
               external_global_unk_imag=this->eqn_number(external_unknown_imag);
               jacobian(local_eqn_real,external_unknown_imag)
                +=d_gamma[ipt][external_global_unk_imag].real()*test[l]*r*W;
              }
            } // End of loop over external data    
          }// End of flag
        }// end of local_eqn_real
       
       if(local_eqn_imag >= 0)
        {
         //Add the gamma contribution in this int_point to the res
         residuals[local_eqn_imag] +=gamma[ipt].imag()*test[l]*r*W; 
         
         // Calculate the jacobian
         //-----------------------
         if(flag)
          {
           //Loop over the shape functions again
           for(unsigned l2=0;l2<n_node;l2++)
            { 
             // Add the contribution of the local data
             local_unknown_real = this->nodal_local_eqn(
              l2,this->U_index_fourier_decomposed_helmholtz.real());
             
             local_unknown_imag = this->nodal_local_eqn(
              l2,this->U_index_fourier_decomposed_helmholtz.imag());

             //If at a non-zero degree of freedom add in the entry
             if(local_unknown_real >= 0)
              {
               // Add the contribution
               global_unk_real=this->eqn_number(local_unknown_real);
               jacobian(local_eqn_imag,local_unknown_real)
                +=d_gamma[ipt][global_unk_real].imag()*test[l]*r*W;
              }
             if(local_unknown_imag >= 0)
              {
               // Add the  contribution
               global_unk_imag=this->eqn_number(local_unknown_imag);
               jacobian(local_eqn_imag,local_unknown_imag)
                +=d_gamma[ipt][global_unk_imag].imag()*test[l]*r*W; 
              }
            } // End of loop over nodes l2  
           
           // Add the contribution of the external data
           unsigned n_ext_data=this->nexternal_data();

           //Loop over the shape functions again
           for(unsigned l2=0;l2<n_ext_data;l2++)
            { 
             // Add the contribution of the local data
             external_unknown_real = this->external_local_eqn(
              l2,this->U_index_fourier_decomposed_helmholtz.real());
             
             external_unknown_imag = this->external_local_eqn(
              l2,this->U_index_fourier_decomposed_helmholtz.imag());
             
             //If at a non-zero degree of freedom add in the entry
             if(external_unknown_real >= 0)
              {
               // Add the contribution
               external_global_unk_real=this->eqn_number(external_unknown_real);
               jacobian(local_eqn_imag,external_unknown_real)
                +=d_gamma[ipt][external_global_unk_real].imag()*test[l]*r*W;
              }
             if(external_unknown_imag >= 0)
              {
               // Add the contribution
               external_global_unk_imag=
                this->eqn_number(external_unknown_imag); 
               jacobian(local_eqn_imag,external_unknown_imag)
                +=d_gamma[ipt][external_global_unk_imag].imag()*test[l]*r*W;
              }
            } // End of loop over external data    
          }// End of flag
        } // end of local_eqn_imag   
      }// end of loop over the nodes
    } //End of loop over int_pt
  } 
  
   
  /// \short Pointer to mesh of all DtN boundary condition elements
  /// (needed to get access to gamma values)
  FourierDecomposedHelmholtzDtNMesh<ELEMENT>* Outer_boundary_mesh_pt;
  
};


////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////



//===========start_compute_gamma_contribution==================
/// \short compute the contribution of the element 
/// to the Gamma integral and its derivates w.r.t 
/// to global unknows; the function takes wavenumber n
/// (from the computation of the gamma integral, not the one 
/// from the Fourier decomposition of the Helmholtz equations!) 
/// and polar angle theta as input.  
//==============================================================  
template<class ELEMENT>
void FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>::
compute_gamma_contribution(
 const double& theta,
 const unsigned& n, 
 std::complex<double>& gamma_con,
 std::map<unsigned,std::complex<double> >& d_gamma_con)
{  
 //Parameters
 int n_fourier_helmholtz=
  dynamic_cast<ELEMENT*>(this->bulk_element_pt())->fourier_wavenumber();
 
 // define the imaginary number
 const std::complex<double> I(0.0,1.0);
  
 //Find out how many nodes there are
 const unsigned n_node = this->nnode();
 
 //Set up memory for the shape  functions
 Shape psi(n_node); 
 DShape dpsi(n_node,1);
 
 // initialise the variable
 int local_unknown_real=0, local_unknown_imag=0;
 int global_eqn_real=0,global_eqn_imag=0;
 
 //Set the value of n_intpt
 const unsigned n_intpt=this->integral_pt()->nweight();
 
 //Set the Vector to hold local coordinates
 Vector<double> s(1);
 
 // Initialise
 gamma_con=std::complex<double>(0.0,0.0);
 d_gamma_con.clear();
 
 //Loop over the integration points
 //--------------------------------
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Assign values of s
   for(unsigned i=0;i<1;i++) 
    {
     s[i]=this->integral_pt()->knot(ipt,i);
    }
   
   //Get the integral weight
   double w=this->integral_pt()->weight(ipt);
   
   // Get the shape functions
   this->dshape_local(s,psi,dpsi);
   
   // Eulerian coordinates at Gauss point
   Vector<double> interpolated_x(2,0.0);
   
   // Derivs of Eulerian coordinates w.r.t. local coordinates
   Vector<double> interpolated_dxds(2);
   std::complex<double> interpolated_u(0.0,0.0);
   
   // Assemble x and its derivs
   for(unsigned l=0;l<n_node;l++) 
    {
     //Loop over directions
     for(unsigned i=0;i<2;i++)
      {
       interpolated_x[i]+=this->nodal_position(l,i)*psi[l];
       interpolated_dxds[i]+=this->nodal_position(l,i)*dpsi(l,0);
      }
     
     //Get the nodal value of the helmholtz unknown
     std::complex<double> u_value(
      this->nodal_value(l,this->U_index_fourier_decomposed_helmholtz.real()),
      this->nodal_value(l,this->U_index_fourier_decomposed_helmholtz.imag()));
     
     interpolated_u += u_value*psi(l);
     
    } // End of loop over the nodes

   
   // calculate the integral
   //-----------------------
   // define the variable theta
   double phi=atan2(interpolated_x[0],interpolated_x[1]);
   
   //define dphi_ds=-z'/r
   double dphi_ds=-std::fabs(-interpolated_dxds[1]/interpolated_x[0]);
   
   //define the associated legendre functions
   double p_theta = Legendre_functions_helper::plgndr2
    (n,n_fourier_helmholtz,cos(theta));

   double p_phi = Legendre_functions_helper::plgndr2
    (n,n_fourier_helmholtz,cos(phi));
   
   gamma_con+=
    interpolated_u*p_phi*p_theta*sin(phi)*w*dphi_ds;
   
   // compute the contribution to each node to the map   
   for(unsigned l=0;l<n_node;l++) 
    {
     // Add the contribution of the real local data
     local_unknown_real = this->nodal_local_eqn(
      l,this->U_index_fourier_decomposed_helmholtz.real());
     if (local_unknown_real >= 0)
      {   
       global_eqn_real=this->eqn_number(local_unknown_real);
       d_gamma_con[global_eqn_real]+=
        psi(l)*p_phi*p_theta*sin(phi)*w*dphi_ds;
      }
     
     // Add the contribution of the imag local data
     local_unknown_imag = this->nodal_local_eqn(
      l,this->U_index_fourier_decomposed_helmholtz.imag());
     if (local_unknown_imag >= 0)
      {   
       global_eqn_imag=this->eqn_number(local_unknown_imag);
       d_gamma_con[global_eqn_imag]+=
        I*psi(l)*p_phi*p_theta*sin(phi)*w*dphi_ds;
      }
    }// end of loop over the node
  }//End of loop over integration points    
 
}


/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////


//===========================================================================
/// \short Namespace for checking radius of nodes on (assumed to be circular) 
/// DtN boundary
//===========================================================================
namespace ToleranceForFourierDecomposedHelmholtzOuterBoundary
{
 /// \short Relative tolerance to within radius of points on DtN boundary
 /// are allowed to deviate from specified value
 extern double Tol;
 
}


//===========================================================================
/// \short Namespace for checking radius of nodes on (assumed to be circular) 
/// DtN boundary
//===========================================================================
namespace ToleranceForFourierDecomposedHelmholtzOuterBoundary
{
 /// \short Relative tolerance to within radius of points on DtN boundary
 /// are allowed to deviate from specified value
 double Tol=1.0e-3;

}

/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////


///================================================================
/// Compute and store the gamma integral and derivates 
// /w.r.t global unknows at all integration points
/// of the mesh's constituent elements
//================================================================
template<class ELEMENT>
void FourierDecomposedHelmholtzDtNMesh<ELEMENT>::setup_gamma()
{ 
#ifdef PARANOID
 {
  // Loop over elements e
  unsigned nel=this->nelement();
  for (unsigned e=0;e<nel;e++)
   {
    FiniteElement* fe_pt=finite_element_pt(e);
    unsigned nnod=fe_pt->nnode();
    for (unsigned j=0;j<nnod;j++)
     {
      Node* nod_pt=fe_pt->node_pt(j);
      
      // Extract nodal coordinates from node:
      Vector<double> x(2);
      x[0]=nod_pt->x(0);
      x[1]=nod_pt->x(1);
      
      // Evaluate the radial distance 
      double r=sqrt(x[0]*x[0]+x[1]*x[1]); 
      
      // Check
      if (Outer_radius==0.0)
       {
        throw OomphLibError(
         "Outer radius for DtN BC must not be zero!",
         OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);
       }
      
      if(std::fabs((r-this->Outer_radius)/Outer_radius)
         >ToleranceForFourierDecomposedHelmholtzOuterBoundary::Tol)
       { 
        std::ostringstream error_stream; 
        error_stream 
         << "Node at " << x[0] << " " << x[1] 
         << " has radius " << r << " which does not "
         << " agree with \nspecified outer radius "
         << this->Outer_radius << " within relative tolerance " 
         << ToleranceForFourierDecomposedHelmholtzOuterBoundary::Tol 
         << ".\nYou can adjust the tolerance via\n"
         << "ToleranceForFourierDecomposedHelmholtzOuterBoundary::Tol\n" 
         << "or recompile without PARANOID.\n";
        throw OomphLibError(
         error_stream.str(),
         OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);
       }
     }
   }
 }
#endif

 
 // Get parameters from first element
 FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>* el_pt=
  dynamic_cast<FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>*>
  (this->element_pt(0));    
 double k=sqrt(dynamic_cast<ELEMENT*>(el_pt->bulk_element_pt())->k_squared());  
 int n_fourier_decomposed=
  dynamic_cast<ELEMENT*>(el_pt->bulk_element_pt())->fourier_wavenumber();
 double n_hankel_order_max=double(N_terms)+0.5;
 double n_hankel_order_tmp=0.0;
 
/// Imaginary unit 
 std::complex<double> I(0.0,1.0);
 
// Precompute factors in sum
 Vector<std::complex<double> > h_a(N_terms+1), 
  hp_a(N_terms+1), q(N_terms+1,std::complex<double>(0.0,0.0));
 Vector<double> jv(N_terms+1), yv(N_terms+1), 
  djv(N_terms+1), dyv(N_terms+1);
 
 // Get Bessel functions
 CRBond_Bessel::bessjyv(n_hankel_order_max,
                        Outer_radius*k,
                        n_hankel_order_tmp,
                        &jv[0],&yv[0],
                        &djv[0],&dyv[0]);

 // Assemble Hankel functions and their derivatives
 for(unsigned j=0; j<N_terms+1;j++) 
  {
   h_a[j]=jv[j]+I*yv[j];
   hp_a[j]=djv[j]+I*dyv[j];
  }
 
 // Precompute relevant terms in sum
 for (unsigned i=n_fourier_decomposed;i<N_terms;i++)
  {
   double n_1=Legendre_functions_helper::factorial(i-n_fourier_decomposed);
   double n_2=Legendre_functions_helper::factorial(i+n_fourier_decomposed);

   q[i]=k*sqrt(MathematicalConstants::Pi/(2.0*k*Outer_radius))*
    (hp_a[i]-h_a[i]/(2.0*k*Outer_radius))*
    (2.0*double(i)+1.0)/
    (2.0*sqrt(MathematicalConstants::Pi/(2.0*k*Outer_radius))*h_a[i])*
    (n_1/n_2);
   
  } 
 
 //first loop over elements e
 unsigned nel=this->nelement();
 for (unsigned e=0;e<nel;e++)
  {
   // Get a pointer to element   
   FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>* el_pt=
    dynamic_cast<FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>*>
    (this->element_pt(e));    
   
   //Set the value of n_intpt
   const unsigned n_intpt =el_pt->integral_pt()->nweight();
   
   // initialise gamma integral and its derivatives
   Vector<std::complex<double> > gamma_vector(
    n_intpt,std::complex<double>(0.0,0.0));
   Vector<std::map<unsigned,std::complex<double> > > 
    d_gamma_vector(n_intpt);
   
   //Loop over the integration points
   for(unsigned ipt=0;ipt<n_intpt;ipt++)
    {
     //Allocate and initialise coordinate
     Vector<double> x(el_pt->dim()+1,0.0);
     
     //Set the Vector to hold local coordinates
     Vector<double> s(el_pt->dim(),0.0);
     for(unsigned i=0;i<el_pt->dim();i++) 
      {
       s[i]=el_pt->integral_pt()->knot(ipt,i);
      }
     
     //Get the coordinates of the integration point
     el_pt->interpolated_x(s,x);
     
     // Polar angle
     double theta=atan2(x[0],x[1]);
     
     // Elemental contribution to gamma integral and its derivative
     std::complex<double> gamma_con(0.0,0.0);
     std::map<unsigned,std::complex<double> > d_gamma_con;
     
     // loop over the Fourier terms
     for (unsigned nn=n_fourier_decomposed;nn<N_terms;nn++)
      {
       //Second loop over the elements 
       //to evaluate the complete integral
       for (unsigned ee=0;ee<nel;ee++)
        {
         FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>* eel_pt=
          dynamic_cast<FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>*>
          (this->element_pt(ee));
         
         // contribution of the positive term in the sum
         eel_pt->compute_gamma_contribution(
          theta,nn,gamma_con,d_gamma_con) ;
         
         unsigned n_node=eel_pt->nnode();

         gamma_vector[ipt]+=q[nn]*gamma_con;
         for(unsigned l=0;l<n_node;l++) 
          {
           // Add the contribution of the real local data
           int local_unknown_real=
            eel_pt->nodal_local_eqn(
             l,eel_pt->u_index_fourier_decomposed_helmholtz().real());
           
           if (local_unknown_real >= 0)
            { 
             int global_eqn_real=eel_pt->eqn_number(
              local_unknown_real);  
             d_gamma_vector[ipt][global_eqn_real]+=
              q[nn]*d_gamma_con[global_eqn_real];
            }
           
           // Add the contribution of the imag local data
           int local_unknown_imag=
            eel_pt->nodal_local_eqn(
               l,eel_pt->u_index_fourier_decomposed_helmholtz().imag());
           
           if (local_unknown_imag >= 0)
            {
             int global_eqn_imag=eel_pt->eqn_number(
              local_unknown_imag);   
             d_gamma_vector[ipt][global_eqn_imag]+=
              q[nn]*d_gamma_con[global_eqn_imag];
            }
          }// end of loop over the node
        }// End of second loop over the elements 
      }// End of loop over Fourier terms  
    }// end of loop over integration point
   
   // Store it in map
   Gamma_at_gauss_point[el_pt]=gamma_vector;
   D_Gamma_at_gauss_point[el_pt]=d_gamma_vector;
   
  }// end of first loop over element  
}

//===========================================================================
/// Constructor, takes the pointer to the "bulk" element and the face index.
//===========================================================================
template<class ELEMENT>
FourierDecomposedHelmholtzBCElementBase<ELEMENT>::
FourierDecomposedHelmholtzBCElementBase(FiniteElement* const &bulk_el_pt, 
                        const int &face_index) : 
 FaceGeometry<ELEMENT>(), FaceElement()
 { 
  // Let the bulk element build the FaceElement, i.e. setup the pointers 
  // to its nodes (by referring to the appropriate nodes in the bulk
  // element), etc.
  bulk_el_pt->build_face_element(face_index,this);
  
  //Set up U_index_fourier_decomposedhelmholtz. 
  U_index_fourier_decomposed_helmholtz = std::complex<unsigned>(0,1);
  
  //Cast to the appropriate FourierDecomposedHelmholtzEquation so that we can
  //find the index at which the variable is stored
  //We assume that the dimension of the full problem is the same
  //as the dimension of the node, if this is not the case you will have
  //to write custom elements, sorry
  FourierDecomposedHelmholtzEquations* eqn_pt = 
   dynamic_cast<FourierDecomposedHelmholtzEquations*>(bulk_el_pt);
  if(eqn_pt==0)
   {
    std::string error_string =
     "Bulk element must inherit from FourierDecomposedHelmholtzEquations.";
    throw OomphLibError(
     error_string,
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }
  //Otherwise read out the value
  else
   {  
    //Read the index from the (cast) bulk element
    U_index_fourier_decomposed_helmholtz = 
     eqn_pt->u_index_fourier_decomposed_helmholtz();
   }
 }
}

#endif