#ifndef FILE_INTEGRATOR
#define FILE_INTEGRATOR

/*********************************************************************/
/* File:   integrator.hpp                                            */
/* Author: Joachim Schoeberl                                         */
/* Date:   25. Mar. 2000                                             */
/*********************************************************************/


template <int M>
void FastMat (int n, double * ba, double * pb, double * pc);


/*
   bilinear-form and linear-form integrators
*/


/**
   Base class for linear-form and bilinear-form integrators.
   Provides integration order, restriction to subdomains
 */
class Integrator
{
protected:
  /// define only on some sub-domains
  BitArray definedon;

  /// if >= 0, use this order of integration
  int integration_order;

  /// if >= 0, use this order of integration for all terms
  static int common_integration_order;

public:
  ///
  Integrator() throw ();

  ///
  virtual ~Integrator();

  ///
  virtual bool BoundaryForm () const = 0;

  ///
  bool DefinedOn (int mat) const;
  
  ///
  void SetDefinedOn (const BitArray & adefinedon);

  ///
  static void SetCommonIntegrationOrder (int cio)
  { 
    common_integration_order = cio; 
  }

  ///
  void SetIntegrationOrder (int io)
  {
    integration_order = io; 
  }

  /// should be pure virtual
  virtual int DimElement () const { return -1; }

  ///
  virtual int DimSpace () const { return -1; }

  ///
  virtual string Name () const;
};




/**
   A BilinearFormIntegrator computes the element matrices. Different
   equations are provided by derived classes. A Integrator can be defined
   in the domain or at the boundary.
 */
class BilinearFormIntegrator : public Integrator
{
public:
  typedef double TSCAL;
  ///
  BilinearFormIntegrator () throw ();
  ///
  virtual ~BilinearFormIntegrator ();

  /// components of flux
  virtual int DimFlux () const { return -1; }

  /**
     Assembles element matrix.
     Result is in elmat, memory is allocated by functions on LocalHeap.
   */
  virtual void
  AssembleElementMatrix (const FiniteElement & fel,
			 const ElementTransformation & eltrans, 
			 FlatMatrix<double> & elmat,
			 LocalHeap & locheap) const = 0;

  virtual void
  AssembleElementMatrix (const FiniteElement & fel,
			 const ElementTransformation & eltrans, 
			 FlatMatrix<Complex> & elmat,
			 LocalHeap & locheap) const;

  virtual void
  AssembleElementMatrixDiag (const FiniteElement & fel,
			     const ElementTransformation & eltrans, 
			     FlatVector<double> & diag,
			     LocalHeap & lh) const;

  virtual void
  AssembleLinearizedElementMatrix (const FiniteElement & fel, 
				   const ElementTransformation & eltrans, 
				   FlatVector<double> & elveclin,
				   FlatMatrix<double> & elmat,
				   LocalHeap & locheap) const;

  virtual void
  AssembleLinearizedElementMatrix (const FiniteElement & fel, 
				   const ElementTransformation & eltrans, 
				   FlatVector<Complex> & elveclin,
				   FlatMatrix<Complex> & elmat,
				   LocalHeap & locheap) const;

  virtual void 
  ApplyElementMatrix (const FiniteElement & fel, 
		      const ElementTransformation & eltrans, 
		      const FlatVector<double> & elx, 
		      FlatVector<double> & ely,
		      LocalHeap & locheap) const;

  virtual void 
  ApplyElementMatrix (const FiniteElement & fel, 
		      const ElementTransformation & eltrans, 
		      const FlatVector<Complex> & elx, 
		      FlatVector<Complex> & ely,
		      LocalHeap & locheap) const;

  virtual void 
  ApplyLinearizedElementMatrix (const FiniteElement & fel, 
				const ElementTransformation & eltrans, 
				const FlatVector<double> & ellin,
				const FlatVector<double> & elx, 
				FlatVector<double> & ely,
				LocalHeap & locheap) const;

  virtual void 
  ApplyLinearizedElementMatrix (const FiniteElement & fel, 
				const ElementTransformation & eltrans, 
				const FlatVector<Complex> & ellin, 
				const FlatVector<Complex> & elx, 
				FlatVector<Complex> & ely,
				LocalHeap & locheap) const;



  virtual double Energy (const FiniteElement & fel, 
			 const ElementTransformation & eltrans, 
			 const FlatVector<double> & elx, 
			 LocalHeap & locheap) const;

  virtual double Energy (const FiniteElement & fel, 
			 const ElementTransformation & eltrans, 
			 const FlatVector<Complex> & elx, 
			 LocalHeap & locheap) const;


  ///
  virtual FlatMatrix<double> 
  AssembleMixedElementMatrix (const FiniteElement & fel1, 
			      const FiniteElement & fel2, 
			      const ElementTransformation & eltrans, 
			      LocalHeap & locheap) const;
  /*
  {
    cerr << "AssembleMixedElementMatrix called for base class" << endl;
    return FlatMatrix<TSCAL> (0,0,0);
  }
  */

  ///
  virtual void 
  ApplyMixedElementMatrix (const FiniteElement & fel1, 
			   const FiniteElement & fel2, 
			   const ElementTransformation & eltrans, 
			   const FlatVector<TSCAL> & elx, 
			   FlatVector<TSCAL> & ely,
			   LocalHeap & locheap) const;
  /*
  {
    ely = AssembleMixedElementMatrix (fel1, fel2, eltrans, locheap) * elx;
  }
  */

  virtual void
  CalcFlux (const FiniteElement & fel,
	    const ElementTransformation & eltrans,
	    const IntegrationPoint & ip,
	    const FlatVector<double> & elx, 
	    FlatVector<double> & flux,
	    bool applyd,
	    LocalHeap & lh) const;
  /*
  {
    cerr << "calcflux<double> called for class " 
	 << typeid(*this).name()
	 << endl;
  }
  */

  virtual void
  CalcFlux (const FiniteElement & fel,
	    const ElementTransformation & eltrans,
	    const IntegrationPoint & ip,
	    const FlatVector<Complex> & elx, 
	    FlatVector<Complex> & flux,
	    bool applyd,
	    LocalHeap & lh) const;
  /*
  {
    cerr << "calcflux<Complex> called for base class " 
	 << typeid(*this).name()
	 << endl;
  }
  */

  virtual void
  CalcFlux (const FiniteElement & fel,
	    const BaseSpecificIntegrationPoint & bsip,
	    const FlatVector<double> & elx, 
	    FlatVector<double> & flux,
	    bool applyd,
	    LocalHeap & lh) const;
    /*
  {
    cerr << "calcflux<double> called for class " 
	 << typeid(*this).name()
	 << endl;
  }
    */

  virtual void
  CalcFlux (const FiniteElement & fel,
	    const BaseSpecificIntegrationPoint & bsip,
	    const FlatVector<Complex> & elx, 
	    FlatVector<Complex> & flux,
	    bool applyd,
	    LocalHeap & lh) const;
  /*
  {
    cerr << "calcflux<Complex> called for base class " 
	 << typeid(*this).name()
	 << endl;
  }
  */


  virtual void
  ApplyBTrans (const FiniteElement & fel,
	       // const ElementTransformation & eltrans,
	       // const IntegrationPoint & ip,
	       const BaseSpecificIntegrationPoint & bsip,
	       const FlatVector<double> & elx, 
	       FlatVector<double> & ely,
	       LocalHeap & lh) const;
  /*
  {
    cerr << "ApplyBTrans<double> called for class " 
	 << typeid(*this).name()
	 << endl;
  }
  */

  virtual void
  ApplyBTrans (const FiniteElement & fel,
	       // const ElementTransformation & eltrans,
	       // const IntegrationPoint & ip,
	       const BaseSpecificIntegrationPoint & bsip,
	       const FlatVector<Complex> & elx, 
	       FlatVector<Complex> & ely,
	       LocalHeap & lh) const;
  /*
  {
    cerr << "ApplyBTrans<Complex> called for class " 
	 << typeid(*this).name()
	 << endl;
  }
  */
};




/**
   Element assembling.
   Assembling for bilinear-forms of type $\int (B v) : D (B u) dx$. \\
   Template argument DiffOp provides differential operator, i.e. B matrix,
   (e.g. gradient, strain operator, curl,...) \\
   DmatOp provides d-matrix (e.g. diagonal, anisotropic, plane stress, ...) \\
   FEL is element type to assemble matrix for (NodalFiniteElement, 
   HCurlFiniteElement, FE_Trig1, ...)
 */
template <class DIFFOP, class DMATOP, class FEL = FiniteElement>
class T_BDBIntegrator : public virtual BilinearFormIntegrator
{
protected:
  DMATOP dmatop;
public:
  typedef double TSCAL;  // necessary ?

  enum { DIM_SPACE   = DIFFOP::DIM_SPACE };
  enum { DIM_ELEMENT = DIFFOP::DIM_ELEMENT };
  enum { DIM_DMAT    = DIFFOP::DIM_DMAT };
  enum { DIM         = DIFFOP::DIM };

  /*
  // same in base-class
  virtual void
  AssembleLinearizedElementMatrix (const FiniteElement & fel, 
				   const ElementTransformation & eltrans, 
				   FlatVector<double> & elveclin,
				   FlatMatrix<double> & elmat,
				   LocalHeap & locheap) const 
  {
    AssembleElementMatrix (fel, eltrans, elmat, locheap);
  }
  */

  ///
  T_BDBIntegrator (const DMATOP & admat)
    : dmatop(admat)
  { ; }

  ///
  virtual ~T_BDBIntegrator ()
  { ; }

  ///
  virtual bool BoundaryForm () const
  { return int(DIM_SPACE) > int(DIM_ELEMENT); }

  virtual int DimElement () const
  { return DIM_ELEMENT; }

  virtual int DimSpace () const
  { return DIM_SPACE; }

  virtual int DimFlux () const
  { return DIM_DMAT; }

  DMATOP & DMat () { return dmatop; }
  const DMATOP & DMat () const { return dmatop; }




#ifdef TEXT_BOOK_VERSION

  virtual void
  AssembleElementMatrix (const FiniteElement & bfel,
			 const ElementTransformation & eltrans, 
			 FlatMatrix<double> & 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, double> bmat (ndof * DIM, locheap);
	FlatMatrixFixHeight<DIM_DMAT, double> dbmat (ndof * DIM, locheap);
	Mat<DIM_DMAT,DIM_DMAT> dmat;

	const IntegrationRule & ir = GetIntegrationRule (fel);

	void * heapp = locheap.GetPointer();

	for (int i = 0 ; i < ir.GetNIP(); i++)
	  {
	    SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> 
	      sip(ir[i], eltrans, locheap);

	    DIFFOP::GenerateMatrix (fel, sip, bmat, locheap);
	    dmatop.GenerateMatrix (fel, sip, dmat, locheap);
	    double fac =  
	      fabs (sip.GetJacobiDet()) * sip.IP().Weight();

	    dbmat = fac * (dmat * bmat);
	    
	    if (DMATOP::SYMMETRIC)
	      elmat += Symmetric (Trans (dbmat) * bmat);
	    else
	      elmat += Trans (bmat) * dbmat; 

	    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;
      }
  }
#endif


#define BLOCK_VERSION
#ifdef BLOCK_VERSION
  ///
  virtual void
  AssembleElementMatrix (const FiniteElement & bfel,
			 const ElementTransformation & eltrans, 
			 FlatMatrix<double> & 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;
	
	enum { BLOCK = 24 / DIM_DMAT + 1 };

	void * heapp = locheap.GetPointer();

	FlatMatrixFixHeight<DIM_DMAT, double> bmat (ndof * DIM, locheap);
	FlatMatrixFixHeight<DIM_DMAT, double> dbmat (ndof * DIM, locheap);

	FlatMatrixFixHeight<DIM_DMAT*BLOCK, double> bbmat (ndof * DIM, locheap);
	FlatMatrixFixHeight<DIM_DMAT*BLOCK, double> bdbmat (ndof * DIM, locheap);
	Mat<DIM_DMAT,DIM_DMAT> dmat;

	const IntegrationRule & ir = GetIntegrationRule (fel);

	int i = 0;
	for (int i1 = 0; i1 < ir.GetNIP() / BLOCK; i1++)
	  {
	    for (int i2 = 0; i2 < BLOCK; i++, i2++)
	      {
		void * heapp2 = locheap.GetPointer();

		SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> 
		  sip(ir[i], eltrans, locheap);

		DIFFOP::GenerateMatrix (fel, sip, bmat, locheap);
		dmatop.GenerateMatrix (fel, sip, dmat, locheap);
		double fac =  
		  fabs (sip.GetJacobiDet()) * sip.IP().Weight();
		
		dbmat = fac * (dmat * bmat);

		for (int l = 0; l < ndof*DIM; l++)
		  for (int k = 0; k < DIM_DMAT; k++)
		    {
		      bbmat(i2*DIM_DMAT+k, l) = bmat(k,l);
		      bdbmat(i2*DIM_DMAT+k, l) = dbmat(k,l);
		    }

		locheap.CleanUp (heapp2);
	      }


#ifndef USE_BLAS
	    /*
	    if (DMATOP::SYMMETRIC)
	      elmat += Symmetric (Trans (bdbmat) * bbmat);
	    else
	      elmat += Trans (bbmat) * bdbmat; 
	    */

	    if (DMATOP::SYMMETRIC)
	      FastMat<DIM_DMAT*BLOCK> (elmat.Height(), &bdbmat(0,0), &bbmat(0,0), &elmat(0,0));
	    else
	      elmat += Trans (bbmat) * bdbmat; 

#else

	    cblas_dgemm (CblasRowMajor, CblasNoTrans, CblasTrans,
			 ndof*DIM, ndof*DIM, DIM_DMAT*BLOCK, 1.0,
			 &bdbmat(0,0), DIM_DMAT*BLOCK, 
			 &bbmat(0,0), DIM_DMAT*BLOCK,
			 1.0, &elmat(0,0), ndof*DIM);

#endif
	  } 

	for ( ; i < ir.GetNIP(); i++)
	  {
	    void * heapp2 = locheap.GetPointer();

	    SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> 
	      sip(ir[i], eltrans, locheap);

	    DIFFOP::GenerateMatrix (fel, sip, bmat, locheap);
	    dmatop.GenerateMatrix (fel, sip, dmat, locheap);
	    double fac =  
	      fabs (sip.GetJacobiDet()) * sip.IP().Weight();

	    dbmat = fac * (dmat * bmat);

	    /*
	    if (DMATOP::SYMMETRIC)
	      elmat += Symmetric (Trans (dbmat) * bmat);
	    else
	      elmat += Trans (bmat) * dbmat; 
	    */

	    if (DMATOP::SYMMETRIC)
	      FastMat<DIM_DMAT> (elmat.Height(), &dbmat(0,0), &bmat(0,0), &elmat(0,0));
	    else
	      elmat += Trans (bmat) * dbmat;

	    locheap.CleanUp (heapp2);
	  } 

	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;
      }
  }


#else




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

    try
      {
	// cout << "start assemble, free = " << locheap.Available() << endl;

	const FEL & fel = dynamic_cast<const FEL&> (bfel);
	int ndof = fel.GetNDof();

	const IntegrationRule & ir = GetIntegrationRule (fel);

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

	FlatMatrix<double> bbmat (ndof * DIM, DIM_DMAT*ir.GetNIP(), locheap);
	FlatMatrix<double> bdbmat (ndof * DIM, DIM_DMAT*ir.GetNIP(), locheap);

	void * heapp = locheap.GetPointer();
	for (int i = 0; i < ir.GetNIP(); i++)
	  {
	    // cout << "assemble, ip loop, free = " << locheap.Available() << endl;

	    SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> 
	      sip(ir[i], eltrans, locheap);

	    DIFFOP::GenerateMatrix (fel, sip, bmat, locheap);
	    dmatop.GenerateMatrix (fel, sip, dmat, locheap);
	    double fac =  
	      fabs (sip.GetJacobiDet()) * sip.IP().Weight();

	    dmat *= fac;
	    dbmat = dmat * bmat;
	    // dbmat = fac * (dmat * bmat);

	    for (int l = 0; l < ndof*DIM; l++)
	      for (int k = 0; k < DIM_DMAT; k++)
		{
		  bbmat(l, i*DIM_DMAT+k) = bmat(k,l);
		  bdbmat(l, i*DIM_DMAT+k) = dbmat(k,l);
		}
	    locheap.CleanUp (heapp);
	  }


	if (DMATOP::SYMMETRIC)
	  elmat += Symmetric ( bdbmat * Trans (bbmat));
	else
	  elmat += bbmat * Trans (bdbmat); 

	/*
	cblas_dgemm (CblasRowMajor, CblasNoTrans, CblasTrans,
		     ndof*DIM, ndof*DIM, bbmat.Width(), 1.0,
		     &bdbmat(0,0), bbmat.Width(), 
		     &bbmat(0,0), bbmat.Width(),
		     1.0, &elmat(0,0), ndof*DIM);
	*/
      } 
    

    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;
      }
  }

#endif







  /*
  ///
  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> 
	      sip(ir[i], eltrans, locheap);

	    DIFFOP::GenerateMatrix (fel, sip, bmat, locheap);
	    
	    dmatop.GenerateMatrix (fel, sip, dmat, locheap);
	   

	    double fac =  
	      fabs (sip.GetJacobiDet()) * sip.IP().Weight();

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

	    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 void
  AssembleElementMatrixDiag (const FiniteElement & bfel,
			     const ElementTransformation & eltrans, 
			     FlatVector<double> & diag,
			     LocalHeap & locheap) const
  {
    try
      {
	const FEL & fel = dynamic_cast<const FEL&> (bfel);
	int ndof = fel.GetNDof();

	diag.AssignMemory (ndof*DIM, locheap);
	diag = 0;
	
	FlatMatrixFixHeight<DIM_DMAT, double> bmat (ndof * DIM, locheap);
	Mat<DIM_DMAT,DIM_DMAT> dmat;

	const IntegrationRule & ir = GetIntegrationRule (fel);

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

	    DIFFOP::GenerateMatrix (fel, sip, bmat, locheap);

	    dmatop.GenerateMatrix (fel, sip, dmat, locheap);

	    double fac =  
	      fabs (sip.GetJacobiDet()) * sip.IP().Weight();

	    for (int j = 0; j < diag.Size(); j++)
	      {
		double sum = 0;
		for (int k = 0; k < DIM_DMAT; k++)
		  for (int l = 0; l < DIM_DMAT; l++)
		    sum += dmat(k,l) * bmat(k,j) * bmat(l,j);
		diag(j) += fac * sum;
	      }

	    locheap.CleanUp (heapp);
	  } 
      }

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


  
  ///
  virtual void 
  ApplyElementMatrix (const FiniteElement & bfel, 
		      const ElementTransformation & eltrans, 
		      const FlatVector<double> & elx, 
		      FlatVector<double> & ely,
		      LocalHeap & lh) const
  {
    const FEL & fel = dynamic_cast<const FEL&> (bfel);

    int ndof = fel.GetNDof ();
    
    ely = 0;
    
    Vec<DIM_DMAT,double> hv1;
    Vec<DIM_DMAT,double> hv2;

    FlatVector<double> hely (ndof*DIM, lh);

    const IntegrationRule & ir = GetIntegrationRule (fel);

    void * heapp = lh.GetPointer();
    for (int i = 0; i < ir.GetNIP(); i++)
      {
	const IntegrationPoint & ip = ir.GetIP(i);
	
	SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> sip (ir[i], eltrans, lh);

	DIFFOP::Apply (fel, sip, elx, hv1, lh);
	dmatop.Apply (fel, sip, hv1, hv2, lh);
	DIFFOP::ApplyTrans (fel, sip, hv2, hely, lh);

	double fac = fabs (sip.GetJacobiDet()) * ip.Weight();
	ely += fac * hely;
	lh.CleanUp (heapp);
      }     
  }

  ///
  virtual void 
  ApplyElementMatrix (const FiniteElement & bfel, 
		      const ElementTransformation & eltrans, 
		      const FlatVector<Complex> & elx, 
		      FlatVector<Complex> & ely,
		      LocalHeap & locheap) const
  {
    const FEL & fel = static_cast<const FEL&> (bfel);
    int ndof = fel.GetNDof ();
    
    ely = 0;
    
    Vec<DIM_DMAT,Complex> hv1;
    Vec<DIM_DMAT,Complex> hv2;
    
    FlatVector<Complex> hely (ndof*DIM, locheap);
    const IntegrationRule & ir = GetIntegrationRule (fel);


    for (int i = 0; i < ir.GetNIP(); i++)
      {
	SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> sip (ir[i], eltrans, locheap);

	DIFFOP::Apply (fel, sip, elx, hv1, locheap);
	dmatop.Apply (fel, sip, hv1, hv2, locheap);
	DIFFOP::ApplyTrans (fel, sip, hv2, hely, locheap);

	double fac = fabs (sip.GetJacobiDet()) * sip.IP().Weight();
	ely += fac * hely;
      }     
  }




  ///
  virtual void 
  ApplyMixedElementMatrix (const FiniteElement & bfel1, 
			   const FiniteElement & bfel2, 
			   const ElementTransformation & eltrans, 
			   const FlatVector<double> & elx, 
			   FlatVector<double> & ely,
			   LocalHeap & lh) const
  {
    const FEL & fel1 = static_cast<const FEL&> (bfel1);
    const FEL & fel2 = static_cast<const FEL&> (bfel2);
    // int ndof1 = fel1.GetNDof ();
    int ndof2 = fel2.GetNDof ();
    
    ely = 0;
    
    Vec<DIM_DMAT,double> hv1;
    Vec<DIM_DMAT,double> hv2;
    
    FlatVector<double> hely (ndof2*DIM, lh);

    const IntegrationRule & ir = GetIntegrationRule (fel2);

    void * heapp = lh.GetPointer();
    for (int i = 0; i < ir.GetNIP(); i++)
      {
	const IntegrationPoint & ip = ir.GetIP(i);
	
	SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> sip (ir[i], eltrans, lh);

	DIFFOP::Apply (fel1, sip, elx, hv1, lh);
	dmatop.Apply (fel1, sip, hv1, hv2, lh);
	DIFFOP::ApplyTrans (fel2, sip, hv2, hely, lh);

	double fac = fabs (sip.GetJacobiDet()) * ip.Weight();
	ely += fac * hely;
	lh.CleanUp (heapp);
      }     
  }





  const IntegrationRule & GetIntegrationRule (const FiniteElement & fel) const
  {
    int order = 2 * fel.Order();

    ELEMENT_TYPE et = fel.ElementType();

    if (et == ET_TET || et == ET_TRIG || et == ET_SEGM)
      order -= 2 * DIFFOP::DIFFORDER;

    if (common_integration_order >= 0)
      order = common_integration_order;

    if (integration_order >= 0)
      order = integration_order;

    return GetIntegrationRules().SelectIntegrationRule (et, order);
  }





  virtual void
  CalcFlux (const FiniteElement & bfel,
	    const ElementTransformation & eltrans,
	    const IntegrationPoint & ip,
	    const FlatVector<double> & elx, 
	    FlatVector<double> & flux,
	    bool applyd,
	    LocalHeap & lh) const
  {
    SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> sip (ip, eltrans, lh);
    CalcFlux (bfel, sip, elx, flux, applyd, lh);
  }


  virtual void
  CalcFlux (const FiniteElement & bfel,
	    const ElementTransformation & eltrans,
	    const IntegrationPoint & ip,
	    const FlatVector<Complex> & elx, 
	    FlatVector<Complex> & flux,
	    bool applyd,
	    LocalHeap & lh) const
  {
    SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> sip (ip, eltrans, lh);
    CalcFlux (bfel, sip, elx, flux, applyd, lh);
  }


  virtual void
  CalcFlux (const FiniteElement & bfel,
	    const BaseSpecificIntegrationPoint & bsip,
	    const FlatVector<double> & elx, 
	    FlatVector<double> & flux,
	    bool applyd,
	    LocalHeap & lh) const
  {
    const FEL & fel = static_cast<const FEL&> (bfel);
    const SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> sip =
      static_cast<const SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bsip);

    flux.AssignMemory (DIM_DMAT, lh);
    if (applyd)
      {
	Vec<DIM_DMAT,double> hv1;
	DIFFOP::Apply (fel, sip, elx, hv1, lh);
	dmatop.Apply (fel, sip, hv1, flux, lh);
      }
    else
      {
	DIFFOP::Apply (fel, sip, elx, flux, lh);
      }
  }
  


  virtual void
  CalcFlux (const FiniteElement & bfel,
	    const BaseSpecificIntegrationPoint & bsip,
	    const FlatVector<Complex> & elx, 
	    FlatVector<Complex> & flux,
	    bool applyd,
	    LocalHeap & lh) const
  {
    const FEL & fel = static_cast<const FEL&> (bfel);
    const SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> sip =
      static_cast<const SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bsip);
    flux.AssignMemory (DIM_DMAT, lh);
    if (applyd)
      {
	Vec<DIM_DMAT,Complex> hv1;
	DIFFOP::Apply (fel, sip, elx, hv1, lh);
	dmatop.Apply (fel, sip, hv1, flux, lh);
      }
    else
      {
	DIFFOP::Apply (fel, sip, elx, flux, lh);
      }
  }
  







  virtual void
  ApplyBTrans (const FiniteElement & bfel,
	       const BaseSpecificIntegrationPoint & bsip,
	       const FlatVector<double> & elx, 
	       FlatVector<double> & ely,
	       LocalHeap & lh) const
  {
    const FEL & fel = dynamic_cast<const FEL&> (bfel);
    const SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> sip =
      static_cast<const SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bsip);
    DIFFOP::ApplyTrans (fel, sip, elx, ely, lh);
    // ely *= sip.IP().Weight() * fabs (sip.GetJacobiDet());
  }
  

  virtual void
  ApplyBTrans (const FiniteElement & bfel,
	       const BaseSpecificIntegrationPoint & bsip,
	       const FlatVector<Complex> & elx, 
	       FlatVector<Complex> & ely,
	       LocalHeap & lh) const
  {
    const FEL & fel = dynamic_cast<const FEL&> (bfel);
    const SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> sip =
      static_cast<const SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bsip);
    DIFFOP::ApplyTrans (fel, sip, elx, ely, lh);
    // ely *= sip.IP().Weight() * fabs (sip.GetJacobiDet());
  }
  


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

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





template <class DIFFOP, class DMATOP, class FEL = FiniteElement>
class T_NonlinearBDBIntegrator : public T_BDBIntegrator<DIFFOP, DMATOP, FEL>
{
  enum { DIM_SPACE   = DIFFOP::DIM_SPACE };
  enum { DIM_ELEMENT = DIFFOP::DIM_ELEMENT };
  enum { DIM_DMAT    = DIFFOP::DIM_DMAT };
  enum { DIM         = DIFFOP::DIM };

public:
  ///
  T_NonlinearBDBIntegrator (const DMATOP & admat)
    : T_BDBIntegrator<DIFFOP,DMATOP,FEL> (admat)
  { ; }

  ///
  virtual ~T_NonlinearBDBIntegrator ()
  { ; }



  virtual void
  AssembleLinearizedElementMatrix (const FiniteElement & bfel, 
				   const ElementTransformation & eltrans, 
				   FlatVector<double> & elveclin,
				   FlatMatrix<double> & 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, double> bmat (ndof * DIM, locheap);
	FlatMatrixFixHeight<DIM_DMAT, double> dbmat (ndof * DIM, locheap);
	Vec<DIM_DMAT,double> hvlin;

	Mat<DIM_DMAT,DIM_DMAT> dmat;

	const IntegrationRule & ir = GetIntegrationRule (fel);

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

	    DIFFOP::Apply (fel, sip, elveclin, hvlin, locheap);

	    DIFFOP::GenerateMatrix (fel, sip, bmat, locheap);

	    this->dmatop . GenerateLinearizedMatrix (fel, sip, hvlin, dmat, locheap);

	    double fac =  
	      fabs (sip.GetJacobiDet()) * sip.IP().Weight();

	    dbmat = dmat * bmat;
	    elmat += fac * (Trans (bmat) * dbmat); 
	    
	    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 void
  AssembleLinearizedElementMatrix (const FiniteElement & bfel, 
				   const ElementTransformation & eltrans, 
				   FlatVector<Complex> & elveclin,
				   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);
	Vec<DIM_DMAT,Complex> hvlin;

	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> 
	      sip(ir[i], eltrans, locheap);

	    DIFFOP::Apply (fel, sip, elveclin, hvlin, locheap);

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

	    double fac =  
	      fabs (sip.GetJacobiDet()) * sip.IP().Weight();

	    dbmat = dmat * bmat;
	    elmat += fac * (Trans (bmat) * dbmat); 
	    
	    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 void 
  ApplyLinearizedElementMatrix (const FiniteElement & bfel, 
				const ElementTransformation & eltrans, 
				const FlatVector<double> & ellin, 
				const FlatVector<double> & elx, 
				FlatVector<double> & ely,
				LocalHeap & locheap) const
  {
    const FEL & fel = static_cast<const FEL&> (bfel);
    int ndof = fel.GetNDof ();
    
    ely = 0;
    
    Vec<DIM_DMAT,double> hvlin;
    Vec<DIM_DMAT,double> hvx;
    Vec<DIM_DMAT,double> hvy;
    
    FlatVector<double> hely (ndof*DIM, locheap);

    /*    
    int order = IntegrationOrder (fel);
    const IntegrationRule & ir = 
      GetIntegrationRules().SelectIntegrationRule (fel.ElementType(), order);
    */
    const IntegrationRule & ir = GetIntegrationRule (fel);


    for (int i = 0; i < ir.GetNIP(); i++)
      {
	const IntegrationPoint & ip = ir.GetIP(i);
	
	SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> sip (ir[i], eltrans, locheap);

	DIFFOP::Apply (fel, sip, ellin, hvlin, locheap);
	DIFFOP::Apply (fel, sip, elx, hvx, locheap);
	this->dmatop.ApplyLinearized (fel, sip, hvlin, hvx, hvy, locheap);
	DIFFOP::ApplyTrans (fel, sip, hvy, hely, locheap);

	double fac = fabs (sip.GetJacobiDet()) * ip.Weight();
	ely += fac * hely;
      }     
  }



 ///
  virtual void 
  ApplyLinearizedElementMatrix (const FiniteElement & bfel, 
				const ElementTransformation & eltrans, 
				const FlatVector<Complex> & ellin, 
				const FlatVector<Complex> & elx, 
				FlatVector<Complex> & ely,
				LocalHeap & locheap) const
  {
    const FEL & fel = static_cast<const FEL&> (bfel);
    int ndof = fel.GetNDof ();
    
    ely = 0;
    
    Vec<DIM_DMAT,Complex> hvlin;
    Vec<DIM_DMAT,Complex> hvx;
    Vec<DIM_DMAT,Complex> hvy;
    
    FlatVector<Complex> hely (ndof*DIM, locheap);
    
    /*
    int order = IntegrationOrder (fel);
    const IntegrationRule & ir = 
      GetIntegrationRules().SelectIntegrationRule (fel.ElementType(), order);
    */
    const IntegrationRule & ir = GetIntegrationRule (fel);



    for (int i = 0; i < ir.GetNIP(); i++)
      {
	SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> sip (ir[i], eltrans, locheap);

	DIFFOP::Apply (fel, sip, ellin, hvlin, locheap);
	DIFFOP::Apply (fel, sip, elx, hvx, locheap);
	this->dmatop.ApplyLinearized (fel, sip, hvlin, hvx, hvy, locheap);
	DIFFOP::ApplyTrans (fel, sip, hvy, hely, locheap);
	
	double fac = fabs (sip.GetJacobiDet()) * sip.IP().Weight();
	ely += fac * hely;
      }     
  }

};




// returns coef * Identity-matrix
class RegularizationBilinearFormIntegrator : public BilinearFormIntegrator
{
  CoefficientFunction * coef;
public:
  RegularizationBilinearFormIntegrator (CoefficientFunction * acoef)
    : coef(acoef)
  { ; }

  virtual bool BoundaryForm () const { return 0; }
  virtual int DimFlux () const { return 1; }
  virtual int DimElement () const { return 3; }
  virtual int DimSpace () const { return 3; }

  virtual void
  AssembleElementMatrix (const FiniteElement & fel,
			 const ElementTransformation & eltrans, 
			 FlatMatrix<double> & elmat,
			 LocalHeap & locheap) const
  {
    int ndof = fel.GetNDof();

    elmat.AssignMemory (ndof, ndof, locheap);
    elmat = 0;
    SpecificIntegrationPoint<3,3> sip (GetIntegrationRules().SelectIntegrationRule (fel.ElementType(), 0)[0], 
				       eltrans, locheap);
    double val = Evaluate (*coef, sip);
    elmat = 0;
    for (int i = 0; i < ndof; i++)
      elmat(i,i) = val;
  }

  virtual string Name () const { return "Regularization"; }
};





class BlockBilinearFormIntegrator : public BilinearFormIntegrator
{
  BilinearFormIntegrator & bfi;
  int dim;
  int comp;
public:
  BlockBilinearFormIntegrator (BilinearFormIntegrator & abfi, int adim, int acomp);
  BlockBilinearFormIntegrator (BilinearFormIntegrator & abfi, int adim);
  virtual ~BlockBilinearFormIntegrator ();

  virtual bool BoundaryForm () const
  { return bfi.BoundaryForm(); }

  virtual int DimFlux () const 
  { return (comp == -1) ? dim * bfi.DimFlux() : bfi.DimFlux(); }

  const BilinearFormIntegrator & Block () const { return bfi; }

  virtual void
  AssembleElementMatrix (const FiniteElement & bfel, 
			 const ElementTransformation & eltrans, 
			 FlatMatrix<double> & elmat,
			 LocalHeap & locheap) const;

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

  virtual void 
  ApplyElementMatrix (const FiniteElement & bfel, 
		      const ElementTransformation & eltrans, 
		      const FlatVector<double> & elx, 
		      FlatVector<double> & ely,
		      LocalHeap & locheap) const;

  virtual void 
  ApplyElementMatrix (const FiniteElement & bfel, 
		      const ElementTransformation & eltrans, 
		      const FlatVector<Complex> & elx, 
		      FlatVector<Complex> & ely,
		      LocalHeap & locheap) const;

  virtual void
  AssembleLinearizedElementMatrix (const FiniteElement & bfel,
				   const ElementTransformation & eltrans,
				   FlatVector<double> & elveclin,
				   FlatMatrix<double> & elmat,
				   LocalHeap & locheap) const;
  virtual void
  AssembleLinearizedElementMatrix (const FiniteElement & bfel, 
				   const ElementTransformation & eltrans, 
				   FlatVector<Complex> & elveclin,
				   FlatMatrix<Complex> & elmat,
				   LocalHeap & locheap) const;


  virtual void
  CalcFlux (const FiniteElement & fel,
	    const ElementTransformation & eltrans,
	    const IntegrationPoint & ip,
	    const FlatVector<double> & elx, 
	    FlatVector<double> & flux,
	    bool applyd,
	    LocalHeap & lh) const;

  virtual void
  CalcFlux (const FiniteElement & fel,
	    const ElementTransformation & eltrans,
	    const IntegrationPoint & ip,
	    const FlatVector<Complex> & elx, 
	    FlatVector<Complex> & flux,
	    bool applyd,
	    LocalHeap & lh) const;



  virtual void
  CalcFlux (const FiniteElement & fel,
	    const BaseSpecificIntegrationPoint & bsip,
	    const FlatVector<double> & elx, 
	    FlatVector<double> & flux,
	    bool applyd,
	    LocalHeap & lh) const;

  virtual void
  CalcFlux (const FiniteElement & fel,
	    const BaseSpecificIntegrationPoint & bsip,
	    const FlatVector<Complex> & elx, 
	    FlatVector<Complex> & flux,
	    bool applyd,
	    LocalHeap & lh) const;



  virtual void
  ApplyBTrans (const FiniteElement & bfel,
	       const BaseSpecificIntegrationPoint & bsip,
	       const FlatVector<double> & elx, 
	       FlatVector<double> & ely,
	       LocalHeap & lh) const;

  virtual void
  ApplyBTrans (const FiniteElement & bfel,
	       const BaseSpecificIntegrationPoint & bsip,
	       const FlatVector<Complex> & elx, 
	       FlatVector<Complex> & ely,
	       LocalHeap & lh) const;

  virtual double Energy (const FiniteElement & fel, 
			 const ElementTransformation & eltrans, 
			 const FlatVector<double> & elx, 
			 LocalHeap & locheap) const;

  virtual string Name () const;
};




class NormalBilinearFormIntegrator : public BilinearFormIntegrator
{
  const BilinearFormIntegrator & bfi;
public:
  NormalBilinearFormIntegrator (const BilinearFormIntegrator & abfi);
  virtual bool BoundaryForm () const
  { return bfi.BoundaryForm(); }

  virtual void
  AssembleElementMatrix (const FiniteElement & bfel, 
			 const ElementTransformation & eltrans, 
			 FlatMatrix<double> & elmat,
			 LocalHeap & locheap) const;
};







class ComplexBilinearFormIntegrator : public BilinearFormIntegrator
{
  const BilinearFormIntegrator & bfi;
  Complex factor;
public:
  ComplexBilinearFormIntegrator (const BilinearFormIntegrator & abfi,
				 Complex afactor);

  virtual bool BoundaryForm () const
  { return bfi.BoundaryForm(); }

  virtual int DimFlux () const
  { return bfi.DimFlux(); }
  virtual int DimElement () const
  { return bfi.DimElement(); }
  virtual int DimSpace () const
  { return bfi.DimSpace(); }



  virtual void
  AssembleElementMatrix (const FiniteElement & fel, 
			 const ElementTransformation & eltrans, 
			 FlatMatrix<double> & elmat,
			 LocalHeap & locheap) const;

  virtual void
  AssembleElementMatrix (const FiniteElement & fel, 
			 const ElementTransformation & eltrans, 
			 FlatMatrix<Complex> & elmat,
			 LocalHeap & locheap) const;

  virtual void 
  ApplyElementMatrix (const FiniteElement & fel, 
		      const ElementTransformation & eltrans, 
		      const FlatVector<Complex> & elx, 
		      FlatVector<Complex> & ely,
		      LocalHeap & locheap) const;

  virtual void
  CalcFlux (const FiniteElement & fel,
	    const ElementTransformation & eltrans,
	    const IntegrationPoint & ip,
	    const FlatVector<Complex> & elx, 
	    FlatVector<Complex> & flux,
	    bool applyd,
	    LocalHeap & lh) const;

  virtual string Name () const;
};






class CompoundBilinearFormIntegrator : public BilinearFormIntegrator
{
  const BilinearFormIntegrator & bfi;
  int comp;
public:
  CompoundBilinearFormIntegrator (const BilinearFormIntegrator & abfi, int acomp);
  
  virtual bool BoundaryForm () const
  { return bfi.BoundaryForm(); }

  virtual int DimFlux () const 
  { return bfi.DimFlux(); }
  virtual int DimElement () const
  { return bfi.DimElement(); }
  virtual int DimSpace () const
  { return bfi.DimSpace(); }


  virtual void
  AssembleElementMatrix (const FiniteElement & bfel, 
			 const ElementTransformation & eltrans, 
			 FlatMatrix<double> & elmat,
			 LocalHeap & locheap) const;

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


  virtual void
  AssembleLinearizedElementMatrix (const FiniteElement & fel, 
				   const ElementTransformation & eltrans, 
				   FlatVector<double> & elveclin,
				   FlatMatrix<double> & elmat,
				   LocalHeap & locheap) const;

  virtual void
  AssembleLinearizedElementMatrix (const FiniteElement & fel, 
				   const ElementTransformation & eltrans, 
				   FlatVector<Complex> & elveclin,
				   FlatMatrix<Complex> & elmat,
				   LocalHeap & locheap) const;

  virtual void
  ApplyElementMatrix (const FiniteElement & bfel, 
		      const ElementTransformation & eltrans, 
		      const FlatVector<double> & elx,
		      FlatVector<double> & ely,
		      LocalHeap & locheap) const;

  virtual void
  ApplyElementMatrix (const FiniteElement & bfel, 
		      const ElementTransformation & eltrans, 
		      const FlatVector<Complex> & elx,
		      FlatVector<Complex> & ely,
		      LocalHeap & locheap) const;

  virtual void
  ApplyLinearizedElementMatrix (const FiniteElement & bfel, 
				const ElementTransformation & eltrans, 
				const FlatVector<double> & ellin,
				const FlatVector<double> & elx,
				FlatVector<double> & ely,
				LocalHeap & locheap) const;

  virtual void
  ApplyLinearizedElementMatrix (const FiniteElement & bfel, 
				const ElementTransformation & eltrans, 
				const FlatVector<Complex> & ellin,
				const FlatVector<Complex> & elx,
				FlatVector<Complex> & ely,
				LocalHeap & locheap) const;

  virtual void
  CalcFlux (const FiniteElement & bfel,
	    const ElementTransformation & eltrans,
	    const IntegrationPoint & ip,
	    const FlatVector<double> & elx, 
	    FlatVector<double> & flux,
	    bool applyd,
	    LocalHeap & lh) const;

  virtual void
  CalcFlux (const FiniteElement & bfel,
	    const ElementTransformation & eltrans,
	    const IntegrationPoint & ip,
	    const FlatVector<Complex> & elx, 
	    FlatVector<Complex> & flux,
	    bool applyd,
	    LocalHeap & lh) const;

  virtual string Name () const;
};











/**
   Mixed element matrix assembling.
   Assembling for mixed bilinear-forms of type $\int (B_2 v) : D (B_1 u) dx$.
 */

template <int SPACEDIM = 2, class SCAL = double>
class B1DB2Integrator : public BilinearFormIntegrator
{
public:
  ///
  B1DB2Integrator ();
  ///
  virtual ~B1DB2Integrator ();
  ///
  /*
  virtual BaseMatrix &
  AssembleMixedElementMatrix (const FiniteElement & fel1, 
			      const FiniteElement & fel2, 
			      const ElementTransformation & eltrans, 
			      LocalHeap & locheap) const;
  */
  ///
  virtual int GetDimension1 () const = 0;
  ///
  virtual int GetDimension2 () const = 0;
  ///
  virtual int GetDimensionD1 () const = 0;
  ///
  virtual int GetDimensionD2 () const = 0;
  ///
  virtual int DiffOrder1 () const { return 0; }
  ///
  virtual int DiffOrder2 () const { return 0; }

  ///
  virtual void GenerateB1Matrix (const FiniteElement & fel,
				 const SpecificIntegrationPoint<> & sip,
				 FlatMatrix<> & bmat,
				 LocalHeap & locheap) const = 0;

  ///
  virtual void GenerateB2Matrix (const FiniteElement & fel,
				 const SpecificIntegrationPoint<> & sip,
				 FlatMatrix<> & bmat,
				 LocalHeap & locheap) const = 0;

  ///
  virtual void GenerateDMatrix (const FiniteElement & fel,
				const SpecificIntegrationPoint<> & sip,
				FlatMatrix<SCAL> & dmat,
				LocalHeap & locheap) const = 0;

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








/*
   Mixed element matrix assembling.
   Element assembling for Bilinear-forms of type $\int (B_2 v) : D (B_1 u) dx$.

template <class SCAL = double>
class B1DB2BoundaryIntegrator : public S_BilinearFormIntegrator<SCAL>
{
public:
  ///
  B1DB2BoundaryIntegrator ();
  ///
  virtual ~B1DB2BoundaryIntegrator ();

  ///

  virtual BaseMatrix & AssembleMixedElementMatrix
  (const FiniteElement & fel1, 
   const FiniteElement & fel2, 
   const ElementTransformation & eltrans, 
   LocalHeap & locheap) const;
  
  ///
  virtual void 
  ApplyMixedElementMatrix (const FiniteElement & fel1, 
			   const FiniteElement & fel2, 
			   const ElementTransformation & eltrans, 
			   const BaseVector & elx, 
			   BaseVector & ely,
			   LocalHeap & locheap) const;
  ///
  virtual bool BoundaryForm () const { return 1; }

  ///
  virtual int GetDimension1 () const = 0;
  ///
  virtual int GetDimension2 () const = 0;
  ///
  virtual int GetDimensionD1 () const = 0;
  ///
  virtual int GetDimensionD2 () const = 0;
  ///
  virtual int DiffOrder1 () const { return 0; }
  ///
  virtual int DiffOrder2 () const { return 0; }

  ///
  virtual void GenerateB1Matrix (const FiniteElement & fel,
				 const SpecificIntegrationPoint<> & sip,
				 FlatMatrix<> & bmat,
				 LocalHeap & locheap) const = 0;

  ///
  virtual void GenerateB2Matrix (const FiniteElement & fel,
				 const SpecificIntegrationPoint<> & sip,
				 FlatMatrix<> & bmat,
				 LocalHeap & locheap) const = 0;

  ///
  virtual void GenerateDMatrix (const FiniteElement & fel,
				const SpecificIntegrationPoint<> & sip,
				FlatMatrix<SCAL> & dmat,
				LocalHeap & locheap) const = 0;


  ///
  virtual void ApplyB1Matrix (const FiniteElement & fel,
			      const SpecificIntegrationPoint<> & sip,
			      const FlatVector<SCAL> & x, 
			      FlatVector<SCAL> & y,
			      LocalHeap & locheap) const;

  ///
  virtual void ApplyB2MatrixTrans (const FiniteElement & fel,
				   const SpecificIntegrationPoint<> & sip,
				   const FlatVector<SCAL> & x, 
				   FlatVector<SCAL> & y,
				   LocalHeap & locheap) const;

  ///
  virtual void ApplyDMatrix (const FiniteElement & fel,
			     const SpecificIntegrationPoint<> & sip,
			     const FlatVector<SCAL> & x, 
			     FlatVector<SCAL> & y,
			     LocalHeap & locheap) const;


  ///
  virtual int Lumping () const
    { return 0; }
  ///
  virtual const char * Name () const { return "B1DB2 boundary integrator"; }
};
*/








template <int DIM_SPACE>
class DirichletPenaltyIntegrator : public BilinearFormIntegrator
{
  CoefficientFunction * penalty;
public:
  DirichletPenaltyIntegrator (CoefficientFunction * apenalty)
    : penalty(apenalty) { ; }

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

  virtual bool BoundaryForm () const
  { return 1; }

  virtual string Name () const
  { return "DirichletPenalty"; }

  virtual void
  AssembleElementMatrix (const FiniteElement & fel, 
			 const ElementTransformation & eltrans, 
			 FlatMatrix<double> & elmat,
			 LocalHeap & locheap) const
  {
    int ndof = fel.GetNDof();

    elmat.AssignMemory (ndof, ndof, locheap);
    
    const IntegrationRule & ir = GetIntegrationRules().SelectIntegrationRule (fel.ElementType(), 0);
    SpecificIntegrationPoint<DIM_SPACE-1,DIM_SPACE> sip (ir[0], eltrans, locheap);

    double val = Evaluate (*penalty, sip);

    elmat = 0;
    for (int i = 0; i < ndof; i++)
      elmat(i,i) = val;
  }


  virtual void 
  ApplyElementMatrix (const FiniteElement & fel, 
		      const ElementTransformation & eltrans, 
		      const FlatVector<Complex> & elx, 
		      FlatVector<Complex> & ely,
		      LocalHeap & locheap) const
  {
    const IntegrationRule & ir = GetIntegrationRules().SelectIntegrationRule (fel.ElementType(), 0);
    SpecificIntegrationPoint<DIM_SPACE-1,DIM_SPACE> sip (ir[0], eltrans, locheap);
    double val = Evaluate (*penalty, sip);

    ely = val * elx;
  }
};





















/**
   Integrator for element vector.
 */
class LinearFormIntegrator : public Integrator

{
public:
  ///
  LinearFormIntegrator () { ; }
  ///
  virtual ~LinearFormIntegrator () { ; }


  virtual void 
  AssembleElementVector (const FiniteElement & fel,
			 const ElementTransformation & eltrans, 
			 FlatVector<double> & elvec,
			 LocalHeap & locheap) const = 0;

  virtual void 
  AssembleElementVector (const FiniteElement & fel,
			 const ElementTransformation & eltrans, 
			 FlatVector<Complex> & elvec,
			 LocalHeap & locheap) const 
  {
    FlatVector<double> rvec;
    AssembleElementVector (fel, eltrans, rvec, locheap);
    elvec.AssignMemory (rvec.Size(), locheap);
    elvec = rvec;
  }

  virtual void
  AssembleElementVectorIndependent (const FiniteElement & gfel,
				    const ElementTransformation & seltrans, 
				    const IntegrationPoint & s_ip,
				    const ElementTransformation & geltrans, 
				    const IntegrationPoint & g_ip,
				    FlatVector<double> & elvec,
				    LocalHeap & locheap) const
  {
    ;
  }
  

  virtual void
  AssembleElementVectorIndependent (const FiniteElement & gfel,
				    const ElementTransformation & seltrans, 
				    const IntegrationPoint & s_ip,
				    const ElementTransformation & geltrans, 
				    const IntegrationPoint & g_ip,
				    FlatVector<Complex> & elvec,
				    LocalHeap & locheap) const
  {
    FlatVector<double> rvec;
    AssembleElementVectorIndependent (gfel, seltrans, s_ip, geltrans, g_ip,
				      rvec, locheap);
    elvec.AssignMemory (rvec.Size(), locheap);
    elvec = rvec;
  }


};


























/**
   Element assembling.
   Assembling for linear-forms of type $\int D (B u) dx$.
 */
template <class DIFFOP, class DVecOp, class FEL = FiniteElement>
class T_BIntegrator : public LinearFormIntegrator
{
  DVecOp dvecop;
public:
  enum { DIM_SPACE = DIFFOP::DIM_SPACE };
  enum { DIM_ELEMENT = DIFFOP::DIM_ELEMENT };
  enum { DIM_DMAT = DIFFOP::DIM_DMAT };
  enum { DIM = DIFFOP::DIM };
  
  ///
  T_BIntegrator (const DVecOp & advec)
    : dvecop(advec)
  { ; }

  ///
  virtual ~T_BIntegrator ()
  { ; }


  ///
  virtual bool BoundaryForm () const
  { return int(DIM_SPACE) > int(DIM_ELEMENT); }

  virtual int DimElement () const
  { return DIM_ELEMENT; }

  virtual int DimSpace () const
  { return DIM_SPACE; }



  ///
  virtual void
  AssembleElementVector (const FiniteElement & bfel,
			 const ElementTransformation & eltrans, 
			 FlatVector<double> & elvec,
			 LocalHeap & locheap) const
  {
   
    try
      {
	const FEL & fel = dynamic_cast<const FEL&> (bfel);
	int ndof = fel.GetNDof();
	
	elvec.AssignMemory (ndof * DIM, locheap);
	elvec = 0;

	FlatVector<double> hv(ndof * DIM, locheap);
	Vec<DIM_DMAT> dvec;
	
	int order = IntegrationOrder (fel);
	
	const IntegrationRule & ir = (Lumping() && fel.Order() == 1 && 0) ?
	  GetIntegrationRules().SelectLumpingIntegrationRule (fel.ElementType()) :
	  GetIntegrationRules().SelectIntegrationRule (fel.ElementType(), order);
	
	void * heapp = locheap.GetPointer();
	for (int i = 0; i < ir.GetNIP(); i++)
	  {
	    SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> sip(ir[i], eltrans, locheap);
	    dvecop.GenerateVector (fel, sip, dvec, locheap);
	    DIFFOP::ApplyTrans (fel, sip, dvec, hv, locheap);
	    double fac = fabs (sip.GetJacobiDet()) * sip.IP().Weight();

	    elvec += fac * hv;
	    locheap.CleanUp (heapp);
	  }
	 
      }

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


  virtual void
  AssembleElementVectorIndependent (const FiniteElement & gfel, 
				    const ElementTransformation & seltrans, 
				    const IntegrationPoint & s_ip,
				    const ElementTransformation & geltrans, 
				    const IntegrationPoint & g_ip,
				    FlatVector<double> & elvec,
				    LocalHeap & locheap) const
  {
    const FEL & fel = dynamic_cast<const FEL&> (gfel);
    int ndof = fel.GetNDof();
    
    elvec.AssignMemory (ndof * DIM, locheap);
    elvec = 0;

    Vec<DIM_DMAT> dvec;
	
    enum { HDIM = (DIM_SPACE > 1) ? DIM_SPACE-1 : 1 };
    SpecificIntegrationPoint<HDIM,DIM_SPACE> s_sip(s_ip, seltrans, locheap);
    SpecificIntegrationPoint<DIM_SPACE,DIM_SPACE> g_sip(g_ip, geltrans, locheap);

    dvecop.GenerateVector (fel, s_sip, dvec, locheap);
    DIFFOP::ApplyTrans (fel, g_sip, dvec, elvec, locheap);
  }  




  int IntegrationOrder (const FEL & fel) const
  {
    // int order = fel.Order()+2;    // low order case
    int order = 2*fel.Order()+1;  // high order case
    
    ELEMENT_TYPE et = fel.ElementType();

    if (et == ET_TET || et == ET_TRIG || et == ET_SEGM)
      order -= DIFFOP::DIFFORDER;

    /*    
    if (common_integration_order >= 0)
      order = common_integration_order;
    */
    if (integration_order >= 0)
      order = integration_order;

    return order;
  }

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

  ///
  virtual int Lumping () const
    { return 0; }

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







class BlockLinearFormIntegrator : public LinearFormIntegrator
{
  const LinearFormIntegrator & lfi;
  int dim;
  int comp;
public:
  BlockLinearFormIntegrator (const LinearFormIntegrator & alfi, int adim, int acomp);

  virtual bool BoundaryForm () const
  { return lfi.BoundaryForm(); }


  virtual void 
  AssembleElementVector (const FiniteElement & bfel, 
			 const ElementTransformation & eltrans, 
			 FlatVector<double> & elvec,
			 LocalHeap & locheap) const;
};





class ComplexLinearFormIntegrator : public LinearFormIntegrator
{
  const LinearFormIntegrator & lfi;
  Complex factor;
public:
  ComplexLinearFormIntegrator (const LinearFormIntegrator & alfi, 
			       Complex afactor)
    : lfi(alfi), factor(afactor)
  { ; }

  virtual bool BoundaryForm () const
  { return lfi.BoundaryForm(); }

  virtual void
  AssembleElementVector (const FiniteElement & fel, 
			 const ElementTransformation & eltrans, 
			 FlatVector<double> & elvec,
			 LocalHeap & locheap) const
  {
    throw Exception ("ComplexLinearFormIntegrator: cannot assemble double vector");
  }

  virtual void
  AssembleElementVector (const FiniteElement & fel, 
			 const ElementTransformation & eltrans, 
			 FlatVector<Complex> & elvec,
			 LocalHeap & locheap) const
  {
    FlatVector<Complex> rvec;
    lfi.AssembleElementVector (fel, eltrans, rvec, locheap);
    elvec.AssignMemory (rvec.Size(), locheap);
    elvec = factor * rvec;
  }  


  virtual void
  AssembleElementVectorIndependent (const FiniteElement & gfel, 
				    const ElementTransformation & seltrans, 
				    const IntegrationPoint & s_ip,
				    const ElementTransformation & geltrans, 
				    const IntegrationPoint & g_ip,
				    FlatVector<double> & elvec,
				    LocalHeap & locheap) const
  {
    throw Exception ("ComplexLinearFormIntegrator: cannot assemble double vector");
  }
  

  virtual void
  AssembleElementVectorIndependent (const FiniteElement & gfel, 
				    const ElementTransformation & seltrans, 
				    const IntegrationPoint & s_ip,
				    const ElementTransformation & geltrans, 
				    const IntegrationPoint & g_ip,
				    FlatVector<Complex> & elvec,
				    LocalHeap & locheap) const
  { 
    FlatVector<double> rvec;
    lfi.AssembleElementVectorIndependent (gfel, seltrans, s_ip, geltrans, g_ip,
					  rvec, locheap);
    elvec.AssignMemory (rvec.Size(), locheap);
    elvec = factor * rvec;
  }




  virtual string Name () const
  {
    return
      string ("ComplexIntegrator (") + lfi.Name() +
      string (")");
  }

};




class CompoundLinearFormIntegrator : public LinearFormIntegrator
{
  const LinearFormIntegrator & lfi;
  int comp;
public:
  CompoundLinearFormIntegrator (const LinearFormIntegrator & alfi, int acomp)
    : lfi(alfi), comp(acomp) { ; }

  virtual bool BoundaryForm () const
  { return lfi.BoundaryForm(); }


  virtual void 
  AssembleElementVector (const FiniteElement & bfel, 
			 const ElementTransformation & eltrans, 
			 FlatVector<double> & elvec,
			 LocalHeap & locheap) const
  {
    const CompoundFiniteElement & fel =
      dynamic_cast<const CompoundFiniteElement&> (bfel);

    int i;
    FlatVector<double> vec1;
    lfi.AssembleElementVector (fel[comp], eltrans, vec1, locheap);
    
    elvec.AssignMemory (fel.GetNDof(), locheap);
    elvec = 0;

    int base = 0;
    for (i = 0; i < comp; i++)
      base += fel[i].GetNDof();

    for (i = 0; i < vec1.Size(); i++)
      elvec(base+i) = vec1(i);
  }  


  virtual void 
  AssembleElementVector (const FiniteElement & bfel, 
			 const ElementTransformation & eltrans, 
			 FlatVector<Complex> & elvec,
			 LocalHeap & locheap) const
  {
    const CompoundFiniteElement & fel =
      dynamic_cast<const CompoundFiniteElement&> (bfel);

    int i;
    FlatVector<Complex> vec1;
    lfi.AssembleElementVector (fel[comp], eltrans, vec1, locheap);
    
    elvec.AssignMemory (fel.GetNDof(), locheap);
    elvec = 0;

    int base = 0;
    for (i = 0; i < comp; i++)
      base += fel[i].GetNDof();

    for (i = 0; i < vec1.Size(); i++)
      elvec(base+i) = vec1(i);
  }  


  virtual string Name () const
  {
    return
      string ("CompoundIntegrator (") + lfi.Name() +
      string (")");
  }
};










class Integrators
{
public:
  class IntegratorInfo
  {
  public:
    string name;
    int spacedim;
    int numcoeffs;
    Integrator* (*creator)(ARRAY<CoefficientFunction*> &);
    
    IntegratorInfo (const string & aname,
		    int aspacedim,
		    int anumcoffs,		      
		    Integrator* (*acreator)(ARRAY<CoefficientFunction*>&));
  };
  

  ARRAY<IntegratorInfo*> bfis;
  ARRAY<IntegratorInfo*> lfis;
  
public:
  Integrators();
  ~Integrators();  
  void AddBFIntegrator (const string & aname, int aspacedim, int anumcoeffs,
		      Integrator* (*acreator)(ARRAY<CoefficientFunction*>&));
  void AddLFIntegrator (const string & aname, int aspacedim, int anumcoeffs,
		      Integrator* (*acreator)(ARRAY<CoefficientFunction*>&));
  
  const ARRAY<IntegratorInfo*> & GetBFIs() const { return bfis; }
  const IntegratorInfo * GetBFI(const string & name, int dim) const;
  BilinearFormIntegrator * CreateBFI(const string & name, int dim,
				     ARRAY<CoefficientFunction*> & coeffs) const;

  const ARRAY<IntegratorInfo*> & GetLFIs() const { return lfis; }
  const IntegratorInfo * GetLFI(const string & name, int dim) const;

  void Print (ostream & ost) const;
};

 
extern Integrators & GetIntegrators ();

#endif