pml_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: 1282 $
//LIC//
//LIC// $LastChangedDate: 2017-01-16 08:27:53 +0000 (Mon, 16 Jan 2017) $
//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//====================================================================
#ifndef OOMPH_PML_ELEMENT_HEADER
#define OOMPH_PML_ELEMENT_HEADER

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

#include "pml_mapping_functions.h"
#include "elements.h"
#include "complex_matrices.h"

#include <complex>

namespace oomph
{

// typedef used to avoid inaccessible error in Rip struct
typedef FiniteElement* FiniteElementPt;
typedef GeneralisedElement GeneralisedElementTypedef;

//=======================================================================
/// Policy class for defining the equivalent quad element for any given
/// element. Used for constructing a structured PML around an unstructured
/// mesh. Needs to be instantiated for each bulk element class
//=======================================================================
template<class ELEMENT>
class EquivalentQElement
{
};

//==============================================================
/// Base class for elements with pml capabilities
//==============================================================
template<unsigned DIM>
class PMLElementBase : public virtual FiniteElement
{

public:

  /// Constructor
  PMLElementBase() : Pml_is_enabled(false) {}

  /// Virtual destructor
  virtual ~PMLElementBase(){}

  /// \short Disable pml
  virtual void disable_pml()
  {
    Pml_is_enabled=false;
  }

  /// \short Enable pml
  virtual void enable_pml()
  {
    Pml_is_enabled=true;
  }

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DenseComplexMatrix& jacobian) = 0;

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DenseComplexMatrix& jacobian_matrix,
    Vector<std::complex<double> >& transformed_x) = 0;

  /// \short Pure virtual function in which we have to specify the
  /// values to be pinned (and set to zero) on the outer edge of
  /// the pml layer. This is usually all of the nodal values
  /// (values 0 and 1 (real and imag part) for Helmholtz; 
  /// values 0,1,2 and 3 (real and imag part of x- and y-displacement
  /// for 2D time-harmonic linear elasticity; etc.). Vector
  /// must be resized internally!
  virtual void values_to_be_pinned_on_outer_pml_boundary(
    Vector<unsigned>& values_to_pin)=0;

  virtual double wavenumber() const
  {
    throw OomphLibError("PML elements must implement wavenumber()",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }

  /// \short Computes laplace matrix (JJ^T)^-1|J| and determinant
  /// |J| for jacobian J. Laplace matrix is so named because it
  /// appears as the factor when the Laplace operator is transformed 
  static void compute_laplace_matrix_and_det(
    const DenseComplexMatrix& jacobian,
    DenseComplexMatrix& laplace_matrix,
    std::complex<double>& jacobian_det);

  /// \short Computes laplace matrix (JJ^T)^-1|J| and determinant
  /// |J| for diagonal matrix jacobian J. Laplace matrix is so named because it
  /// appears as the factor when the Laplace operator is transformed 
  static void compute_laplace_matrix_and_det(
    const DiagonalComplexMatrix& jacobian,
    DiagonalComplexMatrix& laplace_matrix,
    std::complex<double>& jacobian_det);

  /// \short Computes laplace matrix (JJ^T)^-1|J| and determinant
  /// |J| for diagonal matrix jacobian J. Laplace matrix is so named because it
  /// appears as the factor when the Laplace operator is transformed 
  static void compute_laplace_matrix_and_det(
    const DiagonalComplexMatrix& jacobian,
    DenseComplexMatrix& laplace_matrix,
    std::complex<double>& jacobian_det);

  /// \short Computes the inverse J^-1 and determinant |J| of a Jacobian J
  static void compute_jacobian_inverse_and_det(
    const DiagonalComplexMatrix& jacobian,
    DiagonalComplexMatrix& jacobian_inverse,
    std::complex<double>& jacobian_det);

protected:

  /// Boolean indicating if element is used in pml mode
  bool Pml_is_enabled;

};



template<unsigned DIM>
class AxisAlignedPMLElement : public PMLElementBase<DIM>
{
public:
  
  /// Constructor
  AxisAlignedPMLElement() : 
    Pml_direction_active(DIM,false),
    Pml_inner_boundary(DIM,0.0),
    Pml_outer_boundary(DIM,0.0)
  {}

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DiagonalComplexMatrix& jacobian);
    
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DenseComplexMatrix& jacobian)
  {
    throw OomphLibError("AxisAlignedPMLElements produce diagonal matrices,"
                        "use the diagonal version of this function instead",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DiagonalComplexMatrix& jacobian_matrix,
    Vector<std::complex<double> >& transformed_x);

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DenseComplexMatrix& jacobian_matrix,
    Vector<std::complex<double> >& transformed_x)
  {
    throw OomphLibError("AxisAlignedPMLElements produce diagonal matrices,"
                        "use the diagonal version of this function instead",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }

  /// \short Disable pml. Ensures the PML-ification in all directions
  /// has been deactivated
  virtual void disable_pml()
  {
    // Disable base class attributes
    PMLElementBase<DIM>::disable_pml();

    // Loop over the entries in Pml_direction_active and deactivate the
    // PML in this direction
    for (unsigned direction=0;direction<DIM;direction++)
    {
      // Disable the PML in this direction
      Pml_direction_active[direction]=false;
    }
  } // End of disable_pml
 
  /// \short Enable pml. Specify the coordinate direction along which 
  /// pml boundary is constant, as well as the coordinate
  /// along the dimension for the interface between the physical and artificial
  /// domains and the coordinate for the outer boundary.
  /// All of these are used to adjust the perfectly matched layer
  /// mechanism. Needs to be called separately for each pml-ified direction
  /// (if needed -- e.g. in corner elements)
  virtual void enable_pml(const int& direction,
                          const double& pml_inner_boundary, 
                          const double& pml_outer_boundary)
  {
    PMLElementBase<DIM>::enable_pml();
    Pml_direction_active[direction] = true;
    Pml_inner_boundary[direction] = pml_inner_boundary;
    Pml_outer_boundary[direction] = pml_outer_boundary;
  }

  virtual void enable_pml(const int& direction,
                          const double& pml_inner_boundary, 
                          const double& pml_outer_boundary,
                          UniaxialPMLMapping* pml_mapping_pt)
  {
    enable_pml(direction, pml_inner_boundary, pml_outer_boundary);
    Pml_mapping_pt = pml_mapping_pt;
  }
  
  /// Return a pointer to the PML Mapping object
  UniaxialPMLMapping* &pml_mapping_pt() {return Pml_mapping_pt;}

  /// Return a pointer to the PML Mapping object (const version)
  UniaxialPMLMapping* const &pml_mapping_pt() const {return Pml_mapping_pt;}

  /// Vector which points from the inner boundary to x
  inline double nu(const Vector<double>& x, const unsigned& i) const
  {
    return x[i] - this->Pml_inner_boundary[i];
  }

  /// Vector which points from the inner boundary to x
  void nu(const Vector<double>& x, Vector<double> nu_vec) const
  {
    for(unsigned i=0; i<DIM; i++)
    {
      nu_vec[i] = nu(x, i);
    }
  }

  /// Vector which points from the inner boundary to x
  inline double delta(const unsigned& i) const
  {
    return this->Pml_outer_boundary[i] - this->Pml_inner_boundary[i];
  }

  /// Vector which points from the inner boundary to x
  void delta(const Vector<double>& x, Vector<double> delta) const
  {
    for(unsigned i=0; i<DIM; i++)
    {
      delta[i] = this->Pml_outer_boundary[i] - this->Pml_inner_boundary[i];
    }
  }

protected:
  
  /// \short Coordinate direction along which pml boundary is constant; 
  /// alternatively: coordinate direction in which coordinate stretching
  /// is performed.
  std::vector<bool> Pml_direction_active;

  /// \short Coordinate of inner pml boundary 
  /// (Storage is provided for any coordinate
  /// direction; only the entries for "active" directions is used.)
  Vector<double> Pml_inner_boundary;

  /// \short Coordinate of outer pml boundary 
  /// (Storage is provided for any coordinate
  /// direction; only the entries for "active" directions is used.)
  Vector<double> Pml_outer_boundary;
  
  UniaxialPMLMapping* Pml_mapping_pt;
};

class AnnularPMLElementBase : public PMLElementBase<2>
{

public:

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DenseComplexMatrix& jacobian) = 0;

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DenseComplexMatrix& jacobian_matrix,
    Vector<std::complex<double> >& transformed_x) = 0;

  //
  void radial_to_cartesian_jacobian(
    const double& r,
    const double& theta,
    const std::complex<double>& transformed_r,
    const std::complex<double>& dtransformed_r_dr,
    DenseComplexMatrix& cartesian_jacobian) const;

  //
  void radial_to_cartesian_jacobian(
    const double& r,
    const double& theta,
    const std::complex<double>& rt,
    const DenseComplexMatrix& radial_jacobian,
    DenseComplexMatrix& cartesian_jacobian) const;

  /// \short get the Jacobian of the mapping from transformed r to the radial
  /// through-the-pml coordinate. Also returns the transformed_r at this point too
  /// as it is required for the transformation to cartesian coordinates
  virtual void radial_transformation_jacobian(
    const double& nu,
    const double& delta,
    std::complex<double>& transformed_r,
    std::complex<double>& dtransformed_r_dr) const = 0;
  
  //
  virtual double radius(const Vector<double>& s, const Vector<double>& x) = 0;

  //
  virtual double nu(const Vector<double>& s, const Vector<double>& x) = 0;

  //
  virtual double theta(const Vector<double>& s, const Vector<double>& x) = 0;
  
  virtual double delta() = 0;
};


class AnnularFromCartesianPMLElement : public AnnularPMLElementBase
{
public:
  
  /// Constructor
  AnnularFromCartesianPMLElement() :
    Pml_inner_radius(0.0),
    Pml_outer_radius(0.0),
    Origin(2, 0.0),
    Pml_mapping_pt(0)
    {}

  /// Enable pml. Specify the inner and outer radius of PML and the origin
  virtual void enable_pml(
    const double& inner_pml_radius, 
    const double& outer_pml_radius,
    const Vector<double>& origin,
    UniaxialPMLMapping* pml_mapping_pt)
  {
    PMLElementBase<2>::enable_pml();
    Pml_inner_radius = inner_pml_radius;
    Pml_outer_radius = outer_pml_radius;
    Origin = origin;
    Pml_mapping_pt = pml_mapping_pt;
  }

  double radius(const Vector<double>& s, const Vector<double>& x)
  {
    return sqrt(std::pow(x[0]-Origin[0],2) + std::pow(x[1]-Origin[1],2));
  }

  /// Coordinate which goes through the PML
  double nu(const Vector<double>& s, const Vector<double>& x)
  {
    return nu(radius(s, x));
  }

  /// Coordinate which goes through the PML
  double nu(const double& r)
  {
    return r - Pml_inner_radius;
  }

  /// Thickness of PML
  double delta() { return Pml_outer_radius - Pml_inner_radius; }

  /// Angle relative to origin
  double theta(const Vector<double>& x)
  {
    return atan2(x[1]-Origin[1], x[0]-Origin[0]);
  }

  /// \short get the Jacobian of the mapping from transformed r to the radial
  /// through-the-pml coordinate. Also returns the transformed_r at this point too
  /// as it is required for the transformation to cartesian coordinates
  virtual void radial_transformation_jacobian(
    const double& nu,
    const double& delta,
    std::complex<double>& transformed_r,
    std::complex<double>& dtransformed_r_dr)
  {
    const double k = this->wavenumber();
    transformed_r = Pml_inner_radius
                    + Pml_mapping_pt->transformed_nu(nu, delta, k);
    dtransformed_r_dr = Pml_mapping_pt->dtransformed_nu_dnu(nu, delta, k);
  }

protected:
  
  /// \short Coordinate direction along which pml boundary is constant; 
  /// alternatively: coordinate direction in which coordinate stretching
  /// is performed.
  double Pml_inner_radius;

  /// \short Coordinate of inner pml boundary 
  /// (Storage is provided for any coordinate
  /// direction; only the entries for "active" directions is used.)
  double Pml_outer_radius;

  /// \short Coordinate of outer pml boundary 
  /// (Storage is provided for any coordinate
  /// direction; only the entries for "active" directions is used.)
  Vector<double> Origin;
  
  // A uniaxial PML mapping used to transform the radial coordinate
  UniaxialPMLMapping* Pml_mapping_pt;
};


class AnnularFromLocalCoordinatePMLElement : public AnnularPMLElementBase
{
public:
  
  /// Constructor
  AnnularFromLocalCoordinatePMLElement() :
    Nu_s_index(0),
    Nu_at_s_min(0.0),
    Nu_at_s_max(0.0),
    Theta_s_index(0),
    Theta_at_s_min(0.0),
    Theta_at_s_max(0.0),
    Pml_inner_radius(0.0),
    Pml_outer_radius(0.0)
  {}
    
  /// Constructor
  AnnularFromLocalCoordinatePMLElement(
    unsigned nu_s_index,
    double nu_at_s_min,
    double nu_at_s_max,
    unsigned theta_s_index,
    double theta_at_s_min,
    double theta_at_s_max,
    double pml_inner_radius,
    double pml_outer_radius) :
    Nu_s_index(nu_s_index),
    Nu_at_s_min(nu_at_s_min),
    Nu_at_s_max(nu_at_s_max),
    Theta_s_index(theta_s_index),
    Theta_at_s_min(theta_at_s_min),
    Theta_at_s_max(theta_at_s_max),
    Pml_inner_radius(pml_inner_radius),
    Pml_outer_radius(pml_outer_radius)
  {}

  /// Enable pml. Specify the inner and outer radius of PML and the origin
  virtual void enable_pml(const double& inner_pml_radius, 
      const double& outer_pml_radius, const Vector<double>& origin)
  {
    PMLElementBase<2>::enable_pml();
  }

  virtual double s_min();

  virtual double s_max();

  /// Coordinate which goes through the PML in terms of local coordinate
  double nu(const Vector<double>& s, const Vector<double>& x)
  {
    return (s[Nu_s_index]-s_min())/(s_max()-s_min())*(Nu_at_s_max-Nu_at_s_min)
           + Nu_at_s_min;
  }

  /// Thickness of PML
  double delta() { return Pml_outer_radius - Pml_inner_radius; }

  /// Angle relative to origin
  double theta(const Vector<double>& s, const Vector<double>& x)
  {
    return (s[Theta_s_index]-s_min())/(s_max()-s_min())
              *(Theta_at_s_max - Theta_at_s_min)
           + Theta_at_s_min;
  }

  /// \short get the Jacobian of the mapping from transformed r to the radial
  /// through-the-pml coordinate. Also returns the transformed_r at this point too
  /// as it is required for the transformation to cartesian coordinates
  virtual void radial_transformation_jacobian(
    const double& nu,
    const double& delta,
    std::complex<double>& transformed_r,
    std::complex<double>& dtransformed_r_dr)
  {
    const double k = this->wavenumber();
    transformed_r = Pml_inner_radius
                    + Pml_mapping_pt->transformed_nu(nu, delta, k);
    dtransformed_r_dr = Pml_mapping_pt->dtransformed_nu_dnu(nu, delta, k);
  }

protected:

  /// \short Index of local coordinate which nu (through the PML coordinate) is
  /// parametrised by
  unsigned Nu_s_index;

  /// Value of nu at the minimum value of local coordinate
  double Nu_at_s_min;

  /// Value of nu at the maximum value of local coordinate
  double Nu_at_s_max;

  /// Index of local coordinate which theta is parametrised by
  unsigned Theta_s_index;

  /// Value of theta at the minimum value of local coordinate
  double Theta_at_s_min;

  /// Value of theta at the maximum value of local coordinate
  double Theta_at_s_max;
  
  /// \short Value of radius at inner PML boundary
  double Pml_inner_radius;

  /// \short Value of radius at outer PML boundary
  double Pml_outer_radius;

  // A uniaxial PML mapping used to transform the radial coordinate
  UniaxialPMLMapping* Pml_mapping_pt;
};


class Conformal2DPMLElement : public PMLElementBase<2>
{
public:
  
  /// Constructor
  Conformal2DPMLElement() :
    Nu_s_index(0),
    Nu_at_s_min(0.0),
    Nu_at_s_max(0.0),
    Pml_mapping_pt(0)
  {}
    
  /// Constructor
  Conformal2DPMLElement(
    const unsigned nu_s_index,
    const double nu_at_s_min,
    const double nu_at_s_max) :
    Nu_s_index(nu_s_index),
    Nu_at_s_min(nu_at_s_min),
    Nu_at_s_max(nu_at_s_max),
    Pml_mapping_pt(0)
  {}

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DenseComplexMatrix& jacobian);

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DenseComplexMatrix& jacobian_matrix,
    Vector<std::complex<double> >& transformed_x);

  /// Enable pml. Set the pointer to a PML mapping function
  virtual void enable_pml(UniaxialPMLMapping* pml_mapping_pt)
  {
    PMLElementBase::enable_pml();
    pml_mapping_pt = Pml_mapping_pt;
  }

  /// Coordinate which goes through the PML in terms of local coordinate
  double nu(const Vector<double>& s)
  {
    return (s[Nu_s_index]-s_min())/(s_max()-s_min())*(Nu_at_s_max-Nu_at_s_min)
           + Nu_at_s_min;
  }

  //
  void assemble_conformal_jacobian(
    const Vector<double>& p,
    const Vector<double>& dx_inner_dacross,
    const Vector<double>& dp_dacross,
    const double& nu,
    const std::complex<double>& tnu,
    const std::complex<double>& dtnu_dnu,
    DenseComplexMatrix& mapping_jacobian);

  //
  unsigned nu_s_index() const { return Nu_s_index; }
  
  //
  unsigned across_s_index() const { return (Nu_s_index+1) % 2; }

protected:

  // Extract various vectors which characterise the PML and its geometry
  void get_pml_properties(const Vector<double>& s,
                          Vector<double>& x_inner,
                          Vector<double>& x_outer,
                          Vector<double>& p,
                          Vector<double>& dx_inner_dacross,
                          Vector<double>& dp_dacross,
                          double& delta);

  /// \short Index of local coordinate which nu (through the PML coordinate) is
  /// parametrised by
  unsigned Nu_s_index;

  /// Value of nu at the minimum value of local coordinate
  double Nu_at_s_min;

  /// Value of nu at the maximum value of local coordinate
  double Nu_at_s_max;
  
  UniaxialPMLMapping* Pml_mapping_pt;
};

class TangentiallyVaryingConformal2DPMLElement : public Conformal2DPMLElement
{
public:
  
  /// Constructor
  TangentiallyVaryingConformal2DPMLElement() : Conformal2DPMLElement()
  {}

  /// Constructor
  TangentiallyVaryingConformal2DPMLElement(
    const unsigned nu_s_index,
    const double nu_at_s_min,
    const double nu_at_s_max) :
    Conformal2DPMLElement(nu_s_index, nu_at_s_min, nu_at_s_max)
  {}

  /// Enable pml. Set the pointer to a PML mapping function
  virtual void enable_pml(TangentiallyVaryingConformalPMLMapping* pml_mapping_pt)
  {
    PMLElementBase<2>::enable_pml();
    Tangentially_varying_pml_mapping_pt = pml_mapping_pt;
  }

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DenseComplexMatrix& jacobian);

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DenseComplexMatrix& jacobian_matrix,
    Vector<std::complex<double> >& transformed_x);

  //
  void assemble_conformal_jacobian(
    const Vector<double>& p,
    const Vector<double>& dx_inner_dacross,
    const Vector<double>& dp_dacross,
    const double& nu,
    const std::complex<double>& tnu,
    const std::complex<double>& dtnu_dnu,
    const std::complex<double>& dtnu_dacross,
    DenseComplexMatrix& mapping_jacobian);

protected:
  
  // Extract various vectors which characterise the PML and its geometry
  void get_pml_properties(const double& s_across_the_pml,
                          Vector<double>& x_inner,
                          Vector<double>& x_outer,
                          Vector<double>& p,
                          Vector<double>& dx_inner_dacross,
                          Vector<double>& dp_dacross,
                          double& delta,
                          unsigned& across_the_pml_index,
                          unsigned& through_the_pml_index);

   TangentiallyVaryingConformalPMLMapping* Tangentially_varying_pml_mapping_pt;
};


class TangentiallyDiscontinuousConformal2DPMLElement :
  public TangentiallyVaryingConformal2DPMLElement
{
public:
  
  /// Constructor
  TangentiallyDiscontinuousConformal2DPMLElement() :
    TangentiallyVaryingConformal2DPMLElement()
  {}

  /// Constructor
  TangentiallyDiscontinuousConformal2DPMLElement(
    const unsigned nu_s_index,
    const double nu_at_s_min,
    const double nu_at_s_max) :
    TangentiallyVaryingConformal2DPMLElement(nu_s_index, nu_at_s_min,
      nu_at_s_max)
  {}

  /// Enable pml. Set the pointer to a PML mapping function
  virtual void enable_pml(TangentiallyDiscontinuousConformalPMLMapping* pml_mapping_pt)
  {
    PMLElementBase<2>::enable_pml();
    Tangentially_discontinuous_pml_mapping_pt = pml_mapping_pt;
  }

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DenseComplexMatrix& jacobian);

  ///
  virtual void pml_transformation_jacobian(
    const unsigned& ipt,
    const Vector<double>& s,
    const Vector<double>& x,
    DenseComplexMatrix& jacobian_matrix,
    Vector<std::complex<double> >& transformed_x);

  // Use this mode to hijack the assembly process to find rips
  static bool Find_rips_mode;

  /// Information about a rip in two dimensions
  struct TwoDimRip{
    /// Position of onset of rip in local and global coordinates
    FiniteElementPt ripped_element_pt;
    double s_across;
    double s_through;
    double x;
    double y;
    std::complex<double> tnu;
    // Factor which describes behaviour near rip through the PML
    // - \frac{partial^2 u}{\partial \tilde\nu^2} / U at rip
    std::complex<double> alpha;
    // Factor which describes behaviour near rip across the PML
    // - i \frac{partial^2 u}{\partial \tilde\nu^2}
    //     / (U'(1-\bar\nu) -\frac{\partial u}{\partial \zeta}) at rip
    std::complex<double> beta;

    // Constructor with all parameters
    TwoDimRip(FiniteElementPt Ripped_element_pt, const double& S_across,
      const double& S_through, const double& X, const double& Y,
      const std::complex<double> Tnu, const std::complex<double> Alpha,
      const std::complex<double> Beta):
      ripped_element_pt(Ripped_element_pt), s_across(S_across),
      s_through(S_through), x(X), y(Y), tnu(Tnu), alpha(Alpha), beta(Beta)
    { }
  };

  // Holds information about all known rips, which we will use to inform our
  // integration scheme
  static Vector<TwoDimRip> rips;

  /// Look for rips in the transformation and store in rips vector
  void find_rips();

  /// Compute distance from rip to PML element edge
  double  distance_from_element_edge( const Vector<double>& x_inner,
    const Vector<double>& x_outer, const Vector<double>& p,
    const Vector<double>& r);

  /// \short Compute n'J^{-1}, which when applied to the untransformed
  /// grad operator gives us the transformed normal derivative
  void compute_transformed_normal_derivative(
    const DenseComplexMatrix& J,
    const Vector<double>& n,
    Vector<std::complex<double> >& result);

  // You must overload this to do the normal things in your specific element
  virtual void fill_in_contribution_to_jacobian_no_rips(
    Vector<double> &residuals,
    DenseMatrix<double> &jacobian, bool flag) = 0;

  // You must overload this and handle the rips in your specific element
  virtual void fill_in_contribution_to_jacobian_with_rips(
    Vector<double> &residuals,
    DenseMatrix<double> &jacobian, bool flag) = 0;

  /// Add the element's contribution to its residual vector (wrapper)
  void fill_in_contribution_to_residuals(Vector<double> &residuals)
  {
    if(this->Pml_is_enabled && !Disable_rip_integration)
    {
      // Use this mode to hijack the assembly process to find rips
      if(Find_rips_mode)
      {
        find_rips();
        return;
      }

      //Call the generic residuals function with flag set to 0
      //using a dummy matrix argument
      fill_in_contribution_to_jacobian_with_rips(residuals,
        GeneralisedElementTypedef::Dummy_matrix, 0);
    }
    else
    {
      //Call the generic residuals function with flag set to 0
      //using a dummy matrix argument
      fill_in_contribution_to_jacobian_no_rips(
        residuals, GeneralisedElementTypedef::Dummy_matrix, 0);
    }
  }

  /// \short Add the element's contribution to its residual vector and
  /// element Jacobian matrix (wrapper)
  void fill_in_contribution_to_jacobian(
    Vector<double> &residuals,
    DenseMatrix<double> &jacobian)
  {
    if(this->Pml_is_enabled && !Disable_rip_integration)
    {
      // Use this mode to hijack the assembly process to find rips
      if(Find_rips_mode)
      {
        find_rips();
        return;
      }
      // Call the generic routine with the flag set to 1
      fill_in_contribution_to_jacobian_with_rips(
        residuals, jacobian, 1);
    }
    else
    {
      //Call the generic routine with the flag set to 1
      fill_in_contribution_to_jacobian_no_rips(residuals, jacobian,1);
    }
  }
protected:
  // Pretend that there is a rip at some point in the element to test that it has
  // no effect. This is overridden by disable rip integration
  bool Force_rip_integration;

  // Do not do the integration process which find and treats rips separately
  bool Disable_rip_integration;
  
  TangentiallyDiscontinuousConformalPMLMapping* Tangentially_discontinuous_pml_mapping_pt;
};

}

#endif