is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.
| Gauss_Hexahedron.cpp (getdp-3.4.0-source.tgz) | : | Gauss_Hexahedron.cpp (getdp-3.5.0-source.tgz) |
|---|
|
| // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege | | // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege |
| // | | // |
| // See the LICENSE.txt file for license information. Please report all | | // See the LICENSE.txt file for license information. Please report all |
| // 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_Hexahedron.h" | | #include "Gauss_Hexahedron.h" |
| #include "Message.h" | | #include "Message.h" |
| #include "MallocUtils.h" | | #include "MallocUtils.h" |
| | |
| /* Gauss integration over a hexahedron */ | | /* Gauss integration over a hexahedron */ |
| | |
|
| void Gauss_Hexahedron(int Nbr_Points, int Num, | | void Gauss_Hexahedron(int Nbr_Points, int Num, double *u, double *v, double *w, |
| double *u, double *v, double *w, double *wght) | | double *wght) |
| { | | { |
|
| switch (Nbr_Points) { | | switch(Nbr_Points) { |
| | | case 6: |
| case 6 : | | *u = xhex6[Num]; |
| *u = xhex6 [Num] ; *v = yhex6 [Num] ; *w = zhex6 [Num] ; | | *v = yhex6[Num]; |
| *wght = phex6 [Num] ; break ; | | *w = zhex6[Num]; |
| | *wght = phex6[Num]; |
| | break; |
| | |
|
| case 14 : | | case 14: |
| *u = xhex14 [Num] ; *v = yhex14 [Num] ; *w = zhex14 [Num] ; | | *u = xhex14[Num]; |
| *wght = phex14 [Num] ; break ; | | *v = yhex14[Num]; |
| | *w = zhex14[Num]; |
| | *wght = phex14[Num]; |
| | break; |
| | |
|
| case 34 : | | case 34: |
| *u = xhex34[Num] ; *v = yhex34[Num] ; *w = zhex34[Num] ; | | *u = xhex34[Num]; |
| *wght = phex34[Num] ; break ; | | *v = yhex34[Num]; |
| | *w = zhex34[Num]; |
| | *wght = phex34[Num]; |
| | break; |
| | |
|
| case 77 : | | case 77: |
| *u = xhex77[Num] ; *v = yhex77[Num] ; *w = zhex77[Num] ; | | *u = xhex77[Num]; |
| *wght = phex77[Num] ; break ; | | *v = yhex77[Num]; |
| | *w = zhex77[Num]; |
| | *wght = phex77[Num]; |
| | break; |
| | |
|
| default : | | default: |
| Message::Error("Wrong number of Gauss points for Hexahedron: " | | Message::Error("Wrong number of Gauss points for Hexahedron: " |
| "valid choices: 6, 14, 34, 77"); | | "valid choices: 6, 14, 34, 77"); |
| break; | | break; |
|
| | | |
| } | | } |
| } | | } |
| | |
| /* Gauss-Legendre scheme to integrate over a hexahedron */ | | /* Gauss-Legendre scheme to integrate over a hexahedron */ |
| | |
| static int glhex[MAX_LINE_POINTS] = {-1}; | | static int glhex[MAX_LINE_POINTS] = {-1}; |
| static double *glxhex[MAX_LINE_POINTS], *glyhex[MAX_LINE_POINTS]; | | static double *glxhex[MAX_LINE_POINTS], *glyhex[MAX_LINE_POINTS]; |
| static double *glzhex[MAX_LINE_POINTS], *glphex[MAX_LINE_POINTS]; | | static double *glzhex[MAX_LINE_POINTS], *glphex[MAX_LINE_POINTS]; |
| | |
|
| void GaussLegendre_Hexahedron(int Nbr_Points, int Num, | | void GaussLegendre_Hexahedron(int Nbr_Points, int Num, double *u, double *v, |
| double *u, double *v, double *w, double *wght) | | double *w, double *wght) |
| { | | { |
| int i, j, k, index = 0, nb; | | int i, j, k, index = 0, nb; |
| double pt1, pt2, pt3, wt1, wt2, wt3, dum; | | double pt1, pt2, pt3, wt1, wt2, wt3, dum; |
| | |
|
| nb = (int)cbrt((double)Nbr_Points); | | nb = (int)(cbrt((double)Nbr_Points) + 0.5); |
| | |
|
| if(nb * nb * nb != Nbr_Points || nb > MAX_LINE_POINTS){ | | if(nb * nb * nb != Nbr_Points || nb > MAX_LINE_POINTS) { |
| Message::Error("Number of points should be n^3 with n in [1,%d]", MAX_LINE_P | | Message::Error("Number of points should be n^3 with n in [1,%d]", |
| OINTS) ; | | MAX_LINE_POINTS); |
| return; | | return; |
| } | | } |
| | |
|
| if(glhex[0] < 0) for(i = 0; i < MAX_LINE_POINTS; i++) glhex[i] = 0 ; | | if(glhex[0] < 0) |
| | for(i = 0; i < MAX_LINE_POINTS; i++) glhex[i] = 0; |
| | |
|
| if(!glhex[nb - 1]){ | | if(!glhex[nb - 1]) { |
| Message::Info("Computing GaussLegendre %dx%dx%d for Hexahedron", nb, nb, nb) | | Message::Info("Computing GaussLegendre %dx%dx%d for Hexahedron", nb, nb, |
| ; | | nb); |
| glxhex[nb - 1] = (double*)Malloc(Nbr_Points * sizeof(double)); | | glxhex[nb - 1] = (double *)Malloc(Nbr_Points * sizeof(double)); |
| glyhex[nb - 1] = (double*)Malloc(Nbr_Points * sizeof(double)); | | glyhex[nb - 1] = (double *)Malloc(Nbr_Points * sizeof(double)); |
| glzhex[nb - 1] = (double*)Malloc(Nbr_Points * sizeof(double)); | | glzhex[nb - 1] = (double *)Malloc(Nbr_Points * sizeof(double)); |
| glphex[nb - 1] = (double*)Malloc(Nbr_Points * sizeof(double)); | | glphex[nb - 1] = (double *)Malloc(Nbr_Points * sizeof(double)); |
| for(i = 0; 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); |
| for(j = 0; 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); |
| for(k = 0; 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); |
| glxhex[nb - 1][index] = pt1; | | glxhex[nb - 1][index] = pt1; |
| glyhex[nb - 1][index] = pt2; | | glyhex[nb - 1][index] = pt2; |
| glzhex[nb - 1][index] = pt3; | | glzhex[nb - 1][index] = pt3; |
| glphex[nb - 1][index++] = wt1 * wt2 * wt3; | | glphex[nb - 1][index++] = wt1 * wt2 * wt3; |
| } | | } |
| } | | } |
| } | | } |
| glhex[nb - 1] = 1; | | glhex[nb - 1] = 1; |
| } | | } |
| | |
|
| *u = glxhex[nb - 1][Num] ; *v = glyhex[nb - 1][Num] ; *w = glzhex[nb - 1][Num] | | *u = glxhex[nb - 1][Num]; |
| ; | | *v = glyhex[nb - 1][Num]; |
| *wght = glphex[nb - 1][Num] ; | | *w = glzhex[nb - 1][Num]; |
| | *wght = glphex[nb - 1][Num]; |
| } | | } |
| | | |
| End of changes. 15 change blocks. |
|---|
| 44 lines changed or deleted | | 56 lines changed or added |
|---|
"Fossies" - the Fresh Open Source Software Archive 