/*********************************************************************/
/* File:   intrule.cpp                                               */
/* Author: Joachim Schoeberl                                         */
/* Date:   26. Mar. 2000                                             */
/*********************************************************************/

/* 
   Finite Element Integration Rules
*/


#include <fem.hpp>
  
namespace ngfem
{
  using namespace ngfem;
  
  IntegrationPoint :: 
  IntegrationPoint (double api[3], double aw)
  {
    glob_nr = -1;
    nr = -1;
    pi[0] = api[0];
    pi[1] = api[1];
    pi[2] = api[2];
    weight = aw;
  }
  
  IntegrationPoint :: 
  IntegrationPoint (double p1, double p2, double p3, double aw)
  {
    glob_nr = -1;
    nr = -1;
    pi[0] = p1;
    pi[1] = p2;
    pi[2] = p3;
    weight = aw;
  }
  
  IntegrationPoint :: 
  IntegrationPoint (const FlatVector<double> & ap, double aw)
  {
    glob_nr = -1;
    nr = -1;
    pi[0] = (ap.Size() >= 1) ? ap(0) : 0;
    pi[1] = (ap.Size() >= 2) ? ap(1) : 0;
    pi[2] = (ap.Size() >= 3) ? ap(2) : 0;
    weight = aw;
  }
  
  ostream & operator<< (ostream & ost, const IntegrationPoint & ip)
  {
    ost << "IP globnr " << ip.glob_nr << " locnr = " << ip.nr << ": (" 
	<< ip.pi[0] << ", " << ip.pi[1] << ", " << ip.pi[2] 
	<< "), weight = " << ip.weight << endl;
    return ost;
  }
  





  
  template <int S, int R, typename SCAL>
  SpecificIntegrationPoint<S,R,SCAL> :: 
  SpecificIntegrationPoint (const IntegrationPoint & aip,
			    const ElementTransformation & aeltrans,
			    LocalHeap & lh)
    : DimSpecificIntegrationPoint<R> (aip, aeltrans)
  {
    // this->ipnr = aip.Nr();  

    /*    
    this->eltrans.CalcJacobian (this->ip, dxdxi, lh);
    this->eltrans.CalcPoint (this->ip, this->point, lh);
    */

    this->eltrans.CalcPointJacobian (this->ip, this->point, dxdxi, lh);

    if (S == R)
      {
	det = Det (dxdxi);
	if (det < 0 && 0)
	  {
	    cout << "A,det<0" << endl;
	    (*testout) << "A,det<0" << endl;
	  }
	dxidx = Inv (dxdxi);
      }
    else
      {
        // this->eltrans.CalcNormalVector (ip, normalvec, lh);
	// normalvec = eltrans.NVMatrix() * eltrans.GetElement().GetShape (ip, lh);

	if (R == 3)
	  {

	    det = sqrt ( sqr (dxdxi(1,0) * dxdxi(2,1) -
			      dxdxi(2,0) * dxdxi(1,1)) +
			 sqr (dxdxi(0,0) * dxdxi(2,1) -
			      dxdxi(2,0) * dxdxi(0,1)) +
			 sqr (dxdxi(1,0) * dxdxi(0,1) -
			      dxdxi(0,0) * dxdxi(1,1)) );

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

	  }
	else
	  {

	    det = sqrt ( sqr (dxdxi(0,0)) + sqr (dxdxi(1,0)));

	    normalvec(0) = -dxdxi(1,0) / det;
	    normalvec(1) = dxdxi(0,0) / det;
	  }
	


	Mat<S,S> ata, iata;
	
	ata = Trans (dxdxi) * dxdxi;
	iata = Inv (ata);
	dxidx = iata * Trans (dxdxi);
      }
  }


#ifdef OLD
  template <int S, int R, typename SCAL>
  SpecificIntegrationPoint<S,R,SCAL> :: 
  SpecificIntegrationPoint (const IntegrationRule & ir, int aipnr,
			    const class ElementTransformation & aeltrans,
			    LocalHeap & lh)
    : DimSpecificIntegrationPoint<R> (ir.GetIP(aipnr), aeltrans)
  {
    this->ipnr = aipnr;
    
    /*    
	  this->eltrans.CalcJacobian (this->ip, dxdxi, lh);
	  this->eltrans.CalcPoint (this->ip, this->point, lh);
    */
    this->eltrans.CalcPointJacobian (this->ip, this->point, dxdxi, lh);

    //cout << " point " << point << endl; 
    if (S == R)
      {
	det = Det (dxdxi);

	if (det < 0 && 0)
	  {
	    cout << "B,det<0" << endl;
	    (*testout) << "B.det = " << det << endl;
	    /*
	    (*testout) << "p = " << point << endl
		       << "dxdxi = " << FlatMatrix<> (dxdxi) << endl;
	    */
	    // (*testout) << "pmat = " << eltrans.PointMatrix() << endl;
	  }
	dxidx = Inv (dxdxi);
      }
    else
      {

	this->eltrans.CalcNormalVector (this->ip, normalvec, lh);
	// normalvec = eltrans.NVMatrix() * eltrans.GetElement().GetShape (ip, lh);
	
	if (R == 3)
	  {
	    det = sqrt ( sqr (dxdxi(1,0) * dxdxi(2,1) -
			      dxdxi(2,0) * dxdxi(1,1)) +
			 sqr (dxdxi(0,0) * dxdxi(2,1) -
			      dxdxi(2,0) * dxdxi(0,1)) +
			 sqr (dxdxi(1,0) * dxdxi(0,1) -
			      dxdxi(0,0) * dxdxi(1,1)) );
	  }
	else
	  {
	    det = sqrt ( sqr (dxdxi(0,0)) + sqr (dxdxi(1,0)));
	  }
	
	Mat<S,S> ata, iata;
	
	ata = Trans (dxdxi) * dxdxi;
	if (Det (ata) != 0)
	  iata = Inv (ata);
	else
	  iata = 0;
	dxidx = iata * Trans (dxdxi);

      }
      // (*testout) << "det = " << det << endl;
  }
#endif  

  template class SpecificIntegrationPoint<1,1>;
  template class SpecificIntegrationPoint<2,2>;
  template class SpecificIntegrationPoint<3,3>;
  template class SpecificIntegrationPoint<1,2>;
  template class SpecificIntegrationPoint<2,3>;
  
  







  IntegrationRule :: IntegrationRule ()
  {
    ;
  }
  
  IntegrationRule :: IntegrationRule (int nips)
  {
    ipts.SetAllocSize (nips);
  }
  
  IntegrationRule :: ~IntegrationRule ()
  {
    for (int i = 0; i < ipts.Size(); i++)
      delete ipts[i];
  }
  
  
  void IntegrationRule :: 
  AddIntegrationPoint (IntegrationPoint * ip)
  {
    ipts.Append (ip);
  }

  
  // computes Gaussean integration formula on (0,1) with n points
  // in: Numerical algs in C (or, was it the Fortran book ?)
  void ComputeGaussRule (int n, 
			 ARRAY<double> & xi, ARRAY<double> & wi)
  {
    //  cout << "compute gauss rule, n = " << n << endl;
    xi.SetSize (n);
    wi.SetSize (n);
    
    int m = (n+1)/2;
    double p1, p2, p3;
    double pp, z, z1;
    for (int i = 1; i <= m; i++)
      {
	z = cos ( M_PI * (i - 0.25) / (n + 0.5));
	
	while(1)
	  {
	    p1 = 1;
	    p2 = 0;
	    for (int j = 1; j <= n; j++)
	      {
		p3 = p2;
		p2 = p1;
		p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j;
	      }
	    // p1 is legendre polynomial
	    
	    pp = n * (z*p1-p2) / (z*z - 1);
	    z1 = z;
	    z = z1-p1/pp;
	    
	    if (fabs (z - z1) < 1e-14) break;
	  }
	
	xi[i-1] = 0.5 * (1 - z);
	xi[n-i] = 0.5 * (1 + z);
	wi[i-1] = wi[n-i] = 1.0 / ( (1  - z * z) * pp * pp);
      }
    //  (*testout) << "Gauss points with n = " << n << ":" << endl;
    //  for (i = 1; i <= n; i++)
    //    (*testout) << xi.Elem(i) << ",  w= " << wi.Elem(i) << endl;
  }



  //  IntegrationRules * intrules = 0;
  
  IntegrationRules :: IntegrationRules ()
  {
    int i, p;
    IntegrationRule * rule;    
    IntegrationPoint * ip;


    ARRAY<double> xi;
    ARRAY<double> wi;

    /*
    int ordersegm = 200;
    int order2d = 45;
    int tetorder = 40;
    int order3d = 30;

    int ordersegm = 50;
    int order2d = 45;
    int tetorder = 40;
    int order3d = 30;
    */

    /*
    segmentrules.SetSize(ordersegm);
    for (i = 1; i <= ordersegm; i++)
      {
	ComputeGaussRule (i/2+1, xi, wi);

	rule = new IntegrationRule (xi.Size());
	double xip[3] = { 0, 0, 0 };
	for (j = 0; j < xi.Size(); j++)
	  {
	    xip[0] = xi[j];
	    ip = new IntegrationPoint (xip, wi[j]);
	    ip->SetGlobNr (segmentpoints.Append (ip)-1);
	    rule->AddIntegrationPoint (ip);
	  }
	segmentrules[i-1] = rule;
      }
    */
    for (p = 0; p < 20; p++)
      GenerateIntegrationRule (ET_SEGM, p);

    static double qf_segm_lumping_points[][3] = 
      { 
	{ 0 },
	{ 1 },
      };
    static double qf_segm_lumping_weights[] = 
      { 0.5, 0.5 } ;
    for (i = 0; i < 2; i++)
      {
	ip = new IntegrationPoint (qf_segm_lumping_points[i],
				   qf_segm_lumping_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (segmentpoints.Append (ip)-1);
	segmentlumping.AddIntegrationPoint (ip);
      }

    trigrules.SetSize (7);

    static double qf_trig_order1_points[][3] = 
      {
	{ 1.0/3.0, 1.0/3.0 },
      };
    
    static double qf_trig_order1_weights[] = 
      {  0.5} ;


    
    rule = new IntegrationRule (1);
    ip = new IntegrationPoint (qf_trig_order1_points[0],
			       qf_trig_order1_weights[0]);
    ip->SetNr (0);
    ip->SetGlobNr (trigpoints.Append (ip)-1);
    rule->AddIntegrationPoint (ip);
    trigrules[0] = rule;


    rule = new IntegrationRule (1);
    ip = new IntegrationPoint (qf_trig_order1_points[0],
			       qf_trig_order1_weights[0]);
    ip->SetNr (0);
    ip->SetGlobNr (trigpoints.Append (ip)-1);
    rule->AddIntegrationPoint (ip);
    trigrules[1] = rule;


    static double qf_tria_order2_points[][3] = 
      {
	{ 0,   0.5 },
	{ 0.5, 0,  },
	{ 0.5, 0.5 }
      };
   
    static double qf_tria_order2_weights[] = 
      {
	1.0/6.0, 1.0/6.0 , 1.0/6.0
      };
    
    rule = new IntegrationRule (3);
    for (i = 0; i < 3; i++)
      {
	ip = new IntegrationPoint (qf_tria_order2_points[i],
				   qf_tria_order2_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (trigpoints.Append (ip)-1);
	rule->AddIntegrationPoint (ip);
      }

    trigrules[2] = rule;


    static double qf_trig_order4_points[][3] = 
      {
	{ 0.816847572980459, 0.091576213509771, },
	{ 0.091576213509771, 0.816847572980459, },
	{ 0.091576213509771, 0.091576213509771, },
	{ 0.108103018168070, 0.445948490915965, },
	{ 0.445948490915965, 0.108103018168070, },
	{ 0.445948490915965, 0.445948490915965 }
      };
    
    
    static double qf_trig_order4_weights[] = 
      {
	0.054975871827661, 0.054975871827661, 0.054975871827661,
	0.111690794839005, 0.111690794839005, 0.111690794839005
      };

    rule = new IntegrationRule (6);
    for (i = 0; i < 6; i++)
      {
	ip = new IntegrationPoint (qf_trig_order4_points[i],
				   qf_trig_order4_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (trigpoints.Append (ip)-1);
	rule->AddIntegrationPoint (ip);
      }

    trigrules[3] = rule;


    rule = new IntegrationRule (6);
    for (i = 0; i < 6; i++)
      {
	ip = new IntegrationPoint (qf_trig_order4_points[i],
				   qf_trig_order4_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (trigpoints.Append (ip)-1);
	rule->AddIntegrationPoint (ip);
      }

    trigrules[4] = rule;




    static double qf_trig_order6_points[][3] = 
      {
	{ 0.873821971016996, 0.063089014491502, },
	{ 0.063089014491502, 0.873821971016996, },
	{ 0.063089014491502, 0.063089014491502, },
	{ 0.501426509658179, 0.249286745170910, },
	{ 0.249286745170910, 0.501426509658179, },
	{ 0.249286745170910, 0.249286745170910, },

	{ 0.636502499121399, 0.310352451033785, },
	{ 0.310352451033785, 0.053145049844816, },
	{ 0.053145049844816, 0.636502499121399, },
	{ 0.636502499121399, 0.053145049844816, },
	{ 0.310352451033785, 0.636502499121399, },
	{ 0.053145049844816, 0.310352451033785, }
      };

    static double qf_trig_order6_weights[] = 
      {
	0.025422453185103, 0.025422453185103, 0.025422453185103,
	0.058393137863189, 0.058393137863189, 0.058393137863189,

	0.041425537809187, 0.041425537809187, 0.041425537809187,
	0.041425537809187, 0.041425537809187, 0.041425537809187 
      };
  
    rule = new IntegrationRule (12);
    for (i = 0; i < 12; i++)
      {
	ip = new IntegrationPoint (qf_trig_order6_points[i],
				   qf_trig_order6_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (trigpoints.Append (ip)-1);
	rule->AddIntegrationPoint (ip);
      }

    trigrules[5] = rule;

    rule = new IntegrationRule (12);
    for (i = 0; i < 12; i++)
      {
	ip = new IntegrationPoint (qf_trig_order6_points[i],
				   qf_trig_order6_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (trigpoints.Append (ip)-1);
	rule->AddIntegrationPoint (ip);
      }

    trigrules[6] = rule;


    for (p = 7; p <= 10; p++)
      GenerateIntegrationRule (ET_TRIG, p);

    for (p = 0; p <= 10; p++)
      GenerateIntegrationRule (ET_QUAD, p);


    /*

    // trig rules by degenerated tensor product rules:
    trigrules.SetSize(order2d+1);
    for (p = 7; p < trigrules.Size(); p++)
      {
	const IntegrationRule & segmrule = *segmentrules[p+1];
	IntegrationRule * trigrule = new IntegrationRule(segmrule.GetNIP()*segmrule.GetNIP());
	
	double point[3], weight;
	
	for (i = 0; i < segmrule.GetNIP(); i++)
	  for (j = 0; j < segmrule.GetNIP(); j++)
	    {
	      const IntegrationPoint & ipsegmi = segmrule.GetIP(i);
	      const IntegrationPoint & ipsegmj = segmrule.GetIP(j);

	      point[0] = ipsegmi.Point()[0];
	      point[1] = ipsegmj.Point()[0]*(1-point[0]);
	      point[2] = 0;

	      weight = ipsegmi.Weight() *
		ipsegmj.Weight() * (1-point[0]);
	      
	      ip = new IntegrationPoint (point, weight);
	      if (p <= 20)
		ip->SetGlobNr (trigpoints.Append (ip)-1);
	      trigrule->AddIntegrationPoint (ip);
	    }
	trigrules[p] = trigrule;
      }

    */












    static double qf_trig_lumping_points[][3] = 
      {
	{ 1, 0 },
	{ 0, 1, },
	{ 0, 0, }
      };

    static double qf_trig_lumping_weights[] = 
      {
	1.0/6.0, 1.0/6.0 , 1.0/6.0
      };
    
    for (i = 0; i < 3; i++)
      {
	ip = new IntegrationPoint (qf_trig_lumping_points[i],
				   qf_trig_lumping_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (trigpoints.Append (ip)-1);
	triglumping.AddIntegrationPoint (ip);
      }



    static double qf_trig_lumping2_points[][3] = 
      {
	{ 1, 0 },
	{ 0, 1, },
	{ 0, 0, },
	{ 0, 0.5 },
	{ 0.5, 0 },
	{ 0.5, 0.5 }
      };
       
    static double qf_trig_lumping2_weights[] = 
      {
	1.0/12.0, 1.0/12.0 , 1.0/12.0, 
	1.0/12.0, 1.0/12.0 , 1.0/12.0 
      };
    
    for (i = 0; i < 6; i++)
      {
	ip = new IntegrationPoint (qf_trig_lumping2_points[i],
				   qf_trig_lumping2_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (trigpoints.Append (ip)-1);
	triglumping2.AddIntegrationPoint (ip);
      }





    trignodalrules.SetSize(12);
    for (p = 1; p <= trignodalrules.Size(); p++)
      {
	rule = new IntegrationRule ((p+1)*(p+2)/2);
	int nelp = (p*p+3*p+2) / 2;
      
	int lami[3];
	double xi[3];
	int i = 0;
	for (lami[0] = 0; lami[0] <= p; lami[0]++)
	  for (lami[1] = 0; lami[1] <= p-lami[0]; lami[1]++)
	    {
	      lami[2] = p - lami[0] - lami[1];
	      
	      for (int n = 0; n < 3; n++)
		xi[n] = double(lami[n]) / p;
	    
	      ip = new IntegrationPoint (xi, 1.0 / (2.0 * nelp));
	      ip->SetNr (i); i++;
	      ip->SetGlobNr (trigpoints.Append (ip)-1);
	      rule->AddIntegrationPoint (ip);
	    }
	trignodalrules[p-1] = rule;
      }
  




    /*
    quadrules.SetSize(order2d);
    for (k = 0; k < quadrules.Size(); k++)
      {
	const IntegrationRule & segmrule = *segmentrules[k];
	IntegrationRule * quadrule = new IntegrationRule(segmrule.GetNIP()*segmrule.GetNIP());
	
	double point[3], weight;
	
	for (i = 0; i < segmrule.GetNIP(); i++)
	  for (j = 0; j < segmrule.GetNIP(); j++)
	    {
	      const IntegrationPoint & ipsegm1 = segmrule.GetIP(i);
	      const IntegrationPoint & ipsegm2 = segmrule.GetIP(j);
	      point[0] = ipsegm1.Point()[0];
	      point[1] = ipsegm2.Point()[0];
	      point[2] = 0;
	      weight = ipsegm1.Weight() * ipsegm2.Weight();
	      
	      ip = new IntegrationPoint (point, weight);
	      if (p <= 20)
		ip->SetGlobNr (quadpoints.Append (ip)-1);
	      quadrule->AddIntegrationPoint (ip);
	    }
	quadrules[k] = quadrule;
      }
    */  



    {
      static double qf_quad_lumping_points[][4] = 
	{
	  { 0, 0 },
	  { 1, 0 },
	  { 1, 1 },
	  { 0, 1 }
	};
      
      static double qf_quad_lumping_weights[] = 
	{
	  0.25, 0.25, 0.25, 0.25, 
	};
    
      for (i = 0; i < 4; i++)
	{
	  ip = new IntegrationPoint (qf_quad_lumping_points[i],
				     qf_quad_lumping_weights[i]);
	  ip->SetNr (i);
	  ip->SetGlobNr (quadpoints.Append (ip)-1);
	  quadlumping.AddIntegrationPoint (ip);
	}


    /*
      double point[3], weight;
    
      for (i = 0; i < segmentlumping.GetNIP(); i++)
	for (j = 0;j < segmentlumping.GetNIP(); j++)
	  {
	    const IntegrationPoint & ipsegm1 = segmentlumping.GetIP(i);
	    const IntegrationPoint & ipsegm2 = segmentlumping.GetIP(j);
	    point[0] = ipsegm1.Point()[0];
	    point[1] = ipsegm2.Point()[0];
	    point[2] = 0;
	    weight = ipsegm1.Weight() * ipsegm2.Weight();
	      
	    ip = new IntegrationPoint (point, weight);
	    ip->SetGlobNr (quadpoints.Append (ip)-1);
	    quadlumping.AddIntegrationPoint (ip);
	  }  
    */

    }














   


    tetrules.SetSize(5);

    static double qf_tetra_order1_points[][3] = 
      { 
	{ 0.25, 0.25, 0.25 },
      };
    
    static double qf_tetra_order1_weights[] = 
      {
	1.0/6.0
      };
    
    rule = new IntegrationRule (1);

    ip = new IntegrationPoint (qf_tetra_order1_points[0],
			       qf_tetra_order1_weights[0]);
    ip->SetNr (0);
    ip->SetGlobNr (tetpoints.Append (ip)-1);
    rule->AddIntegrationPoint (ip);

    tetrules[0] = rule;
						     
    

    static double qf_tetra_order2_points[][3] = 
      {
	{ 0.585410196624969, 0.138196601125011, 0.138196601125011 },
	{ 0.138196601125011, 0.585410196624969, 0.138196601125011 },
	{ 0.138196601125011, 0.138196601125011, 0.585410196624969 },
	{ 0.138196601125011, 0.138196601125011, 0.138196601125011 }
      };
    
    static double qf_tetra_order2_weights[] = 
      { 1.0/24.0, 1.0/24.0, 1.0/24.0, 1.0/24.0 };
    

    rule = new IntegrationRule (4);
    for (i = 0; i < 4; i++)
      {
	ip = new IntegrationPoint (qf_tetra_order2_points[i],
				   qf_tetra_order2_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (tetpoints.Append (ip)-1);
	rule->AddIntegrationPoint (ip);
      }

    tetrules[1] = rule;

    rule = new IntegrationRule (4);
    for (i = 0; i < 4; i++)
      {
	ip = new IntegrationPoint (qf_tetra_order2_points[i],
				   qf_tetra_order2_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (tetpoints.Append (ip)-1);
	rule->AddIntegrationPoint (ip);
      }

    tetrules[2] = rule;




    static double qf_tetra_order5_points[][3] = 
      {
	{ 0.454496295874350,   0.454496295874350,   0.045503704125650 },
	{ 0.454496295874350,   0.045503704125650,   0.454496295874350 },
	{ 0.454496295874350,   0.045503704125650,   0.045503704125650 },
	{ 0.045503704125650,   0.454496295874350,   0.454496295874350 },
	{ 0.045503704125650,   0.454496295874350,   0.045503704125650 },
	{ 0.045503704125650,   0.045503704125650,   0.454496295874350 },
	{ 0.310885919263301,   0.310885919263301,   0.310885919263301 },   
	{ 0.067342242210098,   0.310885919263301,   0.310885919263301 }, 
	{ 0.310885919263301,   0.067342242210098,   0.310885919263301 },  
	{ 0.310885919263301,   0.310885919263301,   0.067342242210098 },
      
	{ 0.092735250310891,   0.092735250310891,   0.092735250310891 },
	{ 0.721794249067326,   0.092735250310891,   0.092735250310891 },
	{ 0.092735250310891,   0.721794249067326,   0.092735250310891 },
	{ 0.092735250310891,   0.092735250310891,   0.721794249067326 }
      };

    static double qf_tetra_order5_weights[] = 
      {
	0.007091003462847, 0.007091003462847, 0.007091003462847,
	0.007091003462847, 0.007091003462847, 0.007091003462847,
	0.018781320953003, 0.018781320953003, 0.018781320953003, 0.018781320953003,
	0.012248840519394, 0.012248840519394, 0.012248840519394, 0.012248840519394
      };

    rule = new IntegrationRule (14);
    for (i = 0; i < 14; i++)
      {
	ip = new IntegrationPoint (qf_tetra_order5_points[i],
				   qf_tetra_order5_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (tetpoints.Append (ip)-1);
	rule->AddIntegrationPoint (ip);
      }

    tetrules[3] = rule;

    rule = new IntegrationRule (14);
    for (i = 0; i < 14; i++)
      {
	ip = new IntegrationPoint (qf_tetra_order5_points[i],
				   qf_tetra_order5_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (tetpoints.Append (ip)-1);
	rule->AddIntegrationPoint (ip);
      }

    tetrules[4] = rule;

    for (p = 6; p <= 8; p++)
      GenerateIntegrationRule (ET_TET, p);

    /*
    // tet-rules by degenerated tensor product rules:
    tetrules.SetSize(tetorder);
    for (p = 6; p <= tetrules.Size(); p++)
      {
	const IntegrationRule & segmrule = *segmentrules[p+1];
	IntegrationRule * tetrule = 
	  new IntegrationRule(segmrule.GetNIP()* segmrule.GetNIP()* segmrule.GetNIP());

	
	double point[3], weight;
	
	for (i = 0; i < segmrule.GetNIP(); i++)
	  for (j = 0; j < segmrule.GetNIP(); j++)
	    for (k = 0; k < segmrule.GetNIP(); k++)
	      {
		const IntegrationPoint & ipsegmi = segmrule.GetIP(i);
		const IntegrationPoint & ipsegmj = segmrule.GetIP(j);
		const IntegrationPoint & ipsegmk = segmrule.GetIP(k);

		point[0] = ipsegmi.Point()[0];
		point[1] = ipsegmj.Point()[0]*(1-point[0]);
		point[2] = ipsegmk.Point()[0]*(1-point[0]-point[1]);

		weight = ipsegmi.Weight() *
		  ipsegmj.Weight() * (1-point[0]) *
		  ipsegmk.Weight() * (1-point[0]-point[1]);
	      
		ip = new IntegrationPoint (point, weight);
// 		int found = 0;
// 		for (int m = 0; m < tetpoints.Size(); m++)
// 		  if (*ip == *tetpoints[m])
// 		    {
// 		      found = 1;
// 		      delete ip;
// 		      ip = tetpoints[m];
// 		    }
// 		if (!found)
		if (p <= 10)
		  ip->SetGlobNr (tetpoints.Append (ip)-1);
		tetrule->AddIntegrationPoint (ip);
	      }
	tetrules[p-1] = tetrule;
      }
*/


    tetnodalrules.SetSize(8);
    for (p = 3; p <= tetnodalrules.Size(); p++)
      {
	int nelp = (p*p*p + 6 * p * p + 11 * p + 6) / 6;
	rule = new IntegrationRule (nelp);
      
	int lami[4];
	double xi[4];
	int i = 0;
	for (lami[0] = 0; lami[0] <= p; lami[0]++)
	  for (lami[1] = 0; lami[1] <= p-lami[0]; lami[1]++)
	    for (lami[2] = 0; lami[2] <= p-lami[0]-lami[1]; lami[2]++)
	      {
		lami[3] = p - lami[0] - lami[1] - lami[2];
	      
		int n;
		for (n = 0; n < 4; n++)
		  xi[n] = double(lami[n]) / p;

		ip = new IntegrationPoint (xi, 1.0 / (6.0 * nelp));
		ip->SetNr (i); i++;
		ip->SetGlobNr (tetpoints.Append (ip)-1);
		rule->AddIntegrationPoint (ip);
	      }
	tetnodalrules[p-1] = rule;
      }

    double tet1pts[][3] = 
      { { 1, 0, 0 },
	{ 0, 1, 0 },
	{ 0, 0, 1 },
	{ 0, 0, 0 } };
      
    rule = new IntegrationRule (4);
    for (i = 0; i < 4; i++)
      {
	ip = new IntegrationPoint (tet1pts[i], 1.0 / (6.0 * 4));
	ip->SetNr (i);
	ip->SetGlobNr (tetpoints.Append (ip)-1);
	rule->AddIntegrationPoint (ip);
      }

    tetnodalrules[0] = rule;

    rule = new IntegrationRule (4);
    for (i = 0; i < 4; i++)
      {
	ip = new IntegrationPoint (tet1pts[i], 1.0 / (6.0 * 4));
	ip->SetNr (i);
	ip->SetGlobNr (tetpoints.Append (ip)-1);
	rule->AddIntegrationPoint (ip);
      }

    tetnodalrules[1] = rule;
  




    static double qf_tet_lumping_points[][3] = 
      {
	{ 0, 0, 0,},
	{ 1, 0, 0,},
	{ 0, 1, 0,},
	{ 0, 0, 1,},
      };
       
    static double qf_tet_lumping_weights[] = 
      {
	1.0/24.0, 1.0/24.0, 1.0/24.0, 1.0/24.0, 
      };
    
    for (i = 0; i < 4; i++)
      {
	ip = new IntegrationPoint (qf_tet_lumping_points[i],
				   qf_tet_lumping_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (tetpoints.Append (ip)-1);
	tetlumping.AddIntegrationPoint (ip);
      }




    for (p = 0; p <= 4; p++)
      GenerateIntegrationRule (ET_PRISM, p);



    /*
    prismrules.SetSize(order3d);
    for (k = 0; k < prismrules.Size(); k++)
      {
	const IntegrationRule & segmrule = *segmentrules[k];
	const IntegrationRule & trigrule = *trigrules[k];

	IntegrationRule * prismrule = 
	  new IntegrationRule(segmrule.GetNIP()*trigrule.GetNIP());
      
	double point[3], weight;
      
	for (i = 0; i < segmrule.GetNIP(); i++)
	  for (j = 0; j < trigrule.GetNIP(); j++)
	    {
	      const IntegrationPoint & ipsegm = segmrule.GetIP(i);
	      const IntegrationPoint & iptrig = trigrule.GetIP(j);
	      point[0] = iptrig.Point()[0];
	      point[1] = iptrig.Point()[1];
	      point[2] = ipsegm.Point()[0];
	      weight = iptrig.Weight() * ipsegm.Weight();
	      
	      ip = new IntegrationPoint (point, weight);
	      if (k <= 10)
		ip->SetGlobNr (prismpoints.Append (ip)-1);
	      prismrule->AddIntegrationPoint (ip);
	    }
	prismrules[k] = prismrule;
      }
    */

    static double qf_prismfacemidpoint[][3] = 
      {
	{ 0.5, 0, 0.5 },
	{ 0.5, 0.5, 0.5 },
	{ 0, 0.5, 0.5 },
      };

    for (i = 0; i < 3; i++)
      {
	ip = new IntegrationPoint (qf_prismfacemidpoint[i], 0.5);
	ip->SetNr (i);
	ip->SetGlobNr (prismpoints.Append (ip)-1);
	prismfacemidpoint.AddIntegrationPoint (ip);
      }



    static double qf_prism_lumping_points[][3] = 
      {
	{ 1, 0, 0 },
	{ 0, 1, 0 },
	{ 0, 0, 0 },
	{ 1, 0, 1 },
	{ 0, 1, 1 },
	{ 0, 0, 1 }
      };
       
    static double qf_prism_lumping_weights[] = 
      {
	1.0/12.0, 1.0/12.0, 1.0/12.0, 1.0/12.0, 
	1.0/12.0, 1.0/12.0, 
      };
    
    for (i = 0; i < 6; i++)
      {
	ip = new IntegrationPoint (qf_prism_lumping_points[i],
				   qf_prism_lumping_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (prismpoints.Append (ip)-1);
	prismlumping.AddIntegrationPoint (ip);
      }









    /*

    static double qf_pyramidz_order1_points_weight[][2] = 
      {
	{ 0.75,  1.0/3.0 }
      };

    static double qf_pyramidz_order3_points_weight[][2] = 
      {
	{ 0.455848155988775, 0.100785882079825 },
	{ 0.877485177344559, 0.232547451253508 }
      };
    static double qf_pyramidz_order5_points_weight[][2] = 
      {
	{ 0.294997790111502, 0.029950703008581 },
	{ 0.652996233961648, 0.146246269259866 },
	{ 0.927005975926850, 0.157136361064887 }
      };



    pyramidrules.SetSize(order3d);
    for (k = 0; k < pyramidrules.Size(); k++)
      {
	const IntegrationRule & quadrule = *quadrules[k];
	double (*zrule) [2];
	int nz;
	switch (k+1)
	  {
	  case 1: 
	    zrule = qf_pyramidz_order1_points_weight; 
	    nz = 1;
	    break;
	  case 2: case 3:
	    zrule = qf_pyramidz_order3_points_weight; 
	    nz = 2;
	    break;
	    // case 4: case 5:
	  default:
	    zrule = qf_pyramidz_order5_points_weight; 
	    nz = 3;
	    break;
	  }

	IntegrationRule * pyramidrule = new IntegrationRule();
	
	double point[3], weight;
	
	for (i = 0; i < quadrule.GetNIP(); i++)
	  for (j = 0; j < nz; j++)
	    {
	      const IntegrationPoint & ipquad = quadrule.GetIP(i);
	      point[0] = zrule[j][0] * ipquad.Point()[0];
	      point[1] = zrule[j][0] * ipquad.Point()[1];
	      point[2] = 1-zrule[j][0];
	      weight = zrule[j][1] * ipquad.Weight();
	      
	      ip = new IntegrationPoint (point, weight);
	      ip->SetGlobNr (pyramidpoints.Append (ip)-1);
	      pyramidrule->AddIntegrationPoint (ip);
	    }
	pyramidrules[k] = pyramidrule;
      }
    */    


    for (p = 0; p <= 4; p++)
      GenerateIntegrationRule (ET_PYRAMID, p);
    /*
    pyramidrules.SetSize(order3d);
    for (k = 0; k < pyramidrules.Size(); k++)
      {
	const IntegrationRule & quadrule = *quadrules[k];
	const IntegrationRule & segrule = *segmentrules[k];
	IntegrationRule * pyramidrule = 
	  new IntegrationRule(quadrule.GetNIP()*segrule.GetNIP());
	
	double point[3], weight;
	
	for (i = 0; i < quadrule.GetNIP(); i++)
	  for (j = 0; j < segrule.GetNIP(); j++)
	    {
	      const IntegrationPoint & ipquad = quadrule.GetIP(i);
	      const IntegrationPoint & ipseg = segrule.GetIP(j);
	      point[0] = (1-ipseg(0)) * ipquad(0);
	      point[1] = (1-ipseg(0)) * ipquad(1);
	      point[2] = ipseg(0);
	      weight = ipseg.Weight() * sqr (1-ipseg(0)) * ipquad.Weight();
	      
	      ip = new IntegrationPoint (point, weight);
	      if (k <= 10)
		ip->SetGlobNr (pyramidpoints.Append (ip)-1);
	      pyramidrule->AddIntegrationPoint (ip);
	    }
	pyramidrules[k] = pyramidrule;
      }
    */




    static double qf_pyramid_lumping_points[][3] = 
      {
	{ 0, 0, 0 },
	{ 1, 0, 0 },
	{ 1, 1, 0 },
	{ 0, 1, 0 },
	{ 0, 0, 1-1e-14 },
      };
       
    static double qf_pyramid_lumping_weights[] = 
      {
	1.0/15.0,
	1.0/15.0,
	1.0/15.0,
	1.0/15.0,
	1.0/15.0,
      };    // not optimal, JS !!!
    
    for (i = 0; i < 5; i++)
      {
	ip = new IntegrationPoint (qf_pyramid_lumping_points[i],
				   qf_pyramid_lumping_weights[i]);
	ip->SetNr (i);
	ip->SetGlobNr (pyramidpoints.Append (ip)-1);
	pyramidlumping.AddIntegrationPoint (ip);
      }



    for (p = 0; p <= 4; p++)
      GenerateIntegrationRule (ET_HEX, p);


    /*
    hexrules.SetSize(order3d);
    for (k = 0; k < hexrules.Size(); k++)
      {
	const IntegrationRule & segmrule = *segmentrules[k];

	IntegrationRule * hexrule = 
	  new IntegrationRule(segmrule.GetNIP()*segmrule.GetNIP()*segmrule.GetNIP());
	
	double point[3], weight;
	int l;
	for (i = 0; i < segmrule.GetNIP(); i++)
	  for (j = 0; j < segmrule.GetNIP(); j++)
	    for (l = 0; l < segmrule.GetNIP(); l++)
	      {
		const IntegrationPoint & ipsegm1 = segmrule.GetIP(i);
		const IntegrationPoint & ipsegm2 = segmrule.GetIP(j);
		const IntegrationPoint & ipsegm3 = segmrule.GetIP(l);

		point[0] = ipsegm1.Point()[0];
		point[1] = ipsegm2.Point()[0];
		point[2] = ipsegm3.Point()[0];
		weight = ipsegm1.Weight() * ipsegm2.Weight() * ipsegm3.Weight();
	      
		ip = new IntegrationPoint (point, weight);
		int found = 0;
		for (int m = 0; m < hexpoints.Size(); m++)
		  if (*ip == *hexpoints[m])
		    {
		      found = 1;
		      delete ip;
		      ip = hexpoints[m];
		    }
		if (!found)
		  ip->SetGlobNr (hexpoints.Append (ip)-1);
	      

		hexrule->AddIntegrationPoint (ip);
	      }
	hexrules[k] = hexrule;
      }
    */


    /*
    cout << "Check trig intrule:";
    for (int order = 0; order < 10; order++)
      {
	cout << "order = " << order << endl;
	const IntegrationRule & rule = *trigrules[order];
	for (int ox = 0; ox <= 4; ox++)
	  for (int oy = 0; ox+oy <= 4; oy++)
	    {
	      double sum = 0;
	      for (int j = 0; j < rule.GetNIP(); j++)
		{
		  const IntegrationPoint & ip = rule[j];
		  sum += ip.Weight() * pow (ip(0), ox) * pow (ip(1), oy);
		}
	      cout << "\\int x^" << ox << " y^" << oy << " = " << sum << endl;
	    }
      }
    */      

  }




  IntegrationRules :: ~IntegrationRules ()
  {
    for (int i = 0; i < segmentrules.Size(); i++)
      delete segmentrules[i];

    for (int i = 0; i < trigrules.Size(); i++)
      delete trigrules[i];

    for (int i = 0; i < trignodalrules.Size(); i++)
      delete trignodalrules[i];

    for (int i = 0; i < quadrules.Size(); i++)
      delete quadrules[i];

    for (int i = 0; i < tetrules.Size(); i++)
      delete tetrules[i];

    for (int i = 0; i < tetnodalrules.Size(); i++)
      delete tetnodalrules[i];

    for (int i = 0; i < prismrules.Size(); i++)
      delete prismrules[i];

    for (int i = 0; i < pyramidrules.Size(); i++)
      delete pyramidrules[i];

    for (int i = 0; i < hexrules.Size(); i++)
      delete hexrules[i];
  }


  const ARRAY<IntegrationPoint*> & IntegrationRules ::
  GetIntegrationPoints (ELEMENT_TYPE eltyp) const
  {
    switch (eltyp)
      {
      case ET_SEGM:
	return segmentpoints;
      case ET_TRIG:
	return trigpoints;
      case ET_QUAD:
	return quadpoints;
      case ET_TET:
	return tetpoints;
      case ET_PRISM:
	return prismpoints;
      case ET_PYRAMID:
	return pyramidpoints;
      case ET_HEX:
	return hexpoints;
      }

    stringstream str;
    str<< "no intpoints available for element type " 
       << ElementTopology::GetElementName(eltyp) << endl;
    throw Exception (str.str());
  }



  const IntegrationRule & IntegrationRules :: 
  SelectIntegrationRule (ELEMENT_TYPE eltyp, int order) const
  {
    const ARRAY<IntegrationRule*> * ira;
  
    switch (eltyp)
      {
      case ET_SEGM:
	ira = &segmentrules; break;
      case ET_TRIG:
	ira = &trigrules; break;
      case ET_QUAD:
	ira = &quadrules; break;
      case ET_TET:
	ira = &tetrules; break;
      case ET_PYRAMID:
	ira = &pyramidrules; break;
      case ET_PRISM:
	ira = &prismrules; break;
      case ET_HEX:
	ira = &hexrules; break;
      default:
	{
	  stringstream str;
	  str << "no integration rules for element " << int(eltyp) << endl;
	  throw Exception (str.str());
	}
      }

    if (order < 0) 
      { order = 0; }

    if (order >= ira->Size() || (*ira)[order] == 0)
      {
	return const_cast<IntegrationRules&> (*this).
	  GenerateIntegrationRule (eltyp, order);
	/*
	stringstream str;
	str << "Integration rule of order " << order 
	    << " not available for element type " 
	    << ElementTopology::GetElementName(eltyp) << endl;
	throw Exception (str.str());
	*/
      }

    return *((*ira)[order]);
  }
 



  const IntegrationRule & IntegrationRules :: 
  GenerateIntegrationRule (ELEMENT_TYPE eltyp, int order)
  {
    ARRAY<IntegrationRule*> * ira;
  
    switch (eltyp)
      {
      case ET_SEGM:
	ira = &segmentrules; break;
      case ET_TRIG:
	ira = &trigrules; break;
      case ET_QUAD:
	ira = &quadrules; break;
      case ET_TET:
	ira = &tetrules; break;
      case ET_PYRAMID:
	ira = &pyramidrules; break;
      case ET_PRISM:
	ira = &prismrules; break;
      case ET_HEX:
	ira = &hexrules; break;
      default:
	{
	  stringstream str;
	  str << "no integration rules for element " << int(eltyp) << endl;
	  throw Exception (str.str());
	}
      }

    if (ira -> Size() < order+1)
      {
	int oldsize = ira -> Size();
	ira -> SetSize (order+1);
	for (int i = oldsize; i < order+1; i++)
	  (*ira)[i] = 0;
      }

    if ( (*ira)[order] == 0)
      {
	switch (eltyp)
	  {
	  case ET_SEGM:
	    {
	      ARRAY<double> xi, wi;
	      ComputeGaussRule (order/2+1, xi, wi);
	      IntegrationRule * rule = new IntegrationRule (xi.Size());
	      double xip[3] = { 0, 0, 0 };
	      for (int j = 0; j < xi.Size(); j++)
		{
		  xip[0] = xi[j];
		  IntegrationPoint * ip = new IntegrationPoint (xip, wi[j]);
		  ip->SetNr (j);
		  if (order <= 20)
		    ip->SetGlobNr (segmentpoints.Append (ip)-1);
		  rule->AddIntegrationPoint (ip);
		}
	      segmentrules[order] = rule;	      
	      break;
	    }

	  case ET_TRIG:
	    {
	      const IntegrationRule & segmrule = SelectIntegrationRule (ET_SEGM, order+1);
	      IntegrationRule * trigrule = new IntegrationRule(segmrule.GetNIP()*segmrule.GetNIP());
	
	      double point[3], weight;
	      
	      int ii = 0;
	      for (int i = 0; i < segmrule.GetNIP(); i++)
		for (int j = 0; j < segmrule.GetNIP(); j++)
		  {
		    const IntegrationPoint & ipsegmi = segmrule.GetIP(i);
		    const IntegrationPoint & ipsegmj = segmrule.GetIP(j);
		    
		    point[0] = ipsegmi.Point()[0];
		    point[1] = ipsegmj.Point()[0]*(1-point[0]);
		    point[2] = 0;
		    
		    weight = ipsegmi.Weight() *
		      ipsegmj.Weight() * (1-point[0]);
		    
		    IntegrationPoint * ip = new IntegrationPoint (point, weight);
		    ip->SetNr (ii); ii++;
		    if (order <= 10)
		      ip->SetGlobNr (trigpoints.Append (ip)-1);
		    trigrule->AddIntegrationPoint (ip);
		  }
	      trigrules[order] = trigrule;
	      break;
	    }


	  case ET_QUAD:
	    {
	      const IntegrationRule & segmrule = SelectIntegrationRule (ET_SEGM, order);
	      IntegrationRule * quadrule = new IntegrationRule(segmrule.GetNIP()*segmrule.GetNIP());
	      
	      double point[3], weight;

	      int ii = 0;
	      for (int i = 0; i < segmrule.GetNIP(); i++)
		for (int j = 0; j < segmrule.GetNIP(); j++)
		  {
		    const IntegrationPoint & ipsegm1 = segmrule.GetIP(i);
		    const IntegrationPoint & ipsegm2 = segmrule.GetIP(j);
		    point[0] = ipsegm1.Point()[0];
		    point[1] = ipsegm2.Point()[0];
		    point[2] = 0;
		    weight = ipsegm1.Weight() * ipsegm2.Weight();
		    
		    IntegrationPoint * ip = new IntegrationPoint (point, weight);
		    ip->SetNr (ii); ii++;
		    if (order <= 10)
		      ip->SetGlobNr (quadpoints.Append (ip)-1);
		    quadrule->AddIntegrationPoint (ip);
		  }
	      quadrules[order] = quadrule;
	      break;
	    }
  

	  case ET_TET:
	    {
	      // tet-rules by degenerated tensor product rules:

	      const IntegrationRule & segmrule = SelectIntegrationRule (ET_SEGM, order+2);
	      IntegrationRule * tetrule = 
		new IntegrationRule(segmrule.GetNIP()* segmrule.GetNIP()* segmrule.GetNIP());
	
	      double point[3], weight;
	      
	      int ii = 0;
	      for (int i = 0; i < segmrule.GetNIP(); i++)
		for (int j = 0; j < segmrule.GetNIP(); j++)
		  for (int k = 0; k < segmrule.GetNIP(); k++)
		    {
		      const IntegrationPoint & ipsegmi = segmrule.GetIP(i);
		      const IntegrationPoint & ipsegmj = segmrule.GetIP(j);
		      const IntegrationPoint & ipsegmk = segmrule.GetIP(k);

		      point[0] = ipsegmi.Point()[0];
		      point[1] = ipsegmj.Point()[0]*(1-point[0]);
		      point[2] = ipsegmk.Point()[0]*(1-point[0]-point[1]);
		      
		      weight = ipsegmi.Weight() *
			ipsegmj.Weight() * (1-point[0]) *
			ipsegmk.Weight() * (1-point[0]-point[1]);
	      
		      IntegrationPoint * ip = new IntegrationPoint (point, weight);
		      ip->SetNr (ii); ii++;
		      if (order <= 6)
			ip->SetGlobNr (tetpoints.Append (ip)-1);
		      tetrule->AddIntegrationPoint (ip);
		    }
	      tetrules[order] = tetrule;
	      break;
	    }

	  case ET_HEX:
	    {
	      const IntegrationRule & segmrule = SelectIntegrationRule (ET_SEGM, order);

	      IntegrationRule * hexrule = 
		new IntegrationRule(segmrule.GetNIP()*segmrule.GetNIP()*segmrule.GetNIP());
	
	      double point[3], weight;
	      int ii = 0;
	      for (int i = 0; i < segmrule.GetNIP(); i++)
		for (int j = 0; j < segmrule.GetNIP(); j++)
		  for (int l = 0; l < segmrule.GetNIP(); l++)
		    {
		      const IntegrationPoint & ipsegm1 = segmrule.GetIP(i);
		      const IntegrationPoint & ipsegm2 = segmrule.GetIP(j);
		      const IntegrationPoint & ipsegm3 = segmrule.GetIP(l);
		      
		      point[0] = ipsegm1.Point()[0];
		      point[1] = ipsegm2.Point()[0];
		      point[2] = ipsegm3.Point()[0];
		      weight = ipsegm1.Weight() * ipsegm2.Weight() * ipsegm3.Weight();
		      
		      IntegrationPoint * ip = 
			new IntegrationPoint (point, weight);
		      
		      ip->SetNr (ii); ii++;
		      if (order <= 6)
			ip->SetGlobNr (hexpoints.Append (ip)-1);
		      hexrule->AddIntegrationPoint (ip);
		    }
	      hexrules[order] = hexrule;
	      break;
	    }


	  case ET_PRISM:
	    {
	      const IntegrationRule & segmrule = SelectIntegrationRule (ET_SEGM, order);
	      const IntegrationRule & trigrule = SelectIntegrationRule (ET_TRIG, order);

	      IntegrationRule * prismrule = 
		new IntegrationRule(segmrule.GetNIP()*trigrule.GetNIP());
      
	      double point[3], weight;
	      
	      int ii = 0;
	      for (int i = 0; i < segmrule.GetNIP(); i++)
		for (int j = 0; j < trigrule.GetNIP(); j++)
		  {
		    const IntegrationPoint & ipsegm = segmrule.GetIP(i);
		    const IntegrationPoint & iptrig = trigrule.GetIP(j);
		    point[0] = iptrig.Point()[0];
		    point[1] = iptrig.Point()[1];
		    point[2] = ipsegm.Point()[0];
		    weight = iptrig.Weight() * ipsegm.Weight();
	      
		    IntegrationPoint * ip = 
		      new IntegrationPoint (point, weight);
		    ip->SetNr (ii); ii++;
		    if (order <= 6)
		      ip->SetGlobNr (prismpoints.Append (ip)-1);
		    prismrule->AddIntegrationPoint (ip);
		  }
	      prismrules[order] = prismrule;
	      break;
	    }


	  case ET_PYRAMID:
	    {
	      const IntegrationRule & quadrule = SelectIntegrationRule (ET_QUAD, order);
	      const IntegrationRule & segrule = SelectIntegrationRule (ET_SEGM, order);
	      IntegrationRule * pyramidrule = 
		new IntegrationRule(quadrule.GetNIP()*segrule.GetNIP());
	
	      double point[3], weight;
	      
	      int ii = 0;
	      for (int i = 0; i < quadrule.GetNIP(); i++)
		for (int j = 0; j < segrule.GetNIP(); j++)
		  {
		    const IntegrationPoint & ipquad = quadrule.GetIP(i);
		    const IntegrationPoint & ipseg = segrule.GetIP(j);
		    point[0] = (1-ipseg(0)) * ipquad(0);
		    point[1] = (1-ipseg(0)) * ipquad(1);
		    point[2] = ipseg(0);
		    weight = ipseg.Weight() * sqr (1-ipseg(0)) * ipquad.Weight();
		    
		    IntegrationPoint * ip = new IntegrationPoint (point, weight);
		    ip->SetNr (ii); ii++;
		    if (order <= 6)
		      ip->SetGlobNr (pyramidpoints.Append (ip)-1);
		    pyramidrule->AddIntegrationPoint (ip);
		  }
	      pyramidrules[order] = pyramidrule;
	      break;
	    }
	  }
      }

    if ( (*ira)[order] == 0)
      {
	stringstream str;
	str << "could not generate Integration rule of order " << order 
	    << " for element type " 
	    << ElementTopology::GetElementName(eltyp) << endl;
	throw Exception (str.str());
      }
    /*
    else
      {
	cout << "could generate ir of order " << order << " for element " 
	     << ElementTopology::GetElementName(eltyp) << endl;
      }
    */
    return *(*ira)[order];
  }










  const IntegrationRule & IntegrationRules :: 
  SelectLumpingIntegrationRule (ELEMENT_TYPE eltyp) const
  {
    const IntegrationRule * ir = &triglumping;
    //  cerr << "get lumping rule" << endl;

    switch (eltyp)
      {
      case ET_SEGM:
	ir = &segmentlumping; break;
      case ET_TRIG:
	ir = &triglumping; break;
      case ET_QUAD:
	ir = &quadlumping; break;
      case ET_TET:
	ir = &tetlumping; break;
      case ET_PRISM:
	ir = &prismlumping; break;
      case ET_PYRAMID:      
	ir = &pyramidlumping; break;
      default:
	{
	  cout << "no lumping integration rules for element " << int(eltyp) << endl;
	  ir = &SelectIntegrationRule (eltyp, 1);
	  //	ir = &triglumping;
	}
      }

    return *ir;
  }
  



  const IntegrationRule & IntegrationRules :: 
  SelectNodalIntegrationRule (ELEMENT_TYPE eltyp, int order) const
  {
    const IntegrationRule * ir;
    ir = NULL;

    switch (eltyp)
      {
      case ET_TET:
	ir = tetnodalrules[order-1]; break;
      case ET_TRIG:
	if (order == 1)
	  ir = &triglumping;
	else if (order == 2)
	  ir = &triglumping2;
	else ir = trignodalrules[order-1]; break;
	break;
      case ET_PYRAMID:
	ir = &pyramidlumping; break;
      case ET_PRISM:
	ir = &prismlumping; break;
      case ET_QUAD:
	ir = &quadlumping; break;
      default:
	{
	  cout << "no nodal integration rules for element " << int(eltyp) << endl;
	  ir = &triglumping;
	}
      }
    if (!ir)
      {
	cout << "no nodal integration rules for element " << int(eltyp)
	     << ", order " << order << endl;
	ir = &triglumping;
      }
    return *ir;
  }
  


  int Integrator :: common_integration_order = -1;

  const IntegrationRules & GetIntegrationRules ()
  {
    static IntegrationRules intrules;
    return intrules;
  }
}