Index: issm/trunk/src/m/classes/@pentaelem/CreateKMatrix.m
===================================================================
--- issm/trunk/src/m/classes/@pentaelem/CreateKMatrix.m	(revision 327)
+++ issm/trunk/src/m/classes/@pentaelem/CreateKMatrix.m	(revision 328)
@@ -19,8 +19,4 @@
 	Ke=CreateKMatrixVert(pentaelem,grids,materials,inputs);
 
-elseif strcmpi(analysis_type,'diagnostic_basevert'),
-
-	Ke=CreateKMatrixBaseVert(pentaelem,grids,materials,inputs);
-
 elseif strcmpi(analysis_type,'prognostic'),
 
@@ -456,60 +452,4 @@
 		end
 	end%ends friction
-
-end %end function
-
-function Ke=CreateKMatrixBaseVert(pentaelem,grids,materials,inputs)
-
-	%We are dealing with a collapsed formulation on the lower triangle of the penta, 
-	%only on the bedrock.
-	if pentaelem.onbed~=1,
-		Ke=elemmatrix(0);
-		return;
-	end
-
-	%some variables
-	numgrids=3;
-	DOFPERGRID=1;
-	numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.
-
-	%Create elementary stiffness matrix 
-	Ke=elemmatrix(numdof);
-
-	%Get all element grid data:
-	xyz_list=getgriddata(pentaelem,grids); 
-
-	%Just keep the first 3 grids
-	xyz_list=xyz_list(1:3,:);
-
-	%Build linear indices for elementary stiffness matrix.
-	for i=1:numgrids,
-		doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
-		dof=doflist(3); %z degree of freedom
-		Ke.row_indices(i)=dof;
-	end
-
-	% Get gaussian points and weights.
-	[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
-
-	for ig=1:num_gauss2D,
-
-		%Pick up the gaussian point and its weight:
-		gauss_weight=gauss_weights(ig);
-		gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
-		
-		%Get the Jacobian determinant
-		Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
-
-		%Get L matrix if viscous basal drag present:
-		L=GetL(triaelem,gauss_coord,DOFPERGRID);
-   
-		DL_scalar=gauss_weight*Jdettria;
-
-		%Do the triple product tL*D*L. 
-		Ke_gg_vert_gaussian=L'*DL_scalar*L;
-		
-		%Add Ke_gg_drag_gaussian to Ke
-		Ke.terms=Ke.terms+Ke_gg_vert_gaussian;
-	end
 
 end %end function
Index: issm/trunk/src/m/classes/@pentaelem/CreatePVector.m
===================================================================
--- issm/trunk/src/m/classes/@pentaelem/CreatePVector.m	(revision 327)
+++ issm/trunk/src/m/classes/@pentaelem/CreatePVector.m	(revision 328)
@@ -17,8 +17,4 @@
 
 	pe=CreatePVectorVert(pentaelem,grids,materials,inputs);
-
-elseif strcmpi(analysis_type,'diagnostic_basevert'),
-
-	pe=CreatePVectorBaseVert(pentaelem,grids,materials,inputs);
 
 elseif strcmpi(analysis_type,'prognostic'),
@@ -390,147 +386,4 @@
 
 
-function pe=CreatePVectorBaseVert(pentaelem,grids,materials,inputs);
-
-	%We are dealing with a collapsed formulation on the lower triangle of the penta, 
-	%only on the bedrock.
-	if pentaelem.onbed~=1,
-		pe=elemvector(0);
-		return;
-	end
-
-	%some variables
-	numgrids=3;
-	DOFPERGRID=1;
-	numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.
-
-	%recover material parameters
-	matice=materials(pentaelem.matid).material;
-	matpar=materials(end).constants;
-	rho_ice=matpar.rho_ice;
-	rho_water=matpar.rho_water;
-	di=rho_ice/rho_water;
-
-	%Create elementary stiffness matrix 
-	pe=elemvector(numdof);
-
-	%Get all element grid data:
-	xyz_list=getgriddata(pentaelem,grids); 
-
-	%Just keep the first 3 grids
-	xyz_list=xyz_list(1:3,:);
-
-	%recover extra inputs
-	[bed_param bed_is_present]=recover_input(inputs,'bed');
-	[melting_param melting_is_present]=recover_input(inputs,'melting');
-	[accumulation_param accumulation_is_present]=recover_input(inputs,'accumulation');
-	[velocity_param velocity_is_present]=recover_input(inputs,'velocity');
-	[velocity_average_param velocity_average_is_present]=recover_input(inputs,'velocity_average');
-	[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
-	if (~velocity_is_present | ~velocity_average_is_present),
-		error('CreatePVector error message: no input velocity!');
-	end
-
-	%initialize velocity array
-	vxvy_list=zeros(numgrids,2);
-	vxvy_average_list=zeros(numgrids,2);
-	bed_list=zeros(numgrids,1);
-	melting_list=zeros(numgrids,1);
-	accumulation_list=zeros(numgrids,1);
-	thickness_list=zeros(numgrids,1);
-
-	%Build linear indices for elementary stiffness matrix.
-	for i=1:numgrids,
-		doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
-		for j=1:2,
-			dof=doflist(j);
-			vxvy_list(i,j)=velocity_param(dof);
-			vxvy_average_list(i,j)=velocity_average_param(dof);
-		end
-		if(bed_is_present) bed_list(i)=bed_param(doflist(1));end;
-		if(melting_is_present) melting_list(i)=melting_param(doflist(1));end;
-		if(accumulation_is_present) accumulation_list(i)=accumulation_param(doflist(1));end;
-		if(thickness_is_present) thickness_list(i)=thickness_param(doflist(1));end;
-		dof=doflist(3); %z degree of freedom
-		pe.row_indices(i)=dof;
-	end
-
-	if bed_is_present, bed=bed_list(1:3); else bed=pentaelem.b(1:3); end
-	if thickness_is_present, thickness=thickness_list(1:3); else thickness=pentaelem.h(1:3); end
-	if melting_is_present, melting=melting_list(1:3); else melting=pentaelem.melting(1:3); end
-	if accumulation_is_present, accumulation=accumulation_list(1:3); else accumulation=pentaelem.accumulation(1:3); end
-
-	% Get gaussian points and weights.
-	[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
-
-	%for icesheets
-	for ig=1:num_gauss2D,
-	
-		%Pick up the gaussian point and its weight:
-		gauss_weight=gauss_weights(ig);
-		gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
-		
-		%Get melting at gaussian point.
-		melting_gauss=GetParameterValue(triaelem,melting,gauss_coord);
-		
-		%Get velocity at gaussian point.
-		vx=GetParameterValue(triaelem,vxvy_list(:,1),gauss_coord);
-		vy=GetParameterValue(triaelem,vxvy_list(:,2),gauss_coord);
-
-		%Get bed slope
-		db=GetParameterDerivativeValue(triaelem,bed,xyz_list,gauss_coord);
-		dbdx=db(1);
-		dbdy=db(2);
-
-		%Get the Jacobian determinant
-		Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
-
-		%Get L: 
-		L=GetL(triaelem,gauss_coord,DOFPERGRID);
-
-		%Build gaussian vector 
-		pe_g_basevert_gaussian=Jdettria*gauss_weight*(vx*dbdx+vy*dbdy-melting_gauss)*L';
-
-		%Add Ke_gg_drag_gaussian to Ke
-		pe.terms=pe.terms+pe_g_basevert_gaussian;
-	end
-
-	%for iceshelves, vertical velocity at the base is determined using the hydrostatic equilibrium -> we add a corrective term 
-	%to the icesheet term.
-	if pentaelem.shelf,
-
-		%build the vector thickness*velocity_average
-		HU=thickness.*vxvy_average_list(:,1);
-		HV=thickness.*vxvy_average_list(:,2);
-
-		for ig=1:num_gauss2D,
-	
-			%Pick up the gaussian point and its weight:
-			gauss_weight=gauss_weights(ig);
-			gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
-			
-			%get div(H*Vel)
-			dHU=GetParameterDerivativeValue(triaelem,HU,xyz_list,gauss_coord);
-			dHV=GetParameterDerivativeValue(triaelem,HV,xyz_list,gauss_coord);
-			divHVel=dHU(1)+dHV(2);
-
-			%Get melting and accumulation at gaussian point.
-			melting_gauss=GetParameterValue(triaelem,melting,gauss_coord);
-			accumulation_gauss=GetParameterValue(triaelem,accumulation,gauss_coord);
-
-			%Get L: 
-			L=GetL(triaelem,gauss_coord,DOFPERGRID);
-
-			%Get the Jacobian determinant
-			Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
-
-			%Build gaussian vector 
-			pe_g_basevert_iceshelf_gaussian=Jdettria*gauss_weight*(-di*(-divHVel+accumulation_gauss-melting_gauss))*L';
-			
-			%Add Ke_gg_drag_gaussian to Ke
-			pe.terms=pe.terms+pe_g_basevert_iceshelf_gaussian;
-		end
-	end
-end %end function
-
 function pe=CreatePVectorVert(pentaelem,grids,materials,inputs);
 
@@ -593,4 +446,131 @@
 	    end %for ig=1:num_area_gauss,
 	end %for ig=1:num_vert_gauss,
+
+	if pentaelem.onbed,
+
+		numgrids2=3;
+
+		%Just keep the first 3 grids
+		xyz_list=xyz_list(1:3,:);
+
+		%recover extra inputs
+		[bed_param bed_is_present]=recover_input(inputs,'bed');
+		[melting_param melting_is_present]=recover_input(inputs,'melting');
+		[accumulation_param accumulation_is_present]=recover_input(inputs,'accumulation');
+		[velocity_param velocity_is_present]=recover_input(inputs,'velocity');
+		[velocity_average_param velocity_average_is_present]=recover_input(inputs,'velocity_average');
+		[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
+		if (~velocity_is_present | ~velocity_average_is_present),
+			error('CreatePVector error message: no input velocity!');
+		end
+
+		%initialize velocity array
+		vxvy_list=zeros(numgrids2,2);
+		vxvy_average_list=zeros(numgrids2,2);
+		bed_list=zeros(numgrids2,1);
+		melting_list=zeros(numgrids2,1);
+		accumulation_list=zeros(numgrids2,1);
+		thickness_list=zeros(numgrids2,1);
+
+		%Build linear indices for elementary stiffness matrix.
+		for i=1:numgrids2,
+			doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
+			for j=1:2,
+				dof=doflist(j);
+				vxvy_list(i,j)=velocity_param(dof);
+				vxvy_average_list(i,j)=velocity_average_param(dof);
+			end
+			if(bed_is_present) bed_list(i)=bed_param(doflist(1));end;
+			if(melting_is_present) melting_list(i)=melting_param(doflist(1));end;
+			if(accumulation_is_present) accumulation_list(i)=accumulation_param(doflist(1));end;
+			if(thickness_is_present) thickness_list(i)=thickness_param(doflist(1));end;
+		end
+
+		if bed_is_present, bed=bed_list(1:3); else bed=pentaelem.b(1:3); end
+		if thickness_is_present, thickness=thickness_list(1:3); else thickness=pentaelem.h(1:3); end
+		if melting_is_present, melting=melting_list(1:3); else melting=pentaelem.melting(1:3); end
+		if accumulation_is_present, accumulation=accumulation_list(1:3); else accumulation=pentaelem.accumulation(1:3); end
+
+		% Get gaussian points and weights.
+		[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
+
+		%for icesheets
+		for ig=1:num_gauss2D,
+		
+			%Pick up the gaussian point and its weight:
+			gauss_weight=gauss_weights(ig);
+			gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
+			
+			%Get melting at gaussian point.
+			melting_gauss=GetParameterValue(triaelem,melting,gauss_coord);
+			
+			%Get velocity at gaussian point.
+			vx=GetParameterValue(triaelem,vxvy_list(:,1),gauss_coord);
+			vy=GetParameterValue(triaelem,vxvy_list(:,2),gauss_coord);
+
+			%Get bed slope
+			db=GetParameterDerivativeValue(triaelem,bed,xyz_list,gauss_coord);
+			dbdx=db(1);
+			dbdy=db(2);
+
+			%Get the Jacobian determinant
+			Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
+
+			%Get L: 
+			L=GetL(triaelem,gauss_coord,NDOF1);
+
+			%Build gaussian vector 
+			pe_g_basevert_gaussian=-Jdettria*gauss_weight*(vx*dbdx+vy*dbdy-melting_gauss)*L';
+
+			%Add Ke_gg_drag_gaussian to Ke
+			pe.terms(1:3)=pe.terms(1:3)+pe_g_basevert_gaussian;
+		end
+
+		%for iceshelves, vertical velocity at the base is determined using the hydrostatic equilibrium -> we add a corrective term 
+		%to the icesheet term.
+		if pentaelem.shelf,
+
+			%recover material parameters
+			matice=materials(pentaelem.matid).material;
+			matpar=materials(end).constants;
+			rho_ice=matpar.rho_ice;
+			rho_water=matpar.rho_water;
+			di=rho_ice/rho_water;
+
+			%build the vector thickness*velocity_average
+			HU=thickness.*vxvy_average_list(:,1);
+			HV=thickness.*vxvy_average_list(:,2);
+
+			for ig=1:num_gauss2D,
+		
+				%Pick up the gaussian point and its weight:
+				gauss_weight=gauss_weights(ig);
+				gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
+				
+				%get div(H*Vel)
+				dHU=GetParameterDerivativeValue(triaelem,HU,xyz_list,gauss_coord);
+				dHV=GetParameterDerivativeValue(triaelem,HV,xyz_list,gauss_coord);
+				divHVel=dHU(1)+dHV(2);
+
+				%Get melting and accumulation at gaussian point.
+				melting_gauss=GetParameterValue(triaelem,melting,gauss_coord);
+				accumulation_gauss=GetParameterValue(triaelem,accumulation,gauss_coord);
+
+				%Get L: 
+				L=GetL(triaelem,gauss_coord,NDOF1);
+
+				%Get the Jacobian determinant
+				Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
+
+				%Build gaussian vector 
+				pe_g_basevert_iceshelf_gaussian=-Jdettria*gauss_weight*(-di*(-divHVel+accumulation_gauss-melting_gauss))*L';
+				
+				%Add Ke_gg_drag_gaussian to Ke
+				pe.terms(1:3)=pe.terms(1:3)+pe_g_basevert_iceshelf_gaussian;
+			end
+		end
+
+
+	end
 
 end %end function
