#ifndef FILE_BEM
#define FILE_BEM

/*********************************************************************/
/* File:   bem.hh                                                    */
/* Author: Joachim Schoeberl                                         */
/* Date:   5. Dec. 2001                                              */
/*********************************************************************/

/** 
    BEM 
    Bilinear-form for boundary element methods
*/

class BEMElement : public SpecialElement
{
protected:
  ///fe spaces for Dirichlet data
  const FESpace & fes;
  /// fe space for Neumann data
  const FESpace & fes_neumann;

  /// node renumbering:
  ARRAY<int> bem2fem;
  /// node renumbering, reverse:
  ARRAY<int> fem2bem;
  /// number of fem dofs
  int dimfem;
  /// number of Dirichlet - BEM dofs
  int dimd;
  /// number of Neumann - BEM dofs
  int dimn;
public:
  BEMElement (const FESpace & afes,
	      const FESpace & afespacen);

  /*
  virtual void Assemble (FlatMatrix<double> & elmat,
			 LocalHeap & lh) const;
  virtual void Assemble (FlatMatrix<Complex> & elmat,
			 LocalHeap & lh) const;
  */

  virtual void GetDofNrs (ARRAY<int> & dnums) const;

  virtual void Update ();
};



///
template <class SCAL>
class S_BEMElement : public BEMElement
{
public:
  typedef SCAL TSCAL;
  S_BEMElement (const FESpace & afes,
		const FESpace & afes_neumann);


  virtual void Assemble (FlatMatrix<SCAL> & elmat,
			 LocalHeap & lh) const;


  virtual void AssembleIntegralOps (FlatMatrix<SCAL> & singlelayer, 
				    FlatMatrix<SCAL> & doublelayer, 
				    FlatMatrix<SCAL> & hypersingular, 
				    FlatMatrix<SCAL> & duality,
				    LocalHeap & lh) const = 0;
};



template <class SCAL>
class S_LaplaceBEMElement : public S_BEMElement<SCAL>
{
public:
  typedef SCAL TSCAL;
  S_LaplaceBEMElement (const FESpace & afes,
		       const FESpace & afes_neumann)
    : S_BEMElement<SCAL> (afes, afes_neumann) { ; }


  virtual void AssembleIntegralOps (FlatMatrix<SCAL> & singlelayer, 
				    FlatMatrix<SCAL> & doublelayer, 
				    FlatMatrix<SCAL> & hypersingular, 
				    FlatMatrix<SCAL> & duality,
				    LocalHeap & lh) const;


  virtual void AssembleElementMatrix (const NodalFiniteElement & feld1, 
				      const NodalFiniteElement & feld2, 
				      const NodalFiniteElement & feln1, 
				      const NodalFiniteElement & feln2, 
				      const ElementTransformation & eltrans1, 
				      const ElementTransformation & eltrans2, 
				      FlatMatrix<SCAL> & elsingle, 
				      FlatMatrix<SCAL> & eldouble, 
				      FlatMatrix<SCAL> & elduality, 
				      FlatMatrix<SCAL> & elhyper, 
				      LocalHeap & locheap) const;

  double Fundamental (double * x) const;
  double DFundamental (double * x, double * nx) const;
  double DDFundamental (double * x, double * nx, double * ny) const;
};



///
class Helmholtz_BEMElement : public S_BEMElement<Complex>
{
  double k0;
  S_LaplaceBEMElement<Complex> laplace;
public:
  typedef Complex TSCAL;
  Helmholtz_BEMElement (const FESpace & afes,
			const FESpace & afes_neumann,
			double ak);
  

  virtual void Update ();
  
  virtual void AssembleIntegralOps (FlatMatrix<Complex> & singlelayer, 
				    FlatMatrix<Complex> & doublelayer, 
				    FlatMatrix<Complex> & hypersingular, 
				    FlatMatrix<Complex> & duality,
				    LocalHeap & lh) const;
  
  virtual void AssembleElementMatrix (const NodalFiniteElement & feld1, 
				      const NodalFiniteElement & feld2, 
				      const NodalFiniteElement & feln1, 
				      const NodalFiniteElement & feln2, 
				      const ElementTransformation & eltrans1, 
				      const ElementTransformation & eltrans2, 
				      FlatMatrix<Complex> & elsingle, 
				      FlatMatrix<Complex> & eldouble, 
				      FlatMatrix<Complex> & elduality, 
				      FlatMatrix<Complex> & elhyper, 
				      LocalHeap & locheap) const;  

  
  Complex Fundamental (double * x) const;
  Complex DFundamental (double * x, double * nx) const;
  Complex DDFundamental (double * x, double * nx, double * ny) const;
};


#endif