/*
 $Id: integ.cc,v 1.2 1996/10/08 08:17:39 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "integ.h"
#include "utils.h"
#include "mzibutil.h"


// static const double OneSixth = 1./6.;		// volume of a tetrahedron
// static const double OneThird = 1./3.;
// static const double TwoThird = 2./3.;
//-------------------------------------------------------------------------

IntegForm:: IntegForm() : ID(0) { }

//-------------------------------------------------------------------------

void IntegForm:: Info()
{
    cout << "\n\t\t  formula          "  << ID->formName;
    cout << "\n\t\t  spaceDim =       "  << ID->dim;
    cout << "\n\t\t  no of i-Points = "  << ID->noOfIPoints;
    cout << "\n\t\t  degree =         "  << ID->degree << "\n";
}
//-------------------------------------------------------------------------


void IntegForm:: missingFormula(const int dim, const int degree) const
{
    cout << "\n*** GaussForm/LumpedForm: no formula available with degree = " 
         << degree << "in space dimension = " << dim << "\n";
    cout.flush(); abort();
}
//-------------------------------------------------------------------------

void IntegForm:: missingItem(const char* item, const char* formName0)
{
    cout << "\n*** Formula " << formName0 << " without " << item << "\n";
    cout.flush(); abort(); 
}
//-------------------------------------------------------------------------


void IntegForm:: Coordinates(int ip, Vector<Real>& u)
{
    for (int i=1; i<=ID->dim; ++i)  u[i] = ID->point(ip,i);
} 
//-------------------------------------------------------------------------

void IntegForm:: readFormula(const char* formName0, Bool completeMatch) 
{
   int spaceDim, degree, noOfIPoints; 

   Parser parser("integ.dat", "end", CommentFlag);

   if (parser.findWord(formName0, completeMatch) != parser.ok) 
   {
	cout << "\n*** GaussForm/LumpedForm: no formula available with name = " 
	     << formName0  << "\n";
	cout.flush(); abort();
   }

   if (parser.findWord("DIMENSION") != parser.ok) missingItem("DIMENSION",
							      formName0);
   if (parser.findWord("=")  != parser.ok) missingItem("=",formName0);
   parser.getValue(&spaceDim);

   if (parser.findWord("NO_OF_POINTS") != parser.ok) missingItem("NO_OF_POINTS",
								 formName0);
   if (parser.findWord("=")  != parser.ok) missingItem("=",formName0);
   parser.getValue(&noOfIPoints);

   if (parser.findWord("DEGREE") != parser.ok) missingItem("DEGREE",formName0);
   if (parser.findWord("=")  != parser.ok) missingItem("=",formName0);
   parser.getValue(&degree);


   if (parser.findWord("POINTS")  != parser.ok) missingItem("POINTS",formName0);
   if (parser.findWord("WEIGHTS") != parser.ok) missingItem("WEIGHTS",formName0);

   ID = new IntegData(degree,noOfIPoints,spaceDim);

   for (int i=1; i<=noOfIPoints; i++)
   {
       for (int k=1; k<=spaceDim; ++k) parser.getValue(&(ID->point(i,k)));
       parser.getValue(&(ID->weight[i]));
   }

   ID->formName = new char[strlen(formName0)+1];
   strcpy(ID->formName, formName0);
}
//-------------------------------------------------------------------------


GaussForm:: GaussForm(int spaceDim, int degree, const char* formName0)
{
    char name[128];
    
    if (formName0) readFormula(formName0);      // formula specified by user
    else
    {
	if (spaceDim == 1)
	{
	    int noOfIPoints0 = (int)((degree+1+1)/2.0);
	    ID = new GLIntegData1d(degree,noOfIPoints0);
	}	    
	else if (spaceDim == 2)
	{ 
	    switch(degree)
	    {
	            case 1: readFormula("INT_2D_T_1_1");  break;
		    case 2: readFormula("INT_2D_T_2_3A"); break;
		    case 3: readFormula("INT_2D_T_3_4");  break;
		    case 5: readFormula("INT_2D_T_5_7");  break;
		    case 7: readFormula("INT_2D_T_7_13"); break;
		    case 9: readFormula("INT_2D_T_9_19"); break;
		    default:  
		      {  sprintf(name,"INT_2D_T_%1d", degree); 
			 readFormula(name,False);
			 break;
		     }
		  } 
	}
	else if (spaceDim == 3)
	{
	    switch(degree)
	    {
	      	    case 1: readFormula("INT_3D_T_1_1");  break;
		    case 2: readFormula("INT_3D_T_2_4");  break;
		    case 3: readFormula("INT_3D_T_3_5");  break;
		    case 5: readFormula("INT_3D_T_5_15"); break;
		    case 7: readFormula("INT_3D_T_7_35"); break;
		    case 9: readFormula("INT_3D_T_9_70"); break;
		    default:  
		      {  sprintf(name,"INT_3D_T_%1d", degree); 
			 readFormula(name,False);
			 break;
		     }
		  } 
	}
	else missingFormula(spaceDim,degree);
    }   
}
//-------------------------------------------------------------------------


LumpedForm::LumpedForm(int spaceDim, int degree, const char* formName0)
{
    char name[128];

    if (formName0) readFormula(formName0); 	// formula specified by user
    else			
    {
	switch (spaceDim)
	{
		case 1: 
			switch(degree)
			{
			  case 1: readFormula("LumpedLineOrder1");  break;
			  case 2: readFormula("LumpedLineOrder2");  break;
                          default: 
				  {  sprintf(name,"LumpedLineOrder%1d",degree); 
				     readFormula(name,False);
				     break;
				  }
			} 
			break;

		case 2: 
			switch(degree)
			{
			  case 1: readFormula("LumpedTriOrder1");  break;
			  case 2: readFormula("LumpedTriOrder2");  break;
                          default: 
				  {  sprintf(name,"LumpedLineOrder%1d",degree); 
				     readFormula(name,False);
				     break;
				  }
			}
                        break;

		case 3: 
			switch(degree)
			{
			  case 1: readFormula("LumpedTetraOrder1");  break;
			  case 2: readFormula("LumpedTetraOrder2");  break;
                          default: 
				  {  sprintf(name,"LumpedLineOrder%1d",degree); 
				     readFormula(name,False);
				     break;
				  }
			} 
                        break;

		default: missingFormula(spaceDim,degree); break;
	}
   }   
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

IntegData:: IntegData(int degree0, int noOfIPoints0, int dim0) : 

		formName(0), noOfIPoints(noOfIPoints0), 
		degree(degree0), dim(dim0),
                point(noOfIPoints0,dim0), weight(noOfIPoints0)
{ }
//-------------------------------------------------------------------------


GLIntegData1d:: GLIntegData1d(int degree0, int noOfIPoints0) 

	: IntegData(degree0, noOfIPoints0, 1) 
{ 
    int i;
    noOfIPoints = (int)((degree+1+1)/2.0);
    
    Vector<Real> x(noOfIPoints);

    GaussLegendreData(x, weight,noOfIPoints);

    FORALL(weight,i) point(i,1) = x[i];
    degree   = 2*noOfIPoints-1;
    formName = new char[strlen("Gauss-Legendre")+1];
    strcpy(formName, "Gauss-Legendre");
}
//-------------------------------------------------------------------------


void GLIntegData1d:: GaussLegendreData(Vector<Real>& x, Vector<Real>& w, int n)
{
    //   gauss-legendre integration:
    //   calculate points and weights in geometrical range x1...x2 for order n 

    //   Order is 2*n - 1
    //   taken from Numerical Recipes ( = gauleg(...))

    double x1 = 0.0; double x2 = 1.0;
    int m,j,i;
    double z1,z,xm,xl,pp,p3,p2,p1;
    double EPS = 3.0e-11;

    m=(n+1)/2;
    xm=0.5*(x2+x1);
    xl=0.5*(x2-x1);
    for (i=1;i<=m;i++)  
    {
	z=cos(3.141592653589793*(i-0.25)/(n+0.5));
	do {
	    p1=1.0;
	    p2=0.0;
	    for (j=1;j<=n;j++) {
		p3=p2;
		p2=p1;
		p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
	    }
	    pp=n*(z*p1-p2)/(z*z-1.0);
	    z1=z;
	    z=z1-p1/pp;
	} while (fabs(z-z1) > EPS);
	x[i]=xm-xl*z;
	x[n+1-i]=xm+xl*z;
	w[i]=2.0*xl/((1.0-z*z)*pp*pp);
	w[n+1-i]=w[i];
    }
}
//-------------------------------------------------------------------------