Changeset 24946


Ignore:
Timestamp:
06/01/20 14:22:09 (5 years ago)
Author:
Eric.Larour
Message:

CHG: further optimization in the geodetic module, using the geometry module.

Location:
issm/trunk-jpl/src/c/classes
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r24940 r24946  
    380380                virtual void          SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0;
    381381                virtual void          SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea)=0;
    382                 virtual void          SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0;
     382                virtual void          SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz)=0;
    383383                virtual void          SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0;
    384                 virtual void          SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz)=0;
     384                virtual void          SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz)=0;
    385385                #endif
    386386

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r24940 r24946  
    215215                void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    216216                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
    217                 void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
     217                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");};
    218218                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");};
    219219                void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
    220                 void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");};
     220                void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
    221221                #endif
    222222

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r24940 r24946  
    174174                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
    175175                void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    176                 void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
     176                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");};
    177177                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");};
    178178                void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
    179                 void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");};
     179                void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
    180180                IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    181181#endif

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r24940 r24946  
    180180                void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    181181                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
    182                 void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
     182                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");};
    183183                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");};
    184184                void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
    185                 void    SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");};
     185                void    SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
    186186                IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    187187#endif

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r24944 r24946  
    55785578}
    55795579/*}}}*/
    5580 void    Tria::SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){ /*{{{*/
    5581 
     5580void    Tria::SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){ /*{{{*/
    55825581        /*diverse:*/
    55835582        int gsize;
    55845583        bool spherical=true;
    55855584        IssmDouble llr_list[NUMVERTICES][3];
     5585        IssmDouble xyz_list[NUMVERTICES][3];
    55865586        IssmDouble area,eartharea;
    55875587        IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
    … …  
    55905590        IssmDouble lati,longi,ri;
    55915591        IssmDouble constant;
     5592        IssmDouble x_element,y_element,z_element,x,y,z,dx,dy,dz,N_azim,E_azim;
    55925593
    55935594        IssmDouble* G=NULL;
     5595        IssmDouble* GU=NULL;
     5596        IssmDouble* GN=NULL;
     5597        IssmDouble* GE=NULL;
    55945598        IssmDouble* G_elastic_precomputed=NULL;
    55955599        IssmDouble* G_rigid_precomputed=NULL;
     5600        IssmDouble* U_elastic_precomputed=NULL;
     5601        IssmDouble* H_elastic_precomputed=NULL;
    55965602
    55975603        /*elastic green function:*/
    55985604        IssmDouble* indices=NULL;
     5605        int index;
    55995606        int         M;
    56005607
    … …  
    56025609        bool computerigid = true;
    56035610        bool computeelastic = true;
     5611        int  horiz;
    56045612
    56055613        /*recover parameters: */
    … …  
    56095617        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
    56105618        this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
     5619        this->parameters->FindParam(&horiz,SealevelriseHorizEnum);
    56115620
    56125621        /*recover precomputed green function kernels:*/
    … …  
    56165625        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
    56175626        parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
    5618        
     5627
     5628        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter);
     5629        parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M);
     5630
     5631        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
     5632        parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M);
     5633
    56195634        /*allocate indices:*/
    56205635        indices=xNew<IssmDouble>(gsize);
    56215636        G=xNewZeroInit<IssmDouble>(gsize);
     5637        GU=xNewZeroInit<IssmDouble>(gsize);
     5638        if(horiz){
     5639                GN=xNewZeroInit<IssmDouble>(gsize);
     5640                GE=xNewZeroInit<IssmDouble>(gsize);
     5641        }
    56225642       
    56235643        /*compute area:*/
    56245644        area=GetAreaSpherical();
     5645
     5646        /*figure out gravity center of our element (Cartesian): */
     5647        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     5648        x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
     5649        y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
     5650        z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
    56255651
    56265652        /* Where is the centroid of this element?:{{{*/
    … …  
    56745700                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    56755701                indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
    5676 
     5702                index=reCast<int,IssmDouble>(indices[i]);
     5703               
    56775704                /*Rigid earth gravitational perturbation: */
    56785705                if(computerigid){
    5679                         G[i]+=G_rigid_precomputed[reCast<int,IssmDouble>(indices[i])];
     5706                        G[i]+=G_rigid_precomputed[index];
    56805707                }
    56815708                if(computeelastic){
    5682                         G[i]+=G_elastic_precomputed[reCast<int,IssmDouble>(indices[i])];
     5709                        G[i]+=G_elastic_precomputed[index];
    56835710                }
    56845711                G[i]=G[i]*constant;
     5712
     5713                /*Elastic components:*/
     5714                GU[i] =  constant * U_elastic_precomputed[index];
     5715                if(horiz){
     5716                        /*Compute azimuths, both north and east components: */
     5717                        x = xx[i]; y = yy[i]; z = zz[i];
     5718                        if(latitude[i]==90){
     5719                                x=1e-12; y=1e-12;
     5720                        }
     5721                        if(latitude[i]==-90){
     5722                                x=1e-12; y=1e-12;
     5723                        }
     5724                        dx = x_element-x; dy = y_element-y; dz = z_element-z;
     5725                        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);
     5726                        E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
     5727
     5728                        GN[i] = constant*H_elastic_precomputed[index]*N_azim;
     5729                        GE[i] = constant*H_elastic_precomputed[index]*E_azim;
     5730                }
    56855731        }
    56865732
    … …  
    56885734    this->inputs2->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
    56895735    this->inputs2->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
     5736    this->inputs2->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
     5737    if(horiz){
     5738                this->inputs2->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize);
     5739                this->inputs2->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize);
     5740        }
    56905741        this->inputs2->SetDoubleInput(AreaEnum,this->lid,area);
    56915742
    … …  
    56935744        xDelete(indices);
    56945745        xDelete(G);
     5746        xDelete(GU);
     5747        if(horiz){
     5748                xDelete(GN);
     5749                xDelete(GE);
     5750        }
    56955751
    56965752        return;
    … …  
    60346090}
    60356091/*}}}*/
    6036 void    Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){ /*{{{*/
     6092void    Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){ /*{{{*/
    60376093
    60386094        /*diverse:*/
    6039         int gsize,dummy;
    6040         bool spherical=true;
    6041         IssmDouble llr_list[NUMVERTICES][3];
    6042         IssmDouble xyz_list[NUMVERTICES][3];
    6043         IssmDouble area,eartharea;
     6095        int gsize;
    60446096        IssmDouble I, S;                //change in relative ice thickness and sea level
    6045         IssmDouble late,longe,re;
    6046         IssmDouble lati,longi,ri;
    6047         IssmDouble rho_ice,rho_water,rho_earth;
    6048         IssmDouble minlong=400;
    6049         IssmDouble maxlong=-20;
     6097        IssmDouble rho_ice,rho_water;
     6098        int horiz;
    60506099
    60516100        /*precomputed elastic green functions:*/
    6052         IssmDouble* U_elastic_precomputed = NULL;
    6053         IssmDouble* H_elastic_precomputed = NULL;
    6054         int         M;
    6055         IssmDouble* indices=NULL;
    6056         int         index;
    6057 
    6058         /*computation of Green functions:*/
    6059         IssmDouble* U_elastic= NULL;
    6060         IssmDouble* N_elastic= NULL;
    6061         IssmDouble* E_elastic= NULL;
    6062         DoubleVecParam* U_parameter = NULL;
    6063         DoubleVecParam* H_parameter = NULL;
    6064         IssmDouble* U_values=NULL;
    6065         IssmDouble* N_values=NULL;
    6066         IssmDouble* E_values=NULL;
    6067 
    6068         /*optimization:*/
    6069         bool store_green_functions=false;
     6101        IssmDouble* GU=NULL;
     6102        IssmDouble* GN=NULL;
     6103        IssmDouble* GE=NULL;
    60706104
    60716105        /*computational flags:*/
    6072         bool computerigid = true;
    60736106        bool computeelastic= true;
    60746107
    … …  
    60796112        if(masks->isfullyfloating[this->lid])return;
    60806113
    6081         /*recover computational flags: */
    6082         this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
     6114        /*recover parameters:*/
    60836115        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
     6116        this->parameters->FindParam(&horiz,SealevelriseHorizEnum);
    60846117
    60856118        /*early return if elastic not requested:*/
    … …  
    60896122        rho_ice=FindParam(MaterialsRhoIceEnum);
    60906123        rho_water=FindParam(MaterialsRhoSeawaterEnum);
    6091         rho_earth=FindParam(MaterialsEarthDensityEnum);
    6092 
    6093         /*how many dofs are we working with here? */
    6094         this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
    6095 
    6096         /*recover earth area: */
    6097         this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
    6098 
    6099         /*retrieve indices:*/
    6100         if(computerigid){this->inputs2->GetArrayPtr(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);}
    6101 
    6102         /*Get area of element: precomputed in the sealevelrise_core_geometry:*/
    6103         this->GetInput2Value(&area,AreaEnum);
    6104 
    6105         /*figure out gravity center of our element (Cartesian): */
    6106         IssmDouble x_element, y_element, z_element;
    6107         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6108         x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
    6109         y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
    6110         z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
    61116124
    61126125        /*recover elastic Green's functions for displacement:*/
    6113         U_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(U_parameter);
    6114         U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M);
     6126        this->inputs2->GetArrayPtr(SealevelriseGUEnum,this->lid,&GU,&gsize);
    61156127        if(horiz){
    6116                 H_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(H_parameter);
    6117                 H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M);
     6128                this->inputs2->GetArrayPtr(SealevelriseGUEnum,this->lid,&GU,&gsize);
     6129                this->inputs2->GetArrayPtr(SealevelriseGUEnum,this->lid,&GU,&gsize);
    61186130        }
    61196131
    … …  
    61266138        deltathickness_input->GetInputAverage(&I);
    61276139
    6128         /*initialize: */
    6129         U_elastic=xNewZeroInit<IssmDouble>(gsize);
    6130         if(horiz){
    6131                 N_elastic=xNewZeroInit<IssmDouble>(gsize);
    6132                 E_elastic=xNewZeroInit<IssmDouble>(gsize);
    6133         }
    6134 
    6135         //int* indices=xNew<int>(gsize);
    6136         //U_values=xNewZeroInit<IssmDouble>(gsize);
    6137         if(horiz){
    6138                 //N_values=xNewZeroInit<IssmDouble>(gsize);
    6139                 //E_values=xNewZeroInit<IssmDouble>(gsize);
    6140         }
    6141         IssmDouble alpha;
    6142         IssmDouble delPhi,delLambda;
    6143         IssmDouble dx, dy, dz, x, y, z;
    6144         IssmDouble N_azim, E_azim;
    6145 
    6146         /*we are going to use the result of these masks a lot: */
    6147         bool isiceonlyinelement=masks->isiceonly[this->lid];
    6148         bool isoceaninelement=masks->isoceanin[this->lid];
    6149 
    6150         for(int i=0;i<gsize;i++){
    6151 
    6152                 /*Compute alpha angle between centroid and current vertex: */
    6153                 /*lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    6154                 delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
    6155                 alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));*/
    6156 
    6157                 index=reCast<int,IssmDouble>(indices[i]);
    6158 
    6159                 /*Compute azimuths, both north and east components: */
    6160                 x = xx[i]; y = yy[i]; z = zz[i];
    6161                 if(latitude[i]==90){
    6162                         x=1e-12; y=1e-12;
    6163                 }
    6164                 if(latitude[i]==-90){
    6165                         x=1e-12; y=1e-12;
    6166                 }
    6167                 dx = x_element-x; dy = y_element-y; dz = z_element-z;
    6168                 if(horiz){
    6169                         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);
    6170                         E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    6171                 }
    6172 
    6173                 /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
    6174                 U_elastic[i] += U_elastic_precomputed[index];
    6175                 if(horiz){
    6176                         N_elastic[i] += H_elastic_precomputed[index]*N_azim;
    6177                         E_elastic[i] += H_elastic_precomputed[index]*E_azim;
    6178                 }
    6179 
    6180                 /*Add all components to the pUp solution vectors:*/
    6181                 if(isiceonlyinelement){
    6182                         Up[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i];
     6140        if (masks->isiceonly[this->lid]){
     6141                for(int i=0;i<gsize;i++){
     6142                        Up[i]+=rho_ice*I*GU[i];
    61836143                        if(horiz){
    6184                                 North[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i];
    6185                                 East[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i];
     6144                                North[i]+=rho_ice*I*GN[i];
     6145                                East[i]+=rho_ice*I*GE[i];
    61866146                        }
    61876147                }
    6188                 else if(isoceaninelement){
    6189                         Up[i]+=3*rho_water/rho_earth*area/eartharea*S*U_elastic[i];
     6148        }
     6149        else if(masks->isoceanin[this->lid]){
     6150                for(int i=0;i<gsize;i++){
     6151                        Up[i]+=rho_water*S*GU[i];
    61906152                        if(horiz){
    6191                                 North[i]+=3*rho_water/rho_earth*area/eartharea*S*N_elastic[i];
    6192                                 East[i]+=3*rho_water/rho_earth*area/eartharea*S*E_elastic[i];
     6153                                North[i]+=rho_water*S*GN[i];
     6154                                East[i]+=rho_water*S*GE[i];
    61936155                        }
    61946156                }
    6195         }
    6196 
    6197         /*free ressources:*/
    6198         xDelete<IssmDouble>(U_values);
    6199         xDelete<IssmDouble>(U_elastic);
    6200         if(horiz){
    6201                 xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
    6202                 xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);
    62036157        }
    62046158

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r24940 r24946  
    166166                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
    167167                void    SetSealevelMasks(SealevelMasks* masks);
    168                 void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius);
     168                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
    169169                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea);
    170170                void    SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea);
    171171                void    SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea);
    172172                void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius);
    173                 void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz);
     173                void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz);
    174174                #endif
    175175                /*}}}*/

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r24944 r24946  
    46734673#endif
    46744674#ifdef _HAVE_SEALEVELRISE_
    4675 void FemModel::SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/
     4675void FemModel::SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz) { /*{{{*/
    46764676
    46774677        /*Run sealevelrie geometry routine in elements:*/
    46784678        for(int i=0;i<elements->Size();i++){
    46794679                Element*   element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4680                 element->SealevelriseGeometry(latitude,longitude,radius);
     4680                element->SealevelriseGeometry(latitude,longitude,radius,xx,yy,zz);
    46814681        }
    46824682
    … …  
    48544854}
    48554855/*}}}*/
    4856 void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
     4856void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop){/*{{{*/
    48574857
    48584858        /*serialized vectors:*/
    … …  
    48634863        IssmDouble* East  = NULL;
    48644864        int* indices = NULL;
    4865         int         gsize;
    4866 
     4865        int  gsize;
     4866        int  horiz;
     4867       
     4868        /*retrieve parameters:*/
     4869        this->parameters->FindParam(&horiz,SealevelriseHorizEnum);
    48674870
    48684871        /*Serialize vectors from previous iteration:*/
    … …  
    48804883        for(int i=0;i<elements->Size();i++){
    48814884                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4882                 element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz,horiz);
     4885                element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz);
    48834886        }
    48844887

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r24940 r24946  
    164164                #ifdef _HAVE_SEALEVELRISE_
    165165                void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
    166                 void SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
     166                void SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
    167167                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old,  SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
    168168                void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
    169                 void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
     169                void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop);
    170170                IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg,SealevelMasks* masks, IssmDouble oceanarea);
    171171                #endif
Note: See TracChangeset for help on using the changeset viewer.