"Fossies" - the Fresh Open Source Software Archive

Source code changes of the file "Kernel/Gauss_Line.cpp" betweengetdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

Gauss_Line.cpp  (getdp-3.4.0-source.tgz):Gauss_Line.cpp  (getdp-3.5.0-source.tgz)
// GetDP - Copyright (C) 1997-202 P. Dular and C. Geuzaine, University of Liege // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege
// //
// issues on https://gitlab.onelab.info/getdp/getdp/issues. // issues on https://gitlab.onelab.info/getdp/getdp/issues.
#include <math.h> #include <math.h>
#include "Gauss.h" #include "Gauss.h"
#include "Gauss_Line.h" #include "Gauss_Line.h"
#include "Message.h" #include "Message.h"
#include "MallocUtils.h" #include "MallocUtils.h"
/* Gauss integration over a line */ /* Gauss integration over a line */
static int gll[MAX_LINE_POINTS] = {-1}; static int gll[MAX_LINE_POINTS] = {-1};
static double *glxl[MAX_LINE_POINTS], *glpl[MAX_LINE_POINTS]; static double *glxl[MAX_LINE_POINTS], *glpl[MAX_LINE_POINTS];
void Gauss_Line(int Nbr_Points, int Num, void Gauss_Line(int Nbr_Points, int Num, double *u, double *v, double *w,
double *u, double *v, double *w, double *wght) double *wght)
{ {
int i; int i;
{ switch(Nbr_Points) {
case 1:
case *u = *v = *w = *wght = *u = lx1[Num];
case *u = *v = *w = *wght = *v = 0.;
case *u = *v = *w = *wght = *w = 0.;
case *u = *v = *w = *wght = *wght = lp1[Num];
case *u = *v = *w = *wght = break;
case *u = *v = *w = *wght = case 2:
case *u = *v = *w = *wght = *u = lx2[Num];
case *u = *v = *w = *wght = *v = 0.;
case *u = *v = *w = *wght = *w = 0.;
case *u = *v = *w = *wght = *wght = lp2[Num];
case *u = *v = *w = *wght = break;
case *u = *v = *w = *wght = case 3:
case *u = *v = *w = *wght = *u = lx3[Num];
case *u = *v = *w = *wght = *v = 0.;
case *u = *v = *w = *wght = *w = 0.;
case *u = *v = *w = *wght = *wght = lp3[Num];
case *u = *v = *w = *wght = break;
case *u = *v = *w = *wght = case 4:
case *u = *v = *w = *wght = *u = lx4[Num];
case *u = *v = *w = *wght = *v = 0.;
*w = 0.;
if(Nbr_Points <= *wght = lp4[Num];
if(gll[0] < 0) for(i = 0; i < MAX_LINE_POINTS; i++) gll[i] = break;
if(!gll[Nbr_Points - case 5:
Message::Info("Computing GaussLegendre %d for Line", Nbr_Points); *u = lx5[Num];
glxl[Nbr_Points - 1] = * sizeof(double)); *v = 0.;
glpl[Nbr_Points - 1] = * sizeof(double)); *w = 0.;
GaussLegendre(-1., 1., glxl[Nbr_Points - 1], glpl[Nbr_Points - 1], *wght = lp5[Num];
break;
gll[Nbr_Points - 1] = 1; case 6:
*u = lx6[Num];
*v = 0.;
*w = 0.;
*wght = lp6[Num];
break;
case 7:
*u = lx7[Num];
*v = 0.;
*w = 0.;
*wght = lp7[Num];
break;
case 8:
*u = lx8[Num];
*v = 0.;
*w = 0.;
*wght = lp8[Num];
break;
case 9:
*u = lx9[Num];
*v = 0.;
*w = 0.;
*wght = lp9[Num];
break;
case 10:
*u = lx10[Num];
*v = 0.;
*w = 0.;
*wght = lp10[Num];
break;
case 11:
*u = lx11[Num];
*v = 0.;
*w = 0.;
*wght = lp11[Num];
break;
case 12:
*u = lx12[Num];
*v = 0.;
*w = 0.;
*wght = lp12[Num];
break;
case 13:
*u = lx13[Num];
*v = 0.;
*w = 0.;
*wght = lp13[Num];
break;
case 14:
*u = lx14[Num];
*v = 0.;
*w = 0.;
*wght = lp14[Num];
break;
case 15:
*u = lx15[Num];
*v = 0.;
*w = 0.;
*wght = lp15[Num];
break;
case 16:
*u = lx16[Num];
*v = 0.;
*w = 0.;
*wght = lp16[Num];
break;
case 17:
*u = lx17[Num];
*v = 0.;
*w = 0.;
*wght = lp17[Num];
break;
case 18:
*u = lx18[Num];
*v = 0.;
*w = 0.;
*wght = lp18[Num];
break;
case 19:
*u = lx19[Num];
*v = 0.;
*w = 0.;
*wght = lp19[Num];
break;
case 20:
*u = lx20[Num];
*v = 0.;
*w = 0.;
*wght = lp20[Num];
break;
default:
if(Nbr_Points <= MAX_LINE_POINTS) {
if(gll[0] < 0)
for(i = 0; i < MAX_LINE_POINTS; i++) gll[i] = 0;
if(!gll[Nbr_Points - 1]) {
Message::Info("Computing GaussLegendre %d for Line", Nbr_Points);
glxl[Nbr_Points - 1] = (double *)Malloc(Nbr_Points * sizeof(double));
glpl[Nbr_Points - 1] = (double *)Malloc(Nbr_Points * sizeof(double));
GaussLegendre(-1., 1., glxl[Nbr_Points - 1], glpl[Nbr_Points - 1],
Nbr_Points);
gll[Nbr_Points - 1] = 1;
} }
*u = glxl[Nbr_Points - *v = *w = *wght = glpl[Nbr_Points - *u = glxl[Nbr_Points - 1][Num];
*v = *w = 0.;
*wght = glpl[Nbr_Points - 1][Num];
} }
else else
Message::Error("Maximum number of integration points exceeded (%d > %d)", Message::Error("Maximum number of integration points exceeded (%d > %d)",
Nbr_Points, Nbr_Points, MAX_LINE_POINTS);
break;
} }
} }
#define EPS 3.0e-11 #define EPS 3.0e-11
void GaussLegendre(double x1, double x2, double x[], double w[], int n) void GaussLegendre(double x1, double x2, double x[], double w[], int n)
{ {
int m, j, i; int m, j, i;
double z1, z, xm, xl, pp, p3, p2, p1; double z1, z, xm, xl, pp, p3, p2, p1;
m = (n + 1) / 2; m = (n + 1) / 2;
xm = 0.5 * (x2 + x1); xm = 0.5 * (x2 + x1);
xl = 0.5 * (x2 - x1); xl = 0.5 * (x2 - x1);
for(i = 1; i <= m; i++) { for(i = 1; i <= m; i++) {
z = cos(3.141592654 * (i - 0.25) / (n + 0.5)); z = cos(3.141592654 * (i - 0.25) / (n + 0.5));
do { do {
p1 = 1.0; p1 = 1.0;
p2 = 0.0; p2 = 0.0;
= 1; j <= n; j++) { for(j = 1; j <= n; j++) {
p3 = p2; p3 = p2;
p2 = p1; p2 = p1;
p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / j; p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / j;
} }
pp = n * (z * p1 - p2) / (z * z - 1.0); pp = n * (z * p1 - p2) / (z * z - 1.0);
z1 = z; z1 = z;
z = z1 - p1 / pp; z = z1 - p1 / pp;
} while(fabs(z - z1) > EPS); } while(fabs(z - z1) > EPS);
x[i - 1] = xm - xl * z; x[i - 1] = xm - xl * z;
x[n - i] = xm + xl * z; x[n - i] = xm + xl * z;
w[i - 1] = 2.0 * xl((1.0 - z * z) * pp * pp); w[i - 1] = 2.0 * xl / ((1.0 - z * z) * pp * pp);
w[n - i] = w[i - 1]; w[n - i] = w[i - 1];
} }
} }
End of changes. 10 change blocks.
47 lines changed or deleted 148 lines changed or added