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
TetShapeAndDrvL(int p,double par[3],double N[],double dN[][3])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
TetQuadLag(double * par,double * N,double dN[][3])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