Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 23066) @@ -1153,5 +1153,5 @@ /*Clean up and return*/ xDelete(responses); - + }/*}}}*/ void AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector* gradient,int control_index){/*{{{*/ @@ -2112,5 +2112,4 @@ IssmDouble *xyz_list= NULL; - /*Fetch number of vertices for this finite element*/ int numvertices = basalelement->GetNumberOfVertices(); @@ -2131,6 +2130,4 @@ Input* adjointx_input = basalelement->GetInput(AdjointxEnum); _assert_(adjointx_input); Input* adjointy_input = basalelement->GetInput(AdjointyEnum); _assert_(adjointy_input); - - IssmDouble q_exp; @@ -2182,9 +2179,9 @@ Chi = vmag/(pow(C_param,n)*pow(Neff,n)*As); Gamma = (Chi/(1.+alpha*pow(Chi,q_exp))); - + Uder =Neff*C_param/(vmag*vmag*n) * (Gamma-alpha*q_exp*pow(Chi,q_exp-1.)*Gamma*Gamma* pow(Gamma,(1.-n)/n) - n* pow(Gamma,1./n)); - + /*Build gradient vector (actually -dJ/dD): */ for(int i=0;iNewGauss(2); Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 23066) @@ -126,5 +126,5 @@ /*Add input*/ element->AddInput(DamageFEnum,f,element->GetElementType()); - + /*Clean up and return*/ xDelete(f); @@ -170,5 +170,5 @@ for (int i=0;iGaussNode(element->GetElementType(),i); - + eps_xx_input->GetInputValue(&eps_xx,gauss); eps_xy_input->GetInputValue(&eps_xy,gauss); @@ -177,5 +177,5 @@ n_input->GetInputValue(&n,gauss); damage_input->GetInputValue(&damage,gauss); - + /*Calculate principal effective strain rates*/ eps1=(eps_xx+eps_yy)/2.+sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2)); @@ -195,5 +195,5 @@ /*Add input*/ element->AddInput(DamageFEnum,f,element->GetElementType()); - + /*Clean up and return*/ xDelete(f); @@ -262,5 +262,5 @@ for (int i=0;iGaussNode(element->GetElementType(),i); - + damage_input->GetInputValue(&damage,gauss); tau_xx_input->GetInputValue(&tau_xx,gauss); @@ -314,5 +314,5 @@ /*Add input*/ element->AddInput(DamageFEnum,f,element->GetElementType()); - + /*Clean up and return*/ xDelete(f); @@ -375,5 +375,5 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); - + vx_input->GetInputValue(&vx,gauss); vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); @@ -538,5 +538,4 @@ } - /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(2); Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 23066) @@ -103,5 +103,5 @@ int frictionlaw,basalforcing_model,materialstype; int FrictionCoupling; - + /*Now, is the model 3d? otherwise, do nothing: */ if(iomodel->domaintype==Domain2DhorizontalEnum)return; @@ -189,5 +189,5 @@ _error_("not supported"); } - + /*Friction law variables*/ switch(frictionlaw){ @@ -805,5 +805,5 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); - + /*viscous dissipation*/ element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input); @@ -823,5 +823,5 @@ pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss); d2pressure=0.; // for linear elements, 2nd derivative is zero - + d1H_d1P=0.; for(i=0;i<3;i++) d1H_d1P+=d1enthalpypicard[i]*d1pressure[i]; @@ -1056,5 +1056,5 @@ int* basalnodeindices=NULL; Element* element= NULL; - + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); @@ -1388,5 +1388,5 @@ gauss->GaussNode(element->GetElementType(),indices[i]); gaussup->GaussNode(element->GetElementType(),indicesup[i]); - + enthalpy_input->GetInputValue(&enthalpy,gauss); enthalpy_input->GetInputValue(&enthalpyup,gaussup); Index: /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 23066) @@ -40,5 +40,5 @@ IssmDouble* love_h=NULL; IssmDouble* love_l=NULL; - + IssmDouble* U_elastic = NULL; IssmDouble* U_elastic_local = NULL; @@ -69,5 +69,5 @@ U_elastic=xNew(M); H_elastic=xNew(M); - + /*compute combined legendre + love number (elastic green function:*/ m=DetermineLocalSize(M,IssmComm::GetComm()); @@ -96,5 +96,5 @@ deltalove_U = (love_h[n]-love_h[nl-1]); - + /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */ if(n==0){ @@ -199,5 +199,5 @@ /*Default, do nothing*/ return; - + }/*}}}*/ void EsaAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 23066) @@ -145,5 +145,5 @@ GetB(B,workelement,xyz_list,gauss,dim); GetBprime(Bprime,workelement,xyz_list,gauss,dim); - + D_scalar=gauss->weight*Jdet; @@ -174,5 +174,5 @@ else for(i=0;iFindParam(&extvar_enum, ExtrapolationVariableEnum); - + Input* active_input=element->GetInput(IceMaskNodeActivationEnum); _assert_(active_input); Input* extvar_input=element->GetInput(extvar_enum); _assert_(extvar_input); Index: /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp (revision 23066) @@ -199,6 +199,4 @@ } - - return; }/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 23066) @@ -260,5 +260,4 @@ basis,1,numnodes,0, &Ke->values[0],1); - /*Transfer EPL part*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 23066) @@ -213,5 +213,4 @@ return NULL; } - /*Intermediaries */ @@ -504,5 +503,4 @@ } - /*Fetch number of nodes for this finite element*/ int numnodes = basalelement->GetNumberOfNodes(); Index: /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 23066) @@ -295,5 +295,5 @@ n_input->GetInputValue(&n,gauss); A=pow(B,-n); - + /*Compute beta term*/ if(gap0.); - + for(int i=0;ivalues[i]+=Jdet*gauss->weight* ( @@ -398,5 +398,5 @@ /*Calculate effective pressure*/ eff_pressure[i] = rho_ice*g*thickness[i] - rho_water*g*(values[i]-bed[i]); - + if(xIsNan(values[i])) _error_("NaN found in solution vector"); if(xIsInf(values[i])) _error_("Inf found in solution vector"); @@ -444,5 +444,5 @@ reynolds_input->GetInputAverage(&reynolds); gap_input->GetInputAverage(&gap); - + /*Compute conductivity*/ IssmDouble conductivity = pow(gap,3)*g/(12.*NU*(1+OMEGA*reynolds)); @@ -453,5 +453,4 @@ }/*}}}*/ void HydrologyShaktiAnalysis::UpdateGapHeight(FemModel* femmodel){/*{{{*/ - for(int j=0;jelements->Size();j++){ @@ -547,5 +546,5 @@ IssmDouble pressure_water = rho_water*g*(head-bed); if(pressure_water>pressure_ice) pressure_water = pressure_ice; - + /* Compute change in sensible heat due to changes in pressure melting point*/ dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]); @@ -581,5 +580,5 @@ /*Add new gap as an input*/ element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum); - + /*Divide by connectivity, add basal flux as an input*/ q = q/totalweights; Index: /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 23066) @@ -379,5 +379,4 @@ }/*}}}*/ - /*Needed changes to switch to the Johnson formulation*//*{{{*/ /*All the changes are to be done in the velocity computation. Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 23066) @@ -43,5 +43,5 @@ } } - + iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum); @@ -177,5 +177,5 @@ /*Calving threshold*/ - + /*Fetch number of nodes and dof for this finite element*/ int numnodes = basalelement->GetNumberOfNodes(); @@ -326,5 +326,5 @@ } break; - + case CalvingLevermannEnum: calvingratex_input->GetInputValue(&c[0],gauss); @@ -357,5 +357,5 @@ } break; - + case CalvingHabEnum: lsf_slopex_input->GetInputValue(&dlsf[0],gauss); @@ -378,5 +378,5 @@ } break; - + case CalvingCrevasseDepthEnum: lsf_slopex_input->GetInputValue(&dlsf[0],gauss); @@ -385,5 +385,5 @@ if(groundedice<0) meltingrate = 0.; - + norm_dlsf=0.; for(i=0;iIsOnBase()) return NULL; Element* basalelement = element->SpawnBasalElement(); @@ -525,5 +525,5 @@ IssmDouble lsf; IssmDouble* xyz_list = NULL; - + /*Fetch number of nodes and dof for this finite element*/ int numnodes = basalelement->GetNumberOfNodes(); @@ -532,5 +532,5 @@ ElementVector* pe = basalelement->NewElementVector(); basalelement->FindParam(&dt,TimesteppingTimeStepEnum); - + if(dt!=0.){ /*Initialize basis vector*/ @@ -623,5 +623,5 @@ // d=|a x b|/|b| // with a=q-s0, b=s1-s0 - + /* Intermediaries */ const int dim=2; @@ -634,5 +634,5 @@ b[i]=s1[i]-s0[i]; } - + norm_b=0.; for(i=0;iparameters->FindParam(&calvinglaw,CalvingLawEnum); @@ -703,5 +703,5 @@ } } - + if(calvinglaw==CalvingHabEnum){ @@ -709,9 +709,9 @@ femmodel->elements->InputDuplicate(MaskIceLevelsetEnum, DistanceToCalvingfrontEnum); femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum); - + /*Loop over all elements of this partition*/ for(int i=0;ielements->Size();i++){ Element* element = xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - + rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); @@ -744,11 +744,11 @@ } } - + if(calvinglaw==CalvingCrevasseDepthEnum){ - + int nflipped,local_nflipped; Vector* vec_constraint_nodes = NULL; IssmDouble* constraint_nodes = NULL; - + /*Get the DistanceToCalvingfront*/ femmodel->elements->InputDuplicate(MaskIceLevelsetEnum,DistanceToCalvingfrontEnum); @@ -842,5 +842,5 @@ gauss->GaussNode(element->GetElementType(),in); Node* node=element->GetNode(in); - + if(constraint_nodes[node->Sid()]>0.){ node->ApplyConstraint(0,+1.); Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 23066) @@ -165,5 +165,5 @@ iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum); iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum); - + /*Initialize cumdeltalthickness input*/ InputUpdateFromConstantx(elements,0.,SealevelriseCumDeltathicknessEnum); @@ -203,5 +203,5 @@ iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum); } - + }/*}}}*/ void MasstransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ @@ -221,6 +221,4 @@ if(numoutputs)parameters->AddObject(new StringArrayParam(MasstransportRequestedOutputsEnum,requestedoutputs,numoutputs)); iomodel->DeleteData(&requestedoutputs,numoutputs,"md.masstransport.requested_outputs"); - - }/*}}}*/ @@ -805,5 +803,5 @@ xDelete(sealevel); xDelete(bed); - + xDelete(doflist); if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 23066) @@ -69,5 +69,5 @@ IssmDouble* love_k=NULL; IssmDouble* love_l=NULL; - + bool elastic=false; IssmDouble* G_elastic = NULL; @@ -132,5 +132,4 @@ #endif - /*compute combined legendre + love number (elastic green function:*/ m=DetermineLocalSize(M,IssmComm::GetComm()); @@ -168,5 +167,5 @@ deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]); deltalove_U = (love_h[n]-love_h[nl-1]); - + /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */ if(n==0){ @@ -230,5 +229,5 @@ xDelete(H_elastic_local); } - + /*Transitions: */ iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions"); @@ -250,5 +249,4 @@ iomodel->DeleteData(&requestedoutputs,numoutputs,"md.slr.requested_outputs"); - }/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23066) @@ -19,8 +19,8 @@ }/*}}}*/ void SmbAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ - + int smb_model; bool isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled; - + /*Update elements: */ int counter=0; @@ -32,8 +32,8 @@ } } - + /*Figure out smb model: */ iomodel->FindConstant(&smb_model,"md.smb.model"); - + switch(smb_model){ case SMBforcingEnum: @@ -142,6 +142,4 @@ _error_("Surface mass balance model "<AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum)); - + iomodel->FindConstant(&smb_model,"md.smb.model"); iomodel->FindConstant(&interp,"md.timestepping.interp_forcings"); - + switch(smb_model){ case SMBforcingEnum: @@ -198,5 +196,5 @@ parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,M)); iomodel->DeleteData(temp,"md.smb.Pfac"); - + iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2); parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,M)); @@ -266,5 +264,5 @@ /*Figure out smb model: */ femmodel->parameters->FindParam(&smb_model,SmbEnum); - + /*branch to correct module*/ switch(smb_model){ Index: /issm/trunk-jpl/src/c/analyses/SmoothAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmoothAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/SmoothAnalysis.cpp (revision 23066) @@ -153,5 +153,4 @@ } - /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(2); @@ -161,5 +160,4 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); - switch(input_enum){ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23066) @@ -2736,5 +2736,4 @@ } - /*Clean-up*/ xDelete(dbasis); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 23066) @@ -29,5 +29,4 @@ else iscoupling = false; - /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/ @@ -159,5 +158,4 @@ ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixBase(Element* element){/*{{{*/ - if(!element->IsOnBase()) return NULL; @@ -199,5 +197,4 @@ }/*}}}*/ ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/ - if(!element->IsOnSurface()) return NULL; Index: /issm/trunk-jpl/src/c/bamg/BamgQuadtree.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgQuadtree.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/bamg/BamgQuadtree.cpp (revision 23066) @@ -285,5 +285,4 @@ int yiv = b->v[k]->i.y; - int h0 = Norm(xi2,xiv,yi2,yiv); Index: /issm/trunk-jpl/src/c/bamg/Mesh.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/Mesh.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/bamg/Mesh.cpp (revision 23066) @@ -501,5 +501,5 @@ int connectivitymax_1=0; int verbose=0; - + /*Check needed pointer*/ _assert_(bamgmesh); @@ -1034,5 +1034,5 @@ // s0 tt2 s1 //------------------------------- - + /*Intermediaries*/ Triangle* tt[3]; //the three triangles @@ -1158,8 +1158,8 @@ int verbose=0; double lminaniso = 1./ (Max(hminaniso*hminaniso,1e-100)); - + //Get options if(bamgopts) verbose=bamgopts->verbose; - + //display info if (verbose > 1) _printf_(" BoundAnisotropy by " << anisomax << "\n"); @@ -1854,5 +1854,5 @@ /*Check pointer*/ _assert_(bamgopts); - + /*Recover options*/ verbose=bamgopts->verbose; @@ -2728,5 +2728,5 @@ if (verbose>4) _printf_(" Prime Number = "<4) _printf_(" k0 = "< (BAMG v1.01, Metric.cpp/MaxSubDivision)*/ - + /*Intermediaries*/ int verbose=0; @@ -3740,8 +3740,8 @@ void Mesh::SmoothingVertex(BamgOpts* bamgopts,int nbiter,double omega ) { /*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/SmoothingVertex)*/ - + /*Intermediaries*/ int verbose=0; - + /*Get options*/ if(bamgopts) verbose=bamgopts->verbose; @@ -3785,5 +3785,5 @@ void Mesh::SmoothMetric(BamgOpts* bamgopts,double raisonmax) { /*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Metric.cpp/SmoothMetric)*/ - + /*Intermediaries*/ int verbose=0; @@ -4087,5 +4087,5 @@ int i; Metric M1(1); - + /*Get options*/ if(bamgopts) verbose=bamgopts->verbose; Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 23066) @@ -101,5 +101,5 @@ /*Mesh refinement methods*/ void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int** pdatalist,double** pxylist,int** pelementslist){/*{{{*/ - + /*IMPORTANT! pelementslist are in Matlab indexing*/ /*NEOPZ works only in C indexing*/ @@ -118,8 +118,8 @@ /*Intermediaries*/ bool verbose=VerboseSolution(); - + /*Execute refinement*/ this->RefineMeshOneLevel(verbose,gl_distance,if_distance,deviatoricerror,thicknesserror); - + /*Get new geometric mesh in ISSM data structure*/ this->GetMesh(this->previousmesh,pdatalist,pxylist,pelementslist); @@ -130,5 +130,5 @@ /*}}}*/ void AdaptiveMeshRefinement::RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror){/*{{{*/ - + /*Intermediaries*/ int nelem =-1; @@ -204,5 +204,5 @@ } /*}}}*/ - + /*Now, detele the special elements{{{*/ if(this->refinement_type==1) this->DeleteSpecialElements(verbose,gmesh); @@ -216,10 +216,10 @@ } /*}}}*/ - + /*Unrefinement process: loop over level of refinements{{{*/ if(verbose) _printf_("\tunrefinement process...\n"); if(verbose) _printf_("\ttotal: "); count=0; - + nelem=gmesh->NElements();//must keep here for(int i=0;iBuildConnectivity(); /*}}}*/ - + /*Refinement process: loop over level of refinements{{{*/ if(verbose) _printf_("\trefinement process (level max = "<level_max<<")\n"); @@ -320,12 +320,12 @@ * */ if(!geoel) _error_("geoel is NULL!\n"); - + /*Output*/ int type=0; int count=0; - + /*Intermediaries*/ TPZVec sons; - + /*Loop over neighboors (sides 3, 4 and 5)*/ for(int j=3;j<6;j++){ @@ -335,9 +335,9 @@ if(sons.size()>4) count++; //if neighbour's level is > element level+1 } - + /*Verify and return*/ if(count>1) type=2; else type=count; - + return type; } @@ -391,12 +391,12 @@ /*}}}*/ void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/ - + /*Now, insert special elements to avoid hanging nodes*/ if(verbose) _printf_("\trefining to avoid hanging nodes (total: "); - + /*Intermediaries*/ int nelem=-1; int count= 1; - + while(count>0){ nelem=gmesh->NElements();//must keep here @@ -468,11 +468,11 @@ this->index2sid.clear(); this->index2sid.resize(gmesh->NElements()); this->sid2index.clear(); - + /*Fill in the vertex_index2sid vector with non usual index value*/ for(int i=0;iNNodes();i++) vertex_index2sid[i]=-1; - + /*Fill in the this->index2sid vector with non usual index value*/ for(int i=0;iNElements();i++) this->index2sid[i]=-1; - + /*Get elements without sons and fill in the vertex_index2sid with used vertices (indexes) */ sid=0; @@ -511,5 +511,5 @@ newmeshXY[2*sid+1] = coords[1]; // Y } - + for(int i=0;isid2index.size();i++){//over the sid (fill the ISSM elements) for(int j=0;jGetNumberOfNodes();j++) { @@ -533,5 +533,5 @@ detJ=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya); - + /*verify and swap, if necessary*/ if(detJ<0) { @@ -545,5 +545,5 @@ *pxy = newmeshXY; *pelements = newelements; - + /*Cleanup*/ xDelete(vertex_index2sid); @@ -554,5 +554,5 @@ /* IMPORTANT! elements come in Matlab indexing NEOPZ works only in C indexing*/ - + if(nvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n"); if(nelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n"); @@ -561,5 +561,5 @@ /*Verify and creating initial mesh*/ if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!"); - + this->fathermesh = new TPZGeoMesh(); this->fathermesh->NodeVec().Resize(nvertices); @@ -576,5 +576,5 @@ this->fathermesh->NodeVec()[i].SetNodeId(i); } - + /*Generate the elements*/ long index; @@ -604,8 +604,8 @@ /*}}}*/ TPZGeoMesh* AdaptiveMeshRefinement::CreateRefPatternMesh(TPZGeoMesh* gmesh){/*{{{*/ - + TPZGeoMesh *newgmesh = new TPZGeoMesh(); newgmesh->CleanUp(); - + int nnodes = gmesh->NNodes(); int nelem = gmesh->NElements(); @@ -617,9 +617,9 @@ newgmesh->NodeVec().Resize(nnodes); for(int i=0;iNodeVec()[i] = gmesh->NodeVec()[i]; - + //elements for(int i=0;iElement(i); - + if(!geoel){ index=newgmesh->ElementVec().AllocateNewElement(); @@ -627,13 +627,13 @@ continue; } - + TPZManVector elem(3,0); for(int j=0;j<3;j++) elem[j] = geoel->NodeIndex(j); - + newgmesh->CreateGeoElement(ETriangle,elem,mat,index,reftype); newgmesh->ElementVec()[index]->SetId(geoel->Id()); - + TPZGeoElRefPattern* newgeoel = dynamic_cast*>(newgmesh->ElementVec()[index]); - + //old neighbourhood const int nsides = TPZGeoTriangle::NSides; @@ -654,5 +654,5 @@ } } - + //inserting in new element for(int s = 0; s < nsides; s++){ @@ -668,15 +668,15 @@ tempEl->SetNeighbour(byside, tempSide); } - + long fatherindex = geoel->FatherIndex(); if(fatherindex>-1) newgeoel->SetFather(fatherindex); - + if(!geoel->HasSubElement()) continue; - + int nsons = geoel->NSubElements(); TPZAutoPointer ref = gRefDBase.GetUniformRefPattern(ETriangle); newgeoel->SetRefPattern(ref); - + for(int j=0;jSubElement(j); @@ -690,5 +690,5 @@ /*Now, build connectivities*/ newgmesh->BuildConnectivity(); - + return newgmesh; } @@ -705,5 +705,5 @@ /*}}}*/ void AdaptiveMeshRefinement::PrintGMeshVTK(TPZGeoMesh* gmesh,std::ofstream &file,bool matColor){/*{{{*/ - + file.clear(); long nelements = gmesh->NElements(); @@ -731,8 +731,8 @@ MElementType elt = gel->Type(); int elNnodes = MElementType_NNodes(elt); - + size += (1+elNnodes); connectivity << elNnodes; - + for(int t = 0; t < elNnodes; t++) { @@ -743,13 +743,13 @@ } node << std::endl; - + actualNode++; connectivity << " " << actualNode; } connectivity << std::endl; - + int elType = this->GetVTK_ElType(gel); type << elType << std::endl; - + if(matColor == true) { @@ -760,5 +760,5 @@ material << gel->Index() << std::endl; } - + nVALIDelements++; } @@ -790,7 +790,7 @@ /*}}}*/ int AdaptiveMeshRefinement::GetVTK_ElType(TPZGeoEl * gel){/*{{{*/ - + MElementType pzElType = gel->Type(); - + int elType = -1; switch (pzElType) @@ -849,7 +849,6 @@ DebugStop(); } - + return elType; } /*}}}*/ - Index: /issm/trunk-jpl/src/c/classes/AmrBamg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 23066) @@ -34,5 +34,5 @@ this->deviatoricerror_groupthreshold = -1; this->deviatoricerror_maximum = -1; - + /*Geometry and mesh as NULL*/ this->geometry = NULL; @@ -137,5 +137,5 @@ this->options->hmaxVerticesSize[0] = 0; this->options->hmaxVerticesSize[1] = 0; - + /*verify if the metric will be reseted or not*/ if(this->keepmetric==0){ @@ -155,13 +155,13 @@ IssmDouble *xylist= xNew(nbv*2); int* elementslist = xNew(nbt*3); - + datalist[0] = nbv; datalist[1] = nbt; - + for(int i=0;iVertices[i*3+0]; xylist[2*i+1] = meshout->Vertices[i*3+1]; } - + for(int i=0;i(meshout->Triangles[4*i+0]); @@ -174,5 +174,5 @@ *pxylist = xylist; *pelementslist = elementslist; - + /*Cleanup and return*/ delete geomout; @@ -181,5 +181,5 @@ if(!this->options) _error_("AmrBamg->options is NULL!"); - + this->options->hmin = hmin_in; this->options->hmax = hmax_in; Index: /issm/trunk-jpl/src/c/classes/Cfdragcoeffabsgrad.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Cfdragcoeffabsgrad.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Cfdragcoeffabsgrad.cpp (revision 23066) @@ -23,5 +23,5 @@ #include "../classes/gauss/Gauss.h" /*}}}*/ - + /*Cfdragcoeffabsgrad constructors, destructors :*/ Cfdragcoeffabsgrad::Cfdragcoeffabsgrad(){/*{{{*/ @@ -40,5 +40,5 @@ this->definitionenum=in_definitionenum; - + this->name = xNew(strlen(in_name)+1); xMemCpy(this->name,in_name,strlen(in_name)+1); @@ -46,5 +46,5 @@ this->weights_enum=in_weights_enum; this->timepassedflag=in_timepassedflag; - + this->misfit=0; this->lock=0; @@ -102,5 +102,5 @@ /*diverse: */ IssmDouble time; - + /*recover time parameters: */ int i; @@ -116,5 +116,5 @@ ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); J=J_sum; - + timepassedflag = true; return J; @@ -183,3 +183,2 @@ return Jelem; }/*}}}*/ - Index: /issm/trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp (revision 23066) @@ -23,5 +23,5 @@ #include "../classes/gauss/Gauss.h" /*}}}*/ - + /*Cfsurfacelogvel constructors, destructors :*/ Cfsurfacelogvel::Cfsurfacelogvel(){/*{{{*/ @@ -40,5 +40,5 @@ this->definitionenum=in_definitionenum; - + this->name = xNew(strlen(in_name)+1); xMemCpy(this->name,in_name,strlen(in_name)+1); @@ -46,5 +46,5 @@ this->datatime=in_datatime; this->timepassedflag=in_timepassedflag; - + this->misfit=0; this->lock=0; @@ -100,8 +100,8 @@ /*}}}*/ IssmDouble Cfsurfacelogvel::Response(FemModel* femmodel){/*{{{*/ - + /*diverse: */ IssmDouble time; - + /*recover time parameters: */ femmodel->parameters->FindParam(&time,TimeEnum); @@ -112,5 +112,5 @@ IssmDouble J=0.; IssmDouble J_sum=0.; - + if(datatime<=time && !timepassedflag){ for(i=0;ielements->Size();i++){ @@ -122,5 +122,5 @@ ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); J=J_sum; - + timepassedflag = true; return J; @@ -139,5 +139,5 @@ IssmDouble vx,vy,vxobs,vyobs,weight; IssmDouble* xyz_list = NULL; - + /*Get basal element*/ if(!element->IsOnSurface()) return 0.; @@ -160,5 +160,5 @@ /* Get node coordinates*/ topelement->GetVerticesCoordinates(&xyz_list); - + /*Get model values*/ Input* vx_input =topelement->GetInput(VxEnum); _assert_(vx_input); @@ -174,5 +174,4 @@ if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do"); datasetinput = (DatasetInput*)tempinput; - /* Start looping on the number of gaussian points: */ @@ -222,3 +221,2 @@ return Jelem; }/*}}}*/ - Index: /issm/trunk-jpl/src/c/classes/Cfsurfacesquare.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Cfsurfacesquare.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Cfsurfacesquare.cpp (revision 23066) @@ -23,5 +23,5 @@ #include "../classes/gauss/Gauss.h" /*}}}*/ - + /*Cfsurfacesquare constructors, destructors :*/ Cfsurfacesquare::Cfsurfacesquare(){/*{{{*/ @@ -43,5 +43,5 @@ this->definitionenum=in_definitionenum; - + this->name = xNew(strlen(in_name)+1); xMemCpy(this->name,in_name,strlen(in_name)+1); @@ -52,5 +52,5 @@ this->datatime=in_datatime; this->timepassedflag=in_timepassedflag; - + this->misfit=0; this->lock=0; @@ -111,5 +111,5 @@ /*diverse: */ IssmDouble time; - + /*recover time parameters: */ femmodel->parameters->FindParam(&time,TimeEnum); @@ -130,5 +130,5 @@ ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); J=J_sum; - + timepassedflag = true; return J; @@ -172,5 +172,5 @@ /*Get input if it already exists*/ Input* tempinput = topelement->GetInput(definitionenum); - + /*Cast it to a Datasetinput*/ if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do"); @@ -214,3 +214,2 @@ return Jelem; }/*}}}*/ - Index: /issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp (revision 23066) @@ -37,5 +37,5 @@ /*Object virtual functions definitions:*/ Object* SpcStatic::copy() {/*{{{*/ - + SpcStatic* spcstat = new SpcStatic(*this); Index: /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 23066) @@ -20,5 +20,5 @@ int world_rank; ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&world_rank); - + /*Build an femmodel if you are a slave, using the corresponding communicator:*/ if(world_rank!=0){ @@ -33,5 +33,5 @@ int world_rank; ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&world_rank); - + if(world_rank!=0){ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23066) @@ -1105,5 +1105,4 @@ void Element::GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug){/*{{{*/ - /*Get number of nodes for this element*/ int numnodes = this->GetNumberOfNodes(); @@ -1121,5 +1120,4 @@ if(ug[nodes[i]->GetDof(0,GsetEnum)] > input_max) input_max = ug[nodes[i]->GetDof(0,GsetEnum)]; } - /*Second loop to reassign min and max with local extrema*/ @@ -1389,5 +1387,5 @@ _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); } - + /*Clean up*/ xDelete(doflist); @@ -1476,5 +1474,4 @@ parameters->FindParam(&num_controls,InversionNumControlParametersEnum); parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); - /*Get number of vertices*/ @@ -1808,5 +1805,4 @@ } - /*Branch on type of vector: nodal or elementary: */ if(vector_type==1){ //nodal vector @@ -2067,5 +2063,5 @@ this->GetInputListOnVertices(deepwaterel,BasalforcingsDeepwaterElevationEnum); this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum); - + for(int i=0;iupperwaterel[i]) values[i]=0; @@ -2870,5 +2866,4 @@ if (!IsOnSurface()) return; - /*Retrieve material properties and parameters:{{{ */ rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum); @@ -3072,5 +3067,4 @@ 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()); - /*Distribution of absorbed short wave radation with depth:*/ if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid()); Index: /issm/trunk-jpl/src/c/classes/Elements/ElementHook.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/ElementHook.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Elements/ElementHook.cpp (revision 23066) @@ -104,5 +104,5 @@ if(marshall_direction==MARSHALLING_BACKWARD){ - + if (!hnodes_null)this->hnodes = new Hook*[numanalyses]; else this->hnodes=NULL; @@ -140,5 +140,5 @@ _printf_(" ElementHook DeepEcho:\n"); _printf_(" numanalyses : "<< this->numanalyses <<"\n"); - + _printf_(" hnodes:\n"); if(hnodes){ @@ -149,9 +149,9 @@ } else _printf_(" hnodes = NULL\n"); - + _printf_(" hvertices:\n"); if(hvertices) hvertices->DeepEcho(); else _printf_(" hvertices = NULL\n"); - + _printf_(" hmaterial:\n"); if(hmaterial) hmaterial->DeepEcho(); @@ -170,8 +170,8 @@ /*}}}*/ void ElementHook::Echo(){/*{{{*/ - + _printf_(" ElementHook Echo:\n"); _printf_(" numanalyses : "<< this->numanalyses <<"\n"); - + _printf_(" hnodes:\n"); if(hnodes){ @@ -181,9 +181,9 @@ } else _printf_(" hnodes = NULL\n"); - + _printf_(" hvertices:\n"); if(hvertices) hvertices->Echo(); else _printf_(" hvertices = NULL\n"); - + _printf_(" hmaterial:\n"); if(hmaterial) hmaterial->Echo(); Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23066) @@ -286,5 +286,4 @@ IssmDouble calvingratey[NUMVERTICES]; IssmDouble calvingrate[NUMVERTICES]; - /* Get node coordinates and dof list: */ @@ -1048,5 +1047,5 @@ /*}}}*/ void Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/ - + /* Intermediaries */ const int dim=3; @@ -1181,5 +1180,4 @@ if(!IsInputEnum(control_enum)) _error_("Enum "<inputs->GetInput(control_enum); _assert_(input); - /*Cast to Controlinput*/ @@ -2595,5 +2593,5 @@ Penta* penta=this; for(;;){ - + IssmDouble xyz_list[NUMVERTICES][3]; /* Get node coordinates and dof list: */ @@ -2604,5 +2602,5 @@ Jdet[1]=(xyz_list[4][2]-xyz_list[1][2])*0.5; Jdet[2]=(xyz_list[5][2]-xyz_list[2][2])*0.5; - + /*Retrieve all inputs we will need*/ Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); @@ -2615,5 +2613,5 @@ Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - + /* Start looping on the number of 2D vertices: */ for(int ig=0;ig<3;ig++){ @@ -2644,8 +2642,8 @@ delete gauss; } - + /*Stop if we have reached the surface/base*/ if(penta->IsOnSurface()) break; - + /*get upper Penta*/ penta=penta->GetUpperPenta(); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 23066) @@ -104,5 +104,4 @@ return seg; - } /*}}}*/ @@ -141,5 +140,5 @@ /*}}}*/ void Seg::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/ - + /* Intermediaries */ int nrfrontnodes,index; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23066) @@ -106,5 +106,5 @@ } else tria->nodes = NULL; - + tria->vertices = (Vertex**)this->hvertices->deliverp(); tria->material = (Material*)this->hmaterial->delivers(); @@ -115,7 +115,7 @@ /*}}}*/ void Tria::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ - + MARSHALLING_ENUM(TriaEnum); - + /*Call parent classes: */ ElementHook::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); @@ -135,5 +135,5 @@ /*Call inputs method*/ _assert_(this->inputs); - + int domaintype; parameters->FindParam(&domaintype,DomainTypeEnum); @@ -185,5 +185,5 @@ int tria_vertex_ids[3]; - + for(int k=0;k<3;k++){ tria_vertex_ids[k]=reCast(iomodel->elements[3*this->Sid()+k]); //ids for vertices are in the elements array from Matlab @@ -328,5 +328,5 @@ /*}}}*/ void Tria::CalvingCrevasseDepth(){/*{{{*/ - + IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble calvingrate[NUMVERTICES]; @@ -340,8 +340,8 @@ /* Get node coordinates and dof list: */ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - + /*Get the critical fraction of thickness surface and basal crevasses penetrate for calving onset*/ this->parameters->FindParam(&critical_fraction,CalvingCrevasseDepthEnum); - + IssmDouble rho_ice = this->GetMaterialParameter(MaterialsRhoIceEnum); IssmDouble rho_seawater = this->GetMaterialParameter(MaterialsRhoSeawaterEnum); @@ -361,10 +361,10 @@ Input* s_xy_input = inputs->GetInput(DeviatoricStressxyEnum); _assert_(s_xy_input); Input* s_yy_input = inputs->GetInput(DeviatoricStressyyEnum); _assert_(s_yy_input); - + /*Loop over all elements of this partition*/ GaussTria* gauss=new GaussTria(); for (int iv=0;ivGaussVertex(iv); - + H_input->GetInputValue(&thickness,gauss); bed_input->GetInputValue(&bed,gauss); @@ -378,5 +378,5 @@ s_xy_input->GetInputValue(&s_xy,gauss); s_yy_input->GetInputValue(&s_yy,gauss); - + vel=sqrt(vx*vx+vy*vy)+1.e-14; @@ -384,5 +384,5 @@ s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2)); if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;} - + Ho = thickness - (rho_seawater/rho_ice) * (-bed); if(Ho<0.) Ho=0.; @@ -399,5 +399,5 @@ // water_height = surface_crevasse[iv]; //} - + /*basal crevasse*/ //basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (rheology_B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho); @@ -405,11 +405,11 @@ if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.; if (bed>0.) basal_crevasse[iv] = 0.; - + H_surf = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height - critical_fraction*float_depth; H_surfbasal = (surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv])-(critical_fraction*thickness); - + crevasse_depth[iv] = max(H_surf,H_surfbasal); } - + this->inputs->AddInput(new TriaInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1Enum)); this->inputs->AddInput(new TriaInput(BasalCrevasseEnum,&basal_crevasse[0],P1Enum)); @@ -461,5 +461,5 @@ else calvingrate[iv]=0.; - + calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14); calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14); @@ -562,5 +562,5 @@ IssmDouble vorticity_xy[NUMVERTICES]; GaussTria* gauss=NULL; - + /* Get node coordinates and dof list: */ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -569,5 +569,5 @@ Input* vx_input=this->GetInput(EsaXmotionEnum); _assert_(vx_input); Input* vy_input=this->GetInput(EsaYmotionEnum); _assert_(vy_input); - + /* Start looping on the number of vertices: */ gauss=new GaussTria(); @@ -773,5 +773,4 @@ } - }/*}}}*/ void Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/ @@ -1400,5 +1399,5 @@ /*}}}*/ void Tria::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/ - + /* Intermediaries */ int i, dir,nrfrontnodes; @@ -1489,5 +1488,5 @@ }/*}}}*/ void Tria::GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/ - + /* GetLevelsetIntersection computes: * 1. indices of element, sorted in [iceverts, noiceverts] in counterclockwise fashion, @@ -1565,5 +1564,5 @@ /*}}}*/ void Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/ - + /*Computeportion of the element that has a positive levelset*/ @@ -1679,9 +1678,9 @@ parameters->FindParam(&M,NULL,ControlInputSizeMEnum); - + /*Cast to Controlinput*/ if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput"); ControlInput* controlinput = xDynamicCast(input); - + if(strcmp(data,"value")==0){ input = controlinput->values; @@ -1700,5 +1699,5 @@ } /*Check what input we are dealing with*/ - + switch(input->ObjectEnum()){ case TriaInputEnum: @@ -1717,5 +1716,5 @@ break; } - + case TransientInputEnum: { @@ -1982,5 +1981,5 @@ base_input->GetInputAverage(&bed); bed_input->GetInputAverage(&bathymetry); - + /*Return: */ return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.)); @@ -2026,6 +2025,4 @@ if(control_analysis && ad_analysis) iomodel->FindConstant(&num_control_type,"md.autodiff.num_independent_objects"); if(control_analysis && ad_analysis) iomodel->FindConstant(&num_responses,"md.autodiff.num_dependent_objects"); - - /*Recover vertices ids needed to initialize inputs*/ @@ -2394,5 +2391,4 @@ IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/ - /*intermediary: */ IssmDouble* values=NULL; @@ -2407,5 +2403,5 @@ IssmDouble fraction1,fraction2; bool mainlynegative=true; - + /*Output:*/ volume=0; @@ -2416,5 +2412,5 @@ /*Retrieve inputs required:*/ thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); - + /*Retrieve material parameters: */ rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); @@ -2425,5 +2421,5 @@ values[i]=levelset[this->vertices[i]->Sid()]; } - + /*Ok, use the level set values to figure out where we put our gaussian points:*/ this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values); @@ -2877,5 +2873,5 @@ dist_cf_input->GetInputValue(&dist_cf,gauss); delete gauss; - + /*Ensure values are positive for floating ice*/ dist_gl = fabs(dist_gl); @@ -3197,5 +3193,4 @@ } - } /*}}}*/ @@ -3278,5 +3273,4 @@ int idlist[NUMVERTICES],control_init; - /*Get Domain type*/ int domaintype; @@ -3298,5 +3292,5 @@ /*Get out if this is not an element input*/ if(!IsInputEnum(control_enum)) return; - + Input* input = (Input*)this->inputs->GetInput(control_enum); _assert_(input); if(input->ObjectEnum()!=ControlInputEnum){ @@ -3330,5 +3324,4 @@ IssmDouble values[NUMVERTICES]; int vertexpidlist[NUMVERTICES],control_init; - /*Get Domain type*/ @@ -4167,5 +4160,5 @@ IssmDouble* H_elastic_precomputed = NULL; int M, hemi; - + /*computation of Green functions:*/ IssmDouble* U_elastic= NULL; @@ -4174,5 +4167,5 @@ IssmDouble* X_elastic= NULL; IssmDouble* Y_elastic= NULL; - + /*optimization:*/ bool store_green_functions=false; @@ -4182,5 +4175,5 @@ if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!"); deltathickness_input->GetInputAverage(&I); - + /*early return if we are not on the (ice) loading point: */ if(I==0) return; @@ -4195,5 +4188,5 @@ /*which hemisphere? for north-south, east-west components*/ this->parameters->FindParam(&hemi,EsaHemisphereEnum); - + /*compute area of element:*/ area=GetArea(); @@ -4254,5 +4247,5 @@ Y_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*Y_elastic[i]; X_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*X_elastic[i]; - + /*North-south, East-west components */ if (hemi == -1) { @@ -4277,5 +4270,5 @@ pX->SetValues(gsize,indices,X_values,ADD_VAL); pY->SetValues(gsize,indices,Y_values,ADD_VAL); - + /*free ressources:*/ xDelete(indices); @@ -4307,10 +4300,10 @@ IssmDouble* H_elastic_precomputed = NULL; int M; - + /*computation of Green functions:*/ IssmDouble* U_elastic= NULL; IssmDouble* N_elastic= NULL; IssmDouble* E_elastic= NULL; - + /*optimization:*/ bool store_green_functions=false; @@ -4320,5 +4313,5 @@ if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!"); deltathickness_input->GetInputAverage(&I); - + /*early return if we are not on the (ice) loading point: */ if(I==0) return; @@ -4419,5 +4412,5 @@ 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); E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); - + /*Elastic component (from Eq 17 in Adhikari et al, GMD 2015): */ int index=reCast(alpha/PI*(M-1)); @@ -4434,5 +4427,5 @@ pNorth->SetValues(gsize,indices,N_values,ADD_VAL); pEast->SetValues(gsize,indices,E_values,ADD_VAL); - + /*free ressources:*/ xDelete(indices); @@ -4455,5 +4448,5 @@ if(IsWaterInElement()){ - + IssmDouble area; @@ -4528,11 +4521,11 @@ if(IsWaterInElement()){ IssmDouble rho_water, S; - + /*recover material parameters: */ rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum); - + /*From Sg_old, recover water sea level rise:*/ S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; - + /* Perturbation terms for moment of inertia (moi_list): * computed analytically (see Wu & Peltier, eqs 10 & 32) @@ -4546,18 +4539,18 @@ else if(this->inputs->Max(MaskIceLevelsetEnum)<0){ IssmDouble rho_ice, I; - + /*recover material parameters: */ rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); - + /*Compute ice thickness change: */ Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); deltathickness_input->GetInputAverage(&I); - + dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea; dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea; dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea; } - + return; }/*}}}*/ @@ -4592,5 +4585,5 @@ bool computeelastic= true; bool scaleoceanarea= false; - + /*early return if we are not on an ice cap:*/ if(!(this->inputs->Max(MaskIceLevelsetEnum)<=0)){ @@ -4606,10 +4599,10 @@ return; } - + /*If we are here, we are on ice that is fully grounded or half-way to floating: */ if ((this->inputs->Min(MaskGroundediceLevelsetEnum))<0){ notfullygrounded=true; //used later on. } - + /*Inform mask: */ constant=1; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum)); @@ -4636,8 +4629,8 @@ /* Where is the centroid of this element?:{{{*/ - + /*retrieve coordinates: */ ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical); - + IssmDouble minlong=400; IssmDouble maxlong=-20; @@ -4653,10 +4646,10 @@ if (llr_list[2][1]==0)llr_list[2][1]=360; } - + // correction at the north pole if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; - + //correction at the south pole if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; @@ -4684,5 +4677,5 @@ area*=phi; } - + /*Compute ice thickness change: */ Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); @@ -4696,5 +4689,5 @@ int point1; IssmDouble fraction1,fraction2; - + /*Recover portion of element that is grounded*/ this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); @@ -4750,10 +4743,10 @@ } pSgi->SetValues(gsize,indices,values,ADD_VAL); - + /*free ressources:*/ xDelete(values); xDelete(indices); } - + /*Assign output pointer:*/ _assert_(!xIsNan(eustatic)); @@ -4780,5 +4773,5 @@ IssmDouble* G_elastic_precomputed = NULL; int M; - + /*computation of Green functions:*/ IssmDouble* G_elastic= NULL; @@ -4857,7 +4850,7 @@ longe=longe/180*PI; /*}}}*/ - + if(computeelastic){ - + /*recover elastic green function:*/ DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); @@ -4897,5 +4890,5 @@ } } - + pSgo->SetValues(gsize,indices,values,ADD_VAL); @@ -4930,5 +4923,5 @@ IssmDouble* H_elastic_precomputed = NULL; int M; - + /*computation of Green functions:*/ IssmDouble* U_elastic= NULL; @@ -4957,5 +4950,5 @@ this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); - + /*early return if elastic not requested:*/ if(!computeelastic) return; @@ -5026,10 +5019,10 @@ /*From Sg, recover water sea level rise:*/ S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; - + /*Compute ice thickness change: */ Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); deltathickness_input->GetInputAverage(&I); - + /*initialize: */ U_elastic=xNewZeroInit(gsize); @@ -5073,5 +5066,5 @@ E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); } - + /*Elastic component (from Eq 17 in Adhikari et al, GMD 2015): */ int index=reCast(alpha/PI*(M-1)); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23066) @@ -242,5 +242,4 @@ char *lockfilename = NULL; bool waitonlock = false; - /*Write lock file if requested: */ @@ -4393,5 +4392,4 @@ xDelete(RSLg_old); - } /*}}}*/ @@ -4807,5 +4805,4 @@ size_t* poutputbuffersize; - /*Before we delete the profiler, report statistics for this run: */ profiler->Stop(TOTAL); //final tagging @@ -4868,5 +4865,4 @@ }/*}}}*/ #endif - #if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_) @@ -5246,5 +5242,5 @@ newy[i] = newxylist[2*i+1]; } - + /*Assign the pointers*/ (*pnewelementslist) = newelementslist; //Matlab indexing Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 23066) @@ -157,5 +157,5 @@ bool autodiff=false; bool iscontrol=false; - + /*First, keep track of the file handle: */ this->fid=iomodel_handle; @@ -422,5 +422,5 @@ this->FetchData(&autodiff,"md.autodiff.isautodiff"); this->FetchData(&iscontrol,"md.inversion.iscontrol"); - + if(trace || (autodiff && !iscontrol)){ @@ -506,5 +506,5 @@ char** stringarray=*pstringarray; - + if(numstrings){ for(int i=0;ifid)!=1) _error_("could not read integer "); - + /*Convert codes to Enums if needed*/ if(strcmp(record_name,"md.smb.model")==0) integer = IoCodeToEnumSMB(integer); @@ -963,5 +963,5 @@ /*recover my_rank:*/ int my_rank=IssmComm::GetRank(); - + /*Set file pointer to beginning of the data: */ fid=this->SetFilePointerToData(&code,NULL,data_name); @@ -1134,5 +1134,5 @@ } } - + /*output: */ int M,N; @@ -1825,5 +1825,5 @@ /*check we are indeed finding a string, not something else: */ if(codes[i]!=4)_error_("expecting a string for \""<(codes); xDelete(file_positions); - + /*Assign output pointers: */ *pstrings=strings; @@ -1875,5 +1875,5 @@ /*recover my_rank:*/ int my_rank=IssmComm::GetRank(); - + /*Get file pointers to beginning of the data (multiple instances of it): */ file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name); @@ -1890,5 +1890,5 @@ if(code!=2)_error_("expecting an integer for \""<(file_positions); @@ -1928,5 +1928,5 @@ /*recover my_rank:*/ int my_rank=IssmComm::GetRank(); - + /*Get file pointers to beginning of the data (multiple instances of it): */ file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name); @@ -1943,5 +1943,5 @@ if(code!=3)_error_("expecting a double for \""<(file_positions); @@ -1985,5 +1985,5 @@ /*recover my_rank:*/ int my_rank=IssmComm::GetRank(); - + /*Get file pointers to beginning of the data (multiple instances of it): */ file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name); @@ -2014,5 +2014,4 @@ } ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - /*Now allocate matrix: */ @@ -2039,6 +2038,5 @@ else matrix=NULL; - - + /*Assign: */ mdims[i]=M; @@ -2047,5 +2045,5 @@ } } - + /*Free ressources:*/ xDelete(file_positions); @@ -2089,5 +2087,5 @@ /*recover my_rank:*/ int my_rank=IssmComm::GetRank(); - + /*Get file pointers to beginning of the data (multiple instances of it): */ file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name); @@ -2118,5 +2116,4 @@ } ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - /*Now allocate matrix: */ @@ -2144,6 +2141,5 @@ else integer_matrix=NULL; - - + /*Assign: */ mdims[i]=M; @@ -2152,5 +2148,5 @@ } } - + /*Free ressources:*/ xDelete(file_positions); @@ -2376,5 +2372,5 @@ vector_types = xNew(num_instances); } - + /*Reset FILE* position to the beginning of the file, and start again, this time saving the data information * as we find it: */ @@ -2419,10 +2415,10 @@ vector_types[counter] = vector_type; fgetpos(fid,file_positions+counter); - + /*backup and skip over the record, as we have more work to do: */ if(5<=record_code && record_code<=7) fseek(fid,-sizeof(int),SEEK_CUR); fseek(fid,-sizeof(int),SEEK_CUR); fseek(fid,-sizeof(int),SEEK_CUR); - + /*increment counter: */ counter++; Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23066) @@ -196,5 +196,5 @@ if(vmag==0. && (s-1.)<0.) alpha_complement=0.; else alpha_complement=pow(Neff,r)*pow(vmag,(s-1)); - + _assert_(!xIsNan(alpha_complement)); _assert_(!xIsInf(alpha_complement)); @@ -270,5 +270,4 @@ r=drag_q/drag_p; s=1./drag_p; - /*Get effective pressure*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23066) @@ -178,5 +178,5 @@ switch(analysis_type){ - + case HydrologyShaktiAnalysisEnum: pe = this->CreatePVectorHydrologyShakti(); @@ -395,8 +395,8 @@ IssmDouble moulin_load,dt; ElementVector* pe=new ElementVector(&node,1,this->parameters); - + this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum); parameters->FindParam(&dt,TimesteppingTimeStepEnum); - + pe->values[0]=moulin_load*dt; /*Clean up and return*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp (revision 23066) @@ -30,5 +30,4 @@ /*}}}*/ Neumannflux::Neumannflux(int neumannflux_id,int i,IoModel* iomodel,int* segments,int in_analysis_type){/*{{{*/ - /*Some sanity checks*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 23066) @@ -109,5 +109,5 @@ if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook(); - + MARSHALLING_ENUM(MatestarEnum); MARSHALLING(mid); Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 23066) @@ -135,5 +135,5 @@ //if(helement) helement->DeepEcho(); //else _printf_(" helement = NULL\n"); - + _printf_(" element:\n"); _printf_(" note: element not printed to avoid recursion.\n"); @@ -143,10 +143,10 @@ /*}}}*/ void Matice::Echo(void){/*{{{*/ - + _printf_("Matice:\n"); _printf_(" mid: " << mid << "\n"); _printf_(" isdamaged: " << isdamaged << "\n"); _printf_(" isenhanced: " << isenhanced << "\n"); - + /*helement and element Echo were commented to avoid recursion.*/ /*Example: element->Echo calls matice->Echo which calls element->Echo etc*/ @@ -155,5 +155,5 @@ //if(helement) helement->Echo(); //else _printf_(" helement = NULL\n"); - + _printf_(" element:\n"); _printf_(" note: element not printed to avoid recursion.\n"); @@ -167,5 +167,5 @@ if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook(); - + MARSHALLING_ENUM(MaticeEnum); MARSHALLING(mid); @@ -541,5 +541,4 @@ IssmDouble exx,eyy,exy,exz,eyz; - if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && (epsilon[3]==0) && (epsilon[4]==0)){ Index: /issm/trunk-jpl/src/c/classes/Materials/Matlitho.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matlitho.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Materials/Matlitho.cpp (revision 23066) @@ -37,11 +37,11 @@ this->radius=xNew(this->numlayers+1); xMemCpy(this->radius, iomodel->Data("md.materials.radius"),this->numlayers+1); - + this->viscosity=xNew(this->numlayers); xMemCpy(this->viscosity, iomodel->Data("md.materials.viscosity"),this->numlayers); - + this->lame_lambda=xNew(this->numlayers); xMemCpy(this->lame_lambda, iomodel->Data("md.materials.lame_lambda"),this->numlayers); - + this->lame_mu=xNew(this->numlayers); xMemCpy(this->lame_mu, iomodel->Data("md.materials.lame_mu"),this->numlayers); @@ -61,15 +61,15 @@ this->issolid=xNew(this->numlayers); xMemCpy(this->issolid, iomodel->Data("md.materials.issolid"),this->numlayers); - + /*isburgersd= xNew(this->numlayers); this->isburgers=xNew(this->numlayers); xMemCpy(isburgersd, iomodel->Data("md.materials.isburgers"),this->numlayers); for (int i=0;inumlayers;i++)this->isburgers[i]=reCast(isburgersd[i]); - + issolidd= xNew(this->numlayers); this->issolid=xNew(this->numlayers); xMemCpy(issolidd, iomodel->Data("md.materials.issolid"),this->numlayers); for (int i=0;inumlayers;i++)this->issolid[i]=reCast(issolidd[i]);*/ - + /*free ressources: */ xDelete(isburgersd); @@ -78,5 +78,5 @@ /*}}}*/ Matlitho::~Matlitho(){/*{{{*/ - + xDelete(radius); xDelete(viscosity); Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 23066) @@ -59,5 +59,5 @@ mantle_shear_modulus = 0; mantle_density = 0; - + earth_density = 0; @@ -338,5 +338,5 @@ matpar->mantle_shear_modulus=this->mantle_shear_modulus; matpar->mantle_density=this->mantle_density; - + matpar->earth_density=this->earth_density; @@ -439,5 +439,5 @@ MARSHALLING(mantle_shear_modulus); MARSHALLING(mantle_density); - + //slr: MARSHALLING(earth_density); Index: /issm/trunk-jpl/src/c/classes/Misfit.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Misfit.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Misfit.cpp (revision 23066) @@ -23,5 +23,5 @@ #include "../classes/gauss/Gauss.h" /*}}}*/ - + /*Misfit constructors, destructors :*/ Misfit::Misfit(){/*{{{*/ @@ -42,5 +42,5 @@ this->definitionenum=in_definitionenum; - + this->name = xNew(strlen(in_name)+1); xMemCpy(this->name,in_name,strlen(in_name)+1); @@ -48,10 +48,10 @@ this->timeinterpolation = xNew(strlen(in_timeinterpolation)+1); xMemCpy(this->timeinterpolation,in_timeinterpolation,strlen(in_timeinterpolation)+1); - + this->model_enum=in_model_enum; this->observation_enum=in_observation_enum; this->weights_enum=in_weights_enum; this->local=in_local; - + this->misfit=0; this->lock=0; @@ -111,9 +111,9 @@ /*}}}*/ IssmDouble Misfit::Response(FemModel* femmodel){/*{{{*/ - + /*diverse: */ IssmDouble time,starttime,finaltime; IssmDouble dt; - + /*recover time parameters: */ femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); @@ -130,5 +130,4 @@ IssmDouble all_area_t; - /*If we are locked, return time average: */ if(this->lock)return misfit/(time-starttime); @@ -144,5 +143,5 @@ area_t=all_area_t; misfit_t=all_misfit_t; - + /*Divide by surface area if not nill!: */ if (area_t!=0) misfit_t=misfit_t/area_t; @@ -164,9 +163,9 @@ IssmDouble* weights= NULL; int msize,osize,wsize; - + /*Are we transient?:*/ if (time==0){ IssmDouble misfit_t=0.; - + /*get global vectors: */ GetVectorFromInputsx(&model,&msize,femmodel,model_enum); @@ -190,5 +189,5 @@ } else{ - + IssmDouble misfit_t=0.; IssmDouble all_misfit_t=0.; @@ -201,5 +200,5 @@ GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize); GetVectorFromInputsx(&weights,&wsize,femmodel,weights_enum); _assert_(wsize==msize); - + int count=0; for (int i=0;ilock=1; - + /*Free ressources:*/ xDelete(model); @@ -227,7 +226,7 @@ } /*}}}*/ else{ /*global computation: {{{ */ - + IssmDouble model, observation; - + /*If we are locked, return time average: */ if(this->lock)return misfit/(time-starttime); @@ -239,11 +238,11 @@ Input* input = element->GetInput(observation_enum); _assert_(input); input->GetInputAverage(&observation); - + /*Add this time's contribution to curent misfit: */ misfit+=dt*(model-observation); - + /*Do we lock? i.e. are we at final_time? :*/ if(time==finaltime)this->lock=1; - + /*What we return is the value of misfit / time: */ return misfit/(time-starttime); Index: /issm/trunk-jpl/src/c/classes/Nodalvalue.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Nodalvalue.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Nodalvalue.cpp (revision 23066) @@ -87,5 +87,5 @@ /*}}}*/ IssmDouble Nodalvalue::Response(FemModel* femmodel){/*{{{*/ - + /*output:*/ IssmDouble value; Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 23066) @@ -498,5 +498,5 @@ /*}}}*/ void Node::SetApproximation(int in_approximation){/*{{{*/ - + this->approximation = in_approximation; } Index: /issm/trunk-jpl/src/c/classes/Numberedcostfunction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Numberedcostfunction.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Numberedcostfunction.cpp (revision 23066) @@ -8,5 +8,4 @@ #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" #endif - /*Headers:*/ @@ -30,5 +29,5 @@ /*}}}*/ - + /*Numberedcostfunction constructors, destructors :*/ Numberedcostfunction::Numberedcostfunction(){/*{{{*/ @@ -110,5 +109,5 @@ /*}}}*/ IssmDouble Numberedcostfunction::Response(FemModel* femmodel){/*{{{*/ - + _assert_(number_cost_functions>0 && number_cost_functions<1e3); @@ -157,3 +156,2 @@ } /*}}}*/ - Index: /issm/trunk-jpl/src/c/classes/Params/BoolParam.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/BoolParam.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Params/BoolParam.cpp (revision 23066) @@ -63,3 +63,2 @@ } /*}}}*/ - Index: /issm/trunk-jpl/src/c/classes/Params/DataSetParam.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DataSetParam.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Params/DataSetParam.cpp (revision 23066) @@ -55,5 +55,5 @@ if(marshall_direction==MARSHALLING_BACKWARD)value=new DataSet(); - + MARSHALLING_ENUM(DataSetParamEnum); MARSHALLING(enum_type); Index: /issm/trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.cpp (revision 23066) @@ -233,3 +233,2 @@ } /*}}}*/ - Index: /issm/trunk-jpl/src/c/classes/Params/DoubleMatParam.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DoubleMatParam.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleMatParam.cpp (revision 23066) @@ -116,5 +116,5 @@ /*DoubleMatParam specific routines:*/ void DoubleMatParam::GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM,int* pN){/*{{{*/ - + /*Assign output pointers:*/ if(pM) *pM=M; Index: /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp (revision 23066) @@ -116,5 +116,5 @@ /*DoubleVecParam specific routines:*/ void DoubleVecParam::GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM){/*{{{*/ - + /*Assign output pointers:*/ if(pM) *pM=M; Index: /issm/trunk-jpl/src/c/classes/Params/IntParam.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/IntParam.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Params/IntParam.cpp (revision 23066) @@ -64,3 +64,2 @@ } /*}}}*/ - Index: /issm/trunk-jpl/src/c/classes/Params/StringArrayParam.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/StringArrayParam.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Params/StringArrayParam.cpp (revision 23066) @@ -83,5 +83,5 @@ MARSHALLING_ENUM(StringArrayParamEnum); - + MARSHALLING(enum_type); MARSHALLING(numstrings); Index: /issm/trunk-jpl/src/c/classes/Params/StringParam.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/StringParam.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Params/StringParam.cpp (revision 23066) @@ -55,5 +55,5 @@ if(marshall_direction==MARSHALLING_FORWARD || marshall_direction == MARSHALLING_SIZE)size=strlen(value)+1; - + MARSHALLING_ENUM(StringParamEnum); MARSHALLING(enum_type); Index: /issm/trunk-jpl/src/c/classes/Profiler.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Profiler.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/Profiler.cpp (revision 23066) @@ -100,5 +100,4 @@ if(this->running[tag]) _error_("Tag "<time[tag]; @@ -143,5 +142,4 @@ _assert_(tagrunning[tag]) _error_("Tag "<running[tag]) _error_("Tag "<numgauss = order; GaussLegendreLinear(&pcoords1,&pweights,order); - + this->coords1=xNew(numgauss); this->weights=xNew(numgauss); Index: /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp (revision 23066) @@ -741,3 +741,2 @@ xDelete(counter); }/*}}}*/ - Index: /issm/trunk-jpl/src/c/cores/ad_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/ad_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/ad_core.cpp (revision 23066) @@ -46,5 +46,5 @@ /*First, stop tracing: */ trace_off(); - + /*Print tape statistics so that user can kill this run if something is off already:*/ if(VerboseAutodiff()){ /*{{{*/ @@ -93,5 +93,5 @@ femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum); - + /*if no dependents, no point in running a driver: */ if(!(num_dependents*num_independents)) return; @@ -101,5 +101,5 @@ num_dependents=0; num_independents=0; } - + /*retrieve state variable: */ femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum); @@ -117,5 +117,4 @@ femmodel->parameters->FindParam(&driver,AutodiffDriverEnum); - if (strcmp(driver,"fos_forward")==0){ /*{{{*/ @@ -140,5 +139,4 @@ anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx; #endif - /*call driver: */ @@ -184,5 +182,4 @@ #endif // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; - /*seed matrix: */ @@ -244,5 +241,4 @@ #endif - /*call driver: */ fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac ); @@ -284,5 +280,4 @@ anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; #endif - /*seed matrix: */ @@ -317,5 +312,4 @@ else _error_("driver: " << driver << " not yet supported!"); - if(VerboseAutodiff())_printf0_(" end AD core\n"); Index: /issm/trunk-jpl/src/c/cores/adgradient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/adgradient_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/adgradient_core.cpp (revision 23066) @@ -39,5 +39,5 @@ IssmPDouble *aWeightVector=NULL; IssmPDouble *weightVectorTimesJac=NULL; - + /*output: */ IssmPDouble *totalgradient=NULL; @@ -53,5 +53,5 @@ /*First, stop tracing: */ trace_off(); - + if(VerboseAutodiff()){ /*{{{*/ size_t tape_stats[15]; @@ -96,9 +96,9 @@ delete [] sstats; } /*}}}*/ - + /*retrieve parameters: */ femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum); - + /*if no dependents, no point in running a driver: */ if(!(num_dependents*num_independents)) return; @@ -109,8 +109,8 @@ num_dependents=0; num_independents=0; } - + /*retrieve state variable: */ femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum); - + /* driver argument */ xp=xNew(num_independents); @@ -131,5 +131,5 @@ /*Initialize outputs: */ totalgradient=xNewZeroInit(num_independents); - + for(aDepIndex=0;aDepIndex(num_dependents); if (my_rank==0) aWeightVector[aDepIndex]=1.0; - + /*initialize output gradient: */ weightVectorTimesJac=xNew(num_independents); @@ -162,10 +162,10 @@ xDelete(aWeightVector); } - + /*add totalgradient to results*/ femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,totalgradient,num_independents,1,1,0.0)); if(VerboseAutodiff())_printf0_(" end ad core\n"); - + /* delete the allocated space for the parameters and free ressources:{{{*/ xDelete(anEDF_for_solverx_p->dp_x); Index: /issm/trunk-jpl/src/c/cores/bmb_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/bmb_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/bmb_core.cpp (revision 23066) @@ -11,5 +11,4 @@ void bmb_core(FemModel* femmodel){ - /*First, get BMB model from parameters*/ Index: /issm/trunk-jpl/src/c/cores/controlad_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controlad_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/controlad_core.cpp (revision 23066) @@ -105,5 +105,5 @@ default: _printf0_(" Unknown end condition\n"); } - + /*Save results:*/ femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,G,n,1,0,0.0)); @@ -130,5 +130,5 @@ IssmPDouble* X=NULL; int i; - + /*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 *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 @@ -143,5 +143,4 @@ femmodel=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, restartfilename,IssmComm::GetComm(), femmodel->solution_type,NULL); - /*Get initial guess:*/ femmodel->parameters->FindParam(&Xd,&intn,AutodiffXpEnum); @@ -175,5 +174,5 @@ IssmDouble pfd; int i; - + /*Recover Femmodel*/ femmodel = (FemModel*)dzs; @@ -195,5 +194,5 @@ femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast(pfd); _printf0_("f(x) = "<(Jlist); xDelete(G2); - + xDelete(rootpath); xDelete(inputfilename); @@ -251,5 +250,5 @@ IssmDouble pfd; int i; - + /*Recover Femmodel*/ femmodel = (FemModel*)dzs; @@ -271,5 +270,5 @@ femmodelad=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, restartfilename,IssmComm::GetComm(), femmodel->solution_type,X); femmodel=femmodelad; //We can do this, because femmodel is being called from outside, not by reference, so we won't erase it - + /*Recover some parameters*/ femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); @@ -284,5 +283,5 @@ femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast(pfd); _printf0_("f(x) = "<(Jlist); xDelete(G2); - xDelete(rootpath); xDelete(inputfilename); Index: /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 23066) @@ -126,5 +126,5 @@ SetControlInputsFromVectorx(femmodel,aX); xDelete(aX); - + /*Compute solution (forward)*/ void (*solutioncore)(FemModel*)=NULL; @@ -392,5 +392,5 @@ long niter = long(maxsteps); /*Maximum number of iterations*/ long nsim = long(maxiter);/*Maximum number of function calls*/ - + /*Get initial guess*/ Vector *Xpetsc = NULL; @@ -401,5 +401,5 @@ delete Xpetsc; //_assert_(intn==numberofvertices*num_controls); - + /*Get problem dimension and initialize gradient and initial guess*/ long n = long(intn); @@ -415,5 +415,5 @@ N_add+=N[c]; } - + /*Allocate m1qn3 working arrays (see documentation)*/ long m = 100; @@ -471,5 +471,5 @@ N_add +=N[c]; } - + /*Set X as our new control*/ IssmDouble* aX=xNew(intn); @@ -483,7 +483,6 @@ ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG); SetControlInputsFromVectorx(femmodel,aX); - + xDelete(aX); - if (solution_type == TransientSolutionEnum){ @@ -506,5 +505,5 @@ femmodel->results->AddObject(G_output); femmodel->results->AddObject(X_output); - + offset += N[i]*numberofvertices; } Index: /issm/trunk-jpl/src/c/cores/dakota_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/dakota_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/dakota_core.cpp (revision 23066) @@ -235,5 +235,4 @@ void dakota_core(FemModel* femmodel){ /*{{{*/ - int my_rank; char *dakota_input_file = NULL; Index: /issm/trunk-jpl/src/c/cores/damage_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/damage_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/damage_core.cpp (revision 23066) @@ -19,5 +19,4 @@ char **requested_outputs = NULL; - //first recover parameters common to all solutions femmodel->parameters->FindParam(&save_results,SaveResultsEnum); Index: /issm/trunk-jpl/src/c/cores/esa_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/esa_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/esa_core.cpp (revision 23066) @@ -23,5 +23,5 @@ int numoutputs = 0; char **requested_outputs = NULL; - + /*additional parameters: */ int gsize; @@ -80,5 +80,5 @@ U_x = new Vector(gsize); U_y = new Vector(gsize); - + /*call the geodetic main modlule:*/ if(domaintype==Domain3DsurfaceEnum){ @@ -95,5 +95,5 @@ InputUpdateFromVectorx(femmodel,U_north,EsaNmotionEnum,VertexSIdEnum); // north motion InputUpdateFromVectorx(femmodel,U_east,EsaEmotionEnum,VertexSIdEnum); // east motion - + if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); @@ -101,5 +101,5 @@ femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); } - + if(solution_type==EsaSolutionEnum)femmodel->RequestedDependentsx(); @@ -115,3 +115,2 @@ } /*}}}*/ - Index: /issm/trunk-jpl/src/c/cores/gia_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/gia_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/gia_core.cpp (revision 23066) @@ -53,5 +53,5 @@ femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2); } - + xDelete(x); xDelete(y); Index: /issm/trunk-jpl/src/c/cores/love_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/love_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/love_core.cpp (revision 23066) @@ -28,10 +28,10 @@ /*parameters: */ bool save_results; - + if(VerboseSolution()) _printf0_(" computing LOVE numbers\n"); /*Recover some parameters: */ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - + /*recover love number parameters: */ femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum); @@ -45,5 +45,5 @@ femmodel->parameters->FindParam(&love_kernels,LoveKernelsEnum); femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum); - + /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/ Matlitho* matlitho=NULL; Index: /issm/trunk-jpl/src/c/cores/masstransport_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 23066) @@ -31,5 +31,5 @@ femmodel->parameters->FindParam(&stabilization,MasstransportStabilizationEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum); - + if(VerboseSolution()) _printf0_(" computing mass transport\n"); Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 23066) @@ -13,15 +13,14 @@ void sealevelrise_core(FemModel* femmodel){ /*{{{*/ - /*Parameters, variables:*/ bool save_results; bool isslr=0; int solution_type; - + /*Retrieve parameters:*/ femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - + /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/ if(solution_type==SealevelriseSolutionEnum)isslr=1; @@ -41,5 +40,5 @@ /*Run steric core for sure:*/ steric_core(femmodel); - + /*Save results: */ if(save_results){ @@ -58,5 +57,4 @@ /*}}}*/ void geodetic_core(FemModel* femmodel){ /*{{{*/ - /*variables:*/ @@ -88,5 +86,5 @@ /*Should we even be here?:*/ femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); if(!geodetic)return; - + /*Verbose: */ if(VerboseSolution()) _printf0_(" computing geodetic sea level rise\n"); @@ -155,5 +153,5 @@ /*recover N_esa = U_esa + RSLg:*/ N_esa=U_esa->Duplicate(); U_esa->Copy(N_esa); N_esa->AXPY(RSLg,1); - + /*transform these values into rates (as we only run this once each frequency turn:*/ N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency)); @@ -162,5 +160,5 @@ U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency)); RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency)); - + /*get some results into elements:{{{*/ InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum); @@ -182,5 +180,4 @@ } - if(iscoupler){ /*transfer sea level back to ice caps:*/ @@ -228,5 +225,5 @@ Vector *U_gia_rate= NULL; Vector *N_gia_rate= NULL; - + /*parameters: */ bool isslr=0; @@ -243,5 +240,5 @@ /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/ if(solution_type==SealevelriseSolutionEnum)isslr=1; - + /*Should we be here?:*/ if(!isslr)return; @@ -291,5 +288,5 @@ /*}}}*/ Vector* sealevelrise_core_eustatic(FemModel* femmodel){ /*{{{*/ - + /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/ @@ -318,5 +315,5 @@ /*Figure out size of g-set deflection vector and allocate solution vector: */ gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); - + /*Initialize:*/ RSLgi = new Vector(gsize); @@ -331,8 +328,7 @@ * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/ RSLgi->Shift(-eustatic-RSLgi_oceanaverage); - + /*save eustatic value for results: */ femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelRSLEustaticEnum,-eustatic)); - /*clean up and return:*/ @@ -346,5 +342,4 @@ /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5. non eustatic core of the SLR solution */ - Vector *RSLg = NULL; @@ -380,5 +375,5 @@ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum); - + /*computational flag: */ femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum); @@ -389,5 +384,5 @@ /*Figure out size of g-set deflection vector and allocate solution vector: */ gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); - + /*Initialize:*/ RSLg = new Vector(gsize); @@ -400,5 +395,5 @@ count=1; converged=false; - + /*Start loop: */ for(;;){ @@ -413,5 +408,5 @@ /*call the non eustatic module: */ femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, latitude, longitude, radius,verboseconvolution,loop); - + /*assemble solution vector: */ RSLgo->Assemble(); @@ -423,5 +418,5 @@ femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz,latitude,longitude,radius); RSLgo_rot->Assemble(); - + /*save changes in inertia tensor as results: */ femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz)); @@ -431,8 +426,8 @@ RSLgo->AXPY(RSLgo_rot,1); } - + /*we need to average RSLgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */ RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo); - + /*RSLg is the sum of the eustatic term, and the ocean terms: */ RSLg_eustatic->Copy(RSLg); RSLg->AXPY(RSLgo,1); @@ -455,8 +450,8 @@ break; } - + /*some minor verbosing adjustment:*/ if(count>1)verboseconvolution=false; - + } if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n"); @@ -474,5 +469,5 @@ Vector *U_north_esa = NULL; Vector *U_east_esa = NULL; - + /*parameters: */ int configuration_type; @@ -507,5 +502,5 @@ VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); - + /*call the elastic main modlule:*/ femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,latitude,longitude,radius,xx,yy,zz,loop,horiz); @@ -532,5 +527,5 @@ Vector *U_gia = NULL; Vector *N_gia = NULL; - + /*parameters:*/ int frequency; @@ -570,5 +565,5 @@ Vector* forcingglobal=NULL; int* nvs=NULL; - + /*transition vectors:*/ IssmDouble** transitions=NULL; @@ -577,5 +572,5 @@ int* transitions_n=NULL; int nv; - + /*communicators:*/ ISSM_MPI_Comm tocomm; @@ -591,5 +586,5 @@ femmodel->parameters->FindParam(&nummodels,NumModelsEnum); my_rank=IssmComm::GetRank(); - + /*retrieve the inter communicators that will be used to send data from each ice cap to the earth: */ if(modelid==earthid){ @@ -620,5 +615,5 @@ ISSM_MPI_Recv(forcings[i], nvs[i], ISSM_MPI_DOUBLE, 0,i, fromcomms[i], &status); } - + } else{ @@ -631,5 +626,5 @@ /*On the earth model, consolidate all the forcings into one, and update the elements dataset accordingly: {{{*/ if(modelid==earthid){ - + /*Out of all the delta thicknesses, build one delta thickness vector made of all the ice cap contributions. *First, build the global delta thickness vector in the earth model: */ @@ -652,5 +647,4 @@ /*build index to plug values: */ int* index=xNew(M); for(int i=0;i(transition[i])-1; //matlab indexing! - /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular @@ -663,5 +657,5 @@ /*Assemble vector:*/ forcingglobal->Assemble(); - + /*Plug into elements:*/ InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum); @@ -696,5 +690,5 @@ IssmDouble* forcing=NULL; IssmDouble* forcingglobal=NULL; - + /*transition vectors:*/ IssmDouble** transitions=NULL; @@ -703,5 +697,5 @@ int* transitions_n=NULL; int nv; - + /*communicators:*/ ISSM_MPI_Comm fromcomm; @@ -718,5 +712,5 @@ femmodel->parameters->FindParam(&nummodels,NumModelsEnum); my_rank=IssmComm::GetRank(); - + /*retrieve the inter communicators that will be used to send data from earth to ice caps:*/ if(modelid==earthid){ @@ -733,5 +727,4 @@ } - /*Retrieve sea-level on earth model: */ if(modelid==earthid){ @@ -742,10 +735,10 @@ /*Send the forcing to the ice caps:{{{*/ if(my_rank==0){ - + if(modelid==earthid){ - + /*Retrieve transition vectors, used to figure out global forcing contribution to each ice cap's own elements: */ femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum); - + if(ntransitions!=earthid)_error_("TransferSeaLevel error message: number of transition vectors is not equal to the number of icecaps!"); @@ -804,5 +797,5 @@ Vector *cumdeltathickness = NULL; int nv; - + if(VerboseSolution()) _printf0_(" computing earth mass transport\n"); @@ -840,5 +833,5 @@ } /*}}}*/ void slrconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ - + bool converged=true; IssmDouble ndS,nS; @@ -848,12 +841,11 @@ dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0); ndS=dRSLg->Norm(NORM_TWO); - + if (xIsNan(ndS)) _error_("convergence criterion is NaN!"); - + if(!xIsNan(eps_rel)){ nS=RSLg_old->Norm(NORM_TWO); if (xIsNan(nS)) _error_("convergence criterion is NaN!"); } - //clean up Index: /issm/trunk-jpl/src/c/cores/smb_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/smb_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/smb_core.cpp (revision 23066) @@ -29,7 +29,7 @@ femmodel->parameters->FindParam(&numoutputs,SmbNumRequestedOutputsEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SmbRequestedOutputsEnum); - + if(VerboseSolution()) _printf0_(" computing smb \n"); - + analysis = new SmbAnalysis(); analysis->Core(femmodel); Index: /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 23066) @@ -22,5 +22,4 @@ char **requested_outputs = NULL; Analysis *analysis = NULL; - /* recover parameters:*/ @@ -36,5 +35,5 @@ femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,StressbalanceRequestedOutputsEnum); - + if(VerboseSolution()) _printf0_(" computing new velocity\n"); @@ -79,10 +78,9 @@ } - if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); } - + if(solution_type==StressbalanceSolutionEnum && !control_analysis)femmodel->RequestedDependentsx(); Index: /issm/trunk-jpl/src/c/datastructures/DataSet.cpp =================================================================== --- /issm/trunk-jpl/src/c/datastructures/DataSet.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/datastructures/DataSet.cpp (revision 23066) @@ -92,5 +92,5 @@ /*Specific methods*/ void DataSet::Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ - + vector::iterator obj; int obj_size=0; Index: /issm/trunk-jpl/src/c/main/esmfbinders.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/esmfbinders.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/main/esmfbinders.cpp (revision 23066) @@ -2,5 +2,4 @@ * \brief: ESMF binders for ISSM. Binders developed initially for the GEOS-5 framework. */ - #include "./issm.h" @@ -36,5 +35,5 @@ /*Some specific code here for the binding: */ femmodel->parameters->SetParam(SMBgcmEnum,SmbEnum); //bypass SMB model, will be provided by GCM! - + /*Restart file: */ femmodel->Restart(); @@ -113,9 +112,9 @@ IssmDouble surface; - + /*Recover surface from the ISSM element: */ Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); surface_input->GetInputAverage(&surface); - + *(issmoutputs+f*numberofelements+i) = surface; @@ -147,5 +146,5 @@ /*Output results: */ OutputResultsx(femmodel); - + /*Check point: */ femmodel->CheckPoint(); Index: /issm/trunk-jpl/src/c/main/issm_dakota.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/issm_dakota.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/main/issm_dakota.cpp (revision 23066) @@ -16,5 +16,4 @@ int main(int argc,char **argv){ - #if defined(_HAVE_DAKOTA_) && _DAKOTA_MAJOR_ >= 6 @@ -29,9 +28,9 @@ /*Initialize MPI: */ ISSM_MPI_Init(&argc, &argv); // initialize MPI - + /*Recover file name for dakota input file:*/ dakota_input_file=xNew((strlen(argv[2])+strlen(argv[3])+strlen(".qmu.in")+2)); sprintf(dakota_input_file,"%s/%s%s",argv[2],argv[3],".qmu.in"); - + dakota_output_file=xNew((strlen(argv[2])+strlen(argv[3])+strlen(".qmu.out")+2)); sprintf(dakota_output_file,"%s/%s%s",argv[2],argv[3],".qmu.out"); @@ -57,5 +56,5 @@ Dakota::abort_handler(-1); } - + Dakota::ProblemDescDB& problem_db = env.problem_description_db(); Dakota::ModelLIter ml_iter; Index: /issm/trunk-jpl/src/c/main/issm_ocean.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/issm_ocean.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/main/issm_ocean.cpp (revision 23066) @@ -22,5 +22,5 @@ /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/ worldcomm=EnvironmentInit(argc,argv); - + /*What is my rank?:*/ ISSM_MPI_Comm_rank(worldcomm,&my_rank); @@ -39,5 +39,5 @@ FemModel *femmodel = new FemModel(argc,argv,modelcomm); - + /*Now that the models are initialized, keep communicator information in the parameters datasets of each model: */ femmodel->parameters->AddObject(new GenericParam(worldcomm,WorldCommEnum)); Index: /issm/trunk-jpl/src/c/main/issm_slr.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/issm_slr.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/main/issm_slr.cpp (revision 23066) @@ -28,5 +28,5 @@ /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/ worldcomm=EnvironmentInit(argc,argv); - + /*What is my rank?:*/ ISSM_MPI_Comm_rank(worldcomm,&my_rank); @@ -40,9 +40,9 @@ for(int i=0;i(strlen(argv[5+3*i])+1); xMemCpy(string,argv[5+3*i],strlen(argv[5+3*i])+1); dirnames[i]=string; - + string=xNew(strlen(argv[5+3*i+1])+1); xMemCpy(string,argv[5+3*i+1],strlen(argv[5+3*i+1])+1); @@ -94,5 +94,5 @@ /*Initialize femmodel from arguments provided command line: */ FemModel *femmodel = new FemModel(4,arguments,modelcomm); - + /*Now that the models are initialized, keep communicator information in the parameters datasets of each model: */ femmodel->parameters->AddObject(new GenericParam(worldcomm,WorldCommEnum)); Index: /issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp (revision 23066) @@ -78,5 +78,5 @@ if(ppf) pf =new Vector(flocalsize,fsize); } - + /*Free ressources: */ xDelete(toolkittype); Index: /issm/trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp (revision 23066) @@ -39,5 +39,5 @@ IssmDouble d0=1.e+10; IssmDouble xi,yi; - + //recover vertex position: xi=x[i]; yi=y[i]; Index: /issm/trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp (revision 23066) @@ -52,5 +52,5 @@ double x2,y2; double mind; - + for(int i=i0;iparameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum); femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); - + IssmDouble* boxareas=xNewZeroInit(num_basins*maxbox); @@ -122,5 +122,5 @@ ISSM_MPI_Allreduce(boxareas,sumareas,num_basins*maxbox,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); //if(sumareas[0]==0){_error_("No elements in box 0, basal meltrates will be 0. Consider decreasing md.basalforcings.maxboxcount or refining your mesh!");} - + /*Update parameters to keep track of the new areas in future calculations*/ femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoBoxAreaEnum,sumareas,maxbox*num_basins)); Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 23066) @@ -67,3 +67,2 @@ } /*}}}*/ - Index: /issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp (revision 23066) @@ -27,5 +27,4 @@ int size = 0; for(int i=0;iFindParam(&num_controls,InversionNumControlParametersEnum); parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); - /*2. Allocate vector*/ Index: /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 23066) @@ -84,5 +84,5 @@ /*this one is special: we don't specify the type, but let the nature of the inputs dictace. * P0 -> ElementSIdEnum, P1 ->VertexSIdEnum: */ - + /*We go find the input of the first element, and query its interpolation type: */ Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(0)); Index: /issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp (revision 23066) @@ -15,5 +15,4 @@ Element *element = NULL; - /*retrieve parameters: */ parameters->FindParam(&migration_style,GroundinglineMigrationEnum); @@ -21,5 +20,5 @@ if(migration_style==NoneEnum) return; - + if(VerboseModule()) _printf0_(" Migrating grounding line\n"); Index: /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp (revision 23066) @@ -166,8 +166,8 @@ /*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): */ - + /*if (el==318 && contouri==9){ _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]); - + _printf_("Bool" << (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum)); }*/ @@ -199,5 +199,4 @@ segments_dataset->AddObject(new Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1])); - /*This case is impossible: not quite! */ //_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]); @@ -226,5 +225,5 @@ } else if( (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum) ){ - + /*if (el==318 && contouri==9){ _printf_("hello" << " NodeInElement 0 " << (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) << " NodeInElement 1 " << (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1]))); @@ -251,5 +250,5 @@ IsIdenticalNode(xnodes[1],ynodes[1],xsegment[0],ysegment[0],tolerance) || IsIdenticalNode(xnodes[2],ynodes[2],xsegment[0],ysegment[0],tolerance)){ - + /*ok, segments[0] is common to one of our vertices: */ coord1=0; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 23066) @@ -69,9 +69,9 @@ int num_independent_objects,M; char** names = NULL; - + /*this is done somewhere else*/ parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_independent_objects",InversionNumControlParametersEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_dependent_objects",InversionNumCostFunctionsEnum)); - + /*Step1: create controls (independents)*/ iomodel->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects"); @@ -96,8 +96,8 @@ xDelete(cm_responses); parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,costfunc_enums,num_costfunc)); - + iomodel->FetchData(&control_scaling_factors,NULL,NULL,"md.autodiff.independent_scaling_factors"); parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_independent_objects)); - + /*cleanup*/ for(int i=0;iFindConstant(&control_analysis,"md.inversion.iscontrol"); @@ -29,5 +27,5 @@ /*Now, return if no control*/ if(!control_analysis) return; - + /*Process controls and convert from string to enums*/ iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters"); @@ -48,5 +46,5 @@ iomodel->FetchData(3,"md.inversion.cost_functions_coefficients","md.inversion.min_parameters","md.inversion.max_parameters"); - + /*Fetch Observations */ iomodel->FindConstant(&domaintype,"md.mesh.domain_type"); @@ -154,5 +152,4 @@ iomodel->FetchData(&types,NULL,NULL,"md.autodiff.independent_object_types"); - M_all = xNew(num_independent_objects); @@ -162,5 +159,5 @@ iomodel->FetchData(&independents_fullmax,&M_par,&N_par,"md.autodiff.independent_max_parameters"); iomodel->FetchData(&control_sizes,NULL,NULL,"md.autodiff.independent_control_sizes"); - + int* start_point = NULL; start_point = xNew(num_independent_objects); @@ -180,5 +177,5 @@ IssmDouble* independents_min = NULL; IssmDouble* independents_max = NULL; - + FieldAndEnumFromCode(&input_enum,&iofieldname,names[i]); @@ -187,5 +184,5 @@ _assert_(independent); _assert_(N==control_sizes[i]); - + independents_min = NULL; independents_min = xNew(M*N); independents_max = NULL; independents_max = xNew(M*N); @@ -206,5 +203,5 @@ xDelete(independents_min); xDelete(independents_max); - + } else{ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 23066) @@ -35,5 +35,5 @@ /*Create elements*/ if(control_analysis && !adolc_analysis)iomodel->FetchData(2,"md.inversion.min_parameters","md.inversion.max_parameters"); - + switch(iomodel->meshelementtype){ case TriaEnum: @@ -123,5 +123,5 @@ break; case MaterialsEnum: - + //we have several types of materials. Retrieve this info first: iomodel->FetchData(&nature,&nnat,&dummy,"md.materials.nature"); @@ -235,5 +235,5 @@ else iomodel->FetchDataToInput(elements,"md.mesh.scale_factor",MeshScaleFactorEnum,1.); if (isoceancoupling) iomodel->FetchData(2,"md.mesh.lat","md.mesh.long"); - + CreateNumberNodeToElementConnectivity(iomodel,solution_type); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 23066) @@ -11,5 +11,5 @@ int i,j; - + DataSet* output_definitions = NULL; int* output_definition_enums = NULL; @@ -65,5 +65,5 @@ else if (output_definition_enums[i]==MisfitEnum){ /*Deal with misfits: {{{*/ - + /*misfit variables: */ int nummisfits; @@ -115,5 +115,4 @@ _error_("misfit weight size not supported yet"); - /*First create a misfit object for that specific string (misfit_model_string_s[j]):*/ 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]))); @@ -159,5 +158,5 @@ else if (output_definition_enums[i]==CfsurfacesquareEnum){ /*Deal with cfsurfacesquare: {{{*/ - + /*cfsurfacesquare variables: */ int num_cfsurfacesquares; @@ -214,5 +213,5 @@ Element* element=xDynamicCast(elements->GetObjectByOffset(k)); - + 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); 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); @@ -251,5 +250,5 @@ else if (output_definition_enums[i]==CfdragcoeffabsgradEnum){ /*Deal with cfdragcoeffabsgrad: {{{*/ - + /*cfdragcoeffabsgrad variables: */ int num_cfdragcoeffabsgrads; @@ -286,5 +285,5 @@ Element* element=xDynamicCast(elements->GetObjectByOffset(k)); - + 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); @@ -313,5 +312,5 @@ else if (output_definition_enums[i]==CfsurfacelogvelEnum){ /*Deal with cfsurfacelogvel: {{{*/ - + /*cfsurfacelogvel variables: */ int num_cfsurfacelogvels; @@ -369,5 +368,5 @@ Element* element=xDynamicCast(elements->GetObjectByOffset(k)); - + 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); 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); @@ -409,5 +408,5 @@ else if (output_definition_enums[i]==NodalvalueEnum){ /*Deal with nodal values: {{{*/ - + /*nodal value variables: */ int numnodalvalues; @@ -428,5 +427,5 @@ 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. } - + /*Free ressources:*/ for(j=0;jAddObject(new Numberedcostfunction(ncf_name_s[j],StringToEnumx(ncf_definitionstring_s[j]),num_cost_functions,cost_function_enums)); } - + /*Free data: */ iomodel->DeleteData(2,"md.numberedcostfunction.name","md.numberedcostfunction.definitionstring"); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23066) @@ -23,5 +23,5 @@ char* fieldname = NULL; IssmDouble time; - + /*parameters for mass flux:*/ int mass_flux_num_profiles = 0; @@ -62,5 +62,4 @@ parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum)); parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1)); - {/*This is specific to ice...*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 23066) @@ -55,5 +55,4 @@ delete analysis; - /* Update counters, because we have created more nodes, loads and * constraints, and ids for objects created in next call to CreateDataSets Index: /issm/trunk-jpl/src/c/modules/NodalValuex/NodalValuex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/NodalValuex/NodalValuex.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/NodalValuex/NodalValuex.cpp (revision 23066) @@ -21,5 +21,5 @@ cpu_found=-1; found=0; - + for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); Index: /issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp (revision 23066) @@ -22,5 +22,5 @@ /*This is the object that we have been chasing for. compute the response and return: */ IssmDouble return_value=definition->Response(femmodel); - + /*cleanup: */ xDelete(name); @@ -31,5 +31,5 @@ xDelete(name); } - + /*If we are here, did not find the definition for this response, not good!: */ _error_("Could not find the response for output definition " << output_string << " because could not find the definition itself!"); @@ -44,5 +44,5 @@ /*Now, go through the output definitions, and retrieve the object which corresponds to our requested response, output_enum: */ for(int i=0;iSize();i++){ - + //Definition* definition=xDynamicCast(output_definitions->GetObjectByOffset(i)); Definition* definition=dynamic_cast(output_definitions->GetObjectByOffset(i)); @@ -53,10 +53,10 @@ /*This is the object that we have been chasing for. compute the response and return: */ IssmDouble return_value=definition->Response(femmodel); - + /*return:*/ return return_value; } } - + /*If we are here, did not find the definition for this response, not good!: */ _error_("Could not find the response for output definition " << EnumToStringx(output_enum) Index: /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 23066) @@ -92,5 +92,5 @@ /* Intermediaries */ int numvertices = element->GetNumberOfVertices(); - + if(element->IsIceInElement()){ for(int i = 0;i(control_type); } Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23066) @@ -107,5 +107,5 @@ xDelete(dzT); xDelete(dzB); - + //---------NEED TO IMPLEMENT A PROPER GRID STRECHING ALGORITHM------------ @@ -202,5 +202,5 @@ gdn=*pgdn; gsp=*pgsp; - + /*only when aIdx = 1 or 2 do we run grainGrowth: */ if(aIdx!=1 & aIdx!=2){ @@ -220,5 +220,4 @@ //set maximum water content by mass to 9 percent (Brun, 1980) for(int i=0;i9.0-Wtol) lwc[i]=9.0; - /* Calculate temperature gradiant across grid cells @@ -244,5 +243,5 @@ // take absolute value of temperature gradient for(int i=0;i(gsz); @@ -348,5 +347,5 @@ //// Inputs // aIdx = albedo method to use - + // Method 0 // aValue = direct input value for albedo, override all changes to albedo @@ -397,5 +396,5 @@ /*Recover pointers: */ a=*pa; - + //some constants: const IssmDouble dSnow = 300; // density of fresh snow [kg m-3] @@ -513,5 +512,5 @@ /* ENGLACIAL THERMODYNAMICS*/ - + /* Description: computes new temperature profile accounting for energy absorption and @@ -537,5 +536,5 @@ // T: grid cell temperature [k] // EC: evaporation/condensation [kg] - + /*intermediary: */ IssmDouble* K = NULL; @@ -580,5 +579,5 @@ IssmDouble dlw=0.0; IssmDouble dT_dlw=0.0; - + /*outputs:*/ IssmDouble EC=0.0; @@ -597,5 +596,5 @@ //initialize Evaporation - Condenstation EC = 0.0; - + // check if all SW applied to surface or distributed throught subsurface // swIdx = length(swf) > 1 @@ -609,5 +608,5 @@ // if V = 0 goes to infinity therfore if V = 0 change if(V<0.01-Dtol)V=0.01; - + // Bulk-transfer coefficient for turbulent fluxes An = pow(0.4,2) / pow(log(Tz/z0),2); // Bulk-transfer coefficient @@ -627,14 +626,14 @@ // THERMAL DIFFUSION COEFFICIENTS - + // A discretization scheme which truncates the Taylor-Series expansion // after the 3rd term is used. See Patankar 1980, Ch. 3&4 - + // discretized heat equation: - + // Tp = (Au*Tu° + Ad*Td° + (Ap-Au-Ad)Tp° + S) / Ap - + // where neighbor coefficients Au, Ap, & Ad are - + // Au = [dz_u/2KP + dz_p/2KE]^-1 // Ad = [dz_d/2KP + dz_d/2KD]^-1 @@ -649,5 +648,5 @@ KD = xNew(m); KP = xNew(m); - + KU[0] = UNDEF; KD[m-1] = UNDEF; @@ -661,5 +660,5 @@ dzU[0]=UNDEF; dzD[m-1]=UNDEF; - + for(int i=1;i(m); @@ -703,16 +702,16 @@ Np[i]= 1.0 - Nu[i] - Nd[i]; } - + // specify boundary conditions: constant flux at bottom Nu[m-1] = 0.0; Np[m-1] = 1.0; - + // zero flux at surface Np[0] = 1.0 - Nd[0]; - + // Create neighbor arrays for diffusion calculations instead of a tridiagonal matrix Nu[0] = 0.0; Nd[m-1] = 0.0; - + /* RADIATIVE FLUXES*/ @@ -720,5 +719,5 @@ sw = xNew(m); for(int i=0;i(m); @@ -746,5 +745,5 @@ // store initial temperature //T_init = T; - + // calculate temperature of snow surface (Ts) // when incoming SW radition is allowed to penetrate the surface, @@ -754,7 +753,7 @@ Ts = (T[0] + T[1])/2.0; Ts = fmin(CtoK,Ts); // don't allow Ts to exceed 273.15 K (0 degC) - + //TURBULENT HEAT FLUX - + // Monin-Obukhov Stability Correction // Reference: @@ -764,7 +763,7 @@ // calculate the Bulk Richardson Number (Ri) Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0)); - + // calculate Monin-Obukhov stability factors 'coefM' and 'coefH' - + // do not allow Ri to exceed 0.19 Ri = fmin(Ri, 0.19); @@ -778,9 +777,9 @@ coefM =pow (1.0-18*Ri,-0.25); } - + // calculate heat/wind 'coef_H' stability factor if (Ri <= -0.03+Ttol) coefH = coefM/1.3; else coefH = coefM; - + //// Sensible Heat // calculate the sensible heat flux [W m-2](Patterson, 1998) @@ -801,5 +800,5 @@ else{ L = LS; // latent heat of sublimation - + // for an ice surface Murphy and Koop, 2005 [Equation 7] eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts); @@ -823,11 +822,11 @@ ulw = - (SB * pow(Ts,4.0)* teValue) * dt ; dT_ulw = ulw / TCs; - + // new grid point temperature - + //SW penetrates surface for(int j=0;j(Td); - /*Assign output pointers:*/ *pEC=EC; @@ -901,5 +899,5 @@ // swf = absorbed shortwave radiation [W m-2] // - + /*outputs: */ IssmDouble* swf=NULL; @@ -910,8 +908,7 @@ swf=xNewZeroInit(m); - // SHORTWAVE FUNCTION if (swIdx == 0) {// all sw radation is absorbed in by the top grid cell - + // calculate surface shortwave radiation fluxes [W m-2] swf[0] = (1.0 - as) * dsw; @@ -948,5 +945,4 @@ } - // spectral albedos: // 0.3 - 0.8um @@ -988,5 +984,4 @@ B2_cum[i+1]=cum2; } - // flux across grid cell boundaries @@ -1015,6 +1010,5 @@ xDelete(Qs1); xDelete(Qs2); - - + } else{ //function of grid cell density @@ -1144,5 +1138,4 @@ if (P > 0.0+Dtol){ - if (T_air <= CtoK+Ttol){ // if snow @@ -1210,5 +1203,5 @@ mass_diff = mass - massinit - P; - + #ifndef _HAVE_ADOLC_ //avoid round operation. only check in forward mode. mass_diff = round(mass_diff * 100.0)/100.0; @@ -1253,5 +1246,5 @@ IssmDouble* M=NULL; int* D=NULL; - + IssmDouble sumM=0.0; IssmDouble sumER=0.0; @@ -1265,5 +1258,5 @@ IssmDouble X=0.0; IssmDouble Wi=0.0; - + IssmDouble Ztot=0.0; IssmDouble T_bot=0.0; @@ -1279,5 +1272,5 @@ IssmDouble EW_bot=0.0; bool top=false; - + IssmDouble* Zcum=NULL; IssmDouble* dzMin2=NULL; @@ -1286,5 +1279,5 @@ int X1=0; int X2=0; - + int D_size = 0; @@ -1303,5 +1296,5 @@ int n=*pn; IssmDouble* R=NULL; - + if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" melt module\n"); @@ -1335,5 +1328,5 @@ // for(int i=0;i(D); - + // calculate new grid lengths for(int i=0;i(n); dzMin2=xNew(n); - + Zcum[0]=dz[0]; // Compute a cumulative depth vector - + for (int i=1;i(D); - + // check if any of the top 10 cell depths are too large X=0; @@ -1657,5 +1648,5 @@ } } - + int j=0; while(j<=X){ @@ -1686,5 +1677,5 @@ // calculate total model depth Ztot = cellsum(dz,n); - + if (Ztot < zMin-Dtol){ // printf("Total depth < zMin %f \n", Ztot); @@ -1692,5 +1683,5 @@ mAdd = m[n-1]+W[n-1]; addE = T[n-1]*m[n-1]*CI + W[n-1]*(LF+CtoK*CI); - + // add a grid cell of the same size and temperature to the bottom dz_bot=dz[n-1]; @@ -1705,7 +1696,7 @@ EI_bot=EI[n-1]; EW_bot=EW[n-1]; - + dz_add=dz_bot; - + newcell(&dz,dz_bot,top,n); newcell(&T,T_bot,top,n); @@ -1727,9 +1718,9 @@ addE = -(T[n-1]*m[n-1]*CI) - (W[n-1]*(LF+CtoK*CI)); dz_add=-(dz[n-1]); - + // remove a grid cell from the bottom D_size=n-1; D=xNew(D_size); - + for(int i=0;i(m); @@ -1784,5 +1775,5 @@ xDelete(Zcum); xDelete(dzMin2); - + /*Assign output pointers:*/ *pM=sumM; @@ -1862,5 +1853,5 @@ IssmDouble* dz=NULL; IssmDouble* d=NULL; - + if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" densification module\n"); @@ -1871,5 +1862,5 @@ // initial mass IssmDouble* mass_init = xNew(m);for(int i=0;i(m-1); @@ -1880,10 +1871,9 @@ obp[0]=0.0; for(int i=1;i 550 [kg m-3] - - + for(int i=0;i division by zero. Abort"); - + if(waterfraction<=w0) Dret=D0; Index: /issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 23066) @@ -66,5 +66,5 @@ /*Get norm of epsprime*/ epsprime_norm = sqrt(epsprime[0]*epsprime[0] + epsprime[1]*epsprime[1] + epsprime[2]*epsprime[2]); - + /*Assign output pointers*/ *pepsprime_norm=epsprime_norm; Index: /issm/trunk-jpl/src/c/shared/Elements/Paterson.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/Paterson.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/shared/Elements/Paterson.cpp (revision 23066) @@ -15,5 +15,4 @@ /*Switch to celsius from Kelvin: */ T=temperature-273.15; - if(T<=-45.0){ Index: /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 23066) @@ -83,9 +83,9 @@ for (iqj = 0; iqj < 12; iqj++){ imonth = ismon[iqj]; - + /*********compute lapse rate ****************/ st=(s-s0t)/1000.; rtlaps=TdiffTime*rlapslgm + (1.-TdiffTime)*rlaps; // lapse rate - + /*********compute Surface temperature *******/ monthlytemperatures[imonth]=monthlytemperatures[imonth] - rtlaps *max(st,sealevTime*0.001); @@ -99,10 +99,10 @@ else { pd = 0.;} - + /******exp des/elev precip reduction*******/ sp=(s-s0p)/1000.-deselcut; // deselev effect is wrt chng in topo if (sp>0.0){q = exp(-desfac*sp);} else {q = 1.0;} - + qmt= qmt + monthlyprec[imonth]*sconv; //*sconv to convert in m of ice equivalent per month qmpt= q*monthlyprec[imonth]*sconv; @@ -114,5 +114,5 @@ // gaussian=T_m, so ndd=-(Tsurf-pdd) if (iqj>5 && iqj<9){ Tsum=Tsum+tstar;} - + if (tstar >= siglim) {pdd = pdd + tstar*deltm;} else if (tstar> -siglim){ @@ -121,5 +121,5 @@ frzndd = frzndd - (tstar-pddsig)*deltm;} else{frzndd = frzndd - tstar*deltm; } - + /*Assign output pointer*/ *(monthlytemperatures+imonth) = monthlytemperatures[imonth]; @@ -150,5 +150,5 @@ } } - + /***** determine PDD factors *****/ if(Tsum< -1.) { @@ -166,10 +166,10 @@ snwmf=0.95*snwmf; smf=0.95*smf; - + /***** compute PDD ablation and refreezing *****/ pddt = pdd *365; snwm = snwmf*pddt; // snow that could have been melted in a year hmx2 = min(h,dfrz); // refreeze active layer max depth: dfrz - + if(snwm < saccu) { water=prect-saccu + snwm; //water=rain + snowmelt Index: /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 23066) @@ -63,5 +63,5 @@ dtemp = &dtemp_static[0]; } - + /* perform the matrix triple product in the order that minimizes the number of multiplies and the temporary space used, noting that @@ -335,5 +335,5 @@ Matrix2x2Determinant(&det,A); if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon"); - + /*Multiplication is faster than divsion, so we multiply by the reciprocal*/ det_reciprocal = 1./det; @@ -427,5 +427,4 @@ *pvy = vy; - }/*}}}*/ @@ -577,14 +576,14 @@ void newcell(IssmDouble** pcell, IssmDouble newvalue, bool top, int m){ /*{{{*/ - + IssmDouble* cell=NULL; IssmDouble* dummy=NULL; - + /*recover pointer: */ cell=*pcell; - + /*reallocate:*/ dummy=xNew(m+1); - + /*copy data:*/ if(top){ @@ -596,8 +595,8 @@ for(int i=0;i(cell); cell=dummy; - + /*assign output pointer:*/ *pcell=cell; @@ -615,5 +614,5 @@ /*input: */ IssmDouble* cell=*pcell; - + /*output: */ IssmDouble* newcell=xNew(nind); @@ -622,5 +621,5 @@ newcell[i]=cell[indices[i]]; } - + /*free allocation:*/ xDelete(cell); @@ -633,5 +632,5 @@ /*input: */ IssmDouble* cell=*pcell; - + /*output: */ IssmDouble* newcell=xNew(m+1); @@ -641,5 +640,5 @@ newcell[i+1]=scale* cell[i]; for(int j=i+2;j(cell); @@ -662,5 +661,5 @@ } va_end ( arguments ); - + _printf_("Echo of cell: \n"); for(int i=0;i(fxmin)) _error_("Function evaluation returned NaN"); - + /*Get f(xmax)*/ _printf0_(" x = "<(X); Index: /issm/trunk-jpl/src/c/shared/Numerics/legendre.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Numerics/legendre.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/shared/Numerics/legendre.cpp (revision 23066) @@ -121,9 +121,9 @@ for i = 2 : n - + v(:,i+1) = ( ( 2 * i - 1 ) * x .* v(:,i) ... - ( i - 1 ) * v(:,i-1) ) ... / ( i ); - + end }}} */ @@ -240,5 +240,5 @@ polynomials of order 0 through N. }}}*/ - + int i,j; @@ -254,5 +254,5 @@ for ( i = 0; i < m; i++ ) { for ( j = 2; j <= n; j++ ) { - + v[j+i*(n+1)] = ( ( IssmDouble ) ( 2 * j - 1 ) * x[i] * v[(j-1)+i*(n+1)] - ( IssmDouble ) ( j - 1 ) * v[(j-2)+i*(n+1)] ) Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (revision 23066) @@ -185,5 +185,5 @@ /*Initial guess*/ VecZeroEntries(udot); - + /*Richardson's iterations*/ for(int i=0;i<5;i++){ @@ -402,5 +402,5 @@ CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type); delete uf; - + /*Go solve lower order solution*/ femmodel->profiler->Start(SOLVER); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp (revision 23066) @@ -20,5 +20,5 @@ int configuration_type; IssmDouble solver_residue_threshold; - + /*solver convergence test: */ IssmDouble nKUF; @@ -39,7 +39,7 @@ Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); femmodel->profiler->Stop(SOLVER); - + /*Check that the solver converged nicely: */ - + //compute KUF = KU - F = K*U - F KU=uf->Duplicate(); Kff->MatMult(uf,KU); @@ -53,5 +53,4 @@ if(!xIsNan(solver_residue_threshold) && solver_residue>solver_residue_threshold)_error_(" solver residue too high!: norm(KU-F)/norm(F)=" << solver_residue << "\n"); - //clean up delete KU; delete KUF; @@ -63,5 +62,5 @@ // } //#endif - + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys; InputUpdateFromSolutionx(femmodel,ug); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp (revision 23066) @@ -61,5 +61,5 @@ count=1; - + for(;;){ delete tf_old;tf_old=tf; Index: /issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp (revision 23065) +++ /issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp (revision 23066) @@ -178,5 +178,4 @@ #endif - recvcounts=xNew(num_procs); displs=xNew(num_procs); @@ -262,5 +261,4 @@ rhs=xNew(pf_M); #endif - recvcounts=xNew(num_procs);