Changeset 5722
- Timestamp:
- 09/09/10 14:17:56 (15 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r5719 r5722 508 508 int dofp[1]={3}; 509 509 double Jdet2d; 510 Tria* tria=NULL;511 510 512 511 /*Gauss*/ 513 int num_gauss,ig; 514 double* first_gauss_area_coord = NULL; 515 double* second_gauss_area_coord = NULL; 516 double* third_gauss_area_coord = NULL; 517 double* gauss_weights = NULL; 518 double gauss_weight; 519 double gauss_coord[4]; 512 int ig; 513 GaussPenta* gauss; 520 514 521 515 /*Output*/ … … 529 523 bool onbed; 530 524 int approximation; 531 Input* pressure_input=NULL;532 Input* vx_input=NULL;533 Input* vy_input=NULL;534 Input* vz_input=NULL;535 525 536 526 /*parameters: */ … … 569 559 } 570 560 571 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */572 GaussLegendreTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2);573 574 561 /*Retrieve all inputs we will be needing: */ 575 pressure_input=inputs->GetInput(PressureEnum);576 vx_input=inputs->GetInput(VxEnum);577 vy_input=inputs->GetInput(VyEnum);578 vz_input=inputs->GetInput(VzEnum);562 Input* pressure_input=inputs->GetInput(PressureEnum); 563 Input* vx_input=inputs->GetInput(VxEnum); 564 Input* vy_input=inputs->GetInput(VyEnum); 565 Input* vz_input=inputs->GetInput(VzEnum); 579 566 580 567 /* Start looping on the number of gaussian points: */ 568 gauss=new GaussPenta(0,1,2,2); 581 569 for (ig=0; ig<num_gauss; ig++){ 582 570 583 /*Pick up the gaussian point: */ 584 gauss_weight=*(gauss_weights+ig); 585 gauss_coord[0]=*(first_gauss_area_coord+ig); 586 gauss_coord[1]=*(second_gauss_area_coord+ig); 587 gauss_coord[2]=*(third_gauss_area_coord+ig); 588 gauss_coord[3]=-1.0; //take the base 589 590 /*Compute strain rate viscosity and pressure: */ 591 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); 592 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); 593 pressure_input->GetParameterValue(&pressure, &gauss_coord[0]); 594 595 /*Compute Stress*/ 596 sigma_xx=2*viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure 597 sigma_yy=2*viscosity*epsilon[1]-pressure*stokesreconditioning; 598 sigma_zz=2*viscosity*epsilon[2]-pressure*stokesreconditioning; 599 sigma_xy=2*viscosity*epsilon[3]; 600 sigma_xz=2*viscosity*epsilon[4]; 601 sigma_yz=2*viscosity*epsilon[5]; 602 603 /*Get normal vector to the bed */ 604 BedNormal(&bed_normal[0],xyz_list_tria); 605 606 /*basalforce*/ 607 basalforce[0] += sigma_xx*bed_normal[0] + sigma_xy*bed_normal[1] + sigma_xz*bed_normal[2]; 608 basalforce[1] += sigma_xy*bed_normal[0] + sigma_yy*bed_normal[1] + sigma_yz*bed_normal[2]; 609 basalforce[2] += sigma_xz*bed_normal[0] + sigma_yz*bed_normal[1] + sigma_zz*bed_normal[2]; 610 611 /*Get the Jacobian determinant */ 612 tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0],gauss_coord); 613 value+=sigma_zz*Jdet2d*gauss_weight; 614 surface+=Jdet2d*gauss_weight; 571 gauss->GaussPoint(ig); 572 573 /*Compute strain rate viscosity and pressure: */ 574 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 575 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); 576 pressure_input->GetParameterValue(&pressure, &gauss_coord[0]); 577 578 /*Compute Stress*/ 579 sigma_xx=2*viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure 580 sigma_yy=2*viscosity*epsilon[1]-pressure*stokesreconditioning; 581 sigma_zz=2*viscosity*epsilon[2]-pressure*stokesreconditioning; 582 sigma_xy=2*viscosity*epsilon[3]; 583 sigma_xz=2*viscosity*epsilon[4]; 584 sigma_yz=2*viscosity*epsilon[5]; 585 586 /*Get normal vector to the bed */ 587 BedNormal(&bed_normal[0],xyz_list_tria); 588 589 /*basalforce*/ 590 basalforce[0] += sigma_xx*bed_normal[0] + sigma_xy*bed_normal[1] + sigma_xz*bed_normal[2]; 591 basalforce[1] += sigma_xy*bed_normal[0] + sigma_yy*bed_normal[1] + sigma_yz*bed_normal[2]; 592 basalforce[2] += sigma_xz*bed_normal[0] + sigma_yz*bed_normal[1] + sigma_zz*bed_normal[2]; 593 594 /*Get the Jacobian determinant */ 595 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); 596 value+=sigma_zz*Jdet2d*gauss_weight; 597 surface+=Jdet2d*gauss_weight; 615 598 } 616 599 value=value/surface; … … 2735 2718 const int numdof2d=NUMVERTICES2D*NDOF4; 2736 2719 2737 /*Collapsed formulation: */2738 Tria* tria=NULL;2739 2740 2720 /*Grid data: */ 2741 2721 double xyz_list[NUMVERTICES][3]; … … 2761 2741 int ig; 2762 2742 GaussPenta *gauss=NULL; 2763 GaussTria *gauss_tria=NULL;2764 2765 2743 2766 2744 /* specific gaussian point: */ … … 2786 2764 bool shelf; 2787 2765 2788 /*inputs: */2789 Input* vx_input=NULL;2790 Input* vy_input=NULL;2791 Input* vz_input=NULL;2792 2793 2766 /*retrive parameters: */ 2794 2767 parameters->FindParam(&analysis_type,AnalysisTypeEnum); … … 2808 2781 GetDofList(&doflist,StokesApproximationEnum); 2809 2782 2810 /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore2811 get tria gaussian points as well as segment gaussian points. For tria gaussian2812 points, order of integration is 2, because we need to integrate the product tB*D*B'2813 which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian2814 points, same deal, which yields 3 gaussian points.*/2815 2816 2783 /*Retrieve all inputs we will be needing: */ 2817 vx_input=inputs->GetInput(VxEnum);2818 vy_input=inputs->GetInput(VyEnum);2819 vz_input=inputs->GetInput(VzEnum);2784 Input* vx_input=inputs->GetInput(VxEnum); 2785 Input* vy_input=inputs->GetInput(VyEnum); 2786 Input* vz_input=inputs->GetInput(VzEnum); 2820 2787 2821 2788 /* Start looping on the number of gaussian points: */ … … 2853 2820 for(i=0;i<27;i++) for(j=0;j<27;j++) Ke_temp[i][j]+=Ke_gaussian[i][j]; 2854 2821 } 2822 delete gauss; //gauss of previous loop 2855 2823 2856 2824 if((onbed==1) && (shelf==0)){ … … 2866 2834 2867 2835 /* Start looping on the number of gaussian points: */ 2868 delete gauss; //gauss of previous loop 2869 gauss=new GaussPenta(); 2870 gauss->GaussFaceTria(0,1,2,2); 2871 gauss_tria=new GaussTria(); 2836 gauss=new GaussPenta(0,1,2,2); 2872 2837 for (ig=gauss->begin();ig<gauss->end();ig++){ 2873 2838 2874 2839 gauss->GaussPoint(ig); 2875 gauss->SynchronizeGaussTria(gauss_tria);2876 2840 2877 2841 /*Get the Jacobian determinant */ 2878 tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_tria);2842 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); 2879 2843 2880 2844 /*Get L matrix if viscous basal drag present: */ … … 2920 2884 /*Free ressources:*/ 2921 2885 delete friction; 2886 delete gauss; 2922 2887 2923 2888 } //if ( (onbed==1) && (shelf==0)) … … 2930 2895 2931 2896 /*Free ressources:*/ 2932 delete gauss;2933 delete gauss_tria;2934 2897 xfree((void**)&doflist); 2935 2936 return;2937 2898 } 2938 2899 /*}}}*/ … … 3742 3703 int ig; 3743 3704 GaussPenta *gauss=NULL; 3744 GaussTria *gauss_tria=NULL;3745 3705 3746 3706 double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ … … 3854 3814 for(i=0;i<27;i++) for(j=0;j<3;j++) Ke_temp[i][j]+=Ke_gaussian[i][j]; 3855 3815 } 3816 delete gauss; 3856 3817 3857 3818 /*Deal with 2d friction at the bedrock interface: */ … … 3865 3826 3866 3827 /* Start looping on the number of gauss 2d (nodes on the bedrock) */ 3867 delete gauss; //gauss of previous loop 3868 gauss=new GaussPenta(); 3869 gauss->GaussFaceTria(0,1,2,2); 3870 gauss_tria=new GaussTria(); 3828 gauss=new GaussPenta(0,1,2,2) 3871 3829 for (ig=gauss->begin();ig<gauss->end();ig++){ 3872 3830 3873 3831 gauss->GaussPoint(ig); 3874 gauss->SynchronizeGaussTria(gauss_tria);3875 3832 3876 3833 /*Get the Jacobian determinant */ 3877 tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_tria);3834 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); 3878 3835 3879 3836 /* Get bed at gaussian point */ … … 3895 3852 } 3896 3853 } 3854 /*Free ressources:*/ 3855 delete gauss; 3897 3856 } //if ( (onbed==1) && (shelf==1)) 3898 3857 … … 3905 3864 /*Free ressources:*/ 3906 3865 delete gauss; 3907 delete gauss_tria;3908 3866 xfree((void**)&doflist); 3909 3867 -
issm/trunk/src/c/objects/Elements/PentaRef.cpp
r5677 r5722 1677 1677 } 1678 1678 /*}}}*/ 1679 /*FUNCTION PentaRef::GetTriaJacobianDeterminant{{{1*/ 1680 void PentaRef::GetTriaJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss){ 1681 /*The Jacobian determinant is constant over the element, discard the gaussian points. 1682 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 1683 1684 double x1,x2,x3,y1,y2,y3,z1,z2,z3; 1685 1686 x1=*(xyz_list+3*0+0); 1687 y1=*(xyz_list+3*0+1); 1688 z1=*(xyz_list+3*0+2); 1689 x2=*(xyz_list+3*1+0); 1690 y2=*(xyz_list+3*1+1); 1691 z2=*(xyz_list+3*1+2); 1692 x3=*(xyz_list+3*2+0); 1693 y3=*(xyz_list+3*2+1); 1694 z3=*(xyz_list+3*2+2); 1695 1696 *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5); 1697 if(*Jdet<0) ISSMERROR("negative jacobian determinant!"); 1698 1699 } 1700 /*}}}*/ 1679 1701 /*FUNCTION PentaRef::GetJacobianInvert {{{1*/ 1680 1702 void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss){ -
issm/trunk/src/c/objects/Elements/PentaRef.h
r5677 r5722 58 58 void GetJacobian(double* J, double* xyz_list,GaussPenta* gauss); 59 59 void GetJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss); 60 void GetTriaJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss); 60 61 void GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss); 61 62 void GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss); -
issm/trunk/src/c/objects/Gauss/GaussPenta.cpp
r5675 r5722 80 80 } 81 81 /*}}}*/ 82 /*FUNCTION GaussPenta::GaussPenta{{{1*/ 83 GaussPenta::GaussPenta(int index1, int index2, int index3, int order){ 84 85 /*Basal Tria*/ 86 if(index1==0 && index2==1 && index3==2){ 87 88 /*Get GaussTria*/ 89 GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order); 90 91 /*compute z coordinate*/ 92 coords4=(double*)xmalloc(numgauss*sizeof(double)); 93 for(int i=0;i<numgauss;i++) coords4[i]=-1.0; 94 } 95 else{ 96 ISSMERROR("Tria not supported yet"); 97 } 98 99 } 100 /*}}}*/ 82 101 /*FUNCTION GaussPenta::~GaussPenta(){{{1*/ 83 102 GaussPenta::~GaussPenta(){ -
issm/trunk/src/c/objects/Gauss/GaussPenta.h
r5675 r5722 34 34 GaussPenta(); 35 35 GaussPenta(int order_horiz,int order_vert); 36 GaussPenta(int index1, int index2, int index3, int order); 36 37 ~GaussPenta(); 37 38 -
issm/trunk/src/c/objects/Loads/Icefront.cpp
r5720 r5722 1354 1354 } 1355 1355 /*}}}*/ 1356 /*FUNCTION Icefront::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in) {{{1*/1357 void Icefront::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in){1358 1359 /*output: */1360 double value;1361 1362 /*Get value on Element 1*/1363 element->GetParameterValue(&value,nodes[0],nodes[1],gauss_coord,input_in);1364 1365 /*Assign output pointers:*/1366 *pvalue=value;1367 }1368 /*}}}*/1369 1356 /*FUNCTION Icefront::GetNormal {{{1*/ 1370 1357 void Icefront:: GetNormal(double* normal,double xyz_list[2][3]){ -
issm/trunk/src/c/objects/Loads/Icefront.h
r5715 r5722 85 85 void QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, 86 86 double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list); 87 void GetParameterValue(double* pvalue, double gauss_coord,Input* input_in);88 87 void GetNormal(double* normal,double xyz_list[2][3]); 89 88 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.