1 #include <stdio.h> 2 3 #ifdef __cplusplus 4 extern "C" { 5 #endif 6 7 void TetQuadLag(double *par, double *N,double dN[][3]); 8 9 int TetShapeAndDrvL(int p,double par[3],double N[],double dN[][3]) { 10 /* 11 static int TetEMAP[6][2]={{0,1},{1,2},{2,0},{0,3},{1,3},{2,3}}; 12 static int TetFMAP[4][3]={{0,1,2},{0,3,1},{1,3,2},{0,2,3}}; 13 */ 14 static int Nshp[3]={4,10,20}; 15 if( p == 2) { 16 TetQuadLag(par,N,dN); 17 } else { 18 fprintf(stderr,"p != 2 Not implemented for lagrange\n"); 19 } 20 return Nshp[p]; 21 } 22 23 void TetQuadLag(double *par, double *N,double dN[][3]) { 24 double r=par[0],s=par[1],t=par[2],t8,t4; 25 t8 = 1.0-r-s-t; 26 N[0] = (2.0*r-1.0)*r; 27 N[1] = (2.0*s-1.0)*s; 28 N[2] = (2.0*t-1.0)*t; 29 N[3] = (1.0-2.0*r-2.0*s-2.0*t)*t8; 30 N[4] = 4.0*r*s; 31 N[5] = 4.0*s*t; 32 N[6] = 4.0*t*r; 33 N[7] = 4.0*r*t8; 34 N[8] = 4.0*s*t8; 35 N[9] = 4.0*t*t8; 36 37 t4 = -3.0+4.0*r+4.0*s+4.0*t; 38 dN[0][0] = 4.0*r-1.0; 39 dN[0][1] = 0.0; 40 dN[0][2] = 0.0; 41 dN[1][0] = 0.0; 42 dN[1][1] = 4.0*s-1.0; 43 dN[1][2] = 0.0; 44 dN[2][0] = 0.0; 45 dN[2][1] = 0.0; 46 dN[2][2] = 4.0*t-1.0; 47 dN[3][0] = t4; 48 dN[3][1] = t4; 49 dN[3][2] = t4; 50 dN[4][0] = 4.0*s; 51 dN[4][1] = 4.0*r; 52 dN[4][2] = 0.0; 53 dN[5][0] = 0.0; 54 dN[5][1] = 4.0*t; 55 dN[5][2] = 4.0*s; 56 dN[6][0] = 4.0*t; 57 dN[6][1] = 0.0; 58 dN[6][2] = 4.0*r; 59 dN[7][0] = 4.0-8.0*r-4.0*s-4.0*t; 60 dN[7][1] = -4.0*r; 61 dN[7][2] = -4.0*r; 62 dN[8][0] = -4.0*s; 63 dN[8][1] = 4.0-4.0*r-8.0*s-4.0*t; 64 dN[8][2] = -4.0*s; 65 dN[9][0] = -4.0*t; 66 dN[9][1] = -4.0*t; 67 dN[9][2] = 4.0-4.0*r-4.0*s-8.0*t; 68 } 69 70 #ifdef __cplusplus 71 } 72 #endif 73 74