Index: /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 21271) +++ /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 21272) @@ -32,6 +32,4 @@ /*Create inputs: */ iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); - iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum); - iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum); iomodel->FetchDataToInput(elements,"md.esa.deltathickness",EsaDeltathicknessEnum); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 21271) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 21272) @@ -297,5 +297,6 @@ #endif #ifdef _HAVE_ESA_ - virtual void EsaGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea)=0; + virtual void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,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 #ifdef _HAVE_SEALEVELRISE_ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 21271) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 21272) @@ -184,5 +184,6 @@ #endif #ifdef _HAVE_ESA_ - void EsaGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,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 #ifdef _HAVE_SEALEVELRISE_ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 21271) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 21272) @@ -169,5 +169,6 @@ #endif #ifdef _HAVE_ESA_ - void EsaGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,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 #ifdef _HAVE_SEALEVELRISE_ Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 21271) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 21272) @@ -176,5 +176,6 @@ #endif #ifdef _HAVE_ESA_ - void EsaGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,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 #ifdef _HAVE_SEALEVELRISE_ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21271) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21272) @@ -3581,5 +3581,111 @@ #endif #ifdef _HAVE_ESA_ -void Tria::EsaGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/ +void Tria::EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* xx,IssmDouble* yy){ /*{{{*/ + + /*diverse:*/ + int gsize; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble area; + IssmDouble earth_radius = 6371012.0; // Earth's radius [m] + IssmDouble I; //ice/water loading + IssmDouble rho_ice, rho_earth; + + /*precomputed elastic green functions:*/ + 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; + + /*Compute ice thickness change: */ + 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; + + /*recover material parameters: */ + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + rho_earth=matpar->GetMaterialParameter(MaterialsEarthDensityEnum); + + /*how many dofs are we working with here? */ + this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + + /*compute area of element:*/ + area=GetArea(); + + /*figure out gravity center of our element (Cartesian): */ + IssmDouble x_element, y_element; + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; + y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; + + /*recover elastic Green's functions for displacement:*/ + DoubleVecParam* U_parameter = static_cast(this->parameters->FindParamObject(EsaUElasticEnum)); _assert_(U_parameter); + DoubleVecParam* H_parameter = static_cast(this->parameters->FindParamObject(EsaHElasticEnum)); _assert_(H_parameter); + U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M); + H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M); + + /*initialize: */ + U_elastic=xNewZeroInit(gsize); + N_elastic=xNewZeroInit(gsize); + E_elastic=xNewZeroInit(gsize); + + int* indices=xNew(gsize); + IssmDouble* U_values=xNewZeroInit(gsize); + IssmDouble* N_values=xNewZeroInit(gsize); + IssmDouble* E_values=xNewZeroInit(gsize); + IssmDouble dx, dy; + IssmDouble dist, alpha, ang; + IssmDouble N_azim, E_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; + + /*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]; + } + + 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); + xDelete(U_elastic); xDelete(N_elastic); xDelete(E_elastic); + + return; +} +/*}}}*/ +void Tria::EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/ /*diverse:*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 21271) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 21272) @@ -145,5 +145,6 @@ #endif #ifdef _HAVE_ESA_ - void EsaGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea); + void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,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 #ifdef _HAVE_SEALEVELRISE_ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21271) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21272) @@ -2300,5 +2300,39 @@ #endif #ifdef _HAVE_ESA_ -void FemModel::EsaGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/ +void FemModel::EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* xx, IssmDouble* yy){/*{{{*/ + + int ns,nsmax; + + /*Go through elements, and add contribution from each element to the deflection vector wg:*/ + ns = elements->Size(); + + /*Figure out max of ns: */ + ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + + /*Call the esa geodetic core: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + element->EsaGeodetic2D(pUp,pNorth,pEast,xx,yy); + } + if(i%100==0){ + pUp->Assemble(); + pNorth->Assemble(); + pEast->Assemble(); + } + } + + /*One last time: */ + pUp->Assemble(); + pNorth->Assemble(); + pEast->Assemble(); + + /*Free ressources:*/ + xDelete(xx); + xDelete(yy); +} +/*}}}*/ +void FemModel::EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/ IssmDouble eartharea=0; @@ -2326,5 +2360,5 @@ if(i(elements->GetObjectByOffset(i)); - element->EsaGeodetic(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz,eartharea); + element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz,eartharea); } if(i%100==0){ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21271) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21272) @@ -114,5 +114,6 @@ #endif #ifdef _HAVE_ESA_ - void EsaGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); + void EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, 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 #ifdef _HAVE_SEALEVELRISE_ Index: /issm/trunk-jpl/src/c/cores/esa_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/esa_core.cpp (revision 21271) +++ /issm/trunk-jpl/src/c/cores/esa_core.cpp (revision 21272) @@ -17,4 +17,5 @@ bool save_results,isesa,iscoupler; int configuration_type; + int domaintype; int solution_type; int numoutputs = 0; @@ -32,4 +33,5 @@ /*Recover some parameters: */ + femmodel->parameters->FindParam(&domaintype,DomainTypeEnum); femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); @@ -38,8 +40,6 @@ femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); - /*first, recover lat,long and radius vectors from vertices: */ + /* recover coordinates of vertices: */ VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); - - /*recover x,y,z vectors from vertices: */ VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); @@ -78,5 +78,10 @@ /*call the geodetic main modlule:*/ - femmodel->EsaGeodetic(U_radial,U_north,U_east,latitude,longitude,radius,xx,yy,zz); + if(domaintype==Domain3DsurfaceEnum){ + femmodel->EsaGeodetic3D(U_radial,U_north,U_east,latitude,longitude,radius,xx,yy,zz); + } + if(domaintype==Domain2DhorizontalEnum){ + femmodel->EsaGeodetic2D(U_radial,U_north,U_east,xx,yy); + } /*get results into elements:*/