1 #include "shapeFuncInternals.h"
2
3 #ifdef __cplusplus
4 extern "C" {
5 #endif
6
7 /* Calculate shape functions and derivative for tet elements */
TetShapeAndDrv(int p,double par[3],double N[],double dN[][3])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
TriShapeAndDrv(int p,double par[2],double N[],double dN[][2])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