Changeset 17330


Ignore:
Timestamp:
02/21/14 08:58:17 (11 years ago)
Author:
seroussi
Message:

BUG: fixed set FS basal boundary conditions for flow line models

Location:
issm/trunk-jpl/src/c
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r17325 r17330  
    28992899        if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base);
    29002900        if(migration_style==SubelementMigration2Enum){
     2901                if(meshtype==Mesh2DverticalEnum) _error_("Subelement Migration 2 not implemented yet for Flowline");
    29012902                gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
    29022903                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r17309 r17330  
    347347        for(int iv=0;iv<numnodes;iv++){
    348348                gauss->GaussNode(this->FiniteElement(),iv);
     349                input->GetInputValue(&pvalue[iv],gauss);
     350        }
     351        delete gauss;
     352}
     353/*}}}*/
     354void Element::GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype){/*{{{*/
     355
     356        _assert_(pvalue);
     357
     358        int    numnodes = this->GetNumberOfNodesVelocity();
     359        Input *input    = this->GetInput(enumtype);
     360        if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
     361
     362        /* Start looping on the number of vertices: */
     363        Gauss* gauss=this->NewGauss();
     364        for(int iv=0;iv<numnodes;iv++){
     365                gauss->GaussNode(this->VelocityInterpolation(),iv);
    349366                input->GetInputValue(&pvalue[iv],gauss);
    350367        }

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r17322 r17330  
    6666                void       GetInputListOnNodes(IssmDouble* pvalue,int enumtype);
    6767                void       GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
     68                void       GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype);
    6869                void       GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
    6970                void       GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r17322 r17330  
    573573
    574574        bool              mainlyfloating = true;
     575        int               meshtype,index1,index2;
    575576        const IssmPDouble epsilon        = 1.e-15;
    576         IssmDouble         phi,s1,s2,area_init,area_grounded;
    577         IssmDouble         gl[NUMVERTICES];
    578         IssmDouble         xyz_bis[3][3];
     577        IssmDouble        phi,s1,s2,area_init,area_grounded;
     578        IssmDouble        gl[NUMVERTICES];
     579        IssmDouble        xyz_bis[3][3];
    579580
    580581        /*Recover parameters and values*/
     582        parameters->FindParam(&meshtype,MeshTypeEnum);
    581583        GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
    582584
    … …  
    586588        if(gl[2]==0.) gl[2]=gl[2]+epsilon;
    587589
    588         /*Check that not all nodes are grounded or floating*/
    589         if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
    590                 phi=1;
    591         }
    592         else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
    593                 phi=0;
    594         }
    595         else{
    596                 /*Figure out if two nodes are floating or grounded*/
    597                 if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=false;
    598 
    599                 if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
    600                         /*Coordinates of point 2: same as initial point 2*/
    601                         xyz_bis[2][0]=*(xyz_list+3*2+0);
    602                         xyz_bis[2][1]=*(xyz_list+3*2+1);
    603                         xyz_bis[2][2]=*(xyz_list+3*2+2);
    604 
    605                         /*Portion of the segments*/
    606                         s1=gl[2]/(gl[2]-gl[1]);
    607                         s2=gl[2]/(gl[2]-gl[0]);
    608 
    609                         /*New point 1*/
    610                         xyz_bis[1][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*2+0));
    611                         xyz_bis[1][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*2+1));
    612                         xyz_bis[1][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*2+2));
    613 
    614                         /*New point 0*/
    615                         xyz_bis[0][0]=*(xyz_list+3*2+0)+s2*(*(xyz_list+3*0+0)-*(xyz_list+3*2+0));
    616                         xyz_bis[0][1]=*(xyz_list+3*2+1)+s2*(*(xyz_list+3*0+1)-*(xyz_list+3*2+1));
    617                         xyz_bis[0][2]=*(xyz_list+3*2+2)+s2*(*(xyz_list+3*0+2)-*(xyz_list+3*2+2));
    618                 }
    619                 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
    620                         /*Coordinates of point 0: same as initial point 2*/
    621                         xyz_bis[0][0]=*(xyz_list+3*0+0);
    622                         xyz_bis[0][1]=*(xyz_list+3*0+1);
    623                         xyz_bis[0][2]=*(xyz_list+3*0+2);
    624 
    625                         /*Portion of the segments*/
    626                         s1=gl[0]/(gl[0]-gl[1]);
    627                         s2=gl[0]/(gl[0]-gl[2]);
    628 
    629                         /*New point 1*/
    630                         xyz_bis[1][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*0+0));
    631                         xyz_bis[1][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*0+1));
    632                         xyz_bis[1][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*0+2));
    633 
    634                         /*New point 2*/
    635                         xyz_bis[2][0]=*(xyz_list+3*0+0)+s2*(*(xyz_list+3*2+0)-*(xyz_list+3*0+0));
    636                         xyz_bis[2][1]=*(xyz_list+3*0+1)+s2*(*(xyz_list+3*2+1)-*(xyz_list+3*0+1));
    637                         xyz_bis[2][2]=*(xyz_list+3*0+2)+s2*(*(xyz_list+3*2+2)-*(xyz_list+3*0+2));
    638                 }
    639                 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
    640                         /*Coordinates of point 1: same as initial point 2*/
    641                         xyz_bis[1][0]=*(xyz_list+3*1+0);
    642                         xyz_bis[1][1]=*(xyz_list+3*1+1);
    643                         xyz_bis[1][2]=*(xyz_list+3*1+2);
    644 
    645                         /*Portion of the segments*/
    646                         s1=gl[1]/(gl[1]-gl[0]);
    647                         s2=gl[1]/(gl[1]-gl[2]);
    648 
    649                         /*New point 0*/
    650                         xyz_bis[0][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*1+0));
    651                         xyz_bis[0][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*1+1));
    652                         xyz_bis[0][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*1+2));
    653 
    654                         /*New point 2*/
    655                         xyz_bis[2][0]=*(xyz_list+3*1+0)+s2*(*(xyz_list+3*2+0)-*(xyz_list+3*1+0));
    656                         xyz_bis[2][1]=*(xyz_list+3*1+1)+s2*(*(xyz_list+3*2+1)-*(xyz_list+3*1+1));
    657                         xyz_bis[2][2]=*(xyz_list+3*1+2)+s2*(*(xyz_list+3*2+2)-*(xyz_list+3*1+2));
    658                 }
    659 
    660                 /*Compute fraction of grounded element*/
    661                 GetJacobianDeterminant(&area_init, xyz_list,NULL);
    662                 GetJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL);
    663                 if(mainlyfloating==true) area_grounded=area_init-area_grounded;
    664                 phi=area_grounded/area_init;
    665         }
     590        if(meshtype==Mesh2DverticalEnum){
     591                this->EdgeOnBedIndices(&index1,&index2);
     592                if(gl[index1]>0 && gl[index2]>0) phi=1; // All grounded
     593                else if(gl[index1]<0 && gl[index2]<0) phi=0; // All floating
     594                else if(gl[index1]<0 && gl[index2]>0){ //index2 grounded
     595                        phi=1./(1.-gl[index1]/gl[index2]);
     596                }
     597                else if(gl[index2]<0 && gl[index1]>0){ //index1 grounded
     598                        phi=1./(1.-gl[index2]/gl[index1]);
     599                }
     600
     601        }
     602        else if(meshtype==Mesh3DEnum){
     603                /*Check that not all nodes are grounded or floating*/
     604                if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
     605                        phi=1;
     606                }
     607                else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
     608                        phi=0;
     609                }
     610                else{
     611                        /*Figure out if two nodes are floating or grounded*/
     612                        if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=false;
     613
     614                        if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     615                                /*Coordinates of point 2: same as initial point 2*/
     616                                xyz_bis[2][0]=*(xyz_list+3*2+0);
     617                                xyz_bis[2][1]=*(xyz_list+3*2+1);
     618                                xyz_bis[2][2]=*(xyz_list+3*2+2);
     619
     620                                /*Portion of the segments*/
     621                                s1=gl[2]/(gl[2]-gl[1]);
     622                                s2=gl[2]/(gl[2]-gl[0]);
     623
     624                                /*New point 1*/
     625                                xyz_bis[1][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*2+0));
     626                                xyz_bis[1][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*2+1));
     627                                xyz_bis[1][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*2+2));
     628
     629                                /*New point 0*/
     630                                xyz_bis[0][0]=*(xyz_list+3*2+0)+s2*(*(xyz_list+3*0+0)-*(xyz_list+3*2+0));
     631                                xyz_bis[0][1]=*(xyz_list+3*2+1)+s2*(*(xyz_list+3*0+1)-*(xyz_list+3*2+1));
     632                                xyz_bis[0][2]=*(xyz_list+3*2+2)+s2*(*(xyz_list+3*0+2)-*(xyz_list+3*2+2));
     633                        }
     634                        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
     635                                /*Coordinates of point 0: same as initial point 2*/
     636                                xyz_bis[0][0]=*(xyz_list+3*0+0);
     637                                xyz_bis[0][1]=*(xyz_list+3*0+1);
     638                                xyz_bis[0][2]=*(xyz_list+3*0+2);
     639
     640                                /*Portion of the segments*/
     641                                s1=gl[0]/(gl[0]-gl[1]);
     642                                s2=gl[0]/(gl[0]-gl[2]);
     643
     644                                /*New point 1*/
     645                                xyz_bis[1][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*0+0));
     646                                xyz_bis[1][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*0+1));
     647                                xyz_bis[1][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*0+2));
     648
     649                                /*New point 2*/
     650                                xyz_bis[2][0]=*(xyz_list+3*0+0)+s2*(*(xyz_list+3*2+0)-*(xyz_list+3*0+0));
     651                                xyz_bis[2][1]=*(xyz_list+3*0+1)+s2*(*(xyz_list+3*2+1)-*(xyz_list+3*0+1));
     652                                xyz_bis[2][2]=*(xyz_list+3*0+2)+s2*(*(xyz_list+3*2+2)-*(xyz_list+3*0+2));
     653                        }
     654                        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
     655                                /*Coordinates of point 1: same as initial point 2*/
     656                                xyz_bis[1][0]=*(xyz_list+3*1+0);
     657                                xyz_bis[1][1]=*(xyz_list+3*1+1);
     658                                xyz_bis[1][2]=*(xyz_list+3*1+2);
     659
     660                                /*Portion of the segments*/
     661                                s1=gl[1]/(gl[1]-gl[0]);
     662                                s2=gl[1]/(gl[1]-gl[2]);
     663
     664                                /*New point 0*/
     665                                xyz_bis[0][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*1+0));
     666                                xyz_bis[0][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*1+1));
     667                                xyz_bis[0][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*1+2));
     668
     669                                /*New point 2*/
     670                                xyz_bis[2][0]=*(xyz_list+3*1+0)+s2*(*(xyz_list+3*2+0)-*(xyz_list+3*1+0));
     671                                xyz_bis[2][1]=*(xyz_list+3*1+1)+s2*(*(xyz_list+3*2+1)-*(xyz_list+3*1+1));
     672                                xyz_bis[2][2]=*(xyz_list+3*1+2)+s2*(*(xyz_list+3*2+2)-*(xyz_list+3*1+2));
     673                        }
     674
     675                        /*Compute fraction of grounded element*/
     676                        GetJacobianDeterminant(&area_init, xyz_list,NULL);
     677                        GetJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL);
     678                        if(mainlyfloating==true) area_grounded=area_init-area_grounded;
     679                        phi=area_grounded/area_init;
     680                }
     681        }
     682        else _error_("mesh type "<<EnumToStringx(meshtype)<<"not supported yet ");
    666683
    667684        if(phi>1 || phi<0) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1");
    … …  
    19391956void  Tria::ResetFSBasalBoundaryCondition(void){
    19401957
    1941         int        approximation;
    1942         int        numindices;
    1943         int       *indices = NULL;
    1944         IssmDouble slope;
    1945         IssmDouble xz_plane[6];
     1958        int numnodes = this->GetNumberOfNodes();
     1959
     1960        int          approximation;
     1961        IssmDouble*  vertexonbed= NULL;
     1962        IssmDouble   slope,groundedice;
     1963        IssmDouble   xz_plane[6];
    19461964
    19471965        /*For FS only: we want the CS to be tangential to the bedrock*/
    19481966        inputs->GetInputValue(&approximation,ApproximationEnum);
    1949         if(IsFloating() || !HasEdgeOnBed() || approximation!=FSApproximationEnum) return;
    1950 
    1951         /*Get number of nodes for velocity only and base*/
    1952         int index = this->EdgeOnBedIndex();
    1953         NodeOnEdgeIndices(&numindices,&indices,index,this->VelocityInterpolation());
    1954 
     1967        if(IsFloating() || !HasNodeOnBed() ||  approximation!=FSApproximationEnum) return;
     1968
     1969        //printf("element number %i \n",this->id);
    19551970        /*Get inputs*/
    1956         Input* slope_input=inputs->GetInput(BedSlopeXEnum); _assert_(slope_input);
     1971        Input* slope_input=inputs->GetInput(BedSlopeXEnum);                             _assert_(slope_input);
     1972        Input* groundedicelevelset_input=inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(groundedicelevelset_input);
     1973        vertexonbed = xNew<IssmDouble>(numnodes);
     1974        this->GetInputListOnNodesVelocity(&vertexonbed[0],MeshVertexonbedEnum);
    19571975
    19581976        /*Loop over basal nodes and update their CS*/
    19591977        GaussTria* gauss = new GaussTria();
    1960         for(int i=0;i<numindices;i++){//FIXME
    1961 
    1962                 gauss->GaussNode(this->VelocityInterpolation(),indices[i]);
    1963                 slope_input->GetInputValue(&slope,gauss);
    1964                 IssmDouble theta = atan(slope);
    1965 
    1966                 /*New X axis                  New Z axis*/
    1967                 xz_plane[0]=cos(theta);       xz_plane[3]=0.; 
    1968                 xz_plane[1]=sin(theta);       xz_plane[4]=0.; 
    1969                 xz_plane[2]=0.;               xz_plane[5]=1.;         
    1970 
    1971                 XZvectorsToCoordinateSystem(&this->nodes[indices[i]]->coord_system[0][0],&xz_plane[0]);
     1978        for(int i=0;i<this->NumberofNodesVelocity();i++){
     1979
     1980                if(vertexonbed[i]==1){
     1981                        gauss->GaussNode(this->VelocityInterpolation(),i);
     1982                        slope_input->GetInputValue(&slope,gauss);
     1983                        groundedicelevelset_input->GetInputValue(&groundedice,gauss);
     1984                        IssmDouble theta = atan(slope);
     1985
     1986                        if(groundedice>0){
     1987                                /*New X axis                  New Z axis*/
     1988                                xz_plane[0]=cos(theta);       xz_plane[3]=0.; 
     1989                                xz_plane[1]=sin(theta);       xz_plane[4]=0.; 
     1990                                xz_plane[2]=0.;               xz_plane[5]=1.;         
     1991                                this->nodes[i]->DofInSSet(1); //vy
     1992                        }
     1993                        else{
     1994                                /*New X axis                  New Z axis*/
     1995                                xz_plane[0]=1.;               xz_plane[3]=0; 
     1996                                xz_plane[1]=0.;               xz_plane[4]=0; 
     1997                                xz_plane[2]=0;                xz_plane[5]=1.;         
     1998                                this->nodes[i]->DofInFSet(1); //vy
     1999                        }
     2000
     2001                        XZvectorsToCoordinateSystem(&this->nodes[i]->coord_system[0][0],&xz_plane[0]);
     2002                }
    19722003        }
    19732004
    19742005        /*cleanup*/
    1975         xDelete<int>(indices);
     2006        xDelete<IssmDouble>(vertexonbed);
    19762007        delete gauss;
    19772008}
Note: See TracChangeset for help on using the changeset viewer.