/*********************************************************************/
/* File:   postproc.cpp                                              */
/* Author: Joachim Schoeberl                                         */
/* Date:   25. Mar. 2000                                             */
/*********************************************************************/

/* 
   Postprocessing functions
*/

#include <comp.hpp>
#include <fem.hpp>

namespace ngcomp
{
  using namespace ngcomp;

  template <class SCAL>
  void CalcFlux (const MeshAccess & ma, 
		 const S_GridFunction<SCAL> & u,
		 S_GridFunction<SCAL> & flux,
		 const BilinearFormIntegrator & bli,
		 bool applyd, bool add,
		 int domain)
  {
    // cout << "calc flux, type = " << typeid(SCAL).name() << endl;

    const FESpace & fes = u.GetFESpace();
    const FESpace & fesflux = flux.GetFESpace();

    bool bound = bli.BoundaryForm();

    int ne      = bound ? ma.GetNSE() : ma.GetNE();
    int dim     = fes.GetDimension();
    int dimflux = fesflux.GetDimension();
    int i, j;

    ElementTransformation eltrans;

    ARRAY<int> dnums;
    ARRAY<int> dnumsflux;

    LocalHeap lh (1000000);
    
    if (!add) flux.GetVector() = 0.0;

    ARRAY<int> cnti(flux.GetVector().Size());
    cnti = 0;

    for (i = 0; i < ne; i++)
      {
	lh.CleanUp();

	int eldom = 
	  bound ? ma.GetSElIndex(i) : ma.GetElIndex(i);
	
	if ((domain != -1) && (domain != eldom))
	  continue;

	const FiniteElement & fel = 
	  bound ? fes.GetSFE(i, lh) : fes.GetFE (i, lh);

	//const FiniteElement & felflux = 
	//  bound ? fesflux.GetSFE(i, lh) : fesflux.GetFE (i, lh);
 	const NodalFiniteElement & felflux = 
	  dynamic_cast<const NodalFiniteElement&> 
 	  (bound ? fesflux.GetSFE(i, lh) : fesflux.GetFE (i, lh));

	if (bound)
	  {
	    ma.GetSurfaceElementTransformation (i, eltrans, lh);
	    fes.GetSDofNrs (i, dnums);
	    fesflux.GetSDofNrs (i, dnumsflux);
	  }
	else
	  {
	    ma.GetElementTransformation (i, eltrans, lh);
	    fes.GetDofNrs (i, dnums);
	    fesflux.GetDofNrs (i, dnumsflux);
	  }

	FlatVector<SCAL> elu(dnums.Size() * dim, lh);
	FlatVector<SCAL> elflux(dnumsflux.Size() * dimflux, lh);
	FlatVector<SCAL> fluxi(dimflux, lh);


	u.GetElementVector (dnums, elu);
	fes.TransformVec (i, bound, elu, TRANSFORM_SOL);

	//	if (add)
	flux.GetElementVector (dnumsflux, elflux);
	/*
	else
	  elflux = 0;
	*/
	const IntegrationRule & ir = felflux.NodalIntegrationRule();
	//const IntegrationRule & ir = 
	//  GetIntegrationRules().SelectNodalIntegrationRule(fel.ElementType(), 2*felflux.Order());

	for (j = 0; j < ir.GetNIP(); j++)
	  {
	    bli.CalcFlux (fel, eltrans, ir.GetIP(j), elu, fluxi, applyd, lh);
	    for (int k = 0; k < dimflux; k++)
	      elflux(dimflux*j+k) += fluxi(k);
	  }
	flux.SetElementVector (dnumsflux, elflux);
	//	cout << " .... elem " << i << ", val: " << elflux << endl;
	
	for (j = 0; j < dnumsflux.Size(); j++)
	  cnti[dnumsflux[j]]++;
      }

    if (!add)
      {
	FlatVector<SCAL> fluxi(dimflux, lh);
	dnumsflux.SetSize(1);
	for (i = 0; i < cnti.Size(); i++)
	  if (cnti[i])
	    {
	      dnumsflux[0] = i;
	      flux.GetElementVector (dnumsflux, fluxi);
	      fluxi /= cnti[i];
	      flux.SetElementVector (dnumsflux, fluxi);
	    }
      }
  }

  template void CalcFlux<double> (const MeshAccess & ma, 
				  const S_GridFunction<double> & bu,
				  S_GridFunction<double> & bflux,
				  const BilinearFormIntegrator & bli,
				  bool applyd, bool add, int domain);

  template void CalcFlux<Complex> (const MeshAccess & ma, 
				   const S_GridFunction<Complex> & bu,
				   S_GridFunction<Complex> & bflux,
				   const BilinearFormIntegrator & bli,
				   bool applyd, bool add, int domain);








  template <class SCAL>
  void CalcFluxProject (const MeshAccess & ma, 
			const S_GridFunction<SCAL> & u,
			S_GridFunction<SCAL> & flux,
			const BilinearFormIntegrator & bli,
			bool applyd, int domain, LocalHeap & lh)
  {
    ma.PushStatus ("Post-processing");

    const FESpace & fes = u.GetFESpace();
    const FESpace & fesflux = flux.GetFESpace();

    bool bound = bli.BoundaryForm();

    int ne      = bound ? ma.GetNSE() : ma.GetNE();
    int dim     = fes.GetDimension();
    int dimflux = fesflux.GetDimension();

    ElementTransformation eltrans;

    const BilinearFormIntegrator & fluxbli =
      bound ? (*fesflux.GetBoundaryEvaluator()) : (*fesflux.GetEvaluator());

    ARRAY<int> dnums;
    ARRAY<int> dnumsflux;

    ARRAY<int> cnti(fesflux.GetNDof());
    cnti = 0;

    flux.GetVector() = 0.0;

    void * heapp1 = lh.GetPointer();
    for (int i = 0; i < ne; i++)
      {
	ma.SetThreadPercentage ( 100.0*i / ne );

	int eldom = 
	  bound ? ma.GetSElIndex(i) : ma.GetElIndex(i);
	
	if ((domain != -1) && (domain != eldom))
	  continue;

	const FiniteElement & fel = 
	  bound ? fes.GetSFE(i, lh) : fes.GetFE (i, lh);

	const FiniteElement & felflux = 
	  bound ? fesflux.GetSFE(i, lh) : fesflux.GetFE (i, lh);

	if (bound)
	  {
	    ma.GetSurfaceElementTransformation (i, eltrans, lh);
	    fes.GetSDofNrs (i, dnums);
	    fesflux.GetSDofNrs (i, dnumsflux);
	  }
	else
	  {
	    ma.GetElementTransformation (i, eltrans, lh);
	    fes.GetDofNrs (i, dnums);
	    fesflux.GetDofNrs (i, dnumsflux);
	  }

	FlatVector<SCAL> elu(dnums.Size() * dim, lh);
	FlatVector<SCAL> elflux(dnumsflux.Size() * dimflux, lh);
	FlatVector<SCAL> elfluxi(dnumsflux.Size() * dimflux, lh);
	FlatVector<SCAL> fluxi(dimflux, lh);


	u.GetElementVector (dnums, elu);
	fes.TransformVec (i, bound, elu, TRANSFORM_SOL);

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

	elflux = 0;
	for (int j = 0; j < ir.GetNIP(); j++)
	  {
	    void * heapp2 = lh.GetPointer();

	    // bli.CalcFlux (fel, eltrans, ir.GetIP(j), elu, fluxi, applyd, lh);
	    // fluxbli.ApplyBTrans (felflux, eltrans, ir.GetIP(j), fluxi, elfluxi, lh);
	    double fac;
	    (*testout) << "bound =" << bound << endl;
	    if (!bound)
	      {
		if (fel.SpatialDim() == 2)
		  {
		    SpecificIntegrationPoint<2,2> sip (ir.GetIP(j), eltrans, lh);
		    fac = sip.IP().Weight() * fabs (sip.GetJacobiDet());
		    bli.CalcFlux (fel, sip, elu, fluxi, applyd, lh);
		    fluxbli.ApplyBTrans (felflux, sip, fluxi, elfluxi, lh);
		  }
		else
		  {
		    SpecificIntegrationPoint<3,3> sip (ir.GetIP(j), eltrans, lh);
		    fac = sip.IP().Weight() * fabs (sip.GetJacobiDet());
		    bli.CalcFlux (fel, sip, elu, fluxi, applyd, lh);
		    fluxbli.ApplyBTrans (felflux, sip, fluxi, elfluxi, lh);
		  }
	      }
	    else
	      {
		if (fel.SpatialDim() == 2)
		  {
		    SpecificIntegrationPoint<2,3> sip (ir.GetIP(j), eltrans, lh);
		    fac = sip.IP().Weight() * fabs (sip.GetJacobiDet());
		    (*testout) << "elu = " << elu << endl;
		    bli.CalcFlux (fel, sip, elu, fluxi, applyd, lh);
		    (*testout) << "fluxi = " << fluxi << endl;

		    fluxbli.ApplyBTrans (felflux, sip, fluxi, elfluxi, lh);
		  }
		else
		  {
		    SpecificIntegrationPoint<1,2> sip (ir.GetIP(j), eltrans, lh);
		    fac = sip.IP().Weight() * fabs (sip.GetJacobiDet());
		    bli.CalcFlux (fel, sip, elu, fluxi, applyd, lh);
		    fluxbli.ApplyBTrans (felflux, sip, fluxi, elfluxi, lh);
		  }
	      }


	    elflux += fac * elfluxi;
	    lh.CleanUp (heapp2);
	  }




	if (dimflux > 1)
	  {
	    FlatMatrix<SCAL> elmat(dnumsflux.Size(), lh);
	    BlockBilinearFormIntegrator const& fluxbli_ref = dynamic_cast<const BlockBilinearFormIntegrator&> (fluxbli);
	    fluxbli_ref . Block() . AssembleElementMatrix (felflux, eltrans, elmat, lh);
	    CholeskyFactors<SCAL> invelmat(elmat);

	    FlatVector<SCAL> hv1(dnumsflux.Size(), lh);
	    FlatVector<SCAL> hv2(dnumsflux.Size(), lh);
	    for (int j = 0; j < dimflux; j++)
	      {
		for (int k = 0; k < dnumsflux.Size(); k++)
		  hv1(k) = elflux(dimflux*k+j);
		invelmat.Mult (hv1, hv2);
		for (int k = 0; k < dnumsflux.Size(); k++)
		  elfluxi(dimflux*k+j) = hv2(k);
	      }
	  }
	else
	  {
	    FlatMatrix<SCAL> elmat(dnumsflux.Size(), lh);
	    fluxbli.AssembleElementMatrix (felflux, eltrans, elmat, lh);
	    CholeskyFactors<SCAL> invelmat(elmat);
	    invelmat.Mult (elflux, elfluxi);
	  }

	fesflux.TransformVec (i, bound, elfluxi, TRANSFORM_SOL);

	flux.GetElementVector (dnumsflux, elflux);
	elfluxi += elflux;
	flux.SetElementVector (dnumsflux, elfluxi);

	for (int j = 0; j < dnumsflux.Size(); j++)
	  cnti[dnumsflux[j]]++;

	lh.CleanUp(heapp1);
      }

    //    (*testout) << "cnti = " << endl << cnti << endl;
    //    (*testout) << "flux, 1 = " << endl << flux.GetVector() << endl;
    FlatVector<SCAL> fluxi(dimflux, lh);
    dnumsflux.SetSize(1);
    for (int i = 0; i < cnti.Size(); i++)
      if (cnti[i])
	{
	  dnumsflux[0] = i;
	  flux.GetElementVector (dnumsflux, fluxi);
	  fluxi /= double (cnti[i]);
	  flux.SetElementVector (dnumsflux, fluxi);
	}

    //    (*testout) << "flux, 2 = " << endl << flux.GetVector() << endl;
    ma.PopStatus ();
  }

  template void CalcFluxProject<double> (const MeshAccess & ma, 
					 const S_GridFunction<double> & bu,
					 S_GridFunction<double> & bflux,
					 const BilinearFormIntegrator & bli,
					 bool applyd, int domain, LocalHeap & lh);
  
  template void CalcFluxProject<Complex> (const MeshAccess & ma, 
					  const S_GridFunction<Complex> & bu,
					  S_GridFunction<Complex> & bflux,
					  const BilinearFormIntegrator & bli,
					  bool applyd, int domain, LocalHeap & lh);










  template <class SCAL>
  int CalcPointFlux (const MeshAccess & ma, 
		     const GridFunction & bu,
		     const FlatVector<double> & point,
		     FlatVector<SCAL> & flux,
		     const BilinearFormIntegrator & bli,
		     bool applyd,
		     LocalHeap & lh)
  {
    int elnr;
    double lami[3];
    ARRAY<int> dnums;
    ElementTransformation eltrans;

    elnr = Ng_FindElementOfPoint (const_cast<double*>(&point(0)), lami);
    
    if (!elnr) return 0;
    elnr--;

    const S_GridFunction<SCAL> & u = 
      dynamic_cast<const S_GridFunction<SCAL>&> (bu);

    const FESpace & fes = u.GetFESpace();
    const FiniteElement & fel = fes.GetFE (elnr, lh);
    ma.GetElementTransformation (elnr, eltrans, lh);
    fes.GetDofNrs (elnr, dnums);

    FlatVector<SCAL> elu(dnums.Size() * fes.GetDimension(), lh);

    u.GetElementVector (dnums, elu);
    fes.TransformVec (elnr, 0, elu, TRANSFORM_SOL);

    // (*testout) << "lami = " << lami[0] << ", " << lami[1] << ", " << lami[2];

    /*
    if (fel.ElementType() == ET_TRIG)
      {
	double hlam1 = lami[0];
	double hlam2 = lami[1];
	lami[0] = 1-hlam1-hlam2;
	lami[1] = hlam1;
      }

    if (fel.ElementType() == ET_TET)
      {
	double hlam1 = lami[0];
	double hlam2 = lami[1];
	double hlam3 = lami[2];
	lami[0] = 1-hlam1-hlam2-hlam3;
	lami[1] = hlam1;
	lami[2] = hlam2;
      }

    (*testout) << " =?= lami = " << lami[0] << ", " << lami[1] << ", " << lami[2] << endl;
    */

    IntegrationPoint ip(lami, 1);
    //    SpecificIntegrationPoint sip (ip, eltrans, &locheap);
    
    bli.CalcFlux (fel, eltrans, ip, elu, flux, applyd, lh);

    return 1;
  }
  

  template int CalcPointFlux<double> (const MeshAccess & ma, 
				      const GridFunction & u,
				      const FlatVector<double> & point,
				      FlatVector<double> & flux,
				      const BilinearFormIntegrator & bli,
				      bool applyd,
				      LocalHeap & lh);
  
  template int CalcPointFlux<Complex> (const MeshAccess & ma, 
				       const GridFunction & u,
				       const FlatVector<double> & point,
				       FlatVector<Complex> & flux,
				       const BilinearFormIntegrator & bli,
				       bool applyd,
				       LocalHeap & lh);
    



  template <class SCAL>
  void CalcError (const MeshAccess & ma, 
		  const S_GridFunction<SCAL> & u,
		  const S_GridFunction<SCAL> & flux,
		  const BilinearFormIntegrator & bli,
		  FlatVector<double> & err,
		  int domain, LocalHeap & lh)
  {
    ma.PushStatus ("Error estimator");

    const FESpace & fes = u.GetFESpace();
    const FESpace & fesflux = flux.GetFESpace();

    bool bound = bli.BoundaryForm();

    int ne      = bound ? ma.GetNSE() : ma.GetNE();
    int dim     = fes.GetDimension();
    int dimflux = fesflux.GetDimension();

    const BilinearFormIntegrator & fluxbli =
      bound ? (*fesflux.GetBoundaryEvaluator()) : (*fesflux.GetEvaluator());

    ElementTransformation eltrans;

    ARRAY<int> dnums;
    ARRAY<int> dnumsflux;

    double sum = 0;
    for (int i = 0; i < ne; i++)
      {
	cout << "\rerror element " << i << "/" << ne << flush;
	ma.SetThreadPercentage ( 100.0*i / ne );

	lh.CleanUp();

	int eldom = 
	  bound ? ma.GetSElIndex(i) : ma.GetElIndex(i);
	
	if ((domain != -1) && (domain != eldom))
	  continue;

	const FiniteElement & fel = 
	  bound ? fes.GetSFE(i, lh) : fes.GetFE (i, lh);

	const FiniteElement & felflux = 
	  (bound ? fesflux.GetSFE(i, lh) : fesflux.GetFE (i, lh));

	if (bound)
	  {
	    ma.GetSurfaceElementTransformation (i, eltrans, lh);
	    fes.GetSDofNrs (i, dnums);
	    fesflux.GetSDofNrs (i, dnumsflux);
	  }
	else
	  {
	    ma.GetElementTransformation (i, eltrans, lh);
	    fes.GetDofNrs (i, dnums);
	    fesflux.GetDofNrs (i, dnumsflux);
	  }

	FlatVector<SCAL> elu(dnums.Size() * dim, lh);
	FlatVector<SCAL> elflux(dnumsflux.Size() * dimflux, lh);
	FlatVector<SCAL> fluxi(dimflux, lh);
	FlatVector<SCAL> fluxi2(dimflux, lh);


	u.GetElementVector (dnums, elu);
	fes.TransformVec (i, bound, elu, TRANSFORM_SOL);
	flux.GetElementVector (dnumsflux, elflux);
	fesflux.TransformVec (i, bound, elflux, TRANSFORM_SOL);

	double elerr = 0;

	const IntegrationRule & ir = 
	  GetIntegrationRules().SelectIntegrationRule(felflux.ElementType(), 
						      2*felflux.Order());
	double vol = ma.ElementVolume(i);

	for (int j = 0; j < ir.GetNIP(); j++)
	  {
	    void * heapp = lh.GetPointer();

	    //	    bli.CalcFlux (fel, eltrans, ir.GetIP(j), elu, fluxi, 1, lh);
	    //	    fluxbli.CalcFlux (felflux, eltrans, ir.GetIP(j), elflux, fluxi2, 0, lh);


	    if (!bound)
	      {
		if (fel.SpatialDim() == 2)
		  {
		    SpecificIntegrationPoint<2,2> sip (ir.GetIP(j), eltrans, lh);
		    bli.CalcFlux (fel, sip, elu, fluxi, 1, lh);
		    fluxbli.CalcFlux (felflux, sip, elflux, fluxi2, 0, lh);
		  }
		else
		  {
		    SpecificIntegrationPoint<3,3> sip (ir.GetIP(j), eltrans, lh);
		    bli.CalcFlux (fel, sip, elu, fluxi, 1, lh);
		    fluxbli.CalcFlux (felflux, sip, elflux, fluxi2, 0, lh);
		  }
	      }
	    else
	      {
		if (fel.SpatialDim() == 2)
		  {
		    SpecificIntegrationPoint<2,3> sip (ir.GetIP(j), eltrans, lh);
		    bli.CalcFlux (fel, sip, elu, fluxi, 1, lh);
		    fluxbli.CalcFlux (felflux, sip, elflux, fluxi2, 0, lh);
		  }
		else
		  {
		    SpecificIntegrationPoint<1,2> sip (ir.GetIP(j), eltrans, lh);
		    bli.CalcFlux (fel, sip, elu, fluxi, 1, lh);
		    fluxbli.CalcFlux (felflux, sip, elflux, fluxi2, 0, lh);
		  }
	      }





	    //	    (*testout) << "diff: fluxi = " << fluxi << " =?= " << fluxi2 << endl;
	    fluxi -= fluxi2;

	    elerr += ir.GetIP(j).Weight() * vol * L2Norm2 (fluxi);

	    lh.CleanUp (heapp);
	  }

	err(i) += elerr;
	sum += elerr;
      }
    // cout << endl << "sum = " << sum << endl;
    ma.PopStatus ();
  }
  
  template void CalcError<double> (const MeshAccess & ma, 
				   const S_GridFunction<double> & bu,
				   const S_GridFunction<double> & bflux,
				   const BilinearFormIntegrator & bli,
				   FlatVector<double> & err,
				   int domain, LocalHeap & lh);
  
  template void CalcError<Complex> (const MeshAccess & ma, 
				    const S_GridFunction<Complex> & bu,
				    const S_GridFunction<Complex> & bflux,
				    const BilinearFormIntegrator & bli,
				    FlatVector<double> & err,
				    int domain, LocalHeap & lh);    








  template <class SCAL>
  void CalcGradient (const MeshAccess & ma,
		     const FESpace & fesh1,
		     const S_BaseVector<SCAL> & vech1,
		     const FESpace & feshcurl,
		     S_BaseVector<SCAL> & vechcurl)
  {
    cout << "CalcGrad" << endl;
    const NodalFiniteElement * h1fe;
    const HCurlFiniteElementD<2> * hcurlfe2d;
    const HCurlFiniteElementD<3> * hcurlfe3d;

    h1fe = dynamic_cast<const NodalFiniteElement*> (&fesh1.GetFE(ET_TRIG));
    hcurlfe2d = dynamic_cast<const HCurlFiniteElementD<2>*> (&feshcurl.GetFE(ET_TRIG));
    Matrix<> gradtrig(hcurlfe2d->GetNDof(), h1fe->GetNDof());
    ComputeGradientMatrix<2> (*h1fe, *hcurlfe2d, gradtrig);
    (*testout) << "gradtrig = " << gradtrig << endl;

    h1fe = dynamic_cast<const NodalFiniteElement*> (&fesh1.GetFE(ET_TET));
    hcurlfe3d = dynamic_cast<const HCurlFiniteElementD<3>*> (&feshcurl.GetFE(ET_TET));
    Matrix<> gradtet(hcurlfe3d->GetNDof(), h1fe->GetNDof());
    ComputeGradientMatrix<3> (*h1fe, *hcurlfe3d, gradtet);
    (*testout) << "gradtet = " << gradtet << endl;


    int ne = ma.GetNE();
    ARRAY<int> dnumsh1, dnumshcurl;
    LocalHeap lh(100000);
    
    for (int i = 0; i < ne; i++)
      {
	lh.CleanUp();
	fesh1.GetDofNrs (i, dnumsh1);
	feshcurl.GetDofNrs (i, dnumshcurl);

	FlatVector<SCAL> elhcurl(dnumshcurl.Size(), lh);
	FlatVector<SCAL> elh1(dnumsh1.Size(), lh);



	vech1.GetIndirect (dnumsh1, elh1);
	fesh1.TransformVec (i, 0, elh1, TRANSFORM_RHS);

	switch (fesh1.GetFE(i, lh).ElementType())
	  {
	  case ET_TRIG:
	    elhcurl = gradtrig * elh1;
	    break;
	  case ET_TET:
	    elhcurl = gradtet * elh1;
	    break;
	  }

	feshcurl.TransformVec (i, 0, elhcurl, TRANSFORM_RHS);
	vechcurl.SetIndirect (dnumshcurl, elhcurl);
      }
  }
  
  template
  void CalcGradient<double> (const MeshAccess & ma,
			     const FESpace & fesh1,
			     const S_BaseVector<double> & vech1,
			     const FESpace & feshcurl,
			     S_BaseVector<double> & vechcurl);







  
  template <class SCAL>
  void CalcGradientT (const MeshAccess & ma,
		      const FESpace & feshcurl,
		      const S_BaseVector<SCAL> & vechcurl1,
		      const FESpace & fesh1,
		      S_BaseVector<SCAL> & vech1)
  {
    cout << "CalcGrad" << endl;
    const NodalFiniteElement * h1fe;
    const HCurlFiniteElementD<2> * hcurlfe2d;
    const HCurlFiniteElementD<3> * hcurlfe3d;

    h1fe = dynamic_cast<const NodalFiniteElement*> (&fesh1.GetFE(ET_TRIG));
    hcurlfe2d = dynamic_cast<const HCurlFiniteElementD<2>*> (&feshcurl.GetFE(ET_TRIG));
    Matrix<> gradtrig(hcurlfe2d->GetNDof(), h1fe->GetNDof());
    ComputeGradientMatrix<2> (*h1fe, *hcurlfe2d, gradtrig);
    (*testout) << "gradtrig = " << gradtrig << endl;

    h1fe = dynamic_cast<const NodalFiniteElement*> (&fesh1.GetFE(ET_TET));
    hcurlfe3d = dynamic_cast<const HCurlFiniteElementD<3>*> (&feshcurl.GetFE(ET_TET));
    Matrix<> gradtet(hcurlfe3d->GetNDof(), h1fe->GetNDof());
    ComputeGradientMatrix<3> (*h1fe, *hcurlfe3d, gradtet);
    (*testout) << "gradtet = " << gradtet << endl;


    S_BaseVector<SCAL> & vechcurl =
      dynamic_cast<S_BaseVector<SCAL>&> (*vechcurl1.CreateVector());

    int ne = ma.GetNE();
    ARRAY<int> dnumsh1, dnumshcurl;
    LocalHeap lh(100000);
    
    vechcurl = vechcurl1;
    vech1.SetScalar(0); //  = SCAL(0);
    for (int i = 0; i < ne; i++)
      {
	lh.CleanUp();
	fesh1.GetDofNrs (i, dnumsh1);
	feshcurl.GetDofNrs (i, dnumshcurl);

	FlatVector<SCAL> elhcurl(dnumshcurl.Size(), lh);
	FlatVector<SCAL> elh1(dnumsh1.Size(), lh);

	vechcurl.GetIndirect (dnumshcurl, elhcurl);
	feshcurl.TransformVec (i, 0, elhcurl, TRANSFORM_RHS);

	switch (fesh1.GetFE(i, lh).ElementType())
	  {
	  case ET_TRIG:
	    elh1 = Trans (gradtrig) * elhcurl;
	    break;
	  case ET_TET:
	    elh1 = Trans (gradtet) * elhcurl;
	    break;
	  }

	fesh1.TransformVec (i, 0, elh1, TRANSFORM_RHS);
	vech1.AddIndirect (dnumsh1, elh1);

	elhcurl = 0;
	vechcurl.SetIndirect (dnumshcurl, elhcurl);
      }
  }
  
  template
  void CalcGradientT<double> (const MeshAccess & ma,
			      const FESpace & feshcurl,
			      const S_BaseVector<double> & vechcurl,
			      const FESpace & fesh1,
			      S_BaseVector<double> & vech1);


  
}