Index: /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 14806) +++ /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 14807) @@ -477,4 +477,7 @@ StressTensoryzEnum, StressTensorzzEnum, + GiaOutputRatesEnum, + GiaCrossSectionShapeEnum, + GiadWdtEnum, GiaWEnum, /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 14806) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 14807) @@ -1544,5 +1544,5 @@ #endif #ifdef _HAVE_GIA_ -void FemModel::Deflection(Vector* wg,IssmDouble* x, IssmDouble* y){ /*{{{*/ +void FemModel::Deflection(Vector* wg,Vector* dwgdt, IssmDouble* x, IssmDouble* y){ /*{{{*/ int i; @@ -1554,7 +1554,7 @@ for (i=0;iSize();i++){ element=dynamic_cast(elements->GetObjectByOffset(i)); - element->GiaDeflection(wg,x,y); - } -} -/*}}}*/ -#endif + element->GiaDeflection(wg,dwgdt, x,y); + } +} +/*}}}*/ +#endif Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 14806) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 14807) @@ -94,5 +94,5 @@ #endif #ifdef _HAVE_GIA_ - void Deflection(Vector* wg,IssmDouble* x, IssmDouble* y); + void Deflection(Vector* wg,Vector* dwgdt, IssmDouble* x, IssmDouble* y); #endif void SystemMatricesx(Matrix** pKff, Matrix** pKfs, Vector** ppf, Vector** pdf, IssmDouble* pkmax); Index: /issm/trunk-jpl/src/c/classes/GiaDeflectionCoreArgs.h =================================================================== --- /issm/trunk-jpl/src/c/classes/GiaDeflectionCoreArgs.h (revision 14806) +++ /issm/trunk-jpl/src/c/classes/GiaDeflectionCoreArgs.h (revision 14807) @@ -25,4 +25,8 @@ IssmDouble lithosphere_thickness; + /*gia solution parameters: */ + int iedge; + int irate; + /*ice properties: */ IssmDouble rho_ice; Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Element.h (revision 14806) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Element.h (revision 14807) @@ -101,5 +101,5 @@ #ifdef _HAVE_GIA_ - virtual void GiaDeflection(Vector* wg,IssmDouble* x,IssmDouble* y)=0; + virtual void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y)=0; #endif Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 14806) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 14807) @@ -2171,5 +2171,8 @@ name==QmuSurfaceEnum || name==QmuTemperatureEnum || - name==QmuMeltingEnum + name==QmuMeltingEnum || + name==GiaWEnum || + name==GiadWdtEnum + ) { return true; @@ -3498,5 +3501,5 @@ #ifdef _HAVE_GIA_ /*FUNCTION Penta::GiaDeflection {{{*/ -void Penta::GiaDeflection(Vector* wg,IssmDouble* x, IssmDouble* y){ +void Penta::GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y){ _error_("GIA deflection not implemented yet!"); } Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h (revision 14806) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h (revision 14807) @@ -144,5 +144,5 @@ #ifdef _HAVE_GIA_ - void GiaDeflection(Vector* wg,IssmDouble* x,IssmDouble* y); + void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y); #endif Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 14806) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 14807) @@ -2060,5 +2060,7 @@ name==QmuSurfaceEnum || name==QmuTemperatureEnum || - name==QmuMeltingEnum + name==QmuMeltingEnum || + name==GiaWEnum || + name==GiadWdtEnum ){ return true; @@ -3072,5 +3074,5 @@ #ifdef _HAVE_GIA_ /*FUNCTION Tria::GiaDeflection {{{*/ -void Tria::GiaDeflection(Vector* wg,IssmDouble* x, IssmDouble* y){ +void Tria::GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x, IssmDouble* y){ int i; @@ -3087,4 +3089,8 @@ int numtimes; Input* thickness_input=NULL; + + /*gia solution parameters:*/ + int output_rates=0; + int cross_section_shape=0; /*gia material parameters: */ @@ -3102,4 +3108,5 @@ /*output: */ IssmDouble wi; + IssmDouble dwidt; /*arguments to GiaDeflectionCorex: */ @@ -3108,4 +3115,8 @@ /*how many dofs are we working with here? */ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + + /*recover gia solution parameters: */ + this->parameters->FindParam(&output_rates,GiaOutputRatesEnum); + this->parameters->FindParam(&cross_section_shape,GiaCrossSectionShapeEnum); /*what time is it? :*/ @@ -3157,4 +3168,6 @@ arguments.rho_ice=rho_ice; arguments.idisk=this->id; + arguments.irate=output_rates; + arguments.iedge=cross_section_shape; for(i=0;iSetValue(i,wi,ADD_VAL); + dwgdt->SetValue(i,dwidt,ADD_VAL); } Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 14806) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 14807) @@ -142,5 +142,5 @@ #ifdef _HAVE_GIA_ - void GiaDeflection(Vector* wg,IssmDouble* x,IssmDouble* y); + void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y); void InputUpdateFromSolutionGia(IssmDouble* solution); #endif Index: /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 14806) +++ /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 14807) @@ -468,4 +468,7 @@ case StressTensoryzEnum : return "StressTensoryz"; case StressTensorzzEnum : return "StressTensorzz"; + case GiaOutputRatesEnum : return "GiaOutputRates"; + case GiaCrossSectionShapeEnum : return "GiaCrossSectionShape"; + case GiadWdtEnum : return "GiadWdt"; case GiaWEnum : return "GiaW"; case P0Enum : return "P0"; Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 14806) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 14807) @@ -68,17 +68,20 @@ /*}}}*/ -void GiaDeflectionCorex( IssmDouble* pwi, GiaDeflectionCoreArgs* arguments){ +void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble* pdwidt, GiaDeflectionCoreArgs* arguments){ /*output: */ IssmDouble wi=0; + IssmDouble dwidt=0; - /*inputs: */ - int iedge=1; //c iedge ......... = 1 square-edged, = 2 elliptical disc x-section (see naruse.f) + /*inputs: {{{*/ + /*constant: */ int idisk=1; // disk # - int irate=0; // =0 fetch w solution (m) only; =1 dw/dt (mm/yr) only - /*intermediary: */ + /*coming from the model structure, runtime configurable:*/ + /*gia solution parameters: */ + int iedge=0; + int irate=0; - /*inputs: {{{*/ + /*gia inputs: */ IssmDouble ri; //radial distance from center of disk to vertex i IssmDouble re; //radius of disk @@ -121,6 +124,16 @@ rho_ice=arguments->rho_ice; disk_id=arguments->idisk; + irate=arguments->irate; + iedge=arguments->iedge; + + //printf("%g %g %g %i %g %g %g %g %g %g %g %g %i\n", ri,re,current_he,numtimes,currenttime,lithosphere_shear_modulus,lithosphere_density,mantle_shear_modulus,mantle_viscosity, mantle_density,lithosphere_thickness,rho_ice,disk_id); + /*}}}*/ + /*Modify inputs to match naruse code: */ + //from our model, irate comes in with values in [1,2], which maps into [0,1] in naruse: + irate=irate-1; + + /*Prepare block inputs for fortran distme and what0 routines of the naruse code: {{{*/ /*Now, let's set pset from the data that we got in input to GiaDeflectionCorex: */ blockp_.pset[0]=lithosphere_thickness; @@ -147,5 +160,7 @@ /*radial distance of i-th element: */ blockrad_.distrad=ri/1000.0; // in km + /*}}}*/ + /*Call distme driver: */ distme_(&idisk,&iedge); @@ -154,11 +169,17 @@ what0_(&idisk,&iedge); - /*this is the solution: */ - wi = blocks_.aswokm; + /*output solution: */ + if(irate==0){ + wi = blocks_.aswokm; + *pwi=wi; + *pdwidt=0; + } + else if (irate==1){ + dwidt = blocks_.aswokm; + *pdwidt=dwidt; + *pwi=0; + } - /*allocate output pointer: */ - *pwi=wi; - -// printf("wi: %g deflection: %g \n",wi,blocks_.aswokm); + //printf("deflection: %g deflection rate %g\n",wi, dwidt); } Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.h (revision 14806) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.h (revision 14807) @@ -10,5 +10,5 @@ /* local prototypes: */ -void GiaDeflectionCorex( IssmDouble* pwi, GiaDeflectionCoreArgs* arguments); +void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble* pdwidt, GiaDeflectionCoreArgs* arguments); #endif /* _GIADEFLECTIONCOREX_H */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 14806) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 14807) @@ -105,4 +105,6 @@ parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsdelta18oEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIssmbgradientsEnum)); + parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum)); + parameters->AddObject(iomodel->CopyConstantObject(GiaOutputRatesEnum)); iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum); Index: /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 14806) +++ /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 14807) @@ -478,4 +478,7 @@ else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum; + else if (strcmp(name,"GiaOutputRates")==0) return GiaOutputRatesEnum; + else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum; + else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum; else if (strcmp(name,"GiaW")==0) return GiaWEnum; else if (strcmp(name,"P0")==0) return P0Enum; @@ -504,11 +507,11 @@ else if (strcmp(name,"MaxVx")==0) return MaxVxEnum; else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum; - else if (strcmp(name,"MinVy")==0) return MinVyEnum; - else if (strcmp(name,"MaxVy")==0) return MaxVyEnum; - else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"MinVz")==0) return MinVzEnum; + if (strcmp(name,"MinVy")==0) return MinVyEnum; + else if (strcmp(name,"MaxVy")==0) return MaxVyEnum; + else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum; + else if (strcmp(name,"MinVz")==0) return MinVzEnum; else if (strcmp(name,"MaxVz")==0) return MaxVzEnum; else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum; Index: /issm/trunk-jpl/src/c/solutions/gia_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/gia_core.cpp (revision 14806) +++ /issm/trunk-jpl/src/c/solutions/gia_core.cpp (revision 14807) @@ -15,4 +15,5 @@ int i; Vector* wg = NULL; + Vector* dwgdt = NULL; IssmDouble* x = NULL; IssmDouble* y = NULL; @@ -35,4 +36,5 @@ gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); wg = new Vector(gsize); + dwgdt = new Vector(gsize); /*first, recover x and y vectors from vertices: */ @@ -40,14 +42,17 @@ /*call the main module: */ - femmodel->Deflection(wg,x,y); + femmodel->Deflection(wg,dwgdt,x,y); /*assemble vector: */ wg->Assemble(); + dwgdt->Assemble(); InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,wg); + InputUpdateFromVectorx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,dwgdt,GiadWdtEnum,VertexEnum); if(save_results){ if(VerboseSolution()) _pprintLine_(" saving results"); InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GiaWEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GiadWdtEnum); } Index: /issm/trunk-jpl/src/m/classes/gia.m =================================================================== --- /issm/trunk-jpl/src/m/classes/gia.m (revision 14806) +++ /issm/trunk-jpl/src/m/classes/gia.m (revision 14807) @@ -7,4 +7,6 @@ properties (SetAccess=public) lithosphere_thickness = NaN; + output_rates = 0; + cross_section_shape = 0; end methods @@ -18,5 +20,6 @@ end % }}} function obj = setdefaultparameters(obj) % {{{ - + obj.output_rates=1; %don't output rates by default (see irate in GiaDeflectionCorex) + obj.cross_section_shape=1; %square as default (see iedge in GiaDeflectionCorex) end % }}} function md = checkconsistency(obj,md,solution,analyses) % {{{ @@ -24,4 +27,6 @@ if ~ismember(GiaAnalysisEnum(),analyses), return; end md = checkfield(md,'gia.lithosphere_thickness','NaN',1,'size',[md.mesh.numberofvertices 1],'>',0); + md = checkfield(md,'gia.output_rates','numel',[1],'values',[1,2]); + md = checkfield(md,'gia.cross_section_shape','numel',[1],'values',[1,2]); end % }}} @@ -30,8 +35,12 @@ fielddisplay(obj,'lithosphere_thickness','lithosphere thickness[km]'); + fielddisplay(obj,'output_rates','1: fetch w solution (m) (default). 2: fetch dw/dt (mm/yr). See irate in GiaDeflectionCore'); + fielddisplay(obj,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore'); end % }}} function marshall(obj,fid) % {{{ WriteData(fid,'data',obj.lithosphere_thickness,'format','DoubleMat','mattype',1,'enum',GiaLithosphereThicknessEnum()); + WriteData(fid,'object',obj,'fieldname','output_rates','format','Integer'); + WriteData(fid,'object',obj,'fieldname','cross_section_shape','format','Integer'); end % }}} end Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 14806) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 14807) @@ -4519,4 +4519,34 @@ return StringToEnum('StressTensorzz')[0] +def GiaOutputRatesEnum(): + """ + GIAOUTPUTRATESENUM - Enum of GiaOutputRates + + Usage: + macro=GiaOutputRatesEnum() + """ + + return StringToEnum('GiaOutputRates')[0] + +def GiaCrossSectionShapeEnum(): + """ + GIACROSSSECTIONSHAPEENUM - Enum of GiaCrossSectionShape + + Usage: + macro=GiaCrossSectionShapeEnum() + """ + + return StringToEnum('GiaCrossSectionShape')[0] + +def GiadWdtEnum(): + """ + GIADWDTENUM - Enum of GiadWdt + + Usage: + macro=GiadWdtEnum() + """ + + return StringToEnum('GiadWdt')[0] + def GiaWEnum(): """ @@ -5427,4 +5457,4 @@ """ - return 541 - + return 544 + Index: /issm/trunk-jpl/src/m/enum/GiaCrossSectionShapeEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/GiaCrossSectionShapeEnum.m (revision 14807) +++ /issm/trunk-jpl/src/m/enum/GiaCrossSectionShapeEnum.m (revision 14807) @@ -0,0 +1,11 @@ +function macro=GiaCrossSectionShapeEnum() +%GIACROSSSECTIONSHAPEENUM - Enum of GiaCrossSectionShape +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=GiaCrossSectionShapeEnum() + +macro=StringToEnum('GiaCrossSectionShape'); Index: /issm/trunk-jpl/src/m/enum/GiaOutputRatesEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/GiaOutputRatesEnum.m (revision 14807) +++ /issm/trunk-jpl/src/m/enum/GiaOutputRatesEnum.m (revision 14807) @@ -0,0 +1,11 @@ +function macro=GiaOutputRatesEnum() +%GIAOUTPUTRATESENUM - Enum of GiaOutputRates +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=GiaOutputRatesEnum() + +macro=StringToEnum('GiaOutputRates'); Index: /issm/trunk-jpl/src/m/enum/GiadWdtEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/GiadWdtEnum.m (revision 14807) +++ /issm/trunk-jpl/src/m/enum/GiadWdtEnum.m (revision 14807) @@ -0,0 +1,11 @@ +function macro=GiadWdtEnum() +%GIADWDTENUM - Enum of GiadWdt +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=GiadWdtEnum() + +macro=StringToEnum('GiadWdt'); Index: /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 14806) +++ /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 14807) @@ -9,3 +9,3 @@ % macro=MaximumNumberOfEnums() -macro=541; +macro=544;