/*
$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(°ree);
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];
}
}
//-------------------------------------------------------------------------
syntax highlighted by Code2HTML, v. 0.9.1