from numpy import *
import FlagElements as fe

def 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');

	#%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. {{{1
	elementonfloatingice = fe.FlagElements(md, floatingicename)
	elementongroundedice = fe.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 = logical_and(elementonfloatingice,~elementongroundedice)
	elementongroundedice = ~elementonfloatingice

	#the order here is important. we choose vertexongroundedice as default on the grounding line.
	vertexonfloatingice = zeros(md.mesh.numberofvertices,'bool')
	vertexongroundedice = zeros(md.mesh.numberofvertices,'bool')

	pos=argwhere(elementongroundedice==1)
	pos=md.mesh.elements[pos,:]-1
	if pos.size:
		vertexongroundedice[pos]=True

	pos=argwhere(~vertexongroundedice)
	if pos.size:
		vertexonfloatingice[pos]=True;
	#%}}}

	#Return: 
	md.mask.elementonfloatingice = double(elementonfloatingice)
	md.mask.vertexonfloatingice = double(vertexonfloatingice)
	md.mask.elementongroundedice = double(elementongroundedice)
	md.mask.vertexongroundedice = double(vertexongroundedice)
	md.mask.vertexonwater = zeros(md.mesh.numberofvertices)
	md.mask.elementonwater = zeros(md.mesh.numberofelements)
	return md
