Changeset 24336


Ignore:
Timestamp:
11/15/19 20:07:51 (5 years ago)
Author:
Mathieu Morlighem
Message:

CHG: forgot about GIA...

File:
1 edited

Legend:

Unmodified
Added
Removed
  • TabularUnified issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r24335 r24336  
    50275027void       Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x, IssmDouble* y){/*{{{*/
    50285028
    5029         int i;
    5030         int gsize;
    5031         IssmDouble xi,yi,ri,re,area;
    5032         IssmDouble x0,y0;
    50335029        IssmDouble xyz_list[NUMVERTICES][3];
    50345030
    5035         /*thickness averages: */
    5036         IssmDouble* hes=NULL;
    5037         IssmDouble* times=NULL;
    5038         IssmDouble  currenttime;
    5039         int         numtimes;
    5040         Input2* thickness_input=NULL;
    5041 
    50425031        /*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;
    50605033
    50615034        /*output: */
     
    50675040
    50685041        /*how many dofs are we working with here? */
     5042        int gsize;
     5043        IssmDouble yts;
    50695044        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
    50705045        this->parameters->FindParam(&yts,ConstantsYtsEnum);
    50715046
    50725047        /*recover gia solution parameters: */
     5048        int cross_section_shape;
    50735049        this->parameters->FindParam(&cross_section_shape,GiaCrossSectionShapeEnum);
    50745050
    50755051        /*what time is it? :*/
     5052        IssmDouble currenttime;
    50765053        this->parameters->FindParam(&currenttime,TimeEnum);
    50775054
    50785055        /*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,&times,&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,&times,&numtimes,currenttime);
     5069
    50895070
    50905071        /*recover mantle viscosity: */
    5091         mantle_viscosity_input=this->GetInput2(GiaMantleViscosityEnum);
     5072        Input2* mantle_viscosity_input=this->GetInput2(GiaMantleViscosityEnum);
    50925073        if (!mantle_viscosity_input)_error_("mantle viscosity input needed to compute gia deflection!");
    50935074        mantle_viscosity_input->GetInputAverage(&mantle_viscosity);
    50945075
    50955076        /*recover lithosphere thickness: */
    5096         lithosphere_thickness_input=this->GetInput2(GiaLithosphereThicknessEnum);
     5077        Input2* lithosphere_thickness_input=this->GetInput2(GiaLithosphereThicknessEnum);
    50975078        if (!lithosphere_thickness_input)_error_("lithosphere thickness input needed to compute gia deflection!");
    50985079        lithosphere_thickness_input->GetInputAverage(&lithosphere_thickness);
    50995080
    51005081        /*pull area of this Tria: */
    5101         area=this->GetArea();
     5082        IssmDouble area=this->GetArea();
    51025083
    51035084        /*element radius: */
    5104         re=sqrt(area/PI);
     5085        IssmDouble re=sqrt(area/PI);
    51055086
    51065087        /*figure out gravity center of our element: */
    51075088        ::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;
    51105091
    51115092        /*start loading GiaDeflectionCore arguments: */
     
    51265107        arguments.yts=yts;
    51275108
    5128         for(i=0;i<gsize;i++){
     5109        for(int i=0;i<gsize;i++){
    51295110                /*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));
    51325114
    51335115                /*load ri onto arguments for this vertex i: */
     
    51405122                wg->SetValue(i,wi,ADD_VAL);
    51415123                dwgdt->SetValue(i,dwidt,ADD_VAL);
    5142 
    51435124        }
    51445125
Note: See TracChangeset for help on using the changeset viewer.