Changeset 26227


Ignore:
Timestamp:
04/30/21 12:33:54 (4 years ago)
Author:
Eric.Larour
Message:

CHG: lumping all loads into a new class GrdLoads.

Location:
issm/trunk-jpl/src/c
Files:
2 added
9 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r26222 r26227  
    559559        ./cores/sealevelchange_core.cpp \
    560560        ./analyses/SealevelchangeAnalysis.cpp\
     561        ./classes/GrdLoads.cpp\
    561562        ./classes/SealevelGeometry.cpp
    562563

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

    r26222 r26227  
    3131class DatasetInput;
    3232class IoModel;
    33 class SealevelMasks;
    3433class SealevelGeometry;
    3534class Gauss;
     35class GrdLoads;
    3636class ElementVector;
    3737template <class doublematrix> class Matrix;
    … …  
    391391
    392392                virtual void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom)=0;
    393                 virtual void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom)=0;
    394                 virtual void       SealevelchangeShift(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementloads, IssmDouble offset, SealevelGeometry* slgeom)=0;
     393                virtual void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0;
     394                virtual void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom)=0;
    395395                virtual void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
    396396                virtual void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
    397397                virtual void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae)=0;
    398                 virtual void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementiceloads, Vector<IssmDouble>* subelementhydroloads, Vector<IssmDouble>* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom)=0;
    399                 virtual void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
    400                 virtual void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
    401                 virtual void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom)=0;
     398                virtual void       SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom)=0;
     399                virtual void       SealevelchangeConvolution(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
     400                virtual void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
     401                virtual void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0;
    402402                #endif
    403403

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

    r26222 r26227  
    224224                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    225225                void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
    226                 void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");};
    227                 void       SealevelchangeShift(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
     226                void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
     227                void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
    228228                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    229229                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    230230                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};
    231                 void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementiceloads, Vector<IssmDouble>* subelementhydroloads, Vector<IssmDouble>* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");};
    232                 void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    233                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    234                 void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");};
     231                void       SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");};
     232                void       SealevelchangeConvolution(GrdLoads* loads,  Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     233                void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     234                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
    235235                #endif
    236236

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

    r26222 r26227  
    179179                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    180180                void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
    181                 void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");};
    182                 void       SealevelchangeShift(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
     181                void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
     182                void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
    183183                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    184184                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    185185                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};
    186                 void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementiceloads, Vector<IssmDouble>* subelementhydroloads, Vector<IssmDouble>* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");};
    187                 void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    188                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    189                 void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");};
     186                void       SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");};
     187                void       SealevelchangeConvolution(GrdLoads* loads,  Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     188                void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     189                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
    190190#endif
    191191

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

    r26222 r26227  
    186186                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    187187                void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
    188                 void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");};
    189                 void       SealevelchangeShift(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
     188                void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
     189                void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
    190190                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    191191                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    192192                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};
    193                 void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementiceloads, Vector<IssmDouble>* subelementhydroloads, Vector<IssmDouble>* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");};
    194                 void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    195                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    196                 void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");};
     193                void       SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");};
     194                void       SealevelchangeConvolution(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     195                void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     196                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
    197197#endif
    198198

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

    r26223 r26227  
    68066806}
    68076807/*}}}*/
    6808 void       Tria::SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementiceloads, Vector<IssmDouble>* subelementhydroloads, Vector<IssmDouble>* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/
     6808void       Tria::SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/
    68096809
    68106810        /*Inputs:*/
    … …  
    68636863        if(slgeom->issubelement[SLGEOM_ICE][this->lid]){
    68646864                int intj=slgeom->subelementmapping[SLGEOM_ICE][this->lid];
    6865                 subelementiceloads->SetValue(intj,Iavg,INS_VAL);
     6865                loads->vsubloads[SLGEOM_ICE]->SetValue(intj,Iavg,INS_VAL);
    68666866                Iavg=0; //avoid double counting centroid loads and subelement loads!
    68676867        }
    68686868        if(slgeom->issubelement[SLGEOM_WATER][this->lid]){
    68696869                int intj=slgeom->subelementmapping[SLGEOM_WATER][this->lid];
    6870                 subelementhydroloads->SetValue(intj,Wavg,INS_VAL);
     6870                loads->vsubloads[SLGEOM_WATER]->SetValue(intj,Wavg,INS_VAL);
    68716871                Wavg=0;
    68726872        }
    68736873        if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    68746874                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    6875                 subelementbploads->SetValue(intj,BPavg,INS_VAL);
     6875                loads->vsubloads[SLGEOM_OCEAN]->SetValue(intj,BPavg,INS_VAL);
    68766876                BPavg=0;
    68776877        }
    68786878        /*Plug remaining values into centroi load vector:*/
    6879         loads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL);
     6879        loads->vloads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL);
    68806880
    68816881        /*Keep track of barystatic contributions:*/
    … …  
    68846884}
    68856885/*}}}*/
    6886 void       Tria::SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/
     6886void       Tria::SealevelchangeConvolution(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/
    68876887
    68886888        /*sal green function:*/
    68896889        IssmDouble* G=NULL;
    6890         IssmDouble* GsubelIce=NULL;
    6891         IssmDouble* GsubelHydro=NULL;
    6892         IssmDouble* GsubelOcean=NULL;
     6890        IssmDouble* Gsub[SLGEOM_NUMLOADS];
    68936891        IssmDouble SealevelGRD[NUMVERTICES]={0};
    68946892        IssmDouble oceanaverage=0;
    … …  
    69196917        if(sal){
    69206918                this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
    6921                 this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&GsubelIce,&size);
    6922                 this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&GsubelHydro,&size);
    6923                 this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&GsubelOcean,&size);
    6924                
    6925                 if(allsealevelloads){
    6926                         for(int i=0;i<NUMVERTICES;i++) {
    6927                                 for (int e=0;e<nel;e++){
    6928                                         SealevelGRD[i]+=G[i*nel+e]*(allsealevelloads[e]+allloads[e]);
    6929                                 }
    6930                                 nbar=slgeom->nbar[SLGEOM_ICE];
    6931                                 for (int e=0;e<nbar;e++){
    6932                                         SealevelGRD[i]+=GsubelIce[i*nbar+e]*(allsubelementiceloads[e]);
    6933                                 }
    6934                                 nbar=slgeom->nbar[SLGEOM_WATER];
    6935                                 for (int e=0;e<nbar;e++){
    6936                                         SealevelGRD[i]+=GsubelHydro[i*nbar+e]*(allsubelementhydroloads[e]);
    6937                                 }
    6938                                 nbar=slgeom->nbar[SLGEOM_OCEAN];
    6939                                 for (int e=0;e<nbar;e++){
    6940                                         SealevelGRD[i]+=GsubelOcean[i*nbar+e]*(allsubelementbploads[e]+allsubelementsealevelloads[e]);
    6941                                 }
    6942                         }
    6943                 }
    6944                 else{  //this is the initial convolution where only loads are provided
    6945                         for(int i=0;i<NUMVERTICES;i++) {
    6946                                 for (int e=0;e<nel;e++){
    6947                                         SealevelGRD[i]+=G[i*nel+e]*(allloads[e]);
    6948                                 }
    6949                                 nbar=slgeom->nbar[SLGEOM_ICE];
    6950                                 for (int e=0;e<nbar;e++){
    6951                                         SealevelGRD[i]+=GsubelIce[i*nbar+e]*(allsubelementiceloads[e]);
    6952                                 }
    6953                                 nbar=slgeom->nbar[SLGEOM_WATER];
    6954                                 for (int e=0;e<nbar;e++){
    6955                                         SealevelGRD[i]+=GsubelHydro[i*nbar+e]*(allsubelementhydroloads[e]);
    6956                                 }
    6957                                 nbar=slgeom->nbar[SLGEOM_OCEAN];
    6958                                 for (int e=0;e<nbar;e++){
    6959                                         SealevelGRD[i]+=GsubelOcean[i*nbar+e]*(allsubelementbploads[e]);
    6960                                 }
    6961                         }
    6962                 }
     6919                this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&Gsub[SLGEOM_ICE],&size);
     6920                this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size);
     6921                this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
     6922
     6923                this->SealevelchangeGxL(&SealevelGRD[0],G, Gsub, loads, slgeom, nel);
    69636924        }
    69646925
    … …  
    69756936        if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    69766937                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    6977                 subelementsealevelloads->SetValue(intj,oceanaverage*rho_water,INS_VAL);
    6978         }
    6979         else sealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
    6980 
    6981         if(!allsealevelloads){
     6938                loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water,INS_VAL);
     6939        }
     6940        else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
     6941
     6942        if(!loads->sealevelloads){
    69826943                oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
    69836944                if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    … …  
    69906951        return;
    69916952} /*}}}*/
    6992 void       Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/
     6953void       Tria::SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/
    69936954
    69946955        IssmDouble Sealevel[3]={0,0,0};
    … …  
    70036964        IssmDouble* GE=NULL;
    70046965        IssmDouble* GN=NULL;
    7005         IssmDouble* GsubelIce=NULL;
    7006         IssmDouble* GsubelHydro=NULL;
    7007         IssmDouble* GsubelOcean=NULL;
    7008         IssmDouble* GUsubelIce=NULL;
    7009         IssmDouble* GUsubelHydro=NULL;
    7010         IssmDouble* GUsubelOcean=NULL;
    7011         IssmDouble* GNsubelIce=NULL;
    7012         IssmDouble* GNsubelHydro=NULL;
    7013         IssmDouble* GNsubelOcean=NULL;
    7014         IssmDouble* GEsubelIce=NULL;
    7015         IssmDouble* GEsubelHydro=NULL;
    7016         IssmDouble* GEsubelOcean=NULL;
     6966        IssmDouble* Gsub[SLGEOM_NUMLOADS];
     6967        IssmDouble* GUsub[SLGEOM_NUMLOADS];
     6968        IssmDouble* GNsub[SLGEOM_NUMLOADS];
     6969        IssmDouble* GEsub[SLGEOM_NUMLOADS];
    70176970
    70186971        int horiz;
    … …  
    70396992               
    70406993        if(sal){
     6994
    70416995                this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
    7042                 this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&GsubelIce,&size);
    7043                 this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&GsubelHydro,&size);
    7044                 this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&GsubelOcean,&size);
     6996                this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&Gsub[SLGEOM_ICE],&size);
     6997                this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size);
     6998                this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
     6999
    70457000                if(elastic){
    70467001                        this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&size);
    7047                         this->inputs->GetArrayPtr(SealevelchangeGUsubelIceEnum,this->lid,&GUsubelIce,&size);
    7048                         this->inputs->GetArrayPtr(SealevelchangeGUsubelHydroEnum,this->lid,&GUsubelHydro,&size);
    7049                         this->inputs->GetArrayPtr(SealevelchangeGUsubelOceanEnum,this->lid,&GUsubelOcean,&size);
     7002                        this->inputs->GetArrayPtr(SealevelchangeGUsubelIceEnum,this->lid,&GUsub[SLGEOM_ICE],&size);
     7003                        this->inputs->GetArrayPtr(SealevelchangeGUsubelHydroEnum,this->lid,&GUsub[SLGEOM_WATER],&size);
     7004                        this->inputs->GetArrayPtr(SealevelchangeGUsubelOceanEnum,this->lid,&GUsub[SLGEOM_OCEAN],&size);
    70507005                        if(horiz){
     7006                                this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size);
     7007                                this->inputs->GetArrayPtr(SealevelchangeGNsubelIceEnum,this->lid,&GNsub[SLGEOM_ICE],&size);
     7008                                this->inputs->GetArrayPtr(SealevelchangeGNsubelHydroEnum,this->lid,&GNsub[SLGEOM_WATER],&size);
     7009                                this->inputs->GetArrayPtr(SealevelchangeGNsubelOceanEnum,this->lid,&GNsub[SLGEOM_OCEAN],&size);
     7010
    70517011                                this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&size);
    7052                                 this->inputs->GetArrayPtr(SealevelchangeGEsubelIceEnum,this->lid,&GEsubelIce,&size);
    7053                                 this->inputs->GetArrayPtr(SealevelchangeGEsubelOceanEnum,this->lid,&GEsubelOcean,&size);
    7054                                 this->inputs->GetArrayPtr(SealevelchangeGEsubelHydroEnum,this->lid,&GEsubelHydro,&size);
    7055                                 this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size);
    7056                                 this->inputs->GetArrayPtr(SealevelchangeGNsubelIceEnum,this->lid,&GNsubelIce,&size);
    7057                                 this->inputs->GetArrayPtr(SealevelchangeGNsubelHydroEnum,this->lid,&GNsubelHydro,&size);
    7058                                 this->inputs->GetArrayPtr(SealevelchangeGNsubelOceanEnum,this->lid,&GNsubelOcean,&size);
    7059                         }
    7060                 }
    7061 
    7062                 for(int i=0;i<NUMVERTICES;i++) {
    7063                         for (int e=0;e<nel;e++){
    7064                                 SealevelRSL[i]+=G[i*nel+e]*(sealevelloads[e]+loads[e]);
    7065                         }
    7066                         nbar=slgeom->nbar[SLGEOM_ICE];
    7067                         for (int e=0;e<nbar;e++){
    7068                                 SealevelRSL[i]+=GsubelIce[i*nbar+e]*(subelementiceloads[e]);
    7069                         }
    7070                         nbar=slgeom->nbar[SLGEOM_WATER];
    7071                         for (int e=0;e<nbar;e++){
    7072                                 SealevelRSL[i]+=GsubelHydro[i*nbar+e]*(subelementhydroloads[e]);
    7073                         }
    7074                         nbar=slgeom->nbar[SLGEOM_OCEAN];
    7075                         for (int e=0;e<nbar;e++){
    7076                                 SealevelRSL[i]+=GsubelOcean[i*nbar+e]*(subelementbploads[e]+subelementsealevelloads[e]);
    7077                         }
    7078 
    7079                         if(elastic){
    7080                                 for (int e=0;e<nel;e++){
    7081                                         SealevelU[i]+=GU[i*nel+e]*(sealevelloads[e]+loads[e]);
    7082                                 }
    7083                                 nbar=slgeom->nbar[SLGEOM_ICE];
    7084                                 for (int e=0;e<nbar;e++){
    7085                                         SealevelU[i]+=GUsubelIce[i*nbar+e]*(subelementiceloads[e]);
    7086                                 }
    7087                                 nbar=slgeom->nbar[SLGEOM_WATER];
    7088                                 for (int e=0;e<nbar;e++){
    7089                                         SealevelU[i]+=GUsubelHydro[i*nbar+e]*(subelementhydroloads[e]);
    7090                                 }
    7091                                 nbar=slgeom->nbar[SLGEOM_OCEAN];
    7092                                 for (int e=0;e<nbar;e++){
    7093                                         SealevelU[i]+=GUsubelOcean[i*nbar+e]*(subelementbploads[e]+subelementsealevelloads[e]);
    7094                                 }
    7095                         }
    7096                         if(horiz && elastic){
    7097                                 for (int e=0;e<nel;e++){
    7098                                         SealevelN[i]+=GN[i*nel+e]*(sealevelloads[e]+loads[e]);
    7099                                         SealevelE[i]+=GE[i*nel+e]*(sealevelloads[e]+loads[e]);
    7100                                 }
    7101                                 nbar=slgeom->nbar[SLGEOM_ICE];
    7102                                 for (int e=0;e<nbar;e++){
    7103                                         SealevelN[i]+=GNsubelIce[i*nbar+e]*(subelementiceloads[e]);
    7104                                         SealevelE[i]+=GEsubelIce[i*nbar+e]*(subelementiceloads[e]);
    7105                                 }
    7106                                 nbar=slgeom->nbar[SLGEOM_WATER];
    7107                                 for (int e=0;e<nbar;e++){
    7108                                         SealevelN[i]+=GNsubelHydro[i*nbar+e]*(subelementhydroloads[e]);
    7109                                         SealevelE[i]+=GEsubelHydro[i*nbar+e]*(subelementhydroloads[e]);
    7110                                 }
    7111                                 nbar=slgeom->nbar[SLGEOM_OCEAN];
    7112                                 for (int e=0;e<nbar;e++){
    7113                                         SealevelN[i]+=GNsubelOcean[i*nbar+e]*(subelementbploads[e]+subelementsealevelloads[e]);
    7114                                         SealevelE[i]+=GEsubelOcean[i*nbar+e]*(subelementbploads[e]+subelementsealevelloads[e]);
    7115                                 }
     7012                                this->inputs->GetArrayPtr(SealevelchangeGEsubelIceEnum,this->lid,&GEsub[SLGEOM_ICE],&size);
     7013                                this->inputs->GetArrayPtr(SealevelchangeGEsubelHydroEnum,this->lid,&GEsub[SLGEOM_WATER],&size);
     7014                                this->inputs->GetArrayPtr(SealevelchangeGEsubelOceanEnum,this->lid,&GEsub[SLGEOM_OCEAN],&size);
     7015                        }
     7016                }
     7017
     7018                this->SealevelchangeGxL(&SealevelRSL[0],G, Gsub, loads, slgeom, nel);
     7019
     7020                if(elastic){
     7021                        this->SealevelchangeGxL(&SealevelU[0],GU, GUsub, loads, slgeom, nel);
     7022                        if(horiz ){
     7023                                this->SealevelchangeGxL(&SealevelN[0],GN, GNsub, loads, slgeom, nel);
     7024                                this->SealevelchangeGxL(&SealevelE[0],GE, GEsub, loads, slgeom, nel);
    71167025                        }
    71177026                }
    … …  
    71317040
    71327041} /*}}}*/
    7133 void       Tria::SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){/*{{{*/
     7042void       Tria::SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel) { /*{{{*/
     7043
     7044        if(loads->sealevelloads){
     7045                for(int i=0;i<NUMVERTICES;i++) {
     7046                        for (int e=0;e<nel;e++){
     7047                                sealevel[i]+=G[i*nel+e]*(loads->sealevelloads[e]+loads->loads[e]);
     7048                        }
     7049                        for(int l=0;l<SLGEOM_NUMLOADS;l++){
     7050                                int nbar=slgeom->nbar[l];
     7051                                for (int e=0;e<nbar;e++){
     7052                                        sealevel[i]+=Gsub[l][i*nbar+e]*(loads->subloads[l][e]);
     7053                                }
     7054                                if(l==SLGEOM_OCEAN){
     7055                                        for (int e=0;e<nbar;e++){
     7056                                                sealevel[i]+=Gsub[l][i*nbar+e]*(loads->subsealevelloads[e]);
     7057                                        }
     7058                                }
     7059                        }
     7060                }
     7061        }
     7062        else{  //this is the initial convolution where only loads are provided
     7063                for(int i=0;i<NUMVERTICES;i++) {
     7064                        for (int e=0;e<nel;e++){
     7065                                sealevel[i]+=G[i*nel+e]*(loads->loads[e]);
     7066                        }
     7067                        for(int l=0;l<SLGEOM_NUMLOADS;l++){
     7068                                int nbar=slgeom->nbar[l];
     7069                                for (int e=0;e<nbar;e++){
     7070                                        sealevel[i]+=Gsub[l][i*nbar+e]*(loads->subloads[l][e]);
     7071                                }
     7072                        }
     7073                }
     7074        }
     7075
     7076} /*}}}*/
     7077void       Tria::SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){/*{{{*/
    71347078               
    71357079        IssmDouble S=0;
    … …  
    71487092
    71497093        /*recover total load: */
    7150         if(loads) S+=loads[this->Sid()];
    7151         if(sealevelloads) S+=sealevelloads[this->Sid()];
     7094        if(loads->loads) S+=loads->loads[this->Sid()];
     7095        if(loads->sealevelloads) S+=loads->sealevelloads[this->Sid()];
    71527096       
    71537097        /* Perturbation terms for moment of inertia (moi_list):
    … …  
    71617105        return;
    71627106}/*}}}*/
    7163 void       Tria::SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){/*{{{*/
     7107void       Tria::SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* grdloads, SealevelGeometry* slgeom){/*{{{*/
    71647108               
    71657109        IssmDouble  SA=0;
    … …  
    71817125        for(int i=0;i<SLGEOM_NUMLOADS;i++){
    71827126                if(slgeom->issubelement[i][this->lid]){
    7183                         switch(i){
    7184                                 case SLGEOM_ICE:
    7185                                         loads=subelementiceloads;
    7186                                         break;
    7187                                 case SLGEOM_WATER:
    7188                                         loads=subelementhydroloads;
    7189                                         break;
    7190                                 case SLGEOM_OCEAN:
    7191                                         loads=subelementbploads;
    7192                                         sealevelloads=subelementsealevelloads;
    7193                                         break;
    7194                         }
     7127                        loads=grdloads->subloads[i];
     7128                        if(i==SLGEOM_OCEAN) sealevelloads=grdloads->subsealevelloads;
    71957129                        intj=slgeom->subelementmapping[i][this->lid];
    71967130                        late=slgeom->latbarycentre[i][intj]/180*M_PI;
    … …  
    72157149        return;
    72167150}/*}}}*/
    7217 void       Tria::SealevelchangeShift(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
     7151void       Tria::SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
    72187152       
    72197153        if(slgeom->isoceanin[this->lid]){
    72207154                if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    72217155                        int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    7222                         subelementloads->SetValue(intj,offset,ADD_VAL);
    7223                 }
    7224                 else loads->SetValue(this->sid,offset,ADD_VAL);
     7156                        loads->vsubsealevelloads->SetValue(intj,offset,ADD_VAL);
     7157                }
     7158                else loads->vsealevelloads->SetValue(this->sid,offset,ADD_VAL);
    72257159        }
    72267160

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

    r26222 r26227  
    173173                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae);
    174174                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae);
    175                 void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementiceloads, Vector<IssmDouble>* subelementhydroloads, Vector<IssmDouble>* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom);
    176                 void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom);
    177                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom);
    178                 void       SealevelchangeShift(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementloads, IssmDouble offset, SealevelGeometry* slgeom);
    179                 void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom);
    180                 void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom);
     175                void       SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom);
     176                void       SealevelchangeConvolution(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom);
     177                void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom);
     178                void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom);
     179                void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom);
     180                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom);
    181181                #endif
    182182                /*}}}*/
    … …  
    214214                Gauss*         NewGauss(void);
    215215                Gauss*         NewGauss(int order);
    216       Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order);
    217       Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order);
    218       Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,int order);
    219       Gauss*         NewGauss(IssmDouble fraction1,IssmDouble fraction2,int order);
    220       Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
     216        Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order);
     217        Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order);
     218        Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,int order);
     219        Gauss*         NewGauss(IssmDouble fraction1,IssmDouble fraction2,int order);
     220        Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
    221221                Gauss*         NewGaussBase(int order);
    222222                Gauss*         NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
    … …  
    237237                void           SetTemporaryElementType(int element_type_in){_error_("not implemented yet");};
    238238                void           InputServe(Input* input_in);
    239                 Seg*             SpawnSeg(int index1,int index2);
     239                Seg*           SpawnSeg(int index1,int index2);
    240240                IssmDouble     StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
    241241                void           StabilizationParameterAnisotropic(IssmDouble* tau_parameter_ansiotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa){_error_("not implemented yet");};
    242242                void           UpdateConstraintsExtrudeFromBase(void);
    243243                void           UpdateConstraintsExtrudeFromTop(void);
     244                void           SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel);
    244245                /*}}}*/
    245246

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

    r26222 r26227  
    1919#include "./Misfit.h"
    2020#include "./SealevelGeometry.h"
     21#include "./GrdLoads.h"
    2122#include "./BarystaticContributions.h"
    2223#include "./Nodalvalue.h"

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

    r26222 r26227  
    2323void TransferForcing(FemModel* femmodel,int forcingenum);
    2424void TransferSealevel(FemModel* femmodel,int forcingenum);
    25 void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
    26 IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble oceanarea);
    27 void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom);
    28 void ConserveOceanMass(FemModel* femmodel,Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, IssmDouble offset, SealevelGeometry* slgeom);
     25bool slcconvergence(Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
     26IssmDouble  SealevelloadsOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble oceanarea);
     27void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom);
     28void ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom);
    2929void ivins_deformation_core(FemModel* femmodel);
    3030IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel);
    … …  
    3232
    3333/*main cores:*/
    34 void sealevelchange_core(FemModel* femmodel){ /*{{{*/
     34void              sealevelchange_core(FemModel* femmodel){ /*{{{*/
    3535
    3636        SealevelGeometry* slgeom=NULL;
    … …  
    8383}
    8484/*}}}*/
    85 void solidearthexternal_core(FemModel* femmodel){ /*{{{*/
     85void              solidearthexternal_core(FemModel* femmodel){ /*{{{*/
    8686
    8787        /*variables:*/
    … …  
    152152}
    153153/*}}}*/
    154 void couplerinput_core(FemModel* femmodel){  /*{{{*/
     154void              couplerinput_core(FemModel* femmodel){  /*{{{*/
    155155
    156156        /*Be very careful here, everything is well thought through, do not remove
    … …  
    215215       
    216216}; /*}}}*/
    217 void grd_core(FemModel* femmodel, SealevelGeometry* slgeom) { /*{{{*/
     217void              grd_core(FemModel* femmodel, SealevelGeometry* slgeom) { /*{{{*/
    218218
    219219        /*variables:{{{*/
    … …  
    223223        IssmDouble rotationaxismotionvector[3]={0};
    224224       
    225         Vector<IssmDouble>*    loads=NULL;
    226         IssmDouble*            allloads=NULL;
    227         Vector<IssmDouble>*    subelementiceloads=NULL;
    228         IssmDouble*            allsubelementiceloads=NULL;
    229         Vector<IssmDouble>*    subelementhydroloads=NULL;
    230         IssmDouble*            allsubelementhydroloads=NULL;
    231         Vector<IssmDouble>*    subelementbploads=NULL;
    232         IssmDouble*            allsubelementbploads=NULL;
    233         Vector<IssmDouble>*    sealevelloads=NULL;
    234         IssmDouble*            allsealevelloads=NULL;
    235         Vector<IssmDouble>*    subelementsealevelloads=NULL;
    236         IssmDouble*            allsubelementsealevelloads=NULL;
     225        GrdLoads*              loads=NULL;
    237226        Vector<IssmDouble>*    oldsealevelloads=NULL;
    238227        Vector<IssmDouble>*    oceanareas=NULL;
    … …  
    249238        int                  step;
    250239        IssmDouble           time;
    251         bool converged=false;
    252240
    253241        int  modelid,earthid;
    … …  
    301289        femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
    302290
    303         loads=new Vector<IssmDouble>(nel);
    304         subelementiceloads=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_ICE]);
    305         subelementhydroloads=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_WATER]);
    306         subelementbploads=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
    307 
    308         sealevelloads=new Vector<IssmDouble>(nel);
    309         sealevelloads->Set(0);sealevelloads->Assemble();
    310         subelementsealevelloads=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
     291        loads=new GrdLoads(nel,slgeom);
    311292        subelementoceanareas=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
    312293        oceanareas=new Vector<IssmDouble>(nel);
    … …  
    317298        for(Object* & object : femmodel->elements->objects){
    318299                Element* element = xDynamicCast<Element*>(object);
    319                 element->SealevelchangeBarystaticLoads(loads, subelementiceloads, subelementhydroloads, subelementbploads, barycontrib,slgeom);
    320         }
    321         loads->Assemble();
    322         subelementiceloads->Assemble();
    323         subelementhydroloads->Assemble();
    324         subelementbploads->Assemble();
    325        
     300                element->SealevelchangeBarystaticLoads(loads, barycontrib,slgeom);
     301        }
    326302
    327303        //broadcast loads
    328         allloads=loads->ToMPISerial();
    329         allsubelementiceloads=subelementiceloads->ToMPISerial();
    330         allsubelementhydroloads=subelementhydroloads->ToMPISerial();
    331         allsubelementbploads=subelementbploads->ToMPISerial();
     304        loads->BroadcastLoads();
    332305
    333306        //compute rotation axis motion:
    334         RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsubelementiceloads,allsubelementhydroloads,allsubelementbploads, NULL, NULL,slgeom);
     307        RotationAxisMotion(&rotationaxismotionvector[0],femmodel,loads,slgeom);
    335308
    336309        /*skip computation of sea level if requested, which means sea level loads should be zeroed */
    337310        if(!computesealevel){
    338                 allsealevelloads=xNewZeroInit<IssmDouble>(nel);
    339                 allsubelementsealevelloads=xNewZeroInit<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
     311                loads->sealevelloads=xNewZeroInit<IssmDouble>(nel);
     312                loads->subsealevelloads=xNewZeroInit<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
    340313                goto deformation;
    341314        }
    … …  
    344317        for(;;){
    345318
    346                 oldsealevelloads=sealevelloads->Duplicate(); sealevelloads->Copy(oldsealevelloads);
     319                oldsealevelloads=loads->vsealevelloads->Duplicate(); loads->vsealevelloads->Copy(oldsealevelloads);
    347320
    348321                /*convolve load and sealevel loads on oceans:*/
    349322                for(Object* & object : femmodel->elements->objects){
    350323                        Element* element = xDynamicCast<Element*>(object);
    351                         element->SealevelchangeConvolution(sealevelloads, subelementsealevelloads, oceanareas, subelementoceanareas, allsealevelloads, allloads,allsubelementiceloads, allsubelementhydroloads, allsubelementbploads, allsubelementsealevelloads, rotationaxismotionvector,slgeom);
    352                 }
    353 
    354                 sealevelloads->Assemble();
    355                 subelementsealevelloads->Assemble();
    356 
     324                        element->SealevelchangeConvolution(loads, oceanareas, subelementoceanareas, rotationaxismotionvector,slgeom);
     325                }
     326
     327                loads->AssembleSealevelLoads();
     328       
    357329                /*compute ocean areas:*/
    358                 if(!allsealevelloads){ //first time in the loop
     330                if(!loads->sealevelloads){ //first time in the loop
    359331                        oceanareas->Assemble();
    360332                        subelementoceanareas->Assemble();
    … …  
    364336       
    365337                //Conserve ocean mass:
    366                 oceanaverage=SealevelloadsOceanAverage(sealevelloads,subelementsealevelloads, oceanareas,subelementoceanareas, oceanarea);
    367                
    368                 ConserveOceanMass(femmodel,sealevelloads,subelementsealevelloads,barycontrib->Total()/oceanarea - oceanaverage,slgeom);
     338                oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, oceanarea);
     339                ConserveOceanMass(femmodel,loads,barycontrib->Total()/oceanarea - oceanaverage,slgeom);
    369340
    370341                //broadcast sea level loads
    371                 allsealevelloads=sealevelloads->ToMPISerial();
    372                 allsubelementsealevelloads=subelementsealevelloads->ToMPISerial();
     342                loads->BroadcastSealevelLoads();
    373343
    374344                //compute rotation axis motion:
    375                 RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsubelementiceloads,allsubelementhydroloads,allsubelementbploads, allsealevelloads,allsubelementsealevelloads,slgeom);
     345                RotationAxisMotion(&rotationaxismotionvector[0],femmodel,loads, slgeom);
    376346
    377347                //convergence?
    378                 slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs);
    379                 if (converged)break;
     348                if(slcconvergence(loads->vsealevelloads,oldsealevelloads,eps_rel,eps_abs))break;
    380349
    381350                //early return?
    382351                if(iterations>=max_nonlinear_iterations)break;
    383                 else iterations++;
     352                iterations++;
    384353        }
    385354
    … …  
    391360        for(Object* & object : femmodel->elements->objects){
    392361                Element* element = xDynamicCast<Element*>(object);
    393                 element->SealevelchangeDeformationConvolution(allsealevelloads, allsubelementsealevelloads,allloads, allsubelementiceloads, allsubelementhydroloads, allsubelementbploads,rotationaxismotionvector,slgeom);
     362                element->SealevelchangeDeformationConvolution(loads, rotationaxismotionvector,slgeom);
    394363        }
    395364
    … …  
    414383}
    415384/*}}}*/
    416 void dynstr_core(FemModel* femmodel){ /*{{{*/
     385void              dynstr_core(FemModel* femmodel){ /*{{{*/
    417386
    418387        /*variables:*/
    … …  
    480449}
    481450/*}}}*/
    482 void coupleroutput_core(FemModel* femmodel){  /*{{{*/
     451void              coupleroutput_core(FemModel* femmodel){  /*{{{*/
    483452       
    484453        /*parameters:*/
    … …  
    500469        }
    501470}; /*}}}*/
    502 void ivins_deformation_core(FemModel* femmodel){ /*{{{*/
     471void              ivins_deformation_core(FemModel* femmodel){ /*{{{*/
    503472
    504473        int  gsize;
    … …  
    551520}
    552521/*}}}*/
    553 
    554 //Geometry:
    555 void sealevelchange_initialgeometry(FemModel* femmodel) {  /*{{{*/
     522void              sealevelchange_initialgeometry(FemModel* femmodel) {  /*{{{*/
    556523
    557524        /*Geometry core where we compute geometrical kernels and weights:*/
    … …  
    668635}/*}}}*/
    669636
    670 /*Support routines:*/
     637/*subroutines:*/
     638bool slcconvergence(Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
     639
     640        bool converged=true;
     641        IssmDouble ndS,nS;
     642        Vector<IssmDouble> *dRSLg    = NULL;
     643
     644        //compute norm(du) and norm(u) if requested
     645        dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0);
     646        ndS=dRSLg->Norm(NORM_TWO);
     647
     648        if (xIsNan<IssmDouble>(ndS)){
     649                _error_("convergence criterion is NaN (RSL_old=" << RSLg_old->Norm(NORM_TWO) << " RSL=" << RSLg->Norm(NORM_TWO) << ")");
     650        }
     651
     652        if(!xIsNan<IssmDouble>(eps_rel)){
     653                nS=RSLg_old->Norm(NORM_TWO);
     654                if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN! (check the initial RSL)");
     655        }
     656
     657        //clean up
     658        delete dRSLg;
     659
     660        //print
     661        if(!xIsNan<IssmDouble>(eps_rel)){
     662                if((ndS/nS)<eps_rel){
     663                        if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " < " << eps_rel*100 << " %\n");
     664                }
     665                else{
     666                        if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " > " << eps_rel*100 << " %\n");
     667                        converged=false;
     668                }
     669        }
     670        if(!xIsNan<IssmDouble>(eps_abs)){
     671                if(ndS<eps_abs){
     672                        if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)" << ndS << " < " << eps_abs << " \n");
     673                }
     674                else{
     675                        if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)" << ndS << " > " << eps_abs << " \n");
     676                        converged=false;
     677                }
     678        }
     679
     680        /*assign output*/
     681        return converged;
     682
     683} /*}}}*/
     684IssmDouble  SealevelloadsOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* suboceanareas, IssmDouble oceanarea){ /*{{{*/
     685
     686        IssmDouble sealevelloadsaverage;       
     687        IssmDouble subsealevelloadsaverage;     
     688
     689        Vector<IssmDouble>* vsealevelloadsvolume=loads->vsealevelloads->Duplicate();
     690        Vector<IssmDouble>* vsubsealevelloadsvolume=loads->vsubsealevelloads->Duplicate();
     691
     692        vsealevelloadsvolume->PointwiseMult(loads->vsealevelloads,oceanareas);
     693        vsubsealevelloadsvolume->PointwiseMult(loads->vsubsealevelloads,suboceanareas);
     694       
     695        vsealevelloadsvolume->Sum(&sealevelloadsaverage);
     696        vsubsealevelloadsvolume->Sum(&subsealevelloadsaverage);
     697        delete vsealevelloadsvolume;
     698        delete vsubsealevelloadsvolume;
     699       
     700        //return (sealevelloadsaverage+subsealevelloadsaverage)/oceanarea;
     701        return (sealevelloadsaverage)/oceanarea;
     702} /*}}}*/
     703void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom){ /*{{{*/
     704
     705        IssmDouble  moi_list[3]={0,0,0};
     706        IssmDouble  moi_list_sub[3]={0,0,0};
     707        IssmDouble  moi_list_cpu[3]={0,0,0};
     708        IssmDouble*     tide_love_h  = NULL;
     709        IssmDouble*     tide_love_k  = NULL;
     710        IssmDouble*     load_love_k  = NULL;
     711        IssmDouble  tide_love_k2secular;
     712        IssmDouble  moi_e, moi_p;
     713        IssmDouble      m1, m2, m3;
     714        bool rotation=false;
     715       
     716        /*early return?:*/
     717        femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
     718        if(!rotation)return;
     719
     720        /*retrieve parameters: */
     721        femmodel->parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum);
     722        femmodel->parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
     723        femmodel->parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
     724        femmodel->parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum);
     725        femmodel->parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
     726        femmodel->parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
     727
     728        for(Object* & object : femmodel->elements->objects){
     729                Element* element = xDynamicCast<Element*>(object);
     730
     731                element->SealevelchangeMomentOfInertiaCentroid(&moi_list[0],loads,slgeom);
     732
     733                element->SealevelchangeMomentOfInertiaSubElement(&moi_list_sub[0],loads, slgeom);
     734
     735                moi_list_cpu[0] += moi_list[0]+moi_list_sub[0];
     736                moi_list_cpu[1] += moi_list[1]+moi_list_sub[1];
     737                moi_list_cpu[2] += moi_list[2]+moi_list_sub[2];
     738        }
     739
     740
     741        ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     742        ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     743       
     744        ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     745        ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     746       
     747        ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     748        ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     749
     750        /*compute perturbation terms for angular velocity vector: */
     751        m1 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[0];
     752        m2 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[1];
     753        m3 = -(1+load_love_k[2])/moi_p * moi_list[2];   // term associated with fluid number (3-order-of-magnitude smaller) is negelected
     754
     755        /*Assign output pointers:*/
     756        m[0]=m1;
     757        m[1]=m2;
     758        m[2]=m3;
     759} /*}}}*/
     760void        ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
     761
     762        /*Shift sealevel loads by ocean average, only on ocean! :*/
     763        for(Object* & object : femmodel->elements->objects){
     764                Element* element = xDynamicCast<Element*>(object);
     765                element->SealevelchangeShift(loads, offset,slgeom);
     766        }
     767        loads->AssembleSealevelLoads();
     768
     769} /*}}}*/
     770IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel){ /*{{{*/
     771
     772        int* indices=xNew<int>(nel);
     773        for(int i=0;i<nel;i++)indices[i]=i;
     774       
     775        Vector<IssmDouble>* vloadcopy=new Vector<IssmDouble>(nel);
     776        IssmDouble* loadcopy=xNew<IssmDouble>(nel);
     777       
     778        vloadcopy->SetValues(nel,indices,load,INS_VAL);
     779        vloadcopy->Assemble();
     780
     781
     782        if(subload){
     783                for (int i=0;i<femmodel->elements->Size();i++){
     784                        if (slgeom->issubelement[loadtype][i]){
     785                                int se= slgeom->subelementmapping[loadtype][i];
     786                                IssmDouble subloadi=subload[se];
     787                                Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
     788                                vloadcopy->SetValue(element->Sid(),subloadi,ADD_VAL);
     789                        }
     790                }
     791        }
     792        vloadcopy->Assemble();
     793        loadcopy=vloadcopy->ToMPISerial();
     794
     795        return loadcopy;
     796
     797} /*}}}*/
     798
     799/*Coupling routines:*/
    671800void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/
    672801
    … …  
    9141043
    9151044} /*}}}*/
    916 void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
    917 
    918         bool converged=true;
    919         IssmDouble ndS,nS;
    920         Vector<IssmDouble> *dRSLg    = NULL;
    921 
    922         //compute norm(du) and norm(u) if requested
    923         dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0);
    924         ndS=dRSLg->Norm(NORM_TWO);
    925 
    926         if (xIsNan<IssmDouble>(ndS)){
    927                 _error_("convergence criterion is NaN (RSL_old=" << RSLg_old->Norm(NORM_TWO) << " RSL=" << RSLg->Norm(NORM_TWO) << ")");
    928         }
    929 
    930         if(!xIsNan<IssmDouble>(eps_rel)){
    931                 nS=RSLg_old->Norm(NORM_TWO);
    932                 if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN! (check the initial RSL)");
    933         }
    934 
    935         //clean up
    936         delete dRSLg;
    937 
    938         //print
    939         if(!xIsNan<IssmDouble>(eps_rel)){
    940                 if((ndS/nS)<eps_rel){
    941                         if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " < " << eps_rel*100 << " %\n");
    942                 }
    943                 else{
    944                         if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " > " << eps_rel*100 << " %\n");
    945                         converged=false;
    946                 }
    947         }
    948         if(!xIsNan<IssmDouble>(eps_abs)){
    949                 if(ndS<eps_abs){
    950                         if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)" << ndS << " < " << eps_abs << " \n");
    951                 }
    952                 else{
    953                         if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)" << ndS << " > " << eps_abs << " \n");
    954                         converged=false;
    955                 }
    956         }
    957 
    958         /*assign output*/
    959         *pconverged=converged;
    960 
    961 } /*}}}*/
    962 IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble oceanarea){ /*{{{*/
    963 
    964         IssmDouble sealevelloadsaverage;       
    965         IssmDouble subelementsealevelloadsaverage;     
    966 
    967         Vector<IssmDouble>* sealevelloadsvolume=sealevelloads->Duplicate();
    968         Vector<IssmDouble>* subelementsealevelloadsvolume=subelementsealevelloads->Duplicate();
    969 
    970         sealevelloadsvolume->PointwiseMult(sealevelloads,oceanareas);
    971         subelementsealevelloadsvolume->PointwiseMult(subelementsealevelloads,subelementoceanareas);
    972        
    973         sealevelloadsvolume->Sum(&sealevelloadsaverage);
    974         subelementsealevelloadsvolume->Sum(&subelementsealevelloadsaverage);
    975         delete sealevelloadsvolume;
    976         delete subelementsealevelloadsvolume;
    977        
    978         //return (sealevelloadsaverage+subelementsealevelloadsaverage)/oceanarea;
    979         return (sealevelloadsaverage)/oceanarea;
    980 } /*}}}*/
    981 void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){ /*{{{*/
    982 
    983         IssmDouble  moi_list[3]={0,0,0};
    984         IssmDouble  moi_list_sub[3]={0,0,0};
    985         IssmDouble  moi_list_cpu[3]={0,0,0};
    986         IssmDouble*     tide_love_h  = NULL;
    987         IssmDouble*     tide_love_k  = NULL;
    988         IssmDouble*     load_love_k  = NULL;
    989         IssmDouble  tide_love_k2secular;
    990         IssmDouble  moi_e, moi_p;
    991         IssmDouble      m1, m2, m3;
    992         bool rotation=false;
    993        
    994         /*early return?:*/
    995         femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
    996         if(!rotation)return;
    997 
    998         /*retrieve parameters: */
    999         femmodel->parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum);
    1000         femmodel->parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
    1001         femmodel->parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
    1002         femmodel->parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum);
    1003         femmodel->parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
    1004         femmodel->parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
    1005 
    1006         for(Object* & object : femmodel->elements->objects){
    1007                 Element* element = xDynamicCast<Element*>(object);
    1008 
    1009                 element->SealevelchangeMomentOfInertiaCentroid(&moi_list[0],loads,sealevelloads,slgeom);
    1010 
    1011                 element->SealevelchangeMomentOfInertiaSubElement(&moi_list_sub[0],subelementiceloads, subelementhydroloads, subelementbploads, subelementsealevelloads,slgeom);
    1012 
    1013                 moi_list_cpu[0] += moi_list[0]+moi_list_sub[0];
    1014                 moi_list_cpu[1] += moi_list[1]+moi_list_sub[1];
    1015                 moi_list_cpu[2] += moi_list[2]+moi_list_sub[2];
    1016         }
    1017 
    1018 
    1019         ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    1020         ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    1021        
    1022         ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    1023         ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    1024        
    1025         ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    1026         ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    1027 
    1028         /*compute perturbation terms for angular velocity vector: */
    1029         m1 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[0];
    1030         m2 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[1];
    1031         m3 = -(1+load_love_k[2])/moi_p * moi_list[2];   // term associated with fluid number (3-order-of-magnitude smaller) is negelected
    1032 
    1033         /*Assign output pointers:*/
    1034         m[0]=m1;
    1035         m[1]=m2;
    1036         m[2]=m3;
    1037 } /*}}}*/
    1038 void ConserveOceanMass(FemModel* femmodel,Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
    1039 
    1040         /*Shift sealevel loads by ocean average, only on ocean! :*/
    1041         for(Object* & object : femmodel->elements->objects){
    1042                 Element* element = xDynamicCast<Element*>(object);
    1043                 element->SealevelchangeShift(sealevelloads,subelementsealevelloads,offset,slgeom);
    1044         }
    1045         sealevelloads->Assemble();
    1046         subelementsealevelloads->Assemble();
    1047 
    1048 } /*}}}*/
    1049 IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel){ /*{{{*/
    1050 
    1051         int* indices=xNew<int>(nel);
    1052         for(int i=0;i<nel;i++)indices[i]=i;
    1053        
    1054         Vector<IssmDouble>* vloadcopy=new Vector<IssmDouble>(nel);
    1055         IssmDouble* loadcopy=xNew<IssmDouble>(nel);
    1056        
    1057         vloadcopy->SetValues(nel,indices,load,INS_VAL);
    1058         vloadcopy->Assemble();
    1059 
    1060 
    1061         if(subload){
    1062                 for (int i=0;i<femmodel->elements->Size();i++){
    1063                         if (slgeom->issubelement[loadtype][i]){
    1064                                 int se= slgeom->subelementmapping[loadtype][i];
    1065                                 IssmDouble subloadi=subload[se];
    1066                                 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    1067                                 vloadcopy->SetValue(element->Sid(),subloadi,ADD_VAL);
    1068                         }
    1069                 }
    1070         }
    1071         vloadcopy->Assemble();
    1072         loadcopy=vloadcopy->ToMPISerial();
    1073 
    1074         return loadcopy;
    1075 
    1076 } /*}}}*/
    1077 
Note: See TracChangeset for help on using the changeset viewer.