Changeset 24336
- Timestamp:
- 11/15/19 20:07:51 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified issm/trunk-jpl/src/c/classes/Elements/Tria.cpp ¶
r24335 r24336 5027 5027 void Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x, IssmDouble* y){/*{{{*/ 5028 5028 5029 int i;5030 int gsize;5031 IssmDouble xi,yi,ri,re,area;5032 IssmDouble x0,y0;5033 5029 IssmDouble xyz_list[NUMVERTICES][3]; 5034 5030 5035 /*thickness averages: */5036 IssmDouble* hes=NULL;5037 IssmDouble* times=NULL;5038 IssmDouble currenttime;5039 int numtimes;5040 Input2* thickness_input=NULL;5041 5042 5031 /*gia solution parameters:*/ 5043 int cross_section_shape=0; 5044 5045 /*gia material parameters: */ 5046 IssmDouble lithosphere_shear_modulus; 5047 IssmDouble lithosphere_density; 5048 IssmDouble mantle_shear_modulus; 5049 IssmDouble mantle_density; 5050 Input2* mantle_viscosity_input=NULL; 5051 IssmDouble mantle_viscosity; 5052 Input2* lithosphere_thickness_input=NULL; 5053 IssmDouble lithosphere_thickness; 5054 5055 /*ice properties: */ 5056 IssmDouble rho_ice; 5057 5058 /*constants: */ 5059 IssmDouble yts; 5032 IssmDouble lithosphere_thickness,mantle_viscosity; 5060 5033 5061 5034 /*output: */ … … 5067 5040 5068 5041 /*how many dofs are we working with here? */ 5042 int gsize; 5043 IssmDouble yts; 5069 5044 this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); 5070 5045 this->parameters->FindParam(&yts,ConstantsYtsEnum); 5071 5046 5072 5047 /*recover gia solution parameters: */ 5048 int cross_section_shape; 5073 5049 this->parameters->FindParam(&cross_section_shape,GiaCrossSectionShapeEnum); 5074 5050 5075 5051 /*what time is it? :*/ 5052 IssmDouble currenttime; 5076 5053 this->parameters->FindParam(¤ttime,TimeEnum); 5077 5054 5078 5055 /*recover material parameters: */ 5079 lithosphere_shear_modulus=FindParam(MaterialsLithosphereShearModulusEnum); 5080 lithosphere_density=FindParam(MaterialsLithosphereDensityEnum); 5081 mantle_shear_modulus=FindParam(MaterialsMantleShearModulusEnum); 5082 mantle_density=FindParam(MaterialsMantleDensityEnum); 5083 rho_ice=FindParam(MaterialsRhoIceEnum); 5084 5085 /*pull thickness averages: */ 5086 thickness_input=this->GetInput2(ThicknessEnum); 5087 if (!thickness_input)_error_("thickness input needed to compute gia deflection!"); 5088 thickness_input->GetInputAveragesUpToCurrentTime(&hes,×,&numtimes,currenttime); 5056 IssmDouble lithosphere_shear_modulus = FindParam(MaterialsLithosphereShearModulusEnum); 5057 IssmDouble lithosphere_density = FindParam(MaterialsLithosphereDensityEnum); 5058 IssmDouble mantle_shear_modulus = FindParam(MaterialsMantleShearModulusEnum); 5059 IssmDouble mantle_density = FindParam(MaterialsMantleDensityEnum); 5060 IssmDouble rho_ice = FindParam(MaterialsRhoIceEnum); 5061 5062 /*pull thickness averages! */ 5063 TransientInput2* thickness_input=this->inputs2->GetTransientInput(ThicknessEnum); 5064 IssmDouble* hes=NULL; 5065 IssmDouble* times=NULL; 5066 int numtimes; 5067 _error_("Next line not implemented yet"); 5068 //thickness_input->GetInputAveragesUpToCurrentTime(&hes,×,&numtimes,currenttime); 5069 5089 5070 5090 5071 /*recover mantle viscosity: */ 5091 mantle_viscosity_input=this->GetInput2(GiaMantleViscosityEnum);5072 Input2* mantle_viscosity_input=this->GetInput2(GiaMantleViscosityEnum); 5092 5073 if (!mantle_viscosity_input)_error_("mantle viscosity input needed to compute gia deflection!"); 5093 5074 mantle_viscosity_input->GetInputAverage(&mantle_viscosity); 5094 5075 5095 5076 /*recover lithosphere thickness: */ 5096 lithosphere_thickness_input=this->GetInput2(GiaLithosphereThicknessEnum);5077 Input2* lithosphere_thickness_input=this->GetInput2(GiaLithosphereThicknessEnum); 5097 5078 if (!lithosphere_thickness_input)_error_("lithosphere thickness input needed to compute gia deflection!"); 5098 5079 lithosphere_thickness_input->GetInputAverage(&lithosphere_thickness); 5099 5080 5100 5081 /*pull area of this Tria: */ 5101 area=this->GetArea();5082 IssmDouble area=this->GetArea(); 5102 5083 5103 5084 /*element radius: */ 5104 re=sqrt(area/PI);5085 IssmDouble re=sqrt(area/PI); 5105 5086 5106 5087 /*figure out gravity center of our element: */ 5107 5088 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 5108 x0=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;5109 y0=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;5089 IssmDouble x0=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; 5090 IssmDouble y0=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; 5110 5091 5111 5092 /*start loading GiaDeflectionCore arguments: */ … … 5126 5107 arguments.yts=yts; 5127 5108 5128 for(i =0;i<gsize;i++){5109 for(int i=0;i<gsize;i++){ 5129 5110 /*compute distance from the center of the tria to the vertex i: */ 5130 xi=x[i]; yi=y[i]; 5131 ri=sqrt(pow(xi-x0,2)+pow(yi-y0,2)); 5111 IssmDouble xi=x[i]; 5112 IssmDouble yi=y[i]; 5113 IssmDouble ri=sqrt(pow(xi-x0,2)+pow(yi-y0,2)); 5132 5114 5133 5115 /*load ri onto arguments for this vertex i: */ … … 5140 5122 wg->SetValue(i,wi,ADD_VAL); 5141 5123 dwgdt->SetValue(i,dwidt,ADD_VAL); 5142 5143 5124 } 5144 5125
Note:
See TracChangeset
for help on using the changeset viewer.