function [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorDiagnosticHutter(md);
%MODELPROCESSORDIAGNOSTICHUTTER - process model for a Hutter 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]=ModelProcessorDiagnosticHutter(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'),
	pos=find(md.gridonhutter);
	count=[1:length(pos)]';

	elements=struct('element',cell(length(pos),1));
	materials=struct('material',cell(length(pos)+1,1));
	mygrids=zeros(md.numberofgrids,1); %this array determines grid partitioning.
	[elements(count).element]=deal(singelem);

	elements(count)=SetStructureField(elements(count),'element','type','singelem');
	elements(count)=SetStructureField(elements(count),'element','id',count);
	elements(count)=SetStructureField(elements(count),'element','g',pos);
	elements(count)=SetStructureField(elements(count),'element','h',md.thickness(pos));
	elements(count)=SetStructureField(elements(count),'element','k',md.drag(pos));
	elements(count)=SetStructureField(elements(count),'element','matid',count);

	[materials(count).material]=deal(matice);
	materials(count)=SetStructureField(materials(count),'material','id',count);
	materials(count)=SetStructureField(materials(count),'material','B',md.B(pos));
	materials(count)=SetStructureField(materials(count),'material','n',3); %TO BE CHANGED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

else
	pos=find(1-md.gridonsurface);
	grids_el=[pos md.uppergrids(pos)];
	count=1:length(pos);

	elements=struct('element',cell(length(pos),1));
	materials=struct('material',cell(length(pos)+1,1));
	mygrids=zeros(md.numberofgrids,1); %this array determines grid partitioning.
	[elements(count).element]=deal(beamelem);

	elements(count)=SetStructureField(elements(count),'element','type','beamelem');
	elements(count)=SetStructureField(elements(count),'element','id',count');
	elements(count)=SetStructureField(elements(count),'element','g',grids_el(count,1:2));
	elements(count)=SetStructureField(elements(count),'element','h',md.thickness(grids_el));
	elements(count)=SetStructureField(elements(count),'element','s',md.surface(grids_el));
	elements(count)=SetStructureField(elements(count),'element','b',md.bed(grids_el));
	elements(count)=SetStructureField(elements(count),'element','k',md.drag(grids_el));
	elements(count)=SetStructureField(elements(count),'element','onbed',md.gridonbed(grids_el(:,1)));
	elements(count)=SetStructureField(elements(count),'element','onsurface',md.gridonsurface(grids_el(:,2)));
	elements(count)=SetStructureField(elements(count),'element','matid',count');

	[materials(count).material]=deal(matice);
	materials(count)=SetStructureField(materials(count),'material','id',count');
	materials(count)=SetStructureField(materials(count),'material','B',md.B(grids_el)*[1;1]/2);
	materials(count)=SetStructureField(materials(count),'material','n',3); %TO BE CHANGED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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','rho_ice',md.rho_ice);
materials(end)=SetStructureField(materials(end),'constants','beta',md.beta);
materials(end)=SetStructureField(materials(end),'constants','meltingpoint',md.meltingpoint);

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

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','3456');

%Boundary conditions:

loads=struct('load',cell(0,0));

%deal with mpcs for 2d-3d mesh transitions
if strcmpi(md.type,'3d'),
	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
end


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

count=1;
for i=1:md.numberofgrids,
	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;

	end

end


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