Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24914) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24915) @@ -377,7 +377,7 @@ virtual IssmDouble OceanArea(void)=0; virtual void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea)=0; - virtual void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0; - virtual void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0; - virtual void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz)=0; + virtual void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0; + virtual void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0; + virtual void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24914) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24915) @@ -215,7 +215,7 @@ IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; + void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_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 24914) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24915) @@ -173,7 +173,7 @@ #ifdef _HAVE_SEALEVELRISE_ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; + void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; IssmDouble OceanArea(void){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24914) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24915) @@ -179,7 +179,7 @@ #ifdef _HAVE_SEALEVELRISE_ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; + void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; IssmDouble OceanArea(void){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24914) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24915) @@ -5564,5 +5564,5 @@ return; }/*}}}*/ -void Tria::SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ /*Computational flags:*/ @@ -5575,14 +5575,14 @@ /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested *bottom pressure fingerprints:*/ - if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(pSgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); + if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); } //if(!IsIceInElement()){ /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/ - this->SealevelriseEustaticIce(pSgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); + this->SealevelriseEustaticIce(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); //} } /*}}}*/ -void Tria::SealevelriseEustaticIce(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ /*diverse:*/ @@ -5744,10 +5744,7 @@ if(computeelastic | computerigid){ - int* indices=xNew(gsize); - IssmDouble* values=xNew(gsize); IssmDouble alpha; IssmDouble delPhi,delLambda; for(int i=0;iSetValues(gsize,indices,values,ADD_VAL); - - /*free ressources:*/ - xDelete(values); - xDelete(indices); + /*Add all components to the Sgi or Sgo solution vectors:*/ + Sgi[i]+=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic); + + } } @@ -5786,5 +5779,5 @@ } /*}}}*/ -void Tria::SealevelriseEustaticBottomPressure(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ /*diverse:*/ @@ -5898,10 +5891,7 @@ if(computeelastic | computerigid){ - int* indices=xNew(gsize); - IssmDouble* values=xNew(gsize); IssmDouble alpha; IssmDouble delPhi,delLambda; for(int i=0;iSetValues(gsize,indices,values,ADD_VAL); - - /*free ressources:*/ - xDelete(values); - xDelete(indices); + Sgi[i]+=3*rho_water/rho_earth*area/eartharea*I*(G_rigid+G_elastic); + } } @@ -5940,5 +5925,5 @@ } /*}}}*/ -void Tria::SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/ /*diverse:*/ @@ -6046,12 +6031,8 @@ if(computerigid) G_rigid=xNewZeroInit(gsize); - int* indices=xNew(gsize); - IssmDouble* values=xNewZeroInit(gsize); IssmDouble alpha; IssmDouble delPhi,delLambda; for(int i=0;i(alpha/PI*(M-1)); G_elastic[i] += G_elastic_precomputed[index]; - values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i]; - } - } - - pSgo->SetValues(gsize,indices,values,ADD_VAL); - - /*free ressources:*/ - xDelete(values); - xDelete(indices); + //values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i]; + Sgo[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i]; + } + } + /*Free ressources:*/ @@ -6088,5 +6066,5 @@ } /*}}}*/ -void Tria::SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){ /*{{{*/ +void Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){ /*{{{*/ /*diverse:*/ @@ -6217,9 +6195,9 @@ } - int* indices=xNew(gsize); - U_values=xNewZeroInit(gsize); + //int* indices=xNew(gsize); + //U_values=xNewZeroInit(gsize); if(horiz){ - N_values=xNewZeroInit(gsize); - E_values=xNewZeroInit(gsize); + //N_values=xNewZeroInit(gsize); + //E_values=xNewZeroInit(gsize); } IssmDouble alpha; @@ -6230,5 +6208,5 @@ for(int i=0;iSetValues(gsize,indices,U_values,ADD_VAL); - if(horiz){ - pNorth->SetValues(gsize,indices,N_values,ADD_VAL); - pEast->SetValues(gsize,indices,E_values,ADD_VAL); - } /*free ressources:*/ - xDelete(indices); xDelete(U_values); xDelete(U_elastic); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24914) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24915) @@ -166,9 +166,9 @@ IssmDouble OceanAverage(IssmDouble* Sg); void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea); - void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); - void SealevelriseEustaticIce(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); - void SealevelriseEustaticBottomPressure(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); - void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea); - void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz); + void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); + void SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); + void SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); + void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea); + void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz); #endif /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24914) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24915) @@ -4694,11 +4694,16 @@ IssmDouble eartharea = 0.; IssmDouble eartharea_cpu = 0.; - int ns,nsmax; - - /*Go through elements, and add contribution from each element to the deflection vector wg:*/ - ns = elements->Size(); + IssmDouble* RSLgi = NULL; + int* indices = NULL; + int gsize; + + /*Initialize temporary vector that will be used to sum eustatic components on all local elements, prior + * to assembly:*/ + gsize = this->nodes->NumberOfDofs(GsetEnum); + RSLgi=xNewZeroInit(gsize); + indices=xNew(gsize); for (int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); oceanarea_cpu += element->OceanArea(); @@ -4712,25 +4717,14 @@ ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - /*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 sea level rise core: */ - for(int i=0;i(elements->GetObjectByOffset(i)); - element->SealevelriseEustatic(pRSLgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea); - eustatic_cpu+=eustatic_cpu_e; - } - if(i%loop==0)pRSLgi->Assemble(); - } - if(VerboseConvergence())_printf0_("\n"); - - /*One last time: */ + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea); + eustatic_cpu+=eustatic_cpu_e; + } + /*Plug values once and assemble: */ + pRSLgi->SetValues(gsize,indices,RSLgi,ADD_VAL); pRSLgi->Assemble(); - + /*Sum all eustatic components from all cpus:*/ ISSM_MPI_Reduce (&eustatic_cpu,&eustatic,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); @@ -4738,4 +4732,8 @@ _assert_(!xIsNan(eustatic)); + /*Free ressources:*/ + xDelete(indices); + xDelete(RSLgi); + /*Assign output pointers:*/ *peustatic=eustatic; @@ -4751,14 +4749,19 @@ IssmDouble eartharea_cpu=0; - int ns,nsmax; + IssmDouble* RSLgo = NULL; + int* indices = NULL; + int gsize; + + /*Initialize temporary vector that will be used to sum eustatic components on all local elements, prior + * to assembly:*/ + gsize = this->nodes->NumberOfDofs(GsetEnum); + RSLgo=xNewZeroInit(gsize); + indices=xNew(gsize); for (int i=0;iToMPISerial(); - /*Go through elements, and add contribution from each element to the deflection vector wg:*/ - ns = elements->Size(); - /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ - for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); eartharea_cpu += element->GetAreaSpherical(); @@ -4768,21 +4771,17 @@ ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - /*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 sea level rise core: */ - for(int i=0;i(elements->GetObjectByOffset(i)); - element->SealevelriseNonEustatic(pRSLgo,RSLg_old,latitude,longitude,radius,eartharea); - } - if(i%loop==0)pRSLgo->Assemble(); - } - if(verboseconvolution)if(VerboseConvergence())_printf0_("\n"); + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->SealevelriseNonEustatic(RSLgo,RSLg_old,latitude,longitude,radius,eartharea); + } + pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL); + pRSLgo->Assemble(); /*Free ressources:*/ + xDelete(indices); + xDelete(RSLgo); xDelete(RSLg_old); + } /*}}}*/ @@ -4881,14 +4880,24 @@ IssmDouble eartharea_cpu=0; - int ns,nsmax; + IssmDouble* Up = NULL; + IssmDouble* North = NULL; + IssmDouble* East = NULL; + int* indices = NULL; + int gsize; + /*Serialize vectors from previous iteration:*/ RSLg=pRSLg->ToMPISerial(); - /*Go through elements, and add contribution from each element to the deflection vector wg:*/ - ns = elements->Size(); + /*Initialize temporary vector that will be used to sum eustatic components on all local elements, prior + * to assembly:*/ + gsize = this->nodes->NumberOfDofs(GsetEnum); + Up=xNewZeroInit(gsize); + North=xNewZeroInit(gsize); + East=xNewZeroInit(gsize); + indices=xNew(gsize); for (int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); eartharea_cpu += element->GetAreaSpherical(); @@ -4897,33 +4906,24 @@ ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - /*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 sea level rise core: */ - for(int i=0;i(elements->GetObjectByOffset(i)); - element->SealevelriseGeodetic(pUp,pNorth,pEast,RSLg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz); - } - if(i%loop==0){ - pUp->Assemble(); - if (horiz){ - pNorth->Assemble(); - pEast->Assemble(); - } - } - } - - /*One last time: */ + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->SealevelriseGeodetic(Up,North,East,RSLg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz); + } + + pUp->SetValues(gsize,indices,Up,ADD_VAL); pUp->Assemble(); - if(horiz){ + if (horiz){ + pNorth->SetValues(gsize,indices,North,ADD_VAL); pNorth->Assemble(); + pEast->SetValues(gsize,indices,East,ADD_VAL); pEast->Assemble(); } - if(VerboseConvergence())_printf0_("\n"); /*Free ressources:*/ + xDelete(Up); + xDelete(North); + xDelete(East); + xDelete(indices); xDelete(RSLg); }