Index: ../trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 23066) @@ -91,7 +91,7 @@ /* Intermediaries */ int numvertices = element->GetNumberOfVertices(); - + if(element->IsIceInElement()){ for(int i = 0;iSetValue(element->vertices[i]->Sid(),1.,INS_VAL); Index: ../trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp =================================================================== --- ../trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp (revision 23066) @@ -21,7 +21,7 @@ /*This is the object that we have been chasing for. compute the response and return: */ IssmDouble return_value=definition->Response(femmodel); - + /*cleanup: */ xDelete(name); @@ -30,7 +30,7 @@ } 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!"); @@ -43,7 +43,7 @@ /*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)); @@ -52,12 +52,12 @@ /*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) <<" ("<FindParam(&migration_style,GroundinglineMigrationEnum); parameters->FindParam(&analysis_type,AnalysisTypeEnum); if(migration_style==NoneEnum) return; - + if(VerboseModule()) _printf0_(" Migrating grounding line\n"); /*Set toolkit to default*/ Index: ../trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp (revision 23066) @@ -77,7 +77,7 @@ if(pdf) df =new Vector(flocalsize,fsize); if(ppf) pf =new Vector(flocalsize,fsize); } - + /*Free ressources: */ xDelete(toolkittype); Index: ../trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp (revision 23066) @@ -27,7 +27,6 @@ int size = 0; for(int i=0;i(size); @@ -59,7 +58,6 @@ parameters->FindParam(&num_controls,InversionNumControlParametersEnum); parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); - /*2. Allocate vector*/ vector=new Vector(num_controls*vertices->NumberOfVertices()); Index: ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp (revision 23066) @@ -165,10 +165,10 @@ edge3=SegmentIntersect(&gamma1,&gamma2, xel,yel,xsegment,ysegment); /*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)); }*/ @@ -198,7 +198,6 @@ } 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]); /* _error_("error: a line cannot go through 3 different vertices!");*/ @@ -225,7 +224,7 @@ } } 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]))); }*/ @@ -250,7 +249,7 @@ if (IsIdenticalNode(xnodes[0],ynodes[0],xsegment[0],ysegment[0],tolerance) || 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; if(edge1==IntersectEnum){coord2=alpha1;} Index: ../trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp =================================================================== --- ../trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp (revision 23066) @@ -38,7 +38,7 @@ IssmDouble d0=1.e+10; IssmDouble xi,yi; - + //recover vertex position: xi=x[i]; yi=y[i]; Index: ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 23066) @@ -83,7 +83,7 @@ int interpolation_type; /*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)); Input* input=element->GetInput(name); Index: ../trunk-jpl/src/c/modules/NodalValuex/NodalValuex.cpp =================================================================== --- ../trunk-jpl/src/c/modules/NodalValuex/NodalValuex.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/NodalValuex/NodalValuex.cpp (revision 23066) @@ -20,7 +20,7 @@ *element, figure out if they hold the vertex, and the data. If so, return it: */ cpu_found=-1; found=0; - + for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); found=element->NodalValue(&value,index,natureofdataenum); Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 23066) @@ -8,7 +8,6 @@ #include "../../MeshPartitionx/MeshPartitionx.h" #include "../ModelProcessorx.h" - #if !defined(_HAVE_ADOLC_) void UpdateElementsAndMaterialsControl(Elements* elements,Parameters* parameters,Materials* materials, IoModel* iomodel){ /*Intermediary*/ @@ -21,7 +20,6 @@ char **controls = NULL; char **cost_functions = NULL; - /*Fetch parameters: */ iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol"); if(control_analysis) iomodel->FindConstant(&num_controls,"md.inversion.num_control_parameters"); @@ -28,7 +26,7 @@ /*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"); if(num_controls<1) _error_("no controls found"); @@ -47,7 +45,7 @@ } 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"); for(int i=0;iFetchData(&types,NULL,NULL,"md.autodiff.independent_object_types"); - M_all = xNew(num_independent_objects); /*create independent objects, and at the same time, fetch the corresponding independent variables, @@ -161,7 +158,7 @@ iomodel->FetchData(&independents_fullmin,&M_par,&N_par,"md.autodiff.independent_min_parameters"); 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); int counter = 0; @@ -179,7 +176,7 @@ int input_enum; IssmDouble* independents_min = NULL; IssmDouble* independents_max = NULL; - + FieldAndEnumFromCode(&input_enum,&iofieldname,names[i]); /*Fetch required data*/ @@ -186,7 +183,7 @@ iomodel->FetchData(&independent,&M,&N,iofieldname); _assert_(independent); _assert_(N==control_sizes[i]); - + independents_min = NULL; independents_min = xNew(M*N); independents_max = NULL; independents_max = xNew(M*N); for(int m=0;m(independent); xDelete(independents_min); xDelete(independents_max); - + } else{ _error_("Independent cannot be of size " << types[i]); Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 23066) @@ -68,11 +68,11 @@ /*Intermediaries*/ 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"); _assert_(num_independent_objects>0); @@ -95,10 +95,10 @@ } 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;i(names[i]); Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 23066) @@ -10,7 +10,7 @@ void CreateOutputDefinitions(Elements* elements, Parameters* parameters,IoModel* iomodel){ int i,j; - + DataSet* output_definitions = NULL; int* output_definition_enums = NULL; int num_output_definitions; @@ -64,7 +64,7 @@ } else if (output_definition_enums[i]==MisfitEnum){ /*Deal with misfits: {{{*/ - + /*misfit variables: */ int nummisfits; char** misfit_name_s = NULL; @@ -114,7 +114,6 @@ else _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]))); @@ -158,7 +157,7 @@ } else if (output_definition_enums[i]==CfsurfacesquareEnum){ /*Deal with cfsurfacesquare: {{{*/ - + /*cfsurfacesquare variables: */ int num_cfsurfacesquares; char** cfsurfacesquare_name_s = NULL; @@ -213,7 +212,7 @@ for(int k=0;kSize();k++){ 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); @@ -250,7 +249,7 @@ } else if (output_definition_enums[i]==CfdragcoeffabsgradEnum){ /*Deal with cfdragcoeffabsgrad: {{{*/ - + /*cfdragcoeffabsgrad variables: */ int num_cfdragcoeffabsgrads; char** cfdragcoeffabsgrad_name_s = NULL; @@ -285,7 +284,7 @@ for(int k=0;kSize();k++){ 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); } @@ -312,7 +311,7 @@ } else if (output_definition_enums[i]==CfsurfacelogvelEnum){ /*Deal with cfsurfacelogvel: {{{*/ - + /*cfsurfacelogvel variables: */ int num_cfsurfacelogvels; char** cfsurfacelogvel_name = NULL; @@ -368,7 +367,7 @@ for(int k=0;kSize();k++){ 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); element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_weights[j], iomodel,cfsurfacelogvel_weights_M[j],cfsurfacelogvel_weights_N[j],weight_vector_type,StringToEnumx(cfsurfacelogvel_weightstring[j]),7,WeightsSurfaceObservationEnum); @@ -408,7 +407,7 @@ } else if (output_definition_enums[i]==NodalvalueEnum){ /*Deal with nodal values: {{{*/ - + /*nodal value variables: */ int numnodalvalues; char** nodalvalue_name_s = NULL; @@ -427,7 +426,7 @@ /*First create a nodalvalue object for that specific enum (nodalvalue_model_enum_s[j]):*/ 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"); xDelete(cost_function_enums); Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23066) @@ -22,7 +22,7 @@ char** requestedoutputs = NULL; char* fieldname = NULL; IssmDouble time; - + /*parameters for mass flux:*/ int mass_flux_num_profiles = 0; bool qmu_mass_flux_present = false; @@ -61,7 +61,6 @@ parameters->AddObject(iomodel->CopyConstantObject("md.inversion.type",InversionTypeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum)); parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1)); - {/*This is specific to ice...*/ parameters->AddObject(iomodel->CopyConstantObject("md.mesh.elementtype",MeshElementtypeEnum)); Index: ../trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 23066) @@ -54,7 +54,6 @@ analysis->UpdateElements(elements,iomodel,i,analysis_enum); delete analysis; - /* Update counters, because we have created more nodes, loads and * constraints, and ids for objects created in next call to CreateDataSets * will need to start at the end of the updated counters: */ Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 23066) @@ -34,7 +34,7 @@ /*Create elements*/ if(control_analysis && !adolc_analysis)iomodel->FetchData(2,"md.inversion.min_parameters","md.inversion.max_parameters"); - + switch(iomodel->meshelementtype){ case TriaEnum: for(i=0;inumberofelements;i++){ @@ -122,7 +122,7 @@ } break; case MaterialsEnum: - + //we have several types of materials. Retrieve this info first: iomodel->FetchData(&nature,&nnat,&dummy,"md.materials.nature"); @@ -234,7 +234,7 @@ if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->FetchData(3,"md.mesh.lat","md.mesh.long","md.mesh.r"); 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); int lid = 0; Index: ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp (revision 23066) @@ -31,7 +31,6 @@ offset += M[i]*N[i]; } - xDelete(control_type); } else{ Index: ../trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp =================================================================== --- ../trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 23066) @@ -66,4 +66,3 @@ } } /*}}}*/ - Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 23066) @@ -6,7 +6,6 @@ #include "../../shared/shared.h" #include "../../toolkits/toolkits.h" - void SmbGradientsx(FemModel* femmodel){/*{{{*/ // void SurfaceMassBalancex(hd,agd,ni){ @@ -317,7 +316,6 @@ correct for number of seconds in a year [s/yr]*/ smb = smb/yts + anomaly; - /*Update array accordingly*/ smblist[v] = smb; Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23066) @@ -106,7 +106,7 @@ /*Free ressouces:*/ xDelete(dzT); xDelete(dzB); - + //---------NEED TO IMPLEMENT A PROPER GRID STRECHING ALGORITHM------------ /*assign ouput pointers: */ @@ -201,7 +201,7 @@ re=*pre; gdn=*pgdn; gsp=*pgsp; - + /*only when aIdx = 1 or 2 do we run grainGrowth: */ if(aIdx!=1 & aIdx!=2){ /*come out as we came in: */ @@ -220,7 +220,6 @@ //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 * Returns the avereage gradient across the upper and lower grid cell */ @@ -243,7 +242,7 @@ // take absolute value of temperature gradient for(int i=0;i0.0+Gdntol){ @@ -323,7 +322,7 @@ //convert grain size back to effective grain radius: re[i]=gsz[i]/2.0; } - + /*Free ressources:*/ xDelete(gsz); xDelete(dT); @@ -347,7 +346,7 @@ //// Inputs // aIdx = albedo method to use - + // Method 0 // aValue = direct input value for albedo, override all changes to albedo @@ -396,7 +395,7 @@ /*Recover pointers: */ a=*pa; - + //some constants: const IssmDouble dSnow = 300; // density of fresh snow [kg m-3] @@ -512,7 +511,7 @@ void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid) { /*{{{*/ /* ENGLACIAL THERMODYNAMICS*/ - + /* Description: computes new temperature profile accounting for energy absorption and thermal diffusion.*/ @@ -536,7 +535,7 @@ // OUTPUTS // T: grid cell temperature [k] // EC: evaporation/condensation [kg] - + /*intermediary: */ IssmDouble* K = NULL; IssmDouble* KU = NULL; @@ -579,7 +578,7 @@ IssmDouble dT_ulw=0.0; IssmDouble dlw=0.0; IssmDouble dT_dlw=0.0; - + /*outputs:*/ IssmDouble EC=0.0; IssmDouble* T=*pT; @@ -596,7 +595,7 @@ //initialize Evaporation - Condenstation EC = 0.0; - + // check if all SW applied to surface or distributed throught subsurface // swIdx = length(swf) > 1 @@ -608,7 +607,7 @@ // 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 C = An * dAir * V; // shf & lhf common coefficient @@ -626,16 +625,16 @@ for(int i=0;i=dIce-Dtol) K[i] = 9.828 * exp(-5.7E-3*T[i]); // 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 // Ap = d*CI*dz/Dt @@ -648,7 +647,7 @@ KU = xNew(m); KD = xNew(m); KP = xNew(m); - + KU[0] = UNDEF; KD[m-1] = UNDEF; for(int i=1;i(m); dzU[0]=UNDEF; dzD[m-1]=UNDEF; - + for(int i=1;i(m); Nd = xNew(m); @@ -702,24 +701,24 @@ Nd[i] = Ad[i] / Ap[i]; 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*/ // energy supplied by shortwave radiation [J] sw = xNew(m); for(int i=0;i(m); for(int i=0;i(Tu); xDelete(Td); - /*Assign output pointers:*/ *pEC=EC; *pT=T; @@ -900,7 +898,7 @@ // Outputs // swf = absorbed shortwave radiation [W m-2] // - + /*outputs: */ IssmDouble* swf=NULL; @@ -909,10 +907,9 @@ /*Initialize and allocate: */ 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; } @@ -947,7 +944,6 @@ B2_cum[i]=1.0; } - // spectral albedos: // 0.3 - 0.8um IssmDouble a0 = fmin(0.98, 1.0 - 1.58 *pow(gsz[0],0.5)); @@ -988,7 +984,6 @@ B2_cum[i+1]=cum2; } - // flux across grid cell boundaries Qs1 = xNew(m+1); Qs2 = xNew(m+1); @@ -1014,8 +1009,7 @@ xDelete(exp2); xDelete(Qs1); xDelete(Qs2); - - + } else{ //function of grid cell density @@ -1143,7 +1137,6 @@ for(int i=0;i 0.0+Dtol){ - if (T_air <= CtoK+Ttol){ // if snow @@ -1209,7 +1202,7 @@ mass=0; for(int i=0;i 0) _error_("mass not conserved in accumulation function"); @@ -1252,7 +1245,7 @@ IssmDouble* EW=NULL; IssmDouble* M=NULL; int* D=NULL; - + IssmDouble sumM=0.0; IssmDouble sumER=0.0; IssmDouble addE=0.0; @@ -1264,7 +1257,7 @@ IssmDouble dm=0.0; IssmDouble X=0.0; IssmDouble Wi=0.0; - + IssmDouble Ztot=0.0; IssmDouble T_bot=0.0; IssmDouble m_bot=0.0; @@ -1278,7 +1271,7 @@ IssmDouble EI_bot=0.0; IssmDouble EW_bot=0.0; bool top=false; - + IssmDouble* Zcum=NULL; IssmDouble* dzMin2=NULL; IssmDouble zY2=1.025; @@ -1285,7 +1278,7 @@ bool lastCellFlag = false; int X1=0; int X2=0; - + int D_size = 0; /*outputs:*/ @@ -1302,7 +1295,7 @@ IssmDouble* gsp=*pgsp; int n=*pn; IssmDouble* R=NULL; - + if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" melt module\n"); //// INITIALIZATION @@ -1334,7 +1327,7 @@ // new grid point center temperature, T [K] // 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; for(int i=9;i>=0;i--){ @@ -1656,7 +1647,7 @@ break; } } - + int j=0; while(j<=X){ if (dz[j] > dzMin*2.0+Dtol){ @@ -1685,13 +1676,13 @@ // INTERPRETATION // calculate total model depth Ztot = cellsum(dz,n); - + if (Ztot < zMin-Dtol){ // printf("Total depth < zMin %f \n", Ztot); // mass and energy to be added 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]; T_bot=T[n-1]; @@ -1704,9 +1695,9 @@ gsp_bot=gsp[n-1]; 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); newcell(&W,W_bot,top,n); @@ -1726,11 +1717,11 @@ mAdd = -(m[n-1]+W[n-1]); 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); if(EI)xDelete(EI); @@ -1783,7 +1774,7 @@ if(M)xDelete(M); xDelete(Zcum); xDelete(dzMin2); - + /*Assign output pointers:*/ *pM=sumM; *pR=Rsum; @@ -1861,7 +1852,7 @@ /*output: */ IssmDouble* dz=NULL; IssmDouble* d=NULL; - + if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" densification module\n"); /*Recover pointers: */ @@ -1870,7 +1861,7 @@ // initial mass IssmDouble* mass_init = xNew(m);for(int i=0;i(m-1); cumdz[0]=dz[0]; @@ -1879,12 +1870,11 @@ IssmDouble* obp = xNew(m); obp[0]=0.0; for(int i=1;i 550 [kg m-3] - - + for(int i=0;iparameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum); femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); - + IssmDouble* boxareas=xNewZeroInit(num_basins*maxbox); for(int i=0;ielements->Size();i++){ @@ -121,7 +121,7 @@ IssmDouble* sumareas =xNew(num_basins*maxbox); 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: ../trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp (revision 23065) +++ ../trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp (revision 23066) @@ -51,7 +51,7 @@ double x1,y1; double x2,y2; double mind; - + for(int i=i0;i* ys = NULL; int configuration_type; IssmDouble solver_residue_threshold; - + /*solver convergence test: */ IssmDouble nKUF; IssmDouble nF; @@ -38,9 +38,9 @@ femmodel->profiler->Start(SOLVER); 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); KUF=KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.0); @@ -52,7 +52,6 @@ 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; delete Kff; delete pf; delete df; @@ -62,7 +61,7 @@ // ADOLC_DUMP_MACRO(uf->svector->vector[i]); // } //#endif - + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys; InputUpdateFromSolutionx(femmodel,ug); delete ug; Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp =================================================================== --- ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (revision 23065) +++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (revision 23066) @@ -184,7 +184,7 @@ /*Initial guess*/ VecZeroEntries(udot); - + /*Richardson's iterations*/ for(int i=0;i<5;i++){ MatMult(Mc,udot,temp1); @@ -401,7 +401,7 @@ /*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */ 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); SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp =================================================================== --- ../trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp (revision 23065) +++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp (revision 23066) @@ -60,7 +60,7 @@ } count=1; - + for(;;){ delete tf_old;tf_old=tf; Index: ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp =================================================================== --- ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp (revision 23065) +++ ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp (revision 23066) @@ -177,7 +177,6 @@ rhs=xNew(pf_M); #endif - recvcounts=xNew(num_procs); displs=xNew(num_procs); @@ -262,7 +261,6 @@ rhs=xNew(pf_M); #endif - recvcounts=xNew(num_procs); displs=xNew(num_procs); Index: ../trunk-jpl/src/c/cores/controlad_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/controlad_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/controlad_core.cpp (revision 23066) @@ -104,7 +104,7 @@ case 7: _printf0_(" > 0 or <0\n"); break; default: _printf0_(" Unknown end condition\n"); } - + /*Save results:*/ femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,G,n,1,0,0.0)); femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffXpEnum,X,intn,1,0,0.0)); @@ -129,7 +129,7 @@ int intn; 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 in one run of the solution core. So first recover the filenames required for the FemModel constructor, then call a new ad tailored constructor:*/ @@ -142,7 +142,6 @@ femmodel=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, restartfilename,IssmComm::GetComm(), femmodel->solution_type,NULL); - /*Get initial guess:*/ femmodel->parameters->FindParam(&Xd,&intn,AutodiffXpEnum); X=xNew(intn); for(i=0;i(Xd[i]); @@ -174,7 +173,7 @@ FemModel *femmodel = NULL; IssmDouble pfd; int i; - + /*Recover Femmodel*/ femmodel = (FemModel*)dzs; @@ -194,7 +193,7 @@ IssmDouble* Jlist = NULL; femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast(pfd); _printf0_("f(x) = "<(Jlist); xDelete(G2); - + xDelete(rootpath); xDelete(inputfilename); xDelete(outputfilename); @@ -250,7 +249,7 @@ FemModel *femmodelad = NULL; IssmDouble pfd; int i; - + /*Recover Femmodel*/ femmodel = (FemModel*)dzs; @@ -270,7 +269,7 @@ 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); @@ -283,7 +282,7 @@ IssmDouble* Jlist = NULL; femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast(pfd); _printf0_("f(x) = "<(Jlist); xDelete(G2); - xDelete(rootpath); xDelete(inputfilename); xDelete(outputfilename); Index: ../trunk-jpl/src/c/cores/gia_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/gia_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/gia_core.cpp (revision 23066) @@ -52,7 +52,7 @@ int outputs[2] = {GiaWEnum,GiadWdtEnum}; femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2); } - + xDelete(x); xDelete(y); } Index: ../trunk-jpl/src/c/cores/bmb_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/bmb_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/bmb_core.cpp (revision 23066) @@ -11,7 +11,6 @@ void bmb_core(FemModel* femmodel){ - /*First, get BMB model from parameters*/ int basalforcing_model; bool isplume = false; Index: ../trunk-jpl/src/c/cores/dakota_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/dakota_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/dakota_core.cpp (revision 23066) @@ -234,7 +234,6 @@ /*}}}*/ void dakota_core(FemModel* femmodel){ /*{{{*/ - int my_rank; char *dakota_input_file = NULL; char *dakota_output_file = NULL; Index: ../trunk-jpl/src/c/cores/esa_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/esa_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/esa_core.cpp (revision 23066) @@ -22,7 +22,7 @@ int solution_type; int numoutputs = 0; char **requested_outputs = NULL; - + /*additional parameters: */ int gsize; bool spherical=true; @@ -79,7 +79,7 @@ U_east = new Vector(gsize); U_x = new Vector(gsize); U_y = new Vector(gsize); - + /*call the geodetic main modlule:*/ if(domaintype==Domain3DsurfaceEnum){ femmodel->EsaGeodetic3D(U_radial,U_north,U_east,latitude,longitude,radius,xx,yy,zz); @@ -94,13 +94,13 @@ InputUpdateFromVectorx(femmodel,U_radial,EsaUmotionEnum,VertexSIdEnum); // radial displacement 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"); femmodel->parameters->FindParam(&requested_outputs,&numoutputs,EsaRequestedOutputsEnum); femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); } - + if(solution_type==EsaSolutionEnum)femmodel->RequestedDependentsx(); /*Free ressources:*/ @@ -114,4 +114,3 @@ } /*}}}*/ - Index: ../trunk-jpl/src/c/cores/masstransport_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/masstransport_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/masstransport_core.cpp (revision 23066) @@ -30,7 +30,7 @@ femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum); femmodel->parameters->FindParam(&stabilization,MasstransportStabilizationEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum); - + if(VerboseSolution()) _printf0_(" computing mass transport\n"); /*Transport mass or free surface*/ Index: ../trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 23066) @@ -12,17 +12,16 @@ /*cores:*/ 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; @@ -40,7 +39,7 @@ /*Run steric core for sure:*/ steric_core(femmodel); - + /*Save results: */ if(save_results){ int numoutputs = 0; @@ -58,7 +57,6 @@ /*}}}*/ void geodetic_core(FemModel* femmodel){ /*{{{*/ - /*variables:*/ Vector *RSLg = NULL; Vector *RSLg_rate = NULL; @@ -87,7 +85,7 @@ /*Should we even be here?:*/ femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); if(!geodetic)return; - + /*Verbose: */ if(VerboseSolution()) _printf0_(" computing geodetic sea level rise\n"); @@ -154,7 +152,7 @@ /*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)); U_esa_rate=U_esa->Duplicate(); U_esa->Copy(U_esa_rate); U_esa_rate->Scale(1/(dt*frequency)); @@ -161,7 +159,7 @@ N_gia_rate=N_gia->Duplicate(); N_gia->Copy(N_gia_rate); N_gia_rate->Scale(1/(dt*frequency)); 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); InputUpdateFromVectorx(femmodel,N_esa_rate,SealevelNEsaRateEnum,VertexSIdEnum); @@ -181,7 +179,6 @@ } /*}}}*/ } - if(iscoupler){ /*transfer sea level back to ice caps:*/ TransferSealevel(femmodel,SealevelNEsaRateEnum); @@ -227,7 +224,7 @@ Vector *N_esa_rate= NULL; Vector *U_gia_rate= NULL; Vector *N_gia_rate= NULL; - + /*parameters: */ bool isslr=0; int solution_type; @@ -242,7 +239,7 @@ /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/ if(solution_type==SealevelriseSolutionEnum)isslr=1; - + /*Should we be here?:*/ if(!isslr)return; @@ -290,7 +287,7 @@ } /*}}}*/ Vector* sealevelrise_core_eustatic(FemModel* femmodel){ /*{{{*/ - + /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/ Vector *RSLgi = NULL; @@ -317,7 +314,7 @@ /*Figure out size of g-set deflection vector and allocate solution vector: */ gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); - + /*Initialize:*/ RSLgi = new Vector(gsize); @@ -330,11 +327,10 @@ /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the * 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:*/ xDelete(latitude); xDelete(longitude); @@ -346,7 +342,6 @@ /*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; Vector *RSLg_old = NULL; @@ -379,7 +374,7 @@ femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum); femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum); - + /*computational flag: */ femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum); @@ -388,7 +383,7 @@ /*Figure out size of g-set deflection vector and allocate solution vector: */ gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); - + /*Initialize:*/ RSLg = new Vector(gsize); RSLg->Assemble(); @@ -399,7 +394,7 @@ count=1; converged=false; - + /*Start loop: */ for(;;){ @@ -412,7 +407,7 @@ /*call the non eustatic module: */ femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, latitude, longitude, radius,verboseconvolution,loop); - + /*assemble solution vector: */ RSLgo->Assemble(); @@ -422,7 +417,7 @@ RSLgo_rot = new Vector(gsize); RSLgo_rot->Assemble(); 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)); femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz)); @@ -430,10 +425,10 @@ 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); RSLg->Shift(-RSLgo_oceanaverage); @@ -454,10 +449,10 @@ converged=true; break; } - + /*some minor verbosing adjustment:*/ if(count>1)verboseconvolution=false; - + } if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n"); @@ -473,7 +468,7 @@ Vector *U_esa = NULL; Vector *U_north_esa = NULL; Vector *U_east_esa = NULL; - + /*parameters: */ int configuration_type; int gsize; @@ -506,7 +501,7 @@ /*retrieve geometric information: */ 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); @@ -531,7 +526,7 @@ /*variables:*/ Vector *U_gia = NULL; Vector *N_gia = NULL; - + /*parameters:*/ int frequency; IssmDouble dt; @@ -569,7 +564,7 @@ IssmDouble* forcing=NULL; Vector* forcingglobal=NULL; int* nvs=NULL; - + /*transition vectors:*/ IssmDouble** transitions=NULL; int ntransitions; @@ -576,7 +571,7 @@ int* transitions_m=NULL; int* transitions_n=NULL; int nv; - + /*communicators:*/ ISSM_MPI_Comm tocomm; ISSM_MPI_Comm* fromcomms=NULL; @@ -590,7 +585,7 @@ femmodel->parameters->FindParam(&earthid,EarthIdEnum); 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){ GenericParam* parcoms = dynamic_cast*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum)); @@ -619,7 +614,7 @@ forcings[i]=xNew(nvs[i]); ISSM_MPI_Recv(forcings[i], nvs[i], ISSM_MPI_DOUBLE, 0,i, fromcomms[i], &status); } - + } else{ ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, modelid, tocomm); @@ -630,7 +625,7 @@ /*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: */ nv=femmodel->vertices->NumberOfVertices(); @@ -652,7 +647,6 @@ /*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 * ice cap: */ forcingglobal->SetValues(M,index,forcingfromcap,ADD_VAL); @@ -662,7 +656,7 @@ /*Assemble vector:*/ forcingglobal->Assemble(); - + /*Plug into elements:*/ InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum); } @@ -695,7 +689,7 @@ /*forcing being transferred from earth to ice caps: */ IssmDouble* forcing=NULL; IssmDouble* forcingglobal=NULL; - + /*transition vectors:*/ IssmDouble** transitions=NULL; int ntransitions; @@ -702,7 +696,7 @@ int* transitions_m=NULL; int* transitions_n=NULL; int nv; - + /*communicators:*/ ISSM_MPI_Comm fromcomm; ISSM_MPI_Comm* tocomms=NULL; @@ -717,7 +711,7 @@ femmodel->parameters->FindParam(&earthid,EarthIdEnum); 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){ GenericParam* parcoms = dynamic_cast*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum)); @@ -732,7 +726,6 @@ //femmodel->parameters->FindParam((int*)(&fromcomm), IcecapToEarthCommEnum); } - /*Retrieve sea-level on earth model: */ if(modelid==earthid){ nv=femmodel->vertices->NumberOfVertices(); @@ -741,12 +734,12 @@ /*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!"); for(int i=0;i *deltathickness = NULL; Vector *cumdeltathickness = NULL; int nv; - + if(VerboseSolution()) _printf0_(" computing earth mass transport\n"); /*This mass transport module for the Earth is because we might have thickness variations as spcs @@ -839,7 +832,7 @@ } /*}}}*/ void slrconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ - + bool converged=true; IssmDouble ndS,nS; Vector *dRSLg = NULL; @@ -847,15 +840,14 @@ //compute norm(du) and norm(u) if requested 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 delete dRSLg; Index: ../trunk-jpl/src/c/cores/stressbalance_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 23066) @@ -21,7 +21,6 @@ int numoutputs = 0; char **requested_outputs = NULL; Analysis *analysis = NULL; - /* recover parameters:*/ femmodel->parameters->FindParam(&domaintype,DomainTypeEnum); @@ -35,7 +34,7 @@ femmodel->parameters->FindParam(&numoutputs,StressbalanceNumRequestedOutputsEnum); femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,StressbalanceRequestedOutputsEnum); - + if(VerboseSolution()) _printf0_(" computing new velocity\n"); /*Compute slopes if necessary */ @@ -78,12 +77,11 @@ delete analysis; } - if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); } - + if(solution_type==StressbalanceSolutionEnum && !control_analysis)femmodel->RequestedDependentsx(); /*Free ressources:*/ Index: ../trunk-jpl/src/c/cores/adgradient_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/adgradient_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/adgradient_core.cpp (revision 23066) @@ -38,7 +38,7 @@ /*intermediary: */ IssmPDouble *aWeightVector=NULL; IssmPDouble *weightVectorTimesJac=NULL; - + /*output: */ IssmPDouble *totalgradient=NULL; @@ -52,7 +52,7 @@ /*First, stop tracing: */ trace_off(); - + if(VerboseAutodiff()){ /*{{{*/ size_t tape_stats[15]; tapestats(my_rank,tape_stats); //reading of tape statistics @@ -95,11 +95,11 @@ } 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; @@ -108,10 +108,10 @@ if (my_rank!=0){ num_dependents=0; num_independents=0; } - + /*retrieve state variable: */ femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum); - + /* driver argument */ xp=xNew(num_independents); for(i=0;i(num_independents); - + for(aDepIndex=0;aDepIndex(num_dependents); if (my_rank==0) aWeightVector[aDepIndex]=1.0; - + /*initialize output gradient: */ weightVectorTimesJac=xNew(num_independents); @@ -161,12 +161,12 @@ xDelete(weightVectorTimesJac); 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); xDelete(anEDF_for_solverx_p->dp_X); Index: ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 23066) @@ -125,7 +125,7 @@ ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); SetControlInputsFromVectorx(femmodel,aX); xDelete(aX); - + /*Compute solution (forward)*/ void (*solutioncore)(FemModel*)=NULL; CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); @@ -391,7 +391,7 @@ /*Optimization criterions*/ long niter = long(maxsteps); /*Maximum number of iterations*/ long nsim = long(maxiter);/*Maximum number of function calls*/ - + /*Get initial guess*/ Vector *Xpetsc = NULL; @@ -400,7 +400,7 @@ Xpetsc->GetSize(&intn); delete Xpetsc; //_assert_(intn==numberofvertices*num_controls); - + /*Get problem dimension and initialize gradient and initial guess*/ long n = long(intn); G = xNew(n); @@ -414,7 +414,7 @@ } N_add+=N[c]; } - + /*Allocate m1qn3 working arrays (see documentation)*/ long m = 100; long ndz = 4*n+m*(2*n+1); @@ -470,7 +470,7 @@ } N_add +=N[c]; } - + /*Set X as our new control*/ IssmDouble* aX=xNew(intn); IssmDouble* aG=xNew(intn); @@ -482,10 +482,9 @@ ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG); SetControlInputsFromVectorx(femmodel,aX); - + xDelete(aX); - if (solution_type == TransientSolutionEnum){ int step = 1; femmodel->parameters->SetParam(step,StepEnum); @@ -505,7 +504,7 @@ /*Add to results*/ femmodel->results->AddObject(G_output); femmodel->results->AddObject(X_output); - + offset += N[i]*numberofvertices; } } Index: ../trunk-jpl/src/c/cores/love_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/love_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/love_core.cpp (revision 23066) @@ -27,12 +27,12 @@ /*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); femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum); _assert_(nfreq==dummy); @@ -44,7 +44,7 @@ femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum); 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; for (int i=femmodel->materials->Size()-1;i>=0;i--){ Index: ../trunk-jpl/src/c/cores/damage_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/damage_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/damage_core.cpp (revision 23066) @@ -18,7 +18,6 @@ int numoutputs = 0; char **requested_outputs = NULL; - //first recover parameters common to all solutions femmodel->parameters->FindParam(&save_results,SaveResultsEnum); femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); Index: ../trunk-jpl/src/c/cores/smb_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/smb_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/smb_core.cpp (revision 23066) @@ -28,9 +28,9 @@ femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); 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); delete analysis; Index: ../trunk-jpl/src/c/cores/ad_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/ad_core.cpp (revision 23065) +++ ../trunk-jpl/src/c/cores/ad_core.cpp (revision 23066) @@ -45,7 +45,7 @@ /*First, stop tracing: */ trace_off(); - + /*Print tape statistics so that user can kill this run if something is off already:*/ if(VerboseAutodiff()){ /*{{{*/ tapestats(my_rank,tape_stats); //reading of tape statistics @@ -92,7 +92,7 @@ /*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; @@ -100,7 +100,7 @@ if (my_rank!=0){ num_dependents=0; num_independents=0; } - + /*retrieve state variable: */ femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum); @@ -116,7 +116,6 @@ /*Branch according to AD driver: */ femmodel->parameters->FindParam(&driver,AutodiffDriverEnum); - if (strcmp(driver,"fos_forward")==0){ /*{{{*/ int anIndepIndex; @@ -140,7 +139,6 @@ anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx; #endif - /*call driver: */ fos_forward(my_rank,num_dependents,num_independents, 0, xp, tangentDir, theOutput, jacTimesTangentDir ); @@ -184,7 +182,6 @@ #endif // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; - /*seed matrix: */ seed=xNewZeroInit(num_independents,tangentDirNum); @@ -243,7 +240,6 @@ anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF; #endif - /*call driver: */ fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac ); @@ -284,7 +280,6 @@ anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; #endif - /*seed matrix: */ weights=xNewZeroInit(weightNum,num_dependents); @@ -316,7 +311,6 @@ } /*}}}*/ else _error_("driver: " << driver << " not yet supported!"); - if(VerboseAutodiff())_printf0_(" end AD core\n"); /*Free resources: */ Index: ../trunk-jpl/src/c/classes/Materials/Matestar.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 23066) @@ -108,7 +108,7 @@ void Matestar::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook(); - + MARSHALLING_ENUM(MatestarEnum); MARSHALLING(mid); this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); Index: ../trunk-jpl/src/c/classes/Materials/Matlitho.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matlitho.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Materials/Matlitho.cpp (revision 23066) @@ -36,13 +36,13 @@ 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); @@ -60,17 +60,17 @@ 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); xDelete(issolidd); @@ -77,7 +77,7 @@ } /*}}}*/ Matlitho::~Matlitho(){/*{{{*/ - + xDelete(radius); xDelete(viscosity); xDelete(lame_lambda); Index: ../trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 23066) @@ -134,7 +134,7 @@ _printf_(" note: helement not printed to avoid recursion.\n"); //if(helement) helement->DeepEcho(); //else _printf_(" helement = NULL\n"); - + _printf_(" element:\n"); _printf_(" note: element not printed to avoid recursion.\n"); //if(element) element->DeepEcho(); @@ -142,12 +142,12 @@ } /*}}}*/ 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*/ _printf_(" helement:\n"); @@ -154,7 +154,7 @@ _printf_(" note: helement not printed to avoid recursion.\n"); //if(helement) helement->Echo(); //else _printf_(" helement = NULL\n"); - + _printf_(" element:\n"); _printf_(" note: element not printed to avoid recursion.\n"); //if(element) element->Echo(); @@ -166,7 +166,7 @@ void Matice::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook(); - + MARSHALLING_ENUM(MaticeEnum); MARSHALLING(mid); MARSHALLING(isdamaged); @@ -540,7 +540,6 @@ /*input strain rate: */ IssmDouble exx,eyy,exy,exz,eyz; - if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && (epsilon[3]==0) && (epsilon[4]==0)){ mu_prime=0.5*pow(10.,14); Index: ../trunk-jpl/src/c/classes/Numberedcostfunction.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Numberedcostfunction.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Numberedcostfunction.cpp (revision 23066) @@ -8,7 +8,6 @@ #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" #endif - /*Headers:*/ //#include "./Definition.h" //#include "../datastructures/datastructures.h" @@ -29,7 +28,7 @@ #include "../modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h" /*}}}*/ - + /*Numberedcostfunction constructors, destructors :*/ Numberedcostfunction::Numberedcostfunction(){/*{{{*/ @@ -109,7 +108,7 @@ } /*}}}*/ IssmDouble Numberedcostfunction::Response(FemModel* femmodel){/*{{{*/ - + _assert_(number_cost_functions>0 && number_cost_functions<1e3); /*output:*/ @@ -156,4 +155,3 @@ return value_sum; } /*}}}*/ - Index: ../trunk-jpl/src/c/classes/Params/BoolParam.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Params/BoolParam.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Params/BoolParam.cpp (revision 23066) @@ -62,4 +62,3 @@ } /*}}}*/ - Index: ../trunk-jpl/src/c/classes/Params/StringParam.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Params/StringParam.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Params/StringParam.cpp (revision 23066) @@ -54,7 +54,7 @@ int size = 0; if(marshall_direction==MARSHALLING_FORWARD || marshall_direction == MARSHALLING_SIZE)size=strlen(value)+1; - + MARSHALLING_ENUM(StringParamEnum); MARSHALLING(enum_type); MARSHALLING(size); Index: ../trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp (revision 23066) @@ -115,7 +115,7 @@ /*DoubleVecParam specific routines:*/ void DoubleVecParam::GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM){/*{{{*/ - + /*Assign output pointers:*/ if(pM) *pM=M; *pIssmDoublearray=values; Index: ../trunk-jpl/src/c/classes/Params/DoubleMatParam.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Params/DoubleMatParam.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Params/DoubleMatParam.cpp (revision 23066) @@ -115,7 +115,7 @@ /*DoubleMatParam specific routines:*/ void DoubleMatParam::GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM,int* pN){/*{{{*/ - + /*Assign output pointers:*/ if(pM) *pM=M; if(pN) *pN=N; Index: ../trunk-jpl/src/c/classes/Params/DataSetParam.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Params/DataSetParam.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Params/DataSetParam.cpp (revision 23066) @@ -54,7 +54,7 @@ void DataSetParam::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ if(marshall_direction==MARSHALLING_BACKWARD)value=new DataSet(); - + MARSHALLING_ENUM(DataSetParamEnum); MARSHALLING(enum_type); value->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); Index: ../trunk-jpl/src/c/classes/Params/IntParam.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Params/IntParam.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Params/IntParam.cpp (revision 23066) @@ -63,4 +63,3 @@ } /*}}}*/ - Index: ../trunk-jpl/src/c/classes/Params/StringArrayParam.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Params/StringArrayParam.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Params/StringArrayParam.cpp (revision 23066) @@ -82,7 +82,7 @@ } MARSHALLING_ENUM(StringArrayParamEnum); - + MARSHALLING(enum_type); MARSHALLING(numstrings); Index: ../trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.cpp (revision 23066) @@ -232,4 +232,3 @@ } /*}}}*/ - Index: ../trunk-jpl/src/c/classes/AmrBamg.cpp =================================================================== --- ../trunk-jpl/src/c/classes/AmrBamg.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/AmrBamg.cpp (revision 23066) @@ -33,7 +33,7 @@ this->deviatoricerror_threshold = -1; this->deviatoricerror_groupthreshold = -1; this->deviatoricerror_maximum = -1; - + /*Geometry and mesh as NULL*/ this->geometry = NULL; this->fathermesh = NULL; @@ -136,7 +136,7 @@ this->options->hmaxVertices = NULL; this->options->hmaxVerticesSize[0] = 0; this->options->hmaxVerticesSize[1] = 0; - + /*verify if the metric will be reseted or not*/ if(this->keepmetric==0){ if(this->options->metric) xDelete(this->options->metric); @@ -154,15 +154,15 @@ int *datalist = xNew(2); 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]); elementslist[3*i+1] = reCast(meshout->Triangles[4*i+1]); @@ -173,7 +173,7 @@ *pdatalist = datalist; *pxylist = xylist; *pelementslist = elementslist; - + /*Cleanup and return*/ delete geomout; }/*}}}*/ @@ -180,7 +180,7 @@ void AmrBamg::SetBamgOpts(IssmDouble hmin_in,IssmDouble hmax_in,IssmDouble err_in,IssmDouble gradation_in){/*{{{*/ if(!this->options) _error_("AmrBamg->options is NULL!"); - + this->options->hmin = hmin_in; this->options->hmax = hmax_in; this->options->err[0] = err_in; Index: ../trunk-jpl/src/c/classes/kriging/Observations.cpp =================================================================== --- ../trunk-jpl/src/c/classes/kriging/Observations.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/kriging/Observations.cpp (revision 23066) @@ -740,4 +740,3 @@ /*Assign output pointer*/ xDelete(counter); }/*}}}*/ - Index: ../trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp (revision 23066) @@ -22,7 +22,7 @@ #include "../classes/Inputs/Input.h" #include "../classes/gauss/Gauss.h" /*}}}*/ - + /*Cfsurfacelogvel constructors, destructors :*/ Cfsurfacelogvel::Cfsurfacelogvel(){/*{{{*/ @@ -39,13 +39,13 @@ Cfsurfacelogvel::Cfsurfacelogvel(char* in_name, int in_definitionenum, IssmDouble in_datatime, bool in_timepassedflag){/*{{{*/ this->definitionenum=in_definitionenum; - + this->name = xNew(strlen(in_name)+1); xMemCpy(this->name,in_name,strlen(in_name)+1); this->datatime=in_datatime; this->timepassedflag=in_timepassedflag; - + this->misfit=0; this->lock=0; } @@ -99,10 +99,10 @@ } /*}}}*/ IssmDouble Cfsurfacelogvel::Response(FemModel* femmodel){/*{{{*/ - + /*diverse: */ IssmDouble time; - + /*recover time parameters: */ femmodel->parameters->FindParam(&time,TimeEnum); if(time < last_time) timepassedflag = false; @@ -111,7 +111,7 @@ int i; IssmDouble J=0.; IssmDouble J_sum=0.; - + if(datatime<=time && !timepassedflag){ for(i=0;ielements->Size();i++){ Element* element=(Element*)femmodel->elements->GetObjectByOffset(i); @@ -121,7 +121,7 @@ ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); J=J_sum; - + timepassedflag = true; return J; } @@ -138,7 +138,7 @@ IssmDouble misfit,Jdet; IssmDouble vx,vy,vxobs,vyobs,weight; IssmDouble* xyz_list = NULL; - + /*Get basal element*/ if(!element->IsOnSurface()) return 0.; @@ -159,7 +159,7 @@ /* Get node coordinates*/ topelement->GetVerticesCoordinates(&xyz_list); - + /*Get model values*/ Input* vx_input =topelement->GetInput(VxEnum); _assert_(vx_input); Input* vy_input =NULL; @@ -174,7 +174,6 @@ if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do"); datasetinput = (DatasetInput*)tempinput; - /* Start looping on the number of gaussian points: */ Gauss* gauss=topelement->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ @@ -221,4 +220,3 @@ delete gauss; return Jelem; }/*}}}*/ - Index: ../trunk-jpl/src/c/classes/Cfdragcoeffabsgrad.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Cfdragcoeffabsgrad.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Cfdragcoeffabsgrad.cpp (revision 23066) @@ -22,7 +22,7 @@ #include "../classes/Inputs/Input.h" #include "../classes/gauss/Gauss.h" /*}}}*/ - + /*Cfdragcoeffabsgrad constructors, destructors :*/ Cfdragcoeffabsgrad::Cfdragcoeffabsgrad(){/*{{{*/ @@ -39,13 +39,13 @@ Cfdragcoeffabsgrad::Cfdragcoeffabsgrad(char* in_name, int in_definitionenum, int in_weights_enum, bool in_timepassedflag){/*{{{*/ this->definitionenum=in_definitionenum; - + this->name = xNew(strlen(in_name)+1); xMemCpy(this->name,in_name,strlen(in_name)+1); this->weights_enum=in_weights_enum; this->timepassedflag=in_timepassedflag; - + this->misfit=0; this->lock=0; } @@ -101,7 +101,7 @@ IssmDouble Cfdragcoeffabsgrad::Response(FemModel* femmodel){/*{{{*/ /*diverse: */ IssmDouble time; - + /*recover time parameters: */ int i; IssmDouble J=0.; @@ -115,7 +115,7 @@ ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); J=J_sum; - + timepassedflag = true; return J; } @@ -182,4 +182,3 @@ delete gauss; return Jelem; }/*}}}*/ - Index: ../trunk-jpl/src/c/classes/Loads/Neumannflux.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Neumannflux.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Loads/Neumannflux.cpp (revision 23066) @@ -30,7 +30,6 @@ /*}}}*/ Neumannflux::Neumannflux(int neumannflux_id,int i,IoModel* iomodel,int* segments,int in_analysis_type){/*{{{*/ - /*Some sanity checks*/ _assert_(segments); Index: ../trunk-jpl/src/c/classes/Loads/Moulin.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23066) @@ -177,7 +177,7 @@ this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); switch(analysis_type){ - + case HydrologyShaktiAnalysisEnum: pe = this->CreatePVectorHydrologyShakti(); break; @@ -394,10 +394,10 @@ if(!this->node->IsActive()) return NULL; 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*/ return pe; Index: ../trunk-jpl/src/c/classes/Cfsurfacesquare.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Cfsurfacesquare.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Cfsurfacesquare.cpp (revision 23066) @@ -22,7 +22,7 @@ #include "../classes/Inputs/Input.h" #include "../classes/gauss/Gauss.h" /*}}}*/ - + /*Cfsurfacesquare constructors, destructors :*/ Cfsurfacesquare::Cfsurfacesquare(){/*{{{*/ @@ -42,7 +42,7 @@ Cfsurfacesquare::Cfsurfacesquare(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, int in_weights_enum, IssmDouble in_datatime, bool in_timepassedflag){/*{{{*/ this->definitionenum=in_definitionenum; - + this->name = xNew(strlen(in_name)+1); xMemCpy(this->name,in_name,strlen(in_name)+1); @@ -51,7 +51,7 @@ this->weights_enum=in_weights_enum; this->datatime=in_datatime; this->timepassedflag=in_timepassedflag; - + this->misfit=0; this->lock=0; } @@ -110,7 +110,7 @@ IssmDouble Cfsurfacesquare::Response(FemModel* femmodel){/*{{{*/ /*diverse: */ IssmDouble time; - + /*recover time parameters: */ femmodel->parameters->FindParam(&time,TimeEnum); if(time < last_time) timepassedflag = false; @@ -129,7 +129,7 @@ ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); J=J_sum; - + timepassedflag = true; return J; } @@ -171,7 +171,7 @@ /*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"); datasetinput = (DatasetInput*)tempinput; @@ -213,4 +213,3 @@ delete gauss; return Jelem; }/*}}}*/ - Index: ../trunk-jpl/src/c/classes/Profiler.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Profiler.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Profiler.cpp (revision 23066) @@ -99,7 +99,6 @@ _assert_(tagrunning[tag]) _error_("Tag "<time[tag]; #else @@ -143,7 +142,6 @@ _assert_(tagrunning[tag]) _error_("Tag "<running[tag]) _error_("Tag "<fathermesh || !this->previousmesh) _error_("Impossible to execute refinement: fathermesh or previousmesh is NULL!\n"); @@ -117,10 +117,10 @@ /*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); @@ -129,7 +129,7 @@ } /*}}}*/ void AdaptiveMeshRefinement::RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror){/*{{{*/ - + /*Intermediaries*/ int nelem =-1; int side2D = 6; @@ -203,7 +203,7 @@ if(criteria) index.push_back(gmesh->Element(this->specialelementsindex[i])->FatherIndex()); } /*}}}*/ - + /*Now, detele the special elements{{{*/ if(this->refinement_type==1) this->DeleteSpecialElements(verbose,gmesh); else this->specialelementsindex.clear(); @@ -215,12 +215,12 @@ gmesh=this->fathermesh; } /*}}}*/ - + /*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;iElement(i)) continue; @@ -265,7 +265,7 @@ /*Adjust the connectivities before continue*/ gmesh->BuildConnectivity(); /*}}}*/ - + /*Refinement process: loop over level of refinements{{{*/ if(verbose) _printf_("\trefinement process (level max = "<level_max<<")\n"); if(verbose) _printf_("\ttotal: "); @@ -319,14 +319,14 @@ * 2 : uniform refinment * */ 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++){ sons.clear(); @@ -334,11 +334,11 @@ if(sons.size()) count++; //if neighbour was refined 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; } /*}}}*/ @@ -390,14 +390,14 @@ } /*}}}*/ 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 count=0; @@ -467,13 +467,13 @@ long* vertex_index2sid = xNew(ntotalvertices); 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; for(int i=0;iNElements();i++){//over gmesh elements index @@ -510,7 +510,7 @@ newmeshXY[2*sid] = coords[0]; // X 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++) { geoel = gmesh->ElementVec()[this->sid2index[i]]; @@ -532,7 +532,7 @@ xc=newmeshXY[2*c]; yc=newmeshXY[2*c+1]; detJ=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya); - + /*verify and swap, if necessary*/ if(detJ<0) { newelements[i*this->GetNumberOfNodes()+0]=b+1;//a->b @@ -544,7 +544,7 @@ *pdata = newdata; *pxy = newmeshXY; *pelements = newelements; - + /*Cleanup*/ xDelete(vertex_index2sid); } @@ -553,7 +553,7 @@ /* 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"); if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to create initial mesh: refinement type is not defined!\n"); @@ -560,7 +560,7 @@ /*Verify and creating initial mesh*/ if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!"); - + this->fathermesh = new TPZGeoMesh(); this->fathermesh->NodeVec().Resize(nvertices); @@ -575,7 +575,7 @@ this->fathermesh->NodeVec()[i].SetCoord(coord); this->fathermesh->NodeVec()[i].SetNodeId(i); } - + /*Generate the elements*/ long index; const int mat = this->GetElemMaterialID(); @@ -603,10 +603,10 @@ } /*}}}*/ TPZGeoMesh* AdaptiveMeshRefinement::CreateRefPatternMesh(TPZGeoMesh* gmesh){/*{{{*/ - + TPZGeoMesh *newgmesh = new TPZGeoMesh(); newgmesh->CleanUp(); - + int nnodes = gmesh->NNodes(); int nelem = gmesh->NElements(); int mat = this->GetElemMaterialID();; @@ -616,25 +616,25 @@ //nodes 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(); newgmesh->ElementVec()[index] = NULL; 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; TPZVec< std::vector > neighbourhood(nsides); @@ -653,7 +653,7 @@ neighS = neighS.Neighbour(); } } - + //inserting in new element for(int s = 0; s < nsides; s++){ TPZGeoEl * tempEl = newgeoel; @@ -667,17 +667,17 @@ } 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); if(!son){ @@ -689,7 +689,7 @@ /*Now, build connectivities*/ newgmesh->BuildConnectivity(); - + return newgmesh; } /*}}}*/ @@ -704,7 +704,7 @@ } /*}}}*/ void AdaptiveMeshRefinement::PrintGMeshVTK(TPZGeoMesh* gmesh,std::ofstream &file,bool matColor){/*{{{*/ - + file.clear(); long nelements = gmesh->NElements(); TPZGeoEl *gel; @@ -730,10 +730,10 @@ } MElementType elt = gel->Type(); int elNnodes = MElementType_NNodes(elt); - + size += (1+elNnodes); connectivity << elNnodes; - + for(int t = 0; t < elNnodes; t++) { for(int c = 0; c < 3; c++) @@ -742,15 +742,15 @@ node << coord << " "; } node << std::endl; - + actualNode++; connectivity << " " << actualNode; } connectivity << std::endl; - + int elType = this->GetVTK_ElType(gel); type << elType << std::endl; - + if(matColor == true) { material << gel->MaterialId() << std::endl; @@ -759,7 +759,7 @@ { material << gel->Index() << std::endl; } - + nVALIDelements++; } node << std::endl; @@ -789,9 +789,9 @@ } /*}}}*/ int AdaptiveMeshRefinement::GetVTK_ElType(TPZGeoEl * gel){/*{{{*/ - + MElementType pzElType = gel->Type(); - + int elType = -1; switch (pzElType) { @@ -848,8 +848,7 @@ std::cout << "MIGHT BE CURVED ELEMENT (quadratic or quarter point)" << std::endl; DebugStop(); } - + return elType; } /*}}}*/ - Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23066) @@ -105,7 +105,7 @@ for(i=0;inodes[i]) tria->nodes[i]=this->nodes[i]; else tria->nodes[i] = NULL; } else tria->nodes = NULL; - + tria->vertices = (Vertex**)this->hvertices->deliverp(); tria->material = (Material*)this->hmaterial->delivers(); tria->matpar = (Matpar*)this->hmatpar->delivers(); @@ -114,9 +114,9 @@ } /*}}}*/ 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); Element::MarshallElement(pmarshalled_data,pmarshalled_data_size,marshall_direction,this->numanalyses); @@ -134,7 +134,7 @@ /*Call inputs method*/ _assert_(this->inputs); - + int domaintype; parameters->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ @@ -184,7 +184,7 @@ if(N!=num_inputs) _error_("sizes are not consistent"); 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 } @@ -327,7 +327,7 @@ } /*}}}*/ void Tria::CalvingCrevasseDepth(){/*{{{*/ - + IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble calvingrate[NUMVERTICES]; IssmDouble vx,vy,vel; @@ -339,10 +339,10 @@ /* 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); IssmDouble rho_freshwater = this->GetMaterialParameter(MaterialsRhoFreshwaterEnum); @@ -360,12 +360,12 @@ Input* s_xx_input = inputs->GetInput(DeviatoricStressxxEnum); _assert_(s_xx_input); 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); surface_input->GetInputValue(&float_depth,gauss); @@ -377,13 +377,13 @@ s_xx_input->GetInputValue(&s_xx,gauss); s_xy_input->GetInputValue(&s_xy,gauss); s_yy_input->GetInputValue(&s_yy,gauss); - + vel=sqrt(vx*vx+vy*vy)+1.e-14; s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2)); 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.; @@ -398,19 +398,19 @@ //if (surface_crevasse[iv]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)); this->inputs->AddInput(new TriaInput(CrevasseDepthEnum,&crevasse_depth[0],P1Enum)); @@ -460,7 +460,7 @@ } else calvingrate[iv]=0.; - + calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14); calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14); } @@ -561,7 +561,7 @@ IssmDouble strain_xy[NUMVERTICES]; IssmDouble vorticity_xy[NUMVERTICES]; GaussTria* gauss=NULL; - + /* Get node coordinates and dof list: */ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -568,7 +568,7 @@ /*Retrieve all inputs we will be needing: */ 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(); for (int iv=0;ivinputs->GetInput(control_enum); _assert_(input); 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; } @@ -1699,7 +1698,7 @@ _error_("Data " << data << " not supported yet"); } /*Check what input we are dealing with*/ - + switch(input->ObjectEnum()){ case TriaInputEnum: { @@ -1716,7 +1715,7 @@ vector->SetValues(NUMVERTICES,idlist,values,INS_VAL); break; } - + case TransientInputEnum: { TransientInput* transientinput = xDynamicCast(input); @@ -1981,7 +1980,7 @@ surface_input->GetInputAverage(&surface); base_input->GetInputAverage(&bed); bed_input->GetInputAverage(&bathymetry); - + /*Return: */ return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.)); } @@ -2026,8 +2025,6 @@ 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*/ for(i=0;i<3;i++){ tria_vertex_ids[i]=reCast(iomodel->elements[3*index+i]); //ids for vertices are in the elements array from Matlab @@ -2393,7 +2390,6 @@ /*}}}*/ IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/ - /*intermediary: */ IssmDouble* values=NULL; Input* thickness_input=NULL; @@ -2406,7 +2402,7 @@ int point1; IssmDouble fraction1,fraction2; bool mainlynegative=true; - + /*Output:*/ volume=0; @@ -2415,7 +2411,7 @@ /*Retrieve inputs required:*/ thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); - + /*Retrieve material parameters: */ rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); @@ -2424,7 +2420,7 @@ for(int i=0;ivertices[i]->Sid()]; } - + /*Ok, use the level set values to figure out where we put our gaussian points:*/ this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values); Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4); @@ -2876,7 +2872,7 @@ dist_gl_input->GetInputValue(&dist_gl,gauss); dist_cf_input->GetInputValue(&dist_cf,gauss); delete gauss; - + /*Ensure values are positive for floating ice*/ dist_gl = fabs(dist_gl); dist_cf = fabs(dist_cf); @@ -3196,7 +3192,6 @@ } } - } /*}}}*/ void Tria::ResetFSBasalBoundaryCondition(void){/*{{{*/ @@ -3277,7 +3272,6 @@ IssmDouble values[NUMVERTICES]; int idlist[NUMVERTICES],control_init; - /*Get Domain type*/ int domaintype; parameters->FindParam(&domaintype,DomainTypeEnum); @@ -3297,7 +3291,7 @@ /*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){ _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput"); @@ -3330,7 +3324,6 @@ IssmDouble values[NUMVERTICES]; int vertexpidlist[NUMVERTICES],control_init; - /*Get Domain type*/ int domaintype; parameters->FindParam(&domaintype,DomainTypeEnum); @@ -4166,7 +4159,7 @@ IssmDouble* U_elastic_precomputed = NULL; IssmDouble* H_elastic_precomputed = NULL; int M, hemi; - + /*computation of Green functions:*/ IssmDouble* U_elastic= NULL; IssmDouble* N_elastic= NULL; @@ -4173,7 +4166,7 @@ IssmDouble* E_elastic= NULL; IssmDouble* X_elastic= NULL; IssmDouble* Y_elastic= NULL; - + /*optimization:*/ bool store_green_functions=false; @@ -4181,7 +4174,7 @@ Input* deltathickness_input=inputs->GetInput(EsaDeltathicknessEnum); 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; @@ -4194,7 +4187,7 @@ /*which hemisphere? for north-south, east-west components*/ this->parameters->FindParam(&hemi,EsaHemisphereEnum); - + /*compute area of element:*/ area=GetArea(); @@ -4253,7 +4246,7 @@ U_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*U_elastic[i]; 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) { ang2 = PI/2 - atan2(yy[i],xx[i]); @@ -4276,7 +4269,7 @@ pEast->SetValues(gsize,indices,E_values,ADD_VAL); pX->SetValues(gsize,indices,X_values,ADD_VAL); pY->SetValues(gsize,indices,Y_values,ADD_VAL); - + /*free ressources:*/ xDelete(indices); xDelete(U_values); xDelete(N_values); xDelete(E_values); @@ -4306,12 +4299,12 @@ IssmDouble* U_elastic_precomputed = NULL; 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; @@ -4319,7 +4312,7 @@ Input* deltathickness_input=inputs->GetInput(EsaDeltathicknessEnum); 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; @@ -4418,7 +4411,7 @@ dx = x_element-x; dy = y_element-y; dz = z_element-z; 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)); U_elastic[i] += U_elastic_precomputed[index]; @@ -4433,7 +4426,7 @@ pUp->SetValues(gsize,indices,U_values,ADD_VAL); pNorth->SetValues(gsize,indices,N_values,ADD_VAL); pEast->SetValues(gsize,indices,E_values,ADD_VAL); - + /*free ressources:*/ xDelete(indices); xDelete(U_values); xDelete(N_values); xDelete(E_values); @@ -4454,7 +4447,7 @@ IssmDouble Tria::OceanAverage(IssmDouble* Sg){ /*{{{*/ if(IsWaterInElement()){ - + IssmDouble area; /*Compute area of element:*/ @@ -4527,13 +4520,13 @@ 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) * also consistent with my GMD formulation! @@ -4545,20 +4538,20 @@ } 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; }/*}}}*/ void Tria::SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ @@ -4591,7 +4584,7 @@ bool computerigid = true; bool computeelastic= true; bool scaleoceanarea= false; - + /*early return if we are not on an ice cap:*/ if(!(this->inputs->Max(MaskIceLevelsetEnum)<=0)){ constant=0; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum)); @@ -4605,12 +4598,12 @@ *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! 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)); @@ -4635,10 +4628,10 @@ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); /* Where is the centroid of this element?:{{{*/ - + /*retrieve coordinates: */ ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical); - + IssmDouble minlong=400; IssmDouble maxlong=-20; for (int i=0;iGetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias area*=phi; } - + /*Compute ice thickness change: */ Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); @@ -4695,7 +4688,7 @@ bool mainlyfloating = true; int point1; IssmDouble fraction1,fraction2; - + /*Recover portion of element that is grounded*/ this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2); @@ -4749,12 +4742,12 @@ values[i]=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic); } pSgi->SetValues(gsize,indices,values,ADD_VAL); - + /*free ressources:*/ xDelete(values); xDelete(indices); } - + /*Assign output pointer:*/ _assert_(!xIsNan(eustatic)); _assert_(!xIsInf(eustatic)); @@ -4779,7 +4772,7 @@ /*precomputed elastic green functions:*/ IssmDouble* G_elastic_precomputed = NULL; int M; - + /*computation of Green functions:*/ IssmDouble* G_elastic= NULL; IssmDouble* G_rigid= NULL; @@ -4856,9 +4849,9 @@ late=late/180*PI; longe=longe/180*PI; /*}}}*/ - + if(computeelastic){ - + /*recover elastic green function:*/ DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M); @@ -4896,7 +4889,7 @@ values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i]; } } - + pSgo->SetValues(gsize,indices,values,ADD_VAL); /*free ressources:*/ @@ -4929,7 +4922,7 @@ IssmDouble* U_elastic_precomputed = NULL; IssmDouble* H_elastic_precomputed = NULL; int M; - + /*computation of Green functions:*/ IssmDouble* U_elastic= NULL; IssmDouble* N_elastic= NULL; @@ -4956,7 +4949,7 @@ /*recover computational flags: */ this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); - + /*early return if elastic not requested:*/ if(!computeelastic) return; @@ -5025,12 +5018,12 @@ /*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); if(horiz){ @@ -5072,7 +5065,7 @@ 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)); U_elastic[i] += U_elastic_precomputed[index]; Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23066) @@ -286,7 +286,6 @@ IssmDouble calvingratey[NUMVERTICES]; IssmDouble calvingrate[NUMVERTICES]; - /* Get node coordinates and dof list: */ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -1047,7 +1046,7 @@ } /*}}}*/ void Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/ - + /* Intermediaries */ const int dim=3; int i, dir,nrfrontnodes; @@ -1181,7 +1180,6 @@ if(!IsInputEnum(control_enum)) _error_("Enum "<inputs->GetInput(control_enum); _assert_(input); - /*Cast to Controlinput*/ if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput"); ControlInput* controlinput = xDynamicCast(input); @@ -2594,7 +2592,7 @@ Penta* penta=this; for(;;){ - + IssmDouble xyz_list[NUMVERTICES][3]; /* Get node coordinates and dof list: */ ::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES); @@ -2603,7 +2601,7 @@ Jdet[0]=(xyz_list[3][2]-xyz_list[0][2])*0.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); Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); @@ -2614,7 +2612,7 @@ Input* deviayy_input=inputs->GetInput(DeviatoricStressyyEnum); _assert_(deviayy_input); 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++){ GaussPenta* gauss=new GaussPenta(ig,3+ig,11); @@ -2643,10 +2641,10 @@ } delete gauss; } - + /*Stop if we have reached the surface/base*/ if(penta->IsOnSurface()) break; - + /*get upper Penta*/ penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id); Index: ../trunk-jpl/src/c/classes/Elements/ElementHook.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/ElementHook.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Elements/ElementHook.cpp (revision 23066) @@ -103,7 +103,7 @@ MARSHALLING_DYNAMIC(hnodesi_null,bool,numanalyses); if(marshall_direction==MARSHALLING_BACKWARD){ - + if (!hnodes_null)this->hnodes = new Hook*[numanalyses]; else this->hnodes=NULL; this->hvertices = new Hook(); @@ -139,7 +139,7 @@ _printf_(" ElementHook DeepEcho:\n"); _printf_(" numanalyses : "<< this->numanalyses <<"\n"); - + _printf_(" hnodes:\n"); if(hnodes){ for(int i=0;inumanalyses;i++) { @@ -148,11 +148,11 @@ } } else _printf_(" hnodes = NULL\n"); - + _printf_(" hvertices:\n"); if(hvertices) hvertices->DeepEcho(); else _printf_(" hvertices = NULL\n"); - + _printf_(" hmaterial:\n"); if(hmaterial) hmaterial->DeepEcho(); else _printf_(" hmaterial = NULL\n"); @@ -169,10 +169,10 @@ } /*}}}*/ void ElementHook::Echo(){/*{{{*/ - + _printf_(" ElementHook Echo:\n"); _printf_(" numanalyses : "<< this->numanalyses <<"\n"); - + _printf_(" hnodes:\n"); if(hnodes){ for(int i=0;inumanalyses;i++) { @@ -180,11 +180,11 @@ } } else _printf_(" hnodes = NULL\n"); - + _printf_(" hvertices:\n"); if(hvertices) hvertices->Echo(); else _printf_(" hvertices = NULL\n"); - + _printf_(" hmaterial:\n"); if(hmaterial) hmaterial->Echo(); else _printf_(" hmaterial = NULL\n"); Index: ../trunk-jpl/src/c/classes/Elements/Seg.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 23066) @@ -103,7 +103,6 @@ return seg; - } /*}}}*/ void Seg::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ @@ -140,7 +139,7 @@ } /*}}}*/ void Seg::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/ - + /* Intermediaries */ int nrfrontnodes,index; IssmDouble levelset[NUMVERTICES]; Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23066) @@ -1104,7 +1104,6 @@ /*}}}*/ void Element::GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug){/*{{{*/ - /*Get number of nodes for this element*/ int numnodes = this->GetNumberOfNodes(); @@ -1121,7 +1120,6 @@ 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*/ for(int i=0;iGetDof(0,GsetEnum)]>input_min) min[nodes[i]->GetDof(0,GsetEnum)] = input_min; @@ -1388,7 +1386,7 @@ default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); } - + /*Clean up*/ xDelete(doflist); xDelete(values); @@ -1476,7 +1474,6 @@ parameters->FindParam(&num_controls,InversionNumControlParametersEnum); parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); - /*Get number of vertices*/ int numvertices = this->GetNumberOfVertices(); if(isautodiff){ @@ -1807,7 +1804,6 @@ this->inputs->AddInput(datasetinput); } - /*Branch on type of vector: nodal or elementary: */ if(vector_type==1){ //nodal vector @@ -2066,7 +2062,7 @@ this->GetInputListOnVertices(deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum); this->GetInputListOnVertices(deepwaterel,BasalforcingsDeepwaterElevationEnum); this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum); - + for(int i=0;iupperwaterel[i]) values[i]=0; else if (base[i]GetMaterialParameter(MaterialsRhoIceEnum); rho_water = matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum); @@ -3071,7 +3066,6 @@ /*Snow, firn and ice albedo:*/ 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: ../trunk-jpl/src/c/classes/gauss/GaussSeg.cpp =================================================================== --- ../trunk-jpl/src/c/classes/gauss/GaussSeg.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/gauss/GaussSeg.cpp (revision 23066) @@ -30,7 +30,7 @@ /*Get gauss points*/ this->numgauss = order; GaussLegendreLinear(&pcoords1,&pweights,order); - + this->coords1=xNew(numgauss); this->weights=xNew(numgauss); Index: ../trunk-jpl/src/c/classes/Regionaloutput.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Regionaloutput.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Regionaloutput.cpp (revision 23066) @@ -154,4 +154,3 @@ return val_t; } /*}}}*/ - Index: ../trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 23066) @@ -19,7 +19,7 @@ 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){ femmodel_init= new FemModel(argc,argv,evaluation_comm); @@ -32,7 +32,7 @@ int world_rank; ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&world_rank); - + if(world_rank!=0){ /*Wrap up: */ Index: ../trunk-jpl/src/c/classes/Node.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Node.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Node.cpp (revision 23066) @@ -497,7 +497,7 @@ } /*}}}*/ void Node::SetApproximation(int in_approximation){/*{{{*/ - + this->approximation = in_approximation; } /*}}}*/ Index: ../trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp (revision 23066) @@ -36,7 +36,7 @@ /*Object virtual functions definitions:*/ Object* SpcStatic::copy() {/*{{{*/ - + SpcStatic* spcstat = new SpcStatic(*this); spcstat->id=this->id; Index: ../trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/IoModel.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/IoModel.cpp (revision 23066) @@ -156,7 +156,7 @@ bool autodiff=false; bool iscontrol=false; - + /*First, keep track of the file handle: */ this->fid=iomodel_handle; @@ -421,7 +421,7 @@ /*Initialize array detecting whether data[i] is an independent AD mode variable: */ this->FetchData(&autodiff,"md.autodiff.isautodiff"); this->FetchData(&iscontrol,"md.inversion.iscontrol"); - + if(trace || (autodiff && !iscontrol)){ #ifdef _HAVE_ADOLC_ @@ -505,7 +505,7 @@ void IoModel::DeleteData(char*** pstringarray, int numstrings,const char* data_name){/*{{{*/ 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); if(strcmp(record_name,"md.basalforcings.model")==0) integer = IoCodeToEnumBasal(integer); @@ -962,7 +962,7 @@ /*recover my_rank:*/ int my_rank=IssmComm::GetRank(); - + /*Set file pointer to beginning of the data: */ fid=this->SetFilePointerToData(&code,NULL,data_name); @@ -1133,7 +1133,7 @@ return; } } - + /*output: */ int M,N; IssmPDouble *matrix = NULL; @@ -1824,7 +1824,7 @@ if(my_rank==0){ /*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; *pnumstrings=num_instances; @@ -1874,7 +1874,7 @@ /*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); @@ -1889,7 +1889,7 @@ code=codes[i]; if(code!=2)_error_("expecting an integer for \""<(file_positions); xDelete(codes); @@ -1927,7 +1927,7 @@ /*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); @@ -1942,7 +1942,7 @@ code=codes[i]; if(code!=3)_error_("expecting a double for \""<(file_positions); xDelete(codes); @@ -1984,7 +1984,7 @@ /*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,7 +2014,6 @@ } ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - /*Now allocate matrix: */ if(M*N){ pmatrix=xNew(M*N); @@ -2038,8 +2037,7 @@ } else matrix=NULL; - - + /*Assign: */ mdims[i]=M; matrices[i]=matrix; @@ -2046,7 +2044,7 @@ ndims[i]=N; } } - + /*Free ressources:*/ xDelete(file_positions); xDelete(codes); @@ -2088,7 +2086,7 @@ /*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,7 +2116,6 @@ } ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - /*Now allocate matrix: */ if(M*N){ pmatrix=xNew(M*N); @@ -2143,8 +2140,7 @@ } else integer_matrix=NULL; - - + /*Assign: */ mdims[i]=M; matrices[i]=integer_matrix; @@ -2151,7 +2147,7 @@ ndims[i]=N; } } - + /*Free ressources:*/ xDelete(file_positions); xDelete(codes); @@ -2375,7 +2371,7 @@ codes = xNew(num_instances); 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: */ counter=0; @@ -2418,12 +2414,12 @@ codes[counter] = record_code; 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: ../trunk-jpl/src/c/classes/Nodalvalue.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Nodalvalue.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Nodalvalue.cpp (revision 23066) @@ -86,7 +86,7 @@ } /*}}}*/ IssmDouble Nodalvalue::Response(FemModel* femmodel){/*{{{*/ - + /*output:*/ IssmDouble value; Index: ../trunk-jpl/src/c/classes/Misfit.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Misfit.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Misfit.cpp (revision 23066) @@ -22,7 +22,7 @@ #include "../classes/Inputs/Input.h" #include "../classes/gauss/Gauss.h" /*}}}*/ - + /*Misfit constructors, destructors :*/ Misfit::Misfit(){/*{{{*/ @@ -41,18 +41,18 @@ Misfit::Misfit(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, char* in_timeinterpolation, int in_local, int in_weights_enum){/*{{{*/ this->definitionenum=in_definitionenum; - + this->name = xNew(strlen(in_name)+1); xMemCpy(this->name,in_name,strlen(in_name)+1); 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; } @@ -110,11 +110,11 @@ } /*}}}*/ IssmDouble Misfit::Response(FemModel* femmodel){/*{{{*/ - + /*diverse: */ IssmDouble time,starttime,finaltime; IssmDouble dt; - + /*recover time parameters: */ femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum); @@ -129,7 +129,6 @@ IssmDouble area_t=0.; IssmDouble all_area_t; - /*If we are locked, return time average: */ if(this->lock)return misfit/(time-starttime); @@ -143,7 +142,7 @@ ISSM_MPI_Allreduce ( (void*)&area_t,(void*)&all_area_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); 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; @@ -163,11 +162,11 @@ IssmDouble* observation= NULL; 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); GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize); @@ -189,7 +188,7 @@ return misfit; } else{ - + IssmDouble misfit_t=0.; IssmDouble all_misfit_t=0.; @@ -200,7 +199,7 @@ GetVectorFromInputsx(&model,&msize,femmodel,model_enum); 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); xDelete(observation); @@ -226,9 +225,9 @@ } /*}}}*/ else{ /*global computation: {{{ */ - + IssmDouble model, observation; - + /*If we are locked, return time average: */ if(this->lock)return misfit/(time-starttime); @@ -238,13 +237,13 @@ Element* element=(Element*)femmodel->elements->GetObjectByOffset(0); _assert_(element); 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: ../trunk-jpl/src/c/datastructures/DataSet.cpp =================================================================== --- ../trunk-jpl/src/c/datastructures/DataSet.cpp (revision 23065) +++ ../trunk-jpl/src/c/datastructures/DataSet.cpp (revision 23066) @@ -91,7 +91,7 @@ /*Specific methods*/ void DataSet::Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ - + vector::iterator obj; int obj_size=0; int obj_enum=0; Index: ../trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 23066) @@ -593,7 +593,6 @@ xDelete(responses); delete gauss; - }/*}}}*/ void BalancethicknessAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 23066) @@ -260,7 +260,6 @@ basis,1,numnodes,0, &Ke->values[0],1); - /*Transfer EPL part*/ transfer=GetHydrologyKMatrixTransfer(basalelement); D_scalar=dt*transfer*gauss->weight*Jdet; Index: ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 23066) @@ -1152,7 +1152,7 @@ /*Clean up and return*/ xDelete(responses); - + }/*}}}*/ void AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector* gradient,int control_index){/*{{{*/ /*WARNING: We use SSA as an estimate for now*/ @@ -2111,7 +2111,6 @@ IssmDouble vx,vy,lambda,mu; IssmDouble *xyz_list= NULL; - /*Fetch number of vertices for this finite element*/ int numvertices = basalelement->GetNumberOfVertices(); @@ -2131,8 +2130,6 @@ Input* adjointx_input = basalelement->GetInput(AdjointxEnum); _assert_(adjointx_input); Input* adjointy_input = basalelement->GetInput(AdjointyEnum); _assert_(adjointy_input); - - IssmDouble q_exp; IssmDouble C_param; IssmDouble As; @@ -2181,11 +2178,11 @@ vmag = sqrt(vx*vx + vy*vy); 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;iweight*basis[i]; @@ -2318,7 +2315,6 @@ default: _error_("mesh "<NumberofNodesVelocity(); int pnumnodes = element->NumberofNodesPressure(); Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 23066) @@ -42,7 +42,7 @@ counter++; } } - + iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum); iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum); @@ -176,7 +176,7 @@ } /*Calving threshold*/ - + /*Fetch number of nodes and dof for this finite element*/ int numnodes = basalelement->GetNumberOfNodes(); @@ -325,7 +325,7 @@ c[i]=0.; m[i]=0.; } break; - + case CalvingLevermannEnum: calvingratex_input->GetInputValue(&c[0],gauss); if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss); @@ -356,7 +356,7 @@ m[i]=0.; } break; - + case CalvingHabEnum: lsf_slopex_input->GetInputValue(&dlsf[0],gauss); if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); @@ -377,7 +377,7 @@ m[i]=0.; } break; - + case CalvingCrevasseDepthEnum: lsf_slopex_input->GetInputValue(&dlsf[0],gauss); if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); @@ -384,7 +384,7 @@ meltingrate_input->GetInputValue(&meltingrate,gauss); if(groundedice<0) meltingrate = 0.; - + norm_dlsf=0.; for(i=0;iIsOnBase()) return NULL; Element* basalelement = element->SpawnBasalElement(); @@ -524,7 +524,7 @@ IssmDouble Jdet,dt; IssmDouble lsf; IssmDouble* xyz_list = NULL; - + /*Fetch number of nodes and dof for this finite element*/ int numnodes = basalelement->GetNumberOfNodes(); @@ -531,7 +531,7 @@ /*Initialize Element vector*/ ElementVector* pe = basalelement->NewElementVector(); basalelement->FindParam(&dt,TimesteppingTimeStepEnum); - + if(dt!=0.){ /*Initialize basis vector*/ IssmDouble* basis = xNew(numnodes); @@ -622,7 +622,7 @@ // returns distance d of point q to straight going through points s0, s1 // d=|a x b|/|b| // with a=q-s0, b=s1-s0 - + /* Intermediaries */ const int dim=2; int i; @@ -633,7 +633,7 @@ a[i]=q[i]-s0[i]; b[i]=s1[i]-s0[i]; } - + norm_b=0.; for(i=0;iparameters->FindParam(&calvinglaw,CalvingLawEnum); if(calvinglaw==CalvingMinthicknessEnum){ @@ -702,17 +702,17 @@ delete gauss; } } - + if(calvinglaw==CalvingHabEnum){ /*Get the fraction of the flotation thickness at the terminus*/ 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); @@ -743,13 +743,13 @@ delete gauss; } } - + if(calvinglaw==CalvingCrevasseDepthEnum){ - + int nflipped,local_nflipped; Vector* vec_constraint_nodes = NULL; IssmDouble* constraint_nodes = NULL; - + /*Get the DistanceToCalvingfront*/ femmodel->elements->InputDuplicate(MaskIceLevelsetEnum,DistanceToCalvingfrontEnum); femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum); @@ -841,7 +841,7 @@ for(int in=0;inGaussNode(element->GetElementType(),in); Node* node=element->GetNode(in); - + if(constraint_nodes[node->Sid()]>0.){ node->ApplyConstraint(0,+1.); } Index: ../trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 23066) @@ -29,7 +29,6 @@ else iscoupling = false; - /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/ if(!iscoupling){ IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvz",StressbalanceVerticalAnalysisEnum,P1Enum,0); @@ -158,7 +157,6 @@ }/*}}}*/ ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixBase(Element* element){/*{{{*/ - if(!element->IsOnBase()) return NULL; /*Intermediaries*/ @@ -199,7 +197,6 @@ }/*}}}*/ ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/ - if(!element->IsOnSurface()) return NULL; /*Intermediaries*/ Index: ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 23066) @@ -294,7 +294,7 @@ B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); A=pow(B,-n); - + /*Compute beta term*/ if(gap0.); - + for(int i=0;ivalues[i]+=Jdet*gauss->weight* ( meltrate*(1/rho_water-1/rho_ice) @@ -397,7 +397,7 @@ /*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"); } @@ -443,7 +443,7 @@ Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); reynolds_input->GetInputAverage(&reynolds); gap_input->GetInputAverage(&gap); - + /*Compute conductivity*/ IssmDouble conductivity = pow(gap,3)*g/(12.*NU*(1+OMEGA*reynolds)); _assert_(conductivity>0); @@ -453,7 +453,6 @@ }/*}}}*/ void HydrologyShaktiAnalysis::UpdateGapHeight(FemModel* femmodel){/*{{{*/ - for(int j=0;jelements->Size();j++){ Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); UpdateGapHeight(element); @@ -546,7 +545,7 @@ IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.); 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]); dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]); @@ -580,7 +579,7 @@ /*Add new gap as an input*/ element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum); - + /*Divide by connectivity, add basal flux as an input*/ q = q/totalweights; element->AddInput(HydrologyBasalFluxEnum,&q,P0Enum); Index: ../trunk-jpl/src/c/analyses/EsaAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 23066) @@ -39,7 +39,7 @@ int nl; IssmDouble* love_h=NULL; IssmDouble* love_l=NULL; - + IssmDouble* U_elastic = NULL; IssmDouble* U_elastic_local = NULL; IssmDouble* H_elastic = NULL; @@ -68,7 +68,7 @@ M=reCast(180./degacc+1.); U_elastic=xNew(M); H_elastic=xNew(M); - + /*compute combined legendre + love number (elastic green function:*/ m=DetermineLocalSize(M,IssmComm::GetComm()); GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); @@ -95,7 +95,7 @@ IssmDouble deltalove_U; 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){ Pn=1; @@ -198,7 +198,7 @@ void EsaAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ /*Default, do nothing*/ return; - + }/*}}}*/ void EsaAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Default, do nothing*/ Index: ../trunk-jpl/src/c/analyses/SmoothAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/SmoothAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/SmoothAnalysis.cpp (revision 23066) @@ -152,7 +152,6 @@ default: input = element->GetInput(input_enum); } - /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ @@ -161,7 +160,6 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); - switch(input_enum){ case DrivingStressXEnum: case DrivingStressYEnum:{ Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 23066) @@ -164,7 +164,7 @@ } iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum); iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum); - + /*Initialize cumdeltalthickness input*/ InputUpdateFromConstantx(elements,0.,SealevelriseCumDeltathicknessEnum); @@ -202,7 +202,7 @@ iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum); } - + }/*}}}*/ void MasstransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ @@ -220,8 +220,6 @@ parameters->AddObject(new IntParam(MasstransportNumRequestedOutputsEnum,numoutputs)); if(numoutputs)parameters->AddObject(new StringArrayParam(MasstransportRequestedOutputsEnum,requestedoutputs,numoutputs)); iomodel->DeleteData(&requestedoutputs,numoutputs,"md.masstransport.requested_outputs"); - - }/*}}}*/ @@ -804,7 +802,7 @@ xDelete(phi); xDelete(sealevel); xDelete(bed); - + xDelete(doflist); if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; }/*}}}*/ Index: ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 23066) @@ -68,7 +68,7 @@ IssmDouble* love_h=NULL; IssmDouble* love_k=NULL; IssmDouble* love_l=NULL; - + bool elastic=false; IssmDouble* G_elastic = NULL; IssmDouble* G_elastic_local = NULL; @@ -131,7 +131,6 @@ H_elastic=xNew(M); #endif - /*compute combined legendre + love number (elastic green function:*/ m=DetermineLocalSize(M,IssmComm::GetComm()); GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); @@ -167,7 +166,7 @@ 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){ Pn=1; @@ -229,7 +228,7 @@ xDelete(H_elastic); xDelete(H_elastic_local); } - + /*Transitions: */ iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions"); if(transitions){ @@ -249,7 +248,6 @@ if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs)); iomodel->DeleteData(&requestedoutputs,numoutputs,"md.slr.requested_outputs"); - }/*}}}*/ /*Finite Element Analysis*/ Index: ../trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp (revision 23066) @@ -198,8 +198,6 @@ xDelete(lset); } - - return; }/*}}}*/ Index: ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 23066) @@ -378,7 +378,6 @@ return; }/*}}}*/ - /*Needed changes to switch to the Johnson formulation*//*{{{*/ /*All the changes are to be done in the velocity computation. The new velocity needs some new parameter that should be introduce in the hydrologyshreve class: Index: ../trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp (revision 23066) @@ -199,7 +199,6 @@ Ny[i] = -H[i]*Ny[i]/norm; } - /* Start looping on the number of gaussian points: */ Gauss* gauss=basalelement->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 23066) @@ -213,7 +213,6 @@ return NULL; } - /*Intermediaries */ bool active_element,isefficientlayer; IssmDouble D_scalar,Jdet,dt; @@ -503,7 +502,6 @@ default: _error_("mesh "<GetNumberOfNodes(); Index: ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23066) @@ -18,10 +18,10 @@ return 1; }/*}}}*/ 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; for(int i=0;inumberofelements;i++){ @@ -31,10 +31,10 @@ counter++; } } - + /*Figure out smb model: */ iomodel->FindConstant(&smb_model,"md.smb.model"); - + switch(smb_model){ case SMBforcingEnum: iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.); @@ -141,8 +141,6 @@ default: _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: /*Nothing to add to parameters*/ @@ -197,7 +195,7 @@ iomodel->FetchData(&temp,&N,&M,"md.smb.Pfac"); _assert_(N==2); 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)); iomodel->DeleteData(temp,"md.smb.Tdiff"); @@ -265,7 +263,7 @@ /*Figure out smb model: */ femmodel->parameters->FindParam(&smb_model,SmbEnum); - + /*branch to correct module*/ switch(smb_model){ case SMBforcingEnum: Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 23066) @@ -125,7 +125,7 @@ /*Add input*/ element->AddInput(DamageFEnum,f,element->GetElementType()); - + /*Clean up and return*/ xDelete(f); }/*}}}*/ @@ -169,7 +169,7 @@ Gauss* gauss=element->NewGauss(); for (int i=0;iGaussNode(element->GetElementType(),i); - + eps_xx_input->GetInputValue(&eps_xx,gauss); eps_xy_input->GetInputValue(&eps_xy,gauss); eps_yy_input->GetInputValue(&eps_yy,gauss); @@ -176,7 +176,7 @@ B_input->GetInputValue(&B,gauss); 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)); eps2=(eps_xx+eps_yy)/2.-sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2)); @@ -194,7 +194,7 @@ /*Add input*/ element->AddInput(DamageFEnum,f,element->GetElementType()); - + /*Clean up and return*/ xDelete(f); delete gauss; @@ -261,7 +261,7 @@ Gauss* gauss=element->NewGauss(); for (int i=0;iGaussNode(element->GetElementType(),i); - + damage_input->GetInputValue(&damage,gauss); tau_xx_input->GetInputValue(&tau_xx,gauss); tau_xy_input->GetInputValue(&tau_xy,gauss); @@ -313,7 +313,7 @@ } /*Add input*/ element->AddInput(DamageFEnum,f,element->GetElementType()); - + /*Clean up and return*/ xDelete(f); delete gauss; @@ -374,7 +374,7 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); - + vx_input->GetInputValue(&vx,gauss); vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); vy_input->GetInputValue(&vy,gauss); @@ -537,7 +537,6 @@ damaged_input = element->GetInput(DamageDEnum); _assert_(damaged_input); } - /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ Index: ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 23066) @@ -7,7 +7,6 @@ /*Model processing*/ void Balancethickness2Analysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ - int finiteelement = P1Enum; IoModelToConstraintsx(constraints,iomodel,"md.balancethickness.spcthickness",Balancethickness2AnalysisEnum,finiteelement); Index: ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 23066) @@ -102,7 +102,7 @@ bool dakota_analysis,ismovingfront,isenthalpy; int frictionlaw,basalforcing_model,materialstype; int FrictionCoupling; - + /*Now, is the model 3d? otherwise, do nothing: */ if(iomodel->domaintype==Domain2DhorizontalEnum)return; @@ -188,7 +188,7 @@ default: _error_("not supported"); } - + /*Friction law variables*/ switch(frictionlaw){ case 1: @@ -804,7 +804,7 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); - + /*viscous dissipation*/ element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input); @@ -822,7 +822,7 @@ enthalpypicard_input->GetInputDerivativeValue(&d1enthalpypicard[0],xyz_list,gauss); 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]; d1P2=0.; @@ -1055,7 +1055,7 @@ IssmDouble dt; int* basalnodeindices=NULL; Element* element= NULL; - + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); /*depth-integrate the drained water fraction */ @@ -1387,7 +1387,7 @@ for(int i=0;iGaussNode(element->GetElementType(),indices[i]); gaussup->GaussNode(element->GetElementType(),indicesup[i]); - + enthalpy_input->GetInputValue(&enthalpy,gauss); enthalpy_input->GetInputValue(&enthalpyup,gaussup); pressure_input->GetInputValue(&pressure,gauss); Index: ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 23066) @@ -144,7 +144,7 @@ workelement->JacobianDeterminant(&Jdet,xyz_list,gauss); GetB(B,workelement,xyz_list,gauss,dim); GetBprime(Bprime,workelement,xyz_list,gauss,dim); - + D_scalar=gauss->weight*Jdet; if(extrapolatebydiffusion){ @@ -173,7 +173,7 @@ 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: ../trunk-jpl/src/c/main/issm_slr.cpp =================================================================== --- ../trunk-jpl/src/c/main/issm_slr.cpp (revision 23065) +++ ../trunk-jpl/src/c/main/issm_slr.cpp (revision 23066) @@ -27,7 +27,7 @@ /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/ worldcomm=EnvironmentInit(argc,argv); - + /*What is my rank?:*/ ISSM_MPI_Comm_rank(worldcomm,&my_rank); @@ -39,11 +39,11 @@ rankzeros=xNew(nummodels); 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); modelnames[i]=string; @@ -93,7 +93,7 @@ /*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)); femmodel->parameters->AddObject(new IntParam(NumModelsEnum,nummodels)); Index: ../trunk-jpl/src/c/main/issm_dakota.cpp =================================================================== --- ../trunk-jpl/src/c/main/issm_dakota.cpp (revision 23065) +++ ../trunk-jpl/src/c/main/issm_dakota.cpp (revision 23066) @@ -15,7 +15,6 @@ int main(int argc,char **argv){ - #if defined(_HAVE_DAKOTA_) && _DAKOTA_MAJOR_ >= 6 bool parallel=true; @@ -28,11 +27,11 @@ /*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"); @@ -56,7 +55,7 @@ << std::endl; Dakota::abort_handler(-1); } - + Dakota::ProblemDescDB& problem_db = env.problem_description_db(); Dakota::ModelLIter ml_iter; size_t model_index = problem_db.get_db_model_node(); // for restoration Index: ../trunk-jpl/src/c/main/issm_ocean.cpp =================================================================== --- ../trunk-jpl/src/c/main/issm_ocean.cpp (revision 23065) +++ ../trunk-jpl/src/c/main/issm_ocean.cpp (revision 23066) @@ -21,7 +21,7 @@ /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/ worldcomm=EnvironmentInit(argc,argv); - + /*What is my rank?:*/ ISSM_MPI_Comm_rank(worldcomm,&my_rank); ISSM_MPI_Comm_size(worldcomm,&my_size); @@ -38,7 +38,7 @@ ISSM_MPI_Intercomm_create( modelcomm, 0, worldcomm, my_local_size, 0, &tomitgcmcomm); 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)); femmodel->parameters->AddObject(new GenericParam(tomitgcmcomm,ToMITgcmCommEnum)); Index: ../trunk-jpl/src/c/main/esmfbinders.cpp =================================================================== --- ../trunk-jpl/src/c/main/esmfbinders.cpp (revision 23065) +++ ../trunk-jpl/src/c/main/esmfbinders.cpp (revision 23066) @@ -2,7 +2,6 @@ * \brief: ESMF binders for ISSM. Binders developed initially for the GEOS-5 framework. */ - #include "./issm.h" /*GCM specific declarations:*/ @@ -35,7 +34,7 @@ /*Some specific code here for the binding: */ femmodel->parameters->SetParam(SMBgcmEnum,SmbEnum); //bypass SMB model, will be provided by GCM! - + /*Restart file: */ femmodel->Restart(); @@ -112,11 +111,11 @@ { 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; } @@ -146,7 +145,7 @@ /*Output results: */ OutputResultsx(femmodel); - + /*Check point: */ femmodel->CheckPoint(); Index: ../trunk-jpl/src/c/bamg/BamgQuadtree.cpp =================================================================== --- ../trunk-jpl/src/c/bamg/BamgQuadtree.cpp (revision 23065) +++ ../trunk-jpl/src/c/bamg/BamgQuadtree.cpp (revision 23066) @@ -284,7 +284,6 @@ int xiv = b->v[k]->i.x; int yiv = b->v[k]->i.y; - int h0 = Norm(xi2,xiv,yi2,yiv); /*is it smaller than previous value*/ Index: ../trunk-jpl/src/c/bamg/Mesh.cpp =================================================================== --- ../trunk-jpl/src/c/bamg/Mesh.cpp (revision 23065) +++ ../trunk-jpl/src/c/bamg/Mesh.cpp (revision 23066) @@ -500,7 +500,7 @@ int* connectivitysize_1=NULL; int connectivitymax_1=0; int verbose=0; - + /*Check needed pointer*/ _assert_(bamgmesh); @@ -1033,7 +1033,7 @@ // ---------------- ! // s0 tt2 s1 //------------------------------- - + /*Intermediaries*/ Triangle* tt[3]; //the three triangles long long det3local[3]; //three determinants (integer) @@ -1157,10 +1157,10 @@ 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"); @@ -1853,7 +1853,7 @@ /*Check pointer*/ _assert_(bamgopts); - + /*Recover options*/ verbose=bamgopts->verbose; NbJacobi=bamgopts->nbjacobi; @@ -2727,7 +2727,7 @@ long k0=this->RandomNumber(nbv); if (verbose>4) _printf_(" Prime Number = "<4) _printf_(" k0 = "< (BAMG v1.01, Metric.cpp/MaxSubDivision)*/ - + /*Intermediaries*/ int verbose=0; @@ -3739,10 +3739,10 @@ /*}}}*/ 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; @@ -3784,7 +3784,7 @@ /*}}}*/ void Mesh::SmoothMetric(BamgOpts* bamgopts,double raisonmax) { /*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Metric.cpp/SmoothMetric)*/ - + /*Intermediaries*/ int verbose=0; @@ -4086,7 +4086,7 @@ int verbose=0; int i; Metric M1(1); - + /*Get options*/ if(bamgopts) verbose=bamgopts->verbose; Index: ../trunk-jpl/src/c/shared/Elements/Paterson.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Elements/Paterson.cpp (revision 23065) +++ ../trunk-jpl/src/c/shared/Elements/Paterson.cpp (revision 23066) @@ -15,7 +15,6 @@ /*Switch to celsius from Kelvin: */ T=temperature-273.15; - if(T<=-45.0){ B=1.e+8*(-0.000292866376675*pow(T+50.,3)+ 0.011672640664130*pow(T+50.,2) -0.325004442485481*(T+50.)+ 6.524779401948101); } Index: ../trunk-jpl/src/c/shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp (revision 23065) +++ ../trunk-jpl/src/c/shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp (revision 23066) @@ -15,7 +15,6 @@ IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12]; IssmDouble tdiffh; - for (int imonth = 0; imonth<12; imonth++){ tdiffh = TdiffTime*( TemperaturesLgm[imonth] - TemperaturesPresentday[imonth] ); monthlytemperaturestmp[imonth] = tdiffh + TemperaturesPresentday[imonth] ; @@ -28,4 +27,3 @@ } // printf(" tempera %f\n",monthlytemperaturestmp[1]); } - Index: ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 23065) +++ ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 23066) @@ -65,7 +65,7 @@ /*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: ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp (revision 23065) +++ ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp (revision 23066) @@ -10,18 +10,18 @@ bool isPrecipScaled, IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday, IssmDouble* PrecipitationReconstructed,IssmDouble* TemperatureReconstructed, IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout){ - + IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12]; IssmDouble deltaTemp; - + /* Constants */ // dpermil = 2.4;/*degrees C per mil*/ - + /*Create Delta Temp to be applied to monthly temps and used in precip scaling*/ deltaTemp = dpermil * (d018+34.83); - + for (int imonth = 0; imonth<12; imonth++){ - + if (isTemperatureScaled)monthlytemperaturestmp[imonth] = TemperaturePresentday[imonth] + deltaTemp; else{ monthlytemperaturestmp[imonth] = TemperatureReconstructed[imonth]; @@ -30,7 +30,7 @@ if (isPrecipScaled)monthlyprectmp[imonth] = PrecipitationPresentday[imonth]*exp((f/dpermil)*deltaTemp); else monthlyprectmp[imonth] = PrecipitationReconstructed[imonth]; - + /*Assign output pointer*/ *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth]; *(monthlyprecout+imonth) = monthlyprectmp[imonth]; Index: ../trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp (revision 23065) +++ ../trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp (revision 23066) @@ -18,7 +18,7 @@ /*get drainage function value*/ if((w0==w1)||(w1==w2)||(w0==w2)) _error_("Error: equal ordinates in DrainageFunctionWaterfraction -> division by zero. Abort"); - + if(waterfraction<=w0) Dret=D0; else if((waterfraction>w0) && (waterfraction<=w1)) Index: ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 23065) +++ ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 23066) @@ -82,11 +82,11 @@ // seasonal loop 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); tstar = monthlytemperatures[imonth]; @@ -98,12 +98,12 @@ if (tstar >= -siglimc){ pd = pds[reCast(tstar/DT + siglim0c)];}} 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; qmp= qmp + qmpt; @@ -113,7 +113,7 @@ // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of // 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){ pddsig=pdds[reCast(tstar/DT + siglim0)]; @@ -120,7 +120,7 @@ pdd = pdd + pddsig*deltm; frzndd = frzndd - (tstar-pddsig)*deltm;} else{frzndd = frzndd - tstar*deltm; } - + /*Assign output pointer*/ *(monthlytemperatures+imonth) = monthlytemperatures[imonth]; *(monthlyprec+imonth) = monthlyprec[imonth]; @@ -149,7 +149,7 @@ icefac=pddicefac; } } - + /***** determine PDD factors *****/ if(Tsum< -1.) { snwmf=(2.65+snowfac-pddsnowfac0)*0.001; // ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd @@ -165,12 +165,12 @@ } 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 // l 2.2= capillary factor Index: ../trunk-jpl/src/c/shared/Numerics/legendre.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Numerics/legendre.cpp (revision 23065) +++ ../trunk-jpl/src/c/shared/Numerics/legendre.cpp (revision 23066) @@ -120,11 +120,11 @@ v(1:m,2) = x; for i = 2 : n - + v(:,i+1) = ( ( 2 * i - 1 ) * x .* v(:,i) ... - ( i - 1 ) * v(:,i-1) ) ... / ( i ); - + end }}} */ @@ -239,7 +239,7 @@ Output, double P_POLYNOMIAL_VALUE[M*(N+1)], the values of the Legendre polynomials of order 0 through N. }}}*/ - + int i,j; if(n<0) return NULL; @@ -253,7 +253,7 @@ 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)] ) / ( IssmDouble ) ( j ); Index: ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp (revision 23065) +++ ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp (revision 23066) @@ -67,7 +67,7 @@ /*Get current Gradient at xmin=0*/ _printf0_(" x = "<(fxmin)) _error_("Function evaluation returned NaN"); - + /*Get f(xmax)*/ _printf0_(" x = "<(G); J[n]=fxbest; } - + /*return*/ xDelete(X); *pJ=J; Index: ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 23065) +++ ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 23066) @@ -62,7 +62,7 @@ else{ 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 (a*b)*c requires ac(b+d) multiplies and ac IssmDoubles, and a*(b*c) @@ -334,7 +334,7 @@ /*Compute determinant*/ 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; @@ -426,7 +426,6 @@ *pvx = vx; *pvy = vy; - }/*}}}*/ void Matrix3x3Determinant(IssmDouble* Adet,IssmDouble* A){/*{{{*/ @@ -576,16 +575,16 @@ }/*}}}*/ 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){ dummy[0]=newvalue; @@ -595,10 +594,10 @@ dummy[m]=newvalue; for(int i=0;i(cell); cell=dummy; - + /*assign output pointer:*/ *pcell=cell; } /*}}}*/ @@ -614,7 +613,7 @@ /*input: */ IssmDouble* cell=*pcell; - + /*output: */ IssmDouble* newcell=xNew(nind); @@ -621,7 +620,7 @@ for(int i=0;i(cell); @@ -632,7 +631,7 @@ /*input: */ IssmDouble* cell=*pcell; - + /*output: */ IssmDouble* newcell=xNew(m+1); @@ -640,7 +639,7 @@ newcell[i]=scale*cell[i]; newcell[i+1]=scale* cell[i]; for(int j=i+2;j(cell); @@ -661,7 +660,7 @@ celllist[x]= va_arg ( arguments, IssmDouble*); } va_end ( arguments ); - + _printf_("Echo of cell: \n"); for(int i=0;ilithosphere_density=this->lithosphere_density; matpar->mantle_shear_modulus=this->mantle_shear_modulus; matpar->mantle_density=this->mantle_density; - + matpar->earth_density=this->earth_density; return matpar; @@ -438,7 +438,7 @@ MARSHALLING(lithosphere_density); MARSHALLING(mantle_shear_modulus); MARSHALLING(mantle_density); - + //slr: MARSHALLING(earth_density); } Index: ../trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23066) @@ -195,7 +195,7 @@ /*Check to prevent dividing by zero if vmag==0*/ 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,7 +270,6 @@ r=drag_q/drag_p; s=1./drag_p; - /*Get effective pressure*/ IssmDouble Neff = EffectivePressure(gauss); Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23065) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23066) @@ -242,7 +242,6 @@ char *lockfilename = NULL; bool waitonlock = false; - /*Write lock file if requested: */ this->parameters->FindParam(&waitonlock,SettingsWaitonlockEnum); this->parameters->FindParam(&lockfilename,LockFileNameEnum); @@ -4392,7 +4391,6 @@ /*Free ressources:*/ xDelete(RSLg_old); - } /*}}}*/ void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ @@ -4806,7 +4804,6 @@ char** poutputbuffer; size_t* poutputbuffersize; - /*Before we delete the profiler, report statistics for this run: */ profiler->Stop(TOTAL); //final tagging _printf0_("\n"); @@ -4868,7 +4865,6 @@ }/*}}}*/ #endif - #if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_) void FemModel::ReMeshBamg(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist){/*{{{*/ @@ -5245,7 +5241,7 @@ newx[i] = newxylist[2*i]; newy[i] = newxylist[2*i+1]; } - + /*Assign the pointers*/ (*pnewelementslist) = newelementslist; //Matlab indexing (*pnewx) = newx; Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23065) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23066) @@ -2735,7 +2735,6 @@ } } - /*Clean-up*/ xDelete(dbasis); }/*}}}*/