1 /*-------------------------------------------------------------------------
2 Scientific Computation Research Center, RPI, Troy NY
3 (C) Copyright 1995, RPI-SCOREC
4
5 Project : shapeFuntions
6 Author(s): Saikat Dey
7 Creation : Oct., 95
8 Modifi. :
9 Function :
10 evaluate derivatives of entity blend functions over other
11 entities.
12 -------------------------------------------------------------------------*/
13
14 #include <math.h>
15 #include "shapeFuncInternals.h"
16
17 #ifdef __cplusplus
18 extern "C" {
19 #endif
20
V_blendOnEntityDrv(int index,int etype,double * L,double mdrv[3])21 int V_blendOnEntityDrv(int index, int etype, double *L, double mdrv[3]) {
22 /* derivative of blend a vertex mode on a mesh entity */
23 if( etype == Sedge ) {
24 if(index == 0) /* N = 0.5(1-zi) */
25 mdrv[0] = -0.5;
26 else /* N = 0.5(1+zi) */
27 mdrv[0] = +0.5;
28 return 1 ;
29 } else if( etype == Stri ) {
30 mdrv[0] = mdrv[1] = 0.0;
31 if( index == 2 )
32 mdrv[0] = mdrv[1] = -1.0;
33 else
34 mdrv[index] = 1.0;
35 return 1 ;
36 } else if ( etype == Stet ) {
37 mdrv[0] = mdrv[1] = mdrv[2] = 0.0;
38 if( index == 3 )
39 mdrv[0] = mdrv[1] = mdrv[2] = -1.0;
40 else
41 mdrv[index] = 1.0;
42 return 1;
43 } else
44 return 0;
45 }
46
E_blendOnFaceDrv(int index[],int etype,double * L,double bdrv[2])47 int E_blendOnFaceDrv(int index[], int etype, double *L, double bdrv[2]) {
48 if( etype == Stri )
49 return F_edgeBlendTriDrv(index, L, bdrv);
50 else
51 return 0;
52 }
53
F_edgeBlendTriDrv(int index[2],double * L,double drv[])54 int F_edgeBlendTriDrv(int index[2], double *L, double drv[]) {
55 double r,s;
56
57 drv[0] = drv[1] = 0.0;
58 r = L[0] ; s = L[1] ;
59
60 /* figure out which edge we are dealing with */
61 /* v0=0, V1=1 */
62 if( (index[0]==0 && index[1]==1) || (index[0]==1 && index[1]==0) ) {
63 drv[0] = -2.0*s;
64 drv[1] = -2.0*r;
65 /* v0=1, V1=2 */
66 } else if( (index[0]==1 && index[1]==2) || (index[0]==2 && index[1]==1) ) {
67 drv[0] = 2.0*s;
68 drv[1] = -2.0+2.0*r+4.0*s;
69 /* v0=2, V1=0 */
70 } else if( (index[0]==2 && index[1]==0) || (index[0]==0 && index[1]==2) ) {
71 drv[0] = 4.0*r-2.0+2.0*s;
72 drv[1] = 2.0*r;
73 } else
74 return 0;
75
76 return 1 ;
77 }
78
E_blendOnRegionDrv(int index[],int etype,double * L,double bdrv[3])79 int E_blendOnRegionDrv(int index[], int etype, double *L, double bdrv[3]) {
80 if( etype == 4 )
81 return R_edgeBlendTetDrv(index, L, bdrv);
82 else
83 return 0;
84 }
85
F_blendOnRegionDrv(int index[],int etype,double * L,double drv[3])86 int F_blendOnRegionDrv(int index[], int etype, double *L, double drv[3]) {
87 /* blend a face mode on a region */
88 drv[0] = drv[1] = drv[2] = 0.0;
89 if(etype==Stet )
90 return 1;
91 else
92 return 0;
93 }
94
R_edgeBlendTetDrv(int * index,double * L,double drv[])95 int R_edgeBlendTetDrv(int *index, double *L, double drv[]) {
96 double r=L[0],s=L[1],t=L[2];
97
98 /* the blend is given by -2*L[i]*L[j] */
99 /* v0=0, V1=1 */
100 if( (index[0]==0 && index[1]==1) || (index[0]==1 && index[1]==0) ) {
101 drv[0] = -2.0*s;
102 drv[1] = -2.0*r;
103 drv[2] = 0.0;
104 /* v0=0, V1=3 */
105 } else if( (index[0]==0 && index[1]==3) || (index[0]==3 && index[1]==0) ) {
106 drv[0] = -2.0+4.0*r+2.0*s+2.0*t;
107 drv[1] = 2.0*r;
108 drv[2] = 2.0*r;
109 /* v0=1, V1=2 */
110 } else if( (index[0]==1 && index[1]==2) || (index[0]==2 && index[1]==1) ) {
111 drv[0] = 0.0;
112 drv[1] = -2.0*t;
113 drv[2] = -2.0*s;
114 /* v0=1, V1=3 */
115 } else if( (index[0]==1 && index[1]==3) || (index[0]==3 && index[1]==1) ) {
116 drv[0] = 2.0*s;
117 drv[1] = -2.0+2.0*r+4.0*s+2.0*t;
118 drv[2] = 2.0*s;
119 /* v0=2, V1=0 */
120 } else if( (index[0]==2 && index[1]==0) || (index[0]==0 && index[1]==2) ) {
121 drv[0] = -2.0*t;
122 drv[1] = 0.0;
123 drv[2] = -2.0*r;
124 /* v0=2, V1=3 */
125 } else if( (index[0]==2 && index[1]==3) || (index[0]==3 && index[1]==2) ) {
126 drv[0] = 2.0*t;
127 drv[1] = 2.0*t;
128 drv[2] = -2.0+2.0*r+2.0*s+4.0*t;
129 } else
130 return 0;
131
132 return 1;
133 }
134
R_faceBlendTetDrv(int * index,double * L,double drv[])135 int R_faceBlendTetDrv(int *index, double *L, double drv[]) {
136 int isum = index[0]+index[1]+index[2];
137
138 drv[0]=drv[1]=drv[2]=0;
139 if( isum == 3) /* r*s*t */ {
140 drv[0] = L[1]*L[2] ;
141 drv[1] = L[0]*L[2] ;
142 drv[2] = L[1]*L[0] ;
143 } else if( isum == 5) /* r*t*(1-r-s-t) */ {
144 drv[0] = L[2]*(1.0-2.0*L[0]-L[1]-L[2]);
145 drv[1] = -L[0]*L[2] ;
146 drv[2] = L[0]*(1.0-L[0]-L[1]-2.0*L[2]);
147 } else if( isum == 4) /* r*s*(1-r-s-t) */ {
148 drv[0] = L[1]*(1.0-2.0*L[0]-L[1]-L[2]);
149 drv[1] = L[0]*(1.0-L[0]-2.0*L[1]-L[2]);
150 drv[2] = -L[0]*L[1] ;
151 } else if( isum == 6) /* s*t*(1-r-s-t) */ {
152 drv[0] = -L[1]*L[2] ;
153 drv[1] = L[2]*(1.0-L[0]-2.0*L[1]-L[2]);
154 drv[2] = L[1]*(1.0-L[0]-L[1]-2.0*L[2]);
155 } else
156 return 0;
157 return 1;
158 }
159 #ifdef __cplusplus
160 }
161 #endif
162