Changeset 2986 for issm/trunk/src/m/classes/public/geography2.m

Timestamp:

02/08/10 14:12:12 (15 years ago)

Author:

seroussi

Message:

added possibility to use a vector to specify land (instead of an expfile) in geography2

File:

• issm/trunk/src/m/classes/public/geography2.m (modified) (3 diffs)

issm/trunk/src/m/classes/public/geography2.m

@@ -1 +1 @@
 function md=geography2(md,landname,iceshelfname,icesheetname)
+%GEOGRAPHY2 - establish land, ice sheet and ice shelf areas in a domains.
+%
+%   Usage:
+%      md=geography2(md,landname,iceshelfname,icesheetname)
+%
+%   Examples:
+%      md=geography2(md,'LandName.exp','Iceshelves.exp','Islands.exp');
 
 %Get assigned fields
@@ -7 +14 @@
 
 %recover elements and grids on land.
-[gridonland,elementonland]=ContourToMesh(elements,x,y,expread(landname,1),'element and node',2);
+if ischar(landname),
+        [gridonland,elementonland]=ContourToMesh(elements,x,y,expread(landname,1),'element and node',2);
+elseif isfloat(landname),
+        if size(landname,1)~=md.numberofelements,
+                error('Landname for area must be of same size as number of elements in model');
+        end
+        elementonland=landname;
+        gridonland=zeros(md.numberofgrids,1);
+        gridonland(md.elements(find(elementonland),:))=1;
+else
+        error('Invalid area option option');
+end
 
 %any element with 3 grids on land should be on land:
@@ -42 +60 @@
 
 %recover arrays of ice shelf grids and elements, and ice sheet grids and elements.
-if strcmp(iceshelfname,''), %no iceshelf contour file, we are dealing with a pure ice sheet.
-        gridoniceshelf=zeros(md.numberofgrids,1);
-        elementoniceshelf=zeros(md.numberofelements,1);
-elseif strcmp(iceshelfname,'all'), %we are dealing with a pure ice shelf.
-        gridoniceshelf=ones(md.numberofgrids,1);
-        elementoniceshelf=ones(md.numberofelements,1);
-else
-        [gridoniceshelf,elementoniceshelf]=ContourToMesh(elements,x,y,expread(iceshelfname,1),'element and node',2);
-end
+elementoniceshelf=FlagElements(md,iceshelfname);
+elementonicesheet=FlagElements(md,icesheetname);
 
-if strcmp(icesheetname,''), %no icesheet contour file, we are dealing with a pure ice shelf.
-        gridonicesheet=zeros(md.numberofgrids,1);
-        elementonicesheet=zeros(md.numberofelements,1);
-else
-        [gridonicesheet,elementonicesheet]=ContourToMesh(elements,x,y,expread(icesheetname,1),'element and node',2);
-end
-
+%Because icesheet grids and elements can be included into an iceshelf, we need to update. Remember, all the previous 
+%arrays come from domain outlines that can intersect one another: 
+gridoniceshelf=zeros(md.numberofgrids,1);
+gridonicesheet=zeros(md.numberofgrids,1);
+elementoniceshelf=double((elementoniceshelf & ~elementonicesheet));
+elementonicesheet=double(~elementoniceshelf);
+gridoniceshelf(md.elements(find(elementoniceshelf),:))=1;
+gridonicesheet(md.elements(find(elementonicesheet),:))=1;
 
 %now correct, so that none of the iceshelf and icesheet elements and nodes are in the water.