Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 10140) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 10141) @@ -1745,13 +1745,13 @@ const int numdof2d = NDOF1*NUMVERTICES2D; - int i,dummy,hydroadjustment; + int i,hydroadjustment; int* doflist = NULL; - double values[numdof]; double rho_ice,rho_water; - double bed[numdof]; - double surface[numdof]; - double *bed_ptr = NULL; - double *thickness_ptr = NULL; - double *surface_ptr = NULL; + double newthickness[numdof]; + double newbed[numdof]; + double newsurface[numdof]; + double oldbed[NUMVERTICES]; + double oldsurface[NUMVERTICES]; + double oldthickness[NUMVERTICES]; Penta *penta = NULL; @@ -1764,18 +1764,15 @@ /*Use the dof list to index into the solution vector and extrude it */ for(i=0;iGetInput(BedEnum); _assert_(bed_input); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); - bed_input->GetValuesPtr(&bed_ptr,&dummy); - thickness_input->GetValuesPtr(&thickness_ptr,&dummy); - surface_input->GetValuesPtr(&surface_ptr,&dummy); + GetInputListOnVertices(&oldbed[0],BedEnum); + GetInputListOnVertices(&oldsurface[0],SurfaceEnum); + GetInputListOnVertices(&oldthickness[0],ThicknessEnum); /*Fing PrognosticHydrostaticAdjustment to figure out how to update the geometry:*/ @@ -1789,15 +1786,15 @@ /*If shelf: hydrostatic equilibrium*/ if (this->nodes[i]->IsOnSheet()){ - surface[i]=bed_ptr[i]+values[i]; //surface=oldbed+newthickness - bed[i]=bed_ptr[i]; //bed does not change + newsurface[i]=oldbed[i]+newthickness[i]; //surface = oldbed + newthickness + newbed[i]=oldbed[i]; //same bed: do nothing } else{ //so it is an ice shelf if(hydroadjustment==AbsoluteEnum){ - surface[i]=values[i]*(1-rho_ice/rho_water); - bed[i]=values[i]*(-rho_ice/rho_water); + newsurface[i]=newthickness[i]*(1-rho_ice/rho_water); + newbed[i]=newthickness[i]*(-rho_ice/rho_water); } else if(hydroadjustment==IncrementalEnum){ - surface[i]=surface_ptr[i]+(1.0-rho_ice/rho_water)*(values[i]-thickness_ptr[i]); //surface = oldsurface + (1-di) * dH - bed[i]=bed_ptr[i]-rho_ice/rho_water*(values[i]-thickness_ptr[i]); //bed = oldbed + di * dH + newsurface[i]=oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH + newbed[i]=oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed = oldbed + di * dH } else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToStringx(hydroadjustment)); @@ -1809,7 +1806,7 @@ for(;;){ /*Add input to the element: */ - penta->inputs->AddInput(new PentaVertexInput(ThicknessEnum,values)); - penta->inputs->AddInput(new PentaVertexInput(SurfaceEnum,surface)); - penta->inputs->AddInput(new PentaVertexInput(BedEnum,bed)); + penta->inputs->AddInput(new PentaVertexInput(ThicknessEnum,newthickness)); + penta->inputs->AddInput(new PentaVertexInput(SurfaceEnum,newsurface)); + penta->inputs->AddInput(new PentaVertexInput(BedEnum,newbed)); /*Stop if we have reached the surface*/ @@ -7097,5 +7094,5 @@ const int numdof=NDOF2*NUMVERTICES; - int i,dummy; + int i; double rho_ice,g; double values[numdof]; @@ -7108,5 +7105,4 @@ double xyz_list[NUMVERTICES][3]; int *doflist = NULL; - double *vz_ptr = NULL; Penta *penta = NULL; @@ -7141,11 +7137,5 @@ /*Now Compute vel*/ - Input* vz_input=inputs->GetInput(VzEnum); - if (vz_input){ - if (vz_input->ObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as Vz is of type %s",EnumToStringx(vz_input->ObjectEnum())); - vz_input->GetValuesPtr(&vz_ptr,&dummy); - for(i=0;iGetInput(VzEnum); - if (vz_input){ - if (vz_input->ObjectEnum()!=PentaVertexInputEnum){ - _error_("Cannot compute Vel as Vz is of type %s",EnumToStringx(vz_input->ObjectEnum())); - } - vz_input->GetValuesPtr(&vz_ptr,&dummy); - for(i=0;iObjectEnum())); } - vzmacayeal_input->GetValuesPtr(&vzmacayeal_ptr,&dummy); - for(i=0;iObjectEnum())); } - vzpattyn_input->GetValuesPtr(&vzpattyn_ptr,&dummy); - for(i=0;iGetInput(VzEnum); - if (vz_input){ - if (vz_input->ObjectEnum()!=PentaVertexInputEnum){ - _error_("Cannot compute Vel as Vz is of type %s",EnumToStringx(vz_input->ObjectEnum())); - } - vz_input->GetValuesPtr(&vz_ptr,&dummy); - for(i=0;iinputs->AddInput(new PentaVertexInput(VxEnum,vx)); this->inputs->AddInput(new PentaVertexInput(VyEnum,vy)); - this->inputs->AddInput(new TriaVertexInput(VelEnum,vel)); + this->inputs->AddInput(new PentaVertexInput(VelEnum,vel)); this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure)); @@ -7612,5 +7572,5 @@ const int numdof=NDOF1*NUMVERTICES; - int i,dummy; + int i; int approximation; double rho_ice,g; @@ -7627,7 +7587,4 @@ double xyz_list[NUMVERTICES][3]; int* doflist = NULL; - double* vx_ptr = NULL; - double* vy_ptr = NULL; - double* vzstokes_ptr = NULL; @@ -7652,19 +7609,6 @@ /*Get Vx and Vy*/ - Input* vx_input=inputs->GetInput(VxEnum); - if (vx_input){ - if (vx_input->ObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as Vx is of type %s",EnumToStringx(vx_input->ObjectEnum())); - vx_input->GetValuesPtr(&vx_ptr,&dummy); - for(i=0;iGetInput(VyEnum); - if (vy_input){ - if (vy_input->ObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as Vy is of type %s",EnumToStringx(vy_input->ObjectEnum())); - vy_input->GetValuesPtr(&vy_ptr,&dummy); - for(i=0;iGetInput(VzStokesEnum); if (vzstokes_input){ - if (vzstokes_input->ObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as VzStokes is of type %s",EnumToStringx(vy_input->ObjectEnum())); - vzstokes_input->GetValuesPtr(&vzstokes_ptr,&dummy); - for(i=0;iObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as VzStokes is of type %s",EnumToStringx(vzstokes_input->ObjectEnum())); + GetInputListOnVertices(&vzstokes[0],VzStokesEnum); } else _error_("Cannot compute Vz as VzStokes in not present in PattynStokes element"); @@ -7685,7 +7628,6 @@ Input* vzstokes_input=inputs->GetInput(VzStokesEnum); if (vzstokes_input){ - if (vzstokes_input->ObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as VzStokes is of type %s",EnumToStringx(vy_input->ObjectEnum())); - vzstokes_input->GetValuesPtr(&vzstokes_ptr,&dummy); - for(i=0;iObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as VzStokes is of type %s",EnumToStringx(vzstokes_input->ObjectEnum())); + GetInputListOnVertices(&vzstokes[0],VzStokesEnum); } else _error_("Cannot compute Vz as VzStokes in not present in MacAyealStokes element"); @@ -7722,5 +7664,4 @@ this->inputs->AddInput(new PentaVertexInput(VzMacAyealEnum,vzmacayeal)); } - this->inputs->AddInput(new PentaVertexInput(VzEnum,vz)); this->inputs->AddInput(new PentaVertexInput(VelEnum,vel)); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 10140) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 10141) @@ -333,7 +333,6 @@ /*Surface and bed are updated. Update inputs:*/ - surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i]; - bed_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=b[i]; - + this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,&s[0])); + this->inputs->AddInput(new TriaVertexInput(BedEnum,&b[0])); } /*}}}*/ @@ -1651,13 +1650,13 @@ const int numdof = NDOF1*NUMVERTICES; - int i,dummy,hydroadjustment; + int i,hydroadjustment; int* doflist=NULL; double rho_ice,rho_water; - double values[numdof]; - double bed[numdof]; - double surface[numdof]; - double *bed_ptr = NULL; - double *thickness_ptr = NULL; - double *surface_ptr = NULL; + double newthickness[numdof]; + double newbed[numdof]; + double newsurface[numdof]; + double oldbed[NUMVERTICES]; + double oldsurface[NUMVERTICES]; + double oldthickness[NUMVERTICES]; /*Get dof list: */ @@ -1666,22 +1665,17 @@ /*Use the dof list to index into the solution vector: */ for(i=0;iGetInput(BedEnum); _assert_(bed_input); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); - bed_input->GetValuesPtr(&bed_ptr,&dummy); - thickness_input->GetValuesPtr(&thickness_ptr,&dummy); - surface_input->GetValuesPtr(&surface_ptr,&dummy); + GetInputListOnVertices(&oldbed[0],BedEnum); + GetInputListOnVertices(&oldsurface[0],SurfaceEnum); + GetInputListOnVertices(&oldthickness[0],ThicknessEnum); /*Fing PrognosticHydrostaticAdjustment to figure out how to update the geometry:*/ this->parameters->FindParam(&hydroadjustment,PrognosticHydrostaticAdjustmentEnum); - - /*recover material parameters: */ rho_ice=matpar->GetRhoIce(); rho_water=matpar->GetRhoWater(); @@ -1690,16 +1684,16 @@ /*If shelf: hydrostatic equilibrium*/ if (this->nodes[i]->IsOnSheet()){ - surface[i]=bed_ptr[i]+values[i]; //surface=oldbed+newthickness - bed[i]=bed_ptr[i]; //do nothing + newsurface[i]=oldbed[i]+newthickness[i]; //surface = oldbed + newthickness + newbed[i]=oldbed[i]; //same bed: do nothing } else{ //this is an ice shelf if(hydroadjustment==AbsoluteEnum){ - surface[i]=values[i]*(1-rho_ice/rho_water); - bed[i]=values[i]*(-rho_ice/rho_water); + newsurface[i]=newthickness[i]*(1-rho_ice/rho_water); + newbed[i]=newthickness[i]*(-rho_ice/rho_water); } else if(hydroadjustment==IncrementalEnum){ - surface[i]=surface_ptr[i]+(1.0-rho_ice/rho_water)*(values[i]-thickness_ptr[i]); //surface = oldsurface + (1-di) * dH - bed[i]=bed_ptr[i]-rho_ice/rho_water*(values[i]-thickness_ptr[i]); //bed = oldbed + di * dH + newsurface[i]=oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH + newbed[i]=oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed = oldbed + di * dH } else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToStringx(hydroadjustment)); @@ -1708,7 +1702,7 @@ /*Add input to the element: */ - this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,values)); - this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,surface)); - this->inputs->AddInput(new TriaVertexInput(BedEnum,bed)); + this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,newthickness)); + this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,newsurface)); + this->inputs->AddInput(new TriaVertexInput(BedEnum,newbed)); /*Free ressources:*/ @@ -1983,5 +1977,4 @@ void Tria::MigrateGroundingLine(void){ - double *values = NULL; double h[3],s[3],b[3],ba[3]; @@ -1992,10 +1985,5 @@ bool elementonshelf = false; - /*Recover info at the vertices: */ - Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input); - Input* bed_input =inputs->GetInput(BedEnum); _assert_(bed_input); - if((surface_input->ObjectEnum()!=TriaVertexInputEnum) | (bed_input->ObjectEnum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!"); - GetInputListOnVertices(&h[0],ThicknessEnum); GetInputListOnVertices(&s[0],SurfaceEnum); @@ -2020,9 +2008,7 @@ b[i]=ba[i]; s[i]=b[i]+h[i]; - if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i grounding %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]); } else{ /*do nothing, we are still floating.*/ - if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i still floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]); } } @@ -2035,10 +2021,7 @@ s[i]=(1-density)*h[i]; b[i]=-density*h[i]; - printf("%i\n",nodes[i]->Sid()+1); - if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]); } else{ /*do nothing, we are still grounded.*/ - if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i still grounded %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]); } } @@ -2046,11 +2029,8 @@ /*Surface and bed are updated. Update inputs:*/ - surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i]; - bed_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=b[i]; + this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,&s[0])); + this->inputs->AddInput(new TriaVertexInput(BedEnum,&b[0])); - for(i=0;i<3;i++){ - isonshelf[i]=nodes[i]->IsOnShelf(); - if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i second shelf status %i\n",this->Id(),nodes[i]->Sid()+1,isonshelf[i]); - } + for(i=0;i<3;i++) isonshelf[i]=nodes[i]->IsOnShelf(); } /*}}}*/ @@ -2207,4 +2187,64 @@ void Tria::ShelfSync(void){ + double melting[NUMVERTICES]; + double h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES]; + double bed_hydro; + double rho_water,rho_ice,density; + int i; + bool elementonshelf = false; + + /*melting rate at the grounding line: */ + double yts; + int swap; + double gl_melting_rate; + + /*recover parameters: */ + parameters->FindParam(&yts,ConstantsYtsEnum); + parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum); + + /*Recover info at the vertices: */ + GetInputListOnVertices(&h[0],ThicknessEnum); + GetInputListOnVertices(&s[0],SurfaceEnum); + GetInputListOnVertices(&b[0],BedEnum); + GetInputListOnVertices(&ba[0],BathymetryEnum); + + /*material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + density=rho_ice/rho_water; + + /*go through vertices, and update inputs, considering them to be TriaVertex type: */ + for(i=0;iinputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false)); + nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true)); + } + else{ + nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true)); + nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false)); + } + } + + /*Now, update shelf status of element. An element can only be on shelf if all its nodes are on shelf: */ + swap=0; + elementonshelf=false; + for(i=0;iIsOnShelf()){ + elementonshelf=true; + break; + } + } + if(!this->IsOnShelf() && elementonshelf==true)swap=1; + this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,elementonshelf)); + + /*If this element just became ungrounded, set its basal melting rate at 50 m/yr:*/ + if(swap){ + for(i=0;iinputs->AddInput(new TriaVertexInput(BasalforcingsMeltingRateEnum,&melting[0])); + } +} +/*}}}*/ +/*FUNCTION Tria::SoftMigration{{{1*/ +void Tria::SoftMigration(double* sheet_ungrounding){ double *values = NULL; @@ -2215,79 +2255,5 @@ bool elementonshelf = false; - /*melting rate at the grounding line: */ - double yts; - int swap; - double gl_melting_rate; - - /*recover parameters: */ - parameters->FindParam(&yts,ConstantsYtsEnum); - parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum); - /*Recover info at the vertices: */ - Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input); - Input* bed_input =inputs->GetInput(BedEnum); _assert_(bed_input); - if((surface_input->ObjectEnum()!=TriaVertexInputEnum) | (bed_input->ObjectEnum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!"); - - GetInputListOnVertices(&h[0],ThicknessEnum); - GetInputListOnVertices(&s[0],SurfaceEnum); - GetInputListOnVertices(&b[0],BedEnum); - GetInputListOnVertices(&ba[0],BathymetryEnum); - - /*material parameters: */ - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - density=rho_ice/rho_water; - - /*go through vertices, and update inputs, considering them to be TriaVertex type: */ - for(i=0;i<3;i++){ - if(b[i]==ba[i]){ - - nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false)); - nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true)); - } - else{ - nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true)); - nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false)); - - } - } - - /*Now, update shelf status of element. An element can only be on shelf if all its nodes are on shelf: */ - swap=0; - elementonshelf=false; - for(i=0;i<3;i++){ - if(nodes[i]->IsOnShelf()){ - elementonshelf=true; - break; - } - } - if(!this->IsOnShelf() && elementonshelf==true)swap=1; - this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,elementonshelf)); - - /*If this element just became ungrounded, set its basal melting rate at 50 m/yr:*/ - if(swap){ - Input* basal_melting_rate_input =inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_rate_input); - basal_melting_rate_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=gl_melting_rate/yts; - } - - -} -/*}}}*/ -/*FUNCTION Tria::SoftMigration{{{1*/ -void Tria::SoftMigration(double* sheet_ungrounding){ - - - double *values = NULL; - double h[3],s[3],b[3],ba[3]; - double bed_hydro; - double rho_water,rho_ice,density; - int i; - bool elementonshelf = false; - - /*Recover info at the vertices: */ - Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input); - Input* bed_input =inputs->GetInput(BedEnum); _assert_(bed_input); - if((surface_input->ObjectEnum()!=TriaVertexInputEnum) | (bed_input->ObjectEnum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!"); - GetInputListOnVertices(&h[0],ThicknessEnum); GetInputListOnVertices(&s[0],SurfaceEnum); @@ -2333,8 +2299,7 @@ /*Surface and bed are updated. Update inputs:*/ - surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i]; - bed_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=b[i]; - - } + this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,&s[0])); + this->inputs->AddInput(new TriaVertexInput(BedEnum,&b[0])); +} /*}}}*/ /*FUNCTION Tria::SurfaceArea {{{1*/ @@ -3244,9 +3209,9 @@ const int numdof=NDOF2*NUMVERTICES; - int i,dummy; - int* doflist=NULL; - double vx,vy; - double values[numdof]; - GaussTria* gauss=NULL; + int i; + double vx,vy; + double values[numdof]; + int *doflist = NULL; + GaussTria *gauss = NULL; /*Get dof list: */ @@ -3284,5 +3249,4 @@ int i; - int dummy; int* doflist=NULL; double rho_ice,g; @@ -3345,5 +3309,4 @@ int i; - int dummy; int* doflist=NULL; double rho_ice,g; @@ -3355,6 +3318,4 @@ double pressure[NUMVERTICES]; double thickness[NUMVERTICES]; - double* vz_ptr=NULL; - Input* vz_input=NULL; /*Get dof list: */ @@ -3374,18 +3335,6 @@ } - /*Get Vz*/ - vz_input=inputs->GetInput(VzEnum); - if (vz_input){ - if (vz_input->ObjectEnum()!=TriaVertexInputEnum){ - _error_("Cannot compute Vel as Vz is of type %s",EnumToStringx(vz_input->ObjectEnum())); - } - vz_input->GetValuesPtr(&vz_ptr,&dummy); - for(i=0;i