![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
- Timestamp:
- 06/07/17 10:50:54 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp
r21330 r21759 48 48 IssmDouble vx,vy,vz,vmag; 49 49 IssmDouble dvx[3],dvy[3],dvz[3],dvmag[3]; 50 IssmDouble eps o,epsprime;50 IssmDouble epseff,epsprime; 51 51 int dim; 52 52 IssmDouble *xyz_list = NULL; … … 101 101 dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]); 102 102 } 103 EstarStrainrateQuantities(&eps o,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]);104 lambdas[iv]=EstarLambdaS(eps o,epsprime);103 EstarStrainrateQuantities(&epseff,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]); 104 lambdas[iv]=EstarLambdaS(epseff,epsprime); 105 105 } 106 106 … … 527 527 IssmDouble* tmp=xNew<IssmDouble>(numvertices); 528 528 IssmDouble Delta18oTime; 529 IssmDouble dpermil ;529 IssmDouble dpermil,f; 530 530 IssmDouble time,yts,time_yr,month; 531 531 this->parameters->FindParam(&time,TimeEnum); 532 532 this->parameters->FindParam(&yts,ConstantsYtsEnum); 533 this->parameters->FindParam(&f,SmbFEnum); 533 534 time_yr=floor(time/yts)*yts; 534 535 … … 556 557 /*Compute the temperature and precipitation*/ 557 558 for(int iv=0;iv<numvertices;iv++){ 558 ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil, 559 ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,f, 559 560 &PrecipitationsPresentday[iv*12], &TemperaturesPresentday[iv*12], 560 561 &monthlytemperatures[iv*12], &monthlyprec[iv*12]); … … 652 653 653 654 /*Intermediaries*/ 655 int materialstype; 654 656 IssmDouble dmudB; 655 657 IssmDouble epsilon3d[6];/* epsilon=[exx,eyy,exy,exy,exz,eyz]; */ … … 669 671 } 670 672 /*Get viscosity*/ 671 material->GetViscosity_B(&dmudB,eps_eff); 673 materialstype=this->material->ObjectEnum(); 674 switch(materialstype){ 675 case MaticeEnum: 676 material->GetViscosity_B(&dmudB,eps_eff); 677 break; 678 case MatestarEnum: 679 material->ViscosityBFS(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,vz_input); 680 break; 681 default: _error_("not supported"); 682 } 672 683 673 684 /*Assign output pointer*/ … … 679 690 680 691 /*Intermediaries*/ 692 int materialstype; 681 693 IssmDouble dmudB; 682 694 IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exy,exz,eyz]; */ … … 696 708 } 697 709 /*Get viscosity*/ 698 material->GetViscosity_B(&dmudB,eps_eff); 710 materialstype=this->material->ObjectEnum(); 711 switch(materialstype){ 712 case MaticeEnum: 713 material->GetViscosity_B(&dmudB,eps_eff); 714 break; 715 case MatestarEnum: 716 material->ViscosityBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input); 717 break; 718 default: _error_("not supported"); 719 } 699 720 700 721 /*Assign output pointer*/ … … 706 727 707 728 /*Intermediaries*/ 729 int materialstype; 708 730 IssmDouble dmudB; 709 731 IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ … … 723 745 724 746 /*Get viscosity*/ 725 material->GetViscosity_B(&dmudB,eps_eff); 747 materialstype=this->material->ObjectEnum(); 748 switch(materialstype){ 749 case MaticeEnum: 750 material->GetViscosity_B(&dmudB,eps_eff); 751 break; 752 case MatestarEnum: 753 material->ViscosityBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input); 754 break; 755 default: _error_("not supported"); 756 } 726 757 727 758 /*Assign output pointer*/ … … 1064 1095 if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element"); 1065 1096 input->GetInputValue(pvalue,gauss); 1097 1098 }/*}}}*/ 1099 void Element::GetInputsInterpolations(Vector<IssmDouble>* interpolations){/*{{{*/ 1100 1101 int interpolation; 1102 1103 /*Go through all inputs and assign interpolation in vector*/ 1104 _assert_(this->inputs); 1105 for(int i=0;i<this->inputs->Size();i++){ 1106 Input* input=xDynamicCast<Input*>(this->inputs->GetObjectByOffset(i)); 1107 switch(input->ObjectEnum()){ 1108 case BoolInputEnum: 1109 case DoubleInputEnum: 1110 case IntInputEnum: 1111 interpolations->SetValue(input->InstanceEnum(),reCast<IssmDouble>(input->ObjectEnum()),INS_VAL); 1112 break; 1113 case TriaInputEnum: 1114 interpolation = input->GetResultInterpolation(); 1115 interpolations->SetValue(input->InstanceEnum(),interpolation,INS_VAL); 1116 break; 1117 default: 1118 _error_("Input "<<EnumToStringx(input->ObjectEnum())<<" not supported yet"); 1119 } 1120 } 1066 1121 1067 1122 }/*}}}*/ … … 1161 1216 int numnodes = this->GetNumberOfNodes(); 1162 1217 int* doflist = NULL; 1218 IssmDouble value; 1163 1219 IssmDouble* values = NULL; 1220 Input* input = NULL; 1164 1221 1165 1222 switch(type){ 1166 case VertexPIdEnum: 1167 doflist = xNew<int>(numvertices); 1168 values = xNew<IssmDouble>(numvertices); 1169 /*Fill in values*/ 1170 this->GetVertexPidList(doflist); 1171 this->GetInputListOnVertices(values,input_enum); 1172 vector->SetValues(numvertices,doflist,values,INS_VAL); 1173 break; 1174 case VertexSIdEnum: 1175 doflist = xNew<int>(numvertices); 1176 values = xNew<IssmDouble>(numvertices); 1177 /*Fill in values*/ 1178 this->GetVerticesSidList(doflist); 1179 this->GetInputListOnVertices(values,input_enum); 1180 vector->SetValues(numvertices,doflist,values,INS_VAL); 1181 break; 1182 case NodesEnum: 1183 doflist = xNew<int>(numnodes); 1184 values = xNew<IssmDouble>(numnodes); 1185 /*Fill in values*/ 1186 this->GetInputListOnNodes(values,input_enum); 1187 this->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); 1188 vector->SetValues(numnodes,doflist,values,INS_VAL); 1189 break; 1190 case NodeSIdEnum: 1191 doflist = xNew<int>(numnodes); 1192 values = xNew<IssmDouble>(numnodes); 1193 /*Fill in values*/ 1194 this->GetNodesSidList(doflist); 1195 this->GetInputListOnNodes(values,input_enum); 1196 vector->SetValues(numnodes,doflist,values,INS_VAL); 1197 break; 1198 default: 1199 _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); 1200 } 1201 1223 case ElementSIdEnum: 1224 input=inputs->GetInput(input_enum); _assert_(input); 1225 input->GetInputAverage(&value); 1226 vector->SetValue(this->sid,value,INS_VAL); 1227 break; 1228 case VertexPIdEnum: 1229 doflist = xNew<int>(numvertices); 1230 values = xNew<IssmDouble>(numvertices); 1231 /*Fill in values*/ 1232 this->GetVertexPidList(doflist); 1233 this->GetInputListOnVertices(values,input_enum); 1234 vector->SetValues(numvertices,doflist,values,INS_VAL); 1235 break; 1236 case VertexSIdEnum: 1237 doflist = xNew<int>(numvertices); 1238 values = xNew<IssmDouble>(numvertices); 1239 /*Fill in values*/ 1240 this->GetVerticesSidList(doflist); 1241 this->GetInputListOnVertices(values,input_enum); 1242 vector->SetValues(numvertices,doflist,values,INS_VAL); 1243 break; 1244 case NodesEnum: 1245 doflist = xNew<int>(numnodes); 1246 values = xNew<IssmDouble>(numnodes); 1247 /*Fill in values*/ 1248 this->GetInputListOnNodes(values,input_enum); 1249 this->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); 1250 vector->SetValues(numnodes,doflist,values,INS_VAL); 1251 break; 1252 case NodeSIdEnum: 1253 doflist = xNew<int>(numnodes); 1254 values = xNew<IssmDouble>(numnodes); 1255 /*Fill in values*/ 1256 this->GetNodesSidList(doflist); 1257 this->GetInputListOnNodes(values,input_enum); 1258 vector->SetValues(numnodes,doflist,values,INS_VAL); 1259 break; 1260 default: 1261 _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); 1262 } 1263 1202 1264 /*Clean up*/ 1203 1265 xDelete<int>(doflist); … … 1569 1631 name==OldGradientEnum || 1570 1632 name==ConvergedEnum || 1633 name==MaterialsRheologyEEnum || 1634 name==MaterialsRheologyEbarEnum || 1571 1635 name==MaterialsRheologyBEnum || 1572 1636 name==MaterialsRheologyBbarEnum || 1573 1637 name==MaterialsRheologyNEnum || 1574 name==MaterialsRheologyKoEnum ||1575 name==MaterialsRheologyKobarEnum ||1576 1638 name==MaterialsRheologyEcEnum || 1577 1639 name==MaterialsRheologyEcbarEnum || … … 1586 1648 name==SealevelriseDeltathicknessEnum || 1587 1649 name==SealevelriseCumDeltathicknessEnum || 1650 name==EsaUmotionEnum || 1651 name==EsaNmotionEnum || 1652 name==EsaEmotionEnum || 1653 name==EsaStrainratexxEnum || 1654 name==EsaStrainratexyEnum || 1655 name==EsaStrainrateyyEnum || 1656 name==EsaRotationrateEnum || 1657 name==EsaDeltathicknessEnum || 1588 1658 name==GiaWEnum || 1589 1659 name==GiadWdtEnum || … … 1602 1672 name==HydrologyHeadOldEnum || 1603 1673 name==StressbalanceConvergenceNumStepsEnum || 1604 name==MeshVertexonbaseEnum 1674 name==MeshVertexonbaseEnum || 1675 name==FrictionPEnum || 1676 name==FrictionQEnum || 1677 name==FrictionCoefficientcoulombEnum || 1678 name==LoadingforceXEnum || 1679 name==LoadingforceYEnum || 1680 name==VelEnum || 1681 name==VxPicardEnum || 1682 name==VyPicardEnum 1605 1683 1606 1684 ) { … … 2013 2091 desfac, s0t[iv], s0p[iv],rlaps,rlapslgm,TdiffTime,sealevTime, 2014 2092 rho_water,rho_ice); 2015 /*Get yearlytemperatures */2093 /*Get yearlytemperatures */ 2016 2094 for(int month=0;month<12;month++) { 2017 2095 yearlytemperatures[iv]=yearlytemperatures[iv]+(monthlytemperatures[iv*12+month]+273.15)*mavg; // Has to be in Kelvin 2018 2096 } 2019 2097 } … … 2036 2114 switch(this->ObjectEnum()){ 2037 2115 case TriaEnum: 2038 2039 2116 // this->inputs->AddInput(new TriaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); 2117 this->inputs->AddInput(new TriaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum)); 2040 2118 this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&agd[0],P1Enum)); 2041 2119 this->inputs->AddInput(new TriaInput(SmbAccumulationEnum,&accu[0],P1Enum)); … … 2043 2121 break; 2044 2122 case PentaEnum: 2045 if(IsOnSurface()){ 2046 GetInputListOnVertices(&s[0],TemperatureEnum); 2047 yearlytemperatures[0] = s[0]; 2048 yearlytemperatures[1] = s[1]; 2049 yearlytemperatures[2] = s[2]; 2050 this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); 2051 } 2052 this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&agd[0],P1Enum)); 2053 this->inputs->AddInput(new PentaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum)); 2054 this->InputExtrude(TemperaturePDDEnum,-1); 2055 this->InputExtrude(SmbMassBalanceEnum,-1); 2056 break; 2123 if(IsOnSurface()){ 2124 /*Here, we want to change the BC of the thermal model, keep 2125 * the temperatures as they are for the base of the penta and 2126 * yse yearlytemperatures for the top*/ 2127 GetInputListOnVertices(&s[0],TemperatureEnum); 2128 yearlytemperatures[0] = s[0]; 2129 yearlytemperatures[1] = s[1]; 2130 yearlytemperatures[2] = s[2]; 2131 this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); 2132 2133 bool isenthalpy; 2134 this->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); 2135 if(isenthalpy){ 2136 /*Convert that to enthalpy for the enthalpy model*/ 2137 IssmDouble enthalpy[6]; 2138 GetInputListOnVertices(&enthalpy[0],EnthalpyEnum); 2139 ThermalToEnthalpy(&enthalpy[3],yearlytemperatures[3],0.,0.); 2140 ThermalToEnthalpy(&enthalpy[4],yearlytemperatures[4],0.,0.); 2141 ThermalToEnthalpy(&enthalpy[5],yearlytemperatures[5],0.,0.); 2142 this->inputs->AddInput(new PentaInput(EnthalpyEnum,&enthalpy[0],P1Enum)); 2143 } 2144 } 2145 this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&agd[0],P1Enum)); 2146 this->inputs->AddInput(new PentaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum)); 2147 this->InputExtrude(TemperaturePDDEnum,-1); 2148 this->InputExtrude(SmbMassBalanceEnum,-1); 2149 break; 2057 2150 case TetraEnum: 2058 2059 2060 2061 2062 2063 2064 2065 2066 2067 2068 2069 2151 if(IsOnSurface()){ 2152 GetInputListOnVertices(&s[0],TemperatureEnum); 2153 yearlytemperatures[0] = s[0]; 2154 yearlytemperatures[1] = s[1]; 2155 yearlytemperatures[2] = s[2]; 2156 this->inputs->AddInput(new TetraInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); 2157 } 2158 this->inputs->AddInput(new TetraInput(SmbMassBalanceEnum,&agd[0],P1Enum)); 2159 this->inputs->AddInput(new TetraInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum)); 2160 this->InputExtrude(TemperaturePDDEnum,-1); 2161 this->InputExtrude(SmbMassBalanceEnum,-1); 2162 break; 2070 2163 default: _error_("Not implemented yet"); 2071 2164 } … … 2123 2216 case DeviatoricStresszzEnum: 2124 2217 case DeviatoricStresseffectiveEnum: this->ComputeDeviatoricStressTensor(); break; 2218 case EsaStrainratexxEnum: 2219 case EsaStrainratexyEnum: 2220 case EsaStrainrateyyEnum: 2221 case EsaRotationrateEnum: this->ComputeEsaStrainAndVorticity(); break; 2125 2222 case SigmaNNEnum: this->ComputeSigmaNN(); break; 2126 2223 case LambdaSEnum: this->ComputeLambdaS(); break; … … 2662 2759 xDelete<IssmDouble>(a); 2663 2760 xDelete<IssmDouble>(T); 2761 xDelete<IssmDouble>(dzini); 2762 xDelete<IssmDouble>(dini); 2763 xDelete<IssmDouble>(reini); 2764 xDelete<IssmDouble>(gdnini); 2765 xDelete<IssmDouble>(gspini); 2766 xDelete<IssmDouble>(Wini); 2767 xDelete<IssmDouble>(aini); 2768 xDelete<IssmDouble>(Tini); 2769 2664 2770 delete gauss; 2665 2771 /*}}}*/ 2666 2772 } 2667 2773 /*}}}*/ 2774 void Element::StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ 2775 2776 /*Intermediaries*/ 2777 IssmDouble dvx[3]; 2778 IssmDouble dvy[3]; 2779 2780 /*Check that both inputs have been found*/ 2781 if(!vx_input || !vy_input){ 2782 _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n"); 2783 } 2784 2785 /*Get strain rate assuming that epsilon has been allocated*/ 2786 vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); 2787 vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); 2788 epsilon[0] = dvx[0]; // normal strain rate x-direction 2789 epsilon[1] = dvy[1]; // normal strain rate y-direction 2790 epsilon[2] = 0.5*(dvx[1] + dvy[0]); // shear strain rate 2791 epsilon[3] = 0.5*(dvx[1] - dvy[0]); // rotation rate 2792 2793 }/*}}}*/ 2668 2794 void Element::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ 2669 2795 /*Compute the 3d Strain Rate (6 components):
Note:
See TracChangeset
for help on using the changeset viewer.
