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
if strcmpi(md.type,'2d'),
	numberofelements_horiz=length(find(md.elements_type(:,1)==macayealenum()));
else
	numberofelements_horiz=length(find(md.elements_type(:,1)==macayealenum() | md.elements_type(:,1)==pattynenum()));
end
elements=struct('element',cell(numberofelements_horiz,1));
materials=struct('material',cell(numberofelements_horiz+1,1));
mygrids=zeros(md.numberofgrids,1); %this array determines grid partitioning.

%Deal with 2d elements
if strcmpi(md.type,'2d'),
	el3pos=find((md.elements_type(:,1)==macayealenum()) & (element_partitioning==labindex));
	if ~isempty(el3pos),
		pos=[1:length(el3pos)]';
		[elements(pos).element]=deal(triaelem);

		elements(pos)=SetStructureField(elements(pos),'element','type','triaelem');
		elements(pos)=SetStructureField(elements(pos),'element','id',pos);
		elements(pos)=SetStructureField(elements(pos),'element','g',md.elements(el3pos,:));
		elements(pos)=SetStructureField(elements(pos),'element','h',md.thickness(md.elements(el3pos,1:3)));
		elements(pos)=SetStructureField(elements(pos),'element','s',md.surface(md.elements(el3pos,1:3)));
		elements(pos)=SetStructureField(elements(pos),'element','b',md.bed(md.elements(el3pos,1:3)));
		elements(pos)=SetStructureField(elements(pos),'element','friction_type',md.drag_type);
		elements(pos)=SetStructureField(elements(pos),'element','k',md.drag(md.elements(el3pos,1:3)));
		elements(pos)=SetStructureField(elements(pos),'element','p',md.p(el3pos));
		elements(pos)=SetStructureField(elements(pos),'element','q',md.q(el3pos));
		elements(pos)=SetStructureField(elements(pos),'element','shelf',md.elementoniceshelf(el3pos));
		elements(pos)=SetStructureField(elements(pos),'element','meanvel',md.meanvel);
		elements(pos)=SetStructureField(elements(pos),'element','epsvel',md.epsvel);
		elements(pos)=SetStructureField(elements(pos),'element','acceleration',md.acceleration);
		elements(pos)=SetStructureField(elements(pos),'element','matid',pos);


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

		materials(pos)=SetStructureField(materials(pos),'material','id',pos);
		materials(pos)=SetStructureField(materials(pos),'material','B',md.B(md.elements(el3pos,1:3))*[1;1;1]/3);
		materials(pos)=SetStructureField(materials(pos),'material','n',md.n(el3pos));
	end
else

	%3d elements
	%First create penta elements (ie non hutter)
	el6pos=find((md.elements_type(:,1)==macayealenum() | md.elements_type(:,1)==pattynenum()) & (element_partitioning==labindex));

	if ~isempty(el6pos),
		pos=[1:length(el6pos)]'; %Renumber the elements
		[elements(pos).element]=deal(pentaelem);

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

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

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

		%For penta elements where we want to implement MacAyeal's element, we need to collapse 
		%the formulation into trias: 
		pos=find(ismember(el6pos,find(md.elements_type(:,1)==macayealenum())));
		elements(pos)=SetStructureField(elements(pos),'element','collapse',ones(length(pos),1));

	end

end


%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(el3pos,:))=1;
	mygrids(md.elements(el6pos,:))=1;
end

%Deal with acceleration MacAyeal's element. 
if md.acceleration & strcmpi(md.type,'2d'), 
	%build accelerated element.
	elements(end+1).element=acceleratedtriaelem;
	elements(end).element.type='acceleratedtriaelem';
	elements(end).element.x=md.x;
	elements(end).element.y=md.y;
	elements(end).element.index=md.elements;
	elements(end).element.nods=md.numberofgrids;
	elements(end).element.nel=md.numberofelements;
	elements(end).element.thickness=md.thickness;
	elements(end).element.surface=md.surface;
	elements(end).element.bed=md.bed;
	elements(end).element.thickness_el=(md.thickness(md.elements))*[1;1;1]/3;
	elements(end).element.B_bar=(md.B(md.elements))*[1;1;1]/3;
	elements(end).element.glen_coeff=md.n;
	elements(end).element.friction_type=md.drag_type;
	elements(end).element.drag=md.drag;
	elements(end).element.p=md.p;
	elements(end).element.q=md.q;
	elements(end).element.meanvel=md.meanvel;
	elements(end).element.epsvel=md.epsvel;
	if isempty(md.segmentonneumann_diag),
		elements(end).element.index_icefront=[];	
	else
		elements(end).element.index_icefront=md.segmentonneumann_diag(:,1:2);
	end
	elements(end).element.gridoniceshelf=md.gridoniceshelf;
	elements(end).element.elementonicesheet=md.elementonicesheet;
	elements(end).element.matid=1; %point to first matid element.
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 ...
if strcmpi(md.type,'3d')
	grids(pos)=SetStructureField(grids(pos),'grid','gridset','123456');
else
	grids(pos)=SetStructureField(grids(pos),'grid','gridset','3456');
end

%spc macayeal type grids  in 3d meshes
if strcmpi(md.type,'3d'),
	for n=1:length(elements),
		if md.elements_type(el6pos(n),1)==macayealenum(),
			%all dofs associated to this element are frozen, except if the element is on the bed, and acts as a 'fake' pentaelem, 
			%and a true 'macayeal' tria element.
			if elements(n).element.onbed,
				for j=1:3,
					grids(elements(n).element.g(j)).grid.gridset='3456';
				end
			end
		else
			for j=1:length(elements(n).element.g),
				grids(elements(n).element.g(j)).grid.gridset='3456';
			end
		end
	end
end

%Boundary conditions:

%icefront
loads=struct('load',cell(length(md.segmentonneumann_diag),1));

for i=1:size(md.segmentonneumann_diag,1),

	if (element_partitioning(md.segmentonneumann_diag(i,end))~=labindex), %this load does not belong to this cpu element partition.
		continue;
	end

	if strcmpi(md.type,'3d'),
		if md.elements_type(md.segmentonneumann_diag(i,end))==macayealenum(),
			loads(i).load=icefront;
			loads(i).load.eid=find(el6pos==md.segmentonneumann_diag(i,end));
			loads(i).load.g(1)=md.segmentonneumann_diag(i,1);
			loads(i).load.g(2)=md.segmentonneumann_diag(i,2);
			loads(i).load.rho_water=md.rho_water;
			loads(i).load.type='segment';

		elseif md.elements_type(md.segmentonneumann_diag(i,end))==pattynenum(),
			%build a quad ice front for the penta element
			loads(i).load=icefront;
			loads(i).load.eid=find(el6pos==md.segmentonneumann_diag(i,end));
			loads(i).load.g(1)=md.segmentonneumann_diag(i,1);
			loads(i).load.g(2)=md.segmentonneumann_diag(i,2);
			loads(i).load.g(3)=md.segmentonneumann_diag(i,3);
			loads(i).load.g(4)=md.segmentonneumann_diag(i,4);
			loads(i).load.type='quad';

		else
			if ~(md.elements_type(md.segmentonneumann_diag(i,end))==hutterenum()),
				error('diagnostic error message: unsupported  element type');
			end
		end
	else
		if md.elements_type(md.segmentonneumann_diag(i,end))==macayealenum(),
			loads(i).load=icefront;
			loads(i).load.eid=find(el3pos==md.segmentonneumann_diag(i,end));
			loads(i).load.g(1)=md.segmentonneumann_diag(i,1);
			loads(i).load.g(2)=md.segmentonneumann_diag(i,2);
			loads(i).load.rho_water=md.rho_water;
			loads(i).load.type='segment';

		else
			if ~(md.elements_type(md.segmentonneumann_diag(i,end))==hutterenum()),
				error('diagnostic error message: unsupported  element type');
			end
		end
	
	end

end

%Initialize constraints structure
constraints=struct('constraint',cell(2*length(find(md.gridondirichlet_diag)),1));

count=1;
for i=1:md.numberofgrids,
	%deal with spcs
	if md.gridonhutter(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;

	elseif (md.gridonmacayeal(i) | md.gridonpattyn(i)) &  md.gridondirichlet_diag(i),

		%constrain first dof
		constraints(count).constraint=spc;
		constraints(count).constraint.grid=i;
		constraints(count).constraint.dof=1;
		constraints(count).constraint.value=md.dirichletvalues_diag(i,1)/md.yts; %in m/s
		count=count+1;

		%constrain second dof
		constraints(count).constraint=spc;
		constraints(count).constraint.grid=i;
		constraints(count).constraint.dof=2;
		constraints(count).constraint.value=md.dirichletvalues_diag(i,2)/md.yts; %in m/s
		count=count+1;
	end
end

%deal with mpcs
if ~isempty(md.penalties) & ~isnan(md.penalties),
	for i=1:size(md.penalties,1),
		for j=1:(md.numlayers-1),

			%constrain first dof
			constraints(count).constraint=rgb;
			constraints(count).constraint.grid1=md.penalties(i,1);
			constraints(count).constraint.grid2=md.penalties(i,j+1);
			constraints(count).constraint.dof=1;
			count=count+1;
			
			%constrain second dof
			constraints(count).constraint=rgb;
			constraints(count).constraint.grid1=md.penalties(i,1);
			constraints(count).constraint.grid2=md.penalties(i,j+1);
			constraints(count).constraint.dof=2;
			count=count+1;
		end
	end
end

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

end %end function
