Index: /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 22348) +++ /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 22349) @@ -56,4 +56,7 @@ int *transitions_N = NULL; int ntransitions; + + /*some constant parameters: */ + parameters->AddObject(iomodel->CopyConstantObject("md.esa.hemisphere",EsaHemisphereEnum)); /*love numbers: */ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22348) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22349) @@ -1750,4 +1750,6 @@ name==EsaNmotionEnum || name==EsaEmotionEnum || + name==EsaXmotionEnum || + name==EsaYmotionEnum || name==EsaStrainratexxEnum || name==EsaStrainratexyEnum || Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22348) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22349) @@ -310,5 +310,5 @@ #endif #ifdef _HAVE_ESA_ - virtual void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* xx,IssmDouble* yy)=0; + virtual void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast, Vector* pX, Vector* pY,IssmDouble* xx,IssmDouble* yy)=0; virtual void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22348) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22349) @@ -184,5 +184,5 @@ #endif #ifdef _HAVE_ESA_ - void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; + void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_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 22348) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 22349) @@ -167,5 +167,5 @@ #endif #ifdef _HAVE_ESA_ - void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; + void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_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 22348) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 22349) @@ -174,5 +174,5 @@ #endif #ifdef _HAVE_ESA_ - void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; + void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_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 22348) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22349) @@ -539,6 +539,6 @@ /*Retrieve all inputs we will be needing: */ - Input* vx_input=this->GetInput(EsaEmotionEnum); _assert_(vx_input); - Input* vy_input=this->GetInput(EsaNmotionEnum); _assert_(vy_input); + 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: */ @@ -3750,5 +3750,5 @@ #endif #ifdef _HAVE_ESA_ -void Tria::EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* xx,IssmDouble* yy){ /*{{{*/ +void Tria::EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){ /*{{{*/ /*diverse:*/ @@ -3763,5 +3763,5 @@ IssmDouble* U_elastic_precomputed = NULL; IssmDouble* H_elastic_precomputed = NULL; - int M; + int M, hemi; /*computation of Green functions:*/ @@ -3769,4 +3769,6 @@ IssmDouble* N_elastic= NULL; IssmDouble* E_elastic= NULL; + IssmDouble* X_elastic= NULL; + IssmDouble* Y_elastic= NULL; /*optimization:*/ @@ -3788,4 +3790,7 @@ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + /*which hemisphere? for north-south, east-west components*/ + this->parameters->FindParam(&hemi,EsaHemisphereEnum); + /*compute area of element:*/ area=GetArea(); @@ -3807,4 +3812,6 @@ N_elastic=xNewZeroInit(gsize); E_elastic=xNewZeroInit(gsize); + X_elastic=xNewZeroInit(gsize); + Y_elastic=xNewZeroInit(gsize); int* indices=xNew(gsize); @@ -3812,11 +3819,13 @@ IssmDouble* N_values=xNewZeroInit(gsize); IssmDouble* E_values=xNewZeroInit(gsize); - IssmDouble dx, dy; - IssmDouble dist, alpha, ang; - IssmDouble N_azim, E_azim; + IssmDouble* X_values=xNewZeroInit(gsize); + IssmDouble* Y_values=xNewZeroInit(gsize); + IssmDouble dx, dy, dist, alpha, ang, ang2; + IssmDouble N_azim, E_azim, X_azim, Y_azim; for(int i=0;i(alpha/PI*(M-1)); U_elastic[i] += U_elastic_precomputed[index]; - N_elastic[i] += H_elastic_precomputed[index]*N_azim; - E_elastic[i] += H_elastic_precomputed[index]*E_azim; + Y_elastic[i] += H_elastic_precomputed[index]*Y_azim; + X_elastic[i] += H_elastic_precomputed[index]*X_azim; /*Add all components to the pUp solution vectors:*/ U_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*U_elastic[i]; - N_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*N_elastic[i]; - E_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*E_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]); + } + else if (hemi == 1) { + ang2 = PI/2 - atan2(-yy[i],-xx[i]); + } + if (hemi != 0){ + N_azim = Y_azim*cos(ang2) + X_azim*sin(ang2); + E_azim = X_azim*cos(ang2) - Y_azim*sin(ang2); + N_elastic[i] += H_elastic_precomputed[index]*N_azim; + E_elastic[i] += H_elastic_precomputed[index]*E_azim; + N_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*N_elastic[i]; + E_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*E_elastic[i]; + } } @@ -3847,4 +3872,6 @@ pNorth->SetValues(gsize,indices,N_values,ADD_VAL); 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:*/ @@ -3852,4 +3879,6 @@ xDelete(U_values); xDelete(N_values); xDelete(E_values); xDelete(U_elastic); xDelete(N_elastic); xDelete(E_elastic); + xDelete(X_values); xDelete(Y_values); + xDelete(X_elastic); xDelete(Y_elastic); return; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22348) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22349) @@ -143,5 +143,5 @@ #endif #ifdef _HAVE_ESA_ - void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* xx,IssmDouble* yy); + void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast, Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy); void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea); #endif Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22348) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22349) @@ -3976,5 +3976,5 @@ #endif #ifdef _HAVE_ESA_ -void FemModel::EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* xx, IssmDouble* yy){/*{{{*/ +void FemModel::EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pX, Vector* pY, IssmDouble* xx, IssmDouble* yy){/*{{{*/ int ns,nsmax; @@ -3991,5 +3991,5 @@ if(i(elements->GetObjectByOffset(i)); - element->EsaGeodetic2D(pUp,pNorth,pEast,xx,yy); + element->EsaGeodetic2D(pUp,pNorth,pEast,pX,pY,xx,yy); } if(i%100==0){ @@ -3997,4 +3997,6 @@ pNorth->Assemble(); pEast->Assemble(); + pX->Assemble(); + pY->Assemble(); } } @@ -4004,4 +4006,6 @@ pNorth->Assemble(); pEast->Assemble(); + pX->Assemble(); + pY->Assemble(); /*Free ressources:*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22348) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22349) @@ -134,5 +134,5 @@ #endif #ifdef _HAVE_ESA_ - void EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* xx, IssmDouble* yy); + void EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pX, Vector* pY, IssmDouble* xx, IssmDouble* yy); void EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); #endif Index: /issm/trunk-jpl/src/c/cores/esa_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/esa_core.cpp (revision 22348) +++ /issm/trunk-jpl/src/c/cores/esa_core.cpp (revision 22349) @@ -15,4 +15,6 @@ Vector *U_north = NULL; Vector *U_east = NULL; + Vector *U_x = NULL; + Vector *U_y = NULL; bool save_results,isesa,iscoupler; int configuration_type; @@ -76,4 +78,6 @@ U_north = new Vector(gsize); U_east = new Vector(gsize); + U_x = new Vector(gsize); + U_y = new Vector(gsize); /*call the geodetic main modlule:*/ @@ -82,5 +86,7 @@ } if(domaintype==Domain2DhorizontalEnum){ - femmodel->EsaGeodetic2D(U_radial,U_north,U_east,xx,yy); + femmodel->EsaGeodetic2D(U_radial,U_north,U_east,U_x,U_y,xx,yy); + InputUpdateFromVectorx(femmodel,U_x,EsaXmotionEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,U_y,EsaYmotionEnum,VertexSIdEnum); } @@ -102,4 +108,6 @@ delete U_north; delete U_east; + delete U_x; + delete U_y; if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} } Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22348) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22349) @@ -872,4 +872,7 @@ EsaNmotionEnum, EsaEmotionEnum, + EsaXmotionEnum, + EsaYmotionEnum, + EsaHemisphereEnum, EsaStrainratexxEnum, EsaStrainratexyEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22348) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22349) @@ -847,4 +847,7 @@ case EsaNmotionEnum : return "EsaNmotion"; case EsaEmotionEnum : return "EsaEmotion"; + case EsaXmotionEnum : return "EsaXmotion"; + case EsaYmotionEnum : return "EsaYmotion"; + case EsaHemisphereEnum : return "EsaHemisphere"; case EsaStrainratexxEnum : return "EsaStrainratexx"; case EsaStrainratexyEnum : return "EsaStrainratexy"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22348) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22349) @@ -865,4 +865,7 @@ else if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum; else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum; + else if (strcmp(name,"EsaXmotion")==0) return EsaXmotionEnum; + else if (strcmp(name,"EsaYmotion")==0) return EsaYmotionEnum; + else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum; else if (strcmp(name,"EsaStrainratexx")==0) return EsaStrainratexxEnum; else if (strcmp(name,"EsaStrainratexy")==0) return EsaStrainratexyEnum; @@ -872,11 +875,11 @@ else if (strcmp(name,"EsaUElastic")==0) return EsaUElasticEnum; else if (strcmp(name,"EsaHElastic")==0) return EsaHElasticEnum; - else if (strcmp(name,"EsaTransitions")==0) return EsaTransitionsEnum; - else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum; - else if (strcmp(name,"EsaNumRequestedOutputs")==0) return EsaNumRequestedOutputsEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"TransientAmrFrequency")==0) return TransientAmrFrequencyEnum; + if (strcmp(name,"EsaTransitions")==0) return EsaTransitionsEnum; + else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum; + else if (strcmp(name,"EsaNumRequestedOutputs")==0) return EsaNumRequestedOutputsEnum; + else if (strcmp(name,"TransientAmrFrequency")==0) return TransientAmrFrequencyEnum; else if (strcmp(name,"AmrType")==0) return AmrTypeEnum; else if (strcmp(name,"AmrRestart")==0) return AmrRestartEnum; @@ -995,11 +998,11 @@ else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum; - else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum; - else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum; - else if (strcmp(name,"HydrologySommersAnalysis")==0) return HydrologySommersAnalysisEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum; + if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum; + else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum; + else if (strcmp(name,"HydrologySommersAnalysis")==0) return HydrologySommersAnalysisEnum; + else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum; else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum; else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; Index: /issm/trunk-jpl/src/m/classes/esa.m =================================================================== --- /issm/trunk-jpl/src/m/classes/esa.m (revision 22348) +++ /issm/trunk-jpl/src/m/classes/esa.m (revision 22349) @@ -10,4 +10,5 @@ love_l = 0; %ideam degacc = 0; + hemisphere = 0; requested_outputs = {}; transitions = {}; @@ -26,5 +27,8 @@ %numerical discretization accuracy self.degacc=.01; - + + %computational flags: + self.hemisphere=0; + %output default: self.requested_outputs={'default'}; @@ -66,4 +70,5 @@ fielddisplay(self,'love_h','load Love number for radial displacement'); fielddisplay(self,'love_l','load Love number for horizontal displacements'); + fielddisplay(self,'hemisphere','North-south, East-west components of 2-D horiz displacement vector: -1 south, 1 north'); fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'); fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); @@ -75,4 +80,5 @@ WriteData(fid,prefix,'object',self,'fieldname','love_h','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','love_l','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','hemisphere','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray'); @@ -93,4 +99,5 @@ writejs1Darray(fid,[modelname '.esa.love_h'],self.love_h); writejs1Darray(fid,[modelname '.esa.love_l'],self.love_l); + writejsdouble(fid,[modelname '.esa.hemisphere'],self.hemisphere); writejsdouble(fid,[modelname '.esa.degacc'],self.degacc); writejscellstring(fid,[modelname '.esa.requested_outputs'],self.requested_outputs);