![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 22874
- Timestamp:
- 06/25/18 07:32:30 (7 years ago)
- Location:
- issm/trunk-jpl/src/c/analyses
- Files:
-
- 2 edited
-
StressbalanceAnalysis.cpp (modified) (57 diffs)
-
ThermalAnalysis.cpp (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
r22573 r22874 611 611 break; 612 612 case L1L2ApproximationEnum: numdofs =2; break; 613 case HOApproximationEnum: 613 case HOApproximationEnum: 614 614 switch(domaintype){ 615 615 case Domain3DEnum: numdofs=2; break; … … 785 785 iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); 786 786 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); 787 if(FrictionCoupling==1){ 788 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); 787 if(FrictionCoupling==3){ 788 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} 789 else if(FrictionCoupling==4){ 790 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); 789 791 } 790 792 break; … … 798 800 iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum); 799 801 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); 800 if(FrictionCoupling==1){ 801 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); 802 if(FrictionCoupling==3){ 803 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} 804 else if(FrictionCoupling==4){ 805 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); 802 806 } 803 807 break; … … 827 831 iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); 828 832 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); 829 if(FrictionCoupling==1){ 830 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); 833 if(FrictionCoupling==3){ 834 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} 835 else if(FrictionCoupling==4){ 836 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); 837 831 838 } 832 839 break; … … 934 941 solutionsequence_newton(femmodel); 935 942 else 936 solutionsequence_nonlinear(femmodel,conserve_loads); 937 } 938 else if(!isFS && (isSSA || isHO || isL1L2)){ 943 solutionsequence_nonlinear(femmodel,conserve_loads); 944 } 945 else if(!isFS && (isSSA || isHO || isL1L2)){ 939 946 if(VerboseSolution()) _printf0_(" computing velocities\n"); 940 947 … … 943 950 solutionsequence_newton(femmodel); 944 951 else 945 solutionsequence_nonlinear(femmodel,conserve_loads); 952 solutionsequence_nonlinear(femmodel,conserve_loads); 946 953 947 954 if(domaintype==Domain2DverticalEnum && isSSA){ … … 966 973 element->GetInputValue(&approximation,ApproximationEnum); 967 974 switch(approximation){ 968 case FSApproximationEnum: 975 case FSApproximationEnum: 969 976 return CreateDVectorFS(element); 970 977 default: … … 979 986 element->GetInputValue(&approximation,ApproximationEnum); 980 987 switch(approximation){ 981 case SSAApproximationEnum: 988 case SSAApproximationEnum: 982 989 return CreateJacobianMatrixSSA(element); 983 case HOApproximationEnum: 990 case HOApproximationEnum: 984 991 return CreateJacobianMatrixHO(element); 985 case FSApproximationEnum: 992 case FSApproximationEnum: 986 993 return CreateJacobianMatrixFS(element); 987 994 case NoneApproximationEnum: … … 997 1004 case SIAApproximationEnum: 998 1005 return NULL; 999 case SSAApproximationEnum: 1006 case SSAApproximationEnum: 1000 1007 return CreateKMatrixSSA(element); 1001 case L1L2ApproximationEnum: 1008 case L1L2ApproximationEnum: 1002 1009 return CreateKMatrixL1L2(element); 1003 case HOApproximationEnum: 1010 case HOApproximationEnum: 1004 1011 return CreateKMatrixHO(element); 1005 case FSApproximationEnum: 1012 case FSApproximationEnum: 1006 1013 return CreateKMatrixFS(element); 1007 case SSAHOApproximationEnum: 1014 case SSAHOApproximationEnum: 1008 1015 return CreateKMatrixSSAHO(element); 1009 case HOFSApproximationEnum: 1016 case HOFSApproximationEnum: 1010 1017 return CreateKMatrixHOFS(element); 1011 case SSAFSApproximationEnum: 1018 case SSAFSApproximationEnum: 1012 1019 return CreateKMatrixSSAFS(element); 1013 1020 case NoneApproximationEnum: … … 1024 1031 case SIAApproximationEnum: 1025 1032 return NULL; 1026 case SSAApproximationEnum: 1033 case SSAApproximationEnum: 1027 1034 return CreatePVectorSSA(element); 1028 case L1L2ApproximationEnum: 1035 case L1L2ApproximationEnum: 1029 1036 return CreatePVectorL1L2(element); 1030 case HOApproximationEnum: 1037 case HOApproximationEnum: 1031 1038 return CreatePVectorHO(element); 1032 case FSApproximationEnum: 1039 case FSApproximationEnum: 1033 1040 return CreatePVectorFS(element); 1034 case SSAHOApproximationEnum: 1041 case SSAHOApproximationEnum: 1035 1042 return CreatePVectorSSAHO(element); 1036 case HOFSApproximationEnum: 1043 case HOFSApproximationEnum: 1037 1044 return CreatePVectorHOFS(element); 1038 case SSAFSApproximationEnum: 1045 case SSAFSApproximationEnum: 1039 1046 return CreatePVectorSSAFS(element); 1040 1047 case NoneApproximationEnum: … … 1126 1133 InputUpdateFromSolutionFS(solution,element); 1127 1134 return; 1128 case SIAApproximationEnum: 1135 case SIAApproximationEnum: 1129 1136 return; 1130 case SSAApproximationEnum: 1137 case SSAApproximationEnum: 1131 1138 InputUpdateFromSolutionSSA(solution,element); 1132 1139 return; 1133 case HOApproximationEnum: 1140 case HOApproximationEnum: 1134 1141 InputUpdateFromSolutionHO(solution,element); 1135 1142 return; … … 1581 1588 /*Retrieve all inputs and parameters*/ 1582 1589 element->GetVerticesCoordinates(&xyz_list); 1583 Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); 1590 Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); 1584 1591 Input* surface_input =element->GetInput(SurfaceEnum); _assert_(surface_input); 1585 1592 IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); … … 1689 1696 }/*}}}*/ 1690 1697 void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 1691 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 1698 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 1692 1699 * For node i, Bi can be expressed in the actual coordinate system 1693 * by: 1700 * by: 1694 1701 * 2D 1D 1695 1702 * Bi=[ dN/dx 0 ] Bi=[ dN/dx ] 1696 * [ 0 dN/dy ] 1697 * [ 1/2*dN/dy 1/2*dN/dx ] 1703 * [ 0 dN/dy ] 1704 * [ 1/2*dN/dy 1/2*dN/dx ] 1698 1705 * where N is the finiteelement function for node i. 1699 1706 * … … 1729 1736 }/*}}}*/ 1730 1737 void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 1731 /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. 1738 /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. 1732 1739 * For node i, Bi can be expressed in the actual coordinate system 1733 * by: 1740 * by: 1734 1741 * 2D 1D 1735 1742 * Bi=[ N 0 ] Bi = N … … 1766 1773 }/*}}}*/ 1767 1774 void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 1768 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 1775 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 1769 1776 * For node i, Bi' can be expressed in the actual coordinate system 1770 * by: 1777 * by: 1771 1778 * 2D 1D 1772 1779 * Bi_prime=[ 2*dN/dx dN/dy ] Bi_prime=[ 2*dN/dx ] … … 2050 2057 /*Retrieve all inputs and parameters*/ 2051 2058 element->GetVerticesCoordinates(&xyz_list); 2052 Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); 2059 Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); 2053 2060 Input* surface_input =element->GetInput(SurfaceEnum); _assert_(surface_input); 2054 2061 IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); … … 2677 2684 }/*}}}*/ 2678 2685 void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 2679 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 2686 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 2680 2687 * For node i, Bi can be expressed in the actual coordinate system 2681 * by: 2688 * by: 2682 2689 * 3D 2D 2683 2690 * 2684 2691 * Bi=[ dh/dx 0 ] Bi=[ dh/dx] 2685 2692 * [ 0 dh/dy ] [ dh/dy] 2686 * [ 1/2*dh/dy 1/2*dh/dx ] 2687 * [ 1/2*dh/dz 0 ] 2688 * [ 0 1/2*dh/dz ] 2693 * [ 1/2*dh/dy 1/2*dh/dx ] 2694 * [ 1/2*dh/dz 0 ] 2695 * [ 0 1/2*dh/dz ] 2689 2696 * where h is the interpolation function for node i. 2690 2697 * … … 2726 2733 }/*}}}*/ 2727 2734 void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 2728 /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. 2735 /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. 2729 2736 * For node i, Bi can be expressed in the actual coordinate system 2730 * by: 2737 * by: 2731 2738 * 3D 2D 2732 2739 * Bi=[ N 0 ] Bi=N … … 2763 2770 }/*}}}*/ 2764 2771 void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 2765 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 2772 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 2766 2773 * For node i, Bi' can be expressed in the actual coordinate system 2767 * by: 2774 * by: 2768 2775 * 3D 2D 2769 2776 * Bi_prime=[ 2*dN/dx dN/dy ] Bi_prime=[ 2*dN/dx ] 2770 2777 * [ dN/dx 2*dN/dy ] [ dN/dy ] 2771 * [ dN/dy dN/dx ] 2778 * [ dN/dy dN/dx ] 2772 2779 * where hNis the finiteelement function for node i. 2773 2780 * … … 3084 3091 if(dim==2) slope2=slope[0]*slope[0]; 3085 3092 else if(dim==3) slope2=slope[0]*slope[0]+slope[1]*slope[1]; 3086 scalar = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt; 3093 scalar = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt; 3087 3094 for(i=0;i<vnumnodes;i++){ 3088 3095 for(j=0;j<vnumnodes;j++){ … … 3286 3293 &Ke->values[0],1); 3287 3294 3288 if(element->IsOnBase() && 0){ 3295 if(element->IsOnBase() && 0){ 3289 3296 element->FindParam(&rl,AugmentedLagrangianRlambdaEnum); 3290 3297 element->GetVerticesCoordinatesBase(&xyz_list_base); … … 4085 4092 gauss->GaussPoint(ig); 4086 4093 element->JacobianDeterminant(&Jdet,xyz_list,gauss); 4087 4094 4088 4095 pressure_input->GetInputValue(&pressure, gauss); 4089 4096 element->NodalFunctionsDerivativesVelocity(dbasis,xyz_list,gauss); … … 4096 4103 } 4097 4104 4098 if(element->IsOnBase() && 0){ 4105 if(element->IsOnBase() && 0){ 4099 4106 IssmDouble sigmann; 4100 4107 IssmDouble* vbasis = xNew<IssmDouble>(numnodes); … … 4305 4312 }/*}}}*/ 4306 4313 void StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 4307 /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. 4314 /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. 4308 4315 * For node i, Bvi can be expressed in the actual coordinate system 4309 4316 * by: Bvi=[ dphi/dx 0 ] … … 4417 4424 }/*}}}*/ 4418 4425 void StressbalanceAnalysis::GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 4419 /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. 4426 /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. 4420 4427 * For node i, Li can be expressed in the actual coordinate system 4421 * by in 3d 4428 * by in 3d 4422 4429 * Li=[ h 0 0 0 ] 4423 4430 * [ 0 h 0 0 ] … … 4477 4484 }/*}}}*/ 4478 4485 void StressbalanceAnalysis::GetBFSprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 4479 /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 4486 /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 4480 4487 * For node i, Bi' can be expressed in the actual coordinate system 4481 * by: 4488 * by: 4482 4489 * Bvi' = [ dphi/dx 0 ] 4483 4490 * [ 0 dphi/dy ] … … 4590 4597 }/*}}}*/ 4591 4598 void StressbalanceAnalysis::GetBFSprimeUzawa(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 4592 /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2. 4599 /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2. 4593 4600 * For node i, Bi' can be expressed in the actual coordinate system 4594 * by: 4601 * by: 4595 4602 * Bvi' = [ dphi/dx dphi/dy ] 4596 4603 * … … 4626 4633 }/*}}}*/ 4627 4634 void StressbalanceAnalysis::GetBFSprimevel(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 4628 /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 4635 /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 4629 4636 * For node i, Bi' can be expressed in the actual coordinate system 4630 * by: 4637 * by: 4631 4638 * Bvi' = [ dphi/dx 0 ] 4632 4639 * [ 0 dphi/dy ] … … 4689 4696 }/*}}}*/ 4690 4697 void StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 4691 /*Compute B matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi. 4698 /*Compute B matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi. 4692 4699 */ 4693 4700 … … 4711 4718 }/*}}}*/ 4712 4719 void StressbalanceAnalysis::GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 4713 /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. 4720 /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. 4714 4721 * For node i, Bvi can be expressed in the actual coordinate system 4715 4722 * by: Bvi=[ dphi/dx 0 ] … … 4776 4783 void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 4777 4784 /*Compute C matrix. C=[Cp1 Cp2 ...] where: 4778 * Cpi=[phi phi]. 4785 * Cpi=[phi phi]. 4779 4786 */ 4780 4787 … … 4796 4803 }/*}}}*/ 4797 4804 void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 4798 /* Compute C' matrix. C'=[C1' C2' ...] 4805 /* Compute C' matrix. C'=[C1' C2' ...] 4799 4806 * Ci' = [ phi 0 ] 4800 4807 * [ 0 phi ] … … 4930 4937 int tnumnodes = element->GetNumberOfVertices(); //Tensors, P1 DG 4931 4938 IssmDouble* epsxx = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapxx = xNew<IssmDouble>(tnumnodes); 4932 IssmDouble* epsyy = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapyy = xNew<IssmDouble>(tnumnodes); 4933 IssmDouble* epsxy = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapxy = xNew<IssmDouble>(tnumnodes); 4934 IssmDouble* epszz = NULL; IssmDouble* sigmapzz = NULL; 4935 IssmDouble* epsxz = NULL; IssmDouble* sigmapxz = NULL; 4936 IssmDouble* epsyz = NULL; IssmDouble* sigmapyz = NULL; 4939 IssmDouble* epsyy = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapyy = xNew<IssmDouble>(tnumnodes); 4940 IssmDouble* epsxy = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapxy = xNew<IssmDouble>(tnumnodes); 4941 IssmDouble* epszz = NULL; IssmDouble* sigmapzz = NULL; 4942 IssmDouble* epsxz = NULL; IssmDouble* sigmapxz = NULL; 4943 IssmDouble* epsyz = NULL; IssmDouble* sigmapyz = NULL; 4937 4944 if(dim==3){ 4938 4945 epszz = xNew<IssmDouble>(tnumnodes); sigmapzz = xNew<IssmDouble>(tnumnodes); 4939 epsxz = xNew<IssmDouble>(tnumnodes); sigmapxz = xNew<IssmDouble>(tnumnodes); 4940 epsyz = xNew<IssmDouble>(tnumnodes); sigmapyz = xNew<IssmDouble>(tnumnodes); 4946 epsxz = xNew<IssmDouble>(tnumnodes); sigmapxz = xNew<IssmDouble>(tnumnodes); 4947 epsyz = xNew<IssmDouble>(tnumnodes); sigmapyz = xNew<IssmDouble>(tnumnodes); 4941 4948 } 4942 4949 … … 5172 5179 epsxx = dvx[0]; 5173 5180 epsyy = dvy[1]; 5174 epsxy = 0.5*(dvx[1] + dvy[0]); 5181 epsxy = 0.5*(dvx[1] + dvy[0]); 5175 5182 if(dim==3){ 5176 epszz = dvz[2]; 5183 epszz = dvz[2]; 5177 5184 epsxz = 0.5*(dvx[2] + dvz[0]); 5178 5185 epsyz = 0.5*(dvy[2] + dvz[1]); … … 5183 5190 B_input->GetInputValue(&B,gauss); 5184 5191 n_input->GetInputValue(&n,gauss); 5185 coef1 = B*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0 = 2 * B/(2* (1/sqrt(2) )^(n-1)/n ) 5192 coef1 = B*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0 = 2 * B/(2* (1/sqrt(2) )^(n-1)/n ) 5186 5193 coef2 = r; 5187 5194 if(dim==2){ … … 5207 5214 } 5208 5215 else{ 5209 dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + epszz_old*epszz_old 5216 dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + epszz_old*epszz_old 5210 5217 +2.*(epsxy_old*epsxy_old + epsxz_old*epsxz_old + epsyz_old*epsyz_old)); 5211 5218 } … … 5360 5367 epsxx = dvx[0]; 5361 5368 epsyy = dvy[1]; 5362 epsxy = 0.5*(dvx[1] + dvy[0]); 5369 epsxy = 0.5*(dvx[1] + dvy[0]); 5363 5370 if(dim==3){ 5364 epszz = dvz[2]; 5371 epszz = dvz[2]; 5365 5372 epsxz = 0.5*(dvx[2] + dvz[0]); 5366 5373 epsyz = 0.5*(dvy[2] + dvz[1]); … … 5614 5621 for(i=0;i<numdof2dm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftot+j+numdofm]+=Ke_drag[i][j]; 5615 5622 for(i=0;i<numdof2d;i++) for(j=0;j<numdof2dm;j++) Ke->values[(i+numdofm)*numdoftot+j]+=Ke_drag2[i][j]; 5616 5623 5617 5624 /*Transform Coordinate System*/ 5618 5625 element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list); … … 5647 5654 IssmDouble D[4][4]={0.0}; // material matrix, simple scalar matrix. 5648 5655 IssmDouble D2[3][3]={0.0}; // material matrix, simple scalar matrix. 5649 IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix 5650 IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix 5656 IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix 5657 IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix 5651 5658 IssmDouble *xyz_list = NULL; 5652 5659 … … 5719 5726 &Ke_gg2[0][0],1); 5720 5727 5721 } 5728 } 5722 5729 for(i=0;i<numdofs;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j]; 5723 5730 for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[i][j]; … … 5803 5810 5804 5811 DL_scalar=alpha2*gauss->weight*Jdet2d; 5805 for (i=0;i<2;i++) DL[i][i]=DL_scalar; 5812 for (i=0;i<2;i++) DL[i][i]=DL_scalar; 5806 5813 5807 5814 /* Do the triple producte tL*D*L: */ … … 5884 5891 element->JacobianDeterminant(&Jdet, xyz_list,gauss); 5885 5892 this->GetBSSAHO(B, element,xyz_list, gauss); 5886 this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); 5893 this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); 5887 5894 element->material->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input); 5888 5895 … … 5894 5901 Bprime,3,numdofm,0, 5895 5902 Ke_gg,1); 5896 } 5903 } 5897 5904 for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j]; 5898 5905 for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i]; … … 5903 5910 /*Clean-up and return*/ 5904 5911 basaltria->DeleteMaterials(); delete basaltria; 5905 5912 5906 5913 delete gauss; 5907 5914 delete gauss_tria; … … 5941 5948 Ke1->StaticCondensation(3,&indices[0]); 5942 5949 int init = element->FiniteElement(); 5943 element->SetTemporaryElementType(P1Enum); 5950 element->SetTemporaryElementType(P1Enum); 5944 5951 ElementMatrix* Ke2=CreateKMatrixSSA3d(element); 5945 element->SetTemporaryElementType(init); 5952 element->SetTemporaryElementType(init); 5946 5953 ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(element); 5947 5954 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); … … 5987 5994 if(element->IsFloating() || !element->IsOnBase()) return NULL; 5988 5995 5989 /*Build a tria element using the 3 nodes of the base of the penta. Then use 5996 /*Build a tria element using the 3 nodes of the base of the penta. Then use 5990 5997 * the tria functionality to build a friction stiffness matrix on these 3 5991 5998 * nodes: */ … … 6223 6230 element->NodalFunctionsP1(&basis[0],gauss); 6224 6231 element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss); 6225 6232 6226 6233 element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); 6227 6234 vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); … … 6469 6476 }/*}}}*/ 6470 6477 void StressbalanceAnalysis::GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 6471 /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2. 6478 /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2. 6472 6479 * For node i, Bprimei can be expressed in the actual coordinate system 6473 * by: 6480 * by: 6474 6481 * Bprimei=[ 2*dh/dx dh/dy 0 0 ] 6475 6482 * [ dh/dx 2*dh/dy 0 0 ] … … 6518 6525 }/*}}}*/ 6519 6526 void StressbalanceAnalysis::GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 6520 /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. 6527 /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. 6521 6528 * For node i, Bprimei can be expressed in the actual coordinate system 6522 * by: 6529 * by: 6523 6530 * Bprimei=[ dN/dx 0 ] 6524 6531 * [ 0 dN/dy ] … … 6553 6560 }/*}}}*/ 6554 6561 void StressbalanceAnalysis::GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 6555 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 6562 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 6556 6563 * For node i, Bi can be expressed in the actual coordinate system 6557 * by: 6564 * by: 6558 6565 * Bi=[ dh/dx 0 0 0 ] 6559 6566 * [ 0 dh/dy 0 0 ] … … 6611 6618 }/*}}}*/ 6612 6619 void StressbalanceAnalysis::GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 6613 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 6620 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 6614 6621 * For node i, Bi can be expressed in the actual coordinate system 6615 * by: 6622 * by: 6616 6623 * Bi=[ dN/dx 0 ] 6617 6624 * [ 0 dN/dy ] … … 6643 6650 }/*}}}*/ 6644 6651 void StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 6645 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2. 6652 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2. 6646 6653 * For node i, Bi can be expressed in the actual coordinate system 6647 * by: 6654 * by: 6648 6655 * Bi=[ dh/dx 0 ] 6649 6656 * [ 0 dh/dy ] … … 6674 6681 }/*}}}*/ 6675 6682 void StressbalanceAnalysis::GetLprimeFSSSA(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/ 6676 /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 6683 /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 6677 6684 * For node i, Lpi can be expressed in the actual coordinate system 6678 * by: 6685 * by: 6679 6686 * Lpi=[ h 0 ] 6680 6687 * [ 0 h ] … … 6708 6715 }/*}}}*/ 6709 6716 void StressbalanceAnalysis::GetLprimeSSAFS(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/ 6710 /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 6717 /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 6711 6718 * For node i, Lpi can be expressed in the actual coordinate system 6712 * by: 6719 * by: 6713 6720 * Lpi=[ h 0 0 0] 6714 6721 * [ 0 h 0 0] … … 6813 6820 }/*}}}*/ 6814 6821 void StressbalanceAnalysis::GetLFSSSA(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/ 6815 /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. 6822 /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. 6816 6823 * For node i, Li can be expressed in the actual coordinate system 6817 * by: 6824 * by: 6818 6825 * Li=[ h 0 0 ] 6819 6826 * [ 0 h 0 ] … … 6853 6860 void StressbalanceAnalysis::GetLSSAFS(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/ 6854 6861 /* 6855 * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. 6862 * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. 6856 6863 * For node i, Li can be expressed in the actual coordinate system 6857 * by: 6864 * by: 6858 6865 * Li=[ h 0 ] 6859 6866 * [ 0 h ] … … 7158 7165 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]); 7159 7166 7160 /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, 7167 /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, 7161 7168 *so the pressure is just the pressure at the bedrock: */ 7162 7169 rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); -
issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
r22523 r22874 12 12 int finiteelement; 13 13 iomodel->FindConstant(&finiteelement,"md.thermal.fe"); 14 _assert_(finiteelement==P1Enum); 14 _assert_(finiteelement==P1Enum); 15 15 16 16 /*Output*/ … … 94 94 int finiteelement; 95 95 iomodel->FindConstant(&finiteelement,"md.thermal.fe"); 96 _assert_(finiteelement==P1Enum); 97 96 _assert_(finiteelement==P1Enum); 97 98 98 if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); 99 99 ::CreateNodes(nodes,iomodel,ThermalAnalysisEnum,finiteelement); … … 113 113 int finiteelement; 114 114 iomodel->FindConstant(&finiteelement,"md.thermal.fe"); 115 _assert_(finiteelement==P1Enum); 115 _assert_(finiteelement==P1Enum); 116 116 int counter=0; 117 117 for(int i=0;i<iomodel->numberofelements;i++){ … … 189 189 iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); 190 190 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); 191 if (FrictionCoupling==1){ 192 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); 191 if (FrictionCoupling==3){ 192 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} 193 else if(FrictionCoupling==4){ 194 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); 193 195 } 194 196 break; … … 202 204 iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum); 203 205 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); 204 if (FrictionCoupling==1){ 205 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); 206 if (FrictionCoupling==3){ 207 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} 208 else if(FrictionCoupling==4){ 209 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); 206 210 } 207 211 break; … … 231 235 iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); 232 236 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); 233 if (FrictionCoupling==1){ 234 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); 237 if (FrictionCoupling==3){ 238 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} 239 else if(FrictionCoupling==4){ 240 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); 235 241 } 236 242 break; … … 453 459 K[2][0]=h/(2.*vel)*fabs(vz*vx); K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz); 454 460 for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j]; 455 GetBAdvecprime(Bprime,element,xyz_list,gauss); 461 GetBAdvecprime(Bprime,element,xyz_list,gauss); 456 462 457 463 TripleMultiply(Bprime,3,numnodes,1, … … 488 494 }/*}}}*/ 489 495 ElementVector* ThermalAnalysis::CreatePVector(Element* element){/*{{{*/ 490 496 491 497 /* Check if ice in element */ 492 498 if(!element->IsIceInElement()) return NULL; … … 701 707 }/*}}}*/ 702 708 void ThermalAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 703 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 709 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 704 710 * For node i, Bi' can be expressed in the actual coordinate system 705 * by: 711 * by: 706 712 * Bi_advec =[ h ] 707 713 * [ h ] … … 730 736 }/*}}}*/ 731 737 void ThermalAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 732 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 738 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 733 739 * For node i, Bi' can be expressed in the actual coordinate system 734 * by: 740 * by: 735 741 * Biprime_advec=[ dh/dx ] 736 742 * [ dh/dy ] … … 759 765 }/*}}}*/ 760 766 void ThermalAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 761 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 767 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 762 768 * For node i, Bi' can be expressed in the actual coordinate system 763 * by: 769 * by: 764 770 * Bi_conduct=[ dh/dx ] 765 771 * [ dh/dy ]
Note:
See TracChangeset
for help on using the changeset viewer.
