Skip to content
Snippets Groups Projects
Commit 821f27b5 authored by Ruth Sabariego's avatar Ruth Sabariego
Browse files

correction grad basis functions and Gauss points for pyramids

parent 57b5491b
No related branches found
No related tags found
No related merge requests found
......@@ -830,22 +830,24 @@ public:
inline int getNumGaussPoints(){ return 8; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
double u8[8] = {0.3595161057791018,0.09633205020967324,
0.3595161057791018,0.09633205020967324,
0.6920507403468987,0.1854344369976602,
0.6920507403468987,0.1854344369976602};
double v8[8] = {0.3595161057791018,0.3595161057791018,
0.09633205020967324,0.09633205020967324,
0.6920507403468987,0.6920507403468987,
0.1854344369976602,0.1854344369976602};
//Integration points for old_pyramids
double u8[8] = {0.2631840555694285,-0.2631840555694285,
0.2631840555694285,-0.2631840555694285,
0.5066163033492386,-0.5066163033492386,
0.5066163033492386,-0.5066163033492386};
double v8[8] = {0.2631840555694285,0.2631840555694285,
-0.2631840555694285,-0.2631840555694285,
0.5066163033492386,0.5066163033492386,
-0.5066163033492386,-0.5066163033492386};
double w8[8] = {0.544151844011225,0.544151844011225,
0.544151844011225,0.544151844011225,
0.122514822655441,0.122514822655441,
0.122514822655441,0.122514822655441};
double p8[8] = {0.02519647051995625,0.02519647051995625,
0.02519647051995625,0.02519647051995625,
0.058136862813377,0.058136862813377,
0.058136862813377,0.058136862813377};
double p8[8] = {0.100785882079825,0.100785882079825,
0.100785882079825,0.100785882079825,
0.232547451253508,0.232547451253508,
0.232547451253508,0.232547451253508};
if(num < 0 || num > 7) return;
u = u8[num];
v = v8[num];
......@@ -855,6 +857,7 @@ public:
void getShapeFunction(int num, double u, double v, double w, double &s)
{
double r;
if(w != 1. && num != 4) r = u*v*w / (1.-w);
else r = 0.;
switch(num) {
......@@ -869,24 +872,39 @@ public:
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
{
if(w == 1. && num != 4) {
s[0] = 0.25;
switch(num) {
case 0 : s[0] = -0.25 ;
s[1] = -0.25 ;
s[2] = -0.25 ; break ;
case 1 : s[0] = 0.25 ;
s[1] = -0.25 ;
s[2] = -0.25 ; break ;
case 2 : s[0] = 0.25 ;
s[1] = 0.25 ;
s[2] = -0.25;
s[2] = -0.25 ; break ;
case 3 : s[0] = -0.25 ;
s[1] = 0.25 ;
s[2] = -0.25 ; break ;
case 4 : s[0] = 0. ;
s[1] = 0. ;
s[2] = 1. ; break ;
default : s[0] = s[1] = s[2] = 0.; break;
}
}
else{
switch(num) {
case 0 : s[0] = 0.25 * ( -(1.-v) + v*w/(1.-w) ) ;
s[1] = 0.25 * ( -(1.-u) + u*w/(1.-w) ) ;
s[2] = 0.25 * ( -1. + u*v/SQU(1.-w) ) ; break ;
case 1 : s[0] = 0.25 * ( (1.-v) + v*w/(1.-w) ) ;
s[1] = 0.25 * ( -(1.+u) + u*w/(1.-w) ) ;
s[2] = 0.25 * ( -1. + u*v/SQU(1.-w) ) ; break ;
s[2] = 0.25 * ( -1. + u*v/(1.-w) + u*v*w/SQU(1.-w) ) ; break ;
case 1 : s[0] = 0.25 * ( (1.-v) - v*w/(1.-w) ) ;
s[1] = 0.25 * ( -(1.+u) - u*w/(1.-w) ) ;
s[2] = 0.25 * ( -1. - u*v/(1.-w) - u*v*w/SQU(1.-w) ) ; break ;
case 2 : s[0] = 0.25 * ( (1.+v) + v*w/(1.-w) ) ;
s[1] = 0.25 * ( (1.+u) + u*w/(1.-w) ) ;
s[2] = 0.25 * ( -1. + u*v/SQU(1.-w) ) ; break ;
case 3 : s[0] = 0.25 * ( -(1.+v) + v*w/(1.-w) ) ;
s[1] = 0.25 * ( (1.-u) + u*w/(1.-w) ) ;
s[2] = 0.25 * ( -1. + u*v/SQU(1.-w) ) ; break ;
s[2] = 0.25 * ( -1. + u*v/(1.-w) + u*v*w/SQU(1.-w) ) ; break ;
case 3 : s[0] = 0.25 * ( -(1.+v) - v*w/(1.-w) ) ;
s[1] = 0.25 * ( (1.-u) - u*w/(1.-w) ) ;
s[2] = 0.25 * ( -1. - u*v/(1.-w) - u*v*w/SQU(1.-w) ) ; break ;
case 4 : s[0] = 0. ;
s[1] = 0. ;
s[2] = 1. ; break ;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment