Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 6952) +++ /issm/trunk/src/c/Makefile.am (revision 6953) @@ -482,6 +482,4 @@ ./modules/InputArtificialNoisex/InputArtificialNoisex.h\ ./modules/InputArtificialNoisex/InputArtificialNoisex.cpp\ - ./modules/UpdateGeometryx/UpdateGeometryx.h\ - ./modules/UpdateGeometryx/UpdateGeometryx.cpp\ ./modules/TimeAdaptx/TimeAdaptx.h\ ./modules/TimeAdaptx/TimeAdaptx.cpp\ @@ -1053,6 +1051,4 @@ ./modules/InputArtificialNoisex/InputArtificialNoisex.h\ ./modules/InputArtificialNoisex/InputArtificialNoisex.cpp\ - ./modules/UpdateGeometryx/UpdateGeometryx.h\ - ./modules/UpdateGeometryx/UpdateGeometryx.cpp\ ./modules/TimeAdaptx/TimeAdaptx.h\ ./modules/TimeAdaptx/TimeAdaptx.cpp\ Index: /issm/trunk/src/c/modules/modules.h =================================================================== --- /issm/trunk/src/c/modules/modules.h (revision 6952) +++ /issm/trunk/src/c/modules/modules.h (revision 6953) @@ -89,5 +89,4 @@ #include "./TriaSearchx/TriaSearchx.h" #include "./ThicknessAbsMisfitx/ThicknessAbsMisfitx.h" -#include "./UpdateGeometryx/UpdateGeometryx.h" #include "./UpdateVertexPositionsx/UpdateVertexPositionsx.h" #include "./VerticesDofx/VerticesDofx.h" Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 6952) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 6953) @@ -56,5 +56,4 @@ virtual void DeleteResults(void)=0; virtual void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type)=0; - virtual void UpdateGeometry(void)=0; virtual void InputToResult(int enum_type,int step,double time)=0; virtual void ControlInputGetGradient(Vec gradient,int enum_type)=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6952) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6953) @@ -3877,5 +3877,5 @@ } else if (analysis_type==PrognosticAnalysisEnum){ - InputUpdateFromSolutionOneDofCollapsed(solution,ThicknessEnum); + InputUpdateFromSolutionPrognostic(solution); } else if (analysis_type==BalancedthicknessAnalysisEnum){ @@ -4707,4 +4707,94 @@ this->inputs->AddInput(new PentaVertexInput(AdjointyEnum,lambday)); + /*Free ressources:*/ + xfree((void**)&doflist); +} +/*}}}*/ +/*FUNCTION Penta::InputUpdateFromSolutionPrognostic{{{1*/ +void Penta::InputUpdateFromSolutionPrognostic(double* solution){ + + const int numdof = NDOF1*NUMVERTICES; + const int numdof2d = NDOF1*NUMVERTICES2D; + + int i,dummy,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; + Penta *penta = NULL; + + /*If not on bed, return*/ + if (!IsOnBed()) return; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*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); + + /*If shelf: hydrostatic equilibrium*/ + if (this->IsOnShelf()){ + + /*Fing HydrostaticAdjustment to figure out how to update the geometry:*/ + this->parameters->FindParam(&hydroadjustment,HydrostaticAdjustmentEnum); + + /*recover material parameters: */ + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + + if(hydroadjustment==AbsoluteEnum){ + for(i=0;iinputs->AddInput(new PentaVertexInput(ThicknessEnum,values)); + penta->inputs->AddInput(new PentaVertexInput(SurfaceEnum,surface)); + if (penta->IsOnShelf()) penta->inputs->AddInput(new PentaVertexInput(BedEnum,bed)); + + /*Stop if we have reached the surface*/ + if (penta->IsOnSurface()) break; + + /* get upper Penta*/ + penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id); + } + /*Free ressources:*/ xfree((void**)&doflist); @@ -5801,50 +5891,2 @@ } /*}}}*/ -/*FUNCTION Penta::UpdateGeometry{{{1*/ -void Penta::UpdateGeometry(void){ - - /*Intermediaries*/ - int hydroadjustment; - double rho_ice,rho_water; - - /*Correct thickness*/ - this->inputs->ConstrainMin(ThicknessEnum,1.0); - - /*If shelf: hydrostatic equilibrium*/ - if (this->IsOnShelf()){ - - /*Fing HydrostaticAdjustment to figure out how to update the geometry:*/ - this->parameters->FindParam(&hydroadjustment,HydrostaticAdjustmentEnum); - - /*recover material parameters: */ - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - - if(hydroadjustment==AbsoluteEnum){ - /*Create New Surface: s = (1-rho_ice/rho_water) h*/ - this->inputs->DuplicateInput(ThicknessEnum,SurfaceEnum); //1: copy thickness into surface - InputScale(SurfaceEnum,(1-rho_ice/rho_water)); //2: surface = surface * (1-di) - - /*Create New Bed b = -rho_ice/rho_water h*/ - this->inputs->DuplicateInput(ThicknessEnum,BedEnum); //1: copy thickness into bed - InputScale(BedEnum, -rho_ice/rho_water); //2: bed = bed * (-di) - } - else if(hydroadjustment==IncrementalEnum){ - /*The bed is update with di*Thickness*/ - this->inputs->AXPY(SurfaceEnum,1.0-rho_ice/rho_water,ThicknessEnum); //surface = surface + (1-di) * thickness - - /*The surface is update with (1-di)*Thickness*/ - this->inputs->AXPY(BedEnum,rho_ice/rho_water,ThicknessEnum); //bed = bde + di * thickness - } - else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToString(hydroadjustment)); - } - - /*If sheet: surface = bed + thickness*/ - else{ - - /*The bed does not change, update surface only s = b + h*/ - this->inputs->DuplicateInput(BedEnum,SurfaceEnum); //1: copy bed into surface - this->inputs->AXPY(SurfaceEnum,1.0,ThicknessEnum); //2: surface = surface + 1 * thickness - } -} -/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 6952) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 6953) @@ -121,5 +121,4 @@ double SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); - void UpdateGeometry(void); double TimeAdapt(); /*}}}*/ @@ -220,4 +219,5 @@ void InputUpdateFromSolutionDiagnosticVert( double* solutiong); void InputUpdateFromSolutionDiagnosticStokes( double* solutiong); + void InputUpdateFromSolutionPrognostic(double* solutiong); void InputUpdateFromSolutionThermal( double* solutiong); void InputUpdateFromSolutionOneDof(double* solutiong,int enum_type); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6952) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6953) @@ -3696,5 +3696,5 @@ break; case PrognosticAnalysisEnum: - InputUpdateFromSolutionOneDof(solution,ThicknessEnum); + InputUpdateFromSolutionPrognostic(solution); break; case BalancedthicknessAnalysisEnum: @@ -3930,4 +3930,79 @@ /*Add input to the element: */ this->inputs->AddInput(new TriaVertexInput(enum_type,values)); + + /*Free ressources:*/ + xfree((void**)&doflist); +} +/*}}}*/ +/*FUNCTION Tria::InputUpdateFromSolutionPrognostic{{{1*/ +void Tria::InputUpdateFromSolutionPrognostic(double* solution){ + + /*Intermediaries*/ + const int numdof = NDOF1*NUMVERTICES; + + int i,dummy,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; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*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); + + /*If shelf: hydrostatic equilibrium*/ + if (this->IsOnShelf()){ + + /*Fing HydrostaticAdjustment to figure out how to update the geometry:*/ + this->parameters->FindParam(&hydroadjustment,HydrostaticAdjustmentEnum); + + /*recover material parameters: */ + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + + if(hydroadjustment==AbsoluteEnum){ + for(i=0;iinputs->AddInput(new TriaVertexInput(ThicknessEnum,values)); + this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,surface)); + if (this->IsOnShelf()) this->inputs->AddInput(new TriaVertexInput(BedEnum,bed)); /*Free ressources:*/ @@ -5125,50 +5200,2 @@ } /*}}}*/ -/*FUNCTION Tria::UpdateGeometry{{{1*/ -void Tria::UpdateGeometry(void){ - - /*Intermediaries*/ - int hydroadjustment; - double rho_ice,rho_water; - - /*Correct thickness*/ - this->inputs->ConstrainMin(ThicknessEnum,1.0); - - /*If shelf: hydrostatic equilibrium*/ - if (this->IsOnShelf()){ - - /*Fing HydrostaticAdjustment to figure out how to update the geometry:*/ - this->parameters->FindParam(&hydroadjustment,HydrostaticAdjustmentEnum); - - /*recover material parameters: */ - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - - if(hydroadjustment==AbsoluteEnum){ - /*Create New Surface: s = (1-rho_ice/rho_water) h*/ - this->inputs->DuplicateInput(ThicknessEnum,SurfaceEnum); //1: copy thickness into surface - InputScale(SurfaceEnum,(1.0-rho_ice/rho_water)); //2: surface = surface * (1-di) - - /*Create New Bed b = -rho_ice/rho_water h*/ - this->inputs->DuplicateInput(ThicknessEnum,BedEnum); //1: copy thickness into bed - InputScale(BedEnum, -rho_ice/rho_water); //2: bed = bed * (-di) - } - else if(hydroadjustment==IncrementalEnum){ - /*The bed is update with di*Thickness*/ - this->inputs->AXPY(SurfaceEnum,1.0-rho_ice/rho_water,ThicknessEnum); //surface = surface + (1-di) * thickness - - /*The surface is update with (1-di)*Thickness*/ - this->inputs->AXPY(BedEnum,rho_ice/rho_water,ThicknessEnum); //bed = bde + di * thickness - } - else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToString(hydroadjustment)); - } - - /*If sheet: surface = bed + thickness*/ - else{ - - /*The bed does not change, update surface only s = b + h*/ - this->inputs->DuplicateInput(BedEnum,SurfaceEnum); //1: copy bed into surface - this->inputs->AXPY(SurfaceEnum,1.0,ThicknessEnum); //2: surface = surface + 1 * thickness - } -} -/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 6952) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 6953) @@ -124,5 +124,4 @@ double SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); - void UpdateGeometry(void); double TimeAdapt(); /*}}}*/ @@ -179,4 +178,5 @@ void InputUpdateFromSolutionDiagnosticHutter( double* solution); void InputUpdateFromSolutionOneDof(double* solution,int enum_type); + void InputUpdateFromSolutionPrognostic(double* solution); bool IsInput(int name); void SetClone(int* minranks); Index: /issm/trunk/src/c/solutions/transient2d_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient2d_core.cpp (revision 6952) +++ /issm/trunk/src/c/solutions/transient2d_core.cpp (revision 6953) @@ -59,7 +59,4 @@ prognostic_core(femmodel); - _printf_(VerboseSolution(),"%s\n"," updating geometry"); - UpdateGeometryx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); - if(solution_type==Transient2DSolutionEnum && !control_analysis && (step%output_frequency==0)){ _printf_(VerboseSolution(),"%s\n"," saving results\n"); Index: /issm/trunk/src/c/solutions/transient3d_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient3d_core.cpp (revision 6952) +++ /issm/trunk/src/c/solutions/transient3d_core.cpp (revision 6953) @@ -62,7 +62,4 @@ prognostic_core(femmodel); - _printf_(VerboseSolution()," updating geometry\n"); - UpdateGeometryx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); - _printf_(VerboseSolution(),"%s\n"," updating vertices positions"); UpdateVertexPositionsx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);