Changeset 6953


Ignore:
Timestamp:
01/05/11 13:59:17 (14 years ago)
Author:
seroussi
Message:

added Incremental update of ice shelf geometry

Location:
issm/trunk/src/c
Files:
1 deleted
9 edited

Legend:

Unmodified
Added
Removed

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

    r6411 r6953  
    482482                                        ./modules/InputArtificialNoisex/InputArtificialNoisex.h\
    483483                                        ./modules/InputArtificialNoisex/InputArtificialNoisex.cpp\
    484                                         ./modules/UpdateGeometryx/UpdateGeometryx.h\
    485                                         ./modules/UpdateGeometryx/UpdateGeometryx.cpp\
    486484                                        ./modules/TimeAdaptx/TimeAdaptx.h\
    487485                                        ./modules/TimeAdaptx/TimeAdaptx.cpp\
    … …  
    10531051                                        ./modules/InputArtificialNoisex/InputArtificialNoisex.h\
    10541052                                        ./modules/InputArtificialNoisex/InputArtificialNoisex.cpp\
    1055                                         ./modules/UpdateGeometryx/UpdateGeometryx.h\
    1056                                         ./modules/UpdateGeometryx/UpdateGeometryx.cpp\
    10571053                                        ./modules/TimeAdaptx/TimeAdaptx.h\
    10581054                                        ./modules/TimeAdaptx/TimeAdaptx.cpp\

  • issm/trunk/src/c/modules/modules.h

    r6238 r6953  
    8989#include "./TriaSearchx/TriaSearchx.h"
    9090#include "./ThicknessAbsMisfitx/ThicknessAbsMisfitx.h"
    91 #include "./UpdateGeometryx/UpdateGeometryx.h"
    9291#include "./UpdateVertexPositionsx/UpdateVertexPositionsx.h"
    9392#include "./VerticesDofx/VerticesDofx.h"

  • issm/trunk/src/c/objects/Elements/Element.h

    r6238 r6953  
    5656                virtual void   DeleteResults(void)=0;
    5757                virtual void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type)=0;
    58                 virtual void   UpdateGeometry(void)=0;
    5958                virtual void   InputToResult(int enum_type,int step,double time)=0;
    6059                virtual void   ControlInputGetGradient(Vec gradient,int enum_type)=0;

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r6946 r6953  
    38773877        }
    38783878        else if (analysis_type==PrognosticAnalysisEnum){
    3879                 InputUpdateFromSolutionOneDofCollapsed(solution,ThicknessEnum);
     3879                InputUpdateFromSolutionPrognostic(solution);
    38803880        }
    38813881        else if (analysis_type==BalancedthicknessAnalysisEnum){
    … …  
    47074707        this->inputs->AddInput(new PentaVertexInput(AdjointyEnum,lambday));
    47084708
     4709        /*Free ressources:*/
     4710        xfree((void**)&doflist);
     4711}
     4712/*}}}*/
     4713/*FUNCTION Penta::InputUpdateFromSolutionPrognostic{{{1*/
     4714void  Penta::InputUpdateFromSolutionPrognostic(double* solution){
     4715
     4716        const int  numdof   = NDOF1*NUMVERTICES;
     4717        const int  numdof2d = NDOF1*NUMVERTICES2D;
     4718
     4719        int    i,dummy,hydroadjustment;
     4720        int*   doflist = NULL;
     4721        double values[numdof];
     4722        double rho_ice,rho_water;
     4723        double bed[numdof];
     4724        double surface[numdof];
     4725        double *bed_ptr = NULL;
     4726        double *thickness_ptr = NULL;
     4727        double *surface_ptr = NULL;
     4728        Penta  *penta   = NULL;
     4729
     4730        /*If not on bed, return*/
     4731        if (!IsOnBed()) return;
     4732
     4733        /*Get dof list: */
     4734        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     4735
     4736        /*Use the dof list to index into the solution vector and extrude it */
     4737        for(i=0;i<numdof2d;i++){
     4738                values[i]         =solution[doflist[i]];
     4739                values[i+numdof2d]=values[i];
     4740                if(isnan(values[i])) _error_("NaN found in solution vector");
     4741        }
     4742
     4743        /*Get previous bed, thickness and surface*/
     4744        Input* bed_input=inputs->GetInput(BedEnum);             _assert_(bed_input);
     4745        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     4746        Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
     4747        bed_input->GetValuesPtr(&bed_ptr,&dummy);
     4748        thickness_input->GetValuesPtr(&thickness_ptr,&dummy);
     4749        surface_input->GetValuesPtr(&surface_ptr,&dummy);
     4750
     4751        /*If shelf: hydrostatic equilibrium*/
     4752        if (this->IsOnShelf()){
     4753
     4754                /*Fing HydrostaticAdjustment to figure out how to update the geometry:*/
     4755                this->parameters->FindParam(&hydroadjustment,HydrostaticAdjustmentEnum);
     4756
     4757                /*recover material parameters: */
     4758                rho_ice=matpar->GetRhoIce();
     4759                rho_water=matpar->GetRhoWater();
     4760
     4761                if(hydroadjustment==AbsoluteEnum){
     4762                        for(i=0;i<numdof;i++) {
     4763                                surface[i]=values[i]*(1-rho_ice/rho_water);
     4764                                bed[i]=-values[i]*rho_ice/rho_water;
     4765                        }
     4766                }
     4767                else if(hydroadjustment==IncrementalEnum){
     4768                        for(i=0;i<numdof;i++) {
     4769                                surface[i]=surface_ptr[i]+(1.0-rho_ice/rho_water)*(values[i]-thickness_ptr[i]); //surface = oldsurface + (1-di) * dH
     4770                                bed[i]=bed_ptr[i]-rho_ice/rho_water*(values[i]-thickness_ptr[i]); //bed = oldbed + di * dH
     4771                        }
     4772                }
     4773                else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToString(hydroadjustment));
     4774        }
     4775
     4776        /*If sheet: surface = bed + thickness*/
     4777        else{
     4778                /*Now Compute surface only*/
     4779                for(i=0;i<numdof;i++) {
     4780                        surface[i]=bed_ptr[i]+values[i]; //surface=oldbed+newthickness
     4781                }
     4782        }
     4783
     4784        /*Start looping over all elements above current element and update all inputs*/
     4785        penta=this;
     4786        for(;;){
     4787                /*Add input to the element: */
     4788                penta->inputs->AddInput(new PentaVertexInput(ThicknessEnum,values));
     4789                penta->inputs->AddInput(new PentaVertexInput(SurfaceEnum,surface));
     4790                if (penta->IsOnShelf()) penta->inputs->AddInput(new PentaVertexInput(BedEnum,bed));
     4791
     4792                /*Stop if we have reached the surface*/
     4793                if (penta->IsOnSurface()) break;
     4794
     4795                /* get upper Penta*/
     4796                penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
     4797        }
     4798       
    47094799        /*Free ressources:*/
    47104800        xfree((void**)&doflist);
    … …  
    58015891}
    58025892/*}}}*/
    5803 /*FUNCTION Penta::UpdateGeometry{{{1*/
    5804 void  Penta::UpdateGeometry(void){
    5805 
    5806         /*Intermediaries*/
    5807         int hydroadjustment;
    5808         double rho_ice,rho_water;
    5809 
    5810         /*Correct thickness*/
    5811         this->inputs->ConstrainMin(ThicknessEnum,1.0);
    5812 
    5813         /*If shelf: hydrostatic equilibrium*/
    5814         if (this->IsOnShelf()){
    5815 
    5816                 /*Fing HydrostaticAdjustment to figure out how to update the geometry:*/
    5817                 this->parameters->FindParam(&hydroadjustment,HydrostaticAdjustmentEnum);
    5818 
    5819                 /*recover material parameters: */
    5820                 rho_ice=matpar->GetRhoIce();
    5821                 rho_water=matpar->GetRhoWater();
    5822 
    5823                 if(hydroadjustment==AbsoluteEnum){
    5824                         /*Create New Surface: s = (1-rho_ice/rho_water) h*/
    5825                         this->inputs->DuplicateInput(ThicknessEnum,SurfaceEnum);     //1: copy thickness into surface
    5826                         InputScale(SurfaceEnum,(1-rho_ice/rho_water)); //2: surface = surface * (1-di)
    5827 
    5828                         /*Create New Bed b = -rho_ice/rho_water h*/
    5829                         this->inputs->DuplicateInput(ThicknessEnum,BedEnum);         //1: copy thickness into bed
    5830                         InputScale(BedEnum, -rho_ice/rho_water);       //2: bed = bed * (-di)
    5831                 }
    5832                 else if(hydroadjustment==IncrementalEnum){
    5833                         /*The bed is update with di*Thickness*/
    5834                         this->inputs->AXPY(SurfaceEnum,1.0-rho_ice/rho_water,ThicknessEnum);     //surface = surface + (1-di) * thickness
    5835 
    5836                         /*The surface is update with (1-di)*Thickness*/
    5837                         this->inputs->AXPY(BedEnum,rho_ice/rho_water,ThicknessEnum);     //bed = bde + di * thickness
    5838                 }
    5839                 else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToString(hydroadjustment));
    5840         }
    5841 
    5842         /*If sheet: surface = bed + thickness*/
    5843         else{
    5844 
    5845                 /*The bed does not change, update surface only s = b + h*/
    5846                 this->inputs->DuplicateInput(BedEnum,SurfaceEnum);          //1: copy bed into surface
    5847                 this->inputs->AXPY(SurfaceEnum,1.0,ThicknessEnum);     //2: surface = surface + 1 * thickness
    5848         }
    5849 }
    5850 /*}}}*/

  • issm/trunk/src/c/objects/Elements/Penta.h

    r6723 r6953  
    121121                double SurfaceArea(void);
    122122                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
    123                 void   UpdateGeometry(void);
    124123                double TimeAdapt();
    125124                /*}}}*/
    … …  
    220219                void    InputUpdateFromSolutionDiagnosticVert( double* solutiong);
    221220                void    InputUpdateFromSolutionDiagnosticStokes( double* solutiong);
     221                void    InputUpdateFromSolutionPrognostic(double* solutiong);
    222222                void    InputUpdateFromSolutionThermal( double* solutiong);
    223223                void    InputUpdateFromSolutionOneDof(double* solutiong,int enum_type);

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r6946 r6953  
    36963696                        break;
    36973697                case PrognosticAnalysisEnum:
    3698                         InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
     3698                        InputUpdateFromSolutionPrognostic(solution);
    36993699                        break;
    37003700                case BalancedthicknessAnalysisEnum:
    … …  
    39303930        /*Add input to the element: */
    39313931        this->inputs->AddInput(new TriaVertexInput(enum_type,values));
     3932
     3933        /*Free ressources:*/
     3934        xfree((void**)&doflist);
     3935}
     3936/*}}}*/
     3937/*FUNCTION Tria::InputUpdateFromSolutionPrognostic{{{1*/
     3938void  Tria::InputUpdateFromSolutionPrognostic(double* solution){
     3939
     3940        /*Intermediaries*/
     3941        const int numdof = NDOF1*NUMVERTICES;
     3942
     3943        int       i,dummy,hydroadjustment;
     3944        int*      doflist=NULL;
     3945        double    rho_ice,rho_water;
     3946        double    values[numdof];
     3947        double    bed[numdof];
     3948        double    surface[numdof];
     3949        double    *bed_ptr = NULL;
     3950        double    *thickness_ptr = NULL;
     3951        double    *surface_ptr = NULL;
     3952
     3953        /*Get dof list: */
     3954        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     3955
     3956        /*Use the dof list to index into the solution vector: */
     3957        for(i=0;i<numdof;i++){
     3958                values[i]=solution[doflist[i]];
     3959                if(isnan(values[i])) _error_("NaN found in solution vector");
     3960        }
     3961
     3962        /*Get previous bed, thickness and surface*/
     3963        Input* bed_input=inputs->GetInput(BedEnum);             _assert_(bed_input);
     3964        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     3965        Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
     3966        bed_input->GetValuesPtr(&bed_ptr,&dummy);
     3967        thickness_input->GetValuesPtr(&thickness_ptr,&dummy);
     3968        surface_input->GetValuesPtr(&surface_ptr,&dummy);
     3969
     3970        /*If shelf: hydrostatic equilibrium*/
     3971        if (this->IsOnShelf()){
     3972
     3973                /*Fing HydrostaticAdjustment to figure out how to update the geometry:*/
     3974                this->parameters->FindParam(&hydroadjustment,HydrostaticAdjustmentEnum);
     3975
     3976                /*recover material parameters: */
     3977                rho_ice=matpar->GetRhoIce();
     3978                rho_water=matpar->GetRhoWater();
     3979
     3980                if(hydroadjustment==AbsoluteEnum){
     3981                        for(i=0;i<numdof;i++) {
     3982                                surface[i]=values[i]*(1-rho_ice/rho_water);
     3983                                bed[i]=-values[i]*rho_ice/rho_water;
     3984                        }
     3985                }
     3986                else if(hydroadjustment==IncrementalEnum){
     3987                        for(i=0;i<numdof;i++) {
     3988                                surface[i]=surface_ptr[i]+(1.0-rho_ice/rho_water)*(values[i]-thickness_ptr[i]); //surface = oldsurface + (1-di) * dH
     3989                                bed[i]=bed_ptr[i]-rho_ice/rho_water*(values[i]-thickness_ptr[i]); //bed = oldbed + di * dH
     3990                        }
     3991                }
     3992                else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToString(hydroadjustment));
     3993        }
     3994
     3995        /*If sheet: surface = bed + thickness*/
     3996        else{
     3997                /*Now Compute surface only*/
     3998                for(i=0;i<numdof;i++) {
     3999                        surface[i]=bed_ptr[i]+values[i]; //surface=oldbed+newthickness
     4000                }
     4001        }
     4002
     4003        /*Add input to the element: */
     4004        this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,values));
     4005        this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,surface));
     4006        if (this->IsOnShelf()) this->inputs->AddInput(new TriaVertexInput(BedEnum,bed));
    39324007
    39334008        /*Free ressources:*/
    … …  
    51255200}
    51265201/*}}}*/
    5127 /*FUNCTION Tria::UpdateGeometry{{{1*/
    5128 void  Tria::UpdateGeometry(void){
    5129 
    5130         /*Intermediaries*/
    5131         int hydroadjustment;
    5132         double rho_ice,rho_water;
    5133 
    5134         /*Correct thickness*/
    5135         this->inputs->ConstrainMin(ThicknessEnum,1.0);
    5136 
    5137         /*If shelf: hydrostatic equilibrium*/
    5138         if (this->IsOnShelf()){
    5139 
    5140                 /*Fing HydrostaticAdjustment to figure out how to update the geometry:*/
    5141                 this->parameters->FindParam(&hydroadjustment,HydrostaticAdjustmentEnum);
    5142 
    5143                 /*recover material parameters: */
    5144                 rho_ice=matpar->GetRhoIce();
    5145                 rho_water=matpar->GetRhoWater();
    5146 
    5147                 if(hydroadjustment==AbsoluteEnum){
    5148                         /*Create New Surface: s = (1-rho_ice/rho_water) h*/
    5149                         this->inputs->DuplicateInput(ThicknessEnum,SurfaceEnum);     //1: copy thickness into surface
    5150                         InputScale(SurfaceEnum,(1.0-rho_ice/rho_water)); //2: surface = surface * (1-di)
    5151 
    5152                         /*Create New Bed b = -rho_ice/rho_water h*/
    5153                         this->inputs->DuplicateInput(ThicknessEnum,BedEnum);         //1: copy thickness into bed
    5154                         InputScale(BedEnum, -rho_ice/rho_water);       //2: bed = bed * (-di)
    5155                 }
    5156                 else if(hydroadjustment==IncrementalEnum){
    5157                         /*The bed is update with di*Thickness*/
    5158                         this->inputs->AXPY(SurfaceEnum,1.0-rho_ice/rho_water,ThicknessEnum);     //surface = surface + (1-di) * thickness
    5159 
    5160                         /*The surface is update with (1-di)*Thickness*/
    5161                         this->inputs->AXPY(BedEnum,rho_ice/rho_water,ThicknessEnum);     //bed = bde + di * thickness
    5162                 }
    5163                 else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToString(hydroadjustment));
    5164         }
    5165 
    5166         /*If sheet: surface = bed + thickness*/
    5167         else{
    5168                
    5169                 /*The bed does not change, update surface only s = b + h*/
    5170                 this->inputs->DuplicateInput(BedEnum,SurfaceEnum);          //1: copy bed into surface
    5171                 this->inputs->AXPY(SurfaceEnum,1.0,ThicknessEnum);     //2: surface = surface + 1 * thickness
    5172         }
    5173 }
    5174 /*}}}*/

  • issm/trunk/src/c/objects/Elements/Tria.h

    r6412 r6953  
    124124                double SurfaceArea(void);
    125125                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
    126                 void   UpdateGeometry(void);
    127126                double TimeAdapt();
    128127                /*}}}*/
    … …  
    179178                void      InputUpdateFromSolutionDiagnosticHutter( double* solution);
    180179                void      InputUpdateFromSolutionOneDof(double* solution,int enum_type);
     180                void      InputUpdateFromSolutionPrognostic(double* solution);
    181181                bool      IsInput(int name);
    182182                void      SetClone(int* minranks);

  • issm/trunk/src/c/solutions/transient2d_core.cpp

    r6933 r6953  
    5959                prognostic_core(femmodel);
    6060       
    61                 _printf_(VerboseSolution(),"%s\n","   updating geometry");
    62                 UpdateGeometryx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
    63 
    6461                if(solution_type==Transient2DSolutionEnum && !control_analysis && (step%output_frequency==0)){
    6562                        _printf_(VerboseSolution(),"%s\n","   saving results\n");

  • issm/trunk/src/c/solutions/transient3d_core.cpp

    r6933 r6953  
    6262                prognostic_core(femmodel);
    6363       
    64                 _printf_(VerboseSolution(),"   updating geometry\n");
    65                 UpdateGeometryx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
    66                
    6764                _printf_(VerboseSolution(),"%s\n","   updating vertices positions");
    6865                UpdateVertexPositionsx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
Note: See TracChangeset for help on using the changeset viewer.