Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26051) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26052) @@ -366,7 +366,4 @@ virtual void WriteFieldIsovalueSegment(DataSet* segments,int fieldenum,IssmDouble fieldvalue){_error_("not implemented yet");}; - #ifdef _HAVE_GIA_ - virtual void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0; - #endif #ifdef _HAVE_ESA_ virtual void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast, Vector* pX, Vector* pY,IssmDouble* xx,IssmDouble* yy)=0; @@ -385,4 +382,5 @@ virtual void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* mask)=0; virtual void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks)=0; + virtual void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26051) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26052) @@ -207,7 +207,4 @@ #endif - #ifdef _HAVE_GIA_ - void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; - #endif #ifdef _HAVE_ESA_ void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; @@ -224,4 +221,5 @@ void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; + void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26051) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26052) @@ -163,7 +163,4 @@ IssmDouble GetAreaSpherical(void){_error_("not implemented yet!");}; -#ifdef _HAVE_GIA_ - void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; -#endif #ifdef _HAVE_ESA_ void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; @@ -180,4 +177,5 @@ void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; + void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26051) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26052) @@ -169,7 +169,4 @@ void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); -#ifdef _HAVE_GIA_ - void GiaDeflection(Vector* wg,Vector* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; -#endif #ifdef _HAVE_ESA_ void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; @@ -186,4 +183,5 @@ void DeformationFromSurfaceLoads(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; + void GiaDeflection(Vector* wg,Vector* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26051) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26052) @@ -5088,115 +5088,4 @@ /*}}}*/ -#ifdef _HAVE_GIA_ -void Tria::GiaDeflection(Vector* wg,Vector* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/ - - IssmDouble xyz_list[NUMVERTICES][3]; - - /*gia solution parameters:*/ - IssmDouble ice_mask; - - /*output: */ - IssmDouble wi; - IssmDouble dwidt; - - /*arguments to GiaDeflectionCorex: */ - GiaDeflectionCoreArgs arguments; - - /*how many dofs are we working with here? */ - int gsize; - IssmDouble yts; - this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); - this->parameters->FindParam(&yts,ConstantsYtsEnum); - - /*recover gia solution parameters: */ - int cross_section_shape; - this->parameters->FindParam(&cross_section_shape,SolidearthSettingsCrossSectionShapeEnum); - - /*what time is it? :*/ - IssmDouble currenttime; - this->parameters->FindParam(¤ttime,TimeEnum); - - /*recover material parameters: */ - IssmDouble rho_ice = FindParam(MaterialsRhoIceEnum); - - /*recover mantle and lithosphere material properties:*/ - int numlayers=litho->numlayers; - - /*lithosphere is the last layer, mantle is the penultimate layer. Watch out, radius represents the layers - *from center to surface of the Earth:*/ - IssmDouble lithosphere_thickness = litho->radius[numlayers] - litho->radius[numlayers-1]; - IssmDouble lithosphere_shear_modulus = litho->lame_mu[numlayers-1]; - IssmDouble lithosphere_density = litho->density[numlayers-1]; - IssmDouble mantle_shear_modulus = litho->lame_mu[numlayers-2]; - IssmDouble mantle_density = litho->density[numlayers-2]; - IssmDouble mantle_viscosity = litho->viscosity[numlayers-2]; - - /*early return if we are NOT on an icy element:*/ - if(!IsIceInElement()) return; - - /*pull thickness averages! */ - IssmDouble *hes = NULL; - IssmDouble *times = NULL; - int numtimes; - this->GetInputAveragesUpToCurrentTime(TransientAccumulatedDeltaIceThicknessEnum,&hes,×,&numtimes,currenttime); - - if(this->Id()==1){ - _printf_("numtimes: " << numtimes << "\n"); - for (int i=0;iGetArea(); - - /*element radius: */ - IssmDouble re=sqrt(area/PI); - - /*figure out gravity center of our element: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - IssmDouble x0=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; - IssmDouble y0=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; - - /*start loading GiaDeflectionCore arguments: */ - arguments.re=re; - arguments.hes=hes; - arguments.times=times; - arguments.numtimes=numtimes; - arguments.currenttime=currenttime; - arguments.lithosphere_shear_modulus=lithosphere_shear_modulus; - arguments.lithosphere_density=lithosphere_density; - arguments.mantle_shear_modulus=mantle_shear_modulus; - arguments.mantle_viscosity=mantle_viscosity; - arguments.mantle_density=mantle_density; - arguments.lithosphere_thickness=lithosphere_thickness; - arguments.rho_ice=rho_ice; - arguments.idisk=this->id; - arguments.iedge=cross_section_shape; - arguments.yts=yts; - - for(int i=0;iSetValue(i,wi,ADD_VAL); - dwgdt->SetValue(i,dwidt,ADD_VAL); - } - - /*Free ressources: */ - xDelete(hes); - xDelete(times); - - return; -} -/*}}}*/ -#endif #ifdef _HAVE_ESA_ void Tria::EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){ /*{{{*/ @@ -6195,4 +6084,113 @@ } /*}}}*/ +void Tria::GiaDeflection(Vector* wg,Vector* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + + /*gia solution parameters:*/ + IssmDouble ice_mask; + + /*output: */ + IssmDouble wi; + IssmDouble dwidt; + + /*arguments to GiaDeflectionCorex: */ + GiaDeflectionCoreArgs arguments; + + /*how many dofs are we working with here? */ + int gsize; + IssmDouble yts; + this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + + /*recover gia solution parameters: */ + int cross_section_shape; + this->parameters->FindParam(&cross_section_shape,SolidearthSettingsCrossSectionShapeEnum); + + /*what time is it? :*/ + IssmDouble currenttime; + this->parameters->FindParam(¤ttime,TimeEnum); + + /*recover material parameters: */ + IssmDouble rho_ice = FindParam(MaterialsRhoIceEnum); + + /*recover mantle and lithosphere material properties:*/ + int numlayers=litho->numlayers; + + /*lithosphere is the last layer, mantle is the penultimate layer. Watch out, radius represents the layers + *from center to surface of the Earth:*/ + IssmDouble lithosphere_thickness = litho->radius[numlayers] - litho->radius[numlayers-1]; + IssmDouble lithosphere_shear_modulus = litho->lame_mu[numlayers-1]; + IssmDouble lithosphere_density = litho->density[numlayers-1]; + IssmDouble mantle_shear_modulus = litho->lame_mu[numlayers-2]; + IssmDouble mantle_density = litho->density[numlayers-2]; + IssmDouble mantle_viscosity = litho->viscosity[numlayers-2]; + + /*early return if we are NOT on an icy element:*/ + if(!IsIceInElement()) return; + + /*pull thickness averages! */ + IssmDouble *hes = NULL; + IssmDouble *times = NULL; + int numtimes; + this->GetInputAveragesUpToCurrentTime(TransientAccumulatedDeltaIceThicknessEnum,&hes,×,&numtimes,currenttime); + + if(this->Id()==1){ + _printf_("numtimes: " << numtimes << "\n"); + for (int i=0;iGetArea(); + + /*element radius: */ + IssmDouble re=sqrt(area/PI); + + /*figure out gravity center of our element: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + IssmDouble x0=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; + IssmDouble y0=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; + + /*start loading GiaDeflectionCore arguments: */ + arguments.re=re; + arguments.hes=hes; + arguments.times=times; + arguments.numtimes=numtimes; + arguments.currenttime=currenttime; + arguments.lithosphere_shear_modulus=lithosphere_shear_modulus; + arguments.lithosphere_density=lithosphere_density; + arguments.mantle_shear_modulus=mantle_shear_modulus; + arguments.mantle_viscosity=mantle_viscosity; + arguments.mantle_density=mantle_density; + arguments.lithosphere_thickness=lithosphere_thickness; + arguments.rho_ice=rho_ice; + arguments.idisk=this->id; + arguments.iedge=cross_section_shape; + arguments.yts=yts; + + for(int i=0;iSetValue(i,wi,ADD_VAL); + dwgdt->SetValue(i,dwidt,ADD_VAL); + } + + /*Free ressources: */ + xDelete(hes); + xDelete(times); + + return; +} +/*}}}*/ #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26051) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26052) @@ -154,7 +154,4 @@ void WriteFieldIsovalueSegment(DataSet* segments,int fieldenum,IssmDouble fieldvalue); - #ifdef _HAVE_GIA_ - void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y); - #endif #ifdef _HAVE_ESA_ void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast, Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy); @@ -171,4 +168,5 @@ void SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks); void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks); + void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y); #endif /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26051) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26052) @@ -981,5 +981,4 @@ if(profiler->Used(SMBCORE)) _printf0_(" "<TotalTime(SMBCORE) << " sec\n"); if(profiler->Used(GROUNDINGLINECORE)) _printf0_(" "<TotalTime(GROUNDINGLINECORE) << " sec\n"); - if(profiler->Used(GIACORE)) _printf0_(" "<TotalTime(GIACORE) << " sec\n"); if(profiler->Used(ESACORE)) _printf0_(" "<TotalTime(ESACORE) << " sec\n"); if(profiler->Used(SLRCORE)) _printf0_(" "<TotalTime(SLRCORE) << " sec\n"); Index: /issm/trunk-jpl/src/c/classes/Profiler.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Profiler.h (revision 26051) +++ /issm/trunk-jpl/src/c/classes/Profiler.h (revision 26052) @@ -25,17 +25,16 @@ #define SMBCORE 12 /*Profiling SMB */ #define GROUNDINGLINECORE 13 /*Profiling GROUDINGLINE MIGRATION */ -#define GIACORE 14 /*Profiling GIA */ -#define ESACORE 15 /*Profiling ESA */ -#define SLRCORE 16 /*Profiling SLR */ -#define SAMPLINGCORE 17 /*Profiling SAMPLING */ -#define MPISERIAL 18 /*Profiling MPISerial */ -#define SEDLOOP 19 /*Profiling MPISerial */ -#define SEDMatrix 20 /*Profiling MPISerial */ -#define SEDUpdate 21 /*Profiling MPISerial */ -#define EPLLOOP 22 /*Profiling MPISerial */ -#define EPLMasking 23 /*Profiling MPISerial */ -#define EPLMatrices 24 /*Profiling MPISerial */ -#define EPLUpdate 25 /*Profiling MPISerial */ -#define MAXPROFSIZE 26 /*Used to initialize static arrays*/ +#define ESACORE 14 /*Profiling ESA */ +#define SLRCORE 15 /*Profiling SLR */ +#define SAMPLINGCORE 16 /*Profiling SAMPLING */ +#define MPISERIAL 17 /*Profiling MPISerial */ +#define SEDLOOP 18 /*Profiling MPISerial */ +#define SEDMatrix 19 /*Profiling MPISerial */ +#define SEDUpdate 20 /*Profiling MPISerial */ +#define EPLLOOP 21 /*Profiling MPISerial */ +#define EPLMasking 22 /*Profiling MPISerial */ +#define EPLMatrices 23 /*Profiling MPISerial */ +#define EPLUpdate 24 /*Profiling MPISerial */ +#define MAXPROFSIZE 25 /*Used to initialize static arrays*/