![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
source:
issm/oecreview/Archive/22819-23185/ISSM-22873-22874.diff
| Last change on this file was 23186, checked in by , 7 years ago | |
|---|---|
| File size: 33.7 KB | |
-
../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
11 11 int M,N; 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*/ 17 17 IssmDouble *spcvector = NULL; … … 93 93 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); 100 100 iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); … … 112 112 /*Update elements: */ 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++){ 118 118 if(iomodel->my_elements[i]){ … … 188 188 iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); 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; 195 197 case 2: … … 201 203 iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum); 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; 208 212 case 4: … … 230 234 iomodel->FetchDataToInput(elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum); 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; 237 243 case 9: … … 452 458 K[1][0]=h/(2.*vel)*fabs(vy*vx); K[1][1]=h/(2.*vel)*fabs(vy*vy); K[1][2]=h/(2.*vel)*fabs(vy*vz); 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, 458 464 &K[0][0],3,3,0, … … 487 493 return Ke; 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; 493 499 … … 700 706 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 ] 708 714 * [ h ] … … 729 735 xDelete<IssmDouble>(basis); 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 ] 737 743 * [ dh/dz ] … … 758 764 xDelete<IssmDouble>(dbasis); 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 ] 766 772 * [ dh/dz ] -
../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
610 610 } 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; 616 616 case Domain2DverticalEnum: numdofs=1; break; … … 784 784 iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); 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; 791 793 case 2: … … 797 799 iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum); 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; 804 808 case 4: … … 826 830 iomodel->FetchDataToInput(elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum); 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; 833 840 case 9: … … 933 940 else if(newton>0) 934 941 solutionsequence_newton(femmodel); 935 942 else 936 solutionsequence_nonlinear(femmodel,conserve_loads); 943 solutionsequence_nonlinear(femmodel,conserve_loads); 937 944 } 938 else if(!isFS && (isSSA || isHO || isL1L2)){ 945 else if(!isFS && (isSSA || isHO || isL1L2)){ 939 946 if(VerboseSolution()) _printf0_(" computing velocities\n"); 940 947 941 948 femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); … … 942 949 if(newton>0) 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){ 948 955 femmodel->parameters->SetParam(VxEnum,InputToExtrudeEnum); … … 965 972 int approximation; 966 973 element->GetInputValue(&approximation,ApproximationEnum); 967 974 switch(approximation){ 968 case FSApproximationEnum: 975 case FSApproximationEnum: 969 976 return CreateDVectorFS(element); 970 977 default: 971 978 return NULL; //no need for doftypes outside of FS approximation … … 978 985 int approximation; 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: 988 995 return NULL; … … 996 1003 switch(approximation){ 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: 1014 1021 return NULL; … … 1023 1030 switch(approximation){ 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: 1041 1048 return NULL; … … 1125 1132 case FSApproximationEnum: 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; 1136 1143 case L1L2ApproximationEnum: … … 1580 1587 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); 1586 1593 … … 1688 1695 return pe; 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 * 1700 1707 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes) … … 1728 1735 xDelete<IssmDouble>(dbasis); 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 1736 1743 * [ 0 N ] … … 1765 1772 xDelete<IssmDouble>(basis); 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 ] 1773 1780 * [ dN/dx 2*dN/dy ] … … 2049 2056 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); 2055 2062 … … 2676 2683 return pe; 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 * 2691 2698 * We assume B has been allocated already, of size: 5x(NDOF2*numnodes) … … 2725 2732 xDelete<IssmDouble>(dbasis); 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 2733 2740 * [ 0 N ] … … 2762 2769 xDelete<IssmDouble>(basis); 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 * 2774 2781 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes) … … 3083 3090 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); 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++){ 3089 3096 Ke->values[numdof*((i+1)*dim-1)+(j+1)*dim-1] += scalar*vbasis[i]*vbasis[j]; … … 3285 3292 BtBUzawa,pnumdof,numdof,0, 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); 3291 3298 element->NormalBase(&normal[0],xyz_list_base); … … 4084 4091 for(int ig=gauss->begin();ig<gauss->end();ig++){ 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); 4090 4097 … … 4095 4102 } 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); 4101 4108 … … 4304 4311 return pe; 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 ] 4310 4317 * [ 0 dphi/dy ] … … 4416 4423 xDelete<IssmDouble>(pbasis); 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 ] 4424 4431 * in 2d: … … 4476 4483 xDelete<IssmDouble>(vbasis); 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 ] 4484 4491 * [ dphi/dy dphi/dx ] … … 4589 4596 xDelete<IssmDouble>(pbasis); 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 * 4597 4604 * In 3d … … 4625 4632 xDelete<IssmDouble>(vdbasis); 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 ] 4633 4640 * [ dphi/dy dphi/dx ] … … 4688 4695 xDelete<IssmDouble>(vdbasis); 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 4694 4701 /*Fetch number of nodes for this finite element*/ … … 4710 4717 xDelete<IssmDouble>(basis); 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 ] 4716 4723 * [ 0 dphi/dy ] … … 4775 4782 }/*}}}*/ 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 4781 4788 /*Fetch number of nodes for this finite element*/ … … 4795 4802 xDelete<IssmDouble>(basis); 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 ] 4801 4808 * … … 4929 4936 /*Allocate new inputs*/ 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 4943 4950 /*Get d and tau*/ … … 5171 5178 } 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]); 5179 5186 } … … 5182 5189 IssmDouble coef1,coef2,coef3; 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){ 5188 5195 coef3 = sqrt( … … 5206 5213 dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + 2.*epsxy_old*epsxy_old ); 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 } 5212 5219 /*Initial guess cannot be 0 otherwise log(0) - inf*/ … … 5359 5366 } 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]); 5367 5374 } … … 5613 5620 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); 5619 5626 … … 5646 5653 IssmDouble Bprime2[3][numdofs]; 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 5653 5660 /*Find penta on bed as FS must be coupled to the dofs on the bed: */ … … 5718 5725 &Bprime2[0][0],3,numdofs,0, 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]; 5724 5731 … … 5802 5809 this->GetBHOFriction(L,element,3,xyz_list_tria,gauss); 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: */ 5808 5815 TripleMultiply( L,2,numdof,1, … … 5883 5890 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 5889 5896 D_scalar=2*viscosity*gauss->weight*Jdet; … … 5893 5900 &D[0][0],3,3,0, 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]; 5899 5906 … … 5902 5909 5903 5910 /*Clean-up and return*/ 5904 5911 basaltria->DeleteMaterials(); delete basaltria; 5905 5912 5906 5913 delete gauss; 5907 5914 delete gauss_tria; 5908 5915 xDelete<IssmDouble>(B); … … 5940 5947 int indices[3]={18,19,20}; 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); 5948 5955 … … 5986 5993 /*Initialize Element matrix and return if necessary*/ 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: */ 5992 5999 Element* basalelement = element->SpawnBasalElement(); … … 6222 6229 element->JacobianDeterminant(&Jdet, xyz_list,gauss); 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); 6228 6235 … … 6468 6475 return pe; 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 ] 6476 6483 * [ dh/dy dh/dx 0 0 ] … … 6517 6524 } 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 ] 6525 6532 * [ dN/dy dN/dx ] … … 6552 6559 xDelete<IssmDouble>(dbasis); 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 ] 6560 6567 * [ 1/2*dh/dy 1/2*dh/dx 0 0 ] … … 6610 6617 } 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 ] 6618 6625 * [ 1/2*dN/dy 1/2*dN/dx ] … … 6642 6649 xDelete<IssmDouble>(dbasis); 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 ] 6650 6657 * [ 1/2*dh/dy 1/2*dh/dx ] … … 6673 6680 xDelete<IssmDouble>(dbasis); 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 ] 6681 6688 * [ h 0 ] … … 6707 6714 } 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] 6715 6722 * [ 0 0 h 0] … … 6812 6819 } 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 ] 6820 6827 * [ 0 0 h ] … … 6852 6859 }/*}}}*/ 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 ] 6860 6867 * [ h 0 ] … … 7157 7164 element->GetInputListOnNodes(&vz[0],VzEnum,0.); 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); 7163 7170 g = element->GetMaterialParameter(ConstantsGEnum);
Note:
See TracBrowser
for help on using the repository browser.
