sandbox/prouvost/utils/gauss_quadrature.h
Integration, Gauss quadrature and total error computation
At the origin, we where interested in computing accurately the total error by comparing a numerical solution with a known analytical solution. Thus, we implemented some Gauss quadrature rules to compute the integral in \mathcal{L}^p norm of the difference between a scalar field and a known analytical function. Incidentally, the integral of a scalar field (in any \mathcal{L}^p norm) can also be computed.
The Gauss integrations are written in 2D, and some of them (but not all of them) in 3D.
Source for the coefficients of the Gauss quadratures : Krylov V. I., Approximate calculation of integrals, Dover publications, 2005
Usefull functions
We had some odd results with basilisk interpolation function in particular cases. Thus we propose an interpolation function.
#if dimension==2
typedef struct { int x, y;} pseudo_v_i;
#endif
#if dimension==3
typedef struct { int x, y, z;} pseudo_v_i;
#endif
double interpolate_in_cell (Point point, scalar p, coord psi) {
double val;
//shift if required
pseudo_v_i s;
foreach_dimension ()
s.x = -(psi.x < 0);
//redefinition of the relative coordinates
foreach_dimension ()
psi.x -= s.x;
#if dimension == 1
val = p[s.x]*(1. - psi.x)
+ p[s.x+1]*psi.x;
#endif
#if dimension == 2
val = p[s.x,s.y] *(1. - psi.x)*(1. - psi.y)
+ p[s.x+1,s.y] *psi.x*(1. - psi.y)
+ p[s.x,s.y+1] *(1. - psi.x)*psi.y
+ p[s.x+1,s.y+1] *psi.x*psi.y;
#endif
#if dimension == 3
val = p[s.x,s.y,s.z] *(1. - psi.x)*(1. - psi.y)*(1. - psi.z)
+ p[s.x+1,s.y,s.z] *psi.x*(1. - psi.y)*(1. - psi.z)
+ p[s.x,s.y+1,s.z] *(1. - psi.x)*psi.y*(1. - psi.z)
+ p[s.x+1,s.y+1,s.z] *psi.x*psi.y*(1. - psi.z)
+ p[s.x,s.y,s.z+1] *(1. - psi.x)*(1. - psi.y)*psi.z
+ p[s.x+1,s.y,s.z+1] *psi.x*(1. - psi.y)*psi.z
+ p[s.x,s.y+1,s.z+1] *(1. - psi.x)*psi.y*psi.z
+ p[s.x+1,s.y+1,s.z+1] *psi.x*psi.y*psi.z;
#endif
return val;
}
A convenient structure and a convenient function to estimate the interpolated value of a scalar or the difference between the interpolated value of a scalar and a known function at prescribed locations.
struct RealError {
scalar f;
int Lnorm;
double (* func) (double x, double y, double z);
};
double point_error (Point point, struct RealError re, coord psi) {
coord xl;
xl.x = x + Delta*psi.x;
xl.y = y + Delta*psi.y;
xl.z = z + Delta*psi.z;
if (re.func)
return fabs(re.func(xl.x, xl.y, xl.z) - interpolate_in_cell(point, re.f, psi));
else
return fabs(interpolate_in_cell(point, re.f, psi));
}
Some non Gauss integrations
“Child-cell” integrations (integration at the location of the center of child cells).
First, integration of a function.
struct IntegrateFun {
int Lnorm;
double (* func) (double x, double y, double z);
};
double integrate_in_cell (Point point, struct IntegrateFun p, int sublevel) {
int Lnorm = p.Lnorm;
int np = pow(2,sublevel);
double val = 0.;
coord xc;
xc.x = x; xc.y = y; xc.z = z;
#if dimension == 1
for (int i =0; i < np; i++) {
xc.x = x + Delta*(-0.5 + i/((double) np));
val += pow(p.func(xc.x, xc.y, xc.z) ,Lnorm)*dv()/pow(np,dimension);
}
#endif
#if dimension == 2
for (int i =0; i < np; i++)
for (int j =0; j < np; j++) {
xc.x = x + Delta*(-0.5 + i/((double) np));
xc.y = y + Delta*(-0.5 + j/((double) np));
val += pow(p.func(xc.x, xc.y, xc.z) ,Lnorm)*dv()/pow(np,dimension);
}
#endif
#if dimension == 3
for (int i =0; i < np; i++)
for (int j =0; j < np; j++)
for (int k =0; k < np; k++) {
xc.x = x + Delta*(-0.5 + i/((double) np));
xc.y = y + Delta*(-0.5 + j/((double) np));
xc.z = z + Delta*(-0.5 + j/((double) np));
val += pow(p.func(xc.x, xc.y, xc.z) ,Lnorm)*dv()/pow(np,dimension);
}
#endif
return val;
}
double normfromChilds (struct IntegrateFun p, int sublevel)
{
int Lnorm = p.Lnorm;
double val = 0., volume = 0.;
foreach(reduction(+:val) reduction(+:volume))
if (dv() > 0.) {
volume += dv();
val += pow(integrate_in_cell(point, p, sublevel),Lnorm)*dv();
}
return (volume ? pow(val,1./Lnorm) : 0.);
}
“Child”-integration to compute the total error (for scalar (and function)).
double norm_child (struct RealError re)
{
int Lnorm = re.Lnorm;
double val = 0., volume = 0.;
int np = pow(2,dimension);
#if dimension == 1
static coord int_points[2] = {
{-0.25,0.,0.},{0.25,0.,0.}
};
#endif
#if dimension == 2
static coord int_points[4] = {
{-0.25,0.25,0.},{0.25,0.25,0.},{-0.25,-0.25,0.},{0.25,-0.25,0.}
};
#endif
#if dimension == 3
static coord int_points[8] = {
{-0.25,0.25,0.25},{0.25,0.25,0.25},{-0.25,-0.25,0.25},{0.25,-0.25,0.25},
{-0.25,0.25,-0.25},{0.25,0.25,-0.25},{-0.25,-0.25,-0.25},{0.25,-0.25,-0.25}
};
#endif
foreach(reduction(+:val) reduction(+:volume)) {
if (dv() > 0.) {
volume += dv();
for (int i =0; i < np; i++) {
val += pow(point_error(point, re, int_points[i]), Lnorm)*dv()/np;
}
}
}
return (volume ? pow(val/volume,1./Lnorm) : 0.);
}
Gauss-quadratures
3 points Gauss quadrature
Only in 2D for now.
#if dimension == 2
double norm_gauss_3p (struct RealError re)
{
int Lnorm = re.Lnorm;
double val = 0., volume = 0.;
int np = 9;
double p1 = -sqrt(3./5.)/2.;
double p3 = sqrt(3./5.)/2.;
coord int_points[9] = {
{p1,p1,0.},{p1,0.,0.},{p1,p3,0.},
{0.,p1,0.},{0.,0.,0.},{0.,p3,0.},
{p3,p1,0.},{p3,0.,0.},{p3,p3,0.}
};
double coef_points[9] = {
0.3086419753086420,
0.4938271604938272,
0.3086419753086420,
0.4938271604938272,
0.7901234567901235,
0.4938271604938272,
0.3086419753086420,
0.4938271604938272,
0.3086419753086420
};
foreach(reduction(+:val) reduction(+:volume))
if (dv() > 0.) {
volume += dv();
for (int i =0; i < np; i++)
val += pow(point_error(point, re, int_points[i]), Lnorm)*coef_points[i] *dv()/4.;
}
return (volume ? pow(val/volume,1./Lnorm) : 0.);
}
#endif
4 points Gauss quadrature
Only in 2D for now.
#if dimension == 2
double norm_gauss_4p (struct RealError re)
{
int Lnorm = re.Lnorm;
double val = 0., volume = 0.;
int np = 16;
double p1 = sqrt( 3./7. - 2./7.*sqrt(6./5.) ) /2.;
double p2 = -sqrt( 3./7. - 2./7.*sqrt(6./5.) ) /2.;
double p3 = sqrt( 3./7. + 2./7.*sqrt(6./5.) ) /2.;
double p4 = -sqrt( 3./7. + 2./7.*sqrt(6./5.) ) /2.;
coord int_points[16] = {
{p1,p1,0.},{p1,p2,0.},{p1,p3,0.},{p1,p4,0.},
{p2,p1,0.},{p2,p2,0.},{p2,p3,0.},{p2,p4,0.},
{p3,p1,0.},{p3,p2,0.},{p3,p3,0.},{p3,p4,0.},
{p4,p1,0.},{p4,p2,0.},{p4,p3,0.},{p4,p4,0.}
};
double pds12 = (18.+sqrt(30.))/36.;
double pds34 = (18.-sqrt(30.))/36.;
double coef_points[16] = {
pds12*pds12, pds12*pds12, pds12*pds34, pds12*pds34,
pds12*pds12, pds12*pds12, pds12*pds34, pds12*pds34,
pds34*pds12, pds34*pds12, pds34*pds34, pds34*pds34,
pds34*pds12, pds34*pds12, pds34*pds34, pds34*pds34,
};
foreach(reduction(+:val) reduction(+:volume))
if (dv() > 0.) {
volume += dv();
for (int i =0; i < np; i++)
val += pow(point_error(point, re, int_points[i]), Lnorm)*coef_points[i] *dv()/4.;
}
return (volume ? pow(val/volume,1./Lnorm) : 0.);
}
#endif
5 points Gauss quadrature
in 2D and 3D.
double norm_gauss_5p (struct RealError re)
{
int Lnorm = re.Lnorm;
double val = 0., volume = 0.;
int np = 25;
#if dimension == 2
double p1 = 1./3.*sqrt( 5. - 2.*sqrt(10./7.) ) /2.;
double p2 = -1./3.*sqrt( 5. - 2.*sqrt(10./7.) ) /2.;
double p3 = 1./3.*sqrt( 5. + 2.*sqrt(10./7.) ) /2.;
double p4 = -1./3.*sqrt( 5. + 2.*sqrt(10./7.) ) /2.;
double p5 = 0.;
coord int_points[25] = {
{p1,p1,0.},{p1,p2,0.},{p1,p3,0.},{p1,p4,0.},{p1,p5,0.},
{p2,p1,0.},{p2,p2,0.},{p2,p3,0.},{p2,p4,0.},{p2,p5,0.},
{p3,p1,0.},{p3,p2,0.},{p3,p3,0.},{p3,p4,0.},{p3,p5,0.},
{p4,p1,0.},{p4,p2,0.},{p4,p3,0.},{p4,p4,0.},{p4,p5,0.},
{p5,p1,0.},{p5,p2,0.},{p5,p3,0.},{p5,p4,0.},{p5,p5,0.}
};
double pds12 = (322.+13.*sqrt(70.))/900.;
double pds34 = (322.-13.*sqrt(70.))/900.;
double pds5 = 128./225.;
double coef_points[25] = {
pds12*pds12, pds12*pds12, pds12*pds34, pds12*pds34, pds12*pds5,
pds12*pds12, pds12*pds12, pds12*pds34, pds12*pds34, pds12*pds5,
pds34*pds12, pds34*pds12, pds34*pds34, pds34*pds34, pds34*pds5,
pds34*pds12, pds34*pds12, pds34*pds34, pds34*pds34, pds34*pds5,
pds5*pds12, pds5*pds12, pds5*pds34, pds5*pds34, pds5*pds5
};
#endif
#if dimension == 3
double p1 = 1./3.*sqrt( 5. - 2.*sqrt(10./7.) ) /2.;
double p2 = -1./3.*sqrt( 5. - 2.*sqrt(10./7.) ) /2.;
double p3 = 1./3.*sqrt( 5. + 2.*sqrt(10./7.) ) /2.;
double p4 = -1./3.*sqrt( 5. + 2.*sqrt(10./7.) ) /2.;
double p5 = 0.;
coord int_points[125] = {
{p1,p1,p1},{p1,p2,p1},{p1,p3,p1},{p1,p4,p1},{p1,p5,p1},
{p2,p1,p1},{p2,p2,p1},{p2,p3,p1},{p2,p4,p1},{p2,p5,p1},
{p3,p1,p1},{p3,p2,p1},{p3,p3,p1},{p3,p4,p1},{p3,p5,p1},
{p4,p1,p1},{p4,p2,p1},{p4,p3,p1},{p4,p4,p1},{p4,p5,p1},
{p5,p1,p1},{p5,p2,p1},{p5,p3,p1},{p5,p4,p1},{p5,p5,p1},
{p1,p1,p2},{p1,p2,p2},{p1,p3,p2},{p1,p4,p2},{p1,p5,p2},
{p2,p1,p2},{p2,p2,p2},{p2,p3,p2},{p2,p4,p2},{p2,p5,p2},
{p3,p1,p2},{p3,p2,p2},{p3,p3,p2},{p3,p4,p2},{p3,p5,p2},
{p4,p1,p2},{p4,p2,p2},{p4,p3,p2},{p4,p4,p2},{p4,p5,p2},
{p5,p1,p2},{p5,p2,p2},{p5,p3,p2},{p5,p4,p2},{p5,p5,p2},
{p1,p1,p3},{p1,p2,p3},{p1,p3,p3},{p1,p4,p3},{p1,p5,p3},
{p2,p1,p3},{p2,p2,p3},{p2,p3,p3},{p2,p4,p3},{p2,p5,p3},
{p3,p1,p3},{p3,p2,p3},{p3,p3,p3},{p3,p4,p3},{p3,p5,p3},
{p4,p1,p3},{p4,p2,p3},{p4,p3,p3},{p4,p4,p3},{p4,p5,p3},
{p5,p1,p3},{p5,p2,p3},{p5,p3,p3},{p5,p4,p3},{p5,p5,p3},
{p1,p1,p4},{p1,p2,p4},{p1,p3,p4},{p1,p4,p4},{p1,p5,p4},
{p2,p1,p4},{p2,p2,p4},{p2,p3,p4},{p2,p4,p4},{p2,p5,p4},
{p3,p1,p4},{p3,p2,p4},{p3,p3,p4},{p3,p4,p4},{p3,p5,p4},
{p4,p1,p4},{p4,p2,p4},{p4,p3,p4},{p4,p4,p4},{p4,p5,p4},
{p5,p1,p4},{p5,p2,p4},{p5,p3,p4},{p5,p4,p4},{p5,p5,p4},
{p1,p1,p5},{p1,p2,p5},{p1,p3,p5},{p1,p4,p5},{p1,p5,p5},
{p2,p1,p5},{p2,p2,p5},{p2,p3,p5},{p2,p4,p5},{p2,p5,p5},
{p3,p1,p5},{p3,p2,p5},{p3,p3,p5},{p3,p4,p5},{p3,p5,p5},
{p4,p1,p5},{p4,p2,p5},{p4,p3,p5},{p4,p4,p5},{p4,p5,p5},
{p5,p1,p5},{p5,p2,p5},{p5,p3,p5},{p5,p4,p5},{p5,p5,p5}
};
double pds12 = (322.+13.*sqrt(70.))/900.;
double pds34 = (322.-13.*sqrt(70.))/900.;
double pds5 = 128./225.;
double coef_points[125] = {
pds12*pds12*pds12, pds12*pds12*pds12, pds12*pds34*pds12, pds12*pds34*pds12, pds12*pds5*pds12,
pds12*pds12*pds12, pds12*pds12*pds12, pds12*pds34*pds12, pds12*pds34*pds12, pds12*pds5*pds12,
pds34*pds12*pds12, pds34*pds12*pds12, pds34*pds34*pds12, pds34*pds34*pds12, pds34*pds5*pds12,
pds34*pds12*pds12, pds34*pds12*pds12, pds34*pds34*pds12, pds34*pds34*pds12, pds34*pds5*pds12,
pds5*pds12*pds12, pds5*pds12*pds12, pds5*pds34*pds12, pds5*pds34*pds12, pds5*pds5*pds12,
pds12*pds12*pds12, pds12*pds12*pds12, pds12*pds34*pds12, pds12*pds34*pds12, pds12*pds5*pds12,
pds12*pds12*pds12, pds12*pds12*pds12, pds12*pds34*pds12, pds12*pds34*pds12, pds12*pds5*pds12,
pds34*pds12*pds12, pds34*pds12*pds12, pds34*pds34*pds12, pds34*pds34*pds12, pds34*pds5*pds12,
pds34*pds12*pds12, pds34*pds12*pds12, pds34*pds34*pds12, pds34*pds34*pds12, pds34*pds5*pds12,
pds5*pds12*pds12, pds5*pds12*pds12, pds5*pds34*pds12, pds5*pds34*pds12, pds5*pds5*pds12,
pds12*pds12*pds34, pds12*pds12*pds34, pds12*pds34*pds34, pds12*pds34*pds34, pds12*pds5*pds34,
pds12*pds12*pds34, pds12*pds12*pds34, pds12*pds34*pds34, pds12*pds34*pds34, pds12*pds5*pds34,
pds34*pds12*pds34, pds34*pds12*pds34, pds34*pds34*pds34, pds34*pds34*pds34, pds34*pds5*pds34,
pds34*pds12*pds34, pds34*pds12*pds34, pds34*pds34*pds34, pds34*pds34*pds34, pds34*pds5*pds34,
pds5*pds12*pds34, pds5*pds12*pds34, pds5*pds34*pds34, pds5*pds34*pds34, pds5*pds5*pds34,
pds12*pds12*pds34, pds12*pds12*pds34, pds12*pds34*pds34, pds12*pds34*pds34, pds12*pds5*pds34,
pds12*pds12*pds34, pds12*pds12*pds34, pds12*pds34*pds34, pds12*pds34*pds34, pds12*pds5*pds34,
pds34*pds12*pds34, pds34*pds12*pds34, pds34*pds34*pds34, pds34*pds34*pds34, pds34*pds5*pds34,
pds34*pds12*pds34, pds34*pds12*pds34, pds34*pds34*pds34, pds34*pds34*pds34, pds34*pds5*pds34,
pds5*pds12*pds34, pds5*pds12*pds34, pds5*pds34*pds34, pds5*pds34*pds34, pds5*pds5*pds34,
pds12*pds12*pds5, pds12*pds12*pds5, pds12*pds34*pds5, pds12*pds34*pds5, pds12*pds5*pds5,
pds12*pds12*pds5, pds12*pds12*pds5, pds12*pds34*pds5, pds12*pds34*pds5, pds12*pds5*pds5,
pds34*pds12*pds5, pds34*pds12*pds5, pds34*pds34*pds5, pds34*pds34*pds5, pds34*pds5*pds5,
pds34*pds12*pds5, pds34*pds12*pds5, pds34*pds34*pds5, pds34*pds34*pds5, pds34*pds5*pds5,
pds5*pds12*pds5, pds5*pds12*pds5, pds5*pds34*pds5, pds5*pds34*pds5, pds5*pds5*pds5
};
#endif
foreach(reduction(+:val) reduction(+:volume))
if (dv() > 0.) {
volume += dv();
for (int i =0; i < np; i++)
val += pow(point_error(point, re, int_points[i]), Lnorm)*coef_points[i] *dv()/pow(2,dimension); // the sum of the weight coef_points is 2^dimension
}
return (volume ? pow(val/volume,1./Lnorm) : 0.);
}
5 points Gauss quadrature whith local error registration
Only in 2D for now.
#if dimension == 2
double norm_gauss_5p_local (struct RealError re, scalar local_error)
{
int Lnorm = re.Lnorm;
double val = 0., volume = 0.;
int np = 25;
double p1 = 1./3.*sqrt( 5. - 2.*sqrt(10./7.) ) /2.;
double p2 = -1./3.*sqrt( 5. - 2.*sqrt(10./7.) ) /2.;
double p3 = 1./3.*sqrt( 5. + 2.*sqrt(10./7.) ) /2.;
double p4 = -1./3.*sqrt( 5. + 2.*sqrt(10./7.) ) /2.;
double p5 = 0.;
coord int_points[25] = {
{p1,p1,0.},{p1,p2,0.},{p1,p3,0.},{p1,p4,0.},{p1,p5,0.},
{p2,p1,0.},{p2,p2,0.},{p2,p3,0.},{p2,p4,0.},{p2,p5,0.},
{p3,p1,0.},{p3,p2,0.},{p3,p3,0.},{p3,p4,0.},{p3,p5,0.},
{p4,p1,0.},{p4,p2,0.},{p4,p3,0.},{p4,p4,0.},{p4,p5,0.},
{p5,p1,0.},{p5,p2,0.},{p5,p3,0.},{p5,p4,0.},{p5,p5,0.}
};
double pds12 = (322.+13.*sqrt(70.))/900.;
double pds34 = (322.-13.*sqrt(70.))/900.;
double pds5 = 128./225.;
double coef_points[25] = {
pds12*pds12, pds12*pds12, pds12*pds34, pds12*pds34, pds12*pds5,
pds12*pds12, pds12*pds12, pds12*pds34, pds12*pds34, pds12*pds5,
pds34*pds12, pds34*pds12, pds34*pds34, pds34*pds34, pds34*pds5,
pds34*pds12, pds34*pds12, pds34*pds34, pds34*pds34, pds34*pds5,
pds5*pds12, pds5*pds12, pds5*pds34, pds5*pds34, pds5*pds5
};
foreach(reduction(+:val) reduction(+:volume)){
local_error[]=0.;
if (dv() > 0.) {
volume += dv();
for (int i =0; i < np; i++) {
val += pow(point_error(point, re, int_points[i]), Lnorm)*coef_points[i] *dv()/4.;
local_error[] += pow(point_error(point, re, int_points[i]), Lnorm)*coef_points[i] *dv()/4.;
}
}
}
return (volume ? pow(val/volume,1./Lnorm) : 0.);
}
#endif