Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23034) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23035) @@ -2006,6 +2006,6 @@ } else if(migration_style==ContactEnum){ - if(this->inputs->Max(MaskGroundediceLevelsetEnum) > 0.) shelf=false; - else shelf=true; + if(this->inputs->Min(MaskGroundediceLevelsetEnum) < 0.) shelf=true; + else shelf=false; } else if(migration_style==NoneEnum || migration_style==AggressiveMigrationEnum || migration_style==SoftMigrationEnum || migration_style==GroundingOnlyEnum){ //Floating if all nodes are floating @@ -2320,5 +2320,5 @@ for(i=0;iPid()]<10){ + if(migration_style == ContactEnum){ phi[i]=phi_ungrounding[vertices[i]->Pid()]; if(phi[i]>=0.) b[i]=r[i]; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 23034) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 23035) @@ -213,5 +213,5 @@ virtual int FiniteElement(void)=0; virtual IssmDouble FloatingArea(bool scaled)=0; - virtual void FSContactMigration(Vector* vertexgrounded,Vector* vertexfloating)=0; + virtual void FSContactMigration(Vector* vertex_sigmann,Vector* vertex_waterpressure)=0; virtual Element* GetBasalElement(void)=0; virtual int GetElementType(void)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23034) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23035) @@ -741,5 +741,5 @@ } /*}}}*/ -void Penta::FSContactMigration(Vector* vertexgrounded,Vector* vertexfloating){/*{{{*/ +void Penta::FSContactMigration(Vector* vertex_sigmann,Vector* vertex_waterpressure){/*{{{*/ if(!IsOnBase()) return; @@ -747,82 +747,78 @@ int approximation; inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){ - for(int i=0;iSetValue(vertices[i]->Pid(),+9999.,INS_VAL); - vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL); + if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum || approximation==HOFSApproximationEnum){ + _error_("Cannot compute contact condition for non FS elements"); + } + + /*Intermediaries*/ + IssmDouble* xyz_list = NULL; + IssmDouble bed_normal[3],base[NUMVERTICES],bed[NUMVERTICES],surface[NUMVERTICES],phi[NUMVERTICES]; + IssmDouble water_pressure[NUMVERTICES],pressureice[NUMVERTICES],pressure[NUMVERTICES]; + IssmDouble sigmaxx[NUMVERTICES],sigmayy[NUMVERTICES],sigmazz[NUMVERTICES],sigmaxy[NUMVERTICES]; + IssmDouble sigmayz[NUMVERTICES],sigmaxz[NUMVERTICES],sigma_nn[NUMVERTICES]; + IssmDouble viscosity,epsilon[NUMVERTICES]; + GetInputListOnVertices(&base[0],BaseEnum); + GetInputListOnVertices(&bed[0],BedEnum); + GetInputListOnVertices(&surface[0],SurfaceEnum); + GetInputListOnVertices(&pressure[0],PressureEnum); + GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum); + IssmDouble rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble gravity = matpar->GetMaterialParameter(ConstantsGEnum); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will be needing: */ + Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input = inputs->GetInput(VzEnum); _assert_(vz_input); + + /*1. Recover stresses at the base*/ + GaussPenta* gauss=new GaussPenta(); + for (int iv=0;ivGaussVertex(iv); + + /*Compute strain rate viscosity and pressure: */ + this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); + this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input); + /*FIXME: this is for Hongju only*/ + //pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]); + //if (pressure[iv]/pressureice[iv]>1) pressure[iv]=pressureice[iv]; + + /*Compute Stress*/ + sigmaxx[iv]=2*viscosity*epsilon[0]-pressure[iv]; // sigma = nu eps - pressure + sigmayy[iv]=2*viscosity*epsilon[1]-pressure[iv]; + sigmazz[iv]=2*viscosity*epsilon[2]-pressure[iv]; + sigmaxy[iv]=2*viscosity*epsilon[3]; + sigmaxz[iv]=2*viscosity*epsilon[4]; + sigmayz[iv]=2*viscosity*epsilon[5]; + } + + /*2. compute contact condition*/ + for(int i=0;i=0.){ + NormalBase(&bed_normal[0],xyz_list); + 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]); + water_pressure[i]=-gravity*rho_water*base[i]; + vertex_sigmann->SetValue(vertices[i]->Pid(),sigma_nn[i],ADD_VAL); + vertex_waterpressure->SetValue(vertices[i]->Pid(),water_pressure[i],ADD_VAL); } - } - else { - /*Intermediaries*/ - IssmDouble* xyz_list = NULL; - IssmDouble pressure,water_pressure,sigma_nn,viscosity,bed,base; - IssmDouble bed_normal[3]; - IssmDouble epsilon[6]; /* epsilon=[exx eyy ezz exy exz eyz];*/ - IssmDouble surface=0,value=0; - bool grounded; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list); - - /*Retrieve all inputs we will be needing: */ - Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input); - Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input); - Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input); - Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input = inputs->GetInput(VzEnum); _assert_(vz_input); - - /*Create gauss point in the middle of the basal edge*/ - Gauss* gauss=NewGaussBase(1); - gauss->GaussPoint(0); - - if(!IsFloating()){ - /*Check for basal force only if grounded and touching GL*/ - this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); - this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input); - pressure_input->GetInputValue(&pressure, gauss); - base_input->GetInputValue(&base, gauss); _assert_(base<0.); - - /*Compute Stress*/ - IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure; - IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure; - IssmDouble sigma_zz=2.*viscosity*epsilon[2]-pressure; - IssmDouble sigma_xy=2.*viscosity*epsilon[3]; - IssmDouble sigma_xz=2.*viscosity*epsilon[4]; - IssmDouble sigma_yz=2.*viscosity*epsilon[5]; - - /*Get normal vector to the bed */ - NormalBase(&bed_normal[0],xyz_list); - - /*basalforce*/ - 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] - + 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]; - - /*Compute water pressure*/ - IssmDouble rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum); - IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); - IssmDouble gravity = matpar->GetMaterialParameter(ConstantsGEnum); - water_pressure=gravity*rho_water*base; - - /*Compare basal stress to water pressure and determine whether it should ground*/ - if (sigma_nnSetValue(vertices[i]->Pid(),+1.,ADD_VAL); + else vertex_waterpressure->SetValue(vertices[i]->Pid(),+1.,ADD_VAL); } - else{ - /*Check for basal elevation if floating*/ - base_input->GetInputValue(&base, gauss); - bed_input->GetInputValue(&bed, gauss); - if(baseSetValue(vertices[i]->Pid(),+1.,INS_VAL); - else vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL); - } - - /*clean up*/ - delete gauss; - xDelete(xyz_list); - } + } + + /*clean up*/ + delete gauss; + xDelete(xyz_list); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 23034) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 23035) @@ -66,5 +66,5 @@ int FiniteElement(void); IssmDouble FloatingArea(bool scaled); - void FSContactMigration(Vector* vertexgrounded,Vector* vertexfloating); + void FSContactMigration(Vector* vertex_sigmann,Vector* vertex_waterpressure); IssmDouble GetArea3D(void){_error_("not implemented yet!");}; IssmDouble GetAreaSpherical(void){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23034) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23035) @@ -1042,5 +1042,5 @@ } /*}}}*/ -void Tria::FSContactMigration(Vector* vertexgrounded,Vector* vertexfloating){/*{{{*/ +void Tria::FSContactMigration(Vector* vertex_sigmann,Vector* vertex_waterpressure){/*{{{*/ if(!IsOnBase()) return; @@ -1050,80 +1050,68 @@ if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){ - for(int i=0;iSetValue(vertices[i]->Pid(),+9999.,INS_VAL); - vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL); - } - } - else{ - /*Intermediaries*/ - IssmDouble* xyz_list = NULL; - IssmDouble* xyz_list_base = NULL; - IssmDouble pressure,water_pressure,sigma_nn,viscosity,bed,base; - IssmDouble bed_normal[2]; - IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - IssmDouble surface=0,value=0; - bool grounded; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list); - GetVerticesCoordinatesBase(&xyz_list_base); - - /*Retrieve all inputs we will be needing: */ - Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input); - Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input); - Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input); - Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); - - /*Create gauss point in the middle of the basal edge*/ - Gauss* gauss=NewGaussBase(1); - gauss->GaussPoint(0); - - if(!IsFloating()){ - /*Check for basal force only if grounded and touching GL*/ - // if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){ - this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - this->material->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL); - pressure_input->GetInputValue(&pressure, gauss); - base_input->GetInputValue(&base, gauss); - - /*Compute Stress*/ - IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure; - IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure; - IssmDouble sigma_xy=2.*viscosity*epsilon[2]; - - /*Get normal vector to the bed */ - NormalBase(&bed_normal[0],xyz_list_base); - - /*basalforce*/ - 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]; - - /*Compute water pressure*/ - IssmDouble rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum); - IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); - IssmDouble gravity = matpar->GetMaterialParameter(ConstantsGEnum); - water_pressure=gravity*rho_water*base; - - /*Compare basal stress to water pressure and determine whether it should ground*/ - if (sigma_nnGetInputValue(&base, gauss); - bed_input->GetInputValue(&bed, gauss); - if(baseSetValue(vertices[i]->Pid(),+1.,INS_VAL); - else vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL); - } - - /*clean up*/ - delete gauss; - xDelete(xyz_list); - xDelete(xyz_list_base); - } + _error_(" contact contiditon only works for FS elements"); + } + /*Intermediaries*/ + IssmDouble* xyz_list = NULL; + IssmDouble bed_normal[2],base[NUMVERTICES],bed[NUMVERTICES],surface[NUMVERTICES],phi[NUMVERTICES]; + IssmDouble water_pressure[NUMVERTICES],pressureice[NUMVERTICES],pressure[NUMVERTICES]; + IssmDouble sigmaxx[NUMVERTICES],sigmayy[NUMVERTICES],sigmaxy[NUMVERTICES],sigma_nn[NUMVERTICES]; + IssmDouble viscosity,epsilon[NUMVERTICES]; + GetInputListOnVertices(&base[0],BaseEnum); + GetInputListOnVertices(&bed[0],BedEnum); + GetInputListOnVertices(&surface[0],SurfaceEnum); + GetInputListOnVertices(&pressure[0],PressureEnum); + GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum); + IssmDouble rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble gravity = matpar->GetMaterialParameter(ConstantsGEnum); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list); + /*Retrieve all inputs we will be needing: */ + Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); + + /*1. Recover stresses at the base*/ + GaussTria* gauss=new GaussTria(); + for (int iv=0;ivGaussVertex(iv); + + /*Compute strain rate viscosity and pressure: */ + this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + this->material->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL); + /*FIXME: this is for Hongju only*/ + // pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]); + // if (pressure[iv]/pressureice[iv]>1) pressure[iv]=pressureice[iv]; + + /*Compute Stress*/ + sigmaxx[iv]=2*viscosity*epsilon[0]-pressure[iv]; + sigmayy[iv]=2*viscosity*epsilon[1]-pressure[iv]; + sigmaxy[iv]=2*viscosity*epsilon[2]; + } + + /*2. compute contact condition*/ + for(int i=0;i=0.){ + NormalBase(&bed_normal[0],xyz_list); + 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]); + water_pressure[i]=-gravity*rho_water*base[i]; + vertex_sigmann->SetValue(vertices[i]->Pid(),sigma_nn[i],ADD_VAL); + vertex_waterpressure->SetValue(vertices[i]->Pid(),water_pressure[i],ADD_VAL); + } + /*If was floating*/ + else{ + /*Tricky part: + * 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 + * 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*/ + if(base[i]SetValue(vertices[i]->Pid(),+1.,ADD_VAL); + vertex_waterpressure->SetValue(vertices[i]->Pid(),+1.,ADD_VAL); + } + } + + /*clean up*/ + delete gauss; + xDelete(xyz_list); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 23034) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 23035) @@ -75,5 +75,5 @@ int FiniteElement(void); IssmDouble FloatingArea(bool scaled); - void FSContactMigration(Vector* vertexgrounded,Vector* vertexfloating); + void FSContactMigration(Vector* vertex_sigmann,Vector* vertex_waterpressure); Element* GetBasalElement(void){_error_("not implemented yet");}; void GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues); Index: /issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp (revision 23034) +++ /issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp (revision 23035) @@ -59,14 +59,14 @@ IssmDouble* ContactFSLevelset(Elements* elements,Vertices* vertices){ /*{{{*/ - Vector* vertexgrounded = NULL; - Vector* vertexfloating = NULL; - IssmDouble* serial_vertexgrounded = NULL; - IssmDouble* serial_vertexfloating = NULL; + Vector* vertex_sigmann = NULL; + Vector* vertex_waterpressure = NULL; + IssmDouble* serial_vertex_sigmann = NULL; + IssmDouble* serial_vertex_waterpressure = NULL; IssmDouble* phi = NULL; /*Initialize vector with number of vertices*/ int numberofvertices = vertices->NumberOfVertices(); - vertexgrounded = new Vector(numberofvertices); - vertexfloating = new Vector(numberofvertices); + vertex_sigmann = new Vector(numberofvertices); + vertex_waterpressure = new Vector(numberofvertices); phi = xNew(numberofvertices); @@ -74,25 +74,22 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->FSContactMigration(vertexgrounded,vertexfloating); + element->FSContactMigration(vertex_sigmann,vertex_waterpressure); } + /*Assemble vector and serialize */ + vertex_sigmann->Assemble(); + vertex_waterpressure->Assemble(); + serial_vertex_sigmann=vertex_sigmann->ToMPISerial(); + serial_vertex_waterpressure=vertex_waterpressure->ToMPISerial(); - /*Assemble vector and serialize */ - vertexgrounded->Assemble(); - vertexfloating->Assemble(); - serial_vertexgrounded=vertexgrounded->ToMPISerial(); - serial_vertexfloating=vertexfloating->ToMPISerial(); for(int i=0;i10) phi[i]=9999; - else phi[i]=-9999; + if (serial_vertex_waterpressure[i] > serial_vertex_sigmann[i]) phi[i]=-1; + else phi[i]=1; } /*free ressouces and return: */ - delete vertexgrounded; - delete vertexfloating; - xDelete(serial_vertexgrounded); - xDelete(serial_vertexfloating); + delete vertex_sigmann; + delete vertex_waterpressure; + xDelete(serial_vertex_sigmann); + xDelete(serial_vertex_waterpressure); return phi;