Changeset 5726
- Timestamp:
- 09/09/10 15:02:52 (15 years ago)
- Location:
- issm/trunk/src/c/objects/Elements
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r5724 r5726 2740 2740 int ig; 2741 2741 GaussPenta *gauss=NULL; 2742 2743 /* specific gaussian point: */2744 double gauss_weight,area_gauss_weight,vert_gauss_weight;2745 double gauss_coord[4];2746 double gauss_coord_tria[3];2747 int area_order=5;2748 int num_vert_gauss=5;2749 2750 2742 double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ 2751 2743 double viscosity; … … 3416 3408 double pe_g[numdofs]={0.0}; 3417 3409 double xyz_list[NUMVERTICES][3]; 3410 double xyz_list_segment[2][3]; 3418 3411 double z_list[NUMVERTICES]; 3419 3412 double z_segment[2]; 3420 3413 double slope2,constant_part; 3421 double gauss[numdofs][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};3422 3414 int node0,node1; 3423 BeamRef* beam=NULL; 3424 3425 /*gaussian points: */ 3426 int num_gauss; 3427 double* segment_gauss_coord=NULL; 3428 double gauss_coord; 3429 double* gauss_weights=NULL; 3430 double gauss_weight; 3415 GaussPenta* gauss=NULL; 3431 3416 int ig; 3432 3417 double slope[2]; … … 3444 3429 int connectivity[2]; 3445 3430 3446 /*Inputs*/3447 Input* thickness_input;3448 Input* surface_input;3449 Input* slopex_input;3450 Input* slopey_input;3451 3452 3431 /*recover some inputs: */ 3453 3432 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); … … 3469 3448 B=matice->GetB(); 3470 3449 3471 //Initilaize beam3472 beam=new BeamRef();3473 3474 3450 //recover extra inputs 3475 thickness_input=inputs->GetInput(ThicknessEnum);3476 surface_input=inputs->GetInput(SurfaceEnum);3477 slopex_input=inputs->GetInput(SurfaceSlopeXEnum);3478 slopey_input=inputs->GetInput(SurfaceSlopeYEnum);3451 Input* thickness_input=inputs->GetInput(ThicknessEnum); 3452 Input* surface_input=inputs->GetInput(SurfaceEnum); 3453 Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); 3454 Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); 3479 3455 3480 3456 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 3481 3457 for(i=0;i<NUMVERTICES;i++)z_list[i]=xyz_list[i][2]; 3482 3458 3483 //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions3484 num_gauss=3;3485 gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);3486 3487 3459 /*Loop on the three segments*/ 3488 3460 for(i=0;i<3;i++){ 3489 slopex_input->GetParameterValue(&slope[0], &gauss[i][0]);3490 slopey_input->GetParameterValue(&slope[1], &gauss[i][0]);3491 surface_input->GetParameterValue(&surface, &gauss[i][0]);3492 thickness_input->GetParameterValue(&thickness, &gauss[i][0]);3493 3494 //compute slope2 slopex and slopey3495 slope2=pow(slope[0],2)+pow(slope[1],2);3496 3497 //%compute constant_part3498 constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));3499 3500 z_segment[0]=z_list[i];3501 z_segment[1]=z_list[i+3];3502 3461 3503 3462 node0=i; 3504 3463 node1=i+3; 3505 3464 3465 for(j=0;j<3;j++){ 3466 xyz_list_segment[0][j]=xyz_list[node0][j]; 3467 xyz_list_segment[1][j]=xyz_list[node1][j]; 3468 } 3469 3506 3470 connectivity[0]=nodes[node0]->GetConnectivity(); 3507 3471 connectivity[1]=nodes[node1]->GetConnectivity(); 3508 3472 3509 3473 /*Loop on the Gauss points: */ 3510 for(ig=0;ig<num_gauss;ig++){ 3511 3512 //Pick up the gaussian point and its weight: 3513 gauss_weight=gauss_weights[ig]; 3514 gauss_coord=segment_gauss_coord[ig]; 3515 3516 beam->GetParameterValue(&z_g, &z_segment[0],gauss_coord); 3517 3518 //Get Jacobian determinant: 3519 beam->GetJacobianDeterminant(&Jdet, &z_segment[0],gauss_coord); 3474 gauss=new GaussPenta(node0,node1,3); 3475 for(ig=gauss->begin();ig<gauss->end();ig++){ 3476 gauss->GaussPoint(ig); 3477 3478 slopex_input->GetParameterValue(&slope[0],gauss); 3479 slopey_input->GetParameterValue(&slope[1],gauss); 3480 surface_input->GetParameterValue(&surface,gauss); 3481 thickness_input->GetParameterValue(&thickness,gauss); 3482 3483 slope2=pow(slope[0],2)+pow(slope[1],2); 3484 constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2)); 3485 3486 PentaRef::GetParameterValue(&z_g,&z_list[0],gauss); 3487 GetSegmentJacobianDeterminant(&Jdet,&xyz_list_segment[0][0],gauss); 3520 3488 3521 3489 if (onsurface){ 3522 for(j=0;j<NDOF2;j++){ 3523 pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight/(double)connectivity[1]; 3524 } 3490 for(j=0;j<NDOF2;j++) pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight/(double)connectivity[1]; 3525 3491 } 3526 3492 else{//connectivity is too large, should take only half on it 3527 for(j=0;j<NDOF2;j++){ 3528 pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight*2/(double)connectivity[1]; 3529 } 3493 for(j=0;j<NDOF2;j++) pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight*2/(double)connectivity[1]; 3530 3494 } 3531 3495 } 3496 delete gauss; 3532 3497 3533 3498 //Deal with lower surface 3534 3499 if (onbed){ 3535 3536 //compute ub3537 3500 constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thickness; 3538 3501 ub=constant_part*slope[0]; 3539 3502 vb=constant_part*slope[1]; 3540 3503 3541 //Add to pe:3542 3504 pe_g[2*node0]+=ub/(double)connectivity[0]; 3543 3505 pe_g[2*node0+1]+=vb/(double)connectivity[0]; … … 3549 3511 3550 3512 /*Clean up */ 3551 delete beam;3552 xfree((void**)&gauss_weights);3553 xfree((void**)&segment_gauss_coord);3554 3513 xfree((void**)&doflist); 3555 3514 } -
issm/trunk/src/c/objects/Elements/PentaRef.cpp
r5722 r5726 1699 1699 } 1700 1700 /*}}}*/ 1701 /*FUNCTION PentaRef::GetSegmentJacobianDeterminant{{{1*/ 1702 void PentaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss){ 1703 /*The Jacobian determinant is constant over the element, discard the gaussian points. 1704 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 1705 1706 double x1,x2,y1,y2,z1,z2; 1707 1708 x1=*(xyz_list+3*0+0); 1709 y1=*(xyz_list+3*0+1); 1710 z1=*(xyz_list+3*0+2); 1711 x2=*(xyz_list+3*1+0); 1712 y2=*(xyz_list+3*1+1); 1713 z2=*(xyz_list+3*1+2); 1714 1715 *Jdet=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.) + pow(z2-z1,2.)); 1716 if(*Jdet<0) ISSMERROR("negative jacobian determinant!"); 1717 1718 } 1719 /*}}}*/ 1701 1720 /*FUNCTION PentaRef::GetJacobianInvert {{{1*/ 1702 1721 void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss){ -
issm/trunk/src/c/objects/Elements/PentaRef.h
r5722 r5726 59 59 void GetJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss); 60 60 void GetTriaJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss); 61 void GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss); 61 62 void GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss); 62 63 void GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss);
Note:
See TracChangeset
for help on using the changeset viewer.