source: issm/trunk/src/m/classes/public/BasinConstrain.m@ 5024

function md=BasinConstrain(md,domain);
%BASINCONSTRAIN - constrain basin
%
%   Constrain basin using a constraint domain outline, 
%   to dirichlet boundary conditions.
%   constraindomain is an Argus domain outline file enclosing 
%   the geographical area of interest.
%
%   Usage: 
%      md=BasinConstrain(md,constraindomain)
%
%   Example:
%      md=BasinConstrain(md,'DomainOutline.exp');
%      md=BasinConstrain(md,'~Iceshelves.exp');

%now, flag grids and elements outside the domain outline.
if ischar(domain),
        if isempty(domain),
                elementondomain=zeros(md.numberofelements,1);
                gridondomain=zeros(md.numberofgrids,1);
                invert=0;
        elseif strcmpi(domain,'all')
                elementondomain=ones(md.numberofelements,1);
                gridondomain=ones(md.numberofgrids,1);
                invert=0;
        else
                %make sure that we actually don't want the elements outside the domain outline!
                if strcmpi(domain(1),'~'),
                        domain=domain(2:end);
                        invert=1;
                else
                        invert=0;
                end
                %ok, flag elements and nodes
                [gridondomain elementondomain]=ContourToMesh(md.elements(:,1:3),md.x,md.y,domain,'element and node',2);
        end
        if invert,
                gridondomain=~gridondomain;
                elementondomain=~elementondomain;
        end
else
        error('BasinConstrain error message: domain type not supported yet');
end

%list of elements and nodes not on domain
gridnotondomain=find(~gridondomain);
elementnotondomain=find(~elementondomain);

%all elements outside the constraint domain are equivalent to water. all grids outside are spc'd.
md.spcvelocity(gridnotondomain,1:2)=1;
md.spcvelocity(gridnotondomain,4)=md.vx_obs(gridnotondomain);
md.spcvelocity(gridnotondomain,5)=md.vy_obs(gridnotondomain);
md.elementonwater(elementnotondomain)=1;

%now, make sure all elements on water have grids that are spc'd, otherwise, we'll get a singular problem.
pos=find(~md.elementonwater);
numpos=unique(md.elements(pos,:));
grids=setdiff(1:1:md.numberofgrids,numpos);
md.spcvelocity(grids,1:2)=1;
md.spcvelocity(grids,4)=md.vx_obs(grids);
md.spcvelocity(grids,5)=md.vy_obs(grids);

%make sure icefronts that are completely spc'd are taken out:
free_segments=find(sum(md.spcvelocity(md.pressureload(:,1:2),1:2),2)~=2);
md.pressureload=md.pressureload(free_segments,:);