Changeset 10480
- Timestamp:
- 11/04/11 16:41:35 (13 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r10444 r10480 5986 5986 ElementMatrix* Ke1=CreateKMatrixDiagnosticStokesViscous(); 5987 5987 ElementMatrix* Ke2=CreateKMatrixDiagnosticStokesFriction(); 5988 //ElementMatrix* Ke2=CreateKMatrixDiagnosticStokesFriction2(); 5988 5989 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); 5989 5990 … … 6136 6137 6137 6138 /*Transform Coordinate System*/ 6139 if(id==1) printarray(&bed_normal[0],1,3); 6140 if(id==1) printarray(&nodes[0]->coord_system[0][0],3,3); 6141 if(id==1) printf("Before rotation\n"); 6142 if(id==1) printarray(Ke->values,4,Ke->ncols); 6138 6143 TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4); 6144 if(id==1) printf("After rotation\n"); 6145 if(id==1) printarray(Ke->values,4,Ke->ncols); 6146 if(id==1) printf("New:\n"); 6147 if(id==1) Ke=CreateKMatrixDiagnosticStokesFriction2(); 6148 6149 //if(id==56) printarray(Ke->values,1,9); 6150 //if(id==56) printarray(&nodes[0]->coord_system[0][0],3,3); 6151 //if(id==56) _error_("STOP"); 6139 6152 6140 6153 /*Clean up and return*/ 6141 6154 delete gauss; 6142 6155 delete friction; 6156 return Ke; 6157 } 6158 /*}}}*/ 6159 /*FUNCTION Penta::CreateKMatrixDiagnosticStokesFriction2{{{1*/ 6160 ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction2(void){ 6161 6162 /*Constants*/ 6163 const int numdof=NUMVERTICES*NDOF4; 6164 const int numdof2d=NUMVERTICES2D*NDOF4; 6165 6166 /*Intermediaries */ 6167 int i,j,ig; 6168 int analysis_type,approximation; 6169 double alpha2,Jdet2d; 6170 double stokesreconditioning,viscosity; 6171 double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ 6172 double xyz_list[NUMVERTICES][3]; 6173 double xyz_list_tria[NUMVERTICES2D][3]; 6174 double LStokes2[2][numdof2d]; 6175 //double LprimeStokes2[4][numdof2d]; 6176 double DLStokes[2][2]={0.0}; 6177 double Ke_drag_gaussian[numdof2d][numdof2d]; 6178 Friction* friction=NULL; 6179 GaussPenta *gauss=NULL; 6180 6181 /*If on water or not Stokes, skip stiffness: */ 6182 inputs->GetInputValue(&approximation,ApproximationEnum); 6183 if(IsFloating() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL; 6184 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); 6185 6186 /*Retrieve all inputs and parameters*/ 6187 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 6188 parameters->FindParam(&analysis_type,AnalysisTypeEnum); 6189 parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum); 6190 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); 6191 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); 6192 Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); 6193 for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j]; 6194 6195 /*build friction object, used later on: */ 6196 friction=new Friction("3d",inputs,matpar,analysis_type); 6197 6198 /* Start looping on the number of gaussian points: */ 6199 gauss=new GaussPenta(0,1,2,2); 6200 for (ig=gauss->begin();ig<gauss->end();ig++){ 6201 6202 gauss->GaussPoint(ig); 6203 6204 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); 6205 GetLStokes3(&LStokes2[0][0], gauss); 6206 //GetLprimeStokes2(&LprimeStokes2[0][0], &xyz_list[0][0], gauss); 6207 6208 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 6209 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); 6210 6211 friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); 6212 6213 DLStokes[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx 6214 DLStokes[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy 6215 //DLStokes[2][2] = +2*viscosity*gauss->weight*Jdet2d; 6216 //DLStokes[3][3] = -stokesreconditioning*gauss->weight*Jdet2d; 6217 6218 //TripleMultiply( &LStokes2[0][0],4,numdof2d,1, 6219 // &DLStokes[0][0],4,4,0, 6220 // &LprimeStokes2[0][0],4,numdof2d,0, 6221 // &Ke_drag_gaussian[0][0],0); 6222 TripleMultiply( &LStokes2[0][0],2,numdof2d,1, 6223 &DLStokes[0][0],2,2,0, 6224 &LStokes2[0][0],2,numdof2d,0, 6225 &Ke_drag_gaussian[0][0],0); 6226 6227 for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof+j]+=Ke_drag_gaussian[i][j]; 6228 } 6229 6230 /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ 6231 //TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4); 6232 6233 //printf("New:\n"); 6234 //printarray(&Ke_drag_gaussian[0][0],numdof2d,numdof2d); 6235 //printf("Old:\n"); 6236 //Ke=CreateKMatrixDiagnosticStokesFriction(); 6237 //_error_("STOP"); 6238 //if(id==1) printarray(Ke->values,4,Ke->ncols); 6239 //if(id==1) printarray(Ke->values,1,9); 6240 6241 /*Clean up and return*/ 6242 delete gauss; 6243 delete friction; 6244 6245 //return CreateKMatrixDiagnosticStokesFriction(); 6143 6246 return Ke; 6144 6247 } -
issm/trunk/src/c/objects/Elements/Penta.h
r10440 r10480 223 223 ElementMatrix* CreateKMatrixDiagnosticStokesViscous(void); 224 224 ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void); 225 ElementMatrix* CreateKMatrixDiagnosticStokesFriction2(void); 225 226 ElementMatrix* CreateKMatrixDiagnosticVert(void); 226 227 ElementMatrix* CreateKMatrixDiagnosticVertVolume(void); -
issm/trunk/src/c/objects/Elements/PentaRef.cpp
r10135 r10480 651 651 } 652 652 /*}}}*/ 653 /*FUNCTION PentaRef::GetLStokes2{{{1*/ 654 void PentaRef::GetLStokes2(double* LStokes, GaussPenta* gauss){ 655 /* 656 * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. 657 * For node i, Li can be expressed in the actual coordinate system 658 * by: 659 * Li=[ h 0 0 0] 660 * [ 0 h 0 0] 661 * [ 0 0 h h] 662 * [ 0 0 0 0] 663 * where h is the interpolation function for node i. 664 */ 665 666 int i; 667 int num_dof=4; 668 669 double l1l2l3[NUMNODESP1_2d]; 670 671 672 /*Get l1l2l3 in actual coordinate system: */ 673 l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; 674 l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; 675 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; 676 677 /*Build LStokes: */ 678 for (i=0;i<3;i++){ 679 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i]; 680 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.; 681 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.; 682 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.; 683 684 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.; 685 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i]; 686 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.; 687 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.; 688 689 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+0)=0.; 690 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0.; 691 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i]; 692 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0.; 693 694 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+0)=0.; 695 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0.; 696 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i]; 697 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0.; 698 } 699 } 700 /*}}}*/ 701 /*FUNCTION PentaRef::GetLStokes3{{{1*/ 702 void PentaRef::GetLStokes3(double* LStokes, GaussPenta* gauss){ 703 /* 704 * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. 705 * For node i, Li can be expressed in the actual coordinate system 706 * by: 707 * Li=[ h 0 ] 708 * [ 0 h ] 709 * [ 0 0 ] 710 * [ 0 0 ] 711 * where h is the interpolation function for node i. 712 */ 713 714 int i; 715 int num_dof=4; 716 717 double l1l2l3[NUMNODESP1_2d]; 718 719 720 /*Get l1l2l3 in actual coordinate system: */ 721 l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; 722 l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; 723 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; 724 725 /*Build LStokes: */ 726 for (i=0;i<3;i++){ 727 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i]; 728 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.; 729 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.; 730 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.; 731 732 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.; 733 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i]; 734 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.; 735 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.; 736 } 737 } 738 /*}}}*/ 739 /*FUNCTION PentaRef::GetLprimeStokes2{{{1*/ 740 void PentaRef::GetLprimeStokes2(double* LStokes,double* xyz_list,GaussPenta* gauss){ 741 /* 742 * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. 743 * For node i, Li can be expressed in the actual coordinate system 744 * by: 745 * Li=[ h 0 0 0] 746 * [ 0 h 0 0] 747 * [ 0 0 dh\dz 0] 748 * [ 0 0 0 h] 749 * where h is the interpolation function for node i. 750 */ 751 752 int i; 753 int num_dof=4; 754 755 double l1l2l3[NUMNODESP1_2d]; 756 double dh1dh6[3][NUMNODESP1]; 757 758 759 /*Get l1l2l3 in actual coordinate system: */ 760 l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; 761 l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; 762 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; 763 764 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss); 765 766 /*Build LStokes: */ 767 for (i=0;i<3;i++){ 768 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i]; 769 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.; 770 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.; 771 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.; 772 773 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.; 774 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i]; 775 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.; 776 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.; 777 778 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+0)=0.; 779 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0.; 780 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=dh1dh6[2][i]; 781 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0.; 782 783 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+0)=0.; 784 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0.; 785 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=0.; 786 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=l1l2l3[i]; 787 } 788 } 789 /*}}}*/ 653 790 /*FUNCTION PentaRef::GetLprimeStokes {{{1*/ 654 791 void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ … … 658 795 * For node i, Lpi can be expressed in the actual coordinate system 659 796 * by: 660 * Lpi=[ h 0 0 0] 661 * [ 0 h 0 0] 662 * [ h 0 0 0] 663 * [ 0 h 0 0] 664 * [ 0 0 h 0] 665 * [ 0 0 h 0] 666 * [ 0 0 dh/dz 0] 667 * [ 0 0 dh/dz 0] 668 * [ 0 0 dh/dz 0] 669 * [dh/dz 0 dh/dx 0] 670 * [ 0 dh/dz dh/dy 0] 671 * [ 0 0 0 h] 672 * [ 0 0 0 h] 673 * [ 0 0 0 h] 797 * Lpi=[ h 0 0 0]1 798 * [ 0 h 0 0]2 799 * [ h 0 0 0]3 800 * [ 0 h 0 0]4 801 * [ 0 0 h 0]5 802 * [ 0 0 h 0]6 803 * [ 0 0 dh/dz 0]7 804 * [ 0 0 dh/dz 0]8 805 * [ 0 0 dh/dz 0]9 806 * [dh/dz 0 dh/dx 0]0 807 * [ 0 dh/dz dh/dy 0]1 808 * [ 0 0 0 h]2 809 * [ 0 0 0 h]3 810 * [ 0 0 0 h]4 811 * 812 * Li=[ h 0 0 0]1 813 * [ 0 h 0 0]2 814 * [ 0 0 h 0]3 815 * [ 0 0 h 0]4 816 * [ h 0 0 0]5 817 * [ 0 h 0 0]6 818 * [ h 0 0 0]7 819 * [ 0 h 0 0]8 820 * [ 0 0 h 0]9 821 * [ 0 0 h 0]0 822 * [ 0 0 h 0]1 823 * [ h 0 0 0]2 824 * [ 0 h 0 0]3 825 * [ 0 0 h 0]4 674 826 * where h is the interpolation function for node i. 675 827 */ -
issm/trunk/src/c/objects/Elements/PentaRef.h
r10135 r10480 51 51 void GetL(double* L, GaussPenta* gauss,int numdof); 52 52 void GetLStokes(double* LStokes, GaussPenta* gauss); 53 void GetLStokes2(double* LStokes, GaussPenta* gauss); 54 void GetLStokes3(double* LStokes, GaussPenta* gauss); 55 void GetLprimeStokes2(double* LStokes,double* xyz_list,GaussPenta* gauss); 53 56 void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss); 54 57 void GetLMacAyealStokes(double* LMacAyealStokes, GaussPenta* gauss); -
issm/trunk/src/c/shared/Elements/elements.h
r10407 r10480 28 28 for(int i=0;i<lines;i++){ 29 29 printf(" [ "); 30 for(int j=0;j<cols;j++) printf(" % 7.5g ",array[i*cols+j]);30 for(int j=0;j<cols;j++) printf(" %12.7g ",array[i*cols+j]); 31 31 printf(" ]\n"); 32 32 }
Note:
See TracChangeset
for help on using the changeset viewer.