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 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 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 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 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 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 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 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