Changeset 5745
- Timestamp:
- 09/10/10 10:53:39 (15 years ago)
- Location:
- issm/trunk/src/c/objects/Elements
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Element.h
r5743 r5745 32 32 virtual void GetSolutionFromInputs(Vec solution)=0; 33 33 virtual int GetNodeIndex(Node* node)=0; 34 virtual bool GetShelf()=0;35 virtual bool GetOnBed()=0;34 virtual bool IsOnShelf()=0; 35 virtual bool IsOnBed()=0; 36 36 virtual void GetParameterListOnVertices(double* pvalue,int enumtype)=0; 37 37 virtual void GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue)=0; -
issm/trunk/src/c/objects/Elements/Penta.cpp
r5743 r5745 897 897 ISSMERROR("Node provided not found among element nodes"); 898 898 899 }900 /*}}}*/901 /*FUNCTION Penta::GetOnBed {{{1*/902 bool Penta::GetOnBed(){903 904 bool onbed;905 906 inputs->GetParameterValue(&onbed,ElementOnBedEnum);907 908 return onbed;909 }910 /*}}}*/911 /*FUNCTION Penta::GetShelf {{{1*/912 bool Penta::GetShelf(){913 bool onshelf;914 915 /*retrieve inputs :*/916 inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum);917 918 return onshelf;919 899 } 920 900 /*}}}*/ … … 2001 1981 2002 1982 /*If shelf: hydrostatic equilibrium*/ 2003 if (this-> GetShelf()){1983 if (this->IsOnShelf()){ 2004 1984 2005 1985 /*recover material parameters: */ … … 5483 5463 } 5484 5464 /*}}}*/ 5465 /*FUNCTION Penta::IsOnShelf {{{1*/ 5466 bool Penta::IsOnShelf(){ 5467 5468 bool onshelf; 5469 inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum); 5470 return onshelf; 5471 5472 } 5473 /*}}}*/ 5485 5474 /*FUNCTION Penta::IsOnSurface{{{1*/ 5486 5475 bool Penta::IsOnSurface(void){ -
issm/trunk/src/c/objects/Elements/Penta.h
r5743 r5745 79 79 void DeleteResults(void); 80 80 int GetNodeIndex(Node* node); 81 bool GetOnBed();82 bool GetShelf();83 81 void GetSolutionFromInputs(Vec solution); 84 82 void GetVectorFromInputs(Vec vector,int NameEnum); … … 179 177 bool IsOnSurface(void); 180 178 bool IsOnBed(void); 179 bool IsOnShelf(void); 181 180 void ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp); 182 181 void ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp); -
issm/trunk/src/c/objects/Elements/PentaRef.cpp
r5743 r5745 19 19 #include "../../include/include.h" 20 20 /*}}}*/ 21 22 /*Element macros*/ 23 #define NUMNODESP1 6 24 #define NUMNODESP1_2d 3 25 #define NUMNODESMINI 7 26 #define NDOF1 1 27 #define NDOF2 2 28 #define NDOF3 3 29 #define NDOF4 4 21 30 22 31 /*Object constructors and destructor*/ … … 62 71 * where h is the interpolation function for grid i. 63 72 * 64 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) 65 */ 66 67 int i; 68 const int numgrids=6; 69 const int NDOF3=3; 70 const int NDOF2=2; 71 72 double dh1dh6[NDOF3][numgrids]; 73 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) 74 */ 75 76 double dh1dh6[3][NUMNODESP1]; 73 77 74 78 /*Get dh1dh6 in actual coordinate system: */ … … 76 80 77 81 /*Build B: */ 78 for (i =0;i<numgrids;i++){79 *(B+NDOF2* numgrids*0+NDOF2*i)=dh1dh6[0][i];80 *(B+NDOF2* numgrids*0+NDOF2*i+1)=0.0;81 82 *(B+NDOF2* numgrids*1+NDOF2*i)=0.0;83 *(B+NDOF2* numgrids*1+NDOF2*i+1)=dh1dh6[1][i];84 85 *(B+NDOF2* numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];86 *(B+NDOF2* numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];82 for (int i=0;i<NUMNODESP1;i++){ 83 *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dh1dh6[0][i]; 84 *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0; 85 86 *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0; 87 *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dh1dh6[1][i]; 88 89 *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 90 *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 87 91 88 92 } … … 102 106 * where h is the interpolation function for grid i. 103 107 * 104 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) 105 */ 106 107 int i; 108 const int numgrids=6; 109 const int NDOF3=3; 110 const int NDOF2=2; 111 112 double dh1dh6[NDOF3][numgrids]; 108 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) 109 */ 110 111 double dh1dh6[3][NUMNODESP1]; 113 112 114 113 /*Get dh1dh6 in actual coordinate system: */ … … 116 115 117 116 /*Build B: */ 118 for (i =0;i<numgrids;i++){119 *(B+NDOF2* numgrids*0+NDOF2*i)=dh1dh6[0][i];120 *(B+NDOF2* numgrids*0+NDOF2*i+1)=0.0;121 122 *(B+NDOF2* numgrids*1+NDOF2*i)=0.0;123 *(B+NDOF2* numgrids*1+NDOF2*i+1)=dh1dh6[1][i];124 125 *(B+NDOF2* numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];126 *(B+NDOF2* numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];127 128 *(B+NDOF2* numgrids*3+NDOF2*i)=(float).5*dh1dh6[2][i];129 *(B+NDOF2* numgrids*3+NDOF2*i+1)=0.0;130 131 *(B+NDOF2* numgrids*4+NDOF2*i)=0.0;132 *(B+NDOF2* numgrids*4+NDOF2*i+1)=(float).5*dh1dh6[2][i];117 for (int i=0;i<NUMNODESP1;i++){ 118 *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dh1dh6[0][i]; 119 *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0; 120 121 *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0; 122 *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dh1dh6[1][i]; 123 124 *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 125 *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 126 127 *(B+NDOF2*NUMNODESP1*3+NDOF2*i)=(float).5*dh1dh6[2][i]; 128 *(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0; 129 130 *(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0; 131 *(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=(float).5*dh1dh6[2][i]; 133 132 } 134 133 … … 147 146 * where h is the interpolation function for grid i. 148 147 * 149 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) 150 */ 151 152 int i; 153 const int NDOF3=3; 154 const int NDOF2=2; 155 const int numgrids=6; 156 157 double dh1dh6[NDOF3][numgrids]; 148 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) 149 */ 150 double dh1dh6[3][NUMNODESP1]; 158 151 159 152 /*Get dh1dh6 in actual coordinate system: */ … … 161 154 162 155 /*Build BPrime: */ 163 for (i =0;i<numgrids;i++){164 *(B+NDOF2* numgrids*0+NDOF2*i)=2.0*dh1dh6[0][i];165 *(B+NDOF2* numgrids*0+NDOF2*i+1)=dh1dh6[1][i];166 167 *(B+NDOF2* numgrids*1+NDOF2*i)=dh1dh6[0][i];168 *(B+NDOF2* numgrids*1+NDOF2*i+1)=2.0*dh1dh6[1][i];169 170 *(B+NDOF2* numgrids*2+NDOF2*i)=dh1dh6[1][i];171 *(B+NDOF2* numgrids*2+NDOF2*i+1)=dh1dh6[0][i];172 173 *(B+NDOF2* numgrids*3+NDOF2*i)=dh1dh6[2][i];174 *(B+NDOF2* numgrids*3+NDOF2*i+1)=0.0;175 176 *(B+NDOF2* numgrids*4+NDOF2*i)=0.0;177 *(B+NDOF2* numgrids*4+NDOF2*i+1)=dh1dh6[2][i];156 for (int i=0;i<NUMNODESP1;i++){ 157 *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=2.0*dh1dh6[0][i]; 158 *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dh1dh6[1][i]; 159 160 *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=dh1dh6[0][i]; 161 *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2.0*dh1dh6[1][i]; 162 163 *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=dh1dh6[1][i]; 164 *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dh1dh6[0][i]; 165 166 *(B+NDOF2*NUMNODESP1*3+NDOF2*i)=dh1dh6[2][i]; 167 *(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0; 168 169 *(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0; 170 *(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=dh1dh6[2][i]; 178 171 } 179 172 } … … 182 175 void PentaRef::GetBStokes(double* B, double* xyz_list, GaussPenta* gauss){ 183 176 184 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3* DOFPERGRID.177 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4. 185 178 * For grid i, Bi can be expressed in the actual coordinate system 186 179 * by: Bi=[ dh/dx 0 0 0 ] … … 197 190 198 191 int i; 199 const int calculationdof=3; 200 const int numgrids=6; 201 int DOFPERGRID=4; 202 203 double dh1dh7[calculationdof][numgrids+1]; 204 double l1l6[numgrids]; 192 193 double dh1dh7[3][NUMNODESMINI]; 194 double l1l6[NUMNODESP1]; 205 195 206 196 … … 210 200 211 201 /*Build B: */ 212 for (i=0;i< numgrids+1;i++){213 *(B+( DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B[0][DOFPERGRID*i]=dh1dh6[0][i];214 *(B+( DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;215 *(B+( DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;216 *(B+( DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;217 *(B+( DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];218 *(B+( DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;219 *(B+( DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;220 *(B+( DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;221 *(B+( DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];222 *(B+( DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7[1][i];223 *(B+( DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7[0][i];224 *(B+( DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;225 *(B+( DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7[2][i];226 *(B+( DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;227 *(B+( DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7[0][i];228 *(B+( DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;229 *(B+( DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7[2][i];230 *(B+( DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7[1][i];231 *(B+( DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0;232 *(B+( DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0;233 *(B+( DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0;234 *(B+( DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7[0][i];235 *(B+( DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7[1][i];236 *(B+( DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7[2][i];237 } 238 239 for (i=0;i< numgrids;i++){ //last column not for the bubble function240 *(B+( DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;241 *(B+( DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;242 *(B+( DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;243 *(B+( DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;244 *(B+( DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;245 *(B+( DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;246 *(B+( DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=l1l6[i];247 *(B+( DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=0;202 for (i=0;i<NUMNODESMINI;i++){ 203 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i]; 204 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0; 205 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0; 206 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0; 207 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i]; 208 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0; 209 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0; 210 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0; 211 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh7[2][i]; 212 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=(float).5*dh1dh7[1][i]; 213 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=(float).5*dh1dh7[0][i]; 214 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0; 215 *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=(float).5*dh1dh7[2][i]; 216 *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0; 217 *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=(float).5*dh1dh7[0][i]; 218 *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0; 219 *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=(float).5*dh1dh7[2][i]; 220 *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=(float).5*dh1dh7[1][i]; 221 *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=0; 222 *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=0; 223 *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=0; 224 *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=dh1dh7[0][i]; 225 *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=dh1dh7[1][i]; 226 *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=dh1dh7[2][i]; 227 } 228 229 for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function 230 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0; 231 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0; 232 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0; 233 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0; 234 *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0; 235 *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0; 236 *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=l1l6[i]; 237 *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=0; 248 238 } 249 239 … … 269 259 270 260 int i; 271 const int calculationdof=3; 272 const int numgrids=6; 273 int DOFPERGRID=4; 274 275 double dh1dh7[calculationdof][numgrids+1]; 276 double l1l6[numgrids]; 261 double dh1dh7[3][NUMNODESMINI]; 262 double l1l6[NUMNODESP1]; 277 263 278 264 /*Get dh1dh7 in actual coordinate system: */ 279 265 GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss); 280 281 266 GetNodalFunctionsP1(l1l6, gauss); 282 267 283 268 /*B_primeuild B_prime: */ 284 for (i=0;i< numgrids+1;i++){285 *(B_prime+( DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B_prime[0][DOFPERGRID*i]=dh1dh6[0][i];286 *(B_prime+( DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;287 *(B_prime+( DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;288 *(B_prime+( DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;289 *(B_prime+( DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];290 *(B_prime+( DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;291 *(B_prime+( DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;292 *(B_prime+( DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;293 *(B_prime+( DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];294 *(B_prime+( DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=dh1dh7[1][i];295 *(B_prime+( DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=dh1dh7[0][i];296 *(B_prime+( DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;297 *(B_prime+( DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=dh1dh7[2][i];298 *(B_prime+( DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;299 *(B_prime+( DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=dh1dh7[0][i];300 *(B_prime+( DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;301 *(B_prime+( DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=dh1dh7[2][i];302 *(B_prime+( DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=dh1dh7[1][i];303 *(B_prime+( DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=dh1dh7[0][i];304 *(B_prime+( DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=dh1dh7[1][i];305 *(B_prime+( DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=dh1dh7[2][i];306 *(B_prime+( DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=0;307 *(B_prime+( DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=0;308 *(B_prime+( DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=0;309 } 310 311 for (i=0;i< numgrids;i++){ //last column not for the bubble function312 *(B_prime+( DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;313 *(B_prime+( DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;314 *(B_prime+( DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;315 *(B_prime+( DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;316 *(B_prime+( DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;317 *(B_prime+( DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;318 *(B_prime+( DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=0;319 *(B_prime+( DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=l1l6[i];269 for (i=0;i<NUMNODESMINI;i++){ 270 *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B_prime[0][NDOF4*i]=dh1dh6[0][i]; 271 *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0; 272 *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0; 273 *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0; 274 *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i]; 275 *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0; 276 *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0; 277 *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0; 278 *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh7[2][i]; 279 *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=dh1dh7[1][i]; 280 *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=dh1dh7[0][i]; 281 *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0; 282 *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=dh1dh7[2][i]; 283 *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0; 284 *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=dh1dh7[0][i]; 285 *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0; 286 *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=dh1dh7[2][i]; 287 *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=dh1dh7[1][i]; 288 *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=dh1dh7[0][i]; 289 *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=dh1dh7[1][i]; 290 *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=dh1dh7[2][i]; 291 *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=0; 292 *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=0; 293 *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=0; 294 } 295 296 for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function 297 *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0; 298 *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0; 299 *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0; 300 *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0; 301 *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0; 302 *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0; 303 *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=0; 304 *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=l1l6[i]; 320 305 } 321 306 … … 324 309 /*FUNCTION PentaRef::GetBArtdiff {{{1*/ 325 310 void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss){ 326 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5* DOFPERGRID.311 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 327 312 * For grid i, Bi' can be expressed in the actual coordinate system 328 313 * by: … … 331 316 * where h is the interpolation function for grid i. 332 317 * 333 * We assume B has been allocated already, of size: 2x(DOFPERGRID*numgrids) 334 */ 335 336 int i; 337 const int calculationdof=3; 338 const int numgrids=6; 339 int DOFPERGRID=1; 318 * We assume B has been allocated already, of size: 2x(NDOF1*NUMNODESP1) 319 */ 340 320 341 321 /*Same thing in the actual coordinate system: */ 342 double dh1dh6[ calculationdof][numgrids];322 double dh1dh6[3][NUMNODESP1]; 343 323 344 324 /*Get dh1dh6 in actual coordinates system : */ … … 346 326 347 327 /*Build B': */ 348 for (i =0;i<numgrids;i++){349 *(B_artdiff+ DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];350 *(B_artdiff+ DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];328 for (int i=0;i<NUMNODESP1;i++){ 329 *(B_artdiff+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i]; 330 *(B_artdiff+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i]; 351 331 } 352 332 } … … 354 334 /*FUNCTION PentaRef::GetBAdvec{{{1*/ 355 335 void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss){ 356 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5* DOFPERGRID.336 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 357 337 * For grid i, Bi' can be expressed in the actual coordinate system 358 338 * by: … … 362 342 * where h is the interpolation function for grid i. 363 343 * 364 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids) 365 */ 366 367 int i; 368 int calculationdof=3; 369 int numgrids=6; 370 int DOFPERGRID=1; 344 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1) 345 */ 371 346 372 347 /*Same thing in the actual coordinate system: */ … … 377 352 378 353 /*Build B': */ 379 for (i =0;i<numgrids;i++){380 *(B_advec+ DOFPERGRID*numgrids*0+DOFPERGRID*i)=l1l6[i];381 *(B_advec+ DOFPERGRID*numgrids*1+DOFPERGRID*i)=l1l6[i];382 *(B_advec+ DOFPERGRID*numgrids*2+DOFPERGRID*i)=l1l6[i];354 for (int i=0;i<NUMNODESP1;i++){ 355 *(B_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=l1l6[i]; 356 *(B_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=l1l6[i]; 357 *(B_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=l1l6[i]; 383 358 } 384 359 } … … 386 361 /*FUNCTION PentaRef::GetBConduct{{{1*/ 387 362 void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss){ 388 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5* DOFPERGRID.363 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 389 364 * For grid i, Bi' can be expressed in the actual coordinate system 390 365 * by: … … 394 369 * where h is the interpolation function for grid i. 395 370 * 396 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids) 397 */ 398 399 int i; 400 const int calculationdof=3; 401 const int numgrids=6; 402 int DOFPERGRID=1; 371 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1) 372 */ 403 373 404 374 /*Same thing in the actual coordinate system: */ 405 double dh1dh6[ calculationdof][numgrids];375 double dh1dh6[3][NUMNODESP1]; 406 376 407 377 /*Get dh1dh2dh3 in actual coordinates system : */ … … 409 379 410 380 /*Build B': */ 411 for (i =0;i<numgrids;i++){412 *(B_conduct+ DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];413 *(B_conduct+ DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];414 *(B_conduct+ DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i];381 for (int i=0;i<NUMNODESP1;i++){ 382 *(B_conduct+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i]; 383 *(B_conduct+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i]; 384 *(B_conduct+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i]; 415 385 } 416 386 } … … 422 392 423 393 int i; 424 const int NDOF3=3; 425 const int numgrids=6; 426 double dh1dh6[NDOF3][numgrids]; 394 double dh1dh6[3][NUMNODESP1]; 427 395 428 396 /*Get dh1dh6 in actual coordinate system: */ … … 430 398 431 399 /*Build B: */ 432 for (i=0;i< numgrids;i++){400 for (i=0;i<NUMNODESP1;i++){ 433 401 B[i]=dh1dh6[2][i]; 434 402 } … … 438 406 /*FUNCTION PentaRef::GetBprimeAdvec{{{1*/ 439 407 void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss){ 440 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5* DOFPERGRID.408 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 441 409 * For grid i, Bi' can be expressed in the actual coordinate system 442 410 * by: … … 446 414 * where h is the interpolation function for grid i. 447 415 * 448 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids) 449 */ 450 451 int i; 452 const int calculationdof=3; 453 const int numgrids=6; 454 int DOFPERGRID=1; 416 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1) 417 */ 455 418 456 419 /*Same thing in the actual coordinate system: */ 457 double dh1dh6[ calculationdof][numgrids];420 double dh1dh6[3][NUMNODESP1]; 458 421 459 422 /*Get dh1dh2dh3 in actual coordinates system : */ … … 461 424 462 425 /*Build B': */ 463 for (i =0;i<numgrids;i++){464 *(Bprime_advec+ DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];465 *(Bprime_advec+ DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];466 *(Bprime_advec+ DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i];426 for (int i=0;i<NUMNODESP1;i++){ 427 *(Bprime_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i]; 428 *(Bprime_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i]; 429 *(Bprime_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i]; 467 430 } 468 431 } … … 500 463 501 464 int i; 502 const int numgrids2d=3;503 465 int num_dof=4; 504 466 505 double l1l2l3[ numgrids2d];467 double l1l2l3[NUMNODESP1_2d]; 506 468 507 469 … … 513 475 /*Build LStokes: */ 514 476 for (i=0;i<3;i++){ 515 *(LStokes+num_dof* numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];516 *(LStokes+num_dof* numgrids2d*0+num_dof*i+1)=0;517 *(LStokes+num_dof* numgrids2d*0+num_dof*i+2)=0;518 *(LStokes+num_dof* numgrids2d*0+num_dof*i+3)=0;519 *(LStokes+num_dof* numgrids2d*1+num_dof*i)=0;520 *(LStokes+num_dof* numgrids2d*1+num_dof*i+1)=l1l2l3[i];521 *(LStokes+num_dof* numgrids2d*1+num_dof*i+2)=0;522 *(LStokes+num_dof* numgrids2d*1+num_dof*i+3)=0;523 *(LStokes+num_dof* numgrids2d*2+num_dof*i)=0;524 *(LStokes+num_dof* numgrids2d*2+num_dof*i+1)=0;525 *(LStokes+num_dof* numgrids2d*2+num_dof*i+2)=l1l2l3[i];526 *(LStokes+num_dof* numgrids2d*2+num_dof*i+3)=0;527 *(LStokes+num_dof* numgrids2d*3+num_dof*i)=0;528 *(LStokes+num_dof* numgrids2d*3+num_dof*i+1)=0;529 *(LStokes+num_dof* numgrids2d*3+num_dof*i+2)=l1l2l3[i];530 *(LStokes+num_dof* numgrids2d*3+num_dof*i+3)=0;531 *(LStokes+num_dof* numgrids2d*4+num_dof*i)=l1l2l3[i];532 *(LStokes+num_dof* numgrids2d*4+num_dof*i+1)=0;533 *(LStokes+num_dof* numgrids2d*4+num_dof*i+2)=0;534 *(LStokes+num_dof* numgrids2d*4+num_dof*i+3)=0;535 *(LStokes+num_dof* numgrids2d*5+num_dof*i)=0;536 *(LStokes+num_dof* numgrids2d*5+num_dof*i+1)=l1l2l3[i];537 *(LStokes+num_dof* numgrids2d*5+num_dof*i+2)=0;538 *(LStokes+num_dof* numgrids2d*5+num_dof*i+3)=0;539 *(LStokes+num_dof* numgrids2d*6+num_dof*i)=l1l2l3[i];540 *(LStokes+num_dof* numgrids2d*6+num_dof*i+1)=0;541 *(LStokes+num_dof* numgrids2d*6+num_dof*i+2)=0;542 *(LStokes+num_dof* numgrids2d*6+num_dof*i+3)=0;543 *(LStokes+num_dof* numgrids2d*7+num_dof*i)=0;544 *(LStokes+num_dof* numgrids2d*7+num_dof*i+1)=l1l2l3[i];545 *(LStokes+num_dof* numgrids2d*7+num_dof*i+2)=0;546 *(LStokes+num_dof* numgrids2d*7+num_dof*i+3)=0;547 *(LStokes+num_dof* numgrids2d*8+num_dof*i)=0;548 *(LStokes+num_dof* numgrids2d*8+num_dof*i+1)=0;549 *(LStokes+num_dof* numgrids2d*8+num_dof*i+2)=l1l2l3[i];550 *(LStokes+num_dof* numgrids2d*8+num_dof*i+3)=0;551 *(LStokes+num_dof* numgrids2d*9+num_dof*i)=0;552 *(LStokes+num_dof* numgrids2d*9+num_dof*i+1)=0;553 *(LStokes+num_dof* numgrids2d*9+num_dof*i+2)=l1l2l3[i];554 *(LStokes+num_dof* numgrids2d*9+num_dof*i+3)=0;555 *(LStokes+num_dof* numgrids2d*10+num_dof*i)=0;556 *(LStokes+num_dof* numgrids2d*10+num_dof*i+1)=0;557 *(LStokes+num_dof* numgrids2d*10+num_dof*i+2)=l1l2l3[i];558 *(LStokes+num_dof* numgrids2d*10+num_dof*i+3)=0;559 *(LStokes+num_dof* numgrids2d*11+num_dof*i)=l1l2l3[i];560 *(LStokes+num_dof* numgrids2d*11+num_dof*i+1)=0;561 *(LStokes+num_dof* numgrids2d*11+num_dof*i+2)=0;562 *(LStokes+num_dof* numgrids2d*11+num_dof*i+3)=0;563 *(LStokes+num_dof* numgrids2d*12+num_dof*i)=0;564 *(LStokes+num_dof* numgrids2d*12+num_dof*i+1)=l1l2l3[i];565 *(LStokes+num_dof* numgrids2d*12+num_dof*i+2)=0;566 *(LStokes+num_dof* numgrids2d*12+num_dof*i+3)=0;567 *(LStokes+num_dof* numgrids2d*13+num_dof*i)=0;568 *(LStokes+num_dof* numgrids2d*13+num_dof*i+1)=0;569 *(LStokes+num_dof* numgrids2d*13+num_dof*i+2)=l1l2l3[i];570 *(LStokes+num_dof* numgrids2d*13+num_dof*i+3)=0;477 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i]; 478 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0; 479 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0; 480 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0; 481 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0; 482 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i]; 483 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0; 484 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0; 485 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0; 486 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0; 487 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i]; 488 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0; 489 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0; 490 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0; 491 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i]; 492 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0; 493 *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=l1l2l3[i]; 494 *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0; 495 *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=0; 496 *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0; 497 *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0; 498 *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=l1l2l3[i]; 499 *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=0; 500 *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0; 501 *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=l1l2l3[i]; 502 *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0; 503 *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=0; 504 *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=0; 505 *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0; 506 *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=l1l2l3[i]; 507 *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=0; 508 *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=0; 509 *(LStokes+num_dof*NUMNODESP1_2d*8+num_dof*i)=0; 510 *(LStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+1)=0; 511 *(LStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+2)=l1l2l3[i]; 512 *(LStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+3)=0; 513 *(LStokes+num_dof*NUMNODESP1_2d*9+num_dof*i)=0; 514 *(LStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+1)=0; 515 *(LStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+2)=l1l2l3[i]; 516 *(LStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+3)=0; 517 *(LStokes+num_dof*NUMNODESP1_2d*10+num_dof*i)=0; 518 *(LStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+1)=0; 519 *(LStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+2)=l1l2l3[i]; 520 *(LStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+3)=0; 521 *(LStokes+num_dof*NUMNODESP1_2d*11+num_dof*i)=l1l2l3[i]; 522 *(LStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+1)=0; 523 *(LStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+2)=0; 524 *(LStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+3)=0; 525 *(LStokes+num_dof*NUMNODESP1_2d*12+num_dof*i)=0; 526 *(LStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+1)=l1l2l3[i]; 527 *(LStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+2)=0; 528 *(LStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+3)=0; 529 *(LStokes+num_dof*NUMNODESP1_2d*13+num_dof*i)=0; 530 *(LStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+1)=0; 531 *(LStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+2)=l1l2l3[i]; 532 *(LStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+3)=0; 571 533 572 534 } … … 597 559 */ 598 560 int i; 599 const int numgrids2d=3;600 561 int num_dof=4; 601 562 602 double l1l2l3[ numgrids2d];603 double dh1dh6[3][ 6];563 double l1l2l3[NUMNODESP1_2d]; 564 double dh1dh6[3][NUMNODESP1]; 604 565 605 566 /*Get l1l2l3 in actual coordinate system: */ … … 612 573 /*Build LprimeStokes: */ 613 574 for (i=0;i<3;i++){ 614 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i]; 615 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+1)=0; 616 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+2)=0; 617 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+3)=0; 618 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i)=0; 619 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i]; 620 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+2)=0; 621 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+3)=0; 622 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i)=l1l2l3[i]; 623 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+1)=0; 624 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+2)=0; 625 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+3)=0; 626 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i)=0; 627 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+1)=l1l2l3[i]; 628 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+2)=0; 629 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+3)=0; 630 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i)=0; 631 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+1)=0; 632 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+2)=l1l2l3[i]; 633 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+3)=0; 634 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i)=0; 635 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+1)=0; 636 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+2)=l1l2l3[i]; 637 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+3)=0; 638 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i)=0; 639 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+1)=0; 640 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6[2][i]; 641 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+3)=0; 642 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i)=0; 643 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+1)=0; 644 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6[2][i]; 645 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+3)=0; 646 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i)=0; 647 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+1)=0; 648 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6[2][i]; 649 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+3)=0; 650 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6[2][i]; 651 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+1)=0; 652 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6[0][i]; 653 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+3)=0; 654 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i)=0; 655 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6[2][i]; 656 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6[1][i]; 657 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+3)=0; 658 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i)=0; 659 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+1)=0; 660 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+2)=0; 661 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+3)=l1l2l3[i]; 662 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i)=0; 663 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+1)=0; 664 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+2)=0; 665 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+3)=l1l2l3[i]; 666 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i)=0; 667 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+1)=0; 668 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0; 669 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i]; 670 575 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i]; 576 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0; 577 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0; 578 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0; 579 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0; 580 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i]; 581 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0; 582 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0; 583 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i]; 584 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0; 585 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=0; 586 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0; 587 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0; 588 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i]; 589 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=0; 590 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0; 591 *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0; 592 *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0; 593 *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=l1l2l3[i]; 594 *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0; 595 *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0; 596 *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0; 597 *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=l1l2l3[i]; 598 *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0; 599 *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0; 600 *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0; 601 *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=dh1dh6[2][i]; 602 *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=0; 603 *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0; 604 *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0; 605 *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=dh1dh6[2][i]; 606 *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=0; 607 *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i)=0; 608 *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+1)=0; 609 *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+2)=dh1dh6[2][i]; 610 *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+3)=0; 611 *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i)=dh1dh6[2][i]; 612 *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+1)=0; 613 *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+2)=dh1dh6[0][i]; 614 *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+3)=0; 615 *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i)=0; 616 *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+1)=dh1dh6[2][i]; 617 *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+2)=dh1dh6[1][i]; 618 *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+3)=0; 619 *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i)=0; 620 *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+1)=0; 621 *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+2)=0; 622 *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+3)=l1l2l3[i]; 623 *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i)=0; 624 *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+1)=0; 625 *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+2)=0; 626 *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+3)=l1l2l3[i]; 627 *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i)=0; 628 *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+1)=0; 629 *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+2)=0; 630 *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+3)=l1l2l3[i]; 671 631 } 672 632 } … … 675 635 void PentaRef::GetJacobian(double* J, double* xyz_list,GaussPenta* gauss){ 676 636 677 const int NDOF3=3;678 637 int i,j; 679 638 … … 822 781 823 782 int i; 824 const int numgrids = 7; 825 double dh1dh7_ref[3][numgrids]; 783 double dh1dh7_ref[3][NUMNODESMINI]; 826 784 double Jinv[3][3]; 827 785 … … 839 797 */ 840 798 841 for (i=0;i< numgrids;i++){842 *(dh1dh7+ numgrids*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];843 *(dh1dh7+ numgrids*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];844 *(dh1dh7+ numgrids*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];799 for (i=0;i<NUMNODESMINI;i++){ 800 *(dh1dh7+NUMNODESMINI*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i]; 801 *(dh1dh7+NUMNODESMINI*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i]; 802 *(dh1dh7+NUMNODESMINI*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i]; 845 803 } 846 804 … … 852 810 /*This routine returns the values of the nodal functions derivatives (with respect to the 853 811 * natural coordinate system) at the gaussian point. */ 854 855 int numgrids=7; //six plus bubble grids856 857 812 double r=gauss->coord2-gauss->coord1; 858 813 double s=-3.0/SQRT3*(gauss->coord1+gauss->coord2-2.0/3.0); … … 860 815 861 816 /*First nodal function: */ 862 *(dl1dl7+ numgrids*0+0)=-0.5*(1.0-zeta)/2.0;863 *(dl1dl7+ numgrids*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;864 *(dl1dl7+ numgrids*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);817 *(dl1dl7+NUMNODESMINI*0+0)=-0.5*(1.0-zeta)/2.0; 818 *(dl1dl7+NUMNODESMINI*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0; 819 *(dl1dl7+NUMNODESMINI*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD); 865 820 866 821 /*Second nodal function: */ 867 *(dl1dl7+ numgrids*0+1)=0.5*(1.0-zeta)/2.0;868 *(dl1dl7+ numgrids*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;869 *(dl1dl7+ numgrids*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);822 *(dl1dl7+NUMNODESMINI*0+1)=0.5*(1.0-zeta)/2.0; 823 *(dl1dl7+NUMNODESMINI*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0; 824 *(dl1dl7+NUMNODESMINI*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD); 870 825 871 826 /*Third nodal function: */ 872 *(dl1dl7+ numgrids*0+2)=0;873 *(dl1dl7+ numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;874 *(dl1dl7+ numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);827 *(dl1dl7+NUMNODESMINI*0+2)=0; 828 *(dl1dl7+NUMNODESMINI*1+2)=SQRT3/3.0*(1.0-zeta)/2.0; 829 *(dl1dl7+NUMNODESMINI*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD); 875 830 876 831 /*Fourth nodal function: */ 877 *(dl1dl7+ numgrids*0+3)=-0.5*(1.0+zeta)/2.0;878 *(dl1dl7+ numgrids*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;879 *(dl1dl7+ numgrids*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);832 *(dl1dl7+NUMNODESMINI*0+3)=-0.5*(1.0+zeta)/2.0; 833 *(dl1dl7+NUMNODESMINI*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0; 834 *(dl1dl7+NUMNODESMINI*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD); 880 835 881 836 /*Fith nodal function: */ 882 *(dl1dl7+ numgrids*0+4)=0.5*(1.0+zeta)/2.0;883 *(dl1dl7+ numgrids*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;884 *(dl1dl7+ numgrids*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);837 *(dl1dl7+NUMNODESMINI*0+4)=0.5*(1.0+zeta)/2.0; 838 *(dl1dl7+NUMNODESMINI*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0; 839 *(dl1dl7+NUMNODESMINI*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD); 885 840 886 841 /*Sixth nodal function: */ 887 *(dl1dl7+ numgrids*0+5)=0;888 *(dl1dl7+ numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;889 *(dl1dl7+ numgrids*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);842 *(dl1dl7+NUMNODESMINI*0+5)=0; 843 *(dl1dl7+NUMNODESMINI*1+5)=SQRT3/3.0*(1.0+zeta)/2.0; 844 *(dl1dl7+NUMNODESMINI*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD); 890 845 891 846 /*Seventh nodal function: */ 892 *(dl1dl7+ numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);893 *(dl1dl7+ numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));894 *(dl1dl7+ numgrids*2+6)=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);847 *(dl1dl7+NUMNODESMINI*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0); 848 *(dl1dl7+NUMNODESMINI*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0)); 849 *(dl1dl7+NUMNODESMINI*2+6)=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta); 895 850 896 851 } … … 914 869 /*This routine returns the values of the nodal functions derivatives (with respect to the 915 870 * actual coordinate system): */ 916 int i; 917 const int NDOF3 = 3; 918 const int numgrids = 6; 919 double dh1dh6_ref[NDOF3][numgrids]; 871 double dh1dh6_ref[NDOF3][NUMNODESP1]; 920 872 double Jinv[NDOF3][NDOF3]; 921 873 … … 933 885 */ 934 886 935 for (i =0;i<numgrids;i++){936 *(dh1dh6+ numgrids*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];937 *(dh1dh6+ numgrids*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];938 *(dh1dh6+ numgrids*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];887 for (int i=0;i<NUMNODESP1;i++){ 888 *(dh1dh6+NUMNODESP1*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i]; 889 *(dh1dh6+NUMNODESP1*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i]; 890 *(dh1dh6+NUMNODESP1*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i]; 939 891 } 940 892 … … 947 899 * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */ 948 900 949 const int numgrids=6;950 901 double A1,A2,A3,z; 951 902 … … 956 907 957 908 /*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/ 958 *(dl1dl6+ numgrids*0+0)=-0.5*(1.0-z)/2.0;959 *(dl1dl6+ numgrids*1+0)=-0.5/SQRT3*(1.0-z)/2.0;960 *(dl1dl6+ numgrids*2+0)=-0.5*A1;909 *(dl1dl6+NUMNODESP1*0+0)=-0.5*(1.0-z)/2.0; 910 *(dl1dl6+NUMNODESP1*1+0)=-0.5/SQRT3*(1.0-z)/2.0; 911 *(dl1dl6+NUMNODESP1*2+0)=-0.5*A1; 961 912 962 913 /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/ 963 *(dl1dl6+ numgrids*0+1)=0.5*(1.0-z)/2.0;964 *(dl1dl6+ numgrids*1+1)=-0.5/SQRT3*(1.0-z)/2.0;965 *(dl1dl6+ numgrids*2+1)=-0.5*A2;914 *(dl1dl6+NUMNODESP1*0+1)=0.5*(1.0-z)/2.0; 915 *(dl1dl6+NUMNODESP1*1+1)=-0.5/SQRT3*(1.0-z)/2.0; 916 *(dl1dl6+NUMNODESP1*2+1)=-0.5*A2; 966 917 967 918 /*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/ 968 *(dl1dl6+ numgrids*0+2)=0.0;969 *(dl1dl6+ numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0;970 *(dl1dl6+ numgrids*2+2)=-0.5*A3;919 *(dl1dl6+NUMNODESP1*0+2)=0.0; 920 *(dl1dl6+NUMNODESP1*1+2)=1.0/SQRT3*(1.0-z)/2.0; 921 *(dl1dl6+NUMNODESP1*2+2)=-0.5*A3; 971 922 972 923 /*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/ 973 *(dl1dl6+ numgrids*0+3)=-0.5*(1.0+z)/2.0;974 *(dl1dl6+ numgrids*1+3)=-0.5/SQRT3*(1.0+z)/2.0;975 *(dl1dl6+ numgrids*2+3)=0.5*A1;924 *(dl1dl6+NUMNODESP1*0+3)=-0.5*(1.0+z)/2.0; 925 *(dl1dl6+NUMNODESP1*1+3)=-0.5/SQRT3*(1.0+z)/2.0; 926 *(dl1dl6+NUMNODESP1*2+3)=0.5*A1; 976 927 977 928 /*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/ 978 *(dl1dl6+ numgrids*0+4)=0.5*(1.0+z)/2.0;979 *(dl1dl6+ numgrids*1+4)=-0.5/SQRT3*(1.0+z)/2.0;980 *(dl1dl6+ numgrids*2+4)=0.5*A2;929 *(dl1dl6+NUMNODESP1*0+4)=0.5*(1.0+z)/2.0; 930 *(dl1dl6+NUMNODESP1*1+4)=-0.5/SQRT3*(1.0+z)/2.0; 931 *(dl1dl6+NUMNODESP1*2+4)=0.5*A2; 981 932 982 933 /*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/ 983 *(dl1dl6+ numgrids*0+5)=0.0;984 *(dl1dl6+ numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0;985 *(dl1dl6+ numgrids*2+5)=0.5*A3;934 *(dl1dl6+NUMNODESP1*0+5)=0.0; 935 *(dl1dl6+NUMNODESP1*1+5)=1.0/SQRT3*(1.0+z)/2.0; 936 *(dl1dl6+NUMNODESP1*2+5)=0.5*A3; 986 937 } 987 938 /*}}}*/ … … 1010 961 * p is a vector of size 3x1 already allocated. 1011 962 */ 1012 double dh1dh6[3][ 6];963 double dh1dh6[3][NUMNODESP1]; 1013 964 1014 965 /*Get nodal funnctions derivatives in actual coordinate system: */ -
issm/trunk/src/c/objects/Elements/Tria.cpp
r5743 r5745 489 489 490 490 /*build new bed and surface: */ 491 if (this-> GetShelf()){491 if (this->IsOnShelf()){ 492 492 /*hydrostatic equilibrium: */ 493 493 double rho_ice,rho_water,di; … … 819 819 } 820 820 /*}}}*/ 821 /*FUNCTION Tria:: GetOnBed {{{1*/822 bool Tria:: GetOnBed(){821 /*FUNCTION Tria::IsOnBed {{{1*/ 822 bool Tria::IsOnBed(){ 823 823 824 824 bool onbed; … … 828 828 } 829 829 /*}}}*/ 830 /*FUNCTION Tria:: GetShelf {{{1*/831 bool Tria:: GetShelf(){830 /*FUNCTION Tria::IsOnShelf {{{1*/ 831 bool Tria::IsOnShelf(){ 832 832 /*inputs: */ 833 833 bool shelf; … … 2421 2421 2422 2422 /*If shelf: hydrostatic equilibrium*/ 2423 if (this-> GetShelf()){2423 if (this->IsOnShelf()){ 2424 2424 2425 2425 /*recover material parameters: */ -
issm/trunk/src/c/objects/Elements/Tria.h
r5743 r5745 75 75 void CreatePVector(Vec pg); 76 76 int GetNodeIndex(Node* node); 77 bool GetOnBed();78 bool GetShelf();77 bool IsOnBed(); 78 bool IsOnShelf(); 79 79 void GetSolutionFromInputs(Vec solution); 80 80 void GetVectorFromInputs(Vec vector,int NameEnum); -
issm/trunk/src/c/objects/Elements/TriaRef.cpp
r5743 r5745 20 20 /*}}}*/ 21 21 22 /*Element macros*/ 23 #define NUMNODES 3 24 #define NDOF1 1 25 #define NDOF2 2 26 #define NDOF3 3 27 #define NDOF4 4 28 22 29 /*Object constructors and destructor*/ 23 30 /*FUNCTION TriaRef::TriaRef(){{{1*/ … … 27 34 /*}}}*/ 28 35 /*FUNCTION TriaRef::TriaRef(int* types,int nummodels){{{1*/ 36 29 37 TriaRef::TriaRef(const int nummodels){ 30 38 … … 62 70 * where h is the interpolation function for grid i. 63 71 * 64 * We assume B has been allocated already, of size: 3x(NDOF2* numgrids)72 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES) 65 73 */ 66 74 67 75 int i; 68 const int NDOF2=2; 69 const int numgrids=3; 70 71 double dh1dh3[NDOF2][numgrids]; 72 76 double dh1dh3[NDOF2][NUMNODES]; 73 77 74 78 /*Get dh1dh2dh3 in actual coordinate system: */ … … 76 80 77 81 /*Build B: */ 78 for (i=0;i< numgrids;i++){79 *(B+NDOF2* numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];80 *(B+NDOF2* numgrids*0+NDOF2*i+1)=0;81 *(B+NDOF2* numgrids*1+NDOF2*i)=0;82 *(B+NDOF2* numgrids*1+NDOF2*i+1)=dh1dh3[1][i];83 *(B+NDOF2* numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i];84 *(B+NDOF2* numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i];82 for (i=0;i<NUMNODES;i++){ 83 *(B+NDOF2*NUMNODES*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i]; 84 *(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0; 85 *(B+NDOF2*NUMNODES*1+NDOF2*i)=0; 86 *(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dh1dh3[1][i]; 87 *(B+NDOF2*NUMNODES*2+NDOF2*i)=(float).5*dh1dh3[1][i]; 88 *(B+NDOF2*NUMNODES*2+NDOF2*i+1)=(float).5*dh1dh3[0][i]; 85 89 } 86 90 } … … 95 99 */ 96 100 97 const int numgrids=3; 98 double l1l3[numgrids]; 101 double l1l3[NUMNODES]; 99 102 100 103 GetNodalFunctions(&l1l3[0],gauss); … … 115 118 */ 116 119 117 const int numgrids=3; 118 double l1l3[numgrids]; 120 double l1l3[NUMNODES]; 119 121 120 122 GetNodalFunctions(&l1l3[0],gauss); … … 135 137 * where h is the interpolation function for grid i. 136 138 * 137 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids) 138 */ 139 140 int i; 141 const int NDOF1=1; 142 const int numgrids=3; 143 144 double l1l2l3[numgrids]; 145 139 * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODES) 140 */ 141 142 double l1l2l3[NUMNODES]; 146 143 147 144 /*Get dh1dh2dh3 in actual coordinate system: */ … … 149 146 150 147 /*Build B_prog: */ 151 for (i =0;i<numgrids;i++){152 *(B_prog+NDOF1* numgrids*0+NDOF1*i)=l1l2l3[i];153 *(B_prog+NDOF1* numgrids*1+NDOF1*i)=l1l2l3[i];148 for (int i=0;i<NUMNODES;i++){ 149 *(B_prog+NDOF1*NUMNODES*0+NDOF1*i)=l1l2l3[i]; 150 *(B_prog+NDOF1*NUMNODES*1+NDOF1*i)=l1l2l3[i]; 154 151 } 155 152 } … … 166 163 * where h is the interpolation function for grid i. 167 164 * 168 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids) 169 */ 170 171 int i; 172 const int NDOF2=2; 173 const int numgrids=3; 165 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES) 166 */ 174 167 175 168 /*Same thing in the actual coordinate system: */ 176 double dh1dh3[NDOF2][ numgrids];169 double dh1dh3[NDOF2][NUMNODES]; 177 170 178 171 /*Get dh1dh2dh3 in actual coordinates system : */ … … 180 173 181 174 /*Build B': */ 182 for (i =0;i<numgrids;i++){183 *(Bprime+NDOF2* numgrids*0+NDOF2*i)=2*dh1dh3[0][i];184 *(Bprime+NDOF2* numgrids*0+NDOF2*i+1)=dh1dh3[1][i];185 *(Bprime+NDOF2* numgrids*1+NDOF2*i)=dh1dh3[0][i];186 *(Bprime+NDOF2* numgrids*1+NDOF2*i+1)=2*dh1dh3[1][i];187 *(Bprime+NDOF2* numgrids*2+NDOF2*i)=dh1dh3[1][i];188 *(Bprime+NDOF2* numgrids*2+NDOF2*i+1)=dh1dh3[0][i];175 for (int i=0;i<NUMNODES;i++){ 176 *(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=2*dh1dh3[0][i]; 177 *(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=dh1dh3[1][i]; 178 *(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=dh1dh3[0][i]; 179 *(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=2*dh1dh3[1][i]; 180 *(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dh1dh3[1][i]; 181 *(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dh1dh3[0][i]; 189 182 } 190 183 } … … 199 192 * where h is the interpolation function for grid i. 200 193 * 201 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids) 202 */ 203 204 int i; 205 const int NDOF1=1; 206 const int NDOF2=2; 207 const int numgrids=3; 194 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES) 195 */ 208 196 209 197 /*Same thing in the actual coordinate system: */ 210 double dh1dh3[NDOF2][ numgrids];198 double dh1dh3[NDOF2][NUMNODES]; 211 199 212 200 /*Get dh1dh2dh3 in actual coordinates system : */ … … 214 202 215 203 /*Build B': */ 216 for (i =0;i<numgrids;i++){217 *(Bprime_prog+NDOF1* numgrids*0+NDOF1*i)=dh1dh3[0][i];218 *(Bprime_prog+NDOF1* numgrids*1+NDOF1*i)=dh1dh3[1][i];204 for (int i=0;i<NUMNODES;i++){ 205 *(Bprime_prog+NDOF1*NUMNODES*0+NDOF1*i)=dh1dh3[0][i]; 206 *(Bprime_prog+NDOF1*NUMNODES*1+NDOF1*i)=dh1dh3[1][i]; 219 207 } 220 208 } … … 232 220 * where h is the interpolation function for grid i. 233 221 * 234 * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2)222 * We assume L has been allocated already, of size: NUMNODES (numdof=1), or numdofx(numdof*NUMNODES) (numdof=2) 235 223 */ 236 224 237 225 int i; 238 const int NDOF2=2;239 const int numgrids=3;240 226 double l1l2l3[3]; 241 227 … … 245 231 /*Build L: */ 246 232 if(numdof==1){ 247 for (i=0;i< numgrids;i++){233 for (i=0;i<NUMNODES;i++){ 248 234 L[i]=l1l2l3[i]; 249 235 } 250 236 } 251 237 else{ 252 for (i=0;i< numgrids;i++){253 *(L+numdof* numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];254 *(L+numdof* numgrids*0+numdof*i+1)=0;255 *(L+numdof* numgrids*1+numdof*i)=0;256 *(L+numdof* numgrids*1+numdof*i+1)=l1l2l3[i];238 for (i=0;i<NUMNODES;i++){ 239 *(L+numdof*NUMNODES*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i]; 240 *(L+numdof*NUMNODES*0+numdof*i+1)=0; 241 *(L+numdof*NUMNODES*1+numdof*i)=0; 242 *(L+numdof*NUMNODES*1+numdof*i+1)=l1l2l3[i]; 257 243 } 258 244 } … … 263 249 /*The Jacobian is constant over the element, discard the gaussian points. 264 250 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 265 266 const int NDOF2=2;267 const int numgrids=3;268 251 double x1,y1,x2,y2,x3,y3; 269 252 270 x1=*(xyz_list+ numgrids*0+0);271 y1=*(xyz_list+ numgrids*0+1);272 x2=*(xyz_list+ numgrids*1+0);273 y2=*(xyz_list+ numgrids*1+1);274 x3=*(xyz_list+ numgrids*2+0);275 y3=*(xyz_list+ numgrids*2+1);253 x1=*(xyz_list+NUMNODES*0+0); 254 y1=*(xyz_list+NUMNODES*0+1); 255 x2=*(xyz_list+NUMNODES*1+0); 256 y2=*(xyz_list+NUMNODES*1+1); 257 x3=*(xyz_list+NUMNODES*2+0); 258 y3=*(xyz_list+NUMNODES*2+1); 276 259 277 260 … … 388 371 * actual coordinate system): */ 389 372 int i; 390 const int NDOF2 = 2; 391 const int numgrids = 3; 392 double dh1dh3_ref[NDOF2][numgrids]; 373 double dh1dh3_ref[NDOF2][NUMNODES]; 393 374 double Jinv[NDOF2][NDOF2]; 394 375 … … 404 385 * [dhi/dy] [dhi/ds] 405 386 */ 406 for (i=0;i< numgrids;i++){407 dh1dh3[ numgrids*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];408 dh1dh3[ numgrids*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];387 for (i=0;i<NUMNODES;i++){ 388 dh1dh3[NUMNODES*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i]; 389 dh1dh3[NUMNODES*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i]; 409 390 } 410 391 … … 413 394 /*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{1*/ 414 395 void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss){ 415 416 396 /*This routine returns the values of the nodal functions derivatives (with respect to the 417 397 * natural coordinate system) at the gaussian point. */ 418 398 419 const int NDOF2=2;420 const int numgrids=3;421 422 399 /*First nodal function: */ 423 *(dl1dl3+ numgrids*0+0)=-0.5;424 *(dl1dl3+ numgrids*1+0)=-1.0/(2.0*SQRT3);400 *(dl1dl3+NUMNODES*0+0)=-0.5; 401 *(dl1dl3+NUMNODES*1+0)=-1.0/(2.0*SQRT3); 425 402 426 403 /*Second nodal function: */ 427 *(dl1dl3+ numgrids*0+1)=0.5;428 *(dl1dl3+ numgrids*1+1)=-1.0/(2.0*SQRT3);404 *(dl1dl3+NUMNODES*0+1)=0.5; 405 *(dl1dl3+NUMNODES*1+1)=-1.0/(2.0*SQRT3); 429 406 430 407 /*Third nodal function: */ 431 *(dl1dl3+ numgrids*0+2)=0;432 *(dl1dl3+ numgrids*1+2)=1.0/SQRT3;408 *(dl1dl3+NUMNODES*0+2)=0; 409 *(dl1dl3+NUMNODES*1+2)=1.0/SQRT3; 433 410 434 411 }
Note:
See TracChangeset
for help on using the changeset viewer.