source: issm/oecreview/Archive/22819-23185/ISSM-23034-23035.diff

Last change on this file was 23186, checked in by Mathieu Morlighem, 7 years ago

CHG: added Archive/22819-23185
File size: 20.7 KB

  • TabularUnified ../trunk-jpl/src/c/classes/Elements/Penta.h

     
    6565                void           ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    6666                int            FiniteElement(void);
    6767                IssmDouble     FloatingArea(bool scaled);
    68                 void           FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating);
     68                void           FSContactMigration(Vector<IssmDouble>* vertex_sigmann,Vector<IssmDouble>* vertex_waterpressure);
    6969                IssmDouble     GetArea3D(void){_error_("not implemented yet!");};
    7070                IssmDouble     GetAreaSpherical(void){_error_("not implemented yet!");};
    7171                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);

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

     
    212212                virtual void       ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
    213213                virtual int        FiniteElement(void)=0;
    214214                virtual IssmDouble FloatingArea(bool scaled)=0;
    215                 virtual void       FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating)=0;
     215                virtual void       FSContactMigration(Vector<IssmDouble>* vertex_sigmann,Vector<IssmDouble>* vertex_waterpressure)=0;
    216216                virtual Element*   GetBasalElement(void)=0;
    217217                virtual int        GetElementType(void)=0;
    218218                virtual IssmDouble GetHorizontalSurfaceArea(void){_error_("not implemented");};

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

     
    10411041        return floatingarea;
    10421042}
    10431043/*}}}*/
    1044 void       Tria::FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){/*{{{*/
     1044void       Tria::FSContactMigration(Vector<IssmDouble>* vertex_sigmann,Vector<IssmDouble>* vertex_waterpressure){/*{{{*/
    10451045
    10461046        if(!IsOnBase()) return;
    10471047
    … …  
    10491049        inputs->GetInputValue(&approximation,ApproximationEnum);
    10501050
    10511051        if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
    1052                 for(int i=0;i<NUMVERTICES;i++){
    1053                         vertexgrounded->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
    1054                         vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
    1055                 }
     1052                _error_(" contact contiditon only works for FS elements");
    10561053        }
    1057         else{
    1058                 /*Intermediaries*/
    1059                 IssmDouble* xyz_list = NULL;
    1060                 IssmDouble* xyz_list_base = NULL;
    1061                 IssmDouble  pressure,water_pressure,sigma_nn,viscosity,bed,base;
    1062                 IssmDouble  bed_normal[2];
    1063                 IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    1064                 IssmDouble  surface=0,value=0;
    1065                 bool grounded;
     1054        /*Intermediaries*/
     1055        IssmDouble* xyz_list = NULL;
     1056        IssmDouble  bed_normal[2],base[NUMVERTICES],bed[NUMVERTICES],surface[NUMVERTICES],phi[NUMVERTICES];
     1057        IssmDouble  water_pressure[NUMVERTICES],pressureice[NUMVERTICES],pressure[NUMVERTICES];
     1058        IssmDouble  sigmaxx[NUMVERTICES],sigmayy[NUMVERTICES],sigmaxy[NUMVERTICES],sigma_nn[NUMVERTICES];
     1059        IssmDouble  viscosity,epsilon[NUMVERTICES];
     1060        GetInputListOnVertices(&base[0],BaseEnum);
     1061        GetInputListOnVertices(&bed[0],BedEnum);
     1062        GetInputListOnVertices(&surface[0],SurfaceEnum);
     1063        GetInputListOnVertices(&pressure[0],PressureEnum);
     1064        GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
     1065        IssmDouble rho_ice   = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     1066        IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     1067        IssmDouble gravity   = matpar->GetMaterialParameter(ConstantsGEnum);
    10661068
    1067                 /* Get node coordinates and dof list: */
    1068                 GetVerticesCoordinates(&xyz_list);
    1069                 GetVerticesCoordinatesBase(&xyz_list_base);
     1069        /* Get node coordinates and dof list: */
     1070        GetVerticesCoordinates(&xyz_list);
     1071        /*Retrieve all inputs we will be needing: */
     1072        Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
     1073        Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
    10701074
    1071                 /*Retrieve all inputs we will be needing: */
    1072                 Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
    1073                 Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
    1074                 Input* bed_input      = inputs->GetInput(BedEnum);      _assert_(bed_input);
    1075                 Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
    1076                 Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
     1075        /*1. Recover stresses at the base*/
     1076        GaussTria* gauss=new GaussTria();
     1077        for (int iv=0;iv<NUMVERTICES;iv++){
     1078                gauss->GaussVertex(iv);
    10771079
    1078                 /*Create gauss point in the middle of the basal edge*/
    1079                 Gauss* gauss=NewGaussBase(1);
    1080                 gauss->GaussPoint(0);
     1080                /*Compute strain rate viscosity and pressure: */
     1081                this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     1082                this->material->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
     1083                /*FIXME: this is for Hongju only*/
     1084        //      pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]);
     1085        //      if (pressure[iv]/pressureice[iv]>1) pressure[iv]=pressureice[iv];
    10811086
    1082                 if(!IsFloating()){
    1083                         /*Check for basal force only if grounded and touching GL*/
    1084                         //              if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
    1085                         this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    1086                         this->material->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
    1087                         pressure_input->GetInputValue(&pressure, gauss);
    1088                         base_input->GetInputValue(&base, gauss);
     1087                /*Compute Stress*/
     1088                sigmaxx[iv]=2*viscosity*epsilon[0]-pressure[iv];
     1089                sigmayy[iv]=2*viscosity*epsilon[1]-pressure[iv];
     1090                sigmaxy[iv]=2*viscosity*epsilon[2];
     1091        }
    10891092
    1090                         /*Compute Stress*/
    1091                         IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
    1092                         IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
    1093                         IssmDouble sigma_xy=2.*viscosity*epsilon[2];
    1094 
    1095                         /*Get normal vector to the bed */
    1096                         NormalBase(&bed_normal[0],xyz_list_base);
    1097 
    1098                         /*basalforce*/
    1099                         sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1];
    1100 
    1101                         /*Compute water pressure*/
    1102                         IssmDouble rho_ice   = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    1103                         IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    1104                         IssmDouble gravity   = matpar->GetMaterialParameter(ConstantsGEnum);
    1105                         water_pressure=gravity*rho_water*base;
    1106 
    1107                         /*Compare basal stress to water pressure and determine whether it should ground*/
    1108                         if (sigma_nn<water_pressure) grounded=true;
    1109                         else                         grounded=false;
     1093        /*2. compute contact condition*/
     1094        for(int i=0;i<NUMVERTICES;i++){
     1095                /*If was grounded*/
     1096                if (phi[i]>=0.){
     1097                        NormalBase(&bed_normal[0],xyz_list);
     1098                        sigma_nn[i]=-1*(sigmaxx[i]*bed_normal[0]*bed_normal[0] + sigmayy[i]*bed_normal[1]*bed_normal[1]+2*sigmaxy[i]*bed_normal[0]*bed_normal[1]);
     1099                        water_pressure[i]=-gravity*rho_water*base[i];
     1100                        vertex_sigmann->SetValue(vertices[i]->Pid(),sigma_nn[i],ADD_VAL);
     1101                        vertex_waterpressure->SetValue(vertices[i]->Pid(),water_pressure[i],ADD_VAL);
    11101102                }
    1111                 else{
    1112                         /*Check for basal elevation if floating*/
    1113                         base_input->GetInputValue(&base, gauss);
    1114                         bed_input->GetInputValue(&bed, gauss);
    1115                         if(base<bed) grounded=true;
    1116                         else         grounded=false;
     1103                /*If was floating*/
     1104                else{   
     1105                        /*Tricky part:
     1106                         * 1. if base is now touching, we put 1 for sigma_nn and leave water pressure at 0 so that the rest of the module will reground this vertex
     1107                         * 2. if base is still above bed, water pressure is set as 1, sigma_nn is left as 0, so the GL module will keep it afloat*/
     1108                        if(base[i]<bed[i]) vertex_sigmann->SetValue(vertices[i]->Pid(),+1.,ADD_VAL);
     1109                        vertex_waterpressure->SetValue(vertices[i]->Pid(),+1.,ADD_VAL);
    11171110                }
    1118                 for(int i=0;i<NUMVERTICES;i++){
    1119                         if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
    1120                         else         vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
    1121                 }
     1111        }
    11221112
    1123                 /*clean up*/
    1124                 delete gauss;
    1125                 xDelete<IssmDouble>(xyz_list);
    1126                 xDelete<IssmDouble>(xyz_list_base);
    1127         }
     1113        /*clean up*/
     1114        delete gauss;
     1115        xDelete<IssmDouble>(xyz_list);
    11281116}
    11291117/*}}}*/
    11301118IssmDouble Tria::GetArea(void){/*{{{*/

  • TabularUnified ../trunk-jpl/src/c/classes/Elements/Tria.h

     
    7474                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    7575                int         FiniteElement(void);
    7676                IssmDouble  FloatingArea(bool scaled);
    77                 void        FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating);
     77                void        FSContactMigration(Vector<IssmDouble>* vertex_sigmann,Vector<IssmDouble>* vertex_waterpressure);
    7878                Element*    GetBasalElement(void){_error_("not implemented yet");};
    7979                void        GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues);
    8080                void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);

  • TabularUnified ../trunk-jpl/src/c/classes/Elements/Penta.cpp

     
    740740        return floatingarea;
    741741}
    742742/*}}}*/
    743 void       Penta::FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){/*{{{*/
     743void       Penta::FSContactMigration(Vector<IssmDouble>* vertex_sigmann,Vector<IssmDouble>* vertex_waterpressure){/*{{{*/
    744744
    745745        if(!IsOnBase()) return;
    746746
    747747        int approximation;
    748748        inputs->GetInputValue(&approximation,ApproximationEnum);
    749         if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
    750                 for(int i=0;i<NUMVERTICES;i++){
    751                         vertexgrounded->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
    752                         vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
    753                 }
     749        if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum || approximation==HOFSApproximationEnum){
     750                _error_("Cannot compute contact condition for non FS elements");
    754751        }
    755         else {
    756                 /*Intermediaries*/
    757                 IssmDouble* xyz_list = NULL;
    758                 IssmDouble  pressure,water_pressure,sigma_nn,viscosity,bed,base;
    759                 IssmDouble  bed_normal[3];
    760                 IssmDouble  epsilon[6]; /* epsilon=[exx eyy ezz exy exz eyz];*/
    761                 IssmDouble  surface=0,value=0;
    762                 bool grounded;
    763752
    764                 /* Get node coordinates and dof list: */
    765                 GetVerticesCoordinates(&xyz_list);
     753        /*Intermediaries*/
     754        IssmDouble* xyz_list = NULL;
     755        IssmDouble  bed_normal[3],base[NUMVERTICES],bed[NUMVERTICES],surface[NUMVERTICES],phi[NUMVERTICES];
     756        IssmDouble  water_pressure[NUMVERTICES],pressureice[NUMVERTICES],pressure[NUMVERTICES];
     757        IssmDouble  sigmaxx[NUMVERTICES],sigmayy[NUMVERTICES],sigmazz[NUMVERTICES],sigmaxy[NUMVERTICES];
     758        IssmDouble  sigmayz[NUMVERTICES],sigmaxz[NUMVERTICES],sigma_nn[NUMVERTICES];
     759        IssmDouble  viscosity,epsilon[NUMVERTICES];
     760        GetInputListOnVertices(&base[0],BaseEnum);
     761        GetInputListOnVertices(&bed[0],BedEnum);
     762        GetInputListOnVertices(&surface[0],SurfaceEnum);
     763        GetInputListOnVertices(&pressure[0],PressureEnum);
     764        GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
     765        IssmDouble rho_ice   = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     766        IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     767        IssmDouble gravity   = matpar->GetMaterialParameter(ConstantsGEnum);
    766768
    767                 /*Retrieve all inputs we will be needing: */
    768                 Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
    769                 Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
    770                 Input* bed_input      = inputs->GetInput(BedEnum);      _assert_(bed_input);
    771                 Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
    772                 Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
    773                 Input* vz_input       = inputs->GetInput(VzEnum);       _assert_(vz_input);
     769        /* Get node coordinates and dof list: */
     770        GetVerticesCoordinates(&xyz_list);
    774771
    775                 /*Create gauss point in the middle of the basal edge*/
    776                 Gauss* gauss=NewGaussBase(1);
    777                 gauss->GaussPoint(0);
     772        /*Retrieve all inputs we will be needing: */
     773        Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
     774        Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
     775        Input* vz_input = inputs->GetInput(VzEnum); _assert_(vz_input);
    778776
    779                 if(!IsFloating()){
    780                         /*Check for basal force only if grounded and touching GL*/
    781                         this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
    782                         this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
    783                         pressure_input->GetInputValue(&pressure, gauss);
    784                         base_input->GetInputValue(&base, gauss); _assert_(base<0.);
     777        /*1. Recover stresses at the base*/
     778        GaussPenta* gauss=new GaussPenta();
     779        for (int iv=0;iv<NUMVERTICES;iv++){
     780                gauss->GaussVertex(iv);
    785781
    786                         /*Compute Stress*/
    787                         IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
    788                         IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
    789                         IssmDouble sigma_zz=2.*viscosity*epsilon[2]-pressure;
    790                         IssmDouble sigma_xy=2.*viscosity*epsilon[3];
    791                         IssmDouble sigma_xz=2.*viscosity*epsilon[4];
    792                         IssmDouble sigma_yz=2.*viscosity*epsilon[5];
     782                /*Compute strain rate viscosity and pressure: */
     783                this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
     784                this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
     785                /*FIXME: this is for Hongju only*/
     786                //pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]);
     787                //if (pressure[iv]/pressureice[iv]>1)   pressure[iv]=pressureice[iv];
    793788
    794                         /*Get normal vector to the bed */
    795                         NormalBase(&bed_normal[0],xyz_list);
     789                /*Compute Stress*/
     790                sigmaxx[iv]=2*viscosity*epsilon[0]-pressure[iv]; // sigma = nu eps - pressure
     791                sigmayy[iv]=2*viscosity*epsilon[1]-pressure[iv];
     792                sigmazz[iv]=2*viscosity*epsilon[2]-pressure[iv];
     793                sigmaxy[iv]=2*viscosity*epsilon[3];
     794                sigmaxz[iv]=2*viscosity*epsilon[4];
     795                sigmayz[iv]=2*viscosity*epsilon[5];
     796        }
    796797
    797                         /*basalforce*/
    798                         sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + sigma_zz*bed_normal[2]*bed_normal[2]
    799                           + 2.*sigma_xy*bed_normal[0]*bed_normal[1] + 2.*sigma_xz*bed_normal[0]*bed_normal[2] + 2.*sigma_yz*bed_normal[1]*bed_normal[2];
     798        /*2. compute contact condition*/
     799        for(int i=0;i<NUMVERTICES;i++){
    800800
    801                         /*Compute water pressure*/
    802                         IssmDouble rho_ice   = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    803                         IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    804                         IssmDouble gravity   = matpar->GetMaterialParameter(ConstantsGEnum);
    805                         water_pressure=gravity*rho_water*base;
     801                /*If was grounded*/
     802                if (phi[i]>=0.){
     803                        NormalBase(&bed_normal[0],xyz_list);
     804                        sigma_nn[i]=-1*(sigmaxx[i]*bed_normal[0]*bed_normal[0] + sigmayy[i]*bed_normal[1]*bed_normal[1] + sigmazz[i]*bed_normal[2]*bed_normal[2]+2*sigmaxy[i]*bed_normal[0]*bed_normal[1]+2*sigmaxz[i]*bed_normal[0]*bed_normal[2]+2*sigmayz[i]*bed_normal[1]*bed_normal[2]);
     805                        water_pressure[i]=-gravity*rho_water*base[i];
     806                        vertex_sigmann->SetValue(vertices[i]->Pid(),sigma_nn[i],ADD_VAL);
     807                        vertex_waterpressure->SetValue(vertices[i]->Pid(),water_pressure[i],ADD_VAL);
     808                }
    806809
    807                         /*Compare basal stress to water pressure and determine whether it should ground*/
    808                         if (sigma_nn<water_pressure) grounded=true;
    809                         else                         grounded=false;
    810                 }
     810                /*If was floating*/
    811811                else{
    812                         /*Check for basal elevation if floating*/
    813                         base_input->GetInputValue(&base, gauss);
    814                         bed_input->GetInputValue(&bed, gauss);
    815                         if(base<bed) grounded=true;
    816                         else          grounded=false;
     812                        /*Tricky part:
     813                         * 1. if base is now touching, we put 1 for sigma_nn and leave water pressure at 0 so that the rest of the module will reground this vertex
     814                         * 2. if base is still above bed, water pressure is set as 1, sigma_nn is left as 0, so the GL module will keep it afloat*/
     815                        if(base[i]<bed[i]) vertex_sigmann->SetValue(vertices[i]->Pid(),+1.,ADD_VAL);
     816                        else vertex_waterpressure->SetValue(vertices[i]->Pid(),+1.,ADD_VAL);
    817817                }
    818                 for(int i=0;i<NUMVERTICES;i++){
    819                         if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
    820                         else         vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
    821                 }
     818        }
    822819
    823                 /*clean up*/
    824                 delete gauss;
    825                 xDelete<IssmDouble>(xyz_list);
    826         }
     820        /*clean up*/
     821        delete gauss;
     822        xDelete<IssmDouble>(xyz_list);
    827823}
    828824/*}}}*/
    829825void       Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){/*{{{*/

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

     
    20052005                else shelf=false;
    20062006        }
    20072007        else if(migration_style==ContactEnum){
    2008                 if(this->inputs->Max(MaskGroundediceLevelsetEnum) > 0.) shelf=false;
    2009                 else shelf=true;
     2008                if(this->inputs->Min(MaskGroundediceLevelsetEnum) < 0.) shelf=true;
     2009                else shelf=false;
    20102010        }
    20112011        else if(migration_style==NoneEnum || migration_style==AggressiveMigrationEnum || migration_style==SoftMigrationEnum || migration_style==GroundingOnlyEnum){ //Floating if all nodes are floating
    20122012                if(this->inputs->Min(MaskGroundediceLevelsetEnum) > 0.) shelf=false;
    … …  
    23192319        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
    23202320        for(i=0;i<numvertices;i++){
    23212321                /* Contact FS*/
    2322                 if(migration_style == ContactEnum && phi_ungrounding[vertices[i]->Pid()]<10){
     2322                if(migration_style == ContactEnum){
    23232323                        phi[i]=phi_ungrounding[vertices[i]->Pid()];
    23242324                        if(phi[i]>=0.) b[i]=r[i];
    23252325                }

  • TabularUnified ../trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp

     
    5858
    5959IssmDouble*    ContactFSLevelset(Elements* elements,Vertices* vertices){ /*{{{*/
    6060
    61         Vector<IssmDouble>* vertexgrounded = NULL;
    62         Vector<IssmDouble>* vertexfloating = NULL;
    63         IssmDouble*  serial_vertexgrounded = NULL;
    64         IssmDouble*  serial_vertexfloating = NULL;
     61        Vector<IssmDouble>* vertex_sigmann = NULL;
     62        Vector<IssmDouble>* vertex_waterpressure = NULL;
     63        IssmDouble*  serial_vertex_sigmann = NULL;
     64        IssmDouble*  serial_vertex_waterpressure = NULL;
    6565        IssmDouble*  phi                   = NULL;
    6666
    6767        /*Initialize vector with number of vertices*/
    6868        int numberofvertices = vertices->NumberOfVertices();
    69         vertexgrounded = new Vector<IssmDouble>(numberofvertices);
    70         vertexfloating = new Vector<IssmDouble>(numberofvertices);
     69        vertex_sigmann = new Vector<IssmDouble>(numberofvertices);
     70        vertex_waterpressure = new Vector<IssmDouble>(numberofvertices);
    7171        phi            = xNew<IssmDouble>(numberofvertices);
    7272
    7373        /*Fill vector vertices_potentially_floating: */
    7474        for(int i=0;i<elements->Size();i++){
    7575                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    76                 element->FSContactMigration(vertexgrounded,vertexfloating);
     76                element->FSContactMigration(vertex_sigmann,vertex_waterpressure);
    7777        }
     78        /*Assemble vector and serialize */
     79        vertex_sigmann->Assemble();
     80        vertex_waterpressure->Assemble();
     81        serial_vertex_sigmann=vertex_sigmann->ToMPISerial();
     82        serial_vertex_waterpressure=vertex_waterpressure->ToMPISerial();
    7883
    79         /*Assemble vector and serialize */
    80         vertexgrounded->Assemble();
    81         vertexfloating->Assemble();
    82         serial_vertexgrounded=vertexgrounded->ToMPISerial();
    83         serial_vertexfloating=vertexfloating->ToMPISerial();
    8484        for(int i=0;i<numberofvertices;i++){
    85                         if (serial_vertexgrounded[i]==1. && serial_vertexfloating[i]==1.) phi[i]=0.;
    86                         else if (serial_vertexgrounded[i]==1) phi[i]=1;
    87                         else if (serial_vertexfloating[i]==1) phi[i]=-1;
    88                         else if (serial_vertexgrounded[i]>10) phi[i]=9999;
    89                         else phi[i]=-9999;
     85                if (serial_vertex_waterpressure[i] > serial_vertex_sigmann[i]) phi[i]=-1;
     86                else phi[i]=1;
    9087        }
    9188
    9289        /*free ressouces and return: */
    93         delete vertexgrounded;
    94         delete vertexfloating;
    95         xDelete<IssmDouble>(serial_vertexgrounded);
    96         xDelete<IssmDouble>(serial_vertexfloating);
     90        delete vertex_sigmann;
     91        delete vertex_waterpressure;
     92        xDelete<IssmDouble>(serial_vertex_sigmann);
     93        xDelete<IssmDouble>(serial_vertex_waterpressure);
    9794
    9895        return phi;
    9996}
Note: See TracBrowser for help on using the repository browser.