class TetShapesInnerLegendre
{
public:
  template <typename Sx, typename Sy, typename Sz, typename T>
  static int Calc (int n, Sx x, Sy y, Sz z, T & values)
  {
    int ii = 0, i, j, k;
    ArrayMem<Sx, 20> polx(n+1), poly(n+1), polz(n+1);
    
    ScaledLegendrePolynomial (n-4, x, 1-y-z, polx);
    ScaledLegendrePolynomial (n-4, 2*y - (1-z) , (1-z),poly);  //SZ (1-z)-1, poly);
    LegendrePolynomial (n-4, 2*z-1, polz);
    Sx bub = (1-x-y-z) * (1+x-y-z) * y * z ; 
  
    for (i = 0; i <= n-4; i++)
      for (j = 0; j <= n-4-i; j++)
	for (k = 0; k <= n-4-i-j; k++)
	  values[ii++] = bub * polx[i] * poly[j] * polz[k];
    
    return ii;
  }
  template <typename Sx, typename Sy, typename Sz, typename T>
  static void CalcSplitted (int n, Sx x, Sy y, Sz z, T & val1, T & val2, T & val3)
  {
    ArrayMem<Sx, 20> polx(n+1), poly(n+1), polz(n+1);
    
    ScaledLegendrePolynomial (n-4, x, 1-y-z, polx);
    ScaledLegendrePolynomial (n-4, 2*y - (1-z) , (1-z),poly);  //SZ (1-z)-1, poly);
    LegendrePolynomial (n-4, 2*z-1, polz);
    Sx bub1 = (1-x-y-z) * (1+x-y-z);
    Sx bub2 = y; 
    Sx bub3 = z; 
  
    for (int i = 0; i <= n-4; i++)
      {
	val1[i] = bub1 * polx[i];
	val2[i] = bub2 * poly[i]; 
	val3[i] = bub3 * polz[i]; 
	/*
	val1[i] = Mult (bub1, polx[i]);
	val2[i] = Mult (bub2, poly[i]); 
	val3[i] = Mult (bub3, polz[i]); 
	*/
      } 
  }
  
};



class TetShapesInnerJacobi
{
public:
  template <typename Sx, typename Sy, typename Sz, typename T>
  static int Calc (int n, Sx x, Sy y, Sz z, T & values)
  {
    int ii = 0, j, k;
    ArrayMem<Sx, 20> polx(n+1), poly(n+1), polz(n+1);

    Sx bub = y * z * (1-x-y-z) * (1+x-y-z);
    ScaledJacobiPolynomial (n-4, x, (1-y-z), 2, 2, polx);
    
    for (int ix = 0; ix <= n-4; ix++)
      {
	ScaledJacobiPolynomial (n-4, (2*y-1+z),(1-z), 2*ix+5, 2, poly);
	JacobiPolynomial (n-4, 2*z-1, 2*ix+5, 2, polz);
	
	for (j = 0; j <= n-4-ix; j++)
	  for (k = 0; k <= n-4-ix-j; k++)
	    values[ii++] = bub * polx[ix] * poly[j] * polz[k];
      }
    return ii;
  }
};

class TetShapesFaceLegendre
{
public:
  template <typename Sx, typename Sy, typename Sz, typename T>
  static int Calc (int n, Sx x, Sy y, Sz z, T & values)
  {
    int ii = 0, i, j, k;
    ArrayMem<Sx, 20> polx(n+1), poly(n+1);
    
    ScaledLegendrePolynomial (n-3, x, 1-y-z, polx);
    ScaledLegendrePolynomial (n-3, 2*y-(1-z),(1-z), poly);
    Sx bub = (1-x-y-z) * (1+x-y-z)*y; 

    for (i = 0; i <= n-3; i++)
      for (j = 0; j <= n-3-i; j++)
	values[ii++] = bub * polx[i] * poly[j];
 
    return ii;
  }
  template <typename Sx, typename Sy, typename Sz, typename T>
  static void CalcSplitted (int n, Sx x, Sy y, Sz z, T & val1, T & val2)
  {
    /*
    int i;
    ArrayMem<Sx, 20> polx(n+1), poly(n+1);
    
    ScaledLegendrePolynomial (n-3, x, 1-y-z, polx);
    ScaledLegendrePolynomial (n-3, 2*y-(1-z),(1-z), poly);

    Sx bub1 = (1-x-y-z) * (1+x-y-z);
    //    Sx bub1 = Mult ( Minus ( Minus (1-x, y), z),
    //		     Minus ( Minus (1+x, y), z));
    Sx bub2 = y; 

    for (i = 0; i <= n-3; i++)
      {
	val1[i] = bub1 * polx[i];
	val2[i] = bub2 * poly[i]; 
	//	val1[i] = Mult (bub1, polx[i]);
	//	val2[i] = Mult (bub2, poly[i]); 
      } 
    return i;
    */

    ScaledLegendrePolynomial (n-3, x, 1-y-z, val1);
    ScaledLegendrePolynomial (n-3, 2*y-(1-z),(1-z), val2);

    Sx bub1 = (1-x-y-z) * (1+x-y-z);
    Sx bub2 = y; 

    for (int i = 0; i <= n-3; i++)
      {
	val1[i] *= bub1;
	val2[i] *= bub2;
      } 
  }
};




class TetShapesFaceJacobi
{
public:
  template <typename Sx, typename Sy, typename Sz, typename T>
  static int Calc (int n, Sx x, Sy y, Sz z, T & values)
  {
    int ii = 0;
    ArrayMem<Sx, 20> polx(n+1), poly(n+1);

    Sx bub = y * (1-x-y-z) * (1+x-y-z);
    ScaledJacobiPolynomial (n-3, x, 1-y-z, 2, 2, polx);
    
    for (int ix = 0; ix <= n-3; ix++)
      {
	ScaledJacobiPolynomial (n-3, 2*y-1+z, 1-z, 2*ix+5, 2, poly);
	for (int j = 0; j <= n-3-ix; j++)
	  values[ii++] = bub * polx[ix] * poly[j];
      }
    return ii;
  }
};


class TetShapesFaceOpt1
{
public:
  template <typename Sx, typename Sy, typename Sz, typename T>
  static int Calc (int n, Sx x, Sy y, Sz z, T & values)
  {
    int nd = Calc1 (n, x, y, z, values);
    ArrayMem<Sx, 100> hvalues(nd);

    Sx lam1 = 0.5 * (1 + x - y - z);
    Sx lam2 = 0.5 * (1 - x - y - z);
    Sx lam3 = y;
    Sx lam4 = z;

    Sx hlam1 = 0;
    Sx hlam2 = lam2;
    Sx hlam3 = lam3;
    Sx hlam4 = lam4 + lam1;

    Sx hx = hlam1 - hlam2;
    Sx hy = hlam3;
    Sx hz = hlam4;

    Sx frac;
    if (hlam4 < 1e-12)
      frac = 0.0;
    else
      frac = lam4 / hlam4;

    Calc1 (n, hx, hy, hz, hvalues);
    for (int i = 0; i < nd; i++)
      values[i] -= frac * hvalues[i];

    return nd;
  }
  
  template <typename Sx, typename Sy, typename Sz, typename T>
  static int Calc1 (int n, Sx x, Sy y, Sz z, T & values)
  {
    int nd = Calc2 (n, x, y, z, values);
    ArrayMem<Sx, 100> hvalues(nd);

    Sx lam1 = 0.5 * (1 + x - y - z);
    Sx lam2 = 0.5 * (1 - x - y - z);
    Sx lam3 = y;
    Sx lam4 = z;

    Sx hlam1 = lam1;
    Sx hlam2 = 0;
    Sx hlam3 = lam3;
    Sx hlam4 = lam4 + lam2;

    Sx hx = hlam1 - hlam2;
    Sx hy = hlam3;
    Sx hz = hlam4;

    Sx frac;
    if (hlam4 < 1e-12)
      frac = 0.0;
    else
      frac = lam4 / hlam4;

    Calc2 (n, hx, hy, hz, hvalues);
    for (int i = 0; i < nd; i++)
      values[i] -= frac * hvalues[i];
    
    return nd;
  }

  template <typename Sx, typename Sy, typename Sz, typename T>
  static int Calc2 (int n, Sx x, Sy y, Sz z, T & values)
  {
    Sx * hp = &values[0];
    int nd = Calc3 (n, x, y, z, hp);
    ArrayMem<Sx, 100> hvalues(nd);

    Sx lam1 = 0.5 * (1 + x - y - z);
    Sx lam2 = 0.5 * (1 - x - y - z);
    Sx lam3 = y;
    Sx lam4 = z;

    Sx hlam1 = lam1;
    Sx hlam2 = lam2;
    Sx hlam3 = 0;
    Sx hlam4 = lam4 + lam3;

    Sx hx = hlam1 - hlam2;
    Sx hy = hlam3;
    Sx hz = hlam4;

    Sx frac;
    if (hlam4 < 1e-12)
      frac = 0.0;
    else
      frac = lam4 / hlam4;

    hp = &hvalues[0];
    Calc3 (n, hx, hy, hz, hp);
    for (int i = 0; i < nd; i++)
      values[i] -= frac * hvalues[i];
    
    return nd;
  }






  template <typename Sx, typename Sy, typename Sz, typename T>
  static int Calc3 (int n, Sx x, Sy y, Sz z, T & values)
  {
    int ii = 0, i, j, k;
    ArrayMem<Sx, 20> polx(n+1), poly(n+1);

    const IntegrationRule & rule =
      GetIntegrationRules().SelectIntegrationRule (ET_TRIG, n+2);

    for (int ix = 0; ix <= n-3; ix++)
      for (j = 0; j <= n-3-ix; j++)
	values[ii++] = 0;
    for (i = 0; i < rule.GetNIP(); i++)
      {
	ii = 0;
	const IntegrationPoint & ip = rule.GetIP(i);

	Sx hx = x + z * (-1+2*ip(0)+ip(1));
	Sy hy = y + z * ip(1);

	Sx bub = hy * (1-hx-hy) * (1+hx-hy);
	ScaledJacobiPolynomial (n-3, hx, 1-hy, 2, 2, polx);

	Sx fac = 2 * bub * ip.Weight(); //  / (z*z);

	for (int ix = 0; ix <= n-3; ix++)
	  {
	    ScaledJacobiPolynomial (n-3, 2*hy-1, 1, 2*ix+5, 2, poly);
	    for (j = 0; j <= n-3-ix; j++)
	      values[ii++] += fac * polx[ix] * poly[j];
	  }
      }
    return ii;
  }
};




class TetShapesFaceOpt2
{
public:
  template <typename Sx, typename Sy, typename Sz, typename T>
  static int Calc (int n, Sx x, Sy y, Sz z, T & values)
  {
    int ii = 0, i, j, k;
    ArrayMem<Sx, 20> polx(n+1), poly(n+1);

    const IntegrationRule & rule =
      GetIntegrationRules().SelectIntegrationRule (ET_TRIG, n+2);

    for (int ix = 0; ix <= n-3; ix++)
      for (j = 0; j <= n-3-ix; j++)
	values[ii++] = 0;
    for (i = 0; i < rule.GetNIP(); i++)
      {
	ii = 0;
	const IntegrationPoint & ip = rule.GetIP(i);

	Sx hx = x + z * (-1+2*ip(0)+ip(1));
	Sy hy = y + z * ip(1);

	//Sx bub = hy * (1-hx-hy) * (1+hx-hy);
	ScaledJacobiPolynomial (n-3, hx, 1-hy, 2, 2, polx);

	Sx fac = 2  * ip.Weight(); //  / (z*z);

	for (int ix = 0; ix <= n-3; ix++)
	  {
	    ScaledJacobiPolynomial (n-3, 2*hy-1, 1, 2*ix+5, 2, poly);
	    for (j = 0; j <= n-3-ix; j++)
	      values[ii++] += fac * polx[ix] * poly[j];
	  }
      }
      int jj = ii;
      ArrayMem<Sx, 100> hvalues(jj);
      for (int i = 0; i < ii; i++)
          hvalues[i] = values[i];
      ii = 0; jj = 0;
      for (int k = 0; k <= n-3; k++)
      {
         for (int m = 0; m <= n-3-k; m++)
         {
          values[ii++] = y*(1-x-y-z)*(1+x-y-z)*hvalues[jj++];
	 }
      }

    return ii;
  }
};










/*
  /// compute shape
  virtual void CalcShape (const IntegrationPoint & ip,
			  FlatVector<> shape) const
  {
    double lam1 = ip(0);
    double lam2 = 1-ip(0)-ip(1);
    double lam3 = ip(1);
    double bub = lam1*lam2*lam3;

    ArrayMem<double, 50> polx(order-2);
    ArrayMem<double, 50> poly(order-2);

    double fac = 1;
    for (int i = 0; i <= order-3; i++)
      {
	polx[i] *= fac;
	fac *= (1-lam3);
      }

    int ii = 0;
    for (int ix = 0; ix <= order-3; ix++)
      {
	JacobiPolynomial (order-3, lam3-(lam1+lam2), 2*ix+5, 2, poly);
	for (int iy = 0; iy <= order-3-ix; iy++)
	  shape(ii++) = bub*polx[ix] * poly[iy];
      }
  }

*/