Changeset 26149


Ignore:
Timestamp:
03/24/21 12:11:09 (4 years ago)
Author:
Eric.Larour
Message:

CHG: ElementCoordinatesx now also returns areas.
BarystaticContributions bsl components have right sign .
Tria: idem + debugged GeometryOptim G which was not synced with the variable area.

Location:
issm/trunk-jpl/src/c
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r26147 r26149  
    148148
    149149        ice->Sum(&sumice); hydro->Sum(&sumhydro); ocean->Sum(&sumocean);
    150         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcEnum,-this->Total()/oceanarea/rho_water,step,time));
    151         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcIceEnum,-sumice/oceanarea/rho_water,step,time));
    152         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcHydroEnum,-sumice/oceanarea/rho_water,step,time));
    153         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcOceanEnum,-sumocean/oceanarea/rho_water,step,time));
     150        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcEnum,this->Total()/oceanarea/rho_water,step,time));
     151        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcIceEnum,sumice/oceanarea/rho_water,step,time));
     152        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcHydroEnum,sumice/oceanarea/rho_water,step,time));
     153        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcOceanEnum,sumocean/oceanarea/rho_water,step,time));
    154154
    155155        cumice->Sum(&sumice); cumhydro->Sum(&sumhydro); cumocean->Sum(&sumocean);

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

    r26121 r26149  
    239239                virtual void       ElementResponse(IssmDouble* presponse,int response_enum)=0;
    240240                virtual void       ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
    241                 virtual void       ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze,bool spherical=false)=0;
     241                virtual void       ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze, Vector<IssmDouble>* vareae, bool spherical=false)=0;
    242242                virtual int        FiniteElement(void)=0;
    243243                virtual IssmDouble FloatingArea(bool scaled)=0;
    … …  
    392392                virtual void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0;
    393393
    394                 virtual void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze)=0;
     394                virtual void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
    395395                virtual void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0;
    396396                virtual void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks)=0;

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

    r26121 r26149  
    7070                void           CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time,int averaging_method);
    7171                void           ElementResponse(IssmDouble* presponse,int response_enum);
    72                 void           ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze,bool spherical=false){_error_("not implemented yet");};
     72                void           ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze,Vector<IssmDouble>* vareae,bool spherical=false){_error_("not implemented yet");};
    7373                void           ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    7474                int            FiniteElement(void);
    … …  
    229229                void           GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    230230
    231                 void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
     231                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    232232                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    233233                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};

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

    r26121 r26149  
    5454                void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum,int control_interp){_error_("not implemented yet");};
    5555                void        ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};
    56                 void        ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze,bool spherical=false){_error_("not implemented yet");};
     56                void        ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze,Vector<IssmDouble>* vareae,bool spherical=false){_error_("not implemented yet");};
    5757                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
    5858                int         FiniteElement(void);
    … …  
    184184                void        GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    185185
    186                 void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
     186                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    187187                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    188188                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};

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

    r26121 r26149  
    5353                IssmDouble  DragCoefficientAbsGradient(void){_error_("not implemented yet");};
    5454                void        ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};
    55                 void        ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze,bool spherical=false){_error_("not implemented yet");};
     55                void        ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze,Vector<IssmDouble>* vareae,bool spherical=false){_error_("not implemented yet");};
    5656                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    5757                void        FaceOnBaseIndices(int* pindex1,int* pindex2,int* pindex3);
    … …  
    190190                IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    191191                void        GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    192                 void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
     192                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    193193                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    194194                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};

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

    r26126 r26149  
    12861286}
    12871287/*}}}*/
    1288 void       Tria::ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze, bool spherical){/*{{{*/
     1288void       Tria::ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze, Vector<IssmDouble>* vareae, bool spherical){/*{{{*/
    12891289
    12901290        /*Look for x,y,z coordinates:*/
    … …  
    12981298
    12991299        if(!spherical){
     1300                IssmDouble area;
    13001301                vxe->SetValue(this->sid,xe,INS_VAL);
    13011302                vye->SetValue(this->sid,ye,INS_VAL);
    13021303                vze->SetValue(this->sid,ze,INS_VAL);
     1304                area=this->GetAreaSpherical();
     1305                vareae->SetValue(this->sid,area,INS_VAL);
     1306               
     1307                /*in addition, put in in the inputs:*/
     1308                this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
    13031309        }
    13041310        else _error_("spherical coordinates not supported yet!");
    … …  
    62396245
    62406246//new code
    6241 void       Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
     6247void       Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/
    62426248
    62436249        /*diverse:*/
    62446250        int nel;
    6245         bool spherical=true;
    6246         IssmDouble llr_list[NUMVERTICES][3];
    6247         IssmDouble xyz_list[NUMVERTICES][3];
    62486251        IssmDouble area,planetarea,planetradius;
     6252        IssmDouble constant,ratioe;
    62496253        IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
    62506254        IssmDouble rho_earth;
    62516255        IssmDouble lati,longi;
    6252         IssmDouble constant;
    62536256        IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim;
    62546257        int sidlist[NUMVERTICES];
    … …  
    63216324        }
    63226325
    6323         /*compute area and add to inputs:*/
    6324         area=GetAreaSpherical();
    6325         this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
    6326         this->AddInput(SealevelAreaEnum,&area,P0Enum);
     6326        /*early return:*/
    63276327        if(!computerigid)return;
    63286328
    … …  
    63536353        }
    63546354
    6355         /*retrieve coordinates: */
    6356         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6357 
    63586355        /*compute centroids of all elements:*/
    63596356        for(int i=0;i<nel;i++){
    … …  
    63626359        }
    63636360       
    6364         constant=3/rho_earth*area/planetarea;
     6361       
     6362        constant=3/rho_earth/planetarea;
    63656363
    63666364        /*Recover vertex absolute id: */
    … …  
    63686366
    63696367        for(int e=0;e<nel;e++){
     6368                ratioe=constant*areae[e];
    63706369                for (int i=0;i<3;i++){
    63716370                        IssmDouble alpha;
    … …  
    63886387                                G[i*nel+e]+=G_elastic_precomputed[index];
    63896388                        }
    6390                         G[i*nel+e]=G[i*nel+e]*constant;
     6389                        G[i*nel+e]=G[i*nel+e]*ratioe;
    63916390
    63926391                        /*Elastic components:*/
    63936392                        if(computeelastic){
    6394                                 GU[i*nel+e] =  constant * U_elastic_precomputed[index];
     6393                                GU[i*nel+e] =  ratioe * U_elastic_precomputed[index];
    63956394                                if(horiz){
    63966395                                        /*Compute azimuths, both north and east components: */
    … …  
    64066405                                        E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    64076406
    6408                                         GN[i*nel+e] = constant*H_elastic_precomputed[index]*N_azim;
    6409                                         GE[i*nel+e] = constant*H_elastic_precomputed[index]*E_azim;
     6407                                        GN[i*nel+e] = ratioe*H_elastic_precomputed[index]*N_azim;
     6408                                        GE[i*nel+e] = ratioe*H_elastic_precomputed[index]*E_azim;
    64106409                                }
    64116410                        }
    … …  
    64396438                }
    64406439        }
    6441         this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
    6442         this->AddInput(SealevelAreaEnum,&area,P0Enum);
    64436440
    64446441        /*Free allocations:*/
    … …  
    64726469void       Tria::SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){ /*{{{*/
    64736470
     6471        //Compute sea level barystatic loads (ice, water column or bottom pressure, see Farrel and Clarke, Equ. 4)
     6472
    64746473        /*diverse:*/
    64756474        IssmDouble area;
    64766475        IssmDouble phi_ice=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic
    64776476        IssmDouble phi_water=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic
    6478         IssmDouble I,W,BP;  //change in ice thickness, water column or bottom pressure (Farrel and Clarke, Equ. 4)
    6479         bool notfullygrounded=false;
     6477        IssmDouble phi_bp=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic
     6478        IssmDouble I=0; //Do not change this!
     6479        IssmDouble W=0; //Do not change this!
     6480        IssmDouble BP=0; //Do not change this!
     6481
     6482                bool notfullygrounded=false;
    64806483        bool computerigid= false;
    64816484        int  glfraction=1;
    … …  
    65296532        if(!isice  && !ishydro){
    65306533                if(!masks->isoceanin[this->lid]){
     6534                        constant=0;
     6535                        #ifdef _ISSM_DEBUG_
     6536                        this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
     6537                        #endif
    65316538                        return;
    65326539                }
    … …  
    66036610
    66046611                /*Compute barystatic contribution in kg:*/
    6605                 bslcice = rho_ice*area*I;
     6612                bslcice = -rho_ice*phi_ice*area*I;
    66066613                _assert_(!xIsNan<IssmDouble>(bslcice));
    66076614
    … …  
    66126619        /*Deal with water loads if we are on ground:*/
    66136620        if(!masks->isfullyfloating[this->lid]){
     6621
     6622                constant+=10; //1 for water.
     6623                #ifdef _ISSM_DEBUG_
     6624                this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
     6625                #endif
    66146626
    66156627                /*Retrieve water height at vertices: */
    … …  
    66196631
    66206632                /*Compute barystatic component in kg:*/
    6621                 bslchydro = rho_freshwater*area*phi_water*W;
     6633                bslchydro = -rho_freshwater*area*phi_water*W;
    66226634                _assert_(!xIsNan<IssmDouble>(bslchydro));
    66236635
    … …  
    66286640        /*Deal with ocean bottom pressures:*/
    66296641        if(masks->isoceanin[this->lid]){
     6642
     6643                constant+=1; //1 for bottom pressure.
     6644                #ifdef _ISSM_DEBUG_
     6645                this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
     6646                #endif
    66306647
    66316648                /*Retrieve bottom pressure change and average over the element: */
    … …  
    66356652               
    66366653                /*Compute barystatic component in kg:*/
    6637                 bslcbp = rho_water*area*BP;
     6654                bslcbp = -rho_water*area*phi_bp*BP;
    66386655
    66396656                /*convert from m to kg/m^2:*/
    6640                 BP=BP*rho_water;
     6657                BP=BP*rho_water*phi_bp;
    66416658        }
    66426659

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

    r26121 r26149  
    7070                void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum,int control_interp);
    7171                void        CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum);
    72                 void        ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze,bool spherical=false);
     72                void        ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze, Vector<IssmDouble>* vareae, bool spherical=false);
    7373                int         EdgeOnBaseIndex();
    7474                void        EdgeOnBaseIndices(int* pindex1,int* pindex);
    … …  
    174174                void       SetSealevelMasks(SealevelMasks* masks);
    175175               
    176                 void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
     176                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae);
    177177                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks);
    178178                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks);

  • issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp

    r26146 r26149  
    608608        IssmDouble *yye    = NULL;
    609609        IssmDouble *zze    = NULL;
     610        IssmDouble* areae  = NULL;
    610611
    611612        int  horiz;
    … …  
    634635        /*first, recover lat,long and radius vectors from vertices: */
    635636        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
    636         ElementCoordinatesx(&xxe,&yye,&zze,femmodel->elements);
    637637        if(horiz) VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
    638 
     638       
     639        /*first, recover x,y,z and areas from elements: */
     640        ElementCoordinatesx(&xxe,&yye,&zze,&areae,femmodel->elements);
    639641
    640642        /*Run sealevel geometry routine in elements:*/
    641643        for(Object* & object : femmodel->elements->objects){
    642644                Element*   element=xDynamicCast<Element*>(object);
    643                 if(optim) element->SealevelchangeGeometryOptim(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze);
     645                if(optim) element->SealevelchangeGeometryOptim(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze,areae);
    644646                else      element->SealevelchangeGeometry(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze);
    645647        }
    … …  
    654656        xDelete<IssmDouble>(yye);
    655657        xDelete<IssmDouble>(zze);
     658        xDelete<IssmDouble>(areae);
    656659        xDelete<IssmDouble>(latitude);
    657660        xDelete<IssmDouble>(longitude);
    … …  
    12511254} /*}}}*/
    12521255IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea){ /*{{{*/
     1256
    12531257        IssmDouble sealevelloadsaverage;       
     1258
    12541259        Vector<IssmDouble>* sealevelloadsvolume=sealevelloads->Duplicate();
    12551260        sealevelloadsvolume->PointwiseMult(sealevelloads,oceanareas);
    12561261        sealevelloadsvolume->Sum(&sealevelloadsaverage);
    12571262        delete sealevelloadsvolume;
    1258         return sealevelloadsaverage/oceanarea;
     1263       
     1264        return sealevelloadsaverage/oceanarea;-
    12591265} /*}}}*/
    12601266void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealevelloads){ /*{{{*/

  • issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.cpp

    r26099 r26149  
    99#include "../../toolkits/toolkits.h"
    1010
    11 void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze,Elements* elements,bool spherical) {
     11void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze, IssmDouble** pareae, Elements* elements,bool spherical) {
    1212
    1313        /*figure out how many vertices we have: */
    … …  
    1717        Vector<IssmDouble>* vye=new Vector<IssmDouble>(numberofelements);
    1818        Vector<IssmDouble>* vze=new Vector<IssmDouble>(numberofelements);
     19        Vector<IssmDouble>* vareae=new Vector<IssmDouble>(numberofelements);
    1920
    2021        /*march through our elements: */
    2122        for(Object* & object : elements->objects){
    2223                Element* element=(Element*)object;
    23                 element->ElementCoordinates(vxe,vye,vze,spherical);
     24                element->ElementCoordinates(vxe,vye,vze,vareae,spherical);
    2425        }
    2526
    … …  
    2829        vye->Assemble();
    2930        vze->Assemble();
     31        vareae->Assemble();
    3032
    3133        /*serialize: */
    … …  
    3335        IssmDouble* ye=vye->ToMPISerial();
    3436        IssmDouble* ze=vze->ToMPISerial();
     37        IssmDouble* areae=vareae->ToMPISerial();
    3538
    3639        /*Free ressources: */
    … …  
    3841        delete vye;
    3942        delete vze;
     43        delete vareae;
    4044
    4145        /*output: */
    … …  
    4650        if(pze) *pze=ze;
    4751        else xDelete<IssmDouble>(ze);
     52        if(pareae) *pareae=areae;
     53        else xDelete<IssmDouble>(areae);
     54
    4855}

  • issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.h

    r26099 r26149  
    88
    99/* local prototypes: */
    10 void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze,Elements* elements,bool spherical=false);
     10void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze,IssmDouble** pareae, Elements* elements,bool spherical=false);
    1111
    1212#endif  /* _ELEMENT_COORDINATESX_H */
Note: See TracChangeset for help on using the changeset viewer.