Ignore:
Timestamp:
08/07/18 10:22:46 (7 years ago)
Author:
Mathieu Morlighem
Message:

CHG: removing double white lines

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r23059 r23066  
    106106        }
    107107        else tria->nodes = NULL;
    108        
     108
    109109        tria->vertices = (Vertex**)this->hvertices->deliverp();
    110110        tria->material = (Material*)this->hmaterial->delivers();
    … …  
    115115/*}}}*/
    116116void Tria::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
    117        
     117
    118118        MARSHALLING_ENUM(TriaEnum);
    119        
     119
    120120        /*Call parent classes: */
    121121        ElementHook::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    … …  
    135135        /*Call inputs method*/
    136136        _assert_(this->inputs);
    137        
     137
    138138        int domaintype;
    139139        parameters->FindParam(&domaintype,DomainTypeEnum);
    … …  
    185185
    186186        int        tria_vertex_ids[3];
    187        
     187
    188188        for(int k=0;k<3;k++){
    189189                tria_vertex_ids[k]=reCast<int>(iomodel->elements[3*this->Sid()+k]); //ids for vertices are in the elements array from Matlab
    … …  
    328328/*}}}*/
    329329void       Tria::CalvingCrevasseDepth(){/*{{{*/
    330        
     330
    331331        IssmDouble  xyz_list[NUMVERTICES][3];
    332332        IssmDouble  calvingrate[NUMVERTICES];
    … …  
    340340        /* Get node coordinates and dof list: */
    341341        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    342                
     342
    343343        /*Get the critical fraction of thickness surface and basal crevasses penetrate for calving onset*/
    344344        this->parameters->FindParam(&critical_fraction,CalvingCrevasseDepthEnum);
    345                
     345
    346346        IssmDouble rho_ice        = this->GetMaterialParameter(MaterialsRhoIceEnum);
    347347        IssmDouble rho_seawater   = this->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    … …  
    361361        Input*   s_xy_input              = inputs->GetInput(DeviatoricStressxyEnum);     _assert_(s_xy_input);
    362362        Input*   s_yy_input              = inputs->GetInput(DeviatoricStressyyEnum);     _assert_(s_yy_input);
    363        
     363
    364364        /*Loop over all elements of this partition*/
    365365        GaussTria* gauss=new GaussTria();
    366366        for (int iv=0;iv<NUMVERTICES;iv++){
    367367                gauss->GaussVertex(iv);
    368        
     368
    369369                H_input->GetInputValue(&thickness,gauss);
    370370                bed_input->GetInputValue(&bed,gauss);
    … …  
    378378                s_xy_input->GetInputValue(&s_xy,gauss);
    379379                s_yy_input->GetInputValue(&s_yy,gauss);
    380                
     380
    381381                vel=sqrt(vx*vx+vy*vy)+1.e-14;
    382382
    … …  
    384384                s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
    385385                if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}
    386                
     386
    387387                Ho = thickness - (rho_seawater/rho_ice) * (-bed);
    388388                if(Ho<0.)  Ho=0.;
    … …  
    399399                //      water_height = surface_crevasse[iv];
    400400                //}
    401                
     401
    402402                /*basal crevasse*/
    403403                //basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (rheology_B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho);
    … …  
    405405                if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
    406406                if (bed>0.) basal_crevasse[iv] = 0.;
    407        
     407
    408408                H_surf = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height - critical_fraction*float_depth;
    409409                H_surfbasal = (surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv])-(critical_fraction*thickness);
    410                
     410
    411411                crevasse_depth[iv] = max(H_surf,H_surfbasal);
    412412        }
    413        
     413
    414414        this->inputs->AddInput(new TriaInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1Enum));
    415415        this->inputs->AddInput(new TriaInput(BasalCrevasseEnum,&basal_crevasse[0],P1Enum));
    … …  
    461461                else
    462462                        calvingrate[iv]=0.;
    463                
     463
    464464                calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14);
    465465                calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14);
    … …  
    562562        IssmDouble  vorticity_xy[NUMVERTICES];
    563563        GaussTria*  gauss=NULL;
    564        
     564
    565565        /* Get node coordinates and dof list: */
    566566        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    … …  
    569569        Input* vx_input=this->GetInput(EsaXmotionEnum); _assert_(vx_input);
    570570        Input* vy_input=this->GetInput(EsaYmotionEnum); _assert_(vy_input);
    571        
     571
    572572        /* Start looping on the number of vertices: */
    573573        gauss=new GaussTria();
    … …  
    773773        }
    774774
    775                
    776775}/*}}}*/
    777776void       Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
    … …  
    14001399/*}}}*/
    14011400void       Tria::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    1402        
     1401
    14031402        /* Intermediaries */
    14041403        int i, dir,nrfrontnodes;
    … …  
    14891488}/*}}}*/
    14901489void       Tria::GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/
    1491        
     1490
    14921491        /* GetLevelsetIntersection computes:
    14931492         * 1. indices of element, sorted in [iceverts, noiceverts] in counterclockwise fashion,
    … …  
    15651564/*}}}*/
    15661565void       Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/
    1567        
     1566
    15681567        /*Computeportion of the element that has a positive levelset*/
    15691568
    … …  
    16791678
    16801679        parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
    1681        
     1680
    16821681        /*Cast to Controlinput*/
    16831682        if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
    16841683        ControlInput* controlinput = xDynamicCast<ControlInput*>(input);
    1685        
     1684
    16861685        if(strcmp(data,"value")==0){
    16871686                input  = controlinput->values;
    … …  
    17001699        }
    17011700        /*Check what input we are dealing with*/
    1702        
     1701
    17031702        switch(input->ObjectEnum()){
    17041703                case TriaInputEnum:
    … …  
    17171716                                break;
    17181717                          }
    1719                        
     1718
    17201719                case TransientInputEnum:
    17211720                                {
    … …  
    19821981        base_input->GetInputAverage(&bed);
    19831982        bed_input->GetInputAverage(&bathymetry);
    1984        
     1983
    19851984        /*Return: */
    19861985        return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.));
    … …  
    20262025        if(control_analysis && ad_analysis) iomodel->FindConstant(&num_control_type,"md.autodiff.num_independent_objects");
    20272026        if(control_analysis && ad_analysis) iomodel->FindConstant(&num_responses,"md.autodiff.num_dependent_objects");
    2028 
    2029 
    20302027
    20312028        /*Recover vertices ids needed to initialize inputs*/
    … …  
    23942391IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/
    23952392
    2396 
    23972393        /*intermediary: */
    23982394        IssmDouble* values=NULL;
    … …  
    24072403        IssmDouble  fraction1,fraction2;
    24082404        bool        mainlynegative=true;
    2409        
     2405
    24102406        /*Output:*/
    24112407        volume=0;
    … …  
    24162412        /*Retrieve inputs required:*/
    24172413        thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
    2418        
     2414
    24192415        /*Retrieve material parameters: */
    24202416        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    … …  
    24252421                values[i]=levelset[this->vertices[i]->Sid()];
    24262422        }
    2427                
     2423
    24282424        /*Ok, use the level set values to figure out where we put our gaussian points:*/
    24292425        this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values);
    … …  
    28772873        dist_cf_input->GetInputValue(&dist_cf,gauss);
    28782874        delete gauss;
    2879        
     2875
    28802876        /*Ensure values are positive for floating ice*/
    28812877        dist_gl = fabs(dist_gl);
    … …  
    31973193        }
    31983194
    3199 
    32003195}
    32013196/*}}}*/
    … …  
    32783273        int         idlist[NUMVERTICES],control_init;
    32793274
    3280 
    32813275        /*Get Domain type*/
    32823276        int domaintype;
    … …  
    32983292        /*Get out if this is not an element input*/
    32993293        if(!IsInputEnum(control_enum)) return;
    3300        
     3294
    33013295        Input* input     = (Input*)this->inputs->GetInput(control_enum);   _assert_(input);
    33023296        if(input->ObjectEnum()!=ControlInputEnum){
    … …  
    33303324        IssmDouble  values[NUMVERTICES];
    33313325        int         vertexpidlist[NUMVERTICES],control_init;
    3332 
    33333326
    33343327        /*Get Domain type*/
    … …  
    41674160        IssmDouble* H_elastic_precomputed = NULL;
    41684161        int         M, hemi;
    4169        
     4162
    41704163        /*computation of Green functions:*/
    41714164        IssmDouble* U_elastic= NULL;
    … …  
    41744167        IssmDouble* X_elastic= NULL;
    41754168        IssmDouble* Y_elastic= NULL;
    4176        
     4169
    41774170        /*optimization:*/
    41784171        bool store_green_functions=false;
    … …  
    41824175        if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!");
    41834176        deltathickness_input->GetInputAverage(&I);
    4184                
     4177
    41854178        /*early return if we are not on the (ice) loading point: */
    41864179        if(I==0) return;
    … …  
    41954188        /*which hemisphere? for north-south, east-west components*/
    41964189        this->parameters->FindParam(&hemi,EsaHemisphereEnum);
    4197        
     4190
    41984191        /*compute area of element:*/
    41994192        area=GetArea();
    … …  
    42544247                Y_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*Y_elastic[i];
    42554248                X_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*X_elastic[i];
    4256                
     4249
    42574250                /*North-south, East-west components */
    42584251                if (hemi == -1) {
    … …  
    42774270        pX->SetValues(gsize,indices,X_values,ADD_VAL);
    42784271        pY->SetValues(gsize,indices,Y_values,ADD_VAL);
    4279        
     4272
    42804273        /*free ressources:*/
    42814274        xDelete<int>(indices);
    … …  
    43074300        IssmDouble* H_elastic_precomputed = NULL;
    43084301        int         M;
    4309        
     4302
    43104303        /*computation of Green functions:*/
    43114304        IssmDouble* U_elastic= NULL;
    43124305        IssmDouble* N_elastic= NULL;
    43134306        IssmDouble* E_elastic= NULL;
    4314        
     4307
    43154308        /*optimization:*/
    43164309        bool store_green_functions=false;
    … …  
    43204313        if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!");
    43214314        deltathickness_input->GetInputAverage(&I);
    4322                
     4315
    43234316        /*early return if we are not on the (ice) loading point: */
    43244317        if(I==0) return;
    … …  
    44194412                N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    44204413                E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    4421                
     4414
    44224415                /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
    44234416                int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
    … …  
    44344427        pNorth->SetValues(gsize,indices,N_values,ADD_VAL);
    44354428        pEast->SetValues(gsize,indices,E_values,ADD_VAL);
    4436        
     4429
    44374430        /*free ressources:*/
    44384431        xDelete<int>(indices);
    … …  
    44554448
    44564449        if(IsWaterInElement()){
    4457                
     4450
    44584451                IssmDouble area;
    44594452
    … …  
    45284521        if(IsWaterInElement()){
    45294522                IssmDouble rho_water, S;
    4530                
     4523
    45314524                /*recover material parameters: */
    45324525                rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    4533        
     4526
    45344527                /*From Sg_old, recover water sea level rise:*/
    45354528                S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
    4536                
     4529
    45374530                /* Perturbation terms for moment of inertia (moi_list):
    45384531                 * computed analytically (see Wu & Peltier, eqs 10 & 32)
    … …  
    45464539        else if(this->inputs->Max(MaskIceLevelsetEnum)<0){
    45474540                IssmDouble rho_ice, I;
    4548                
     4541
    45494542                /*recover material parameters: */
    45504543                rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    4551        
     4544
    45524545                /*Compute ice thickness change: */
    45534546                Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
    45544547                if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
    45554548                deltathickness_input->GetInputAverage(&I);
    4556                
     4549
    45574550                dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea;
    45584551                dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea;
    45594552                dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea;
    45604553        }
    4561        
     4554
    45624555        return;
    45634556}/*}}}*/
    … …  
    45924585        bool computeelastic= true;
    45934586        bool scaleoceanarea= false;
    4594        
     4587
    45954588        /*early return if we are not on an ice cap:*/
    45964589        if(!(this->inputs->Max(MaskIceLevelsetEnum)<=0)){
    … …  
    46064599                return;
    46074600        }
    4608        
     4601
    46094602        /*If we are here, we are on ice that is fully grounded or half-way to floating: */
    46104603        if ((this->inputs->Min(MaskGroundediceLevelsetEnum))<0){
    46114604                notfullygrounded=true; //used later on.
    46124605        }
    4613                
     4606
    46144607        /*Inform mask: */
    46154608        constant=1; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum));
    … …  
    46364629
    46374630        /* Where is the centroid of this element?:{{{*/
    4638        
     4631
    46394632        /*retrieve coordinates: */
    46404633        ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
    4641        
     4634
    46424635        IssmDouble minlong=400;
    46434636        IssmDouble maxlong=-20;
    … …  
    46534646                if (llr_list[2][1]==0)llr_list[2][1]=360;
    46544647        }
    4655        
     4648
    46564649        // correction at the north pole
    46574650        if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
    46584651        if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
    46594652        if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
    4660                        
     4653
    46614654        //correction at the south pole
    46624655        if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
    … …  
    46844677                area*=phi;
    46854678        }
    4686        
     4679
    46874680        /*Compute ice thickness change: */
    46884681        Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
    … …  
    46964689                int         point1;
    46974690                IssmDouble  fraction1,fraction2;
    4698                
     4691
    46994692                /*Recover portion of element that is grounded*/
    47004693                this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    … …  
    47504743                }
    47514744                pSgi->SetValues(gsize,indices,values,ADD_VAL);
    4752                
     4745
    47534746                /*free ressources:*/
    47544747                xDelete<IssmDouble>(values);
    47554748                xDelete<int>(indices);
    47564749        }
    4757        
     4750
    47584751        /*Assign output pointer:*/
    47594752        _assert_(!xIsNan<IssmDouble>(eustatic));
    … …  
    47804773        IssmDouble* G_elastic_precomputed = NULL;
    47814774        int         M;
    4782        
     4775
    47834776        /*computation of Green functions:*/
    47844777        IssmDouble* G_elastic= NULL;
    … …  
    48574850        longe=longe/180*PI;
    48584851        /*}}}*/
    4859        
     4852
    48604853        if(computeelastic){
    4861        
     4854
    48624855                /*recover elastic green function:*/
    48634856                DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
    … …  
    48974890                }
    48984891        }
    4899        
     4892
    49004893        pSgo->SetValues(gsize,indices,values,ADD_VAL);
    49014894
    … …  
    49304923        IssmDouble* H_elastic_precomputed = NULL;
    49314924        int         M;
    4932        
     4925
    49334926        /*computation of Green functions:*/
    49344927        IssmDouble* U_elastic= NULL;
    … …  
    49574950        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
    49584951        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
    4959        
     4952
    49604953        /*early return if elastic not requested:*/
    49614954        if(!computeelastic) return;
    … …  
    50265019        /*From Sg, recover water sea level rise:*/
    50275020        S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;
    5028        
     5021
    50295022        /*Compute ice thickness change: */
    50305023        Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
    50315024        if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
    50325025        deltathickness_input->GetInputAverage(&I);
    5033                
     5026
    50345027        /*initialize: */
    50355028        U_elastic=xNewZeroInit<IssmDouble>(gsize);
    … …  
    50735066                        E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    50745067                }
    5075                
     5068
    50765069                /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
    50775070                int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
Note: See TracChangeset for help on using the changeset viewer.