/*********************************************************************/
/* File:   pml.cpp                                                   */
/* Author: Joachim Schoeberl                                         */
/* Date:   Nov 3, 2003                                               */
/*********************************************************************/

/*
  Billinear form integrator over complex domain
*/



// #include <fem.hpp>
#include <solve.hpp>
  

extern ngsolve::PDE * pde;




namespace ngfem
{
  using namespace ngfem;

  static Complex alpha(0,1);
  static double pml_r = 1;
  static double pml_x = 1;

  static double pml_xmin[3] = { -1, -1, -1};
  static double pml_xmax[3] = {  1,  1,  1};

  /*
    rect_pml = 0 .... circular pml with radius pml_r
    rect_pml = 1 .... square pml on square (-pml_x, pml_x)^d
    rect_pml = 2 .... rectangular pml on (pml_xmin, pml_xmax) x (pml_ymin, pml_ymax) x (pml_zmin, pml_zmax)
   */
  static int rect_pml = 0;
  

  
  template <>
  SpecificIntegrationPoint<2,2,Complex> :: 
  SpecificIntegrationPoint (const IntegrationPoint & aip,
			    const ElementTransformation & aeltrans,
			    LocalHeap & lh)
    : DimSpecificIntegrationPoint<2,Complex> (aip, aeltrans)
  {
    Mat<2,2> hdxdxi;
    Vec<2> hpoint;

    eltrans.CalcPointJacobian (ip, hpoint, hdxdxi, lh);


    switch  (rect_pml)
      {
      case 0:
	{
	  double abs_x = L2Norm (hpoint);
	  if (abs_x <= pml_r)
	    {
	      point = hpoint;
	      dxdxi = hdxdxi;
	    }
	  else
	    {
	      Complex g = 1.+alpha*(1.0-pml_r/abs_x);
	      point = g * hpoint;
	      Mat<2,2,Complex> trans =
		g * Id<2>() + (pml_r*alpha/(abs_x*abs_x*abs_x)) * (hpoint * Trans(hpoint));
	      dxdxi = trans * hdxdxi;
	      cout << "hi, bin ich hier " << endl;
	    }
	  break;
	}
      case 1:
	{
	  Vec<2> dabs_dpoint;
	  double abs_x;
	  if (fabs (hpoint(0)) > fabs(hpoint(1)))
	    {
	      abs_x = fabs (hpoint(0));
	      dabs_dpoint(0) = (hpoint(0) > 0) ? 1 : -1;
	      dabs_dpoint(1) = 0;
	    }
	  else
	    {
	      abs_x = fabs (hpoint(1));
	      dabs_dpoint(0) = 0;
	      dabs_dpoint(1) = (hpoint(1) > 0) ? 1 : -1;
	    }
	  
	  if (abs_x <= pml_x)   //  || hpoint(1) < 0)
	    {
	      point = hpoint;
	      dxdxi = hdxdxi;
	    }
	  else
	    {
	      Complex g = 1.+alpha*(1.0-pml_x/abs_x);
	      point = g * hpoint;
	      Mat<2,2,Complex> trans =
		g * Id<2>() + (pml_r*alpha/(abs_x*abs_x)) * (hpoint * Trans(dabs_dpoint));
	      dxdxi = trans * hdxdxi;
	    }
	  break;
	}

      case 2:
	{
	  point = hpoint;
	  dxdxi = hdxdxi;

	  for (int i = 0; i < 2; i++)
	    {
	      if (hpoint(i) > pml_xmax[i])
		{
		  point(i) += alpha * (point(i) - pml_xmax[i]);

		  Mat<2,2,Complex> trans = Id<2>();
		  trans(i,i) += alpha;
		  Mat<2,2,Complex> hm = dxdxi;
		  dxdxi = trans * hm;
		}

	      else if (hpoint(i) < pml_xmin[i])
		{
		  point(i) -= alpha * (pml_xmin[i] - point(i));
		  
		  Mat<2,2,Complex> trans = Id<2>();
		  trans(i,i) += alpha;
		  Mat<2,2,Complex> hm = dxdxi;
		  dxdxi = trans * hm;
		}
	    }
	  break;
	}
      default:
        break;
      }
	
    det = Det (dxdxi);
    dxidx = Inv (dxdxi);
  }



  
  template <>
  SpecificIntegrationPoint<3,3,Complex> :: 
  SpecificIntegrationPoint (const IntegrationPoint & aip,
			    const ElementTransformation & aeltrans,
			    LocalHeap & lh)
    : DimSpecificIntegrationPoint<3,Complex> (aip, aeltrans)
  {
    Mat<3,3> hdxdxi;
    Vec<3> hpoint;

    eltrans.CalcPointJacobian (ip, hpoint, hdxdxi, lh);
    // eltrans.CalcJacobian (ip, hdxdxi, lh);
    // eltrans.CalcPoint (ip, hpoint, lh);

    
    double abs_x = L2Norm (hpoint);
    switch (rect_pml)
      {
      case 0:
	{
	  if (abs_x <= pml_r)
	    {
	      point = hpoint;
	      dxdxi = hdxdxi;
	    }
	  else
	    {
	      Complex g = 1.+alpha*(1-pml_r/abs_x);
	      point = g * hpoint;
	      Mat<3,3,Complex> trans =
		g * Id<3>() + (pml_r*alpha/(abs_x*abs_x*abs_x)) * (hpoint * Trans(hpoint));
	      dxdxi = trans * hdxdxi;
	      
	      // dxdxi = g * Id<3>() + (alpha/(abs_x*abs_x*abs_x)) * (hpoint * Trans(hpoint));
	    }
	  break;
	}
      case 1:
	{
	  double abs_x = fabs (hpoint(0));
	  if (abs_x <= 3.0)
	    {
	      point = hpoint;
	      dxdxi = hdxdxi;
	    }
	  else
	    {
	      Complex g = 1.+ alpha * (1.0-3.0/abs_x);  
	      point = g * hpoint;
	    dxdxi = g * Id<3>();
	    dxdxi(0,0) += 3.0*alpha/(abs_x*abs_x*abs_x) * hpoint(0);
	    }
	  break;
	}
      case 2:
	{
	  point = hpoint;
	  dxdxi = hdxdxi;

	  for (int i = 0; i < 3; i++)
	    {
	      if (hpoint(i) > pml_xmax[i])
		{
		  point(i) += alpha * (point(i) - pml_xmax[i]);

		  Mat<2,2,Complex> trans = Id<2>();
		  trans(i,i) += alpha;
		  Mat<2,2,Complex> hm = dxdxi;
		  dxdxi = trans * hm;
		}

	      else if (hpoint(i) < pml_xmin[i])
		{
		  point(i) -= alpha * (pml_xmin[i] - point(i));
		  
		  Mat<2,2,Complex> trans = Id<2>();
		  trans(i,i) += alpha;
		  Mat<2,2,Complex> hm = dxdxi;
		  dxdxi = trans * hm;
		}
	    }
	  break;
	}

      default:
	break;
      }
    
    det = Det (dxdxi);
    dxidx = Inv (dxdxi);
  }



  template class SpecificIntegrationPoint<2,2,Complex>;
  template class SpecificIntegrationPoint<3,3,Complex>;




















  

template <class DIFFOP, class DMATOP, class FEL = FiniteElement>
class PML_BDBIntegrator : public T_BDBIntegrator<DIFFOP,DMATOP,FEL>
{
public:
  
  enum { DIM_SPACE   = DIFFOP::DIM_SPACE };
  enum { DIM_ELEMENT = DIFFOP::DIM_ELEMENT };
  enum { DIM_DMAT    = DIFFOP::DIM_DMAT };
  enum { DIM         = DIFFOP::DIM };
  
  ///
  PML_BDBIntegrator (const DMATOP & admat)
    : T_BDBIntegrator<DIFFOP,DMATOP,FEL> (admat)
  { 
    if (pde->GetConstantTable().Used ("pml_r"))
      pml_r = pde->GetConstant("pml_r");
    
    if (pde->GetConstantTable().Used ("pml_x"))
      {
	pml_x = pde->GetConstant("pml_x");
	rect_pml = 1;
      }


    if (pde->GetConstantTable().Used ("pml_xmin"))
      {
	pml_xmin[0] = pde->GetConstant("pml_xmin");
	rect_pml = 2;
      }
    if (pde->GetConstantTable().Used ("pml_xmax"))
      {
	pml_xmax[0] = pde->GetConstant("pml_xmax");
	rect_pml = 2;
      }

    if (pde->GetConstantTable().Used ("pml_ymin"))
      {
	pml_xmin[1] = pde->GetConstant("pml_ymin");
	rect_pml = 2;
      }
    if (pde->GetConstantTable().Used ("pml_ymax"))
      {
	pml_xmax[1] = pde->GetConstant("pml_ymax");
	rect_pml = 2;
      }

    if (pde->GetConstantTable().Used ("pml_zmin"))
      {
	pml_xmin[2] = pde->GetConstant("pml_zmin");
	rect_pml = 2;
      }
    if (pde->GetConstantTable().Used ("pml_zmax"))
      {
	pml_xmax[2] = pde->GetConstant("pml_zmax");
	rect_pml = 2;
      }


    switch (rect_pml)
      {
      case 0: cout << "circular pml of radius " << pml_r << endl; break;
      case 1: cout << "square pml on +/- " << pml_x << endl; break;
      case 2: cout << "rectangular pml on " 
		   << "(" << pml_xmin[0] << "," << pml_xmax[0] 
		   << ") x (" << pml_xmin[1] << "," << pml_xmax[1] 
		   << ") x (" << pml_xmin[2] << "," << pml_xmax[2] 
		   << endl; 
	break;
      }

    if (pde->GetConstantTable().Used ("pml_alpha"))
      alpha = Complex (0, pde->GetConstant("pml_alpha"));
  }

  ///
  virtual ~PML_BDBIntegrator ()
  { ; }


  ///
  virtual void
  AssembleElementMatrix (const FiniteElement & bfel, 
			 const ElementTransformation & eltrans, 
			 FlatMatrix<double> & elmat,
			 LocalHeap & locheap) const
  {
    throw Exception ("PML cannot generae real matrices");
  }

  ///
  virtual void
  AssembleElementMatrix (const FiniteElement & bfel, 
			 const ElementTransformation & eltrans, 
			 FlatMatrix<Complex> & elmat,
			 LocalHeap & locheap) const
  {
    try
      {
	const FEL & fel = dynamic_cast<const FEL&> (bfel);
	int ndof = fel.GetNDof();

	elmat.AssignMemory (ndof*DIM, ndof*DIM, locheap);
	elmat = 0;
	
	FlatMatrixFixHeight<DIM_DMAT, Complex> bmat (ndof * DIM, locheap);
	FlatMatrixFixHeight<DIM_DMAT, Complex> dbmat (ndof * DIM, locheap);

	Mat<DIM_DMAT,DIM_DMAT, Complex> dmat;

	const IntegrationRule & ir = GetIntegrationRule (fel);

	void * heapp = locheap.GetPointer();
	for (int i = 0; i < ir.GetNIP(); i++)
	  {
	    SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE,Complex> 
	      sip(ir[i], eltrans, locheap);
	    SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE,double> 
	      sip_real(ir[i], eltrans, locheap);

	    DIFFOP::GenerateMatrix (fel, sip, bmat, locheap);
	    this->dmatop.GenerateMatrix (fel, sip_real, dmat, locheap);

	    Complex fac = sip.GetJacobiDet() * sip.IP().Weight();

	    dbmat = dmat * bmat;
	    elmat += fac * (Trans (dbmat) * bmat);

	    locheap.CleanUp (heapp);
	  } 
      }

    catch (Exception & e)
      {
	e.Append ("in AssembleElementMatrix, type = ");
	e.Append (typeid(*this).name());
	e.Append ("\n");
	throw;
      }
    catch (exception & e)
      {
	Exception e2(e.what());
	e2.Append ("\nin AssembleElementMatrix, type = ");
	e2.Append (typeid(*this).name());
	e2.Append ("\n");
	throw e2;
      }
  }


  ///
  virtual int GetDimension () const
  { return DIM; }

  ///
  virtual int Lumping () const
    { return 0; }
  ///
  virtual string Name () const { return "BDB integrator"; }
};




/// 
template <int D, typename FEL = NodalFiniteElement>
class PML_LaplaceIntegrator 
  : public PML_BDBIntegrator<DiffOpGradient<D>, DiagDMat<D>, FEL>
{
public:
  ///
  PML_LaplaceIntegrator (CoefficientFunction * coeff)
    : PML_BDBIntegrator<DiffOpGradient<D>, DiagDMat<D>, FEL> (DiagDMat<D> (coeff))
  {
    ;
  }
  
  static Integrator * Create (ARRAY<CoefficientFunction*> & coeffs)
  {
    return new PML_LaplaceIntegrator (coeffs[0]);
  }

  ///
  virtual string Name () const { return "PML_Laplace"; }
};




/// 
template <int D, typename FEL = NodalFiniteElement>
class PML_MassIntegrator 
  : public PML_BDBIntegrator<DiffOpId<D>, DiagDMat<1>, FEL>
{
public:
  ///
  PML_MassIntegrator (CoefficientFunction * coeff)
    : PML_BDBIntegrator<DiffOpId<D>, DiagDMat<1>, FEL> (DiagDMat<1> (coeff))
  {
    ;
  }
  
  static Integrator * Create (ARRAY<CoefficientFunction*> & coeffs)
  {
    return new PML_MassIntegrator (coeffs[0]);
  }

  ///
  virtual string Name () const { return "PML_Mass"; }
};




///
class PML_ElasticityIntegrator
  : public PML_BDBIntegrator<DiffOpStrain<2>, PlaneStressDMat, NodalFiniteElement>
{
public:
  ///
  PML_ElasticityIntegrator (CoefficientFunction * coefe,
			    CoefficientFunction * coefnu)
    : PML_BDBIntegrator<DiffOpStrain<2>, PlaneStressDMat, NodalFiniteElement> 
  (PlaneStressDMat (coefe, coefnu))
  { ; }
  
  static Integrator * Create (ARRAY<CoefficientFunction*> & coeffs)
  {
    return new PML_ElasticityIntegrator (coeffs[0], coeffs[1]);
  }



  ///
  virtual string Name () const { return "Elasticity"; }
};






///
class PML_CurlCurlEdgeIntegrator
  : public PML_BDBIntegrator<DiffOpCurlEdge<3>, DiagDMat<3>, HCurlFiniteElementD<3> >
{
public:
  ///
  PML_CurlCurlEdgeIntegrator (CoefficientFunction * coef)
    : PML_BDBIntegrator<DiffOpCurlEdge<3>, DiagDMat<3>, HCurlFiniteElementD<3> > 
  (DiagDMat<3> (coef))
  { ; }
  
  static Integrator * Create (ARRAY<CoefficientFunction*> & coeffs)
  {
    return new PML_CurlCurlEdgeIntegrator (coeffs[0]);
  }

  ///
  virtual string Name () const { return "PML_CurlCurlEdge"; }
};


///
class PML_MassEdgeIntegrator
  : public PML_BDBIntegrator<DiffOpIdEdge<3>, DiagDMat<3>, HCurlFiniteElementD<3> >
{
public:
  ///
  PML_MassEdgeIntegrator (CoefficientFunction * coef)
    : PML_BDBIntegrator<DiffOpIdEdge<3>, DiagDMat<3>, HCurlFiniteElementD<3> > 
  (DiagDMat<3> (coef))
  { ; }
  
  static Integrator * Create (ARRAY<CoefficientFunction*> & coeffs)
  {
    return new PML_MassEdgeIntegrator (coeffs[0]);
  }

  ///
  virtual string Name () const { return "PML_Massedge"; }
};


















class PML_ReducedMassIntegrator : public BilinearFormIntegrator
{
  DiagDMat<1> dmatop;
  CoefficientFunction * coef;
public:
  
  enum { DIM_SPACE   = 2 };
  enum { DIM_ELEMENT = 2 };
  enum { DIM_DMAT    = 1 };
  enum { DIM         = 1 };
  

  PML_ReducedMassIntegrator (CoefficientFunction * acoef)
    : BilinearFormIntegrator(), dmatop(acoef), coef(acoef)
  { ; }

  ///
  virtual ~PML_ReducedMassIntegrator ()
  { ; }

  static Integrator * Create (ARRAY<CoefficientFunction*> & coeffs)
  {
    return new PML_ReducedMassIntegrator (coeffs[0]);
  }

  virtual bool BoundaryForm () const
  { return 0; }

  ///
  virtual void
  AssembleElementMatrix (const FiniteElement & bfel, 
			 const ElementTransformation & eltrans, 
			 FlatMatrix<double> & elmat,
			 LocalHeap & locheap) const
  {
    throw Exception ("PML cannot generae real matrices");
  }

  ///
  virtual void
  AssembleElementMatrix (const FiniteElement & bfel, 
			 const ElementTransformation & eltrans, 
			 FlatMatrix<Complex> & elmat,
			 LocalHeap & locheap) const
  {

    FE_Trig2 trig;
    NodalFiniteElement & felq = trig;
    const NodalFiniteElement & fel = dynamic_cast<const NodalFiniteElement&> (bfel);
    int ndof = fel.GetNDof();
    int ndofq = felq.GetNDof();

    elmat.AssignMemory (ndof*DIM, ndof*DIM, locheap);
    elmat = 0;
	
    FlatMatrix<Complex> matb(ndofq, ndof, locheap);
    FlatMatrix<Complex> matc(ndofq, ndofq, locheap);
    FlatMatrix<Complex> matcinv(ndofq, ndofq, locheap);


    const IntegrationRule & ir = 
      GetIntegrationRules().SelectIntegrationRule (fel.ElementType(), 6);
    
    matb = 0;
    matc = 0;
    void * heapp = locheap.GetPointer();
    for (int i = 0; i < ir.GetNIP(); i++)
      {
	SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE,Complex> 
	  sip(ir[i], eltrans, locheap);
	SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE,double> 
	  sip_real(ir[i], eltrans, locheap);
	
	
	const FlatVector<> shape = fel.GetShape (sip.IP(), locheap);
	const FlatVector<> shapeq = felq.GetShape (sip.IP(), locheap);
	

	double val = Evaluate (*coef, sip_real);
	Complex fac = sip.GetJacobiDet() * sip.IP().Weight();
	
	matb += fac * (shapeq * Trans (shape));
	matc += (fac/val) * (shapeq * Trans (shapeq));
	locheap.CleanUp (heapp);
      }

    CalcInverse (matc, matcinv);
    elmat = Trans (matb) * matcinv * matb;

    /*
    try
      {
	const NodalFiniteElement & fel = dynamic_cast<const NodalFiniteElement&> (bfel);
	int ndof = fel.GetNDof();

	elmat.AssignMemory (ndof*DIM, ndof*DIM, locheap);
	elmat = 0;
	
	FlatMatrixFixHeight<DIM_DMAT, Complex> bmat (ndof * DIM, locheap);
	FlatMatrixFixHeight<DIM_DMAT, Complex> dbmat (ndof * DIM, locheap);

	Mat<DIM_DMAT,DIM_DMAT, Complex> dmat;

	const IntegrationRule & ir = 
	  GetIntegrationRules().SelectIntegrationRule (fel.ElementType(), 2);

	void * heapp = locheap.GetPointer();
	for (int i = 0; i < ir.GetNIP(); i++)
	  {
	    SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE,Complex> 
	      sip(ir[i], eltrans, locheap);
	    SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE,double> 
	      sip_real(ir[i], eltrans, locheap);


	    const FlatVector<> shape = fel.GetShape (sip.IP(), locheap);
	    for (int j = 0; j < shape.Height(); j++)
	      bmat(0, j) = shape(j);
	    
	    dmatop.GenerateMatrix (fel, sip_real, dmat, locheap);
	    Complex fac = sip.GetJacobiDet() * sip.IP().Weight();

	    dbmat = dmat * bmat;
	    elmat += fac * (Trans (dbmat) * bmat);

	    locheap.CleanUp (heapp);
	  } 
      }

    catch (Exception & e)
      {
	e.Append ("in AssembleElementMatrix, type = ");
	e.Append (typeid(*this).name());
	e.Append ("\n");
	throw;
      }
    catch (exception & e)
      {
	Exception e2(e.what());
	e2.Append ("\nin AssembleElementMatrix, type = ");
	e2.Append (typeid(*this).name());
	e2.Append ("\n");
	throw e2;
      }
    */
  }


  ///
  virtual int GetDimension () const
  { return 1; }

  ///
  virtual int Lumping () const
    { return 0; }
  ///
  virtual string Name () const { return "PML_ReducedMass integrator"; }
};












class PML_InterpolationMassIntegrator : public BilinearFormIntegrator
{
  DiagDMat<1> dmatop;
  CoefficientFunction * coef;
public:
  
  enum { DIM_SPACE   = 2 };
  enum { DIM_ELEMENT = 2 };
  enum { DIM_DMAT    = 1 };
  enum { DIM         = 1 };
  

  PML_InterpolationMassIntegrator (CoefficientFunction * acoef)
    : BilinearFormIntegrator(), dmatop(acoef), coef(acoef)
  { ; }

  ///
  virtual ~PML_InterpolationMassIntegrator ()
  { ; }

  static Integrator * Create (ARRAY<CoefficientFunction*> & coeffs)
  {
    return new PML_InterpolationMassIntegrator (coeffs[0]);
  }

  virtual bool BoundaryForm () const
  { return 0; }

  ///
  virtual void
  AssembleElementMatrix (const FiniteElement & bfel, 
			 const ElementTransformation & eltrans, 
			 FlatMatrix<double> & elmat,
			 LocalHeap & locheap) const
  {
    throw Exception ("PML cannot generae real matrices");
  }

  ///
  virtual void
  AssembleElementMatrix (const FiniteElement & bfel, 
			 const ElementTransformation & eltrans, 
			 FlatMatrix<Complex> & elmat,
			 LocalHeap & locheap) const
  {

    FE_Trig1 trig;
    NodalFiniteElement & felq = trig;
    const NodalFiniteElement & fel = dynamic_cast<const NodalFiniteElement&> (bfel);
    int ndof = fel.GetNDof();
    int ndofq = felq.GetNDof();

    elmat.AssignMemory (ndof*DIM, ndof*DIM, locheap);
    elmat = 0;

    FlatMatrix<Complex> matm(ndofq, ndofq, locheap);

    FlatMatrix<Complex> matb(ndofq, ndof, locheap);
    FlatMatrix<Complex> matb2(ndofq, ndof, locheap);
    FlatMatrix<Complex> matc(ndofq, ndofq, locheap);
    FlatMatrix<Complex> matcinv(ndofq, ndofq, locheap);

    const IntegrationRule & ir = 
      GetIntegrationRules().SelectIntegrationRule (fel.ElementType(), 6);
    const IntegrationRule & irseg = 
      GetIntegrationRules().SelectIntegrationRule (ET_SEGM, 0);
    
    matm = 0;
    matb = 0;
    matc = 0;

    FlatVector<> shapee (3, locheap);


    void * heapp = locheap.GetPointer();
    for (int i = 0; i < ir.GetNIP(); i++)
      {
	SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE,Complex> 
	  sip(ir[i], eltrans, locheap);
	SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE,double> 
	  sip_real(ir[i], eltrans, locheap);
	
	const FlatVector<> shape = fel.GetShape (sip.IP(), locheap);
	const FlatVector<> shapeq = felq.GetShape (sip.IP(), locheap);

	double val = Evaluate (*coef, sip_real);
	Complex fac = sip.GetJacobiDet() * sip.IP().Weight() * val;
	
	matm += fac * (shapeq * Trans (shapeq));
	locheap.CleanUp (heapp);
      }


    const EDGE * edges = ElementTopology::GetEdges (ET_TRIG);
    const POINT3D * points = ElementTopology::GetVertices (ET_TRIG);
    for (int en = 0; en < 3; en++)
      for (int i = 0; i < irseg.GetNIP(); i++)
	{
	  int v1 = edges[en][0];
	  int v2 = edges[en][1];
	  double lam = irseg[i](0);
	  IntegrationPoint ip(lam*points[v1][0]+(1-lam)*points[v2][0],
			      lam*points[v1][1]+(1-lam)*points[v2][1], 0, 0);
	  
	  SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE,Complex> 
	    sip(ip, eltrans, locheap);
	  SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE,double> 
	    sip_real(ip, eltrans, locheap);
	  
	  const FlatVector<> shape = fel.GetShape (sip.IP(), locheap);
	  const FlatVector<> shapeq = felq.GetShape (sip.IP(), locheap);

	  shapee = 0.0;
	  shapee(en) = 1;

	  Complex fac = irseg[i].Weight();

	  matb += fac * (shapee * Trans (shape));
	  matc += fac * (shapee * Trans (shapeq));
	  locheap.CleanUp (heapp);
	}


    CalcInverse (matc, matcinv);
    matb2 = matcinv * matb;
    elmat = Trans (matb2) * matm * matb2;

    (*testout) << "matb = " << endl << matb << endl
	       << "matc = " << endl << matc << endl
	       << "matb2 = " << endl << matb2 << endl
	       << "matm = " << endl << matm << endl;
  }


  ///
  virtual int GetDimension () const
  { return 1; }

  ///
  virtual int Lumping () const
    { return 0; }
  ///
  virtual string Name () const { return "PML_InterpolationMass integrator"; }
};




















  namespace {
    class Init
    { 
    public: 
      Init ();
    };        
    

    Init::Init()
    {
      GetIntegrators().AddBFIntegrator ("PML_laplace", 2, 1,
					PML_LaplaceIntegrator<2>::Create);
      GetIntegrators().AddBFIntegrator ("PML_mass", 2, 1,
					PML_MassIntegrator<2>::Create);
      GetIntegrators().AddBFIntegrator ("PML_elasticity", 2, 2,
					PML_ElasticityIntegrator::Create);
      GetIntegrators().AddBFIntegrator ("PML_laplace", 3, 1,
					PML_LaplaceIntegrator<3>::Create);
      GetIntegrators().AddBFIntegrator ("PML_mass", 3, 1,
					PML_MassIntegrator<3>::Create);
      GetIntegrators().AddBFIntegrator ("PML_reducedmass", 2, 1,
					PML_ReducedMassIntegrator::Create);
      GetIntegrators().AddBFIntegrator ("PML_interpolationmass", 2, 1,
					PML_InterpolationMassIntegrator::Create);

      GetIntegrators().AddBFIntegrator ("PML_curlcurledge", 3, 1,
					PML_CurlCurlEdgeIntegrator::Create);
      GetIntegrators().AddBFIntegrator ("PML_massedge", 3, 1,
					PML_MassEdgeIntegrator::Create);
    }
    
    Init init;
    
  }
 }