Index: issm/trunk/src/m/parameterization/setflowequation.py
===================================================================
--- issm/trunk/src/m/parameterization/setflowequation.py	(revision 14310)
+++ issm/trunk/src/m/parameterization/setflowequation.py	(revision 16137)
@@ -11,9 +11,9 @@
 
 	   This routine works like plotmodel: it works with an even number of inputs
-	   'hutter','macayeal','pattyn','l1l2','stokes' and 'fill' are the possible options
+	   'SIA','SSA','HO','L1L2','FS' and 'fill' are the possible options
 	   that must be followed by the corresponding exp file or flags list
 	   It can either be a domain file (argus type, .exp extension), or an array of element flags. 
 	   If user wants every element outside the domain to be 
-	   setflowequationd, add '~' to the name of the domain file (ex: '~Pattyn.exp');
+	   setflowequationd, add '~' to the name of the domain file (ex: '~HO.exp');
 	   an empty string '' will be considered as an empty domain
 	   a string 'all' will be considered as the entire domain
@@ -24,6 +24,5 @@
 
 	   Example:
-	      md=setflowequation(md,'pattyn','Pattyn.exp','macayeal',md.mask.elementonfloatingice,'fill','hutter');
-	      md=setflowequation(md,'pattyn','Pattyn.exp',fill','hutter','coupling','tiling');
+	      md=setflowequation(md,'HO','HO.exp',fill','SIA','coupling','tiling');
 	"""
 
@@ -42,87 +41,87 @@
 
 	#recover elements distribution
-	hutterflag   = FlagElements(md,options.getfieldvalue('hutter',''))
-	macayealflag = FlagElements(md,options.getfieldvalue('macayeal',''))
-	pattynflag   = FlagElements(md,options.getfieldvalue('pattyn',''))
-	l1l2flag     = FlagElements(md,options.getfieldvalue('l1l2',''))
-	stokesflag   = FlagElements(md,options.getfieldvalue('stokes',''))
+	SIAflag   = FlagElements(md,options.getfieldvalue('SIA',''))
+	SSAflag = FlagElements(md,options.getfieldvalue('SSA',''))
+	HOflag   = FlagElements(md,options.getfieldvalue('HO',''))
+	L1L2flag     = FlagElements(md,options.getfieldvalue('L1L2',''))
+	FSflag   = FlagElements(md,options.getfieldvalue('FS',''))
 	filltype     = options.getfieldvalue('fill','none')
 
 	#Flag the elements that have not been flagged as filltype
-	if   strcmpi(filltype,'hutter'):
-		hutterflag[numpy.nonzero(numpy.logical_not(logical_or_n(macayealflag,pattynflag)))]=True
-	elif strcmpi(filltype,'macayeal'):
-		macayealflag[numpy.nonzero(numpy.logical_not(logical_or_n(hutterflag,pattynflag,stokesflag)))]=True
-	elif strcmpi(filltype,'pattyn'):
-		pattynflag[numpy.nonzero(numpy.logical_not(logical_or_n(hutterflag,macayealflag,stokesflag)))]=True
+	if   strcmpi(filltype,'SIA'):
+		SIAflag[numpy.nonzero(numpy.logical_not(logical_or_n(SSAflag,HOflag)))]=True
+	elif strcmpi(filltype,'SSA'):
+		SSAflag[numpy.nonzero(numpy.logical_not(logical_or_n(SIAflag,HOflag,FSflag)))]=True
+	elif strcmpi(filltype,'HO'):
+		HOflag[numpy.nonzero(numpy.logical_not(logical_or_n(SIAflag,SSAflag,FSflag)))]=True
 
 	#check that each element has at least one flag
-	if not any(hutterflag+macayealflag+l1l2flag+pattynflag+stokesflag):
-		raise TypeError("elements type not assigned, must be specified")
+	if not any(SIAflag+SSAflag+L1L2flag+HOflag+FSflag):
+		raise TypeError("elements type not assigned, supported models are 'SIA','SSA','HO' and 'FS'")
 
 	#check that each element has only one flag
-	if any(hutterflag+macayealflag+l1l2flag+pattynflag+stokesflag>1):
+	if any(SIAflag+SSAflag+L1L2flag+HOflag+FSflag>1):
 		print "setflowequation warning message: some elements have several types, higher order type is used for them"
-		hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]=False
-		hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]=False
-		macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]=False
-
-	#Check that no pattyn or stokes for 2d mesh
+		SIAflag[numpy.nonzero(numpy.logical_and(SIAflag,SSAflag))]=False
+		SIAflag[numpy.nonzero(numpy.logical_and(SIAflag,HOflag))]=False
+		SSAflag[numpy.nonzero(numpy.logical_and(SSAflag,HOflag))]=False
+
+	#Check that no HO or FS for 2d mesh
 	if md.mesh.dimension==2:
-		if numpy.any(logical_or_n(l1l2flag,stokesflag,pattynflag)):
-			raise TypeError("stokes and pattyn elements not allowed in 2d mesh, extrude it first")
-
-	#Stokes can only be used alone for now:
-	if any(stokesflag) and any(hutterflag):
-		raise TypeError("stokes cannot be used with any other model for now, put stokes everywhere")
+		if numpy.any(logical_or_n(L1L2flag,FSflag,HOflag)):
+			raise TypeError("FS and HO elements not allowed in 2d mesh, extrude it first")
+
+	#FS can only be used alone for now:
+	if any(FSflag) and any(SIAflag):
+		raise TypeError("FS cannot be used with any other model for now, put FS everywhere")
 
 	#Initialize node fields
-	nodeonhutter=numpy.zeros(md.mesh.numberofvertices,bool)
-	nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:]-1]=True
-	nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices,bool)
-	nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
-	nodeonl1l2=numpy.zeros(md.mesh.numberofvertices,bool)
-	nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:]-1]=True
-	nodeonpattyn=numpy.zeros(md.mesh.numberofvertices,bool)
-	nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
-	nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonSIA=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonSIA[md.mesh.elements[numpy.nonzero(SIAflag),:]-1]=True
+	nodeonSSA=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True
+	nodeonL1L2=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonL1L2[md.mesh.elements[numpy.nonzero(L1L2flag),:]-1]=True
+	nodeonHO=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonHO[md.mesh.elements[numpy.nonzero(HOflag),:]-1]=True
+	nodeonFS=numpy.zeros(md.mesh.numberofvertices,bool)
 	noneflag=numpy.zeros(md.mesh.numberofelements,bool)
 
-	#First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal)
-	if any(stokesflag):
-#		fullspcnodes=double((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy)+~isnan(md.diagnostic.spcvz))==3 | (nodeonpattyn & nodeonstokes));         %find all the nodes on the boundary of the domain without icefront
-		fullspcnodes=numpy.logical_or(numpy.logical_not(numpy.isnan(md.diagnostic.spcvx)).astype(int)+ \
-		                              numpy.logical_not(numpy.isnan(md.diagnostic.spcvy)).astype(int)+ \
-		                              numpy.logical_not(numpy.isnan(md.diagnostic.spcvz)).astype(int)==3, \
-		                              numpy.logical_and(nodeonpattyn,nodeonstokes).reshape(-1,1)).astype(int)    #find all the nodes on the boundary of the domain without icefront
+	#First modify FSflag to get rid of elements contrained everywhere (spc + border with HO or SSA)
+	if any(FSflag):
+#		fullspcnodes=double((~isnan(md.stressbalance.spcvx)+~isnan(md.stressbalance.spcvy)+~isnan(md.stressbalance.spcvz))==3 | (nodeonHO & nodeonFS));         %find all the nodes on the boundary of the domain without icefront
+		fullspcnodes=numpy.logical_or(numpy.logical_not(numpy.isnan(md.stressbalance.spcvx)).astype(int)+ \
+		                              numpy.logical_not(numpy.isnan(md.stressbalance.spcvy)).astype(int)+ \
+		                              numpy.logical_not(numpy.isnan(md.stressbalance.spcvz)).astype(int)==3, \
+		                              numpy.logical_and(nodeonHO,nodeonFS).reshape(-1,1)).astype(int)    #find all the nodes on the boundary of the domain without icefront
 #		fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6);         %find all the nodes on the boundary of the domain without icefront
 		fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements-1],axis=1)==6).astype(int)    #find all the nodes on the boundary of the domain without icefront
-		stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=False
-		nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
-
-	#Then complete with NoneApproximation or the other model used if there is no stokes
-	if any(stokesflag): 
-		if   any(pattynflag):    #fill with pattyn
-			pattynflag[numpy.logical_not(stokesflag)]=True
-			nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
-		elif any(macayealflag):    #fill with macayeal
-			macayealflag[numpy.logical_not(stokesflag)]=True
-			nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
+		FSflag[numpy.nonzero(fullspcelems.reshape(-1))]=False
+		nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True
+
+	#Then complete with NoneApproximation or the other model used if there is no FS
+	if any(FSflag): 
+		if   any(HOflag):    #fill with HO
+			HOflag[numpy.logical_not(FSflag)]=True
+			nodeonHO[md.mesh.elements[numpy.nonzero(HOflag),:]-1]=True
+		elif any(SSAflag):    #fill with SSA
+			SSAflag[numpy.logical_not(FSflag)]=True
+			nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True
 		else:    #fill with none 
-			noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=True
-
-	#Now take care of the coupling between MacAyeal and Pattyn
-	md.diagnostic.vertex_pairing=numpy.array([])
-	nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices,bool)
-	nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices,bool)
-	nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices,bool)
-	macayealpattynflag=numpy.zeros(md.mesh.numberofelements,bool)
-	macayealstokesflag=numpy.zeros(md.mesh.numberofelements,bool)
-	pattynstokesflag=numpy.zeros(md.mesh.numberofelements,bool)
+			noneflag[numpy.nonzero(numpy.logical_not(FSflag))]=True
+
+	#Now take care of the coupling between SSA and HO
+	md.stressbalance.vertex_pairing=numpy.array([])
+	nodeonSSAHO=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonHOFS=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonSSAFS=numpy.zeros(md.mesh.numberofvertices,bool)
+	SSAHOflag=numpy.zeros(md.mesh.numberofelements,bool)
+	SSAFSflag=numpy.zeros(md.mesh.numberofelements,bool)
+	HOFSflag=numpy.zeros(md.mesh.numberofelements,bool)
 	if   strcmpi(coupling_method,'penalties'):
-		#Create the border nodes between Pattyn and MacAyeal and extrude them
+		#Create the border nodes between HO and SSA and extrude them
 		numnodes2d=md.mesh.numberofvertices2d
 		numlayers=md.mesh.numberoflayers
-		bordernodes2d=numpy.nonzero(numpy.logical_and(nodeonpattyn[0:numnodes2d],nodeonmacayeal[0:numnodes2d]))[0]+1    #Nodes connected to two different types of elements
+		bordernodes2d=numpy.nonzero(numpy.logical_and(nodeonHO[0:numnodes2d],nodeonSSA[0:numnodes2d]))[0]+1    #Nodes connected to two different types of elements
 
 		#initialize and fill in penalties structure
@@ -131,160 +130,161 @@
 			for	i in xrange(1,numlayers):
 				penalties=numpy.vstack((penalties,numpy.hstack((bordernodes2d.reshape(-1,1),bordernodes2d.reshape(-1,1)+md.mesh.numberofvertices2d*(i)))))
-			md.diagnostic.vertex_pairing=penalties
+			md.stressbalance.vertex_pairing=penalties
 
 	elif strcmpi(coupling_method,'tiling'):
-		if   any(macayealflag) and any(pattynflag):    #coupling macayeal pattyn
+		if   any(SSAflag) and any(HOflag):    #coupling SSA HO
 			#Find node at the border
-			nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=True
-			#Macayeal elements in contact with this layer become MacAyealPattyn elements
-			matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealpattyn)[0])
+			nodeonSSAHO[numpy.nonzero(numpy.logical_and(nodeonSSA,nodeonHO))]=True
+			#SSA elements in contact with this layer become SSAHO elements
+			matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonSSAHO)[0])
 			commonelements=numpy.sum(matrixelements,axis=1)!=0
-			commonelements[numpy.nonzero(pattynflag)]=False    #only one layer: the elements previously in macayeal
-			macayealflag[numpy.nonzero(commonelements)]=False    #these elements are now macayealpattynelements
-			macayealpattynflag[numpy.nonzero(commonelements)]=True
-			nodeonmacayeal[:]=False
-			nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
+			commonelements[numpy.nonzero(HOflag)]=False    #only one layer: the elements previously in SSA
+			SSAflag[numpy.nonzero(commonelements)]=False    #these elements are now SSAHOelements
+			SSAHOflag[numpy.nonzero(commonelements)]=True
+			nodeonSSA[:]=False
+			nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True
 
 			#rule out elements that don't touch the 2 boundaries
-			pos=numpy.nonzero(macayealpattynflag)[0]
+			pos=numpy.nonzero(SSAHOflag)[0]
 			elist=numpy.zeros(numpy.size(pos),dtype=int)
-			elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
-			elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1]  ,axis=1).astype(bool)
+			elist = elist + numpy.sum(nodeonSSA[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
+			elist = elist - numpy.sum(nodeonHO[md.mesh.elements[pos,:]-1]  ,axis=1).astype(bool)
 			pos1=numpy.nonzero(elist==1)[0]
-			macayealflag[pos[pos1]]=True
-			macayealpattynflag[pos[pos1]]=False
+			SSAflag[pos[pos1]]=True
+			SSAHOflag[pos[pos1]]=False
 			pos2=numpy.nonzero(elist==-1)[0]
-			pattynflag[pos[pos2]]=True
-			macayealpattynflag[pos[pos2]]=False
+			HOflag[pos[pos2]]=True
+			SSAHOflag[pos[pos2]]=False
 
 			#Recompute nodes associated to these elements
-			nodeonmacayeal[:]=False
-			nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
-			nodeonpattyn[:]=False
-			nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
-			nodeonmacayealpattyn[:]=False
-			nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:]-1]=True
-
-		elif any(pattynflag) and any(stokesflag):    #coupling pattyn stokes
+			nodeonSSA[:]=False
+			nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True
+			nodeonHO[:]=False
+			nodeonHO[md.mesh.elements[numpy.nonzero(HOflag),:]-1]=True
+			nodeonSSAHO[:]=False
+			nodeonSSAHO[md.mesh.elements[numpy.nonzero(SSAHOflag),:]-1]=True
+
+		elif any(HOflag) and any(FSflag):    #coupling HO FS
 			#Find node at the border
-			nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=True
-			#Stokes elements in contact with this layer become PattynStokes elements
-			matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonpattynstokes)[0])
+			nodeonHOFS[numpy.nonzero(numpy.logical_and(nodeonHO,nodeonFS))]=True
+			#FS elements in contact with this layer become HOFS elements
+			matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonHOFS)[0])
 			commonelements=numpy.sum(matrixelements,axis=1)!=0
-			commonelements[numpy.nonzero(pattynflag)]=False    #only one layer: the elements previously in macayeal
-			stokesflag[numpy.nonzero(commonelements)]=False    #these elements are now macayealpattynelements
-			pattynstokesflag[numpy.nonzero(commonelements)]=True
-			nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool)
-			nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
+			commonelements[numpy.nonzero(HOflag)]=False    #only one layer: the elements previously in SSA
+			FSflag[numpy.nonzero(commonelements)]=False    #these elements are now SSAHOelements
+			HOFSflag[numpy.nonzero(commonelements)]=True
+			nodeonFS=numpy.zeros(md.mesh.numberofvertices,bool)
+			nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True
 
 			#rule out elements that don't touch the 2 boundaries
-			pos=numpy.nonzero(pattynstokesflag)[0]
+			pos=numpy.nonzero(HOFSflag)[0]
 			elist=numpy.zeros(numpy.size(pos),dtype=int)
-			elist = elist + numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
-			elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
+			elist = elist + numpy.sum(nodeonFS[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
+			elist = elist - numpy.sum(nodeonHO[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
 			pos1=numpy.nonzero(elist==1)[0]
-			stokesflag[pos[pos1]]=True
-			pattynstokesflag[pos[pos1]]=False
+			FSflag[pos[pos1]]=True
+			HOFSflag[pos[pos1]]=False
 			pos2=numpy.nonzero(elist==-1)[0]
-			pattynflag[pos[pos2]]=True
-			pattynstokesflag[pos[pos2]]=False
+			HOflag[pos[pos2]]=True
+			HOFSflag[pos[pos2]]=False
 
 			#Recompute nodes associated to these elements
-			nodeonstokes[:]=False
-			nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
-			nodeonpattyn[:]=False
-			nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
-			nodeonpattynstokes[:]=False
-			nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:]-1]=True
-
-		elif any(stokesflag) and any(macayealflag):
+			nodeonFS[:]=False
+			nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True
+			nodeonHO[:]=False
+			nodeonHO[md.mesh.elements[numpy.nonzero(HOflag),:]-1]=True
+			nodeonHOFS[:]=False
+			nodeonHOFS[md.mesh.elements[numpy.nonzero(HOFSflag),:]-1]=True
+
+		elif any(FSflag) and any(SSAflag):
 			#Find node at the border
-			nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=True
-			#Stokes elements in contact with this layer become MacAyealStokes elements
-			matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealstokes)[0])
+			nodeonSSAFS[numpy.nonzero(numpy.logical_and(nodeonSSA,nodeonFS))]=True
+			#FS elements in contact with this layer become SSAFS elements
+			matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonSSAFS)[0])
 			commonelements=numpy.sum(matrixelements,axis=1)!=0
-			commonelements[numpy.nonzero(macayealflag)]=False    #only one layer: the elements previously in macayeal
-			stokesflag[numpy.nonzero(commonelements)]=False    #these elements are now macayealmacayealelements
-			macayealstokesflag[numpy.nonzero(commonelements)]=True
-			nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool)
-			nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
+			commonelements[numpy.nonzero(SSAflag)]=False    #only one layer: the elements previously in SSA
+			FSflag[numpy.nonzero(commonelements)]=False    #these elements are now SSASSAelements
+			SSAFSflag[numpy.nonzero(commonelements)]=True
+			nodeonFS=numpy.zeros(md.mesh.numberofvertices,bool)
+			nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True
 
 			#rule out elements that don't touch the 2 boundaries
-			pos=numpy.nonzero(macayealstokesflag)[0]
+			pos=numpy.nonzero(SSAFSflag)[0]
 			elist=numpy.zeros(numpy.size(pos),dtype=int)
-			elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
-			elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1]  ,axis=1).astype(bool)
+			elist = elist + numpy.sum(nodeonSSA[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
+			elist = elist - numpy.sum(nodeonFS[md.mesh.elements[pos,:]-1]  ,axis=1).astype(bool)
 			pos1=numpy.nonzero(elist==1)[0]
-			macayealflag[pos[pos1]]=True
-			macayealstokesflag[pos[pos1]]=False
+			SSAflag[pos[pos1]]=True
+			SSAFSflag[pos[pos1]]=False
 			pos2=numpy.nonzero(elist==-1)[0]
-			stokesflag[pos[pos2]]=True
-			macayealstokesflag[pos[pos2]]=False
+			FSflag[pos[pos2]]=True
+			SSAFSflag[pos[pos2]]=False
 
 			#Recompute nodes associated to these elements
-			nodeonmacayeal[:]=False
-			nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
-			nodeonstokes[:]=False
-			nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
-			nodeonmacayealstokes[:]=False
-			nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:]-1]=True
-
-		elif any(stokesflag) and any(hutterflag):
+			nodeonSSA[:]=False
+			nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True
+			nodeonFS[:]=False
+			nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True
+			nodeonSSAFS[:]=False
+			nodeonSSAFS[md.mesh.elements[numpy.nonzero(SSAFSflag),:]-1]=True
+
+		elif any(FSflag) and any(SIAflag):
 			raise TypeError("type of coupling not supported yet")
 
-	#Create MacAyealPattynApproximation where needed
+	#Create SSAHOApproximation where needed
 	md.flowequation.element_equation=numpy.zeros(md.mesh.numberofelements,int)
 	md.flowequation.element_equation[numpy.nonzero(noneflag)]=0
-	md.flowequation.element_equation[numpy.nonzero(hutterflag)]=1
-	md.flowequation.element_equation[numpy.nonzero(macayealflag)]=2
-	md.flowequation.element_equation[numpy.nonzero(l1l2flag)]=8
-	md.flowequation.element_equation[numpy.nonzero(pattynflag)]=3
-	md.flowequation.element_equation[numpy.nonzero(stokesflag)]=4
-	md.flowequation.element_equation[numpy.nonzero(macayealpattynflag)]=5
-	md.flowequation.element_equation[numpy.nonzero(macayealstokesflag)]=6
-	md.flowequation.element_equation[numpy.nonzero(pattynstokesflag)]=7
+	md.flowequation.element_equation[numpy.nonzero(SIAflag)]=1
+	md.flowequation.element_equation[numpy.nonzero(SSAflag)]=2
+	md.flowequation.element_equation[numpy.nonzero(L1L2flag)]=8
+	md.flowequation.element_equation[numpy.nonzero(HOflag)]=3
+	md.flowequation.element_equation[numpy.nonzero(FSflag)]=4
+	md.flowequation.element_equation[numpy.nonzero(SSAHOflag)]=5
+	md.flowequation.element_equation[numpy.nonzero(SSAFSflag)]=6
+	md.flowequation.element_equation[numpy.nonzero(HOFSflag)]=7
 
 	#border
-	md.flowequation.borderpattyn=nodeonpattyn
-	md.flowequation.bordermacayeal=nodeonmacayeal
-	md.flowequation.borderstokes=nodeonstokes
+	md.flowequation.borderHO=nodeonHO
+	md.flowequation.borderSSA=nodeonSSA
+	md.flowequation.borderFS=nodeonFS
 
 	#Create vertices_type
 	md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices,int)
-	pos=numpy.nonzero(nodeonmacayeal)
+	pos=numpy.nonzero(nodeonSSA)
 	md.flowequation.vertex_equation[pos]=2
-	pos=numpy.nonzero(nodeonl1l2)
+	pos=numpy.nonzero(nodeonL1L2)
 	md.flowequation.vertex_equation[pos]=8
-	pos=numpy.nonzero(nodeonpattyn)
+	pos=numpy.nonzero(nodeonHO)
 	md.flowequation.vertex_equation[pos]=3
-	pos=numpy.nonzero(nodeonhutter)
+	pos=numpy.nonzero(nodeonSIA)
 	md.flowequation.vertex_equation[pos]=1
-	pos=numpy.nonzero(nodeonmacayealpattyn)
+	pos=numpy.nonzero(nodeonSSAHO)
 	md.flowequation.vertex_equation[pos]=5
-	pos=numpy.nonzero(nodeonstokes)
+	pos=numpy.nonzero(nodeonFS)
 	md.flowequation.vertex_equation[pos]=4
-	if any(stokesflag):
-		pos=numpy.nonzero(numpy.logical_not(nodeonstokes))
-		if not (any(pattynflag) or any(macayealflag)):
+	if any(FSflag):
+		pos=numpy.nonzero(numpy.logical_not(nodeonFS))
+		if not (any(HOflag) or any(SSAflag)):
 			md.flowequation.vertex_equation[pos]=0
-	pos=numpy.nonzero(nodeonpattynstokes)
+	pos=numpy.nonzero(nodeonHOFS)
 	md.flowequation.vertex_equation[pos]=7
-	pos=numpy.nonzero(nodeonmacayealstokes)
+	pos=numpy.nonzero(nodeonSSAFS)
 	md.flowequation.vertex_equation[pos]=6
 
 	#figure out solution types
-	md.flowequation.ishutter=any(md.flowequation.element_equation==1)
-	md.flowequation.ismacayealpattyn=bool(numpy.any(numpy.logical_or(md.flowequation.element_equation==2,md.flowequation.element_equation==3)))
-	md.flowequation.isl1l2=any(md.flowequation.element_equation==8)
-	md.flowequation.isstokes=any(md.flowequation.element_equation==4)
+	md.flowequation.isSIA=any(md.flowequation.element_equation==1)
+	md.flowequation.isSSA=any(md.flowequation.element_equation==2)
+	md.flowequation.isL1L2=any(md.flowequation.element_equation==8)
+	md.flowequation.isHO=any(md.flowequation.element_equation==3)
+	md.flowequation.isFS=any(md.flowequation.element_equation==4)
 
 	return md
 
 	#Check that tiling can work:
-	if any(md.flowequation.bordermacayeal) and any(md.flowequation.borderpattyn) and any(md.flowequation.borderpattyn + md.flowequation.bordermacayeal !=1):
+	if any(md.flowequation.borderSSA) and any(md.flowequation.borderHO) and any(md.flowequation.borderHO + md.flowequation.borderSSA !=1):
 		raise TypeError("error coupling domain too irregular")
-	if any(md.flowequation.bordermacayeal) and any(md.flowequation.borderstokes) and any(md.flowequation.borderstokes + md.flowequation.bordermacayeal !=1):
+	if any(md.flowequation.borderSSA) and any(md.flowequation.borderFS) and any(md.flowequation.borderFS + md.flowequation.borderSSA !=1):
 		raise TypeError("error coupling domain too irregular")
-	if any(md.flowequation.borderstokes) and any(md.flowequation.borderpattyn) and any(md.flowequation.borderpattyn + md.flowequation.borderstokes !=1):
+	if any(md.flowequation.borderFS) and any(md.flowequation.borderHO) and any(md.flowequation.borderHO + md.flowequation.borderFS !=1):
 		raise TypeError("error coupling domain too irregular")
 
