Changeset 23066


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

CHG: removing double white lines

Location:
issm/trunk-jpl/src/c
Files:
112 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

    r22105 r23066  
    11531153        /*Clean up and return*/
    11541154        xDelete<int>(responses);
    1155                          
     1155
    11561156}/*}}}*/
    11571157void           AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    … …  
    21122112        IssmDouble *xyz_list= NULL;
    21132113
    2114 
    21152114        /*Fetch number of vertices for this finite element*/
    21162115        int numvertices = basalelement->GetNumberOfVertices();
    … …  
    21312130        Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);    _assert_(adjointx_input);
    21322131        Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);    _assert_(adjointy_input);
    2133 
    2134 
    21352132
    21362133        IssmDouble  q_exp;
    … …  
    21822179                Chi   = vmag/(pow(C_param,n)*pow(Neff,n)*As);
    21832180                Gamma = (Chi/(1.+alpha*pow(Chi,q_exp)));
    2184                
     2181
    21852182                Uder =Neff*C_param/(vmag*vmag*n) *
    21862183                        (Gamma-alpha*q_exp*pow(Chi,q_exp-1.)*Gamma*Gamma* pow(Gamma,(1.-n)/n) -
    21872184                         n* pow(Gamma,1./n));
    2188                
     2185
    21892186                /*Build gradient vector (actually -dJ/dD): */
    21902187                for(int i=0;i<numvertices;i++){
    … …  
    23182315                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    23192316        }
    2320 
    23212317
    23222318        /*Fetch number of nodes and dof for this finite element*/

  • issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp

    r22526 r23066  
    77/*Model processing*/
    88void Balancethickness2Analysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
    9 
    109
    1110        int finiteelement = P1Enum;

  • issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp

    r20690 r23066  
    594594        delete gauss;
    595595
    596 
    597596}/*}}}*/
    598597void           BalancethicknessAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/

  • issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp

    r20690 r23066  
    200200        }
    201201
    202 
    203202        /* Start  looping on the number of gaussian points: */
    204203        Gauss* gauss=basalelement->NewGauss(2);

  • issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp

    r22608 r23066  
    126126        /*Add input*/
    127127        element->AddInput(DamageFEnum,f,element->GetElementType());
    128        
     128
    129129        /*Clean up and return*/
    130130        xDelete<IssmDouble>(f);
    … …  
    170170        for (int i=0;i<numnodes;i++){
    171171                gauss->GaussNode(element->GetElementType(),i);
    172                
     172
    173173                eps_xx_input->GetInputValue(&eps_xx,gauss);
    174174                eps_xy_input->GetInputValue(&eps_xy,gauss);
    … …  
    177177                n_input->GetInputValue(&n,gauss);
    178178                damage_input->GetInputValue(&damage,gauss);
    179        
     179
    180180                /*Calculate principal effective strain rates*/
    181181                eps1=(eps_xx+eps_yy)/2.+sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2));
    … …  
    195195        /*Add input*/
    196196        element->AddInput(DamageFEnum,f,element->GetElementType());
    197        
     197
    198198        /*Clean up and return*/
    199199        xDelete<IssmDouble>(f);
    … …  
    262262        for (int i=0;i<numnodes;i++){
    263263                gauss->GaussNode(element->GetElementType(),i);
    264                
     264
    265265                damage_input->GetInputValue(&damage,gauss);
    266266                tau_xx_input->GetInputValue(&tau_xx,gauss);
    … …  
    314314        /*Add input*/
    315315        element->AddInput(DamageFEnum,f,element->GetElementType());
    316        
     316
    317317        /*Clean up and return*/
    318318        xDelete<IssmDouble>(f);
    … …  
    375375                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    376376                element->NodalFunctions(basis,gauss);
    377                
     377
    378378                vx_input->GetInputValue(&vx,gauss);
    379379                vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
    … …  
    538538        }
    539539
    540 
    541540        /* Start  looping on the number of gaussian points: */
    542541        Gauss* gauss=element->NewGauss(2);

  • issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

    r22511 r23066  
    103103        int frictionlaw,basalforcing_model,materialstype;
    104104        int FrictionCoupling;
    105        
     105
    106106        /*Now, is the model 3d? otherwise, do nothing: */
    107107        if(iomodel->domaintype==Domain2DhorizontalEnum)return;
    … …  
    189189                        _error_("not supported");
    190190        }
    191        
     191
    192192        /*Friction law variables*/
    193193        switch(frictionlaw){
    … …  
    805805                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    806806                element->NodalFunctions(basis,gauss);
    807                
     807
    808808                /*viscous dissipation*/
    809809                element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
    … …  
    823823                        pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
    824824                        d2pressure=0.; // for linear elements, 2nd derivative is zero
    825                        
     825
    826826                        d1H_d1P=0.;
    827827                        for(i=0;i<3;i++) d1H_d1P+=d1enthalpypicard[i]*d1pressure[i];
    … …  
    10561056        int* basalnodeindices=NULL;
    10571057        Element* element= NULL;
    1058        
     1058
    10591059        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    10601060
    … …  
    13881388                gauss->GaussNode(element->GetElementType(),indices[i]);
    13891389                gaussup->GaussNode(element->GetElementType(),indicesup[i]);
    1390                
     1390
    13911391                enthalpy_input->GetInputValue(&enthalpy,gauss);
    13921392                enthalpy_input->GetInputValue(&enthalpyup,gaussup);

  • issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp

    r22609 r23066  
    4040        IssmDouble* love_h=NULL;
    4141        IssmDouble* love_l=NULL;
    42        
     42
    4343        IssmDouble* U_elastic = NULL;
    4444        IssmDouble* U_elastic_local = NULL;
    … …  
    6969        U_elastic=xNew<IssmDouble>(M);
    7070        H_elastic=xNew<IssmDouble>(M);
    71        
     71
    7272        /*compute combined legendre + love number (elastic green function:*/
    7373        m=DetermineLocalSize(M,IssmComm::GetComm());
    … …  
    9696
    9797                        deltalove_U = (love_h[n]-love_h[nl-1]);
    98        
     98
    9999                        /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
    100100                        if(n==0){
    … …  
    199199        /*Default, do nothing*/
    200200        return;
    201        
     201
    202202}/*}}}*/
    203203void           EsaAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/

  • issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp

    r22461 r23066  
    145145                GetB(B,workelement,xyz_list,gauss,dim);
    146146                GetBprime(Bprime,workelement,xyz_list,gauss,dim);
    147                
     147
    148148                D_scalar=gauss->weight*Jdet;
    149149
    … …  
    174174                        else
    175175                                for(i=0;i<dim;i++)      normal[i]=0.;
    176                        
     176
    177177                        for(row=0;row<dim;row++)
    178178                                for(col=0;col<dim;col++)
    … …  
    337337        /* Get parameters */
    338338        element->FindParam(&extvar_enum, ExtrapolationVariableEnum);
    339        
     339
    340340        Input* active_input=element->GetInput(IceMaskNodeActivationEnum); _assert_(active_input);
    341341        Input* extvar_input=element->GetInput(extvar_enum); _assert_(extvar_input);

  • issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp

    r23037 r23066  
    199199        }
    200200
    201 
    202 
    203201        return;
    204202}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

    r22856 r23066  
    260260                                                basis,1,numnodes,0,
    261261                                                &Ke->values[0],1);
    262 
    263262
    264263                        /*Transfer EPL part*/

  • issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

    r22902 r23066  
    213213                return NULL;
    214214        }
    215 
    216215
    217216        /*Intermediaries */
    … …  
    504503        }
    505504
    506 
    507505        /*Fetch number of nodes for this finite element*/
    508506        int numnodes = basalelement->GetNumberOfNodes();

  • issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp

    r23020 r23066  
    295295                n_input->GetInputValue(&n,gauss);
    296296                A=pow(B,-n);
    297                
     297
    298298                /*Compute beta term*/
    299299                if(gap<br)
    … …  
    324324        meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
    325325                _assert_(meltrate>0.);
    326        
     326
    327327                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
    328328                 (
    … …  
    398398                /*Calculate effective pressure*/
    399399                eff_pressure[i] = rho_ice*g*thickness[i] - rho_water*g*(values[i]-bed[i]);
    400        
     400
    401401                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    402402                if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
    … …  
    444444        reynolds_input->GetInputAverage(&reynolds);
    445445        gap_input->GetInputAverage(&gap);
    446        
     446
    447447        /*Compute conductivity*/
    448448        IssmDouble conductivity = pow(gap,3)*g/(12.*NU*(1+OMEGA*reynolds));
    … …  
    453453}/*}}}*/
    454454void HydrologyShaktiAnalysis::UpdateGapHeight(FemModel* femmodel){/*{{{*/
    455 
    456455
    457456        for(int j=0;j<femmodel->elements->Size();j++){
    … …  
    547546                IssmDouble pressure_water = rho_water*g*(head-bed);
    548547                if(pressure_water>pressure_ice) pressure_water = pressure_ice;
    549      
     548
    550549      /* Compute change in sensible heat due to changes in pressure melting point*/
    551550           dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
    … …  
    581580        /*Add new gap as an input*/
    582581        element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum);
    583  
     582
    584583        /*Divide by connectivity, add basal flux as an input*/
    585584        q = q/totalweights;

  • issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp

    r22284 r23066  
    379379}/*}}}*/
    380380
    381 
    382381/*Needed changes to switch to the Johnson formulation*//*{{{*/
    383382/*All the changes are to be done in the velocity computation.

  • issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

    r22987 r23066  
    4343                }
    4444        }
    45        
     45
    4646        iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    4747        iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
    … …  
    177177
    178178        /*Calving threshold*/
    179        
     179
    180180        /*Fetch number of nodes and dof for this finite element*/
    181181        int numnodes    = basalelement->GetNumberOfNodes();
    … …  
    326326                                 }
    327327                                break;
    328                                
     328
    329329                        case CalvingLevermannEnum:
    330330                                calvingratex_input->GetInputValue(&c[0],gauss);
    … …  
    357357                                 }
    358358                                break;
    359                        
     359
    360360                        case CalvingHabEnum:
    361361                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    … …  
    378378                                 }
    379379                                break;
    380                        
     380
    381381                        case CalvingCrevasseDepthEnum:
    382382                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    … …  
    385385
    386386                                if(groundedice<0) meltingrate = 0.;
    387                                
     387
    388388                                norm_dlsf=0.;
    389389                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    … …  
    516516}/*}}}*/
    517517ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
    518        
     518
    519519        if(!element->IsOnBase()) return NULL;
    520520        Element* basalelement = element->SpawnBasalElement();
    … …  
    525525        IssmDouble  lsf;
    526526        IssmDouble* xyz_list = NULL;
    527        
     527
    528528        /*Fetch number of nodes and dof for this finite element*/
    529529        int numnodes = basalelement->GetNumberOfNodes();
    … …  
    532532        ElementVector* pe = basalelement->NewElementVector();
    533533        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
    534        
     534
    535535        if(dt!=0.){
    536536                /*Initialize basis vector*/
    … …  
    623623        // d=|a x b|/|b|
    624624        // with a=q-s0, b=s1-s0
    625        
     625
    626626        /* Intermediaries */
    627627        const int dim=2;
    … …  
    634634                b[i]=s1[i]-s0[i];
    635635        }
    636        
     636
    637637        norm_b=0.;
    638638        for(i=0;i<dim;i++)
    … …  
    671671        IssmDouble  bed,water_depth;
    672672        IssmDouble  levelset;
    673        
     673
    674674        femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
    675675
    … …  
    703703                }
    704704        }
    705        
     705
    706706        if(calvinglaw==CalvingHabEnum){
    707707
    … …  
    709709                femmodel->elements->InputDuplicate(MaskIceLevelsetEnum, DistanceToCalvingfrontEnum);
    710710                femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
    711                
     711
    712712                /*Loop over all elements of this partition*/
    713713                for(int i=0;i<femmodel->elements->Size();i++){
    714714                        Element* element  = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    715                        
     715
    716716                        rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
    717717                        rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    … …  
    744744                }
    745745        }
    746        
     746
    747747        if(calvinglaw==CalvingCrevasseDepthEnum){
    748        
     748
    749749                int                 nflipped,local_nflipped;
    750750                Vector<IssmDouble>* vec_constraint_nodes = NULL;
    751751                IssmDouble* constraint_nodes = NULL;
    752                
     752
    753753                /*Get the DistanceToCalvingfront*/
    754754                femmodel->elements->InputDuplicate(MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
    … …  
    842842                                gauss->GaussNode(element->GetElementType(),in);
    843843                                Node* node=element->GetNode(in);
    844                                
     844
    845845                                if(constraint_nodes[node->Sid()]>0.){
    846846                                        node->ApplyConstraint(0,+1.);

  • issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

    r23014 r23066  
    165165        iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
    166166        iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
    167        
     167
    168168        /*Initialize cumdeltalthickness input*/
    169169        InputUpdateFromConstantx(elements,0.,SealevelriseCumDeltathicknessEnum);
    … …  
    203203                iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
    204204        }
    205        
     205
    206206}/*}}}*/
    207207void MasstransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    … …  
    221221        if(numoutputs)parameters->AddObject(new StringArrayParam(MasstransportRequestedOutputsEnum,requestedoutputs,numoutputs));
    222222        iomodel->DeleteData(&requestedoutputs,numoutputs,"md.masstransport.requested_outputs");
    223        
    224        
    225223
    226224}/*}}}*/
    … …  
    805803        xDelete<IssmDouble>(sealevel);
    806804        xDelete<IssmDouble>(bed);
    807        
     805
    808806        xDelete<int>(doflist);
    809807        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};

  • issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp

    r22968 r23066  
    6969        IssmDouble* love_k=NULL;
    7070        IssmDouble* love_l=NULL;
    71        
     71
    7272        bool elastic=false;
    7373        IssmDouble* G_elastic = NULL;
    … …  
    132132                #endif
    133133
    134                
    135134                /*compute combined legendre + love number (elastic green function:*/
    136135                m=DetermineLocalSize(M,IssmComm::GetComm());
    … …  
    168167                                deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
    169168                                deltalove_U = (love_h[n]-love_h[nl-1]);
    170                
     169
    171170                                /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
    172171                                if(n==0){
    … …  
    230229                xDelete<IssmDouble>(H_elastic_local);
    231230        }
    232        
     231
    233232        /*Transitions: */
    234233        iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions");
    … …  
    250249        iomodel->DeleteData(&requestedoutputs,numoutputs,"md.slr.requested_outputs");
    251250
    252 
    253251}/*}}}*/
    254252

  • issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

    r22852 r23066  
    1919}/*}}}*/
    2020void SmbAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    21        
     21
    2222        int    smb_model;
    2323        bool   isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled;
    24        
     24
    2525        /*Update elements: */
    2626        int counter=0;
    … …  
    3232                }
    3333        }
    34        
     34
    3535        /*Figure out smb model: */
    3636        iomodel->FindConstant(&smb_model,"md.smb.model");
    37                        
     37
    3838        switch(smb_model){
    3939                case SMBforcingEnum:
    … …  
    142142                        _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
    143143        }
    144        
    145        
    146144
    147145}/*}}}*/
    … …  
    154152        IssmDouble *temp = NULL;
    155153        int         N,M;
    156        
     154
    157155        parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
    158        
     156
    159157        iomodel->FindConstant(&smb_model,"md.smb.model");
    160158        iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
    161        
     159
    162160        switch(smb_model){
    163161                case SMBforcingEnum:
    … …  
    198196                          parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,M));
    199197                          iomodel->DeleteData(temp,"md.smb.Pfac");
    200                        
     198
    201199                          iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2);
    202200                          parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,M));
    … …  
    266264        /*Figure out smb model: */
    267265        femmodel->parameters->FindParam(&smb_model,SmbEnum);
    268        
     266
    269267        /*branch to correct module*/
    270268        switch(smb_model){

  • issm/trunk-jpl/src/c/analyses/SmoothAnalysis.cpp

    r18930 r23066  
    153153        }
    154154
    155 
    156155        /* Start  looping on the number of gaussian points: */
    157156        Gauss* gauss=element->NewGauss(2);
    … …  
    161160                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    162161                element->NodalFunctions(basis,gauss);
    163 
    164162
    165163                switch(input_enum){

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r23014 r23066  
    27362736        }
    27372737
    2738 
    27392738        /*Clean-up*/
    27402739        xDelete<IssmDouble>(dbasis);

  • issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp

    r22576 r23066  
    2929        else
    3030         iscoupling = false;
    31 
    3231
    3332        /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/
    … …  
    159158ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixBase(Element* element){/*{{{*/
    160159
    161 
    162160        if(!element->IsOnBase()) return NULL;
    163161
    … …  
    199197}/*}}}*/
    200198ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/
    201 
    202199
    203200        if(!element->IsOnSurface()) return NULL;

  • issm/trunk-jpl/src/c/bamg/BamgQuadtree.cpp

    r21629 r23066  
    285285                                int yiv = b->v[k]->i.y;
    286286
    287 
    288287                                int h0 = Norm(xi2,xiv,yi2,yiv);
    289288

  • issm/trunk-jpl/src/c/bamg/Mesh.cpp

    r22380 r23066  
    501501                int  connectivitymax_1=0;
    502502                int verbose=0;
    503        
     503
    504504                /*Check needed pointer*/
    505505                _assert_(bamgmesh);
    … …  
    10341034                //   s0       tt2       s1
    10351035                //-------------------------------
    1036                
     1036
    10371037                /*Intermediaries*/
    10381038                Triangle* tt[3];       //the three triangles
    … …  
    11581158                int verbose=0;
    11591159                double lminaniso = 1./ (Max(hminaniso*hminaniso,1e-100));
    1160        
     1160
    11611161                //Get options
    11621162                if(bamgopts) verbose=bamgopts->verbose;
    1163                        
     1163
    11641164                //display info
    11651165                if (verbose > 1)  _printf_("   BoundAnisotropy by " << anisomax << "\n");
    … …  
    18541854                /*Check pointer*/
    18551855                _assert_(bamgopts);
    1856                
     1856
    18571857                /*Recover options*/
    18581858                verbose=bamgopts->verbose;
    … …  
    27282728                if (verbose>4) _printf_("      Prime Number = "<<PrimeNumber<<"\n");
    27292729                if (verbose>4) _printf_("      k0 = "<<k0<<"\n");
    2730                
     2730
    27312731                //Build orderedvertices
    27322732                for (i=0; i<nbv; i++){
    … …  
    29532953        void  Mesh::MaxSubDivision(BamgOpts* bamgopts,double maxsubdiv) {/*{{{*/
    29542954                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MaxSubDivision)*/
    2955                
     2955
    29562956                /*Intermediaries*/
    29572957                int verbose=0;
    … …  
    37403740        void Mesh::SmoothingVertex(BamgOpts* bamgopts,int nbiter,double omega ) { /*{{{*/
    37413741                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SmoothingVertex)*/
    3742                
     3742
    37433743                /*Intermediaries*/
    37443744                int verbose=0;
    3745        
     3745
    37463746                /*Get options*/
    37473747                if(bamgopts) verbose=bamgopts->verbose;
    … …  
    37853785        void Mesh::SmoothMetric(BamgOpts* bamgopts,double raisonmax) { /*{{{*/
    37863786                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/SmoothMetric)*/
    3787                
     3787
    37883788                /*Intermediaries*/
    37893789                int verbose=0;
    … …  
    40874087                int i;
    40884088                Metric M1(1);
    4089        
     4089
    40904090                /*Get options*/
    40914091                if(bamgopts) verbose=bamgopts->verbose;

  • issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp

    r23065 r23066  
    101101/*Mesh refinement methods*/
    102102void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int** pdatalist,double** pxylist,int** pelementslist){/*{{{*/
    103        
     103
    104104        /*IMPORTANT! pelementslist are in Matlab indexing*/
    105105        /*NEOPZ works only in C indexing*/
    … …  
    118118        /*Intermediaries*/
    119119        bool verbose=VerboseSolution();
    120        
     120
    121121        /*Execute refinement*/
    122122        this->RefineMeshOneLevel(verbose,gl_distance,if_distance,deviatoricerror,thicknesserror);
    123    
     123
    124124        /*Get new geometric mesh in ISSM data structure*/
    125125        this->GetMesh(this->previousmesh,pdatalist,pxylist,pelementslist);
    … …  
    130130/*}}}*/
    131131void AdaptiveMeshRefinement::RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror){/*{{{*/
    132        
     132
    133133        /*Intermediaries*/
    134134        int nelem                                                       =-1;
    … …  
    204204        }
    205205        /*}}}*/
    206        
     206
    207207        /*Now, detele the special elements{{{*/
    208208        if(this->refinement_type==1) this->DeleteSpecialElements(verbose,gmesh);
    … …  
    216216        }
    217217        /*}}}*/
    218        
     218
    219219        /*Unrefinement process: loop over level of refinements{{{*/
    220220        if(verbose) _printf_("\tunrefinement process...\n");
    221221        if(verbose) _printf_("\ttotal: ");
    222222        count=0;
    223          
     223
    224224        nelem=gmesh->NElements();//must keep here
    225225        for(int i=0;i<nelem;i++){
    … …  
    266266        gmesh->BuildConnectivity();
    267267        /*}}}*/
    268        
     268
    269269        /*Refinement process: loop over level of refinements{{{*/
    270270        if(verbose) _printf_("\trefinement process (level max = "<<this->level_max<<")\n");
    … …  
    320320         * */
    321321        if(!geoel) _error_("geoel is NULL!\n");
    322        
     322
    323323        /*Output*/
    324324        int type=0;
    325325        int count=0;
    326        
     326
    327327        /*Intermediaries*/
    328328        TPZVec<TPZGeoEl*> sons;
    329        
     329
    330330        /*Loop over neighboors (sides 3, 4 and 5)*/
    331331        for(int j=3;j<6;j++){
    … …  
    335335                if(sons.size()>4) count++; //if neighbour's level is > element level+1
    336336        }
    337        
     337
    338338        /*Verify and return*/
    339339        if(count>1) type=2;
    340340        else type=count;
    341        
     341
    342342        return type;
    343343}
    … …  
    391391/*}}}*/
    392392void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/
    393    
     393
    394394        /*Now, insert special elements to avoid hanging nodes*/
    395395        if(verbose) _printf_("\trefining to avoid hanging nodes (total: ");
    396        
     396
    397397        /*Intermediaries*/
    398398        int nelem=-1;
    399399        int count= 1;
    400        
     400
    401401        while(count>0){
    402402                nelem=gmesh->NElements();//must keep here
    … …  
    468468        this->index2sid.clear(); this->index2sid.resize(gmesh->NElements());
    469469        this->sid2index.clear();
    470        
     470
    471471        /*Fill in the vertex_index2sid vector with non usual index value*/
    472472        for(int i=0;i<gmesh->NNodes();i++) vertex_index2sid[i]=-1;
    473        
     473
    474474        /*Fill in the this->index2sid vector with non usual index value*/
    475475        for(int i=0;i<gmesh->NElements();i++) this->index2sid[i]=-1;
    476        
     476
    477477        /*Get elements without sons and fill in the vertex_index2sid with used vertices (indexes) */
    478478        sid=0;
    … …  
    511511                newmeshXY[2*sid+1]      = coords[1]; // Y
    512512        }
    513                
     513
    514514        for(int i=0;i<this->sid2index.size();i++){//over the sid (fill the ISSM elements)
    515515                for(int j=0;j<this->GetNumberOfNodes();j++) {
    … …  
    533533
    534534                detJ=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya);
    535        
     535
    536536                /*verify and swap, if necessary*/
    537537                if(detJ<0) {
    … …  
    545545        *pxy               = newmeshXY;
    546546        *pelements      = newelements;
    547    
     547
    548548        /*Cleanup*/
    549549        xDelete<long>(vertex_index2sid);
    … …  
    554554        /* IMPORTANT! elements come in Matlab indexing
    555555                NEOPZ works only in C indexing*/
    556        
     556
    557557        if(nvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n");
    558558   if(nelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n");
    … …  
    561561    /*Verify and creating initial mesh*/
    562562   if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!");
    563    
     563
    564564   this->fathermesh = new TPZGeoMesh();
    565565        this->fathermesh->NodeVec().Resize(nvertices);
    … …  
    576576                this->fathermesh->NodeVec()[i].SetNodeId(i);
    577577        }
    578        
     578
    579579        /*Generate the elements*/
    580580        long index;
    … …  
    604604/*}}}*/
    605605TPZGeoMesh* AdaptiveMeshRefinement::CreateRefPatternMesh(TPZGeoMesh* gmesh){/*{{{*/
    606        
     606
    607607        TPZGeoMesh *newgmesh = new TPZGeoMesh();
    608608   newgmesh->CleanUp();
    609    
     609
    610610   int nnodes  = gmesh->NNodes();
    611611        int nelem   = gmesh->NElements();
    … …  
    617617        newgmesh->NodeVec().Resize(nnodes);
    618618   for(int i=0;i<nnodes;i++) newgmesh->NodeVec()[i] = gmesh->NodeVec()[i];
    619    
     619
    620620   //elements
    621621   for(int i=0;i<nelem;i++){
    622622        TPZGeoEl * geoel = gmesh->Element(i);
    623                
     623
    624624                if(!geoel){
    625625                        index=newgmesh->ElementVec().AllocateNewElement();
    … …  
    627627                        continue;
    628628                }
    629      
     629
    630630                TPZManVector<long> elem(3,0);
    631631      for(int j=0;j<3;j++) elem[j] = geoel->NodeIndex(j);
    632      
     632
    633633      newgmesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);
    634634                newgmesh->ElementVec()[index]->SetId(geoel->Id());
    635        
     635
    636636      TPZGeoElRefPattern<TPZGeoTriangle>* newgeoel = dynamic_cast<TPZGeoElRefPattern<TPZGeoTriangle>*>(newgmesh->ElementVec()[index]);
    637        
     637
    638638      //old neighbourhood
    639639      const int nsides = TPZGeoTriangle::NSides;
    … …  
    654654         }
    655655      }
    656        
     656
    657657      //inserting in new element
    658658      for(int s = 0; s < nsides; s++){
    … …  
    668668         tempEl->SetNeighbour(byside, tempSide);
    669669      }
    670        
     670
    671671      long fatherindex = geoel->FatherIndex();
    672672      if(fatherindex>-1) newgeoel->SetFather(fatherindex);
    673        
     673
    674674      if(!geoel->HasSubElement()) continue;
    675        
     675
    676676      int nsons = geoel->NSubElements();
    677677
    678678      TPZAutoPointer<TPZRefPattern> ref = gRefDBase.GetUniformRefPattern(ETriangle);
    679679      newgeoel->SetRefPattern(ref);
    680        
     680
    681681      for(int j=0;j<nsons;j++){
    682682        TPZGeoEl* son = geoel->SubElement(j);
    … …  
    690690        /*Now, build connectivities*/
    691691        newgmesh->BuildConnectivity();
    692    
     692
    693693        return newgmesh;
    694694}
    … …  
    705705/*}}}*/
    706706void AdaptiveMeshRefinement::PrintGMeshVTK(TPZGeoMesh* gmesh,std::ofstream &file,bool matColor){/*{{{*/
    707    
     707
    708708        file.clear();
    709709        long nelements = gmesh->NElements();
    … …  
    731731          MElementType elt = gel->Type();
    732732          int elNnodes = MElementType_NNodes(elt);
    733          
     733
    734734          size += (1+elNnodes);
    735735          connectivity << elNnodes;
    736          
     736
    737737          for(int t = 0; t < elNnodes; t++)
    738738          {
    … …  
    743743                        }
    744744                        node << std::endl;
    745                        
     745
    746746                        actualNode++;
    747747                        connectivity << " " << actualNode;
    748748          }
    749749          connectivity << std::endl;
    750          
     750
    751751          int elType = this->GetVTK_ElType(gel);
    752752          type << elType << std::endl;
    753          
     753
    754754          if(matColor == true)
    755755          {
    … …  
    760760                        material << gel->Index() << std::endl;
    761761          }
    762          
     762
    763763          nVALIDelements++;
    764764        }
    … …  
    790790/*}}}*/
    791791int AdaptiveMeshRefinement::GetVTK_ElType(TPZGeoEl * gel){/*{{{*/
    792    
     792
    793793        MElementType pzElType = gel->Type();
    794    
     794
    795795    int elType = -1;
    796796    switch (pzElType)
    … …  
    849849        DebugStop();
    850850    }
    851    
     851
    852852    return elType;
    853853}
    854854/*}}}*/
    855 

  • issm/trunk-jpl/src/c/classes/AmrBamg.cpp

    r23065 r23066  
    3434        this->deviatoricerror_groupthreshold    = -1;
    3535        this->deviatoricerror_maximum                           = -1;
    36        
     36
    3737        /*Geometry and mesh as NULL*/
    3838        this->geometry                                                                  = NULL;
    … …  
    137137        this->options->hmaxVerticesSize[0] = 0;
    138138        this->options->hmaxVerticesSize[1] = 0;
    139        
     139
    140140        /*verify if the metric will be reseted or not*/
    141141        if(this->keepmetric==0){
    … …  
    155155        IssmDouble *xylist= xNew<IssmDouble>(nbv*2);
    156156        int* elementslist = xNew<int>(nbt*3);
    157        
     157
    158158        datalist[0] = nbv;
    159159        datalist[1] = nbt;
    160        
     160
    161161        for(int i=0;i<nbv;i++){
    162162                xylist[2*i]             = meshout->Vertices[i*3+0];
    163163                xylist[2*i+1]   = meshout->Vertices[i*3+1];
    164164        }
    165        
     165
    166166        for(int i=0;i<nbt;i++){
    167167                elementslist[3*i+0] = reCast<int>(meshout->Triangles[4*i+0]);
    … …  
    174174        *pxylist                        = xylist;
    175175        *pelementslist = elementslist;
    176        
     176
    177177        /*Cleanup and return*/
    178178        delete geomout;
    … …  
    181181
    182182        if(!this->options) _error_("AmrBamg->options is NULL!");
    183        
     183
    184184        this->options->hmin     = hmin_in;
    185185        this->options->hmax     = hmax_in;

  • issm/trunk-jpl/src/c/classes/Cfdragcoeffabsgrad.cpp

    r22971 r23066  
    2323#include "../classes/gauss/Gauss.h"
    2424/*}}}*/
    25                
     25
    2626/*Cfdragcoeffabsgrad constructors, destructors :*/
    2727Cfdragcoeffabsgrad::Cfdragcoeffabsgrad(){/*{{{*/
    … …  
    4040
    4141        this->definitionenum=in_definitionenum;
    42        
     42
    4343        this->name              = xNew<char>(strlen(in_name)+1);
    4444        xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
    … …  
    4646        this->weights_enum=in_weights_enum;
    4747        this->timepassedflag=in_timepassedflag;
    48        
     48
    4949        this->misfit=0;
    5050        this->lock=0;
    … …  
    102102         /*diverse: */
    103103         IssmDouble time;
    104          
     104
    105105         /*recover time parameters: */
    106106                 int i;
    … …  
    116116                 ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    117117                 J=J_sum;
    118                
     118
    119119                 timepassedflag = true;
    120120                 return J;
    … …  
    183183        return Jelem;
    184184}/*}}}*/
    185 

  • issm/trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp

    r22612 r23066  
    2323#include "../classes/gauss/Gauss.h"
    2424/*}}}*/
    25                
     25
    2626/*Cfsurfacelogvel constructors, destructors :*/
    2727Cfsurfacelogvel::Cfsurfacelogvel(){/*{{{*/
    … …  
    4040
    4141        this->definitionenum=in_definitionenum;
    42        
     42
    4343        this->name              = xNew<char>(strlen(in_name)+1);
    4444        xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
    … …  
    4646        this->datatime=in_datatime;
    4747        this->timepassedflag=in_timepassedflag;
    48        
     48
    4949        this->misfit=0;
    5050        this->lock=0;
    … …  
    100100/*}}}*/
    101101IssmDouble Cfsurfacelogvel::Response(FemModel* femmodel){/*{{{*/
    102                  
     102
    103103         /*diverse: */
    104104         IssmDouble time;
    105          
     105
    106106         /*recover time parameters: */
    107107         femmodel->parameters->FindParam(&time,TimeEnum);
    … …  
    112112                 IssmDouble J=0.;
    113113                 IssmDouble J_sum=0.;
    114        
     114
    115115         if(datatime<=time && !timepassedflag){
    116116                 for(i=0;i<femmodel->elements->Size();i++){
    … …  
    122122                 ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    123123                 J=J_sum;
    124                
     124
    125125                 timepassedflag = true;
    126126                 return J;
    … …  
    139139        IssmDouble vx,vy,vxobs,vyobs,weight;
    140140        IssmDouble* xyz_list = NULL;
    141        
     141
    142142        /*Get basal element*/
    143143        if(!element->IsOnSurface()) return 0.;
    … …  
    160160        /* Get node coordinates*/
    161161        topelement->GetVerticesCoordinates(&xyz_list);
    162        
     162
    163163        /*Get model values*/
    164164        Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
    … …  
    174174        if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do");
    175175        datasetinput = (DatasetInput*)tempinput;
    176 
    177176
    178177        /* Start  looping on the number of gaussian points: */
    … …  
    222221        return Jelem;
    223222}/*}}}*/
    224 

  • issm/trunk-jpl/src/c/classes/Cfsurfacesquare.cpp

    r22612 r23066  
    2323#include "../classes/gauss/Gauss.h"
    2424/*}}}*/
    25                
     25
    2626/*Cfsurfacesquare constructors, destructors :*/
    2727Cfsurfacesquare::Cfsurfacesquare(){/*{{{*/
    … …  
    4343
    4444        this->definitionenum=in_definitionenum;
    45        
     45
    4646        this->name              = xNew<char>(strlen(in_name)+1);
    4747        xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
    … …  
    5252        this->datatime=in_datatime;
    5353        this->timepassedflag=in_timepassedflag;
    54        
     54
    5555        this->misfit=0;
    5656        this->lock=0;
    … …  
    111111         /*diverse: */
    112112         IssmDouble time;
    113          
     113
    114114         /*recover time parameters: */
    115115         femmodel->parameters->FindParam(&time,TimeEnum);
    … …  
    130130                 ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    131131                 J=J_sum;
    132                
     132
    133133                 timepassedflag = true;
    134134                 return J;
    … …  
    172172        /*Get input if it already exists*/
    173173        Input*  tempinput = topelement->GetInput(definitionenum);
    174        
     174
    175175        /*Cast it to a Datasetinput*/
    176176        if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do");
    … …  
    214214        return Jelem;
    215215}/*}}}*/
    216 

  • issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp

    r22361 r23066  
    3737/*Object virtual functions definitions:*/
    3838Object* SpcStatic::copy() {/*{{{*/
    39        
     39
    4040        SpcStatic* spcstat = new SpcStatic(*this);
    4141

  • issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp

    r22617 r23066  
    2020                int world_rank;
    2121                ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&world_rank);
    22                
     22
    2323                /*Build an femmodel if you are a slave, using the corresponding communicator:*/
    2424                if(world_rank!=0){
    … …  
    3333                int world_rank;
    3434                ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&world_rank);
    35                
     35
    3636                if(world_rank!=0){
    3737

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

    r23059 r23066  
    11051105void       Element::GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug){/*{{{*/
    11061106
    1107 
    11081107        /*Get number of nodes for this element*/
    11091108        int numnodes = this->GetNumberOfNodes();
    … …  
    11211120                if(ug[nodes[i]->GetDof(0,GsetEnum)] > input_max) input_max = ug[nodes[i]->GetDof(0,GsetEnum)];
    11221121        }
    1123 
    11241122
    11251123        /*Second loop to reassign min and max with local extrema*/
    … …  
    13891387                _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
    13901388        }
    1391        
     1389
    13921390        /*Clean up*/
    13931391        xDelete<int>(doflist);
    … …  
    14761474        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    14771475        parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
    1478 
    14791476
    14801477        /*Get number of vertices*/
    … …  
    18081805        }
    18091806
    1810 
    18111807    /*Branch on type of vector: nodal or elementary: */
    18121808    if(vector_type==1){ //nodal vector
    … …  
    20672063        this->GetInputListOnVertices(deepwaterel,BasalforcingsDeepwaterElevationEnum);
    20682064        this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum);
    2069        
     2065
    20702066        for(int i=0;i<numvertices;i++){
    20712067                if(base[i]>upperwaterel[i])      values[i]=0;
    … …  
    28702866        if (!IsOnSurface()) return;
    28712867
    2872 
    28732868        /*Retrieve material properties and parameters:{{{ */
    28742869        rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    … …  
    30723067                if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
    30733068
    3074 
    30753069                /*Distribution of absorbed short wave radation with depth:*/
    30763070                if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());

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

    r22004 r23066  
    104104
    105105        if(marshall_direction==MARSHALLING_BACKWARD){
    106                
     106
    107107                if (!hnodes_null)this->hnodes = new Hook*[numanalyses];
    108108                else this->hnodes=NULL;
    … …  
    140140        _printf_(" ElementHook DeepEcho:\n");
    141141        _printf_("  numanalyses : "<< this->numanalyses <<"\n");
    142        
     142
    143143        _printf_("  hnodes:\n");
    144144        if(hnodes){
    … …  
    149149        }
    150150        else _printf_("  hnodes = NULL\n");
    151        
     151
    152152        _printf_("  hvertices:\n");
    153153        if(hvertices) hvertices->DeepEcho();
    154154        else _printf_("  hvertices = NULL\n");
    155        
     155
    156156        _printf_("  hmaterial:\n");
    157157        if(hmaterial) hmaterial->DeepEcho();
    … …  
    170170/*}}}*/
    171171void ElementHook::Echo(){/*{{{*/
    172        
     172
    173173        _printf_(" ElementHook Echo:\n");
    174174        _printf_("  numanalyses : "<< this->numanalyses <<"\n");
    175        
     175
    176176        _printf_("  hnodes:\n");
    177177        if(hnodes){
    … …  
    181181        }
    182182        else _printf_("  hnodes = NULL\n");
    183        
     183
    184184        _printf_("  hvertices:\n");
    185185        if(hvertices) hvertices->Echo();
    186186        else _printf_("  hvertices = NULL\n");
    187        
     187
    188188        _printf_("  hmaterial:\n");
    189189        if(hmaterial) hmaterial->Echo();

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

    r23059 r23066  
    286286        IssmDouble  calvingratey[NUMVERTICES];
    287287        IssmDouble  calvingrate[NUMVERTICES];
    288 
    289288
    290289        /* Get node coordinates and dof list: */
    … …  
    10481047/*}}}*/
    10491048void       Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    1050        
     1049
    10511050        /* Intermediaries */
    10521051        const int dim=3;
    … …  
    11811180        if(!IsInputEnum(control_enum)) _error_("Enum "<<EnumToStringx(control_enum)<<" is not in IsInput");
    11821181        Input* input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
    1183 
    11841182
    11851183        /*Cast to Controlinput*/
    … …  
    25952593        Penta* penta=this;
    25962594        for(;;){
    2597        
     2595
    25982596                IssmDouble  xyz_list[NUMVERTICES][3];
    25992597                /* Get node coordinates and dof list: */
    … …  
    26042602                Jdet[1]=(xyz_list[4][2]-xyz_list[1][2])*0.5;
    26052603                Jdet[2]=(xyz_list[5][2]-xyz_list[2][2])*0.5;
    2606        
     2604
    26072605                /*Retrieve all inputs we will need*/
    26082606                Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    … …  
    26152613                Input* surface_input=inputs->GetInput(SurfaceEnum);                                                             _assert_(surface_input);
    26162614                Input* thickness_input=inputs->GetInput(ThicknessEnum);                                                 _assert_(thickness_input);
    2617                
     2615
    26182616                /* Start looping on the number of 2D vertices: */
    26192617                for(int ig=0;ig<3;ig++){
    … …  
    26442642                        delete gauss;
    26452643                }
    2646                        
     2644
    26472645                /*Stop if we have reached the surface/base*/
    26482646                if(penta->IsOnSurface()) break;
    2649                
     2647
    26502648                /*get upper Penta*/
    26512649                penta=penta->GetUpperPenta();

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

    r22647 r23066  
    104104        return seg;
    105105
    106 
    107106}
    108107/*}}}*/
    … …  
    141140/*}}}*/
    142141void       Seg::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    143        
     142
    144143        /* Intermediaries */
    145144        int nrfrontnodes,index;

  • 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));

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r23065 r23066  
    242242        char *lockfilename   = NULL;
    243243        bool  waitonlock     = false;
    244 
    245244
    246245        /*Write lock file if requested: */
    … …  
    43934392        xDelete<IssmDouble>(RSLg_old);
    43944393
    4395 
    43964394}
    43974395/*}}}*/
    … …  
    48074805        size_t* poutputbuffersize;
    48084806
    4809 
    48104807        /*Before we delete the profiler, report statistics for this run: */
    48114808        profiler->Stop(TOTAL);  //final tagging
    … …  
    48684865}/*}}}*/
    48694866#endif
    4870 
    48714867
    48724868#if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_)
    … …  
    52465242      newy[i] = newxylist[2*i+1];
    52475243   }
    5248        
     5244
    52495245        /*Assign the pointers*/
    52505246   (*pnewelementslist)  = newelementslist; //Matlab indexing

  • issm/trunk-jpl/src/c/classes/IoModel.cpp

    r22993 r23066  
    157157        bool autodiff=false;
    158158        bool iscontrol=false;
    159        
     159
    160160        /*First, keep track of the file handle: */
    161161        this->fid=iomodel_handle;
    … …  
    422422        this->FetchData(&autodiff,"md.autodiff.isautodiff");
    423423        this->FetchData(&iscontrol,"md.inversion.iscontrol");
    424        
     424
    425425        if(trace || (autodiff && !iscontrol)){
    426426
    … …  
    506506
    507507        char** stringarray=*pstringarray;
    508        
     508
    509509        if(numstrings){
    510510                for(int i=0;i<numstrings;i++){
    … …  
    597597                                        /*Read the integer and broadcast it to other cpus:*/
    598598                                        if(fread(&integer,sizeof(int),1,this->fid)!=1) _error_("could not read integer ");
    599                                        
     599
    600600                                        /*Convert codes to Enums if needed*/
    601601                                        if(strcmp(record_name,"md.smb.model")==0) integer = IoCodeToEnumSMB(integer);
    … …  
    963963        /*recover my_rank:*/
    964964        int my_rank=IssmComm::GetRank();
    965        
     965
    966966        /*Set file pointer to beginning of the data: */
    967967        fid=this->SetFilePointerToData(&code,NULL,data_name);
    … …  
    11341134                }
    11351135        }
    1136          
     1136
    11371137        /*output: */
    11381138        int          M,N;
    … …  
    18251825                                /*check we are indeed finding a string, not something else: */
    18261826                                if(codes[i]!=4)_error_("expecting a string for \""<<data_name<<"\"");
    1827                
     1827
    18281828                                /*We have to read a string from disk. First read the dimensions of the string, then the string: */
    18291829                                fsetpos(fid,file_positions+i);
    … …  
    18541854        xDelete<int>(codes);
    18551855        xDelete<fpos_t>(file_positions);
    1856        
     1856
    18571857        /*Assign output pointers: */
    18581858        *pstrings=strings;
    … …  
    18751875        /*recover my_rank:*/
    18761876        int my_rank=IssmComm::GetRank();
    1877        
     1877
    18781878        /*Get file pointers to beginning of the data (multiple instances of it): */
    18791879        file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name);
    … …  
    18901890
    18911891                                if(code!=2)_error_("expecting an integer for \""<<data_name<<"\"");
    1892                                
     1892
    18931893                                /*We have to read a integer from disk. First read the dimensions of the integer, then the integer: */
    18941894                                fsetpos(fid,file_positions+i);
    … …  
    19031903                }
    19041904        }
    1905                        
     1905
    19061906        /*Free ressources:*/
    19071907        xDelete<fpos_t>(file_positions);
    … …  
    19281928        /*recover my_rank:*/
    19291929        int my_rank=IssmComm::GetRank();
    1930        
     1930
    19311931        /*Get file pointers to beginning of the data (multiple instances of it): */
    19321932        file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name);
    … …  
    19431943
    19441944                                if(code!=3)_error_("expecting a double for \""<<data_name<<"\"");
    1945                                
     1945
    19461946                                /*We have to read a double from disk: */
    19471947                                fsetpos(fid,file_positions+i);
    … …  
    19561956                }
    19571957        }
    1958                        
     1958
    19591959        /*Free ressources:*/
    19601960        xDelete<fpos_t>(file_positions);
    … …  
    19851985        /*recover my_rank:*/
    19861986        int my_rank=IssmComm::GetRank();
    1987        
     1987
    19881988        /*Get file pointers to beginning of the data (multiple instances of it): */
    19891989        file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name);
    … …  
    20142014                        }
    20152015                        ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    2016 
    20172016
    20182017                        /*Now allocate matrix: */
    … …  
    20392038                        else
    20402039                                matrix=NULL;
    2041                        
    2042                        
     2040
    20432041                        /*Assign: */
    20442042                        mdims[i]=M;
    … …  
    20472045                }
    20482046        }
    2049                        
     2047
    20502048        /*Free ressources:*/
    20512049        xDelete<fpos_t>(file_positions);
    … …  
    20892087        /*recover my_rank:*/
    20902088        int my_rank=IssmComm::GetRank();
    2091        
     2089
    20922090        /*Get file pointers to beginning of the data (multiple instances of it): */
    20932091        file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name);
    … …  
    21182116                        }
    21192117                        ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    2120 
    21212118
    21222119                        /*Now allocate matrix: */
    … …  
    21442141                        else
    21452142                                integer_matrix=NULL;
    2146                        
    2147                        
     2143
    21482144                        /*Assign: */
    21492145                        mdims[i]=M;
    … …  
    21522148                }
    21532149        }
    2154                        
     2150
    21552151        /*Free ressources:*/
    21562152        xDelete<fpos_t>(file_positions);
    … …  
    23762372                        vector_types   = xNew<int>(num_instances);
    23772373                }
    2378        
     2374
    23792375                /*Reset FILE* position to the beginning of the file, and start again, this time saving the data information
    23802376                 * as we find it: */
    … …  
    24192415                                vector_types[counter] = vector_type;
    24202416                                fgetpos(fid,file_positions+counter);
    2421                                
     2417
    24222418                                /*backup and skip over the record, as we have more work to do: */
    24232419                                if(5<=record_code && record_code<=7) fseek(fid,-sizeof(int),SEEK_CUR);
    24242420                                fseek(fid,-sizeof(int),SEEK_CUR);
    24252421                                fseek(fid,-sizeof(int),SEEK_CUR);
    2426                                
     2422
    24272423                                /*increment counter: */
    24282424                                counter++;

  • issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

    r23020 r23066  
    196196        if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
    197197        else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
    198        
     198
    199199        _assert_(!xIsNan<IssmDouble>(alpha_complement));
    200200        _assert_(!xIsInf<IssmDouble>(alpha_complement));
    … …  
    270270        r=drag_q/drag_p;
    271271        s=1./drag_p;
    272 
    273272
    274273        /*Get effective pressure*/

  • issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp

    r23020 r23066  
    178178
    179179        switch(analysis_type){
    180                
     180
    181181        case HydrologyShaktiAnalysisEnum:
    182182                pe = this->CreatePVectorHydrologyShakti();
    … …  
    395395        IssmDouble moulin_load,dt;
    396396        ElementVector* pe=new ElementVector(&node,1,this->parameters);
    397        
     397
    398398        this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
    399399        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    400        
     400
    401401        pe->values[0]=moulin_load*dt;
    402402        /*Clean up and return*/

  • issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp

    r23020 r23066  
    3030/*}}}*/
    3131Neumannflux::Neumannflux(int neumannflux_id,int i,IoModel* iomodel,int* segments,int in_analysis_type){/*{{{*/
    32 
    3332
    3433        /*Some sanity checks*/

  • issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp

    r22306 r23066  
    109109
    110110        if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook();
    111        
     111
    112112        MARSHALLING_ENUM(MatestarEnum);
    113113        MARSHALLING(mid);

  • issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

    r22307 r23066  
    135135        //if(helement) helement->DeepEcho();
    136136        //else _printf_("   helement = NULL\n");
    137        
     137
    138138        _printf_("   element:\n");
    139139        _printf_("     note: element not printed to avoid recursion.\n");
    … …  
    143143/*}}}*/
    144144void      Matice::Echo(void){/*{{{*/
    145        
     145
    146146        _printf_("Matice:\n");
    147147        _printf_("   mid: " << mid << "\n");
    148148        _printf_("   isdamaged: " << isdamaged << "\n");
    149149        _printf_("   isenhanced: " << isenhanced << "\n");
    150        
     150
    151151        /*helement and element Echo were commented to avoid recursion.*/
    152152        /*Example: element->Echo calls matice->Echo which calls element->Echo etc*/
    … …  
    155155        //if(helement) helement->Echo();
    156156        //else _printf_("   helement = NULL\n");
    157        
     157
    158158        _printf_("   element:\n");
    159159        _printf_("     note: element not printed to avoid recursion.\n");
    … …  
    167167
    168168        if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook();
    169        
     169
    170170        MARSHALLING_ENUM(MaticeEnum);
    171171        MARSHALLING(mid);
    … …  
    541541        IssmDouble exx,eyy,exy,exz,eyz;
    542542
    543 
    544543        if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
    545544                                (epsilon[3]==0) && (epsilon[4]==0)){

  • issm/trunk-jpl/src/c/classes/Materials/Matlitho.cpp

    r22004 r23066  
    3737        this->radius=xNew<IssmDouble>(this->numlayers+1);
    3838        xMemCpy<IssmDouble>(this->radius, iomodel->Data("md.materials.radius"),this->numlayers+1);
    39        
     39
    4040        this->viscosity=xNew<IssmDouble>(this->numlayers);
    4141        xMemCpy<IssmDouble>(this->viscosity, iomodel->Data("md.materials.viscosity"),this->numlayers);
    42        
     42
    4343        this->lame_lambda=xNew<IssmDouble>(this->numlayers);
    4444        xMemCpy<IssmDouble>(this->lame_lambda, iomodel->Data("md.materials.lame_lambda"),this->numlayers);
    45        
     45
    4646        this->lame_mu=xNew<IssmDouble>(this->numlayers);
    4747        xMemCpy<IssmDouble>(this->lame_mu, iomodel->Data("md.materials.lame_mu"),this->numlayers);
    … …  
    6161        this->issolid=xNew<IssmDouble>(this->numlayers);
    6262        xMemCpy<IssmDouble>(this->issolid, iomodel->Data("md.materials.issolid"),this->numlayers);
    63        
     63
    6464        /*isburgersd= xNew<IssmDouble>(this->numlayers);
    6565        this->isburgers=xNew<bool>(this->numlayers);
    6666        xMemCpy<IssmDouble>(isburgersd, iomodel->Data("md.materials.isburgers"),this->numlayers);
    6767        for (int i=0;i<this->numlayers;i++)this->isburgers[i]=reCast<bool,IssmDouble>(isburgersd[i]);
    68        
     68
    6969        issolidd= xNew<IssmDouble>(this->numlayers);
    7070        this->issolid=xNew<bool>(this->numlayers);
    7171        xMemCpy<IssmDouble>(issolidd, iomodel->Data("md.materials.issolid"),this->numlayers);
    7272        for (int i=0;i<this->numlayers;i++)this->issolid[i]=reCast<bool,IssmDouble>(issolidd[i]);*/
    73        
     73
    7474        /*free ressources: */
    7575        xDelete<IssmDouble>(isburgersd);
    … …  
    7878/*}}}*/
    7979Matlitho::~Matlitho(){/*{{{*/
    80        
     80
    8181        xDelete<IssmDouble>(radius);
    8282        xDelete<IssmDouble>(viscosity);

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

    r23020 r23066  
    5959        mantle_shear_modulus      = 0;
    6060        mantle_density            = 0;
    61        
     61
    6262        earth_density             = 0;
    6363
    … …  
    338338        matpar->mantle_shear_modulus=this->mantle_shear_modulus;
    339339        matpar->mantle_density=this->mantle_density;
    340        
     340
    341341        matpar->earth_density=this->earth_density;
    342342
    … …  
    439439        MARSHALLING(mantle_shear_modulus);
    440440        MARSHALLING(mantle_density);
    441        
     441
    442442        //slr:
    443443        MARSHALLING(earth_density);

  • issm/trunk-jpl/src/c/classes/Misfit.cpp

    r22438 r23066  
    2323#include "../classes/gauss/Gauss.h"
    2424/*}}}*/
    25                
     25
    2626/*Misfit constructors, destructors :*/
    2727Misfit::Misfit(){/*{{{*/
    … …  
    4242
    4343        this->definitionenum=in_definitionenum;
    44        
     44
    4545        this->name              = xNew<char>(strlen(in_name)+1);
    4646        xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
    … …  
    4848        this->timeinterpolation = xNew<char>(strlen(in_timeinterpolation)+1);
    4949        xMemCpy<char>(this->timeinterpolation,in_timeinterpolation,strlen(in_timeinterpolation)+1);
    50                                
     50
    5151        this->model_enum=in_model_enum;
    5252        this->observation_enum=in_observation_enum;
    5353        this->weights_enum=in_weights_enum;
    5454        this->local=in_local;
    55        
     55
    5656        this->misfit=0;
    5757        this->lock=0;
    … …  
    111111/*}}}*/
    112112IssmDouble Misfit::Response(FemModel* femmodel){/*{{{*/
    113                  
     113
    114114         /*diverse: */
    115115         IssmDouble time,starttime,finaltime;
    116116         IssmDouble dt;
    117          
     117
    118118         /*recover time parameters: */
    119119         femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    … …  
    130130                 IssmDouble all_area_t;
    131131
    132        
    133132                 /*If we are locked, return time average: */
    134133                 if(this->lock)return misfit/(time-starttime);
    … …  
    144143                 area_t=all_area_t;
    145144                 misfit_t=all_misfit_t;
    146                  
     145
    147146                 /*Divide by surface area if not nill!: */
    148147                 if (area_t!=0) misfit_t=misfit_t/area_t;
    … …  
    164163                 IssmDouble* weights= NULL;
    165164                 int msize,osize,wsize;
    166                  
     165
    167166                 /*Are we transient?:*/
    168167                 if (time==0){
    169168                         IssmDouble misfit_t=0.;
    170                          
     169
    171170                         /*get global vectors: */
    172171                         GetVectorFromInputsx(&model,&msize,femmodel,model_enum);
    … …  
    190189                 }
    191190                 else{
    192                          
     191
    193192                         IssmDouble misfit_t=0.;
    194193                         IssmDouble all_misfit_t=0.;
    … …  
    201200                         GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize);
    202201                         GetVectorFromInputsx(&weights,&wsize,femmodel,weights_enum); _assert_(wsize==msize);
    203                          
     202
    204203                         int count=0;
    205204                         for (int i=0;i<msize;i++){
    … …  
    215214                         /*Do we lock? i.e. are we at final_time? :*/
    216215                         if(time==finaltime)this->lock=1;
    217                          
     216
    218217                         /*Free ressources:*/
    219218                         xDelete<IssmDouble>(model);
    … …  
    227226         } /*}}}*/
    228227         else{ /*global computation: {{{ */
    229                  
     228
    230229                 IssmDouble model, observation;
    231                  
     230
    232231                 /*If we are locked, return time average: */
    233232                 if(this->lock)return misfit/(time-starttime);
    … …  
    239238                 Input* input = element->GetInput(observation_enum); _assert_(input);
    240239                 input->GetInputAverage(&observation);
    241                  
     240
    242241                 /*Add this time's contribution to curent misfit: */
    243242                 misfit+=dt*(model-observation);
    244                  
     243
    245244                 /*Do we lock? i.e. are we at final_time? :*/
    246245                 if(time==finaltime)this->lock=1;
    247                  
     246
    248247                 /*What we return is the value of misfit / time: */
    249248                 return misfit/(time-starttime);

  • issm/trunk-jpl/src/c/classes/Nodalvalue.cpp

    r22427 r23066  
    8787/*}}}*/
    8888IssmDouble Nodalvalue::Response(FemModel* femmodel){/*{{{*/
    89        
     89
    9090         /*output:*/
    9191         IssmDouble value;

  • issm/trunk-jpl/src/c/classes/Node.cpp

    r21509 r23066  
    498498/*}}}*/
    499499void  Node::SetApproximation(int in_approximation){/*{{{*/
    500        
     500
    501501        this->approximation = in_approximation;
    502502}

  • issm/trunk-jpl/src/c/classes/Numberedcostfunction.cpp

    r22423 r23066  
    88#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
    99#endif
    10 
    1110
    1211/*Headers:*/
    … …  
    3029
    3130/*}}}*/
    32                
     31
    3332/*Numberedcostfunction constructors, destructors :*/
    3433Numberedcostfunction::Numberedcostfunction(){/*{{{*/
    … …  
    110109/*}}}*/
    111110IssmDouble Numberedcostfunction::Response(FemModel* femmodel){/*{{{*/
    112        
     111
    113112         _assert_(number_cost_functions>0 && number_cost_functions<1e3);
    114113
    … …  
    157156 }
    158157 /*}}}*/
    159 

  • issm/trunk-jpl/src/c/classes/Params/BoolParam.cpp

    r20827 r23066  
    6363}
    6464/*}}}*/
    65                

  • issm/trunk-jpl/src/c/classes/Params/DataSetParam.cpp

    r22588 r23066  
    5555
    5656        if(marshall_direction==MARSHALLING_BACKWARD)value=new DataSet();
    57        
     57
    5858        MARSHALLING_ENUM(DataSetParamEnum);
    5959        MARSHALLING(enum_type);

  • issm/trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.cpp

    r20827 r23066  
    233233}
    234234/*}}}*/
    235                

  • issm/trunk-jpl/src/c/classes/Params/DoubleMatParam.cpp

    r20635 r23066  
    116116/*DoubleMatParam specific routines:*/
    117117void  DoubleMatParam::GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM,int* pN){/*{{{*/
    118        
     118
    119119        /*Assign output pointers:*/
    120120        if(pM) *pM=M;

  • issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp

    r22628 r23066  
    116116/*DoubleVecParam specific routines:*/
    117117void  DoubleVecParam::GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM){/*{{{*/
    118        
     118
    119119        /*Assign output pointers:*/
    120120        if(pM) *pM=M;

  • issm/trunk-jpl/src/c/classes/Params/IntParam.cpp

    r20827 r23066  
    6464}
    6565/*}}}*/
    66 

  • issm/trunk-jpl/src/c/classes/Params/StringArrayParam.cpp

    r21936 r23066  
    8383
    8484        MARSHALLING_ENUM(StringArrayParamEnum);
    85        
     85
    8686        MARSHALLING(enum_type);
    8787        MARSHALLING(numstrings);

  • issm/trunk-jpl/src/c/classes/Params/StringParam.cpp

    r21915 r23066  
    5555
    5656        if(marshall_direction==MARSHALLING_FORWARD || marshall_direction == MARSHALLING_SIZE)size=strlen(value)+1;
    57        
     57
    5858        MARSHALLING_ENUM(StringParamEnum);
    5959        MARSHALLING(enum_type);

  • issm/trunk-jpl/src/c/classes/Profiler.cpp

    r22554 r23066  
    100100        if(this->running[tag]) _error_("Tag "<<tag<<" has not been stopped");
    101101
    102 
    103102        #ifdef _HAVE_MPI_
    104103        return this->time[tag];
    … …  
    143142        _assert_(tag<MAXIMUMSIZE);
    144143        if(this->running[tag]) _error_("Tag "<<tag<<" is already running");
    145 
    146144
    147145        /*If mpisync requested, make sure all the cpus are at the same point in the execution: */
    … …  
    184182        if(!this->running[tag]) _error_("Tag "<<tag<<" is not running");
    185183
    186 
    187184        /*If mpisync requested, make sure all the cpus are at the same point in the execution: */
    188185        if(!dontmpisync){

  • issm/trunk-jpl/src/c/classes/Regionaloutput.cpp

    r22438 r23066  
    155155}
    156156/*}}}*/
    157 

  • issm/trunk-jpl/src/c/classes/gauss/GaussSeg.cpp

    r20827 r23066  
    3131        this->numgauss = order;
    3232        GaussLegendreLinear(&pcoords1,&pweights,order);
    33        
     33
    3434        this->coords1=xNew<IssmDouble>(numgauss);
    3535        this->weights=xNew<IssmDouble>(numgauss);

  • issm/trunk-jpl/src/c/classes/kriging/Observations.cpp

    r21915 r23066  
    741741        xDelete<IssmPDouble>(counter);
    742742}/*}}}*/
    743 

  • issm/trunk-jpl/src/c/cores/ad_core.cpp

    r22610 r23066  
    4646                        /*First, stop tracing: */
    4747                        trace_off();
    48                        
     48
    4949                        /*Print tape statistics so that user can kill this run if something is off already:*/
    5050                        if(VerboseAutodiff()){ /*{{{*/
    … …  
    9393                        femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
    9494                        femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
    95        
     95
    9696                        /*if no dependents, no point in running a driver: */
    9797                        if(!(num_dependents*num_independents)) return;
    … …  
    101101                                num_dependents=0; num_independents=0;
    102102                        }
    103                        
     103
    104104                        /*retrieve state variable: */
    105105                        femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum);
    … …  
    117117                        femmodel->parameters->FindParam(&driver,AutodiffDriverEnum);
    118118
    119 
    120119                        if (strcmp(driver,"fos_forward")==0){ /*{{{*/
    121120
    … …  
    140139                                anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx;
    141140#endif
    142 
    143141
    144142                                /*call driver: */
    … …  
    184182#endif
    185183                                // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
    186 
    187184
    188185                                /*seed matrix: */
    … …  
    244241#endif
    245242
    246 
    247243                                /*call driver: */
    248244                                fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
    … …  
    284280                                anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
    285281                                #endif
    286 
    287282
    288283                                /*seed matrix: */
    … …  
    317312                        else _error_("driver: " << driver << " not yet supported!");
    318313
    319 
    320314                        if(VerboseAutodiff())_printf0_("   end AD core\n");
    321315

  • issm/trunk-jpl/src/c/cores/adgradient_core.cpp

    r18940 r23066  
    3939        IssmPDouble *aWeightVector=NULL;
    4040        IssmPDouble *weightVectorTimesJac=NULL;
    41        
     41
    4242        /*output: */
    4343        IssmPDouble *totalgradient=NULL;
    … …  
    5353                        /*First, stop tracing: */
    5454                        trace_off();
    55                        
     55
    5656                        if(VerboseAutodiff()){ /*{{{*/
    5757                                size_t  tape_stats[15];
    … …  
    9696                        delete [] sstats;
    9797                } /*}}}*/
    98                
     98
    9999                        /*retrieve parameters: */
    100100                        femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
    101101                        femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
    102                        
     102
    103103                        /*if no dependents, no point in running a driver: */
    104104                        if(!(num_dependents*num_independents)) return;
    … …  
    109109                                num_dependents=0; num_independents=0;
    110110                        }
    111                        
     111
    112112                        /*retrieve state variable: */
    113113                        femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum);
    114                        
     114
    115115                        /* driver argument */
    116116                        xp=xNew<double>(num_independents);
    … …  
    131131                        /*Initialize outputs: */
    132132                        totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
    133                        
     133
    134134                        for(aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
    135135
    … …  
    137137                                aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
    138138                                if (my_rank==0) aWeightVector[aDepIndex]=1.0;
    139                                
     139
    140140                                /*initialize output gradient: */
    141141                                weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
    … …  
    162162                                xDelete(aWeightVector);
    163163                        }
    164                
     164
    165165                        /*add totalgradient to results*/
    166166                        femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,totalgradient,num_independents,1,1,0.0));
    167167
    168168                        if(VerboseAutodiff())_printf0_("   end ad core\n");
    169                        
     169
    170170                        /* delete the allocated space for the parameters and free ressources:{{{*/
    171171                        xDelete(anEDF_for_solverx_p->dp_x);

  • issm/trunk-jpl/src/c/cores/bmb_core.cpp

    r22974 r23066  
    1111
    1212void bmb_core(FemModel* femmodel){
    13 
    1413
    1514        /*First, get BMB model from parameters*/

  • issm/trunk-jpl/src/c/cores/controlad_core.cpp

    r20803 r23066  
    105105                default: _printf0_("   Unknown end condition\n");
    106106        }
    107        
     107
    108108        /*Save results:*/
    109109        femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,G,n,1,0,0.0));
    … …  
    130130        IssmPDouble*   X=NULL;
    131131        int            i;
    132        
     132
    133133        /*Now things get complicated. The femmodel we recovered did not initialize an AD trace, so we can't compute gradients with it. We are going to recreate
    134134         *a new femmodel, identical in all aspects to the first one, with trace on though, which will allow us to run the forward mode and get the gradient
    … …  
    143143        femmodel=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, restartfilename,IssmComm::GetComm(), femmodel->solution_type,NULL);
    144144
    145        
    146145        /*Get initial guess:*/
    147146        femmodel->parameters->FindParam(&Xd,&intn,AutodiffXpEnum);
    … …  
    175174        IssmDouble    pfd;
    176175        int            i;
    177        
     176
    178177        /*Recover Femmodel*/
    179178        femmodel  = (FemModel*)dzs;
    … …  
    195194        femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast<IssmPDouble>(pfd);
    196195        _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
    197        
     196
    198197        /*Compute gradient using AD. Gradient is in the results after the ad_core is called*/
    199198        adgradient_core(femmodel);
    … …  
    211210        /*MPI broadcast results:*/
    212211        ISSM_MPI_Bcast(G2,*n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
    213        
     212
    214213        /*Send gradient to m1qn3 core: */
    215214        for(long i=0;i<*n;i++) G[i] = G2[i];
    216        
     215
    217216        /*Constrain X and G*/
    218217        IssmDouble  Gnorm = 0.;
    … …  
    228227        xDelete<IssmDouble>(Jlist);
    229228        xDelete<IssmPDouble>(G2);
    230        
     229
    231230        xDelete<char>(rootpath);
    232231        xDelete<char>(inputfilename);
    … …  
    251250        IssmDouble    pfd;
    252251        int            i;
    253        
     252
    254253        /*Recover Femmodel*/
    255254        femmodel  = (FemModel*)dzs;
    … …  
    271270        femmodelad=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, restartfilename,IssmComm::GetComm(), femmodel->solution_type,X);
    272271        femmodel=femmodelad; //We can do this, because femmodel is being called from outside, not by reference, so we won't erase it
    273        
     272
    274273        /*Recover some parameters*/
    275274        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    … …  
    284283        femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast<IssmPDouble>(pfd);
    285284        _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
    286        
     285
    287286        /*Compute gradient using AD. Gradient is in the results after the ad_core is called*/
    288287        adgradient_core(femmodel);
    … …  
    300299        /*MPI broadcast results:*/
    301300        ISSM_MPI_Bcast(G2,*n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
    302        
     301
    303302        /*Send gradient to m1qn3 core: */
    304303        for(long i=0;i<*n;i++) G[i] = G2[i];
    305        
     304
    306305        /*Recover Gnorm: */
    307306        IssmDouble  Gnorm = 0.;
    … …  
    317316        xDelete<IssmDouble>(Jlist);
    318317        xDelete<IssmPDouble>(G2);
    319        
    320318        xDelete<char>(rootpath);
    321319        xDelete<char>(inputfilename);

  • issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp

    r23052 r23066  
    126126        SetControlInputsFromVectorx(femmodel,aX);
    127127        xDelete<IssmDouble>(aX);
    128        
     128
    129129        /*Compute solution (forward)*/
    130130        void (*solutioncore)(FemModel*)=NULL;
    … …  
    392392        long niter = long(maxsteps); /*Maximum number of iterations*/
    393393        long nsim  = long(maxiter);/*Maximum number of function calls*/
    394        
     394
    395395        /*Get initial guess*/
    396396        Vector<double> *Xpetsc = NULL;
    … …  
    401401        delete Xpetsc;
    402402        //_assert_(intn==numberofvertices*num_controls);
    403  
     403
    404404        /*Get problem dimension and initialize gradient and initial guess*/
    405405        long n = long(intn);
    … …  
    415415                N_add+=N[c];
    416416        }
    417        
     417
    418418        /*Allocate m1qn3 working arrays (see documentation)*/
    419419        long      m   = 100;
    … …  
    471471                N_add +=N[c];
    472472        }
    473        
     473
    474474        /*Set X as our new control*/
    475475        IssmDouble* aX=xNew<IssmDouble>(intn);
    … …  
    483483        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
    484484        SetControlInputsFromVectorx(femmodel,aX);
    485        
     485
    486486        xDelete(aX);
    487 
    488487
    489488        if (solution_type == TransientSolutionEnum){
    … …  
    506505                        femmodel->results->AddObject(G_output);
    507506                        femmodel->results->AddObject(X_output);
    508                        
     507
    509508                        offset += N[i]*numberofvertices;
    510509                }

  • issm/trunk-jpl/src/c/cores/dakota_core.cpp

    r19634 r23066  
    235235void dakota_core(FemModel* femmodel){  /*{{{*/
    236236
    237 
    238237        int                my_rank;
    239238        char              *dakota_input_file  = NULL;

  • issm/trunk-jpl/src/c/cores/damage_core.cpp

    r19527 r23066  
    1919        char   **requested_outputs = NULL;
    2020
    21                        
    2221        //first recover parameters common to all solutions
    2322        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);

  • issm/trunk-jpl/src/c/cores/esa_core.cpp

    r22349 r23066  
    2323        int        numoutputs        = 0;
    2424        char     **requested_outputs = NULL;
    25        
     25
    2626        /*additional parameters: */
    2727        int  gsize;
    … …  
    8080                U_x = new Vector<IssmDouble>(gsize);
    8181                U_y = new Vector<IssmDouble>(gsize);
    82                
     82
    8383                /*call the geodetic main modlule:*/
    8484                if(domaintype==Domain3DsurfaceEnum){
    … …  
    9595                InputUpdateFromVectorx(femmodel,U_north,EsaNmotionEnum,VertexSIdEnum);  // north motion
    9696                InputUpdateFromVectorx(femmodel,U_east,EsaEmotionEnum,VertexSIdEnum);           // east motion
    97                
     97
    9898                if(save_results){
    9999                        if(VerboseSolution()) _printf0_("   saving results\n");
    … …  
    101101                        femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    102102                }
    103                        
     103
    104104                if(solution_type==EsaSolutionEnum)femmodel->RequestedDependentsx();
    105105
    … …  
    115115}
    116116/*}}}*/
    117 

  • issm/trunk-jpl/src/c/cores/gia_core.cpp

    r22009 r23066  
    5353                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
    5454        }
    55        
     55
    5656        xDelete<IssmDouble>(x);
    5757        xDelete<IssmDouble>(y);

  • issm/trunk-jpl/src/c/cores/love_core.cpp

    r22382 r23066  
    2828        /*parameters: */
    2929        bool save_results;
    30        
     30
    3131        if(VerboseSolution()) _printf0_("   computing LOVE numbers\n");
    3232
    3333        /*Recover some parameters: */
    3434        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    35        
     35
    3636        /*recover love number parameters: */
    3737        femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
    … …  
    4545        femmodel->parameters->FindParam(&love_kernels,LoveKernelsEnum);
    4646        femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);
    47        
     47
    4848        /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
    4949        Matlitho* matlitho=NULL;

  • issm/trunk-jpl/src/c/cores/masstransport_core.cpp

    r21077 r23066  
    3131        femmodel->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
    3232        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
    33                        
     33
    3434        if(VerboseSolution()) _printf0_("   computing mass transport\n");
    3535

  • issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp

    r22982 r23066  
    1313void sealevelrise_core(FemModel* femmodel){ /*{{{*/
    1414
    15 
    1615        /*Parameters, variables:*/
    1716        bool save_results;
    1817        bool isslr=0;
    1918        int solution_type;
    20                
     19
    2120        /*Retrieve parameters:*/
    2221        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
    2322        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    2423        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    25        
     24
    2625        /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
    2726        if(solution_type==SealevelriseSolutionEnum)isslr=1;
    … …  
    4140        /*Run steric core for sure:*/
    4241        steric_core(femmodel);
    43        
     42
    4443        /*Save results: */
    4544        if(save_results){
    … …  
    5857/*}}}*/
    5958void geodetic_core(FemModel* femmodel){ /*{{{*/
    60 
    6159
    6260        /*variables:*/
    … …  
    8886        /*Should we even be here?:*/
    8987        femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); if(!geodetic)return;
    90        
     88
    9189        /*Verbose: */
    9290        if(VerboseSolution()) _printf0_("         computing geodetic sea level rise\n");
    … …  
    155153                /*recover N_esa  = U_esa + RSLg:*/
    156154                N_esa=U_esa->Duplicate(); U_esa->Copy(N_esa); N_esa->AXPY(RSLg,1);
    157                
     155
    158156                /*transform these values into rates (as we only run this once each frequency turn:*/
    159157                N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency));
    … …  
    162160                U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency));
    163161                RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency));
    164                
     162
    165163                /*get some results into elements:{{{*/
    166164                InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum);
    … …  
    182180        }
    183181
    184 
    185182        if(iscoupler){
    186183                /*transfer sea level back to ice caps:*/
    … …  
    228225        Vector<IssmDouble> *U_gia_rate= NULL;
    229226        Vector<IssmDouble> *N_gia_rate= NULL;
    230                
     227
    231228        /*parameters: */
    232229        bool isslr=0;
    … …  
    243240        /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
    244241        if(solution_type==SealevelriseSolutionEnum)isslr=1;
    245        
     242
    246243        /*Should we be here?:*/
    247244        if(!isslr)return;
    … …  
    291288/*}}}*/
    292289Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel){ /*{{{*/
    293  
     290
    294291        /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
    295292
    … …  
    318315        /*Figure out size of g-set deflection vector and allocate solution vector: */
    319316        gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
    320        
     317
    321318        /*Initialize:*/
    322319        RSLgi = new Vector<IssmDouble>(gsize);
    … …  
    331328         * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/
    332329        RSLgi->Shift(-eustatic-RSLgi_oceanaverage);
    333        
     330
    334331        /*save eustatic value for results: */
    335332        femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelRSLEustaticEnum,-eustatic));
    336 
    337333
    338334        /*clean up and return:*/
    … …  
    346342        /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
    347343          non eustatic core of the SLR solution */
    348 
    349344
    350345        Vector<IssmDouble> *RSLg    = NULL;
    … …  
    380375        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    381376        femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum);
    382        
     377
    383378        /*computational flag: */
    384379        femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum);
    … …  
    389384        /*Figure out size of g-set deflection vector and allocate solution vector: */
    390385        gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
    391        
     386
    392387        /*Initialize:*/
    393388        RSLg = new Vector<IssmDouble>(gsize);
    … …  
    400395        count=1;
    401396        converged=false;
    402        
     397
    403398        /*Start loop: */
    404399        for(;;){
    … …  
    413408                /*call the non eustatic module: */
    414409                femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, latitude, longitude, radius,verboseconvolution,loop);
    415        
     410
    416411                /*assemble solution vector: */
    417412                RSLgo->Assemble();
    … …  
    423418                        femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz,latitude,longitude,radius);
    424419                        RSLgo_rot->Assemble();
    425                        
     420
    426421                        /*save changes in inertia tensor as results: */
    427422                        femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz));
    … …  
    431426                        RSLgo->AXPY(RSLgo_rot,1);
    432427                }
    433                
     428
    434429                /*we need to average RSLgo over the ocean: RHS term  5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
    435430                RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo);
    436        
     431
    437432                /*RSLg is the sum of the eustatic term, and the ocean terms: */
    438433                RSLg_eustatic->Copy(RSLg); RSLg->AXPY(RSLgo,1);
    … …  
    455450                        break;
    456451                }       
    457                
     452
    458453                /*some minor verbosing adjustment:*/
    459454                if(count>1)verboseconvolution=false;
    460                
     455
    461456        }
    462457        if(VerboseConvergence()) _printf0_("\n              total number of iterations: " << count-1 << "\n");
    … …  
    474469        Vector<IssmDouble> *U_north_esa   = NULL;
    475470        Vector<IssmDouble> *U_east_esa    = NULL;
    476                
     471
    477472        /*parameters: */
    478473        int configuration_type;
    … …  
    507502        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
    508503        VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
    509        
     504
    510505        /*call the elastic main modlule:*/
    511506        femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,latitude,longitude,radius,xx,yy,zz,loop,horiz);
    … …  
    532527        Vector<IssmDouble> *U_gia  = NULL;
    533528        Vector<IssmDouble> *N_gia  = NULL;
    534        
     529
    535530        /*parameters:*/
    536531        int                                     frequency;
    … …  
    570565        Vector<IssmDouble>* forcingglobal=NULL;
    571566        int*         nvs=NULL;
    572        
     567
    573568        /*transition vectors:*/
    574569        IssmDouble** transitions=NULL;
    … …  
    577572        int*         transitions_n=NULL;
    578573        int          nv;
    579        
     574
    580575        /*communicators:*/
    581576        ISSM_MPI_Comm tocomm;
    … …  
    591586        femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
    592587        my_rank=IssmComm::GetRank();
    593        
     588
    594589        /*retrieve the inter communicators that will be used to send data from each ice cap to the earth: */
    595590        if(modelid==earthid){
    … …  
    620615                                ISSM_MPI_Recv(forcings[i], nvs[i], ISSM_MPI_DOUBLE, 0,i, fromcomms[i], &status);
    621616                        }
    622                        
     617
    623618                }
    624619                else{
    … …  
    631626        /*On the earth model, consolidate all the forcings into one, and update the elements dataset accordingly: {{{*/
    632627        if(modelid==earthid){
    633                
     628
    634629                /*Out of all the delta thicknesses, build one delta thickness vector made of all the ice cap contributions.
    635630                 *First, build the global delta thickness vector in the earth model: */
    … …  
    652647                                /*build index to plug values: */
    653648                                int*        index=xNew<int>(M); for(int i=0;i<M;i++)index[i]=reCast<int>(transition[i])-1; //matlab indexing!
    654 
    655649
    656650                                /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular
    … …  
    663657                /*Assemble vector:*/
    664658                forcingglobal->Assemble();
    665                
     659
    666660                /*Plug into elements:*/
    667661                InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum);
    … …  
    696690        IssmDouble*  forcing=NULL;
    697691        IssmDouble*  forcingglobal=NULL;
    698        
     692
    699693        /*transition vectors:*/
    700694        IssmDouble** transitions=NULL;
    … …  
    703697        int*         transitions_n=NULL;
    704698        int          nv;
    705        
     699
    706700        /*communicators:*/
    707701        ISSM_MPI_Comm fromcomm;
    … …  
    718712        femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
    719713        my_rank=IssmComm::GetRank();
    720        
     714
    721715        /*retrieve the inter communicators that will be used to send data from earth to ice caps:*/
    722716        if(modelid==earthid){
    … …  
    733727        }
    734728
    735 
    736729        /*Retrieve sea-level on earth model: */
    737730        if(modelid==earthid){
    … …  
    742735        /*Send the forcing to the ice caps:{{{*/
    743736        if(my_rank==0){
    744                
     737
    745738                if(modelid==earthid){
    746                        
     739
    747740                        /*Retrieve transition vectors, used to figure out global forcing contribution to each ice cap's own elements: */
    748741                        femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
    749                        
     742
    750743                        if(ntransitions!=earthid)_error_("TransferSeaLevel error message: number of transition vectors is not equal to the number of icecaps!");
    751744
    … …  
    804797        Vector<IssmDouble> *cumdeltathickness    = NULL;
    805798        int nv;
    806        
     799
    807800        if(VerboseSolution()) _printf0_("              computing earth mass transport\n");
    808801
    … …  
    840833} /*}}}*/
    841834void slrconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
    842        
     835
    843836        bool converged=true;
    844837        IssmDouble ndS,nS;
    … …  
    848841        dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0);
    849842        ndS=dRSLg->Norm(NORM_TWO);
    850        
     843
    851844        if (xIsNan<IssmDouble>(ndS)) _error_("convergence criterion is NaN!");
    852        
     845
    853846        if(!xIsNan<IssmDouble>(eps_rel)){
    854847                nS=RSLg_old->Norm(NORM_TWO);
    855848                if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN!");
    856849        }
    857 
    858850
    859851        //clean up

  • issm/trunk-jpl/src/c/cores/smb_core.cpp

    r19528 r23066  
    2929        femmodel->parameters->FindParam(&numoutputs,SmbNumRequestedOutputsEnum);
    3030        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SmbRequestedOutputsEnum);
    31                        
     31
    3232        if(VerboseSolution()) _printf0_("   computing smb \n");
    33  
     33
    3434        analysis = new SmbAnalysis();
    3535        analysis->Core(femmodel);

  • issm/trunk-jpl/src/c/cores/stressbalance_core.cpp

    r23030 r23066  
    2222        char     **requested_outputs = NULL;
    2323        Analysis  *analysis          = NULL;
    24                        
    2524
    2625        /* recover parameters:*/
    … …  
    3635        femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum);
    3736        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,StressbalanceRequestedOutputsEnum);
    38        
     37
    3938        if(VerboseSolution()) _printf0_("   computing new velocity\n");
    4039
    … …  
    7978        }
    8079
    81 
    8280        if(save_results){
    8381                if(VerboseSolution()) _printf0_("   saving results\n");
    8482                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    8583        }
    86        
     84
    8785        if(solution_type==StressbalanceSolutionEnum && !control_analysis)femmodel->RequestedDependentsx();
    8886

  • issm/trunk-jpl/src/c/datastructures/DataSet.cpp

    r21701 r23066  
    9292/*Specific methods*/
    9393void  DataSet::Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
    94        
     94
    9595        vector<Object*>::iterator obj;
    9696        int obj_size=0;

  • issm/trunk-jpl/src/c/main/esmfbinders.cpp

    r20972 r23066  
    22 * \brief: ESMF binders for ISSM. Binders developed initially for the GEOS-5 framework.
    33 */
    4 
    54
    65#include "./issm.h"
    … …  
    3635                /*Some specific code here for the binding: */
    3736                femmodel->parameters->SetParam(SMBgcmEnum,SmbEnum); //bypass SMB model, will be provided by GCM!
    38        
     37
    3938                /*Restart file: */
    4039                femmodel->Restart();
    … …  
    113112
    114113                                                IssmDouble surface;
    115                                                
     114
    116115                                                /*Recover surface from the ISSM element: */
    117116                                                Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
    118117                                                surface_input->GetInputAverage(&surface);
    119                        
     118
    120119                                                *(issmoutputs+f*numberofelements+i) = surface;
    121120
    … …  
    147146                /*Output results: */
    148147                OutputResultsx(femmodel);
    149                        
     148
    150149                /*Check point: */
    151150                femmodel->CheckPoint();

  • issm/trunk-jpl/src/c/main/issm_dakota.cpp

    r21876 r23066  
    1616int main(int argc,char **argv){
    1717
    18 
    1918        #if defined(_HAVE_DAKOTA_) && _DAKOTA_MAJOR_ >= 6
    2019
    … …  
    2928        /*Initialize MPI: */
    3029        ISSM_MPI_Init(&argc, &argv); // initialize MPI
    31        
     30
    3231        /*Recover file name for dakota input file:*/
    3332        dakota_input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".qmu.in")+2));
    3433        sprintf(dakota_input_file,"%s/%s%s",argv[2],argv[3],".qmu.in");
    35        
     34
    3635        dakota_output_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".qmu.out")+2));
    3736        sprintf(dakota_output_file,"%s/%s%s",argv[2],argv[3],".qmu.out");
    … …  
    5756                Dakota::abort_handler(-1);
    5857        }
    59        
     58
    6059        Dakota::ProblemDescDB& problem_db = env.problem_description_db();
    6160        Dakota::ModelLIter ml_iter;

  • issm/trunk-jpl/src/c/main/issm_ocean.cpp

    r22656 r23066  
    2222        /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/
    2323        worldcomm=EnvironmentInit(argc,argv);
    24        
     24
    2525        /*What is my rank?:*/
    2626        ISSM_MPI_Comm_rank(worldcomm,&my_rank);
    … …  
    3939
    4040        FemModel *femmodel = new FemModel(argc,argv,modelcomm);
    41        
     41
    4242        /*Now that the models are initialized, keep communicator information in the parameters datasets of each model: */
    4343        femmodel->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(worldcomm,WorldCommEnum));

  • issm/trunk-jpl/src/c/main/issm_slr.cpp

    r22190 r23066  
    2828        /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/
    2929        worldcomm=EnvironmentInit(argc,argv);
    30        
     30
    3131        /*What is my rank?:*/
    3232        ISSM_MPI_Comm_rank(worldcomm,&my_rank);
    … …  
    4040        for(int i=0;i<nummodels;i++){
    4141                char* string=NULL;
    42                
     42
    4343                string=xNew<char>(strlen(argv[5+3*i])+1);
    4444                xMemCpy<char>(string,argv[5+3*i],strlen(argv[5+3*i])+1);
    4545                dirnames[i]=string;
    46                
     46
    4747                string=xNew<char>(strlen(argv[5+3*i+1])+1);
    4848                xMemCpy<char>(string,argv[5+3*i+1],strlen(argv[5+3*i+1])+1);
    … …  
    9494        /*Initialize femmodel from arguments provided command line: */
    9595        FemModel *femmodel = new FemModel(4,arguments,modelcomm);
    96        
     96
    9797        /*Now that the models are initialized, keep communicator information in the parameters datasets of each model: */
    9898        femmodel->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(worldcomm,WorldCommEnum));

  • issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp

    r21889 r23066  
    7878                if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
    7979        }
    80        
     80
    8181        /*Free ressources: */
    8282        xDelete<char>(toolkittype);

  • issm/trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp

    r22379 r23066  
    3939                IssmDouble d0=1.e+10;
    4040                IssmDouble xi,yi;
    41                
     41
    4242                //recover vertex position:
    4343                xi=x[i];  yi=y[i];

  • issm/trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp

    r22458 r23066  
    5252        double x2,y2;
    5353        double mind;
    54        
     54
    5555        for(int i=i0;i<i1;i++){
    5656
    … …  
    9393      return sqrt(pow(x2-x0,2)+pow(y2-y0,2));
    9494   }
    95        
     95
    9696   /*Projection falls on segment*/
    9797        double projx= x1 + t* (x2-x1);

  • issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp

    r23053 r23066  
    107107        femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum);
    108108        femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
    109        
     109
    110110        IssmDouble* boxareas=xNewZeroInit<IssmDouble>(num_basins*maxbox);
    111111
    … …  
    122122        ISSM_MPI_Allreduce(boxareas,sumareas,num_basins*maxbox,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    123123        //if(sumareas[0]==0){_error_("No elements in box 0, basal meltrates will be 0. Consider decreasing md.basalforcings.maxboxcount or refining your mesh!");}
    124        
     124
    125125        /*Update parameters to keep track of the new areas in future calculations*/
    126126        femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoBoxAreaEnum,sumareas,maxbox*num_basins));

  • issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp

    r22953 r23066  
    6767}
    6868/*}}}*/
    69 

  • issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp

    r22828 r23066  
    2727                int size = 0;
    2828                for(int i=0;i<num_controls;i++) size+=M[i]*N[i];
    29 
    3029
    3130                /*2. Allocate vector*/
    … …  
    5958                parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    6059                parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    61 
    6260
    6361                /*2. Allocate vector*/

  • issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp

    r22952 r23066  
    8484        /*this one is special: we don't specify the type, but let the nature of the inputs dictace.
    8585         * P0 -> ElementSIdEnum, P1 ->VertexSIdEnum: */
    86        
     86
    8787        /*We go find the input of the first element, and query its interpolation type: */
    8888        Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(0));

  • issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp

    r23035 r23066  
    1515        Element            *element                          = NULL;
    1616
    17 
    1817        /*retrieve parameters: */
    1918        parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
    … …  
    2120
    2221        if(migration_style==NoneEnum) return;
    23        
     22
    2423        if(VerboseModule()) _printf0_("   Migrating grounding line\n");
    2524

  • issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp

    r21915 r23066  
    166166
    167167        /*edge can be either IntersectEnum (one and only one intersection between the edge and the segment), ColinearEnum (edge and segment are collinear) and SeparateEnum (no intersection): */
    168                
     168
    169169        /*if (el==318 && contouri==9){
    170170                _printf_(edge1 << " " << edge2 << " " << edge3 << " "  << alpha1 << " " << alpha2 << " " << beta1 << " " << beta2 << " " << gamma1 << " " << gamma2 << " " << xsegment[0] << " "  << xsegment[1] << " " << ysegment[0] << " " << ysegment[1] << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << ynodes[2]);
    171        
     171
    172172        _printf_("Bool" << (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum));
    173173        }*/
    … …  
    199199                segments_dataset->AddObject(new  Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
    200200
    201                
    202201                /*This case is impossible: not quite! */
    203202                //_printf_(alpha1 << " " << alpha2 << " " << beta1 << " " << beta2 << " " << gamma1 << " " << gamma2 << " " << xsegment[0] << " "  << xsegment[1] << " " << ysegment[0] << " " << ysegment[1] << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << ynodes[2]);
    … …  
    226225        }
    227226        else if(  (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum)   ){
    228        
     227
    229228                /*if (el==318 && contouri==9){
    230229                        _printf_("hello" <<  " NodeInElement 0 " << (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) <<  " NodeInElement 1 " << (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])));
    … …  
    251250                                IsIdenticalNode(xnodes[1],ynodes[1],xsegment[0],ysegment[0],tolerance) ||
    252251                                IsIdenticalNode(xnodes[2],ynodes[2],xsegment[0],ysegment[0],tolerance)){
    253                                
     252
    254253                                /*ok, segments[0] is common to one of our vertices: */
    255254                                coord1=0;

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp

    r22612 r23066  
    6969                        int            num_independent_objects,M;
    7070                        char**         names                   = NULL;
    71                                
     71
    7272                                /*this is done somewhere else*/
    7373                                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_independent_objects",InversionNumControlParametersEnum));
    7474                           parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_dependent_objects",InversionNumCostFunctionsEnum));
    75                                
     75
    7676                                /*Step1: create controls (independents)*/
    7777                                iomodel->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects");
    … …  
    9696                                xDelete<char*>(cm_responses);
    9797                                parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,costfunc_enums,num_costfunc));
    98                                
     98
    9999                                iomodel->FetchData(&control_scaling_factors,NULL,NULL,"md.autodiff.independent_scaling_factors");
    100100                                parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_independent_objects));
    101        
     101
    102102                                /*cleanup*/
    103103                                for(int i=0;i<num_independent_objects;i++){

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

    r22753 r23066  
    88#include "../../MeshPartitionx/MeshPartitionx.h"
    99#include "../ModelProcessorx.h"
    10 
    1110
    1211#if !defined(_HAVE_ADOLC_)
    … …  
    2221        char     **cost_functions   = NULL;
    2322
    24 
    2523        /*Fetch parameters: */
    2624        iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
    … …  
    2927        /*Now, return if no control*/
    3028        if(!control_analysis) return;
    31        
     29
    3230        /*Process controls and convert from string to enums*/
    3331        iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters");
    … …  
    4846
    4947        iomodel->FetchData(3,"md.inversion.cost_functions_coefficients","md.inversion.min_parameters","md.inversion.max_parameters");
    50        
     48
    5149        /*Fetch Observations */
    5250        iomodel->FindConstant(&domaintype,"md.mesh.domain_type");
    … …  
    154152        iomodel->FetchData(&types,NULL,NULL,"md.autodiff.independent_object_types");
    155153
    156                
    157154        M_all = xNew<int>(num_independent_objects);
    158155
    … …  
    162159        iomodel->FetchData(&independents_fullmax,&M_par,&N_par,"md.autodiff.independent_max_parameters");
    163160        iomodel->FetchData(&control_sizes,NULL,NULL,"md.autodiff.independent_control_sizes");
    164        
     161
    165162        int* start_point = NULL;
    166163        start_point = xNew<int>(num_independent_objects);
    … …  
    180177                        IssmDouble*     independents_min                        = NULL;
    181178                        IssmDouble*        independents_max                     = NULL;
    182        
     179
    183180                        FieldAndEnumFromCode(&input_enum,&iofieldname,names[i]);
    184181
    … …  
    187184                        _assert_(independent);
    188185                        _assert_(N==control_sizes[i]);
    189                
     186
    190187                        independents_min = NULL; independents_min = xNew<IssmDouble>(M*N);
    191188                        independents_max = NULL; independents_max = xNew<IssmDouble>(M*N);
    … …  
    206203                        xDelete<IssmDouble>(independents_min);
    207204                        xDelete<IssmDouble>(independents_max);
    208        
     205
    209206                }
    210207                else{

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

    r23053 r23066  
    3535        /*Create elements*/
    3636        if(control_analysis && !adolc_analysis)iomodel->FetchData(2,"md.inversion.min_parameters","md.inversion.max_parameters");
    37        
     37
    3838        switch(iomodel->meshelementtype){
    3939                case TriaEnum:
    … …  
    123123                        break;
    124124                case MaterialsEnum:
    125        
     125
    126126                        //we have several types of materials. Retrieve this info first:
    127127                        iomodel->FetchData(&nature,&nnat,&dummy,"md.materials.nature");
    … …  
    235235        else iomodel->FetchDataToInput(elements,"md.mesh.scale_factor",MeshScaleFactorEnum,1.);
    236236        if (isoceancoupling) iomodel->FetchData(2,"md.mesh.lat","md.mesh.long");
    237        
     237
    238238        CreateNumberNodeToElementConnectivity(iomodel,solution_type);
    239239

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp

    r22971 r23066  
    1111
    1212        int i,j;
    13        
     13
    1414        DataSet*     output_definitions      = NULL;
    1515        int*         output_definition_enums = NULL;
    … …  
    6565                        else if (output_definition_enums[i]==MisfitEnum){
    6666                                /*Deal with misfits: {{{*/
    67                        
     67
    6868                                /*misfit variables: */
    6969                                int          nummisfits;
    … …  
    115115                                         _error_("misfit weight size not supported yet");
    116116
    117 
    118117                                        /*First create a misfit object for that specific string (misfit_model_string_s[j]):*/
    119118                                        output_definitions->AddObject(new Misfit(misfit_name_s[j],StringToEnumx(misfit_definitionstring_s[j]),StringToEnumx(misfit_model_string_s[j]),StringToEnumx(misfit_observation_string_s[j]),misfit_timeinterpolation_s[j],misfit_local_s[j],StringToEnumx(misfit_weights_string_s[j])));
    … …  
    159158                        else if (output_definition_enums[i]==CfsurfacesquareEnum){
    160159                                /*Deal with cfsurfacesquare: {{{*/
    161                                
     160
    162161                                /*cfsurfacesquare variables: */
    163162                                int          num_cfsurfacesquares;
    … …  
    214213
    215214                                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
    216                                                
     215
    217216                                                element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_observation_s[j], iomodel,cfsurfacesquare_observation_M_s[j],cfsurfacesquare_observation_N_s[j],obs_vector_type,StringToEnumx(cfsurfacesquare_observation_string_s[j]),7,SurfaceObservationEnum);
    218217                                                element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_weights_s[j], iomodel,cfsurfacesquare_weights_M_s[j],cfsurfacesquare_weights_N_s[j],weight_vector_type,StringToEnumx(cfsurfacesquare_weights_string_s[j]),7,WeightsSurfaceObservationEnum);
    … …  
    251250                        else if (output_definition_enums[i]==CfdragcoeffabsgradEnum){
    252251                                /*Deal with cfdragcoeffabsgrad: {{{*/
    253                                
     252
    254253                                /*cfdragcoeffabsgrad variables: */
    255254                                int          num_cfdragcoeffabsgrads;
    … …  
    286285
    287286                                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
    288                                                
     287
    289288                                                element->DatasetInputAdd(StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),cfdragcoeffabsgrad_weights_s[j], iomodel,cfdragcoeffabsgrad_weights_M_s[j],cfdragcoeffabsgrad_weights_N_s[j],weight_vector_type,StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),7,WeightsSurfaceObservationEnum);
    290289
    … …  
    313312                        else if (output_definition_enums[i]==CfsurfacelogvelEnum){
    314313                                /*Deal with cfsurfacelogvel: {{{*/
    315                                
     314
    316315                                /*cfsurfacelogvel variables: */
    317316                                int          num_cfsurfacelogvels;
    … …  
    369368
    370369                                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
    371                                                
     370
    372371                                                element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vxobs[j], iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vxobs_string[j]),7,VxObsEnum);
    373372                                                        element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vyobs[j], iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vyobs_string[j]),7,VyObsEnum);
    … …  
    409408                        else if (output_definition_enums[i]==NodalvalueEnum){
    410409                                /*Deal with nodal values: {{{*/
    411                                
     410
    412411                                /*nodal value variables: */
    413412                                int          numnodalvalues;
    … …  
    428427                                        output_definitions->AddObject(new Nodalvalue(nodalvalue_name_s[j],StringToEnumx(nodalvalue_definitionstrings[j]),StringToEnumx(nodalvalue_modelstrings[j]),nodalvalue_node_s[j]-1)); //-1 because matlab to c indexing.
    429428                                }
    430                                        
     429
    431430                                /*Free ressources:*/
    432431                                for(j=0;j<numnodalvalues;j++){
    … …  
    482481                        else if (output_definition_enums[i]==MassconaxpbyEnum){
    483482                                /*Deal with masscon combinations: {{{*/
    484                                
     483
    485484                                /*masscon variables: */
    486485                                char**       masscon_name_s             = NULL;   
    … …  
    617616                                        output_definitions->AddObject(new Numberedcostfunction(ncf_name_s[j],StringToEnumx(ncf_definitionstring_s[j]),num_cost_functions,cost_function_enums));
    618617                                }
    619                                
     618
    620619                                /*Free data: */
    621620                                iomodel->DeleteData(2,"md.numberedcostfunction.name","md.numberedcostfunction.definitionstring");

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r22993 r23066  
    2323        char*       fieldname = NULL;
    2424        IssmDouble  time;
    25        
     25
    2626        /*parameters for mass flux:*/
    2727        int          mass_flux_num_profiles     = 0;
    … …  
    6262        parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum));
    6363        parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1)); 
    64        
    6564
    6665          {/*This is specific to ice...*/

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp

    r22739 r23066  
    5555                delete analysis;
    5656
    57 
    5857                /* Update counters, because we have created more nodes, loads and
    5958                 * constraints, and ids for objects created in next call to CreateDataSets

  • issm/trunk-jpl/src/c/modules/NodalValuex/NodalValuex.cpp

    r20038 r23066  
    2121        cpu_found=-1;
    2222        found=0;
    23        
     23
    2424        for(int i=0;i<elements->Size();i++){
    2525                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));

  • issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp

    r22386 r23066  
    2222                        /*This is the object that we have been chasing for. compute the response and return: */
    2323                        IssmDouble return_value=definition->Response(femmodel);
    24                
     24
    2525                        /*cleanup: */
    2626                        xDelete<char>(name);
    … …  
    3131                xDelete<char>(name);
    3232        }
    33        
     33
    3434        /*If we are here, did not find the definition for this response, not good!: */
    3535        _error_("Could not find the response for output definition " << output_string << " because could not find the definition itself!");
    … …  
    4444        /*Now, go through the output definitions, and retrieve the object which corresponds to our requested response, output_enum: */
    4545        for(int i=0;i<output_definitions->Size();i++){
    46                
     46
    4747                //Definition* definition=xDynamicCast<Definition*>(output_definitions->GetObjectByOffset(i));
    4848                Definition* definition=dynamic_cast<Definition*>(output_definitions->GetObjectByOffset(i));
    … …  
    5353                        /*This is the object that we have been chasing for. compute the response and return: */
    5454                        IssmDouble return_value=definition->Response(femmodel);
    55                
     55
    5656                        /*return:*/
    5757                        return return_value;
    5858                }
    5959        }
    60        
     60
    6161        /*If we are here, did not find the definition for this response, not good!: */
    6262        _error_("Could not find the response for output definition " << EnumToStringx(output_enum)

  • issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp

    r22004 r23066  
    9292        /* Intermediaries */
    9393        int numvertices = element->GetNumberOfVertices();
    94        
     94
    9595        if(element->IsIceInElement()){
    9696                for(int i = 0;i<numvertices;i++){

  • issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp

    r22739 r23066  
    3232                }
    3333
    34 
    3534                xDelete<int>(control_type);
    3635        }

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

    r22840 r23066  
    107107        xDelete(dzT);
    108108        xDelete(dzB);
    109        
     109
    110110        //---------NEED TO IMPLEMENT A PROPER GRID STRECHING ALGORITHM------------
    111111
    … …  
    202202        gdn=*pgdn;
    203203        gsp=*pgsp;
    204        
     204
    205205        /*only when aIdx = 1 or 2 do we run grainGrowth: */
    206206        if(aIdx!=1 & aIdx!=2){
    … …  
    220220        //set maximum water content by mass to 9 percent (Brun, 1980)
    221221        for(int i=0;i<m;i++)if(lwc[i]>9.0-Wtol) lwc[i]=9.0;
    222 
    223222
    224223        /* Calculate temperature gradiant across grid cells
    … …  
    244243        // take absolute value of temperature gradient
    245244        for(int i=0;i<m;i++)dT[i]=fabs(dT[i]);
    246        
     245
    247246        /*Snow metamorphism. Depends on value of snow dendricity and wetness of the snowpack: */
    248247        for(int i=0;i<m;i++){
    … …  
    324323                re[i]=gsz[i]/2.0;
    325324        }
    326        
     325
    327326        /*Free ressources:*/
    328327        xDelete<IssmDouble>(gsz);
    … …  
    348347        //// Inputs
    349348        // aIdx      = albedo method to use
    350        
     349
    351350        // Method 0
    352351        //  aValue   = direct input value for albedo, override all changes to albedo
    … …  
    397396        /*Recover pointers: */
    398397        a=*pa;
    399        
     398
    400399        //some constants:
    401400        const IssmDouble dSnow = 300;   // density of fresh snow [kg m-3]       
    … …  
    513512
    514513        /* ENGLACIAL THERMODYNAMICS*/
    515          
     514
    516515        /* Description:
    517516           computes new temperature profile accounting for energy absorption and
    … …  
    537536        // T: grid cell temperature [k]
    538537        // EC: evaporation/condensation [kg]
    539        
     538
    540539        /*intermediary: */
    541540        IssmDouble* K = NULL;
    … …  
    580579        IssmDouble dlw=0.0;
    581580        IssmDouble dT_dlw=0.0;
    582        
     581
    583582        /*outputs:*/
    584583        IssmDouble EC=0.0;
    … …  
    597596        //initialize Evaporation - Condenstation
    598597        EC = 0.0;
    599        
     598
    600599        // check if all SW applied to surface or distributed throught subsurface
    601600        // swIdx = length(swf) > 1
    … …  
    609608        // if V = 0 goes to infinity therfore if V = 0 change
    610609        if(V<0.01-Dtol)V=0.01;
    611        
     610
    612611        // Bulk-transfer coefficient for turbulent fluxes
    613612        An =  pow(0.4,2) / pow(log(Tz/z0),2);     // Bulk-transfer coefficient
    … …  
    627626
    628627        // THERMAL DIFFUSION COEFFICIENTS
    629  
     628
    630629        // A discretization scheme which truncates the Taylor-Series expansion
    631630        // after the 3rd term is used. See Patankar 1980, Ch. 3&4
    632  
     631
    633632        // discretized heat equation:
    634  
     633
    635634        //                 Tp = (Au*Tu° + Ad*Td° + (Ap-Au-Ad)Tp° + S) / Ap
    636  
     635
    637636        // where neighbor coefficients Au, Ap, & Ad are
    638  
     637
    639638        //                   Au = [dz_u/2KP + dz_p/2KE]^-1
    640639        //                   Ad = [dz_d/2KP + dz_d/2KD]^-1
    … …  
    649648        KD = xNew<IssmDouble>(m);
    650649        KP = xNew<IssmDouble>(m);
    651        
     650
    652651        KU[0] = UNDEF;
    653652        KD[m-1] = UNDEF;
    … …  
    661660        dzU[0]=UNDEF;
    662661        dzD[m-1]=UNDEF;
    663        
     662
    664663        for(int i=1;i<m;i++) dzU[i]= dz[i-1];
    665664        for(int i=0;i<m-1;i++) dzD[i] = dz[i+1];
    … …  
    693692                Ap[i] = (d[i]*dz[i]*CI)/dt;
    694693        }
    695        
     694
    696695        // create "neighbor" coefficient matrix
    697696        Nu = xNew<IssmDouble>(m);
    … …  
    703702                Np[i]= 1.0 - Nu[i] - Nd[i];
    704703        }
    705        
     704
    706705        // specify boundary conditions: constant flux at bottom
    707706        Nu[m-1] = 0.0;
    708707        Np[m-1] = 1.0;
    709        
     708
    710709        // zero flux at surface
    711710        Np[0] = 1.0 - Nd[0];
    712        
     711
    713712        // Create neighbor arrays for diffusion calculations instead of a tridiagonal matrix
    714713        Nu[0] = 0.0;
    715714        Nd[m-1] = 0.0;
    716        
     715
    717716        /* RADIATIVE FLUXES*/
    718717
    … …  
    720719        sw = xNew<IssmDouble>(m);
    721720        for(int i=0;i<m;i++) sw[i]= swf[i] * dt;
    722        
     721
    723722        // temperature change due to SW
    724723        dT_sw = xNew<IssmDouble>(m);
    … …  
    746745                // store initial temperature
    747746                //T_init = T;
    748    
     747
    749748                // calculate temperature of snow surface (Ts)
    750749                // when incoming SW radition is allowed to penetrate the surface,
    … …  
    754753                Ts = (T[0] + T[1])/2.0;
    755754                Ts = fmin(CtoK,Ts);    // don't allow Ts to exceed 273.15 K (0 degC)
    756                
     755
    757756                //TURBULENT HEAT FLUX
    758    
     757
    759758                // Monin-Obukhov Stability Correction
    760759                // Reference:
    … …  
    764763                // calculate the Bulk Richardson Number (Ri)
    765764                Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
    766                
     765
    767766                // calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
    768    
     767
    769768                // do not allow Ri to exceed 0.19
    770769                Ri = fmin(Ri, 0.19);
    … …  
    778777                        coefM =pow (1.0-18*Ri,-0.25);
    779778                }
    780                
     779
    781780                // calculate heat/wind 'coef_H' stability factor
    782781                if (Ri <= -0.03+Ttol) coefH = coefM/1.3;
    783782                else coefH = coefM;
    784                
     783
    785784                //// Sensible Heat
    786785                // calculate the sensible heat flux [W m-2](Patterson, 1998)
    … …  
    801800                else{
    802801                        L = LS; // latent heat of sublimation
    803                        
     802
    804803                        // for an ice surface Murphy and Koop, 2005 [Equation 7]
    805804                        eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts);
    … …  
    823822                ulw = - (SB * pow(Ts,4.0)* teValue) * dt ;
    824823                dT_ulw = ulw / TCs;
    825                
     824
    826825                // new grid point temperature
    827    
     826
    828827                //SW penetrates surface
    829828                for(int j=0;j<m;j++) T[j] = T[j] + dT_sw[j];
    830829                T[0] = T[0] + dT_dlw + dT_ulw + dT_turb;
    831                
     830
    832831                // temperature diffusion
    833832                for(int j=0;j<m;j++) T0[1+j]=T[j];
    … …  
    835834                for(int j=0;j<m;j++) Td[j] = T0[2+j];
    836835                for(int j=0;j<m;j++) T[j] = (Np[j] * T[j]) + (Nu[j] * Tu[j]) + (Nd[j] * Td[j]);
    837                
     836
    838837                // calculate cumulative evaporation (+)/condensation(-)
    839838                EC = EC + (EC_day/dts)*dt;
    … …  
    873872        xDelete<IssmDouble>(Td);
    874873
    875 
    876874        /*Assign output pointers:*/
    877875        *pEC=EC;
    … …  
    901899        //   swf     = absorbed shortwave radiation [W m-2]
    902900        //
    903        
     901
    904902        /*outputs: */
    905903        IssmDouble* swf=NULL;
    … …  
    910908        swf=xNewZeroInit<IssmDouble>(m);
    911909
    912 
    913910        // SHORTWAVE FUNCTION
    914911        if (swIdx == 0) {// all sw radation is absorbed in by the top grid cell
    915        
     912
    916913                // calculate surface shortwave radiation fluxes [W m-2]
    917914                swf[0] = (1.0 - as) * dsw;
    … …  
    948945                        }
    949946
    950 
    951947                        // spectral albedos:
    952948                        // 0.3 - 0.8um
    … …  
    988984                                B2_cum[i+1]=cum2;
    989985                        }
    990 
    991986
    992987                        // flux across grid cell boundaries
    … …  
    10151010                        xDelete<IssmDouble>(Qs1);
    10161011                        xDelete<IssmDouble>(Qs2);
    1017                        
    1018                        
     1012
    10191013                }
    10201014                else{  //function of grid cell density
    … …  
    11441138
    11451139        if (P > 0.0+Dtol){
    1146                        
    11471140
    11481141                if (T_air <= CtoK+Ttol){ // if snow
    … …  
    12101203
    12111204                mass_diff = mass - massinit - P;
    1212                
     1205
    12131206                #ifndef _HAVE_ADOLC_  //avoid round operation. only check in forward mode.
    12141207                mass_diff = round(mass_diff * 100.0)/100.0;
    … …  
    12531246        IssmDouble* M=NULL;
    12541247        int*       D=NULL;
    1255        
     1248
    12561249        IssmDouble sumM=0.0;
    12571250        IssmDouble sumER=0.0;
    … …  
    12651258        IssmDouble X=0.0;
    12661259        IssmDouble Wi=0.0;
    1267    
     1260
    12681261        IssmDouble Ztot=0.0;
    12691262        IssmDouble T_bot=0.0;
    … …  
    12791272        IssmDouble EW_bot=0.0;
    12801273        bool        top=false;
    1281    
     1274
    12821275        IssmDouble* Zcum=NULL;
    12831276        IssmDouble* dzMin2=NULL;
    … …  
    12861279        int X1=0;
    12871280        int X2=0;
    1288    
     1281
    12891282        int        D_size = 0;
    12901283
    … …  
    13031296        int         n=*pn;
    13041297        IssmDouble* R=NULL;
    1305        
     1298
    13061299        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   melt module\n");
    13071300
    … …  
    13351328        // for(int i=0;i<n;i++) T[i]-=exsT[i];
    13361329        for(int i=0;i<n;i++) T[i]=fmin(T[i],CtoK);
    1337        
     1330
    13381331        // specify irreducible water content saturation [fraction]
    13391332        const IssmDouble Swi = 0.07;                     // assumed constant after Colbeck, 1974
    … …  
    14931486                                T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature
    14941487
    1495 
    14961488                                // check if an ice layer forms
    14971489                                if (fabs(d[i] - dIce) < Dtol){
    … …  
    15051497                }
    15061498
    1507 
    15081499                //// GRID CELL SPACING AND MODEL DEPTH
    15091500                for(int i=0;i<n;i++)if (W[i] < 0.0 - Wtol) _error_("negative pore water generated in melt equations");
    1510                
     1501
    15111502                // delete all cells with zero mass
    15121503                // adjust pore water
    … …  
    15331524                n=D_size;
    15341525                xDelete<int>(D);
    1535        
     1526
    15361527                // calculate new grid lengths
    15371528                for(int i=0;i<n;i++)dz[i] = m[i] / d[i];
    … …  
    15451536        Zcum=xNew<IssmDouble>(n);
    15461537        dzMin2=xNew<IssmDouble>(n);
    1547    
     1538
    15481539        Zcum[0]=dz[0]; // Compute a cumulative depth vector
    1549    
     1540
    15501541        for (int i=1;i<n;i++){
    15511542                Zcum[i]=Zcum[i-1]+dz[i];
    15521543        }
    1553    
     1544
    15541545        for (int i=0;i<n;i++){
    15551546                if (Zcum[i]<=zTop+Dtol){
    … …  
    15691560                }
    15701561        }
    1571    
     1562
    15721563        //Last cell has to be treated separately if has to be merged (no underlying cell so need to merge with overlying cell)
    15731564        if(X==n-1){
    … …  
    15891580                        gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new;
    15901581                        gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new;
    1591            
     1582
    15921583                        // merge with underlying grid cell and delete old cell
    15931584                        dz[i+1] = dz[i] + dz[i+1];                 // combine cell depths
    … …  
    15951586                        W[i+1] = W[i+1] + W[i];                     // combine liquid water
    15961587                        m[i+1] = m_new;                             // combine top masses
    1597            
     1588
    15981589                        // set cell to 99999 for deletion
    15991590                        m[i] = -99999;
    … …  
    16111602                        }
    16121603                }
    1613        
     1604
    16141605                // adjust variables as a linearly weighted function of mass
    16151606                IssmDouble m_new = m[X2] + m[X1];
    … …  
    16191610                gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
    16201611                gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new;
    1621        
     1612
    16221613                // merge with underlying grid cell and delete old cell
    16231614                dz[X1] = dz[X2] + dz[X1];                 // combine cell depths
    … …  
    16251616                W[X1] = W[X1] + W[X2];                     // combine liquid water
    16261617                m[X1] = m_new;                             // combine top masses
    1627        
     1618
    16281619                // set cell to 99999 for deletion
    16291620                m[X2] = -99999;
    … …  
    16481639        n=D_size;
    16491640        xDelete<int>(D);
    1650    
     1641
    16511642        // check if any of the top 10 cell depths are too large
    16521643        X=0;
    … …  
    16571648                }
    16581649        }
    1659        
     1650
    16601651        int j=0;
    16611652        while(j<=X){
    … …  
    16861677        // calculate total model depth
    16871678        Ztot = cellsum(dz,n);
    1688    
     1679
    16891680        if (Ztot < zMin-Dtol){
    16901681                // printf("Total depth < zMin %f \n", Ztot);
    … …  
    16921683                mAdd = m[n-1]+W[n-1];
    16931684                addE = T[n-1]*m[n-1]*CI + W[n-1]*(LF+CtoK*CI);
    1694        
     1685
    16951686                // add a grid cell of the same size and temperature to the bottom
    16961687                dz_bot=dz[n-1];
    … …  
    17051696                EI_bot=EI[n-1];
    17061697                EW_bot=EW[n-1];
    1707        
     1698
    17081699                dz_add=dz_bot;
    1709        
     1700
    17101701                newcell(&dz,dz_bot,top,n);
    17111702                newcell(&T,T_bot,top,n);
    … …  
    17271718                addE = -(T[n-1]*m[n-1]*CI) - (W[n-1]*(LF+CtoK*CI));
    17281719                dz_add=-(dz[n-1]);
    1729        
     1720
    17301721                // remove a grid cell from the bottom
    17311722                D_size=n-1;
    17321723                D=xNew<int>(D_size);
    1733        
     1724
    17341725                for(int i=0;i<n-1;i++) D[i]=i;
    17351726                celldelete(&dz,n,D,D_size);
    … …  
    17681759                        << "dm: " << dm << " dE: " << dE << "\n");
    17691760        #endif
    1770        
     1761
    17711762        /*Free ressources:*/
    17721763        if(m)xDelete<IssmDouble>(m);
    … …  
    17841775        xDelete<IssmDouble>(Zcum);
    17851776        xDelete<IssmDouble>(dzMin2);
    1786    
     1777
    17871778        /*Assign output pointers:*/
    17881779        *pM=sumM;
    … …  
    18621853        IssmDouble* dz=NULL;
    18631854        IssmDouble* d=NULL;
    1864        
     1855
    18651856        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   densification module\n");
    18661857
    … …  
    18711862        // initial mass
    18721863        IssmDouble* mass_init = xNew<IssmDouble>(m);for(int i=0;i<m;i++) mass_init[i]=d[i] * dz[i];
    1873        
     1864
    18741865        /*allocations and initialization of overburden pressure and factor H: */
    18751866        IssmDouble* cumdz = xNew<IssmDouble>(m-1);
    … …  
    18801871        obp[0]=0.0;
    18811872        for(int i=1;i<m;i++)obp[i]=cumdz[i-1]*d[i-1];
    1882        
     1873
    18831874        // calculate new snow/firn density for:
    18841875        //   snow with densities <= 550 [kg m-3]
    18851876        //   snow with densities > 550 [kg m-3]
    1886                
    1887        
     1877
    18881878        for(int i=0;i<m;i++){
    18891879                switch (denIdx){
    … …  
    20041994        IssmDouble lhf=0.0;
    20051995        IssmDouble EC=0.0;
    2006        
     1996
    20071997        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   turbulentFlux module\n");
    20081998
    … …  
    20292019
    20302020        // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H'
    2031        
     2021
    20322022        // do not allow Ri to exceed 0.19
    20332023        Ri = fmin(Ri, 0.19);
    … …  
    20732063        }
    20742064
    2075 
    20762065        // Latent heat flux [W m-2]
    20772066        lhf = C * L * (eAir - eS) * 0.622 / pAir;

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

    r22448 r23066  
    66#include "../../shared/shared.h"
    77#include "../../toolkits/toolkits.h"
    8 
    98
    109void SmbGradientsx(FemModel* femmodel){/*{{{*/
    … …  
    318317                        smb = smb/yts + anomaly;
    319318
    320 
    321319                        /*Update array accordingly*/
    322320                        smblist[v] = smb;

  • issm/trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp

    r22852 r23066  
    1111                        IssmDouble* PrecipitationReconstructed,IssmDouble* TemperatureReconstructed, IssmDouble* monthlytemperaturesout,
    1212                        IssmDouble* monthlyprecout){
    13  
     13
    1414  IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12];
    1515  IssmDouble deltaTemp;
    16  
     16
    1717  /* Constants */
    1818  // dpermil = 2.4;/*degrees C per mil*/
    19  
     19
    2020  /*Create Delta Temp to be applied to monthly temps and used in precip scaling*/
    2121  deltaTemp = dpermil * (d018+34.83);   
    22    
     22
    2323  for (int imonth = 0; imonth<12; imonth++){
    24    
     24
    2525         if (isTemperatureScaled)monthlytemperaturestmp[imonth] = TemperaturePresentday[imonth] + deltaTemp;
    2626         else{
    … …  
    3131         if (isPrecipScaled)monthlyprectmp[imonth] = PrecipitationPresentday[imonth]*exp((f/dpermil)*deltaTemp);
    3232         else monthlyprectmp[imonth] = PrecipitationReconstructed[imonth];
    33          
     33
    3434    /*Assign output pointer*/
    3535    *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth];

  • issm/trunk-jpl/src/c/shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp

    r19268 r23066  
    1616  IssmDouble tdiffh; 
    1717
    18 
    1918  for (int imonth = 0; imonth<12; imonth++){
    2019    tdiffh = TdiffTime*( TemperaturesLgm[imonth] - TemperaturesPresentday[imonth] );
    … …  
    2928  // printf(" tempera %f\n",monthlytemperaturestmp[1]);
    3029}
    31  

  • issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp

    r18063 r23066  
    1919        if((w0==w1)||(w1==w2)||(w0==w2))
    2020                _error_("Error: equal ordinates in DrainageFunctionWaterfraction -> division by zero. Abort");
    21        
     21
    2222        if(waterfraction<=w0)
    2323                Dret=D0;

  • issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp

    r21826 r23066  
    6666        /*Get norm of epsprime*/
    6767        epsprime_norm = sqrt(epsprime[0]*epsprime[0] + epsprime[1]*epsprime[1] + epsprime[2]*epsprime[2]);
    68        
     68
    6969        /*Assign output pointers*/
    7070        *pepsprime_norm=epsprime_norm;

  • issm/trunk-jpl/src/c/shared/Elements/Paterson.cpp

    r17457 r23066  
    1515        /*Switch to celsius from Kelvin: */
    1616        T=temperature-273.15;
    17 
    1817
    1918        if(T<=-45.0){

  • issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp

    r22448 r23066  
    8383  for (iqj = 0; iqj < 12; iqj++){
    8484    imonth =  ismon[iqj];
    85    
     85
    8686    /*********compute lapse rate ****************/
    8787    st=(s-s0t)/1000.;
    8888    rtlaps=TdiffTime*rlapslgm + (1.-TdiffTime)*rlaps; // lapse rate
    89    
     89
    9090    /*********compute Surface temperature *******/
    9191    monthlytemperatures[imonth]=monthlytemperatures[imonth] - rtlaps *max(st,sealevTime*0.001);
    … …  
    9999    else {
    100100      pd = 0.;}
    101    
     101
    102102    /******exp des/elev precip reduction*******/
    103103    sp=(s-s0p)/1000.-deselcut; // deselev effect is wrt chng in topo
    104104    if (sp>0.0){q = exp(-desfac*sp);}
    105105    else {q = 1.0;}
    106    
     106
    107107    qmt= qmt + monthlyprec[imonth]*sconv;  //*sconv to convert in m of ice equivalent per month 
    108108    qmpt= q*monthlyprec[imonth]*sconv;
    … …  
    114114    // gaussian=T_m, so ndd=-(Tsurf-pdd)
    115115    if (iqj>5 && iqj<9){ Tsum=Tsum+tstar;}
    116    
     116
    117117    if (tstar >= siglim) {pdd = pdd + tstar*deltm;}
    118118    else if (tstar> -siglim){
    … …  
    121121      frzndd = frzndd - (tstar-pddsig)*deltm;}
    122122    else{frzndd = frzndd - tstar*deltm; }
    123    
     123
    124124    /*Assign output pointer*/
    125125    *(monthlytemperatures+imonth) = monthlytemperatures[imonth];
    … …  
    150150          }
    151151  }
    152  
     152
    153153  /***** determine PDD factors *****/
    154154  if(Tsum< -1.) {
    … …  
    166166  snwmf=0.95*snwmf;
    167167  smf=0.95*smf;
    168  
     168
    169169  /*****  compute PDD ablation and refreezing *****/
    170170  pddt = pdd *365;
    171171  snwm = snwmf*pddt;       // snow that could have been melted in a year
    172172  hmx2 = min(h,dfrz);   // refreeze active layer max depth: dfrz
    173  
     173
    174174  if(snwm < saccu) {
    175175    water=prect-saccu + snwm; //water=rain + snowmelt

  • issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp

    r21216 r23066  
    6363                dtemp = &dtemp_static[0];
    6464        }
    65        
     65
    6666        /*  perform the matrix triple product in the order that minimizes the
    6767                 number of multiplies and the temporary space used, noting that
    … …  
    335335        Matrix2x2Determinant(&det,A);
    336336        if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
    337        
     337
    338338        /*Multiplication is faster than divsion, so we multiply by the reciprocal*/
    339339        det_reciprocal = 1./det; 
    … …  
    427427        *pvy      = vy;
    428428
    429 
    430429}/*}}}*/
    431430
    … …  
    577576
    578577void newcell(IssmDouble** pcell, IssmDouble newvalue, bool top, int m){  /*{{{*/
    579    
     578
    580579    IssmDouble* cell=NULL;
    581580    IssmDouble* dummy=NULL;
    582    
     581
    583582    /*recover pointer: */
    584583    cell=*pcell;
    585    
     584
    586585    /*reallocate:*/
    587586    dummy=xNew<IssmDouble>(m+1);
    588    
     587
    589588        /*copy data:*/
    590589    if(top){
    … …  
    596595        for(int i=0;i<m;i++)dummy[i]=cell[i];
    597596    }
    598    
     597
    599598    /*delete and reassign: */
    600599    xDelete<IssmDouble>(cell); cell=dummy;
    601    
     600
    602601    /*assign output pointer:*/
    603602    *pcell=cell;
    … …  
    615614        /*input: */
    616615        IssmDouble* cell=*pcell;
    617        
     616
    618617        /*output: */
    619618        IssmDouble* newcell=xNew<IssmDouble>(nind);
    … …  
    622621                newcell[i]=cell[indices[i]];
    623622        }
    624        
     623
    625624        /*free allocation:*/
    626625        xDelete<IssmDouble>(cell);
    … …  
    633632        /*input: */
    634633        IssmDouble* cell=*pcell;
    635        
     634
    636635        /*output: */
    637636        IssmDouble* newcell=xNew<IssmDouble>(m+1);
    … …  
    641640        newcell[i+1]=scale* cell[i];
    642641        for(int j=i+2;j<m+1;j++)newcell[j]=cell[j-1];
    643        
     642
    644643        /*free allocation:*/
    645644        xDelete<IssmDouble>(cell);
    … …  
    662661        }
    663662        va_end ( arguments );                 
    664        
     663
    665664        _printf_("Echo of cell: \n");
    666665        for(int i=0;i<m;i++){

  • issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp

    r18137 r23066  
    6868                _printf0_(" x = "<<setw(9)<<xmin<<" | ");
    6969                fxmin = (*g)(&G,X0,usr); if(xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
    70                
     70
    7171                /*Get f(xmax)*/
    7272                _printf0_(" x = "<<setw(9)<<xmax<<" | ");
    … …  
    244244                J[n]=fxbest;
    245245        }
    246        
     246
    247247        /*return*/
    248248        xDelete<IssmDouble>(X);

  • issm/trunk-jpl/src/c/shared/Numerics/legendre.cpp

    r20033 r23066  
    121121
    122122          for i = 2 : n
    123          
     123
    124124                v(:,i+1) = ( ( 2 * i - 1 ) * x .* v(:,i)   ...
    125125                                        -  (     i - 1 ) *    v(:,i-1) ) ...
    126126                                        /  (     i     );
    127          
     127
    128128          end
    129129          }}}  */
    … …  
    240240                                           polynomials of order 0 through N.
    241241        }}}*/
    242        
     242
    243243        int i,j;
    244244
    … …  
    254254        for ( i = 0; i < m; i++ ) {
    255255                for ( j = 2; j <= n; j++ ) {
    256                        
     256
    257257                        v[j+i*(n+1)] = ( ( IssmDouble ) ( 2 * j - 1 ) * x[i] * v[(j-1)+i*(n+1)]
    258258                                        - ( IssmDouble ) (     j - 1 ) *        v[(j-2)+i*(n+1)] )

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp

    r22551 r23066  
    185185        /*Initial guess*/
    186186        VecZeroEntries(udot);
    187        
     187
    188188        /*Richardson's iterations*/
    189189        for(int i=0;i<5;i++){
    … …  
    402402        CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
    403403        delete uf;
    404        
     404
    405405        /*Go solve lower order solution*/
    406406        femmodel->profiler->Start(SOLVER);

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp

    r22551 r23066  
    2020        int  configuration_type;
    2121        IssmDouble solver_residue_threshold;
    22        
     22
    2323        /*solver convergence test: */
    2424        IssmDouble nKUF;
    … …  
    3939        Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters);
    4040        femmodel->profiler->Stop(SOLVER);
    41        
     41
    4242        /*Check that the solver converged nicely: */
    43                
     43
    4444        //compute KUF = KU - F = K*U - F
    4545        KU=uf->Duplicate(); Kff->MatMult(uf,KU);
    … …  
    5353        if(!xIsNan(solver_residue_threshold) && solver_residue>solver_residue_threshold)_error_("   solver residue too high!: norm(KU-F)/norm(F)=" << solver_residue << "\n");
    5454
    55 
    5655        //clean up
    5756        delete KU; delete KUF;
    … …  
    6362        //        }
    6463        //#endif
    65        
     64
    6665        Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
    6766        InputUpdateFromSolutionx(femmodel,ug);

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp

    r22551 r23066  
    6161
    6262        count=1;
    63        
     63
    6464        for(;;){
    6565                delete tf_old;tf_old=tf;

  • issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp

    r23044 r23066  
    178178#endif
    179179
    180 
    181180        recvcounts=xNew<int>(num_procs);
    182181        displs=xNew<int>(num_procs);
    … …  
    262261        rhs=xNew<IssmDouble>(pf_M);
    263262#endif
    264 
    265263
    266264        recvcounts=xNew<int>(num_procs);
Note: See TracChangeset for help on using the changeset viewer.