"Fossies" - the Fresh Open Source Software Archive

Source code changes of the file "Kernel/Gauss_Tetrahedron.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_Tetrahedron.cpp  (getdp-3.4.0-source.tgz):Gauss_Tetrahedron.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_Tetrahedron.h" #include "Gauss_Tetrahedron.h"
#include "Message.h" #include "Message.h"
#include "MallocUtils.h" #include "MallocUtils.h"
/* Gauss integration over a tetrahedron */ /* Gauss integration over a tetrahedron */
void Gauss_Tetrahedron(int Nbr_Points, int Num, void Gauss_Tetrahedron(int Nbr_Points, int Num, double *u, double *v, double *w,
double *u, double *v, double *w, double *wght) double *wght)
{ {
{ switch(Nbr_Points) {
case 1:
case *u = xtet1[Num];
*u = *v = *w = *v = ytet1[Num];
*wght = *w = ztet1[Num];
*wght = ptet1[Num];
break;
case case 4:
*u = *v = *w = *u = xtet4[Num];
*wght = *v = ytet4[Num];
*w = ztet4[Num];
*wght = ptet4[Num];
break;
case case 5:
*u = *v = *w = *u = xtet5[Num];
*wght = *v = ytet5[Num];
*w = ztet5[Num];
*wght = ptet5[Num];
break;
case case 15:
*u = *v = *w = *u = xtet15[Num];
*wght = *v = ytet15[Num];
*w = ztet15[Num];
*wght = ptet15[Num];
break;
case case 16:
*u = *v = *w = *u = xtet16[Num];
*wght = *v = ytet16[Num];
*w = ztet16[Num];
*wght = ptet16[Num];
break;
case case 17:
*u = *v = *w = *u = xtet17[Num];
*wght = *v = ytet17[Num];
*w = ztet17[Num];
*wght = ptet17[Num];
break;
case case 29:
*u = *v = *w = *u = xtet29[Num];
*wght = *v = ytet29[Num];
*w = ztet29[Num];
*wght = ptet29[Num];
break;
default: default:
Message::Error("Wrong number of Gauss Points for Tetrahedron: " Message::Error("Wrong number of Gauss Points for Tetrahedron: "
"valid choices: 1, 4, 5, 15, 16, 17, 29"); "valid choices: 1, 4, 5, 15, 16, 17, 29");
break; break;
} }
} }
/* Degenerate n1Xn2Xn3 Gauss-Legendre scheme to integrate over a tet */ /* Degenerate n1Xn2Xn3 Gauss-Legendre scheme to integrate over a tet */
static int gltet[MAX_LINE_POINTS] = {-1}; static int gltet[MAX_LINE_POINTS] = {-1};
static double *glxtet[MAX_LINE_POINTS], *glytet[MAX_LINE_POINTS]; static double *glxtet[MAX_LINE_POINTS], *glytet[MAX_LINE_POINTS];
static double *glztet[MAX_LINE_POINTS], *glptet[MAX_LINE_POINTS]; static double *glztet[MAX_LINE_POINTS], *glptet[MAX_LINE_POINTS];
void hexToTet(double eta, double zeta, void hexToTet(double xi, double eta, double zeta, double *r, double *s,
double *r, double *s, double *t, double *J) double *t, double *J)
{ {
double r1,rs1; double r1, rs1;
*r = *r = 0.5e0 * (1.0e0 + xi);
r1 = r1 = 1.0e0 - (*r);
*s = *s = 0.5e0 * (1.0e0 + eta) * r1;
rs1 = rs1 = 1.0e0 - (*r) - (*s);
*t = *t = 0.5e0 * (1.0e0 + zeta) * rs1;
*J = *J = 0.125e0 * r1 * rs1;
} }
void GaussLegendre_Tetrahedron(int Nbr_Points, int Num, void GaussLegendre_Tetrahedron(int Nbr_Points, int Num, double *u, double *v,
double *u, double *v, double *w, double *wght) double *w, double *wght)
{ {
int int i, j, k, index = 0, nb;
double double pt1, pt2, pt3, wt1, wt2, wt3, dJ, dum;
nb = (int)); nb = (int)(cbrt((double)Nbr_Points) + 0.5);
!= Nbr_Points || nb > if(nb * nb * nb != Nbr_Points || nb > MAX_LINE_POINTS) {
Message::Error("Number of points should be n^3 with n in [1,%d]", Message::Error("Number of points should be n^3 with n in [1,%d]",
MAX_LINE_POINTS);
return; return;
} }
if(gltet[0] < 0) i < i++) gltet[i] = if(gltet[0] < 0)
for(i = 0; i < MAX_LINE_POINTS; i++) gltet[i] = 0;
if(!gltet[nb{ if(!gltet[nb - 1]) {
Message::Info("Computing degenerate GaussLegendre %dX%dX%d for Tetrahedron", Message::Info("Computing degenerate GaussLegendre %dX%dX%d for Tetrahedron",
nb, nb, nb); nb, nb, nb);
= glxtet[nb - 1] = (double *)Malloc(Nbr_Points * sizeof(double));
= glytet[nb - 1] = (double *)Malloc(Nbr_Points * sizeof(double));
= glztet[nb - 1] = (double *)Malloc(Nbr_Points * sizeof(double));
= glptet[nb - 1] = (double *)Malloc(Nbr_Points * sizeof(double));
i < nb; i++) { for(i = 0; i < nb; i++) {
Gauss_Line(nb, i, &pt1, &dum, &dum, &wt1); Gauss_Line(nb, i, &pt1, &dum, &dum, &wt1);
j < nb; j++) { for(j = 0; j < nb; j++) {
Gauss_Line(nb, j, &pt2, &dum, &dum, &wt2); Gauss_Line(nb, j, &pt2, &dum, &dum, &wt2);
k < nb; k++) { for(k = 0; k < nb; k++) {
Gauss_Line(nb, k, &pt3, &dum, &dum, &wt3); Gauss_Line(nb, k, &pt3, &dum, &dum, &wt3);
hexToTet(pt1, pt2, pt3, hexToTet(pt1, pt2, pt3, &glxtet[nb - 1][index],
&dJ); &glytet[nb - 1][index], &glztet[nb - 1][index], &dJ);
= glptet[nb - 1][index++] = dJ * wt1 * wt2 * wt3;
} }
} }
} }
gltet[nb1] = 1; gltet[nb - 1] = 1;
} }
*u = *v = *w = *u = glxtet[nb - 1][Num];
*wght = *v = glytet[nb - 1][Num];
*w = glztet[nb - 1][Num];
*wght = glptet[nb - 1][Num];
} }
End of changes. 24 change blocks.
63 lines changed or deleted 86 lines changed or added