1 #include "shapeFuncInternals.h" 2 3 #ifdef __cplusplus 4 extern "C" { 5 #endif 6 7 /* Calculate shape functions and derivative for tet elements */ 8 int TetShapeAndDrv(int p,double par[3],double N[],double dN[][3]) { 9 static int TetEMAP[6][2]={{0,1},{1,2},{2,0},{0,3},{1,3},{2,3}}; 10 static int TetFMAP[4][3]={{0,1,2},{0,3,1},{1,3,2},{0,2,3}}; 11 int NS,is,i,j,ip,nshp=0,(*fpdrv)[2]; 12 double L[4],Le[2],Lf[3],mode,blend,bdrv[3],mdrv[3],epdrv[3][2],mfdrv[2]; 13 double tmp,rst,rsw,stw,rstw,rtw; 14 if(p<1) 15 return nshp; 16 L[0]=par[0]; 17 L[1]=par[1]; 18 L[2]=par[2]; 19 L[3]=1.0e0-par[0]-par[1]-par[2]; 20 /* collect all vertex modes */ 21 for(i=0; i <4; i++) { 22 N[i] = L[i]; 23 if(i==3) 24 dN[i][0]=dN[i][1]=dN[i][2]=-1.0e0; 25 else { 26 for(j=0; j <3; j++) { 27 if(i==j) 28 dN[i][j] = 1.0e0; 29 else 30 dN[i][j] = 0.0e0; 31 } 32 } 33 } 34 nshp=4; 35 if( p > 1 ) { 36 /* collect all edge modes, for all p */ 37 for(i=0; i <6; i++) { 38 Le[0]=L[TetEMAP[i][0]]; 39 Le[1]=L[TetEMAP[i][1]]; 40 blend = R_edgeBlendTet(TetEMAP[i],L); 41 R_edgeBlendTetDrv(TetEMAP[i],L,bdrv); 42 E_parDrv(TetEMAP[i][0],TetEMAP[i][1],Stet,epdrv); 43 for(ip=2; ip <= p; ip++) { 44 mode = E_modeShape(ip,Le); 45 E_modeShapeDrv(ip, Le, mdrv); 46 N[nshp] = blend*mode; 47 dN[nshp][0] = bdrv[0]*mode + 48 blend*(mdrv[0]*epdrv[0][0]+mdrv[1]*epdrv[0][1]); 49 dN[nshp][1] = bdrv[1]*mode + 50 blend*(mdrv[0]*epdrv[1][0]+mdrv[1]*epdrv[1][1]); 51 dN[nshp][2] = bdrv[2]*mode + 52 blend*(mdrv[0]*epdrv[2][0]+mdrv[1]*epdrv[2][1]); 53 nshp++; 54 } 55 } 56 } 57 if( p > 2 ) { 58 /* collect all face modes for all p */ 59 for(i=0; i <4; i++) { 60 Lf[0]=L[TetFMAP[i][0]]; 61 Lf[1]=L[TetFMAP[i][1]]; 62 Lf[2]=L[TetFMAP[i][2]]; 63 blend=Lf[0]*Lf[1]*Lf[2]; 64 R_faceBlendTetDrv(TetFMAP[i],L,bdrv); 65 F_parDrv(TetFMAP[i][0],TetFMAP[i][1],TetFMAP[i][2],Stet,&fpdrv); 66 for(ip=3; ip <= p; ip++) { 67 NS = ip-2; 68 for(is=0; is < NS; is++) { 69 mode = F_modeShapeTri(ip,is,Lf); 70 N[nshp] = blend*mode; 71 F_modeShapeTriDrv(ip,is,Lf,mfdrv); 72 mdrv[0]=mfdrv[0]*(double)fpdrv[0][0]+ mfdrv[1]*(double)fpdrv[0][1]; 73 mdrv[1]=mfdrv[0]*(double)fpdrv[1][0]+ mfdrv[1]*(double)fpdrv[1][1]; 74 mdrv[2]=mfdrv[0]*(double)fpdrv[2][0]+ mfdrv[1]*(double)fpdrv[2][1]; 75 dN[nshp][0] = bdrv[0]*mode + blend*mdrv[0] ; 76 dN[nshp][1] = bdrv[1]*mode + blend*mdrv[1] ; 77 dN[nshp][2] = bdrv[2]*mode + blend*mdrv[2] ; 78 nshp++; 79 } 80 } 81 } 82 } 83 if( p > 3 ) { 84 /* collect all region modes */ 85 for(ip=4; ip <= p; ip++) { 86 NS = (ip-3)*(ip-2)/2 ; 87 blend = L[0]*L[1]*L[2]*L[3]; 88 tmp = 1.0e0-L[0]-L[1]-L[2] ; 89 rsw = L[0]*L[1]*tmp; 90 stw = L[1]*L[2]*tmp ; 91 rtw = L[0]*L[2]*tmp ; 92 rst = L[0]*L[1]*L[2] ; 93 rstw = rst*L[3] ; 94 for(is=0; is < NS; is++) { 95 mode = R_modeShapeTet(ip,is,L); 96 N[nshp] = blend*mode; 97 R_modeShapeTetDrv(ip, is, L, mdrv); 98 dN[nshp][0] = mode*(stw-rst) + rstw*mdrv[0] ; 99 dN[nshp][1] = mode*(rtw-rst) + rstw*mdrv[1] ; 100 dN[nshp][2] = mode*(rsw-rst) + rstw*mdrv[2] ; 101 nshp++; 102 } 103 } 104 } 105 return nshp; 106 } 107 108 /* calculate the shape functions and their derivatives for 109 triangular faces 110 */ 111 112 int TriShapeAndDrv(int p,double par[2],double N[],double dN[][2]){ 113 int i,j,nshp=0; 114 double L[3]; 115 116 if(p > 2) /* not supported */ 117 return nshp; 118 L[0]=par[0]; 119 L[1]=par[1]; 120 L[2]=1.0e0-par[0]-par[1]; 121 122 /* define shape functions for a quadratic triangle */ 123 124 /* collect all vertex modes */ 125 for(i=0; i<3; i++) { 126 N[i] = L[i]; 127 if(i==2) 128 dN[i][0]=dN[i][1]=-1.0e0; 129 else { 130 for(j=0; j<2; j++) { 131 if(i==j) 132 dN[i][j] = 1.0e0; 133 else 134 dN[i][j] = 0.0e0; 135 } 136 } 137 } 138 nshp=3; 139 if( p > 1 ){ 140 /* collect edge modes (only quadratic for now) */ 141 N[3] = -2.0e0*L[0]*L[1]; 142 N[4] = -2.0e0*L[1]*L[2]; 143 N[5] = -2.0e0*L[0]*L[2]; 144 dN[3][0] = -2.0e0*L[1]; 145 dN[3][1] = -2.0e0*L[0]; 146 dN[4][0] = 2.0e0*L[1]; 147 dN[4][1] = -2.0e0+2.0e0*L[0]+4.0e0*L[1]; 148 dN[5][0] = -2.0e0+4.0e0*L[0]+2.0e0*L[1]; 149 dN[5][1] = 2.0e0*L[0]; 150 nshp=6; 151 } 152 return nshp; 153 } 154 155 156 157 158 159 160 #ifdef __cplusplus 161 } 162 #endif 163 164