Changeset 24336

Timestamp:

11/15/19 20:07:51 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: forgot about GIA...

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp ¶

@@ -5027 +5027 @@
 void       Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x, IssmDouble* y){/*{{{*/
 
-        int i;
-        int gsize;
-        IssmDouble xi,yi,ri,re,area;
-        IssmDouble x0,y0;
         IssmDouble xyz_list[NUMVERTICES][3];
 
-        /*thickness averages: */
-        IssmDouble* hes=NULL;
-        IssmDouble* times=NULL;
-        IssmDouble  currenttime;
-        int         numtimes;
-        Input2* thickness_input=NULL;
-
         /*gia solution parameters:*/
-        int cross_section_shape=0;
-
-        /*gia material parameters: */
-        IssmDouble lithosphere_shear_modulus;
-        IssmDouble lithosphere_density;
-        IssmDouble mantle_shear_modulus;
-        IssmDouble mantle_density;
-        Input2* mantle_viscosity_input=NULL;
-        IssmDouble mantle_viscosity;
-        Input2* lithosphere_thickness_input=NULL;
-        IssmDouble lithosphere_thickness;
-
-        /*ice properties: */
-        IssmDouble rho_ice;
-
-        /*constants: */
-        IssmDouble yts;
+        IssmDouble lithosphere_thickness,mantle_viscosity;
 
         /*output: */
@@ -5067 +5040 @@
 
         /*how many dofs are we working with here? */
+        int gsize;
+        IssmDouble yts;
         this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
         this->parameters->FindParam(&yts,ConstantsYtsEnum);
 
         /*recover gia solution parameters: */
+        int cross_section_shape;
         this->parameters->FindParam(&cross_section_shape,GiaCrossSectionShapeEnum);
 
         /*what time is it? :*/
+        IssmDouble currenttime;
         this->parameters->FindParam(&currenttime,TimeEnum);
 
         /*recover material parameters: */
-        lithosphere_shear_modulus=FindParam(MaterialsLithosphereShearModulusEnum);
-        lithosphere_density=FindParam(MaterialsLithosphereDensityEnum);
-        mantle_shear_modulus=FindParam(MaterialsMantleShearModulusEnum);
-        mantle_density=FindParam(MaterialsMantleDensityEnum);
-        rho_ice=FindParam(MaterialsRhoIceEnum);
-
-        /*pull thickness averages: */
-        thickness_input=this->GetInput2(ThicknessEnum); 
-        if (!thickness_input)_error_("thickness input needed to compute gia deflection!");
-        thickness_input->GetInputAveragesUpToCurrentTime(&hes,&times,&numtimes,currenttime);
+        IssmDouble lithosphere_shear_modulus = FindParam(MaterialsLithosphereShearModulusEnum);
+        IssmDouble lithosphere_density       = FindParam(MaterialsLithosphereDensityEnum);
+        IssmDouble mantle_shear_modulus      = FindParam(MaterialsMantleShearModulusEnum);
+        IssmDouble mantle_density            = FindParam(MaterialsMantleDensityEnum);
+        IssmDouble rho_ice                   = FindParam(MaterialsRhoIceEnum);
+
+        /*pull thickness averages! */
+        TransientInput2* thickness_input=this->inputs2->GetTransientInput(ThicknessEnum); 
+        IssmDouble* hes=NULL;
+        IssmDouble* times=NULL;
+        int         numtimes;
+        _error_("Next line not implemented yet");
+        //thickness_input->GetInputAveragesUpToCurrentTime(&hes,&times,&numtimes,currenttime);
+
 
         /*recover mantle viscosity: */
-        mantle_viscosity_input=this->GetInput2(GiaMantleViscosityEnum);
+        Input2* mantle_viscosity_input=this->GetInput2(GiaMantleViscosityEnum);
         if (!mantle_viscosity_input)_error_("mantle viscosity input needed to compute gia deflection!");
         mantle_viscosity_input->GetInputAverage(&mantle_viscosity);
 
         /*recover lithosphere thickness: */
-        lithosphere_thickness_input=this->GetInput2(GiaLithosphereThicknessEnum);
+        Input2* lithosphere_thickness_input=this->GetInput2(GiaLithosphereThicknessEnum);
         if (!lithosphere_thickness_input)_error_("lithosphere thickness input needed to compute gia deflection!");
         lithosphere_thickness_input->GetInputAverage(&lithosphere_thickness);
 
         /*pull area of this Tria: */
-        area=this->GetArea();
+        IssmDouble area=this->GetArea();
 
         /*element radius: */
-        re=sqrt(area/PI);
+        IssmDouble re=sqrt(area/PI);
 
         /*figure out gravity center of our element: */
         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-        x0=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
-        y0=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
+        IssmDouble x0=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
+        IssmDouble y0=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
 
         /*start loading GiaDeflectionCore arguments: */
@@ -5126 +5107 @@
         arguments.yts=yts;
 
-        for(i=0;i<gsize;i++){
+        for(int i=0;i<gsize;i++){
                 /*compute distance from the center of the tria to the vertex i: */
-                xi=x[i]; yi=y[i];
-                ri=sqrt(pow(xi-x0,2)+pow(yi-y0,2));
+                IssmDouble xi=x[i];
+                IssmDouble yi=y[i];
+                IssmDouble ri=sqrt(pow(xi-x0,2)+pow(yi-y0,2));
 
                 /*load ri onto arguments for this vertex i: */
@@ -5140 +5122 @@
                 wg->SetValue(i,wi,ADD_VAL);
                 dwgdt->SetValue(i,dwidt,ADD_VAL);
-
         }