Changeset 15795

Timestamp:

08/12/13 13:44:11 (12 years ago)

Author:

seroussi

Message:

CHG: added spcs on the border of ice front to get coherent boundary conditions

Location:

issm/trunk-jpl/src/m/boundaryconditions

Files:

• SetIceShelfBC.m (modified) (2 diffs)
• SetMarineIceSheetBC.m (modified) (2 diffs)

issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.m

@@ -31 +31 @@
 md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
 md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
-md.stressbalance.spcvx(pos)=0;
-md.stressbalance.spcvy(pos)=0;
-md.stressbalance.spcvz(pos)=0;
 md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
 md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
+
+%Ice front position: 
+md.mask.icelevelset(find(nodeonicefront))=0;
+
+%First find segments that are not completely on the front
+if(md.mesh.dimension==2),
+        numbernodesfront=2;
+else 
+        numbernodesfront=4;
+end
+segmentsfront=md.mask.icelevelset(md.mesh.segments(:,1:numbernodesfront))==0;
+segments=find(sum(segmentsfront,2)~=numbernodesfront);
+%Find all nodes for these segments and spc them
+pos=[md.mesh.segments(segments,1:numbernodesfront)];
+md.stressbalance.spcvx(pos(:))=0;
+md.stressbalance.spcvy(pos(:))=0;
+md.stressbalance.spcvz(pos(:))=0;
 
 %Dirichlet Values
@@ -45 +59 @@
         disp('      boundary conditions for stressbalance model: spc set as zero');
 end
-
-%Ice front position
-md.mask.icelevelset(find(nodeonicefront))=0;
 
 %Create zeros basalforcings and surfaceforcings

issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m

@@ -38 +38 @@
 md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
 md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
-md.stressbalance.spcvx(pos)=0;
-md.stressbalance.spcvy(pos)=0;
-md.stressbalance.spcvz(pos)=0;
 md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
 md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
+
+
+%Position of ice front
+md.mask.icelevelset(find(vertexonicefront))=0;
+
+%First find segments that are not completely on the front
+if(md.mesh.dimension==2),
+        numbernodesfront=2;
+else 
+        numbernodesfront=4;
+end
+segmentsfront=md.mask.icelevelset(md.mesh.segments(:,1:numbernodesfront))==0;
+segments=find(sum(segmentsfront,2)~=numbernodesfront);
+%Find all nodes for these segments and spc them
+pos=[md.mesh.segments(segments,1:numbernodesfront)];
+md.stressbalance.spcvx(pos(:))=0;
+md.stressbalance.spcvy(pos(:))=0;
+md.stressbalance.spcvz(pos(:))=0;
 
 %Dirichlet Values
@@ -56 +71 @@
 pos=find(md.mesh.vertexonboundary);
 md.hydrology.spcwatercolumn(pos,1)=1;
-
-%Position of ice front
-md.mask.icelevelset(find(vertexonicefront))=0;
 
 %Create zeros basalforcings and surfaceforcings