#ifndef FILE_ELEMENTTRANSFORMATION
#define FILE_ELEMENTTRANSFORMATION

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

/*
  Transformation from reference element to actual element
*/


#ifndef NETGEN_ELTRANS

/*
    Transformation from reference element to physical element.
    Uses finite element fel to describe mapping
*/
class ElementTransformation
{
  /// finite element defining transformation.
  const NodalFiniteElement * fel;


  /// matrix with points, dim * np
  FlatMatrix<> pointmat;
  ///
  bool pointmat_ownmem;

  /// normal vectors (only surfelements)
  FlatMatrix<> nvmat;

  /// element number
  int elnr;

  /// material property
  int elindex;
public:
  ///
  ElementTransformation ()
    : pointmat(0,0,0), nvmat(0,0,0), pointmat_ownmem(0)
  { ; }

  ~ElementTransformation ()
  {
    if (pointmat_ownmem) delete [] &pointmat(0,0); 
  }

  ///
  void SetElement (const NodalFiniteElement * afel, int aelnr, int aelindex)
  {
    fel = afel; 
    elnr = aelnr; 
    elindex = aelindex;
  }


  ///
  const NodalFiniteElement & GetElement () const { return *fel; }
  ///
  int GetElementNr () const { return elnr; }
  ///
  int GetElementIndex () const { return elindex; }


  ELEMENT_TYPE GetElementType () const
   { return fel->ElementType(); }
  
  ///
  template <typename T>
  void CalcJacobian (const IntegrationPoint & ip,
		     MatExpr<T> & dxdxi,
		     LocalHeap & lh) const
  {
    dxdxi = pointmat * fel->GetDShape(ip, lh);
  }

  ///
  template <typename T>
  void CalcPoint (const IntegrationPoint & ip, 
		  MatExpr<T> & point,
		  LocalHeap & lh) const
  {
    point = pointmat * fel->GetShape(ip, lh);
  }

  template <typename T1, typename T2>
  void CalcPointJacobian (const IntegrationPoint & ip,
			  T1 & point, T2 & dxdxi,
			  LocalHeap & lh) const
  {
    dxdxi = pointmat * fel->GetDShape(ip, lh);
    point = pointmat * fel->GetShape(ip, lh);
  }

  ///
  const FlatMatrix<> & PointMatrix () const { return pointmat; }
  ///
  FlatMatrix<> & PointMatrix () { return pointmat; }
  ///
  void AllocPointMatrix (int spacedim, int vertices)
  {
    if (pointmat_ownmem) delete [] &pointmat(0,0);
    pointmat.AssignMemory (spacedim, vertices, new double[spacedim*vertices]);
    pointmat_ownmem = 1;
  }
  ///
  const FlatMatrix<> & NVMatrix () const { return nvmat; }

  template <typename T>
  void CalcNormalVector (const IntegrationPoint & ip,
		         T & nv,
		         LocalHeap & lh) const
  {
    for (int i = 0; i < nvmat.Height(); i++)
      nv(i) = nvmat(i,0) ;
  }

  ///
  int SpaceDim () const
  {
    return pointmat.Height(); 
  }
};




#else




/**
    Transformation from reference element to physical element.
    Uses Netgen-meshaccess to perform transformation
*/
class ElementTransformation
{
  /// element number
  int elnr;

  /// material property
  int elindex;
  
  /// elemente in R^dim
  int dim;
  /// is it a boundary element ?
  bool boundary;
  /// geometry of element
  ELEMENT_TYPE eltype;
public:
  ///
  ElementTransformation () { ; }
 
  ///
  void SetElement (bool aboundary, int aelnr, int aelindex)
  {
    boundary = aboundary;
    elnr = aelnr;
    elindex = aelindex;
    dim = Ng_GetDimension();
  }

  ///
  int GetElementNr () const { return elnr; }
  ///
  int GetElementIndex () const { return elindex; }
  ///
  ELEMENT_TYPE GetElementType () const { return ET_TRIG; }
  
  ///
  template <typename T>
  void CalcJacobian (const IntegrationPoint & ip,
		     T & dxdxi,
		     LocalHeap & lh) const
  {
    if (boundary)
      Ng_GetSurfaceElementTransformation (elnr+1, &ip(0),
					  0, &dxdxi(0,0));
					  
    else
      Ng_GetElementTransformation (elnr+1, &ip(0),
				   0, &dxdxi(0,0));
  }

  ///
  template <typename T>
  void CalcPoint (const IntegrationPoint & ip,
		  T & point,
		  LocalHeap & lh) const
  {
    if (boundary)
      Ng_GetSurfaceElementTransformation (elnr+1, &ip(0), &point(0), 0);
    else
      Ng_GetElementTransformation (elnr+1, &ip(0), &point(0), 0);
  }



  template <typename T1, typename T2>
  void CalcPointJacobian (const IntegrationPoint & ip,
			  T1 & point, T2 & dxdxi,
			  LocalHeap & lh) const
  {
    if (boundary)
      Ng_GetSurfaceElementTransformation (elnr+1, &ip(0), &point(0), &dxdxi(0,0));
    else
      Ng_GetElementTransformation (elnr+1, &ip(0), &point(0), &dxdxi(0,0));
  }





  ///
  template <typename T>
  void CalcNormalVector (const IntegrationPoint & ip,
		         T & nv,
		         LocalHeap & lh) const
  {
    //cout << "calc normal vec" << endl;
    //cout << "bound = " << boundary << ", dim = " << dim << endl;
    if (boundary)
      {
	if (dim == 2)
	  {
	    Mat<2,1> dxdxi;
	    Ng_GetSurfaceElementTransformation (elnr+1, &ip(0), 0, &dxdxi(0));
	    //  cout << "dxdxi = " << dxdxi << endl;
	    double len = sqrt (sqr (dxdxi(0,0)) + sqr(dxdxi(1,0)));
	    //nv(0) = dxdxi(1,0) / len;  // original 
	    //nv(1) = -dxdxi(0,0) / len;
	    nv(0) = -dxdxi(1,0) / len; //SZ 
	    nv(1) = dxdxi(0,0) / len;
	    //cout << "nv = " << nv << endl;
	    // (*testout)<<"NormalVector="<<nv<<endl;
	  }
	else
	  {
	    // berechne aus dxdxi

	    Mat<3,2> dxdxi;
	    Ng_GetSurfaceElementTransformation (elnr+1, &ip(0), 0, &dxdxi(0));

	    nv(0) = dxdxi(1,0) * dxdxi(2,1) - dxdxi(2,0) * dxdxi(1,1);
	    nv(1) = dxdxi(2,0) * dxdxi(0,1) - dxdxi(0,0) * dxdxi(2,1);
	    nv(2) = dxdxi(0,0) * dxdxi(1,1) - dxdxi(1,0) * dxdxi(0,1);
	    nv /= L2Norm (nv);
	  }

      }
  }
  
  
  ///
  int SpaceDim () const
  {
    return dim;
  }

};



#endif
#endif