source: issm/trunk-jpl/src/m/parameterization/setmask.m@ 15749

function md=setmask(md,floatingicename,groundedicename)
%SETMASK - establish boundaries between grounded and floating ice.
%
%   By default, ice is considered grounded. The contour floatingicename defines nodes 
%   for which ice is floating. The contour groundedicename defines nodes inside an floatingice, 
%   that are grounded (ie: ice rises, islands, etc ...)
%   All input files are in the Argus format (extension .exp).
%
%   Usage:
%      md=setmask(md,floatingicename,groundedicename)
%
%   Examples:
%      md=setmask(md,'all','');
%      md=setmask(md,'Iceshelves.exp','Islands.exp');

%some checks on list of arguments
if ((nargin~=3) | (nargout~=1)),
        help mask
        error('mask error message');
end

%Get assigned fields
x=md.mesh.x;
y=md.mesh.y;
elements=md.mesh.elements;

%Assign elementonfloatingice, elementongroundedice, vertexongroundedice and vertexonfloatingice. Only change at your own peril! This is synchronized heavily with the GroundingLineMigration module. {{{
elementonfloatingice=FlagElements(md,floatingicename);
elementongroundedice=FlagElements(md,groundedicename);

%Because groundedice nodes and elements can be included into an floatingice, we need to update. Remember, all the previous 
%arrays come from domain outlines that can intersect one another: 
elementonfloatingice=double((elementonfloatingice & ~elementongroundedice));
elementongroundedice=double(~elementonfloatingice);

%the order here is important. we choose vertexongroundedice as default on the grounding line.
vertexonfloatingice=zeros(md.mesh.numberofvertices,1);
vertexongroundedice=zeros(md.mesh.numberofvertices,1);
vertexongroundedice(md.mesh.elements(find(elementongroundedice),:))=1;
vertexonfloatingice(find(~vertexongroundedice))=1;
%}}}

%Return: 
md.mask.elementonfloatingice=elementonfloatingice;
md.mask.vertexonfloatingice=vertexonfloatingice;
md.mask.elementongroundedice=elementongroundedice;
md.mask.vertexongroundedice=vertexongroundedice;
md.mask.icelevelset=ones(md.mesh.numberofvertices,1);