#ifndef FILE_EQUILIBRIUM
#define FILE_EQUILIBRIUM

/*********************************************************************/
/* File:   equilibrium.hpp                                           */
/* Author: Joachim Schoeberl                                         */
/* Date:   5. Jul. 2001                                              */
/*********************************************************************/

/** 
   realizations of equilibrium methods
   ref: R. Stenberg (1988) "A Family of Mixed Finite Elements for
   the Elasticity Prolbem", Numer. Math. 53, pp 513-538
*/


class EquilibriumIntegrator : public BilinearFormIntegrator
{
public:
  ///
  EquilibriumIntegrator ();
  ///
  virtual ~EquilibriumIntegrator ();

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

  ///
  virtual void ComputeInternalVariables (const FiniteElement & fel, 
					 const ElementTransformation & eltrans, 
					 const FlatVector<> & uext, 
					 FlatVector<> & sigma, 
					 FlatVector<> & uint,
					 LocalHeap & locheap) const;

  ///
  virtual void ComputePointValues (const FiniteElement & fel,
				   const ElementTransformation & eltrans, 
				   const IntegrationPoint & sip,
				   const FlatVector<> & sigma, 
				   const FlatVector<> & uint,
				   FlatVector<> & psigma, FlatVector<> & pu, 
				   LocalHeap & locheap ) const = 0;

  ///
  virtual void GetInternalDofs (const FiniteElement & fel,
				BitArray & internal) const = 0;
  
  ///
  virtual void ComputeMatrices (const FiniteElement & fel, 
				const ElementTransformation & eltrans, 
				ngbla::FlatMatrix<> & mata,
				ngbla::FlatMatrix<> & matb,
				ngbla::FlatMatrix<> & matc,
				LocalHeap & locheap) const = 0;
  
  ///
  virtual void ComputeVectors (const FiniteElement & fel, 
			       const ElementTransformation & eltrans, 
			       ngbla::FlatVector<> & vecfsigma,
			       ngbla::FlatVector<> & vecf,
			       LocalHeap & locheap) const = 0;

  virtual bool BoundaryForm() const
  { return 0; }
};


///
template <int DIM>
class ElasticityEquilibriumIntegrator : public EquilibriumIntegrator
{
  typedef CompositeFiniteElement<HDivFiniteElement<DIM>, NodalFiniteElement,NodalFiniteElement> TFEL;

  ///
  CoefficientFunction * coeffe;
  ///
  CoefficientFunction * coeffnu;
  ///
  CoefficientFunction * fx, *fy, *fz;

public:
  ///
  ElasticityEquilibriumIntegrator (CoefficientFunction * acoeffe,
				   CoefficientFunction * acoeffnu,
				   CoefficientFunction * afx = 0,
				   CoefficientFunction * afy = 0,
				   CoefficientFunction * afz = 0);
  ///
  virtual ~ElasticityEquilibriumIntegrator ();
 

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



  ///
  virtual void GetInternalDofs (const FiniteElement & fel,
				BitArray & internal) const;
  
  ///
  virtual void ComputeMatrices (const FiniteElement & fel, 
				const ElementTransformation & eltrans, 
				FlatMatrix<> & mata,
				FlatMatrix<> & matb,
				FlatMatrix<> & matc,
				LocalHeap & locheap ) const;
  
  ///
  virtual void ComputeVectors (const FiniteElement & fel, 
			       const ElementTransformation & eltrans, 
			       FlatVector<> & vsigma,
			       FlatVector<> & vu,
			       LocalHeap & locheap ) const;

  ///
  virtual void ComputePointValues (const FiniteElement & fel,
				   const ElementTransformation & eltrans, 
				   const IntegrationPoint & ip,
				   const FlatVector<> & sigma, const FlatVector<> & uint,
				   FlatVector<> & psigma, FlatVector<> & pu, 
				   LocalHeap & locheap ) const;

  ///
  void GetElement (ELEMENT_TYPE eltype,
		   HDivFiniteElement<DIM> *& sigmael,
		   FiniteElement *& uel,
		   FiniteElement *& gammael) const;
};



///
class ElasticityEquilibriumIntegratorStab : public EquilibriumIntegrator
{
  ///
  CoefficientFunction * coeffe;
  ///
  CoefficientFunction * coeffnu;
  ///
  CoefficientFunction * fx, *fy, *fz;

public:
  ///
  ElasticityEquilibriumIntegratorStab (CoefficientFunction * acoeffe,
				       CoefficientFunction * acoeffnu,
				       CoefficientFunction * afx,
				       CoefficientFunction * afy,
				       CoefficientFunction * afz);
  ///
  virtual ~ElasticityEquilibriumIntegratorStab ();

  ///
  virtual void GetInternalDofs (const FiniteElement & fel,
				BitArray & internal) const;
  
  ///
  virtual void ComputeMatrices (const FiniteElement & fel, 
				const ElementTransformation & eltrans, 
				FlatMatrix<> & mata,
				FlatMatrix<> & matb,
				FlatMatrix<> & matc,
				LocalHeap & locheap ) const;
  
  ///
  virtual void ComputeVectors (const FiniteElement & fel, 
			       const ElementTransformation & eltrans, 
			       FlatVector<> & vsigma,
			       FlatVector<> & vu,
			       LocalHeap & locheap ) const;

  ///
  virtual void ComputePointValues (const FiniteElement & fel,
				   const ElementTransformation & eltrans, 
				   const IntegrationPoint & ip,
				   const FlatVector<> & sigma, 
				   const FlatVector<> & uint,
				   FlatVector<> & psigma, FlatVector<> & pu, 
				   LocalHeap & locheap ) const;

  ///
  void GetElement (ELEMENT_TYPE eltype,
		   HDivFiniteElement<2> *& sigmael,
		   FiniteElement *& uel,
		   FiniteElement *& gammael) const;
};






///
class EquilibriumForceIntegrator : public LinearFormIntegrator
{
  ///
  EquilibriumIntegrator * gequ;

public:
  ///
  EquilibriumForceIntegrator (EquilibriumIntegrator * agequ);
  ///
  virtual ~EquilibriumForceIntegrator ();
  ///
  virtual void
  AssembleElementVector (const FiniteElement & fel, 
			 const ElementTransformation & eltrans, 
			 FlatVector<double> & elvec,
			 LocalHeap & locheap ) const;
};










///
class ScalarEquilibriumIntegrator : public EquilibriumIntegrator
{
  ///
  CoefficientFunction * lam;
  ///
  CoefficientFunction * rho;
  ///
  CoefficientFunction * f;

public:
  ///
  ScalarEquilibriumIntegrator (CoefficientFunction * alam,
			       CoefficientFunction * arho,
			       CoefficientFunction * af);
  ///
  virtual ~ScalarEquilibriumIntegrator ();

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


  ///
  virtual void GetInternalDofs (const FiniteElement & fel,
				BitArray & internal) const;
  
  ///
  virtual void ComputeMatrices (const FiniteElement & fel, 
				const ElementTransformation & eltrans, 
				FlatMatrix<> & mata,
				FlatMatrix<> & matb,
				FlatMatrix<> & matc,
				LocalHeap & locheap ) const;
  
  ///
  virtual void ComputeVectors (const FiniteElement & fel, 
			       const ElementTransformation & eltrans, 
			       FlatVector<> & vsigma,
			       FlatVector<> & vu,
			       LocalHeap & locheap ) const;

  ///
  virtual void ComputePointValues (const FiniteElement & fel,
				   const ElementTransformation & eltrans, 
				   const IntegrationPoint & ip,
				   const FlatVector<> & sigma, 
				   const FlatVector<> & uint,
				   FlatVector<> & psigma, FlatVector<> & pu, 
				   LocalHeap & locheap ) const;

  ///
  void GetElement (ELEMENT_TYPE eltype,
		   HDivFiniteElement<2> *& sigmael,
		   FiniteElement *& uel) const;
};






///
class EquilibriumSourceIntegrator : public LinearFormIntegrator
{
  ///
  EquilibriumIntegrator * gequ;

public:
  ///
  EquilibriumSourceIntegrator (EquilibriumIntegrator * agequ);
  ///
  virtual ~EquilibriumSourceIntegrator ();

  ///
  static Integrator * Create (ARRAY<CoefficientFunction*> & coeffs)
  {
    return new EquilibriumSourceIntegrator 
      ( new ScalarEquilibriumIntegrator (coeffs[0], coeffs[1], coeffs[2]));
  }

  virtual bool BoundaryForm () const { return 0; }
  ///
  virtual void
  AssembleElementVector (const FiniteElement & fel, 
			 const ElementTransformation & eltrans, 
			 FlatVector<double> & elvec,
			 LocalHeap & locheap ) const;
};










#ifdef NONE

// returns the matrix for stresses, reduce other variables
/// by penalty
class EquilibriumIntegratorxxx : public BilinearFormIntegrator
{
  ///
  int peers;
  ///
  CoefficientFunction * coeffe;
  ///
  CoefficientFunction * coeffnu;

public:
  ///
  EquilibriumIntegrator (CoefficientFunction * acoeffe,
			 CoefficientFunction * acoeffnu);
  ///
  virtual ~EquilibriumIntegrator ();
  ///
  virtual BaseMatrix &  
  AssembleElementMatrix (const FiniteElement & fel, 
			 const ElementTransformation & eltrans, 
			 LocalHeap & locheap ) const;

  ///
  virtual void ComputeDisplacement (const ElementTransformation & eltrans, 
				    const FlatVector<> & uf, FlatVector<> & sigma, FlatVector<> & disp,
				    CoefficientFunction * afx,
				    CoefficientFunction * afy,
				    LocalHeap & locheap ) const;

  ///
  virtual void ComputePointValues (const FiniteElement & fel,
				   const ElementTransformation & eltrans, 
				   const IntegrationPoint & ip,
				   const FlatVector<> & sigma, const FlatVector<> & disp,
				   FlatVector<> & psigma, FlatVector<> & pdisp,
				   LocalHeap & locheap ) const;

  ///
  void ComputeMatrices (const FiniteElement & fel, 
			const ElementTransformation & eltrans, 
			FlatMatrix<> & mata,
			FlatMatrix<> & matb1,
			FlatMatrix<> & matb2,
			FlatMatrix<> & matb3,
			LocalHeap & locheap ) const;

  ///
  void ComputeVector (const ElementTransformation & eltrans, 
		      CoefficientFunction * acoeff1,
		      CoefficientFunction * acoeff2,
		      FlatVector<> & vecf,
		      LocalHeap & locheap ) const;

  ///
  void GetElement (ELEMENT_TYPE eltype,
		   HDivFiniteElement<DIM> *& sigmael,
		   FiniteElement *& uel,
		   FiniteElement *& gammael) const;
};

///
class EquilibriumForceIntegrator : public LinearFormIntegrator
{
  ///
  CoefficientFunction * fx;
  ///
  CoefficientFunction * fy;
  ///
  EquilibriumIntegrator equ1;

public:
  ///
  EquilibriumForceIntegrator (CoefficientFunction * acoeffe,
			      CoefficientFunction * acoeffnu);
  ///
  virtual ~EquilibriumForceIntegrator ();
  ///
  virtual BaseVector &
  AssembleElementVector (const FiniteElement & fel, 
			 const ElementTransformation & eltrans, 
			 LocalHeap & locheap ) const;
};
#endif

#endif