%MESH class definition
%
%   Usage:
%      planetmesh=planetmesh();

classdef planetmesh
	properties (SetAccess=public) 
		x                           = NaN;
		y                           = NaN;
		z                           = NaN;
		r                           = NaN;
		theta                       = NaN;
		phi                         = NaN
		elements                    = NaN
		dimension                   = 0;
		numberoflayers              = 0;
		numberofelements            = 0;
		numberofvertices            = 0;

		lat                         = NaN
		long                        = NaN

		vertexconnectivity          = NaN
		elementconnectivity         = NaN
		average_vertex_connectivity = 0;
	end
	methods
		function obj = planetmesh(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				otherwise
					error('constructor not supported');
			end
		end % }}}
		function obj = setdefaultparameters(obj) % {{{

			%the connectivity is the avergaded number of nodes linked to a
			%given node through an edge. This connectivity is used to initially
			%allocate memory to the stiffness matrix. A value of 16 seems to
			%give a good memory/time ration. This value can be checked in
			%trunk/test/Miscellaneous/runme.m
			obj.average_vertex_connectivity=25;
		end % }}}
		function md = checkconsistency(obj,md,solution,analyses) % {{{

			md = checkfield(md,'planetmesh.x','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
			md = checkfield(md,'planetmesh.y','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
			md = checkfield(md,'planetmesh.z','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
			md = checkfield(md,'planetmesh.r','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
			md = checkfield(md,'planetmesh.theta','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
			md = checkfield(md,'planetmesh.phi','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
			md = checkfield(md,'planetmesh.elements','NaN',1,'>',0,'values',1:md.planetmesh.numberofvertices);
			if(md.planetmesh.dimension==2),
				md = checkfield(md,'planetmesh.elements','size',[md.planetmesh.numberofelements 3]);
			else
				md = checkfield(md,'planetmesh.elements','size',[md.planetmesh.numberofelements 6]);
			end
			if any(~ismember(1:md.planetmesh.numberofvertices,sort(unique(md.planetmesh.elements(:)))));
				md = checkmessage(md,'orphan nodes have been found. Check the planetmesh outline');
			end
			md = checkfield(md,'planetmesh.dimension','values',[2 3]);
			md = checkfield(md,'planetmesh.numberoflayers','>=',0);
			md = checkfield(md,'planetmesh.numberofelements','>',0);
			md = checkfield(md,'planetmesh.numberofvertices','>',0);
			%no checks for numberofedges lat long and hemisphere
			if (md.planetmesh.dimension==2),
				md = checkfield(md,'planetmesh.average_vertex_connectivity','>=',9,'message','''planetmesh.average_vertex_connectivity'' should be at least 9 in 2d');
			else
				md = checkfield(md,'planetmesh.average_vertex_connectivity','>=',24,'message','''planetmesh.average_vertex_connectivity'' should be at least 24 in 3d');
			end
			md = checkfield(md,'planetmesh.elementconnectivity','size',[md.planetmesh.numberofelements 3],'NaN',1);

			%Solution specific checks
			switch(solution),
				case PrognosticSolutionEnum(),
					if md.prognostic.stabilization==3,
						md = checkfield(md,'planetmesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d planetmeshes');
					end
				case TransientSolutionEnum(),
					if md.transient.isprognostic & md.prognostic.stabilization==3,
						md = checkfield(md,'planetmesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d planetmeshes');
					end
				case ThermalSolutionEnum(),
					md = checkfield(md,'planetmesh.dimension','values',3,'message','thermal solution only supported on 3d planetmeshes');
			end
		end % }}}
		function disp(obj) % {{{
			disp(sprintf('   Mesh:')); 

			disp(sprintf('\n      Elements and vertices:'));
			fielddisplay(obj,'numberofelements','number of elements');
			fielddisplay(obj,'numberofvertices','number of vertices');
			fielddisplay(obj,'elements','index into (x,y,z), coordinates of the vertices');
			fielddisplay(obj,'x','vertices x coordinate');
			fielddisplay(obj,'y','vertices y coordinate');
			fielddisplay(obj,'z','vertices z coordinate');
			fielddisplay(obj,'r','vertices r coordinate');
			fielddisplay(obj,'theta','vertices theta coordinate');
			fielddisplay(obj,'phi','vertices phi coordinate');

			disp(sprintf('\n      Properties:'));
			fielddisplay(obj,'dimension','planetmesh dimension (2d or 3d)');
			fielddisplay(obj,'numberoflayers','number of extrusion layers');

			fielddisplay(obj,'vertexconnectivity','list of vertices connected to vertex_i');
			fielddisplay(obj,'elementconnectivity','list of vertices connected to element_i');
			fielddisplay(obj,'average_vertex_connectivity','average number of vertices connected to one vertex');

		end % }}}
		function marshall(obj,fid) % {{{
			WriteData(fid,'object',obj,'fieldname','x','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','y','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','z','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','r','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','theta','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','phi','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','elements','format','DoubleMat','mattype',2);
			WriteData(fid,'object',obj,'fieldname','dimension','format','Integer');
			WriteData(fid,'object',obj,'fieldname','numberoflayers','format','Integer');
			WriteData(fid,'object',obj,'fieldname','numberofelements','format','Integer');
			WriteData(fid,'object',obj,'fieldname','numberofvertices','format','Integer');
			WriteData(fid,'object',obj,'fieldname','elementconnectivity','format','DoubleMat','mattype',3);
			WriteData(fid,'object',obj,'fieldname','average_vertex_connectivity','format','Integer');
		end % }}}
	end
end
