Changeset 7089


Ignore:
Timestamp:
01/14/11 08:14:31 (14 years ago)
Author:
Eric.Larour
Message:

Handle grounding line migration.

Location:
issm/trunk/src/c
Files:
18 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Container/Elements.cpp

    r6850 r7089  
    204204}
    205205/*}}}*/
     206/*FUNCTION Elements::NumberOfElements{{{1*/
     207int Elements::NumberOfElements(void){
     208
     209        int local_nelem=0;
     210        int numberofelements;
     211
     212        #ifdef _PARALLEL_
     213        local_nelem=this->Size();
     214        MPI_Allreduce ( (void*)&local_nelem,(void*)&numberofelements,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
     215        #else
     216        numberofelements=this->Size();
     217        #endif
     218
     219        return numberofelements;
     220}
     221/*}}}*/

  • issm/trunk/src/c/Container/Elements.h

    r6411 r7089  
    3333                void ToResults(Results* results,Parameters* parameters,int step, double time);
    3434                Patch* ResultsToPatch(void);
     35                int   NumberOfElements(void);
    3536                /*}}}*/
    3637

  • issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

    r6967 r7089  
    3636        /*Fetch data needed: */
    3737        IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
     38        IoModelFetchData(&iomodel->elementconnectivity,NULL,NULL,iomodel_handle,"elementconnectivity");
    3839        IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements");
    3940        IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements");
    … …  
    4849
    4950                        /*Create and add tria element to elements dataset: */
    50                         if(iomodel->dim==2)elements->AddObject(new Tria(i+1,i,iomodel,nummodels));
    51                         else elements->AddObject(new Penta(i+1,i,iomodel,nummodels));
     51                        if(iomodel->dim==2)elements->AddObject(new Tria(i+1,i,i,iomodel,nummodels));
     52                        else elements->AddObject(new Penta(i+1,i,i,iomodel,nummodels));
    5253
    5354                        /*Create and add material property to materials dataset: */
    … …  
    5859        /*Free data: */
    5960        xfree((void**)&iomodel->elements);
     61        xfree((void**)&iomodel->elementconnectivity);
    6062        xfree((void**)&iomodel->upperelements);
    6163        xfree((void**)&iomodel->lowerelements);

  • issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r6946 r7089  
    7575        parameters->AddObject(new BoolParam(KffEnum,iomodel->kff));
    7676        parameters->AddObject(new BoolParam(IoGatherEnum,iomodel->io_gather));
     77        parameters->AddObject(new BoolParam(GroundingLineMigrationEnum,iomodel->grounding_line_migration));
    7778
    7879        /*Deal with more complex parameters*/

  • issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp

    r6972 r7089  
    3939        IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
    4040        IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
     41        if(iomodel->grounding_line_migration) IoModelFetchData(&iomodel->bathymetry,NULL,NULL,iomodel_handle,"bathymetry");
    4142        if (iomodel->dim==3){
    4243                IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed");

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

    r6953 r7089  
    3232                virtual void   GetSolutionFromInputs(Vec solution)=0;
    3333                virtual int    GetNodeIndex(Node* node)=0;
     34                virtual int    Sid()=0;
    3435                virtual bool   IsOnShelf()=0;
     36                virtual bool   IsNodeOnShelf()=0;
    3537                virtual bool   IsOnBed()=0;
    3638                virtual void   GetParameterListOnVertices(double* pvalue,int enumtype)=0;
    … …  
    7981                virtual bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums)=0;
    8082                virtual void   AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part)=0;
     83                virtual int*   GetHorizontalNeighboorSids(void)=0;
    8184                virtual double TimeAdapt()=0;
     85                virtual void   MigrateGroundingLine()=0;
     86                virtual void   UpdateShelfStatus()=0;
    8287
    8388                /*Implementation: */

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

    r7073 r7089  
    2626/*FUNCTION Penta::Penta(){{{1*/
    2727Penta::Penta(){
     28
     29        int i;
     30
    2831        this->nodes=NULL;
    2932        this->matice=NULL;
    3033        this->matpar=NULL;
    31         this->neighbors=NULL;
     34        this->verticalneighbors=NULL;
    3235        this->inputs=NULL;
    3336        this->parameters=NULL;
    3437        this->results=NULL;
     38        for(i=0;i<3;i++)this->horizontalneighborsids[i]=UNDEF;
    3539}
    3640/*}}}*/
    … …  
    4347/*}}}*/
    4448/*FUNCTION Penta::Penta(int id, int index, IoModel* iomodel,int nummodels) {{{1*/
    45 Penta::Penta(int penta_id, int index, IoModel* iomodel,int nummodels)
     49Penta::Penta(int penta_id, int penta_sid, int index, IoModel* iomodel,int nummodels)
    4650        :PentaRef(nummodels)
    4751        ,PentaHook(nummodels,index+1,iomodel->numberofelements+1) //index+1: matice id, iomodel->numberofelements+1: matpar id
    … …  
    5963        /*id: */
    6064        this->id=penta_id;
     65        this->sid=penta_sid;
    6166
    6267        /*Build neighbors list*/
    … …  
    6671        else                                    penta_elements_ids[0]=(int)(iomodel->lowerelements[index]);
    6772        this->InitHookNeighbors(penta_elements_ids);
     73
     74        /*Build horizontalneighborsids list: */
     75        for (i=0;i<3;i++) this->horizontalneighborsids[i]=(int)iomodel->elementconnectivity[3*index+i]-1;
    6876
    6977        //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
    … …  
    7886        this->matice=NULL;
    7987        this->matpar=NULL;
    80         this->neighbors=NULL;
     88        this->verticalneighbors=NULL;
    8189}
    8290/*}}}*/
    … …  
    104112        /*deal with Penta  copy fields: */
    105113        penta->id=this->id;
     114        penta->sid=this->sid;
    106115        if(this->inputs){
    107116                penta->inputs=(Inputs*)this->inputs->Copy();
    … …  
    124133        penta->matice=(Matice*)penta->hmatice->delivers();
    125134        penta->matpar=(Matpar*)penta->hmatpar->delivers();
    126         penta->neighbors=(Penta**)penta->hneighbors->deliverp();
     135        penta->verticalneighbors=(Penta**)penta->hneighbors->deliverp();
     136
     137        /*neighbors: */
     138        for(i=0;i<3;i++)penta->horizontalneighborsids[i]=this->horizontalneighborsids[i];
    127139
    128140        return penta;
    … …  
    154166        /*marshall Penta data: */
    155167        memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
     168        memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid);
    156169        memcpy(marshalled_dataset,&numanalyses,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses);
    157170
    … …  
    195208        xfree((void**)&marshalled_results);
    196209
     210        /*marshall horizontal neighbors: */
     211        memcpy(marshalled_dataset,horizontalneighborsids,3*sizeof(int));marshalled_dataset+=3*sizeof(int);
     212
    197213        *pmarshalled_dataset=marshalled_dataset;
    198214        return;
    … …  
    211227
    212228        return sizeof(id)
     229                +sizeof(sid)
    213230                +hnodes_size
    214231                +sizeof(numanalyses)
    … …  
    219236                +inputs->MarshallSize()
    220237                +results->MarshallSize()
     238                +3*sizeof(int)
    221239                +sizeof(int); //sizeof(int) for enum type
    222240}
    … …  
    235253         *object data (thanks to DataSet::Demarshall):*/
    236254        memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
     255        memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid);
    237256        memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses);
    238257
    … …  
    260279        matice=NULL;
    261280        matpar=NULL;
    262         neighbors=NULL;
     281        verticalneighbors=NULL;
    263282
    264283        /*demarshall inputs and results: */
    … …  
    268287        /*parameters: may not exist even yet, so let Configure handle it: */
    269288        this->parameters=NULL;
     289
     290        /*neighbors: */
     291        memcpy(&this->horizontalneighborsids,marshalled_dataset,3*sizeof(int));marshalled_dataset+=3*sizeof(int);
    270292
    271293        /*return: */
    … …  
    425447        this->matice=(Matice*)this->hmatice->delivers();
    426448        this->matpar=(Matpar*)this->hmatpar->delivers();
    427         this->neighbors=(Penta**)this->hneighbors->deliverp();
     449        this->verticalneighbors=(Penta**)this->hneighbors->deliverp();
    428450
    429451        /*point parameters to real dataset: */
    … …  
    499521
    500522        /*Checks in debugging {{{2*/
    501         _assert_(this->nodes && this->matice && this->matpar && this->neighbors && this->parameters && this->inputs);
     523        _assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
    502524        /*}}}*/
    503525
    … …  
    17941816
    17951817        /*if debugging mode, check that all pointers exist {{{2*/
    1796         _assert_(this->nodes && this->matice && this->matpar && this->neighbors && this->parameters && this->inputs);
     1818        _assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
    17971819        /*}}}*/
    17981820
    … …  
    28782900        matice->DeepEcho();
    28792901        matpar->DeepEcho();
    2880         printf("   neighbor ids: %i-%i\n",neighbors[0]->Id(),neighbors[1]->Id());
     2902        printf("   neighbor ids: %i-%i\n",verticalneighbors[0]->Id(),verticalneighbors[1]->Id());
    28812903        printf("   parameters\n");
    28822904        parameters->DeepEcho();
    … …  
    28852907        printf("   results\n");
    28862908        results->DeepEcho();
     2909        printf("neighboor sids: \n");
     2910        printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]);
     2911
    28872912        return;
    28882913}
    … …  
    30633088}
    30643089/*}}}*/
     3090/*FUNCTION Penta::GetHorizontalNeighboorSids {{{1*/
     3091int* Penta::GetHorizontalNeighboorSids(){
     3092
     3093        /*return PentaRef field*/
     3094        return &this->horizontalneighborsids[0];
     3095
     3096}
     3097/*}}}*/
    30653098/*FUNCTION Penta::GetLowerElement{{{1*/
    30663099Penta* Penta::GetLowerElement(void){
    … …  
    30683101        Penta* upper_penta=NULL;
    30693102
    3070         upper_penta=(Penta*)neighbors[0]; //first one (0) under, second one (1) above
     3103        upper_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
    30713104
    30723105        return upper_penta;
    … …  
    34973530        Penta* upper_penta=NULL;
    34983531
    3499         upper_penta=(Penta*)neighbors[1]; //first one under, second one above
     3532        upper_penta=(Penta*)verticalneighbors[1]; //first one under, second one above
    35003533
    35013534        return upper_penta;
    … …  
    35333566
    35343567        return z;
     3568}
     3569/*}}}*/
     3570/*FUNCTION Penta::Sid {{{1*/
     3571int    Penta::Sid(){
     3572       
     3573        return sid;
     3574
    35353575}
    35363576/*}}}*/
    … …  
    51915231}
    51925232/*}}}*/
     5233/*FUNCTION Penta::IsNodeOnShelf {{{1*/
     5234bool   Penta::IsNodeOnShelf(){
     5235
     5236        int  i;
     5237        bool shelf=false;
     5238
     5239        for(i=0;i<6;i++){
     5240                if (nodes[i]->IsOnShelf()){
     5241                        shelf=true;
     5242                        break;
     5243                }
     5244        }
     5245        return shelf;
     5246}
     5247/*}}}*/
    51935248/*FUNCTION Penta::IsOnSurface{{{1*/
    51945249bool Penta::IsOnSurface(void){
    … …  
    53025357        /*Assign output pointers:*/
    53035358        *pmaxvz=maxvz;
     5359}
     5360/*}}}*/
     5361/*FUNCTION Penta::MigrateGroundingLine{{{1*/
     5362void  Penta::MigrateGroundingLine(void){
     5363        _error_("not supported yet!");
    53045364}
    53055365/*}}}*/
    … …  
    54175477        int     numrows  = 0;
    54185478        int     numnodes = 0;
     5479        int     temp_numnodes=0;
    54195480
    54205481        /*Go through all the results objects, and update the counters: */
    … …  
    54245485                numrows++;
    54255486                /*now, how many vertices and how many nodal values for this result? :*/
    5426                 numnodes=elementresult->NumberOfNodalValues(); //ask result object.
     5487                temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object.
     5488                if(temp_numnodes>numnodes)numnodes=temp_numnodes;
    54275489        }
    54285490
    … …  
    60716133}
    60726134/*}}}*/
     6135/*FUNCTION Penta::UpdateShelfStatus{{{1*/
     6136void Penta::UpdateShelfStatus(void){
     6137        _error_("Not implemented yet");
     6138}
     6139/*}}}*/

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

    r6987 r7089  
    3232
    3333                int          id;
     34                int          sid;
    3435
    3536                Node       **nodes;        // 6 nodes
    3637                Matice      *matice;       // 1 material ice
    3738                Matpar      *matpar;       // 1 material parameter
    38                 Penta      **neighbors;   // 2 neighbors: first one under, second one above
     39                Penta      **verticalneighbors;   // 2 neighbors: first one under, second one above
     40                int          horizontalneighborsids[3];
    3941
    4042                Parameters  *parameters;   //pointer to solution parameters
    … …  
    4446                /*Penta constructors and destructor: {{{1*/
    4547                Penta();
    46                 Penta(int penta_id,int i, IoModel* iomodel,int nummodels);
     48                Penta(int penta_id,int penta_sid,int i, IoModel* iomodel,int nummodels);
    4749                ~Penta();
    4850                /*}}}*/
    … …  
    8890                void   GradjB(Vec gradient);
    8991                void   GradjDrag(Vec gradient);
     92                int    Sid();
    9093                void   InputControlUpdate(double scalar,bool save_parameter);
    9194                void   InputArtificialNoise(int enum_type,double min, double max);
    … …  
    106109                void   MaxVy(double* pmaxvy, bool process_units);
    107110                void   MaxVz(double* pmaxvz, bool process_units);
     111                void   MigrateGroundingLine();
    108112                void   MinVel(double* pminvel, bool process_units);
    109113                void   MinVx(double* pminvx, bool process_units);
    … …  
    121125                double SurfaceArea(void);
    122126                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
     127                void   UpdateShelfStatus(void);
    123128                double TimeAdapt();
     129                int*   GetHorizontalNeighboorSids(void);
    124130                /*}}}*/
    125131                /*Penta specific routines:{{{1*/
    … …  
    229235                bool      IsOnBed(void);
    230236                bool    IsOnShelf(void);
     237                bool    IsNodeOnShelf(void);
    231238                bool    IsOnWater(void);
    232239                double  MinEdgeLength(double xyz_list[6][3]);

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

    r6971 r7089  
    2626Tria::Tria(){
    2727
     28        int i;
     29
    2830        this->nodes=NULL;
    2931        this->matice=NULL;
    3032        this->matpar=NULL;
     33        for(i=0;i<3;i++)this->horizontalneighborsids[i]=UNDEF;
    3134        this->inputs=NULL;
    3235        this->parameters=NULL;
    … …  
    3538}
    3639/*}}}*/
    37 /*FUNCTION Tria::Tria(int id, int index, IoModel* iomodel,int nummodels){{{1*/
    38 Tria::Tria(int tria_id, int index, IoModel* iomodel,int nummodels)
     40/*FUNCTION Tria::Tria(int id, int sid,int index, IoModel* iomodel,int nummodels){{{1*/
     41Tria::Tria(int tria_id, int tria_sid, int index, IoModel* iomodel,int nummodels)
    3942        :TriaRef(nummodels)
    4043        ,TriaHook(nummodels,index+1,iomodel->numberofelements+1){
    41         /*id: */
    42         this->id=tria_id;
    43 
    44         //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
    45         this->parameters=NULL;
    46 
    47         /*intialize inputs and results: */
    48         this->inputs=new Inputs();
    49         this->results=new Results();
    50 
    51         /*initialize pointers:*/
    52         this->nodes=NULL;
    53         this->matice=NULL;
    54         this->matpar=NULL;
     44               
     45                int i;
     46                /*id: */
     47                this->id=tria_id;
     48                this->sid=tria_sid;
     49
     50                //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
     51                this->parameters=NULL;
     52
     53                /*Build horizontalneighborsids list: */
     54                for (i=0;i<3;i++) this->horizontalneighborsids[i]=(int)iomodel->elementconnectivity[3*index+i]-1;
     55
     56                /*intialize inputs and results: */
     57                this->inputs=new Inputs();
     58                this->results=new Results();
     59
     60                /*initialize pointers:*/
     61                this->nodes=NULL;
     62                this->matice=NULL;
     63                this->matpar=NULL;
    5564
    5665}
    … …  
    8493        /*deal with Tria fields: */
    8594        tria->id=this->id;
     95        tria->sid=this->sid;
    8696        if(this->inputs){
    8797                tria->inputs=(Inputs*)this->inputs->Copy();
    … …  
    104114        tria->matice=(Matice*)tria->hmatice->delivers();
    105115        tria->matpar=(Matpar*)tria->hmatpar->delivers();
     116
     117        /*neighbors: */
     118        for(i=0;i<3;i++)tria->horizontalneighborsids[i]=this->horizontalneighborsids[i];
    106119
    107120        return tria;
    … …  
    133146        /*marshall Tria data: */
    134147        memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
     148        memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid);
    135149        memcpy(marshalled_dataset,&numanalyses,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses);
    136150
    … …  
    174188        xfree((void**)&marshalled_results);
    175189
     190        /*marshall horizontal neighbors: */
     191        memcpy(marshalled_dataset,horizontalneighborsids,3*sizeof(int));marshalled_dataset+=3*sizeof(int);
     192
    176193        *pmarshalled_dataset=marshalled_dataset;
    177194        return;
    … …  
    190207
    191208        return sizeof(id)
     209          +sizeof(sid)
    192210          +hnodes_size
    193211          +sizeof(numanalyses)
    … …  
    197215          +inputs->MarshallSize()
    198216          +results->MarshallSize()
     217          +3*sizeof(int)
    199218          +sizeof(int); //sizeof(int) for enum type
    200219}
    … …  
    213232         *object data (thanks to DataSet::Demarshall):*/
    214233        memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
     234        memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid);
    215235        memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses);
    216236
    … …  
    244264        /*parameters: may not exist even yet, so let Configure handle it: */
    245265        this->parameters=NULL;
     266
     267        /*neighbors: */
     268        memcpy(&this->horizontalneighborsids,marshalled_dataset,3*sizeof(int));marshalled_dataset+=3*sizeof(int);
    246269
    247270        /*return: */
    … …  
    26142637        if (results) results->DeepEcho();
    26152638        else printf("results=NULL\n");
     2639
     2640        printf("neighboor sids: \n");
     2641        printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]);
    26162642       
    26172643        return;
    … …  
    26542680        if (results) results->Echo();
    26552681        else printf("results=NULL\n");
     2682
     2683        printf("neighboor sids: \n");
     2684        printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]);
    26562685}
    26572686/*}}}*/
    … …  
    27152744        /*return TriaRef field*/
    27162745        return this->element_type;
     2746
     2747}
     2748/*}}}*/
     2749/*FUNCTION Tria::GetHorizontalNeighboorSids {{{1*/
     2750int* Tria::GetHorizontalNeighboorSids(){
     2751
     2752        /*return TriaRef field*/
     2753        return &this->horizontalneighborsids[0];
    27172754
    27182755}
    … …  
    33403377}
    33413378/*}}}*/
     3379/*FUNCTION Tria::Sid {{{1*/
     3380int    Tria::Sid(){
     3381       
     3382        return sid;
     3383
     3384}
     3385/*}}}*/
    33423386/*FUNCTION Tria::InputArtificialNoise{{{1*/
    33433387void  Tria::InputArtificialNoise(int enum_type,double min,double max){
    … …  
    35483592                for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_vertex_ids[i]-1];
    35493593                this->inputs->AddInput(new TriaVertexInput(BedEnum,nodeinputs));
     3594        }
     3595        if(iomodel->grounding_line_migration){
     3596                for(i=0;i<3;i++)nodeinputs[i]=iomodel->bathymetry[tria_vertex_ids[i]-1];
     3597                this->inputs->AddInput(new TriaVertexInput(BathymetryEnum,nodeinputs));
    35503598        }
    35513599        if (iomodel->drag_coefficient) {
    … …  
    41724220}
    41734221/*}}}*/
     4222/*FUNCTION Tria::IsNodeOnShelf {{{1*/
     4223bool   Tria::IsNodeOnShelf(){
     4224
     4225        int  i;
     4226        bool shelf=false;
     4227
     4228        for(i=0;i<3;i++){
     4229                if (nodes[i]->IsOnShelf()){
     4230                        shelf=true;
     4231                        break;
     4232                }
     4233        }
     4234        return shelf;
     4235}
     4236/*}}}*/
    41744237/*FUNCTION Tria::IsOnWater {{{1*/
    41754238bool   Tria::IsOnWater(){
    … …  
    43354398}
    43364399/*}}}*/
     4400/*FUNCTION Tria::MigrateGroundingLine{{{1*/
     4401void  Tria::MigrateGroundingLine(void){
     4402
     4403        int     i;
     4404
     4405        double  h[3];
     4406        double  s[3];
     4407        double  b[3];
     4408        double  ba[3];
     4409        double  bed_hydro;
     4410        int     isonshelf[3];
     4411        Input  *surface_input  = NULL;
     4412        Input  *bed_input      = NULL;
     4413        double  rho_water;
     4414        double  rho_ice;
     4415        double  density;
     4416        double *values         = NULL;
     4417        bool    elementonshelf = false;
     4418
     4419
     4420        /*Recover info at the vertices: */
     4421        surface_input   =inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
     4422        bed_input   =inputs->GetInput(BedEnum);   _assert_(bed_input);
     4423        if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
     4424
     4425        GetParameterListOnVertices(&h[0],ThicknessEnum);
     4426        GetParameterListOnVertices(&s[0],SurfaceEnum);
     4427        GetParameterListOnVertices(&b[0],BedEnum);
     4428        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     4429        for(i=0;i<3;i++)isonshelf[i]=nodes[i]->IsOnShelf();
     4430
     4431        /*material parameters: */
     4432        rho_water=matpar->GetRhoWater();
     4433        rho_ice=matpar->GetRhoIce();
     4434        density=rho_ice/rho_water;
     4435
     4436        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
     4437        for(i=0;i<3;i++){
     4438                if (isonshelf[i]){
     4439                        /*This node is on the shelf. See if its bed is going under the bathymetry: */
     4440                        if(b[i]<=ba[i]){ //<= because Neff being 0 when b=ba, drag will be 0 anyway.
     4441                                /*The ice shelf is getting grounded, the thickness is the same, so just update the bed to stick to the bathymetry and elevate the surface accordingly: */
     4442                                b[i]=ba[i];
     4443                                s[i]=b[i]+h[i];
     4444                        }
     4445                        else{
     4446                                /*do nothing, we are still floating.*/
     4447                        }
     4448                }
     4449                else{
     4450                        /*This node is on the sheet, near the grounding line. See if wants to unground. To
     4451                         * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */
     4452                        bed_hydro=-density*h[i];
     4453                        if (bed_hydro>ba[i]){
     4454                                /*We are now floating, bed and surface are determined from hydrostatic equilibrium: */
     4455                                s[i]=(1-density)*h[i];
     4456                                b[i]=-density*h[i];
     4457                        }
     4458                        else{
     4459                                /*do nothing, we are still grounded.*/
     4460                        }
     4461                }
     4462        }
     4463
     4464        /*Surface and bed are updated. Update inputs:*/
     4465        surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i];
     4466        bed_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=b[i];
     4467
     4468}
     4469/*}}}*/
    43374470/*FUNCTION Tria::MinVel{{{1*/
    43384471void  Tria::MinVel(double* pminvel, bool process_units){
    … …  
    44264559        int     numrows     = 0;
    44274560        int     numnodes    = 0;
     4561        int     temp_numnodes    = 0;
    44284562
    44294563        /*Go through all the results objects, and update the counters: */
    … …  
    44334567                numrows++;
    44344568                /*now, how many vertices and how many nodal values for this result? :*/
    4435                 numnodes=elementresult->NumberOfNodalValues(); //ask result object.
     4569                temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object.
     4570                if(temp_numnodes>numnodes)numnodes=temp_numnodes;
    44364571        }
    44374572
    … …  
    51965331}
    51975332/*}}}*/
     5333/*FUNCTION Tria::UpdateShelfStatus{{{1*/
     5334void  Tria::UpdateShelfStatus(void){
     5335
     5336        int     i;
     5337
     5338        double  h[3];
     5339        double  s[3];
     5340        double  b[3];
     5341        double  ba[3];
     5342        Input  *surface_input  = NULL;
     5343        Input  *bed_input      = NULL;
     5344        double  rho_water;
     5345        double  rho_ice;
     5346        double  density;
     5347        bool    elementonshelf = false;
     5348
     5349
     5350        /*Recover info at the vertices: */
     5351        surface_input   =inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
     5352        bed_input   =inputs->GetInput(BedEnum);   _assert_(bed_input);
     5353        if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
     5354
     5355        GetParameterListOnVertices(&h[0],ThicknessEnum);
     5356        GetParameterListOnVertices(&s[0],SurfaceEnum);
     5357        GetParameterListOnVertices(&b[0],BedEnum);
     5358        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     5359
     5360        /*material parameters: */
     5361        rho_water=matpar->GetRhoWater();
     5362        rho_ice=matpar->GetRhoIce();
     5363        density=rho_ice/rho_water;
     5364
     5365        /*go through vertices, and figure out if they are grounded or not, then update their status: */
     5366        for(i=0;i<3;i++){
     5367                if(b[i]<=ba[i]){
     5368                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,false));
     5369                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,true));
     5370                }
     5371                else{
     5372                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,true));
     5373                        nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,false));
     5374                }
     5375        }
     5376
     5377        /*Now, update  shelf status of element. An element can only be on shelf if all its nodes are on shelf: */
     5378        elementonshelf=true;
     5379        for(i=0;i<3;i++){
     5380                if(nodes[i]->IsOnSheet()){
     5381                        elementonshelf=false;
     5382                        break;
     5383                }
     5384        }
     5385    this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,elementonshelf));
     5386}
     5387/*}}}*/

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

    r6953 r7089  
    3030
    3131                int  id;
     32                int  sid;
    3233
    3334                Node   **nodes;    // 3 nodes
    3435                Matice  *matice;   // 1 material ice
    3536                Matpar  *matpar;   // 1 material parameter
     37                int      horizontalneighborsids[3];
    3638
    3739                Parameters *parameters;   //pointer to solution parameters
    … …  
    4143                /*Tria constructors, destructors {{{1*/
    4244                Tria();
    43                 Tria(int tria_id,int i, IoModel* iomodel,int nummodels);
     45                Tria(int tria_id,int tria_sid,int i, IoModel* iomodel,int nummodels);
    4446                ~Tria();
    4547                /*}}}*/
    … …  
    7880                void   CreatePVector(Vec pg, Vec pf);
    7981                int    GetNodeIndex(Node* node);
     82                int    Sid();
    8083                bool   IsOnBed();
    8184                bool   IsOnShelf();
     85                bool   IsNodeOnShelf();
    8286                bool   IsOnWater();
    8387                void   GetSolutionFromInputs(Vec solution);
    … …  
    109113                void   MaxVy(double* pmaxvy, bool process_units);
    110114                void   MaxVz(double* pmaxvz, bool process_units);
     115                void   MigrateGroundingLine();
    111116                void   MinVel(double* pminvel, bool process_units);
    112117                void   MinVx(double* pminvx, bool process_units);
    … …  
    124129                double SurfaceArea(void);
    125130                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
     131                void   UpdateShelfStatus(void);
    126132                double TimeAdapt();
     133                int*   GetHorizontalNeighboorSids(void);
    127134                /*}}}*/
    128135                /*Tria specific routines:{{{1*/

  • issm/trunk/src/c/objects/Inputs/BoolInput.cpp

    r6412 r7089  
    148148/*FUNCTION BoolInput::SpawnResult{{{1*/
    149149ElementResult* BoolInput::SpawnResult(int step, double time){
    150 
    151         _error_(" not supported yet!");
     150       
     151        return new BoolElementResult(this->enum_type,this->value,step,time);
    152152
    153153}

  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp

    r6413 r7089  
    399399
    400400        *pvalues=this->values;
    401         *pnum_values=3;
    402 
    403 }
    404 /*}}}*/
     401        if(pnum_values)*pnum_values=3;
     402
     403}
     404/*}}}*/

  • issm/trunk/src/c/objects/IoModel.cpp

    r6946 r7089  
    3838        xfree((void**)&this->z);
    3939        xfree((void**)&this->elements);
     40        xfree((void**)&this->elementconnectivity);
    4041        xfree((void**)&this->elements_type);
    4142        xfree((void**)&this->vertices_type);
    … …  
    5960        xfree((void**)&this->surface);
    6061        xfree((void**)&this->bed);
     62        xfree((void**)&this->bathymetry);
    6163        xfree((void**)&this->vx_obs);
    6264        xfree((void**)&this->vy_obs);
    … …  
    191193        IoModelFetchData(&this->waitonlock,iomodel_handle,"waitonlock");
    192194        IoModelFetchData(&this->kff,iomodel_handle,"kff");
     195        IoModelFetchData(&this->grounding_line_migration,iomodel_handle,"grounding_line_migration");
    193196
    194197        /*!Get thermal parameters: */
    … …  
    214217                IoModelFetchData(&this->qmu_save_femmodel,iomodel_handle,"qmu_save_femmodel");
    215218        }
     219
     220
    216221        /*i/o: */
    217222        IoModelFetchData(&this->io_gather,iomodel_handle,"io_gather");
    … …  
    242247        this->z=NULL;
    243248        this->elements=NULL;
     249        this->elementconnectivity=NULL;
    244250        this->elements_type=NULL;
    245251        this->vertices_type=NULL;
    … …  
    275281        this->surface=NULL;
    276282        this->bed=NULL;
     283        this->bathymetry=NULL;
    277284        this->elementoniceshelf=NULL;
    278285        this->elementonwater=NULL;

  • issm/trunk/src/c/objects/IoModel.h

    r6946 r7089  
    3232                double* elements_type;
    3333                double* vertices_type;
     34                double* elementconnectivity;
    3435
    3536                /*3d mesh: */
    … …  
    8384                double* surface;
    8485                double* bed;
     86                double* bathymetry;
    8587                double* elementoniceshelf;
    8688                double* elementonwater;
    … …  
    189191                int      qmu_save_femmodel;
    190192
     193                /*grounding line migration: */
     194                int      grounding_line_migration;
     195
    191196                /*exterior partitioning data, to be carried around: */
    192197                bool*   my_elements;

  • issm/trunk/src/c/objects/Loads/Riftfront.cpp

    r6748 r7089  
    8787        this->length=*(iomodel->riftinfo+RIFTINFOSIZE*i+6);
    8888        this->fraction=*(iomodel->riftinfo+RIFTINFOSIZE*i+9);
    89         this->state=*(iomodel->riftinfo+RIFTINFOSIZE*i+11);
     89        this->state=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+11);
    9090
    9191        //intialize inputs, and add as many inputs per element as requested:

  • issm/trunk/src/c/objects/Node.h

    r6412 r7089  
    2121class Node: public Object ,public Update{
    2222
    23         private:
     23        public:
    2424
    2525                int         id;  //unique arbitrary id.
    … …  
    3131                int            analysis_type;
    3232               
    33         public:
    3433
    3534                /*Node constructors, destructors {{{1*/

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

    r6200 r7089  
    5757#include "./ElementResults/DoubleElementResult.h"
    5858#include "./ElementResults/TriaVertexElementResult.h"
    59 #include "./ElementResults/PentaVertexElementResult.h"
     59#include "./ElementResults/PentaVertexElementResult.h"
     60#include "./ElementResults/BoolElementResult.h"
    6061
    6162/*ExternalResults: */

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

    r6962 r7089  
    2222        bool   time_adapt=false;
    2323        int    output_frequency;
     24        bool   grounding_line_migration;
    2425
    2526        /*intermediary: */
    … …  
    3536        femmodel->parameters->FindParam(&output_frequency,OutputFrequencyEnum);
    3637        femmodel->parameters->FindParam(&time_adapt,TimeAdaptEnum);
     38        femmodel->parameters->FindParam(&grounding_line_migration,GroundingLineMigrationEnum);
    3739
    3840        /*initialize: */
    … …  
    5860                _printf_(VerboseSolution(),"%s\n","   computing new thickness");
    5961                prognostic_core(femmodel);
    60        
     62
     63                if (grounding_line_migration){
     64                        _printf_(VerboseSolution(),"%s\n","   computing new grounding line position");
     65                        GroundingLineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     66                }
     67
    6168                if(solution_type==Transient2DSolutionEnum && !control_analysis && (step%output_frequency==0 || time==finaltime)){
    6269                        _printf_(VerboseSolution(),"%s\n","   saving results\n");
    … …  
    6875                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time);
    6976                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time);
     77                        if(grounding_line_migration)InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ElementOnIceShelfEnum,step,time);
    7078
    7179                        /*unload results*/
Note: See TracChangeset for help on using the changeset viewer.