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// Copyright (C) 2006-2026 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
#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
{
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
 
 
  //======================================================================
  /// 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:
    /// 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);
 
    /// 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) = delete;
 
    /// 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&) = delete;*/
 
 
    /// 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);
    }
 
    /// 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);
    }
 
    /// 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);
    }
 
    /// 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());
    }
 
    /// 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);
    }
 
    /// 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:
    /// 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);
    }
 
    /// 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 equal to the shape functions
      test.shallow_copy_from(psi);
      dtest_ds.shallow_copy_from(dpsi_ds);
 
      // Return the value of the jacobian
      return J_eulerian(s);
    }
 
    /// 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)
    {
    }
 
    /// Compute and store the gamma integral at all integration
    /// points of the constituent elements.
    void setup_gamma();
 
    /// 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];
    }
 
    /// 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];
    }
 
    /// The outer radius
    double& outer_radius()
    {
      return Outer_radius;
    }
 
    /// 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;
 
 
    /// Container to store the gamma integral for given Gauss point
    /// and element
    std::map<FiniteElement*, Vector<std::complex<double>>> Gamma_at_gauss_point;
 
 
    /// 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:
    /// 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);
    }
 
 
    /// 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);
    }
 
    /// 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);
 
 
    /// 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;
    }
 
    /// Set mesh of all DtN boundary condition elements
    void set_outer_boundary_mesh_pt(
      FourierDecomposedHelmholtzDtNMesh<ELEMENT>* mesh_pt)
    {
      Outer_boundary_mesh_pt = mesh_pt;
    }
 
 
    /// 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:
    /// 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
    }
 
 
    /// 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==================
  /// 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
  }
 
 
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
 
 
  //===========================================================================
  /// Namespace for checking radius of nodes on (assumed to be circular)
  /// DtN boundary
  //===========================================================================
  namespace ToleranceForFourierDecomposedHelmholtzOuterBoundary
  {
    /// Relative tolerance to within radius of points on DtN boundary
    /// are allowed to deviate from specified value
    extern double Tol;
 
  } // namespace ToleranceForFourierDecomposedHelmholtzOuterBoundary
 
 
  //===========================================================================
  /// Namespace for checking radius of nodes on (assumed to be circular)
  /// DtN boundary
  //===========================================================================
  namespace ToleranceForFourierDecomposedHelmholtzOuterBoundary
  {
    /// Relative tolerance to within radius of points on DtN boundary
    /// are allowed to deviate from specified value
    double Tol = 1.0e-3;
 
  } // namespace ToleranceForFourierDecomposedHelmholtzOuterBoundary
 
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
 
 
  /// ================================================================
  /// 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();
    }
  }
} // namespace oomph
 
#endif