function [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorDiagnosticHoriz(md);
%MODELPROCESSORDIAGNOSTICHORIZ - process model for an horizontal diagnostic solution
%
%   This routine uses all the informations in the model md and put them
%   into different structures (grids, elements, loads, constrained,materials)
%   that will be used to create the stiffness matrix and load vector.
%   After this routine, the model md should not be called until the end of the
%   solution sequence.
%
%   Usage:
%      [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorDiagnosticHoriz(md)

global cluster

if cluster,
	%We are running in parallel, we need to partition the elements 
	element_partitioning=MeshPartition(md,numlabs);
else
	%We are running in serial, all elements belong to the same partition.
	element_partitioning=ones(md.numberofelements,1);
	labindex=1; %older versions of matlab do not include the parallel toolbox labindex variable.
end

%Allocate grids and elements
numberofstokeselements=size(find(md.elements_type(:,2)==stokesenum()),1);
elements=struct('element',cell(numberofstokeselements,1));
materials=struct('material',cell(numberofstokeselements+1,1));
mygrids=zeros(md.numberofgrids,1); %this array determines grid partitioning.

%3d elements
%pos=find(element_partitioning==labindex);
pos=find(md.elements_type(:,2)==stokesenum());
stokesnewnumber=[1:numberofstokeselements]';
[elements(stokesnewnumber).element]=deal(pentaelem);

elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','type','pentaelem');
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','id',stokesnewnumber);
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','g',md.elements(pos,1:6));
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','h',md.thickness(md.elements(pos,1:6)));
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','s',md.surface(md.elements(pos,1:6)));
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','b',md.bed(md.elements(pos,1:6)));
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','friction_type',md.drag_type);
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','k',md.drag(md.elements(pos,1:6)));
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','p',md.p(pos));
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','q',md.q(pos));
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','shelf',md.elementoniceshelf(pos));
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','onbed',md.elementonbed(pos));
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','onsurface',md.elementonsurface(pos));
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','meanvel',md.meanvel);
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','epsvel',md.epsvel);
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','acceleration',0);%acceleration 0 since no acceleration is posssible
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','collapse',0);
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','matid',stokesnewnumber);
elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','reconditioning_number',md.stokesreconditioning);

[materials(stokesnewnumber).material]=deal(matice);

materials(stokesnewnumber)=SetStructureField(materials(stokesnewnumber),'material','id',stokesnewnumber);
materials(stokesnewnumber)=SetStructureField(materials(stokesnewnumber),'material','B',md.B(md.elements(pos,1:6))*[1;1;1;1;1;1]/6);
materials(stokesnewnumber)=SetStructureField(materials(stokesnewnumber),'material','n',md.n(pos));


%Add physical constants in materials
[materials(end).constants]=deal(matpar);
materials(end)=SetStructureField(materials(end),'constants','g',md.g);
materials(end)=SetStructureField(materials(end),'constants','viscosity_overshoot',md.viscosity_overshoot);
materials(end)=SetStructureField(materials(end),'constants','rho_ice',md.rho_ice);
materials(end)=SetStructureField(materials(end),'constants','rho_water',md.rho_water);
materials(end)=SetStructureField(materials(end),'constants','thermalconductivity',md.thermalconductivity);
materials(end)=SetStructureField(materials(end),'constants','heatcapacity',md.heatcapacity);
materials(end)=SetStructureField(materials(end),'constants','latentheat',md.latentheat);
materials(end)=SetStructureField(materials(end),'constants','beta',md.beta);
materials(end)=SetStructureField(materials(end),'constants','meltingpoint',md.meltingpoint);
materials(end)=SetStructureField(materials(end),'constants','mixed_layer_capacity',md.mixed_layer_capacity);
materials(end)=SetStructureField(materials(end),'constants','thermal_exchange_velocity',md.thermal_exchange_velocity);


if cluster, 
	%For elements, the corresponding grids belong to this cpu. Keep track of it. 
	mygrids(md.elements(pos,:))=1;
end

if cluster, 
	%Figure out which grids from the partitioning belong to different element partitions. We'll 
	%call them 'border' grids.
	bordergrids=double(gplus(mygrids)>1);
else
	bordergrids=zeros(md.numberofgrids,1); %no partitioning serially.
end

%Get the grids set up:
grids=struct('grid',cell(md.numberofgrids,1));

pos=[1:md.numberofgrids]';
[grids(pos).grid]=deal(node);
grids(pos)=SetStructureField(grids(pos),'grid','id',pos);
grids(pos)=SetStructureField(grids(pos),'grid','x',md.x(pos));
grids(pos)=SetStructureField(grids(pos),'grid','y',md.y(pos));
grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos));
grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos));
grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos));

%spc degrees of freedom:
%	 for each grid, 6 degrees of freedom are allowed. These dofs are numbered from 1 to 6. The first 3
%    deal with the (x,y,z) velocities, or deformations. The last 3 deal with the (x,y,z) rotations.
%    If a certain degree of freedom i (1<=i<=6) is constrained to the value 0, the number i should be added to the
%    gridset field of a grid.
%    The gridset field holds all the numbers corresponding to the dofs that have been constrained to 0 value. Because
%    we do not know firshand how many dofs have been constrained for a certain grid, we need a flexible way
%    of keeping track of these constraints. Hence gridset is a string array, of no given size, with no given
%    numerical order.
%    Ex: if a grid has 0 values for the x and z deformations, and 0 values for the y rotation, we could add any of the
%    following strings to the gridset: '135', '153', '315', etc ...
grids(pos)=SetStructureField(grids(pos),'grid','gridset','56');

%Build grids on the border of stokes
%Get border grids
stokes_elements=find(md.elements_type(:,2)==stokesenum()); %find the elements on the stokes domain
borderflags=zeros(md.numberofgrids,1); 
borderflags(md.elements(stokes_elements,:))=1; %find all the grids of the elements on stokes domain, ie stokes grids and borderstokes
borderstokes=borderflags-md.gridonstokes; %remove stokes grids from this list

%spc the non stokes grids to 0 for all the dof
for n=1:md.numberofgrids,
	if md.gridonstokes(n)==0, %the grids does not belong to the stokes part of the mesh
		grids(n).grid.gridset='123456'; %the six dofs are frozen
	end
	if borderstokes(n),
		grids(n).grid.gridset='12356';
	end
end

%Boundary conditions:

%icefront
segmentonneumann_diag_stokes=md.segmentonneumann_diag(find(md.elements_type(md.segmentonneumann_diag(:,end),2)),:);
length_segmentonneumann_diag_stokes= size(segmentonneumann_diag_stokes,1);
loads=struct('load',cell(length_segmentonneumann_diag_stokes,1));

for i=1:length_segmentonneumann_diag_stokes,

	if (element_partitioning(segmentonneumann_diag_stokes(i,end))~=labindex), %this load does not belong to this cpu element partition.
		continue;
	end
	
	%build a quad ice front for the penta element
	loads(i).load=icefront;
	loads(i).load.eid=find(segmentonneumann_diag_stokes(i,end)==find(md.elements_type(:,2)==stokesenum())); %elements contain only stokes elements, so we have to renumbered
	loads(i).load.g(1)=segmentonneumann_diag_stokes(i,1);
	loads(i).load.g(2)=segmentonneumann_diag_stokes(i,2);
	loads(i).load.g(3)=segmentonneumann_diag_stokes(i,3);
	loads(i).load.g(4)=segmentonneumann_diag_stokes(i,4);
	loads(i).load.rho_water=md.rho_water;
	loads(i).load.type='quad';
end

%create penalties for grids on the base of icesheet. We must have wb=ub*db/dx+vb*db/dy
count=length(loads)+1;
for i=1:md.numberofgrids,
	if grids(i).grid.onbed & md.gridonicesheet(i)  & md.gridonstokes(i),
		pengridobject=pengrid;
		pengridobject.id=count;
		pengridobject.dof=1;
		pengridobject.grid=i;
		pengridobject.penalty_offset=md.penalty_offset;

		loads(count).load=pengridobject;
		count=count+1;
	end
end

%Single point constraints:
spcs=find(md.gridondirichlet_diag);
constraints=struct('constraint',cell(3*length(spcs),1));

count=1;
for i=1:md.numberofgrids,
	if ~md.gridonstokes(i),

		%constrain first dof
		constraints(count).constraint=spc;
		constraints(count).constraint.grid=i;
		constraints(count).constraint.dof=1;
		constraints(count).constraint.value=0; %this will be change to vx in the solution sequences
		count=count+1;

		%constrain second dof
		constraints(count).constraint=spc;
		constraints(count).constraint.grid=i;
		constraints(count).constraint.dof=2;
		constraints(count).constraint.value=0; %this will be change to vy in the solution sequences
		count=count+1;

		%constrain third dof
		constraints(count).constraint=spc;
		constraints(count).constraint.grid=i;
		constraints(count).constraint.dof=3;
		constraints(count).constraint.value=0; %this will be change to vz in the solution sequences
		count=count+1;
	end
end


%Last thing, return a partitioning vector, and its transpose.
[part,tpart]=PartitioningVector(md,grids);


end %end function
