Index: /issm/trunk/src/m/classes/public/BasinConstrain2.m
===================================================================
--- /issm/trunk/src/m/classes/public/BasinConstrain2.m	(revision 1317)
+++ /issm/trunk/src/m/classes/public/BasinConstrain2.m	(revision 1317)
@@ -0,0 +1,66 @@
+function md=BasinConstrain2(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,expread(domain,1),'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.gridondirichlet_diag(gridnotondomain)=1;
+md.dirichletvalues_diag(gridnotondomain,1)=md.vx_obs(gridnotondomain);
+md.dirichletvalues_diag(gridnotondomain,2)=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.gridondirichlet_diag(grids)=1;
+md.dirichletvalues_diag(grids,1)=md.vx_obs(grids);
+md.dirichletvalues_diag(grids,2)=md.vy_obs(grids);
+
+
+%make sure icefronts that are completely spc'd are taken out:
+free_segments=find(sum(md.gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2);
+md.segmentonneumann_diag=md.segmentonneumann_diag(free_segments,:);
Index: /issm/trunk/src/m/classes/public/BasinConstrainShelf.m
===================================================================
--- /issm/trunk/src/m/classes/public/BasinConstrainShelf.m	(revision 1317)
+++ /issm/trunk/src/m/classes/public/BasinConstrainShelf.m	(revision 1317)
@@ -0,0 +1,79 @@
+function md=BasinConstrainShelf(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,expread(domain,1),'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.gridondirichlet_diag(gridnotondomain)=1;
+md.dirichletvalues_diag(gridnotondomain,1)=md.vx_obs(gridnotondomain);
+md.dirichletvalues_diag(gridnotondomain,2)=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.gridondirichlet_diag(grids)=1;
+md.dirichletvalues_diag(grids,1)=md.vx_obs(grids);
+md.dirichletvalues_diag(grids,2)=md.vy_obs(grids);
+
+
+%make sure any grid with NaN velocity is spc'd:
+pos=find(isnan(md.vel_obs_raw));
+md.gridondirichlet_diag(pos)=1;
+%we spc to the smoothed value, so that control methods don't go berserk trying to figure out what reaction force to apply for the spc to stand.
+md.dirichletvalues_diag(pos,1)=md.vx_obs(pos); 
+md.dirichletvalues_diag(pos,2)=md.vy_obs(pos); 
+
+%iceshelves: any grid on icesheet is spc'd
+pos=find(md.gridonicesheet);
+md.gridondirichlet_diag(pos)=1;
+md.dirichletvalues_diag(pos,1)=md.vx_obs(pos); 
+md.dirichletvalues_diag(pos,2)=md.vy_obs(pos); 
+
+%make sure icefronts that are completely spc'd are taken out:
+free_segments=find(sum(md.gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2);
+md.segmentonneumann_diag=md.segmentonneumann_diag(free_segments,:);