Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26163) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26164) @@ -248,4 +248,6 @@ virtual void GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0; virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0; + virtual void GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl)=0; + virtual IssmDouble GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl)=0; virtual IssmDouble GetIcefrontArea(){_error_("not implemented");}; virtual void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0; @@ -295,4 +297,5 @@ virtual void JacobianDeterminantSurface(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss)=0; virtual void JacobianDeterminantTop(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0; + virtual void LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum)=0; virtual void Marshall(MarshallHandle* marshallhandle)=0; virtual IssmDouble Masscon(IssmDouble* levelset)=0; @@ -394,8 +397,6 @@ virtual void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0; virtual void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0; - virtual void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks)=0; - virtual void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks)=0; - virtual void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks)=0; - virtual void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks)=0; + virtual void SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector)=0; + virtual void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector)=0; virtual void SealevelchangeShift(Vector* loads, IssmDouble offset, SealevelMasks* masks)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26163) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26164) @@ -83,4 +83,6 @@ void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating); IssmDouble GetGroundedPortion(IssmDouble* xyz_list); + void GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){_error_("not implemented yet");}; + IssmDouble GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){_error_("not implemented yet");}; IssmDouble GetIcefrontArea(); void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum); @@ -125,4 +127,5 @@ void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); + void LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum){_error_("not implemented yet");}; IssmDouble Masscon(IssmDouble* levelset){_error_("not implemented yet");}; IssmDouble MassFlux(IssmDouble* segment); @@ -231,8 +234,6 @@ void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector,SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; - void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){_error_("not implemented yet");}; void SealevelchangeShift(Vector* loads, IssmDouble offset, SealevelMasks* masks){_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 26163) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26164) @@ -63,4 +63,6 @@ void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");}; IssmDouble GetGroundedPortion(IssmDouble* xyz_list); + void GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){_error_("not implemented yet");}; + IssmDouble GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){_error_("not implemented yet");}; void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum); Input* GetInput(int enumtype); @@ -97,4 +99,5 @@ void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; + void LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum){_error_("not implemented yet");}; IssmDouble Masscon(IssmDouble* levelset){_error_("not implemented yet");}; IssmDouble MassFlux(IssmDouble* segment){_error_("not implemented yet");}; @@ -186,8 +189,6 @@ void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; - void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){_error_("not implemented yet");}; void SealevelchangeShift(Vector* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26163) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26164) @@ -67,4 +67,6 @@ void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");}; IssmDouble GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");}; + void GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){_error_("not implemented yet");}; + IssmDouble GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){_error_("not implemented yet");}; void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");}; Input* GetInput(int enumtype); @@ -103,4 +105,5 @@ void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); + void LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum){_error_("not implemented yet");}; IssmDouble Masscon(IssmDouble* levelset){_error_("not implemented yet");}; IssmDouble MassFlux(IssmDouble* segment){_error_("not implemented yet");}; @@ -192,8 +195,6 @@ void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; - void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){_error_("not implemented yet");}; void SealevelchangeShift(Vector* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26163) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26164) @@ -1600,5 +1600,5 @@ } /*}}}*/ -void Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){/*{{{*/ +void Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlyfloating){/*{{{*/ /*Computeportion of the element that is grounded*/ @@ -1651,5 +1651,5 @@ *fraction1=f1; *fraction2=f2; - *mainlyfloating=floating; + *pmainlyfloating=floating; } /*}}}*/ @@ -1721,4 +1721,121 @@ _assert_(phi<=1. && phi>=0.); return phi; +} +/*}}}*/ +IssmDouble Tria::GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){/*{{{*/ + + /*Computeportion of the element that is on a levelset: */ + bool mainlyfloating = true; + int domaintype,index1,index2; + const IssmPDouble epsilon = 1.e-15; + IssmDouble phi,s1,s2; + + /*Recover parameters and values*/ + parameters->FindParam(&domaintype,DomainTypeEnum); + + /*Be sure that values are not zero*/ + if(gl[0]==0.) gl[0]=gl[0]+epsilon; + if(gl[1]==0.) gl[1]=gl[1]+epsilon; + if(gl[2]==0.) gl[2]=gl[2]+epsilon; + + if(domaintype==Domain2DverticalEnum){ + this->EdgeOnBaseIndices(&index1,&index2); + if(gl[index1]>0 && gl[index2]>0) phi=1; // All grounded + else if(gl[index1]<0 && gl[index2]<0) phi=0; // All floating + else if(gl[index1]<0 && gl[index2]>0){ //index2 grounded + phi=1./(1.-gl[index1]/gl[index2]); + } + else if(gl[index2]<0 && gl[index1]>0){ //index1 grounded + phi=1./(1.-gl[index2]/gl[index1]); + } + + } + else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){ + /*Check that not all nodes are grounded or floating*/ + if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded + phi=1; + } + else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating + phi=0; + } + else{ + /*Figure out if two nodes are floating or grounded*/ + if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=false; + + if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 + s1=gl[2]/(gl[2]-gl[1]); + s2=gl[2]/(gl[2]-gl[0]); + } + else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 + s1=gl[0]/(gl[0]-gl[1]); + s2=gl[0]/(gl[0]-gl[2]); + } + else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 + s1=gl[1]/(gl[1]-gl[0]); + s2=gl[1]/(gl[1]-gl[2]); + } + else _error_("case not possible"); + if(mainlyfloating){ + phi = (1-s1*s2); + } + else{ + phi = s1*s2; + } + } + } + else _error_("mesh type "<=0.); + return 1-phi; +} +/*}}}*/ +void Tria::GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){/*{{{*/ + + /*Computeportion of the element that is grounded*/ + bool negative=false; + int point; + const IssmPDouble epsilon= 1.e-15; + IssmDouble f1,f2; + + /*Be sure that values are not zero*/ + if(gl[0]==0.) gl[0]=gl[0]+epsilon; + if(gl[1]==0.) gl[1]=gl[1]+epsilon; + if(gl[2]==0.) gl[2]=gl[2]+epsilon; + + /*Check that not all nodes are positive or negative: */ + if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All positive + point=0; + f1=1.; + f2=1.; + } + else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All negative + point=0; + f1=0.; + f2=0.; + } + else{ + if(gl[0]*gl[1]*gl[2]<0) negative=true; + + if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 + point=2; + f1=gl[2]/(gl[2]-gl[0]); + f2=gl[2]/(gl[2]-gl[1]); + } + else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 + point=0; + f1=gl[0]/(gl[0]-gl[1]); + f2=gl[0]/(gl[0]-gl[2]); + } + else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 + point=1; + f1=gl[1]/(gl[1]-gl[2]); + f2=gl[1]/(gl[1]-gl[0]); + } + else _error_("case not possible"); + } + *point1=point; + *fraction1=f1; + *fraction2=f2; + *pmainlynegative=negative; } /*}}}*/ @@ -5545,15 +5662,4 @@ IssmDouble phi=1.0; /*{{{*/ - - - - if (masks->isfullyfloating[this->lid]){ - I=0; - this->AddInput(EffectivePressureEnum,&I,P0Enum); - return; - } - - - if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on. @@ -5605,5 +5711,5 @@ total_weight+=gauss->weight; } - I=I/total_weight; + if(total_weight) I=I/total_weight; delete gauss; } @@ -6551,5 +6657,5 @@ Grotm1[i]= - pre* 0.5*sin(2.*lati)*cos(longi); Grotm2[i]= - pre* 0.5*sin(2.*lati)*sin(longi); - Grotm3[i]= - pre* (1/6.0 - 0.5*cos(2.0*lati)); + Grotm3[i]= - pre* (1.0/6.0 - 0.5*cos(2.0*lati)); } this->AddInput(SealevelGrotm1Enum,&Grotm1[0],P1Enum); @@ -6687,5 +6793,5 @@ /*Deal with ice loads if we are on grounded ice:*/ - if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){ + if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid] && isice){ /*Compute fraction of the element that is grounded: {{{*/ @@ -6734,5 +6840,5 @@ total_weight+=gauss->weight; } - I=I/total_weight; + if(total_weight) I=I/total_weight; delete gauss; } @@ -6748,5 +6854,5 @@ /*Deal with water loads if we are on ground:*/ - if(!masks->isfullyfloating[this->lid]){ + if(!masks->isfullyfloating[this->lid] && ishydro){ constant+=10; //1 for water. @@ -6769,5 +6875,5 @@ /*Deal with ocean bottom pressures:*/ - if(masks->isoceanin[this->lid]){ + if(masks->isoceanin[this->lid] && isocean){ constant+=1; //1 for bottom pressure. @@ -6796,9 +6902,9 @@ } /*}}}*/ -void Tria::SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){ /*{{{*/ +void Tria::SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationvector){ /*{{{*/ /*sal green function:*/ IssmDouble* G=NULL; - IssmDouble SealevelGRD[NUMVERTICES]={0,0,0}; + IssmDouble SealevelGRD[NUMVERTICES]={0}; IssmDouble oceanaverage,oceanarea=0; IssmDouble rho_water; @@ -6828,7 +6934,16 @@ this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size); - for(int i=0;iSetValue(this->sid,oceanaverage*rho_water,INS_VAL); - oceanareas->SetValue(this->sid,oceanarea,INS_VAL); + if(!allsealevelloads)oceanareas->SetValue(this->sid,oceanarea,INS_VAL); return; } /*}}}*/ -void Tria::SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){ /*{{{*/ - - IssmDouble SealevelGRD[3]={0,0,0}; - IssmDouble oceanaverage,oceanarea=0; - IssmDouble rho_water; - int nel; - bool sal = false; - IssmDouble* G=NULL; - int size; - IssmDouble Grotm1[3]; - IssmDouble Grotm2[3]; - IssmDouble Grotm3[3]; - bool rotation= false; - - this->parameters->FindParam(&nel,MeshNumberofelementsEnum); - this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum); - this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); - this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum); - - if(rotation){ - Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum); - Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum); - Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum); - - for(int i=0;iinputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size); - - for(int i=0;iSetValue(this->sid,oceanaverage*rho_water,INS_VAL); - -} /*}}}*/ -void Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){ /*{{{*/ +void Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){ /*{{{*/ IssmDouble Sealevel[3]={0,0,0}; @@ -6958,24 +7033,26 @@ } /*}}}*/ -void Tria::OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/ +void Tria::LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum){ /*{{{*/ IssmDouble phi=1.0; - bool iscoastline=false; IssmDouble area; - IssmDouble Sg_avg=0; //output + IssmDouble average = 0; + IssmDouble total_weight=0; + int point1; + IssmDouble fraction1,fraction2; + IssmDouble levelsetvalues[NUMVERTICES]; + bool fractiongeometryflag=true; - /*Do we have an ocean?:*/ - if(!masks->isoceanin[this->lid]){ - *poceanarea=0; - *poceanaverage=0; - } - - /*Do we have a coastline?:*/ - if(!masks->isfullyfloating[this->lid])iscoastline=true; - - if(iscoastline){ - IssmDouble xyz_list[NUMVERTICES][3]; - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - phi=1.0-this->GetGroundedPortion(&xyz_list[0][0]); + + + /*retrieve value of levelset:*/ + Input *levelset= this->GetInput(levelsetenum); + this->Element::GetInputListOnVertices(&levelsetvalues[0],levelsetenum); + + /*Early return if no vertices are inside the desired area:*/ + if(levelset->GetInputMin()>=0){ + *paverage=0; + *parea=0; + return; } @@ -6983,44 +7060,38 @@ this->Element::GetInputValue(&area,AreaEnum); - /*Average over ocean if there is no coastline, area of the ocean - *is the areaa of the element:*/ - if(!iscoastline){ - - /*Average Sg over vertices:*/ - for(int i=0;iGetInputMax()<=0){ + for(int i=0;iGetFractionArea(&xyz_list[0][0],levelsetvalues); area=phi*area; - - IssmDouble total_weight=0; - bool mainlyfloating = true; - int point1; - IssmDouble fraction1,fraction2; - - /*Recover portion of element that is grounded*/ - this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); - //!mainlyfloating so that the integration is done on the ocean (and not the grounded) part - Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,!mainlyfloating,2); - - /* Start looping on the number of gaussian points and average over these gaussian points: */ + + /*Average over the fraction area only, using the right gaussian points: */ + this->GetFractionGeometry(&point1,&fraction1,&fraction2,&fractiongeometryflag,levelsetvalues); + Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,fractiongeometryflag,2); + total_weight=0; - Sg_avg=0; + average=0; while(gauss->next()){ - IssmDouble Sg_gauss=0; - TriaRef::GetInputValue(&Sg_gauss, Sg, gauss,P1Enum); - Sg_avg+=Sg_gauss*gauss->weight; + IssmDouble field_gauss=0; + TriaRef::GetInputValue(&field_gauss, field_on_localvertices, gauss,P1Enum); + average+=field_gauss*gauss->weight; total_weight+=gauss->weight; } - if(total_weight) Sg_avg=Sg_avg/total_weight; + if(total_weight) average=average/total_weight; + + /*free ressources:*/ delete gauss; - *poceanaverage=Sg_avg; - *poceanarea=area; + *paverage=average; + *parea=area; return; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26163) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26164) @@ -84,4 +84,6 @@ void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating); IssmDouble GetGroundedPortion(IssmDouble* xyz_list); + void GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl); + IssmDouble GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl); IssmDouble GetIcefrontArea(); void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum); @@ -111,4 +113,5 @@ bool IsNodeOnShelfFromFlags(IssmDouble* flags); bool IsZeroLevelset(int levelset_enum); + void LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum); IssmDouble Masscon(IssmDouble* levelset); IssmDouble MassFlux(IssmDouble* segment); @@ -176,8 +179,6 @@ void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae); void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks); - void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks); - void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks); - void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks); - void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks); + void SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector); + void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector); void SealevelchangeShift(Vector* loads, IssmDouble offset, SealevelMasks* masks); #endif Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26163) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26164) @@ -405,5 +405,5 @@ SealevelMasks* masks=NULL; GenericParam* barycontribparam=NULL; - IssmDouble rotationaxismotionvector[3]; + IssmDouble rotationaxismotionvector[3]={0}; Vector* loads=NULL; @@ -468,6 +468,7 @@ femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum); loads=new Vector(nel); + oceanareas=new Vector(nel); sealevelloads=new Vector(nel); - oceanareas=new Vector(nel); + sealevelloads->Set(0);sealevelloads->Assemble(); /*call masks core: */ @@ -494,28 +495,5 @@ } - /*convolve loads:*/ - for(Object* & object : femmodel->elements->objects){ - Element* element = xDynamicCast(object); - element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,rotationaxismotionvector,masks); - } - sealevelloads->Assemble(); - oceanareas->Assemble(); - - //Get ocean area: - oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.); - if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 - - //conserve ocean mass: - oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea); - ConserveOceanMass(femmodel,sealevelloads,barycontrib->Total()/oceanarea - oceanaverage,masks); - - //broadcast sea level loads - allsealevelloads=sealevelloads->ToMPISerial(); - - - //compute rotation axis motion: - RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsealevelloads); - - if(VerboseSolution()) _printf0_(" converging ocean GRD convolutions\n"); + if(VerboseSolution()) _printf0_(" converging GRD convolutions\n"); for(;;){ @@ -525,8 +503,14 @@ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,rotationaxismotionvector,masks); + element->SealevelchangeConvolution(sealevelloads, oceanareas, allsealevelloads, allloads,rotationaxismotionvector); } sealevelloads->Assemble(); - + + /*compute ocean areas:*/ + if(!allsealevelloads){ //first time in the loop + oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.); + if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 + } + //Conserve ocean mass: oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea); @@ -544,5 +528,5 @@ //early return? - if(iterations>max_nonlinear_iterations)break; + if(iterations>=max_nonlinear_iterations)break; else iterations++; } @@ -555,5 +539,5 @@ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeDeformationConvolution(allsealevelloads, allloads, rotationaxismotionvector,masks); + element->SealevelchangeDeformationConvolution(allsealevelloads, allloads, rotationaxismotionvector); } @@ -1227,9 +1211,11 @@ ndS=dRSLg->Norm(NORM_TWO); - if (xIsNan(ndS)) _error_("convergence criterion is NaN!"); + if (xIsNan(ndS)){ + _error_("convergence criterion is NaN (RSL_old=" << RSLg_old->Norm(NORM_TWO) << " RSL=" << RSLg->Norm(NORM_TWO) << ")"); + } if(!xIsNan(eps_rel)){ nS=RSLg_old->Norm(NORM_TWO); - if (xIsNan(nS)) _error_("convergence criterion is NaN!"); + if (xIsNan(nS)) _error_("convergence criterion is NaN! (check the initial RSL)"); }