Index: /issm/trunk/src/m/classes/public/BasinConstrain.m
===================================================================
--- /issm/trunk/src/m/classes/public/BasinConstrain.m	(revision 1234)
+++ /issm/trunk/src/m/classes/public/BasinConstrain.m	(revision 1234)
@@ -0,0 +1,33 @@
+function md=BasinConstrain(md,domain);
+%BASINCONSTRAIN constrain basin, using a constraint domain outline, to dirichlet boundary conditions.
+%
+% usage: md=BasinConstrain(md,constraindomain)
+%          
+%            where constraindomain is an Argus domain outline file enclosing the geographical area of interest.
+%
+% See also: N/A
+
+%prepare barycentre positions.
+xelem=md.x(md.elements)*[1;1;1]/3;
+yelem=md.y(md.elements)*[1;1;1]/3;
+
+%now, flag grids and elements outside the domain outline.
+[gridondomain elementondomain]=ContourToMesh(md.elements,md.x,md.y,expread(domain,1),'element and node',2);
+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);
+grids=unique(md.elements(pos,:));
+
+md.gridondirichlet_diag(grids)=1;
+md.dirichletvalues_diag(grids,1)=md.vx_obs(grids);
+md.dirichletvalues_diag(grids,2)=md.vy_obs(grids);
+