Changeset 14767
- Timestamp:
- 04/26/13 10:07:43 (12 years ago)
- Location:
- issm/trunk-jpl/src/c/classes/objects/Elements
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/objects/Elements/TriaRef.cpp
r14591 r14767 21 21 22 22 /*Element macros*/ 23 #define NUMNODES 3 23 #define NUMNODESP1 3 24 #define NUMNODESP2 6 24 25 25 26 /*Object constructors and destructor*/ … … 65 66 * where h is the interpolation function for node i. 66 67 * 67 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES )68 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1) 68 69 */ 69 70 70 71 int i; 71 IssmDouble dbasis[NDOF2][NUMNODES ];72 IssmDouble dbasis[NDOF2][NUMNODESP1]; 72 73 73 74 /*Get dh1dh2dh3 in actual coordinate system: */ … … 75 76 76 77 /*Build B: */ 77 for (i=0;i<NUMNODES ;i++){78 B[NDOF1*NUMNODES *0+NDOF1*i]=dbasis[0][i];79 B[NDOF1*NUMNODES *1+NDOF1*i]=dbasis[1][i];78 for (i=0;i<NUMNODESP1;i++){ 79 B[NDOF1*NUMNODESP1*0+NDOF1*i]=dbasis[0][i]; 80 B[NDOF1*NUMNODESP1*1+NDOF1*i]=dbasis[1][i]; 80 81 } 81 82 } … … 91 92 * where h is the interpolation function for node i. 92 93 * 93 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES )94 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1) 94 95 */ 95 96 96 97 int i; 97 IssmDouble dbasis[NDOF2][NUMNODES ];98 IssmDouble dbasis[NDOF2][NUMNODESP1]; 98 99 99 100 /*Get dh1dh2dh3 in actual coordinate system: */ … … 101 102 102 103 /*Build B: */ 103 for (i=0;i<NUMNODES ;i++){104 *(B+NDOF2*NUMNODES *0+NDOF2*i)=dbasis[0][i]; //B[0][NDOF2*i]=dbasis[0][i];105 *(B+NDOF2*NUMNODES *0+NDOF2*i+1)=0;106 *(B+NDOF2*NUMNODES *1+NDOF2*i)=0;107 *(B+NDOF2*NUMNODES *1+NDOF2*i+1)=dbasis[1][i];108 *(B+NDOF2*NUMNODES *2+NDOF2*i)=(float).5*dbasis[1][i];109 *(B+NDOF2*NUMNODES *2+NDOF2*i+1)=(float).5*dbasis[0][i];104 for (i=0;i<NUMNODESP1;i++){ 105 *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i]; //B[0][NDOF2*i]=dbasis[0][i]; 106 *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0; 107 *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0; 108 *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i]; 109 *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i]; 110 *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i]; 110 111 } 111 112 } … … 122 123 * where h is the interpolation function for node i. 123 124 * 124 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES )125 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1) 125 126 */ 126 127 127 128 /*Same thing in the actual coordinate system: */ 128 IssmDouble dbasis[NDOF2][NUMNODES ];129 IssmDouble dbasis[NDOF2][NUMNODESP1]; 129 130 130 131 /*Get dh1dh2dh3 in actual coordinates system : */ … … 132 133 133 134 /*Build B': */ 134 for (int i=0;i<NUMNODES ;i++){135 *(B+NDOF2*NUMNODES *0+NDOF2*i)=dbasis[0][i];136 *(B+NDOF2*NUMNODES *0+NDOF2*i+1)=0;137 *(B+NDOF2*NUMNODES *1+NDOF2*i)=0;138 *(B+NDOF2*NUMNODES *1+NDOF2*i+1)=dbasis[1][i];139 *(B+NDOF2*NUMNODES *2+NDOF2*i)=0.5*dbasis[1][i];140 *(B+NDOF2*NUMNODES *2+NDOF2*i+1)=0.5*dbasis[0][i];135 for (int i=0;i<NUMNODESP1;i++){ 136 *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i]; 137 *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0; 138 *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0; 139 *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i]; 140 *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=0.5*dbasis[1][i]; 141 *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=0.5*dbasis[0][i]; 141 142 } 142 143 } … … 151 152 */ 152 153 153 IssmDouble l1l3[NUMNODES ];154 IssmDouble l1l3[NUMNODESP1]; 154 155 155 156 GetNodalFunctions(&l1l3[0],gauss); … … 170 171 */ 171 172 172 IssmDouble l1l3[NUMNODES ];173 IssmDouble l1l3[NUMNODESP1]; 173 174 174 175 GetNodalFunctions(&l1l3[0],gauss); … … 189 190 * where h is the interpolation function for node i. 190 191 * 191 * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODES )192 */ 193 194 IssmDouble basis[NUMNODES ];192 * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODESP1) 193 */ 194 195 IssmDouble basis[NUMNODESP1]; 195 196 196 197 /*Get dh1dh2dh3 in actual coordinate system: */ … … 198 199 199 200 /*Build B_prog: */ 200 for (int i=0;i<NUMNODES ;i++){201 *(B_prog+NDOF1*NUMNODES *0+NDOF1*i)=basis[i];202 *(B_prog+NDOF1*NUMNODES *1+NDOF1*i)=basis[i];201 for (int i=0;i<NUMNODESP1;i++){ 202 *(B_prog+NDOF1*NUMNODESP1*0+NDOF1*i)=basis[i]; 203 *(B_prog+NDOF1*NUMNODESP1*1+NDOF1*i)=basis[i]; 203 204 } 204 205 } … … 215 216 * where h is the interpolation function for node i. 216 217 * 217 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES )218 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODESP1) 218 219 */ 219 220 220 221 /*Same thing in the actual coordinate system: */ 221 IssmDouble dbasis[NDOF2][NUMNODES ];222 IssmDouble dbasis[NDOF2][NUMNODESP1]; 222 223 223 224 /*Get dh1dh2dh3 in actual coordinates system : */ … … 225 226 226 227 /*Build B': */ 227 for (int i=0;i<NUMNODES ;i++){228 *(Bprime+NDOF2*NUMNODES *0+NDOF2*i)=2*dbasis[0][i];229 *(Bprime+NDOF2*NUMNODES *0+NDOF2*i+1)=dbasis[1][i];230 *(Bprime+NDOF2*NUMNODES *1+NDOF2*i)=dbasis[0][i];231 *(Bprime+NDOF2*NUMNODES *1+NDOF2*i+1)=2*dbasis[1][i];232 *(Bprime+NDOF2*NUMNODES *2+NDOF2*i)=dbasis[1][i];233 *(Bprime+NDOF2*NUMNODES *2+NDOF2*i+1)=dbasis[0][i];228 for (int i=0;i<NUMNODESP1;i++){ 229 *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i)=2*dbasis[0][i]; 230 *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dbasis[1][i]; 231 *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i)=dbasis[0][i]; 232 *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2*dbasis[1][i]; 233 *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i]; 234 *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i]; 234 235 } 235 236 } … … 247 248 * where h is the interpolation function for node i. 248 249 * 249 * We assume Bprime has been allocated already, of size: 3x(NDOF2*NUMNODES )250 * We assume Bprime has been allocated already, of size: 3x(NDOF2*NUMNODESP1) 250 251 */ 251 252 252 253 /*Same thing in the actual coordinate system: */ 253 IssmDouble dbasis[NDOF2][NUMNODES ];254 IssmDouble dbasis[NDOF2][NUMNODESP1]; 254 255 255 256 /*Get dh1dh2dh3 in actual coordinates system : */ … … 257 258 258 259 /*Build Bprime: */ 259 for (int i=0;i<NUMNODES ;i++){260 *(Bprime+NDOF2*NUMNODES *0+NDOF2*i)=dbasis[0][i];261 *(Bprime+NDOF2*NUMNODES *0+NDOF2*i+1)=0;262 *(Bprime+NDOF2*NUMNODES *1+NDOF2*i)=0;263 *(Bprime+NDOF2*NUMNODES *1+NDOF2*i+1)=dbasis[1][i];264 *(Bprime+NDOF2*NUMNODES *2+NDOF2*i)=dbasis[1][i];265 *(Bprime+NDOF2*NUMNODES *2+NDOF2*i+1)=dbasis[0][i];266 *(Bprime+NDOF2*NUMNODES *3+NDOF2*i)=dbasis[0][i];267 *(Bprime+NDOF2*NUMNODES *3+NDOF2*i+1)=dbasis[1][i];260 for (int i=0;i<NUMNODESP1;i++){ 261 *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i]; 262 *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0; 263 *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i)=0; 264 *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i]; 265 *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i]; 266 *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i]; 267 *(Bprime+NDOF2*NUMNODESP1*3+NDOF2*i)=dbasis[0][i]; 268 *(Bprime+NDOF2*NUMNODESP1*3+NDOF2*i+1)=dbasis[1][i]; 268 269 } 269 270 } … … 278 279 * where h is the interpolation function for node i. 279 280 * 280 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES )281 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODESP1) 281 282 */ 282 283 283 284 /*Same thing in the actual coordinate system: */ 284 IssmDouble dbasis[NDOF2][NUMNODES ];285 IssmDouble dbasis[NDOF2][NUMNODESP1]; 285 286 286 287 /*Get dh1dh2dh3 in actual coordinates system : */ … … 288 289 289 290 /*Build B': */ 290 for (int i=0;i<NUMNODES ;i++){291 *(Bprime_prog+NDOF1*NUMNODES *0+NDOF1*i)=dbasis[0][i];292 *(Bprime_prog+NDOF1*NUMNODES *1+NDOF1*i)=dbasis[1][i];291 for (int i=0;i<NUMNODESP1;i++){ 292 *(Bprime_prog+NDOF1*NUMNODESP1*0+NDOF1*i)=dbasis[0][i]; 293 *(Bprime_prog+NDOF1*NUMNODESP1*1+NDOF1*i)=dbasis[1][i]; 293 294 } 294 295 } … … 306 307 * where h is the interpolation function for node i. 307 308 * 308 * We assume L has been allocated already, of size: NUMNODES (numdof=1), or numdofx(numdof*NUMNODES) (numdof=2)309 * We assume L has been allocated already, of size: NUMNODESP1 (numdof=1), or numdofx(numdof*NUMNODESP1) (numdof=2) 309 310 */ 310 311 … … 317 318 /*Build L: */ 318 319 if(numdof==1){ 319 for (i=0;i<NUMNODES ;i++){320 for (i=0;i<NUMNODESP1;i++){ 320 321 L[i]=basis[i]; 321 322 } 322 323 } 323 324 else{ 324 for (i=0;i<NUMNODES ;i++){325 *(L+numdof*NUMNODES *0+numdof*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];326 *(L+numdof*NUMNODES *0+numdof*i+1)=0;327 *(L+numdof*NUMNODES *1+numdof*i)=0;328 *(L+numdof*NUMNODES *1+numdof*i+1)=basis[i];325 for (i=0;i<NUMNODESP1;i++){ 326 *(L+numdof*NUMNODESP1*0+numdof*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i]; 327 *(L+numdof*NUMNODESP1*0+numdof*i+1)=0; 328 *(L+numdof*NUMNODESP1*1+numdof*i)=0; 329 *(L+numdof*NUMNODESP1*1+numdof*i+1)=basis[i]; 329 330 } 330 331 } … … 337 338 IssmDouble x1,y1,x2,y2,x3,y3; 338 339 339 x1=*(xyz_list+NUMNODES *0+0);340 y1=*(xyz_list+NUMNODES *0+1);341 x2=*(xyz_list+NUMNODES *1+0);342 y2=*(xyz_list+NUMNODES *1+1);343 x3=*(xyz_list+NUMNODES *2+0);344 y3=*(xyz_list+NUMNODES *2+1);340 x1=*(xyz_list+NUMNODESP1*0+0); 341 y1=*(xyz_list+NUMNODESP1*0+1); 342 x2=*(xyz_list+NUMNODESP1*1+0); 343 y2=*(xyz_list+NUMNODESP1*1+1); 344 x3=*(xyz_list+NUMNODESP1*2+0); 345 y3=*(xyz_list+NUMNODESP1*2+1); 345 346 346 347 *(J+NDOF2*0+0)=0.5*(x2-x1); … … 405 406 } 406 407 /*}}}*/ 408 /*FUNCTION TriaRef::GetNodalFunctionsP2{{{*/ 409 void TriaRef::GetNodalFunctionsP2(IssmDouble** pbasis,GaussTria* gauss){ 410 /*This routine returns the values of the nodal functions at the gaussian point.*/ 411 412 /*Allocate*/ 413 IssmDouble* basis =xNew<IssmDouble>(NUMNODESP2); 414 415 basis[0]=gauss->coord1*(2.*gauss->coord1-1.); 416 basis[1]=gauss->coord2*(2.*gauss->coord2-1.); 417 basis[2]=gauss->coord3*(2.*gauss->coord3-1.); 418 basis[3]=4.*gauss->coord3*gauss->coord2; 419 basis[4]=4.*gauss->coord3*gauss->coord1; 420 basis[5]=4.*gauss->coord1*gauss->coord2; 421 422 /*Assign output pointer*/ 423 *pbasis = basis; 424 425 } 426 /*}}}*/ 407 427 /*FUNCTION TriaRef::GetSegmentNodalFunctions{{{*/ 408 428 void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss,int index1,int index2){ … … 425 445 * actual coordinate system): */ 426 446 int i; 427 IssmDouble dbasis_ref[NDOF2][NUMNODES ];447 IssmDouble dbasis_ref[NDOF2][NUMNODESP1]; 428 448 IssmDouble Jinv[NDOF2][NDOF2]; 429 449 … … 439 459 * [dhi/dy] [dhi/ds] 440 460 */ 441 for (i=0;i<NUMNODES ;i++){442 dbasis[NUMNODES *0+i]=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i];443 dbasis[NUMNODES *1+i]=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i];461 for (i=0;i<NUMNODESP1;i++){ 462 dbasis[NUMNODESP1*0+i]=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i]; 463 dbasis[NUMNODESP1*1+i]=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i]; 444 464 } 445 465 … … 452 472 453 473 /*First nodal function: */ 454 *(dl1dl3+NUMNODES *0+0)=-0.5;455 *(dl1dl3+NUMNODES *1+0)=-1.0/(2.0*SQRT3);474 *(dl1dl3+NUMNODESP1*0+0)=-0.5; 475 *(dl1dl3+NUMNODESP1*1+0)=-1.0/(2.0*SQRT3); 456 476 457 477 /*Second nodal function: */ 458 *(dl1dl3+NUMNODES *0+1)=0.5;459 *(dl1dl3+NUMNODES *1+1)=-1.0/(2.0*SQRT3);478 *(dl1dl3+NUMNODESP1*0+1)=0.5; 479 *(dl1dl3+NUMNODESP1*1+1)=-1.0/(2.0*SQRT3); 460 480 461 481 /*Third nodal function: */ 462 *(dl1dl3+NUMNODES*0+2)=0; 463 *(dl1dl3+NUMNODES*1+2)=1.0/SQRT3; 464 482 *(dl1dl3+NUMNODESP1*0+2)=0; 483 *(dl1dl3+NUMNODESP1*1+2)=1.0/SQRT3; 484 485 } 486 /*}}}*/ 487 /*FUNCTION TriaRef::GetNodalFunctionsDerivativesP2Reference{{{*/ 488 void TriaRef::GetNodalFunctionsDerivativesP2Reference(IssmDouble** pdbasis,GaussTria* gauss){ 489 /*This routine returns the values of the nodal functions derivatives (with respect to the 490 * natural coordinate system) at the gaussian point. */ 491 492 /*Allocate*/ 493 IssmDouble* dbasis =xNew<IssmDouble>(NUMNODESP2*2); 494 495 /*Nodal function 1*/ 496 dbasis[NUMNODESP2*0+0]=-2.*gauss->coord1 + 0.5; 497 dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.; 498 /*Nodal function 2*/ 499 dbasis[NUMNODESP2*0+0]=+2.*gauss->coord2 + 0.5; 500 dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.; 501 /*Nodal function 3*/ 502 dbasis[NUMNODESP2*0+0]=0.; 503 dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.; 504 /*Nodal function 4*/ 505 dbasis[NUMNODESP2*0+0]=+2.*gauss->coord3; 506 dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3; 507 /*Nodal function 5*/ 508 dbasis[NUMNODESP2*0+0]=-2.*gauss->coord3; 509 dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3; 510 /*Nodal function 6*/ 511 dbasis[NUMNODESP2*0+0]=2.*(gauss->coord1-gauss->coord2); 512 dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*(gauss->coord1+gauss->coord2); 513 514 /*Assign output pointer*/ 515 *pdbasis = dbasis; 465 516 } 466 517 /*}}}*/ -
issm/trunk-jpl/src/c/classes/objects/Elements/TriaRef.h
r14591 r14767 35 35 void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss); 36 36 void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTria* gauss); 37 void GetNodalFunctions(IssmDouble* l1l2l3,GaussTria* gauss); 38 void GetSegmentNodalFunctions(IssmDouble* l1l2l3,GaussTria* gauss, int index1,int index2); 37 void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss); 38 void GetNodalFunctionsP2(IssmDouble** pbasis,GaussTria* gauss); 39 void GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss, int index1,int index2); 39 40 void GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2); 40 41 void GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2); 41 void GetNodalFunctionsDerivatives(IssmDouble* l1l2l3,IssmDouble* xyz_list, GaussTria* gauss);42 void GetNodalFunctionsDerivatives(IssmDouble* basis,IssmDouble* xyz_list, GaussTria* gauss); 42 43 void GetNodalFunctionsDerivativesReference(IssmDouble* dl1dl3,GaussTria* gauss); 44 void GetNodalFunctionsDerivativesP2Reference(IssmDouble** pdbasis,GaussTria* gauss); 43 45 void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss); 44 46 void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);
Note:
See TracChangeset
for help on using the changeset viewer.