from numpy import * 
def FlagElements(md,region):
#FLAGELEMENTS - flag the elements in an region
#
#   The region can be given with an exp file, a list of elements.
#
#   Usage: 
#      flag=FlagElements(md,region);
#
#   Example:
#      flag=FlagElements(md,'all');
#      flag=FlagElements(md,'');
#      flag=FlagElements(md,'Domain.exp');
#      flag=FlagElements(md,'~Domain.exp');
#      flag=FlagElements(md,md.mask.elementongroundedice);

	if isinstance(region,basestring):
		if not(region):
			flag=zeros(md.mesh.numberofelements,'bool')
			invert=0;
		elif region=='all':
			flag=ones(md.mesh.numberofelements,'bool')
			invert=0;
		else:
			#make sure that we actually don't want the elements outside the domain outline!
			if region[0]=='~':
				region=region[1:]
				invert=1;
			else:
				invert=0;
			
			#does the region domain outline exist or do we have to look for xlim,ylim in basinzoom?
			if not os.path.isfile(region):
				[xlim,ylim]=basinzoom('basin',region);
				flag_nodes=double(md.mesh.x<xlim(2) & md.mesh.x>xlim(1) &  md.mesh.y<ylim(2) & md.mesh.y>ylim(1));
				flag=prod(flag_nodes(md.mesh.elements),2);
			else:
				#ok, flag elements
				flag=ContourToMesh(md.mesh.elements[:,0:3],md.mesh.x,md.mesh.y,region,'element',1);
		
		if invert:
			flag=~flag;
	
	elif isinstance(region,nparray):
		if len(region)!=md.mesh.numberofelements:
			print FlagElements.__doc__
			print 'Flaglist for region must be of same size as number of elements in model'
			return []
		flag=region;
	else:
		print 'Invalid region option'
		return []

	return flag;
