Index: /issm/trunk-jpl/src/m/classes/flowequation.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/flowequation.m	(revision 18075)
+++ /issm/trunk-jpl/src/m/classes/flowequation.m	(revision 18076)
@@ -144,6 +144,6 @@
 					md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[1:2]);
 				elseif strcmp(domaintype(md.mesh),'2Dvertical')
-					md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[2:4]);
-					md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[2:4]);
+					md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[2,4,5]);
+					md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[2,4,5]);
 				elseif strcmp(domaintype(md.mesh),'3D'),
 					md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[0:8]);
@@ -202,10 +202,11 @@
 			pos=find(data==1); data(pos,end)=SIAApproximationEnum();
 			pos=find(data==2); data(pos,end)=SSAApproximationEnum();
-			pos=find(data==3); data(pos,end)=HOApproximationEnum();
-			pos=find(data==4); data(pos,end)=FSApproximationEnum();
-			pos=find(data==5); data(pos,end)=SSAHOApproximationEnum();
-			pos=find(data==6); data(pos,end)=SSAFSApproximationEnum();
+			pos=find(data==3); data(pos,end)=L1L2ApproximationEnum();
+			pos=find(data==4); data(pos,end)=HOApproximationEnum();
+			pos=find(data==5); data(pos,end)=FSApproximationEnum();
+
+			pos=find(data==6); data(pos,end)=SSAHOApproximationEnum();
 			pos=find(data==7); data(pos,end)=HOFSApproximationEnum();
-			pos=find(data==8); data(pos,end)=L1L2ApproximationEnum();
+			pos=find(data==8); data(pos,end)=SSAFSApproximationEnum();
 			WriteData(fid,'data',data,'enum',FlowequationVertexEquationEnum(),'format','DoubleMat','mattype',1);
 			data=obj.element_equation;
@@ -213,10 +214,11 @@
 			pos=find(data==1); data(pos,end)=SIAApproximationEnum();
 			pos=find(data==2); data(pos,end)=SSAApproximationEnum();
-			pos=find(data==3); data(pos,end)=HOApproximationEnum();
-			pos=find(data==4); data(pos,end)=FSApproximationEnum();
-			pos=find(data==5); data(pos,end)=SSAHOApproximationEnum();
-			pos=find(data==6); data(pos,end)=SSAFSApproximationEnum();
-			pos=find(data==7); data(pos,end)=HOFSApproximationEnum();
-			pos=find(data==8); data(pos,end)=L1L2ApproximationEnum();
+			pos=find(data==3); data(pos,end)=L1L2ApproximationEnum();
+			pos=find(data==4); data(pos,end)=HOApproximationEnum();
+			pos=find(data==5); data(pos,end)=FSApproximationEnum();
+
+			pos=find(data==6); data(pos,end)=SSAHOApproximationEnum();
+			pos=find(data==7); data(pos,end)=SSAFSApproximationEnum();
+			pos=find(data==8); data(pos,end)=HOFSApproximationEnum();
 			WriteData(fid,'data',data,'enum',FlowequationElementEquationEnum(),'format','DoubleMat','mattype',2);
 		end % }}}
Index: /issm/trunk-jpl/src/m/classes/flowequation.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/flowequation.py	(revision 18075)
+++ /issm/trunk-jpl/src/m/classes/flowequation.py	(revision 18076)
@@ -121,10 +121,10 @@
 		data[numpy.nonzero(data==1)]=SIAApproximationEnum()
 		data[numpy.nonzero(data==2)]=SSAApproximationEnum()
-		data[numpy.nonzero(data==3)]=HOApproximationEnum()
-		data[numpy.nonzero(data==4)]=FSApproximationEnum()
-		data[numpy.nonzero(data==5)]=SSAHOApproximationEnum()
-		data[numpy.nonzero(data==6)]=SSAFSApproximationEnum()
+		data[numpy.nonzero(data==3)]=L1L2ApproximationEnum()
+		data[numpy.nonzero(data==4)]=HOApproximationEnum()
+		data[numpy.nonzero(data==5)]=FSApproximationEnum()
+		data[numpy.nonzero(data==6)]=SSAHOApproximationEnum()
 		data[numpy.nonzero(data==7)]=HOFSApproximationEnum()
-		data[numpy.nonzero(data==8)]=L1L2ApproximationEnum()
+		data[numpy.nonzero(data==8)]=SSAFSApproximationEnum()
 		WriteData(fid,'data',data,'enum',FlowequationVertexEquationEnum(),'format','DoubleMat','mattype',1)
 		data=copy.deepcopy(self.element_equation)
@@ -132,10 +132,10 @@
 		data[numpy.nonzero(data==1)]=SIAApproximationEnum()
 		data[numpy.nonzero(data==2)]=SSAApproximationEnum()
-		data[numpy.nonzero(data==3)]=HOApproximationEnum()
-		data[numpy.nonzero(data==4)]=FSApproximationEnum()
-		data[numpy.nonzero(data==5)]=SSAHOApproximationEnum()
-		data[numpy.nonzero(data==6)]=SSAFSApproximationEnum()
-		data[numpy.nonzero(data==7)]=HOFSApproximationEnum()
-		data[numpy.nonzero(data==8)]=L1L2ApproximationEnum()
+		data[numpy.nonzero(data==3)]=L1L2ApproximationEnum()
+		data[numpy.nonzero(data==4)]=HOApproximationEnum()
+		data[numpy.nonzero(data==5)]=FSApproximationEnum()
+		data[numpy.nonzero(data==6)]=SSAHOApproximationEnum()
+		data[numpy.nonzero(data==7)]=SSAFSApproximationEnum()
+		data[numpy.nonzero(data==8)]=HOFSApproximationEnum()
 		WriteData(fid,'data',data,'enum',FlowequationElementEquationEnum(),'format','DoubleMat','mattype',2)
 	# }}}
Index: /issm/trunk-jpl/src/m/classes/inversionvalidation.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/inversionvalidation.m	(revision 18075)
+++ /issm/trunk-jpl/src/m/classes/inversionvalidation.m	(revision 18076)
@@ -53,6 +53,6 @@
 			md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0 1]);
 			md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',...
-				{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy' 'Thickness'});
-			md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:506]);
+				{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy' 'Thickness' 'BalancethicknessApparentMassbalance'});
+			md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:507]);
 			md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
 			md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
@@ -108,5 +108,10 @@
 			WriteData(fid,'object',obj,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts);
 			WriteData(fid,'object',obj,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts);
-			WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1);
+			if(numel(obj.thickness_obs)==md.mesh.numberofelements),
+				mattype=2; 
+			else
+				mattype=1;
+			end
+			WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',mattype);
 
 			%process control parameters
@@ -134,4 +139,5 @@
 			pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum();
 			pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum();
+			pos=find(obj.cost_functions==507); data(pos)=Balancethickness2MisfitEnum();
 			WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
 			WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');
Index: /issm/trunk-jpl/src/m/classes/m1qn3inversion.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/m1qn3inversion.m	(revision 18075)
+++ /issm/trunk-jpl/src/m/classes/m1qn3inversion.m	(revision 18076)
@@ -71,10 +71,10 @@
 			md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0 1]);
 			md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',...
-				{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy' 'Thickness'});
+				{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy' 'Thickness' 'BalancethicknessApparentMassbalance'});
 			md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0);
 			md = checkfield(md,'fieldname','inversion.maxiter','numel',1,'>=',0);
 			md = checkfield(md,'fieldname','inversion.dxmin','numel',1,'>',0);
 			md = checkfield(md,'fieldname','inversion.gttol','numel',1,'>',0);
-			md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:506]);
+			md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:507]);
 			md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
 			md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
@@ -138,5 +138,10 @@
 			WriteData(fid,'object',obj,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts);
 			WriteData(fid,'object',obj,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts);
-			WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1);
+			if(numel(obj.thickness_obs)==md.mesh.numberofelements),
+				mattype=2;
+			else
+				mattype=1;
+			end
+			WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',mattype);
 
 			%process control parameters
@@ -164,4 +169,5 @@
 			pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum();
 			pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum();
+			pos=find(obj.cost_functions==507); data(pos)=Balancethickness2MisfitEnum();
 			WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
 			WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');
Index: /issm/trunk-jpl/src/m/parameterization/setflowequation.m
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/setflowequation.m	(revision 18075)
+++ /issm/trunk-jpl/src/m/parameterization/setflowequation.m	(revision 18076)
@@ -70,7 +70,7 @@
 end
 
-%Check that no L1L2 or HO or FS for 2d mesh
+%Check that no HO or FS for 2d mesh
 if strcmp(domaintype(md.mesh),'2Dhorizontal')
-	if any(L1L2flag | FSflag | HOflag)
+	if any(FSflag | HOflag)
 		error('FS and HO elements not allowed in 2d mesh, extrude it first')
 	end
@@ -239,15 +239,15 @@
 end
 
-%Create MacaAyealHOApproximation where needed
+%Create element equations
 md.flowequation.element_equation=zeros(md.mesh.numberofelements,1);
 md.flowequation.element_equation(find(noneflag))=0;
 md.flowequation.element_equation(find(SIAflag))=1;
 md.flowequation.element_equation(find(SSAflag))=2;
-md.flowequation.element_equation(find(L1L2flag))=8;
-md.flowequation.element_equation(find(HOflag))=3;
-md.flowequation.element_equation(find(FSflag))=4;
-md.flowequation.element_equation(find(SSAHOflag))=5;
-md.flowequation.element_equation(find(SSAFSflag))=6;
-md.flowequation.element_equation(find(HOFSflag))=7;
+md.flowequation.element_equation(find(L1L2flag))=3;
+md.flowequation.element_equation(find(HOflag))=4;
+md.flowequation.element_equation(find(FSflag))=5;
+md.flowequation.element_equation(find(SSAHOflag))=6;
+md.flowequation.element_equation(find(SSAFSflag))=7;
+md.flowequation.element_equation(find(HOFSflag))=8;
 
 %border
@@ -258,16 +258,14 @@
 %Create vertices_type
 md.flowequation.vertex_equation=zeros(md.mesh.numberofvertices,1);
+pos=find(nodeonSIA);
+md.flowequation.vertex_equation(pos)=1;
 pos=find(nodeonSSA);
 md.flowequation.vertex_equation(pos)=2;
 pos=find(nodeonL1L2);
-md.flowequation.vertex_equation(pos)=8;
+md.flowequation.vertex_equation(pos)=3;
 pos=find(nodeonHO);
-md.flowequation.vertex_equation(pos)=3;
-pos=find(nodeonSIA);
-md.flowequation.vertex_equation(pos)=1;
-pos=find(nodeonSSAHO);
+md.flowequation.vertex_equation(pos)=4;
+pos=find(nodeonFS);
 md.flowequation.vertex_equation(pos)=5;
-pos=find(nodeonFS);
-md.flowequation.vertex_equation(pos)=4;
 if any(FSflag),
 	pos=find(~nodeonFS);
@@ -276,15 +274,17 @@
 	end
 end
+pos=find(nodeonSSAHO);
+md.flowequation.vertex_equation(pos)=6;
 pos=find(nodeonHOFS);
 md.flowequation.vertex_equation(pos)=7;
 pos=find(nodeonSSAFS);
-md.flowequation.vertex_equation(pos)=6;
+md.flowequation.vertex_equation(pos)=8;
 
 %figure out solution types
 md.flowequation.isSIA  = double(any(md.flowequation.element_equation == 1));
 md.flowequation.isSSA  = double(any(md.flowequation.element_equation == 2));
-md.flowequation.isHO   = double(any(md.flowequation.element_equation == 3));
-md.flowequation.isL1L2 = double(any(md.flowequation.element_equation == 8));
-md.flowequation.isFS   = double(any(md.flowequation.element_equation == 4));
+md.flowequation.isL1L2 = double(any(md.flowequation.element_equation == 3));
+md.flowequation.isHO   = double(any(md.flowequation.element_equation == 4));
+md.flowequation.isFS   = double(any(md.flowequation.element_equation == 5));
 
 return
Index: /issm/trunk-jpl/src/m/parameterization/setflowequation.py
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/setflowequation.py	(revision 18075)
+++ /issm/trunk-jpl/src/m/parameterization/setflowequation.py	(revision 18076)
@@ -232,10 +232,10 @@
 	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
+	md.flowequation.element_equation[numpy.nonzero(L1L2flag)]=3
+	md.flowequation.element_equation[numpy.nonzero(HOflag)]=4
+	md.flowequation.element_equation[numpy.nonzero(FSflag)]=5
+	md.flowequation.element_equation[numpy.nonzero(SSAHOflag)]=6
+	md.flowequation.element_equation[numpy.nonzero(SSAFSflag)]=7
+	md.flowequation.element_equation[numpy.nonzero(HOFSflag)]=8
 
 	#border
@@ -246,31 +246,31 @@
 	#Create vertices_type
 	md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices,int)
+	pos=numpy.nonzero(nodeonSIA)
+	md.flowequation.vertex_equation[pos]=1
 	pos=numpy.nonzero(nodeonSSA)
 	md.flowequation.vertex_equation[pos]=2
 	pos=numpy.nonzero(nodeonL1L2)
-	md.flowequation.vertex_equation[pos]=8
+	md.flowequation.vertex_equation[pos]=3
 	pos=numpy.nonzero(nodeonHO)
-	md.flowequation.vertex_equation[pos]=3
-	pos=numpy.nonzero(nodeonSIA)
-	md.flowequation.vertex_equation[pos]=1
-	pos=numpy.nonzero(nodeonSSAHO)
+	md.flowequation.vertex_equation[pos]=4
+	pos=numpy.nonzero(nodeonFS)
 	md.flowequation.vertex_equation[pos]=5
-	pos=numpy.nonzero(nodeonFS)
-	md.flowequation.vertex_equation[pos]=4
 	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(nodeonSSAHO)
+	md.flowequation.vertex_equation[pos]=6
 	pos=numpy.nonzero(nodeonHOFS)
 	md.flowequation.vertex_equation[pos]=7
 	pos=numpy.nonzero(nodeonSSAFS)
-	md.flowequation.vertex_equation[pos]=6
+	md.flowequation.vertex_equation[pos]=8
 
 	#figure out solution types
 	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)
+	md.flowequation.isL1L2=any(md.flowequation.element_equation==3)
+	md.flowequation.isHO=any(md.flowequation.element_equation==4)
+	md.flowequation.isFS=any(md.flowequation.element_equation==5)
 
 	return md
Index: /issm/trunk-jpl/src/m/plot/plot_elementstype.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_elementstype.m	(revision 18075)
+++ /issm/trunk-jpl/src/m/plot/plot_elementstype.m	(revision 18076)
@@ -18,46 +18,54 @@
 
 if is2d
+	%None elements
+	posNONE=find(data==0);
+	A=elements(posNONE,1); B=elements(posNONE,2); C=elements(posNONE,3); 
+	p1=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',0,'FaceColor','flat','EdgeColor',edgecolor);
 	%SIA elements
 	posH=find(data==1);
 	A=elements(posH,1); B=elements(posH,2); C=elements(posH,3); 
-	p1=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',1,'FaceColor','flat','EdgeColor',edgecolor);
+	p2=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',1,'FaceColor','flat','EdgeColor',edgecolor);
 	%SSA element
 	posM=find(data==2);
 	A=elements(posM,1); B=elements(posM,2); C=elements(posM,3); 
-	p2=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',2,'FaceColor','flat','EdgeColor',edgecolor);
+	p3=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',2,'FaceColor','flat','EdgeColor',edgecolor);
+	%L1L2 element
+	posM=find(data==3);
+	A=elements(posM,1); B=elements(posM,2); C=elements(posM,3); 
+	p4=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',2,'FaceColor','flat','EdgeColor',edgecolor);
 	%HO element
-	posP=find(data==3);
+	posP=find(data==4);
 	A=elements(posP,1); B=elements(posP,2); C=elements(posP,3); 
-	p3=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',3,'FaceColor','flat','EdgeColor',edgecolor);
-	%SSAHO element
-	posMP=find(data==5);
-	A=elements(posMP,1); B=elements(posMP,2); C=elements(posMP,3); 
-	p5=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',5,'FaceColor','flat','EdgeColor',edgecolor);
+	p5=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',3,'FaceColor','flat','EdgeColor',edgecolor);
 	%FS elements
-	posS=find(data==4);
+	posS=find(data==5);
 	A=elements(posS,1); B=elements(posS,2); C=elements(posS,3); 
 	p6=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',4,'FaceColor','flat','EdgeColor',edgecolor);
-	%SSAFS elements
-	posMS=find(data==6);
-	A=elements(posMS,1); B=elements(posMS,2); C=elements(posMS,3); 
-	p7=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',6,'FaceColor','flat','EdgeColor',edgecolor);
+	%SSAHO element
+	posMP=find(data==6);
+	A=elements(posMP,1); B=elements(posMP,2); C=elements(posMP,3); 
+	p7=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',5,'FaceColor','flat','EdgeColor',edgecolor);
 	%HOFS elements
 	posPS=find(data==7);
 	A=elements(posPS,1); B=elements(posPS,2); C=elements(posPS,3); 
 	p8=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',7,'FaceColor','flat','EdgeColor',edgecolor);
-	%None elements
-	posNONE=find(data==0);
-	A=elements(posNONE,1); B=elements(posNONE,2); C=elements(posNONE,3); 
-	p9=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',0,'FaceColor','flat','EdgeColor',edgecolor);
-
-	legend([p1 p2 p3 p5 p6 p7 p8 p9],...
-		'SIA''s elements','SSA''s elements','HO''s elements',...
-		'SSAHO''s elements','FS''s elements','SSAFS''s elements','HOFS''s elements','None element');
+	%SSAFS elements
+	posMS=find(data==8);
+	A=elements(posMS,1); B=elements(posMS,2); C=elements(posMS,3); 
+	p9=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',6,'FaceColor','flat','EdgeColor',edgecolor);
 
 else
+	%None elements
+	PosNONE=find(data==0);
+	A=elements(PosNONE,1); B=elements(PosNONE,2); C=elements(PosNONE,3); D=elements(PosNONE,4); E=elements(PosNONE,5); F=elements(PosNONE,6);
+	p1=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
 	%SIA elements
 	posH=find(data==1);
 	A=elements(posH,1); B=elements(posH,2); C=elements(posH,3); D=elements(posH,4); E=elements(posH,5); F=elements(posH,6);
-	p1=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 1,'FaceColor','flat','EdgeColor',edgecolor);
+	p2=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 1,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 1,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 1,'FaceColor','flat','EdgeColor',edgecolor);
@@ -67,13 +75,21 @@
 	posM=find(data==2);
 	A=elements(posM,1); B=elements(posM,2); C=elements(posM,3); D=elements(posM,4); E=elements(posM,5); F=elements(posM,6);
-	p2=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 2,'FaceColor','flat','EdgeColor',edgecolor);
+	p3=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 2,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 2,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 2,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 2,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 2,'FaceColor','flat','EdgeColor',edgecolor);
+	%L1L2 elements
+	posP=find(data==3);
+	A=elements(posP,1); B=elements(posP,2); C=elements(posP,3); D=elements(posP,4); E=elements(posP,5); F=elements(posP,6);
+	p4=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
 	%HO elements
 	posP=find(data==3);
 	A=elements(posP,1); B=elements(posP,2); C=elements(posP,3); D=elements(posP,4); E=elements(posP,5); F=elements(posP,6);
-	p3=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
+	p5=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
@@ -83,5 +99,5 @@
 	posS=find(data==4);
 	A=elements(posS,1); B=elements(posS,2); C=elements(posS,3); D=elements(posS,4); E=elements(posS,5); F=elements(posS,6);
-	p4=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 4,'FaceColor','flat','EdgeColor',edgecolor);
+	p6=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 4,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 4,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 4,'FaceColor','flat','EdgeColor',edgecolor);
@@ -91,5 +107,5 @@
 	posP=find(data==5);
 	A=elements(posP,1); B=elements(posP,2); C=elements(posP,3); D=elements(posP,4); E=elements(posP,5); F=elements(posP,6);
-	p5=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 5,'FaceColor','flat','EdgeColor',edgecolor);
+	p7=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 5,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 5,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 5,'FaceColor','flat','EdgeColor',edgecolor);
@@ -99,5 +115,5 @@
 	PosPS=find(data==7);
 	A=elements(PosPS,1); B=elements(PosPS,2); C=elements(PosPS,3); D=elements(PosPS,4); E=elements(PosPS,5); F=elements(PosPS,6);
-	p6=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 7,'FaceColor','flat','EdgeColor',edgecolor);
+	p8=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 7,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 7,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 7,'FaceColor','flat','EdgeColor',edgecolor);
@@ -107,22 +123,14 @@
 	PosMS=find(data==6);
 	A=elements(PosMS,1); B=elements(PosMS,2); C=elements(PosMS,3); D=elements(PosMS,4); E=elements(PosMS,5); F=elements(PosMS,6);
-	p7=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 6,'FaceColor','flat','EdgeColor',edgecolor);
+	p9=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 6,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 6,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 6,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 6,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 6,'FaceColor','flat','EdgeColor',edgecolor);
-	%None elements
-	PosNONE=find(data==0);
-	A=elements(PosNONE,1); B=elements(PosNONE,2); C=elements(PosNONE,3); D=elements(PosNONE,4); E=elements(PosNONE,5); F=elements(PosNONE,6);
-	p8=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
 
-	legend([p1 p2 p3 p4 p5 p6 p7 p8],...
-		'SIA''s elements','SSA''s elements','HO''s elements','FS''s elements',...
-		'SSAHO''s elements','HOFS''s elements','SSAFS''s elements','None elements');
 end
+legend([p1 p2 p3 p4 p5 p6 p7 p8 p9],...
+		'None','SIA','SSA','L1L2','HO',...
+		'SSAHO','FS','SSAFS','HOFS');
 
 %apply options