%MESH3DTETRAS class definition
%
%   Usage:
%      mesh=mesh3dtetras();

classdef mesh3dtetras
	properties (SetAccess=public) 
		x                           = NaN;
		y                           = NaN;
		z                           = NaN
		elements                    = NaN
		numberoflayers              = 0;
		numberofelements            = 0;
		numberofvertices            = 0;

		lat                         = NaN
		long                        = NaN
		hemisphere                  = NaN

		vertexonbase                 = NaN
		vertexonsurface             = NaN
		lowerelements               = NaN
		lowervertex                 = NaN
		upperelements               = NaN
		uppervertex                 = NaN
		vertexonboundary            = NaN

		vertexconnectivity          = NaN
		elementconnectivity         = NaN
		average_vertex_connectivity = 0;

		x2d                         = NaN
		y2d                         = NaN
		elements2d                  = NaN
		numberofvertices2d          = 0;
		numberofelements2d          = 0;

		extractedvertices           = NaN
		extractedelements           = NaN
	end
	methods
        function createxml(obj,fid) % {{{
            fprintf(fid, '<!-- 3D Tetra Mesh -->\n');
            
            % Elements and verticies of the original 2d mesh
            fprintf(fid,'%s\n%s\n%s\n','       <frame key="1" label="Elements and vertices of the orginal 2d mesh">','       <section name="mesh3dprisms" />');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','          <parameter key ="numberofelements2d" type="',class(obj.numberofelements2d),'" default="',convert2str(obj.numberofelements2d),'">','              <help> number of elements </help>','          </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','          <parameter key ="numberofvertices2d" type="',class(obj.numberofvertices2d),'" default="',convert2str(obj.numberofvertices2d),'">','              <help> number of vertices </help>','          </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','          <parameter key ="elements2d" type="',class(obj.elements2d),'" default="',convert2str(obj.elements2d),'">','              <help> vertex indices of the mesh elements </help>','          </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="x2d" type="',class(obj.x2d),'" default="',convert2str(obj.x2d),'">','     <section name="mesh3dprisms" />','     <help> vertices x coordinate [m] </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="y2d" type="',class(obj.y2d),'" default="',convert2str(obj.y2d),'">','     <section name="mesh3dprisms" />','     <help> vertices y coordinate [m] </help>','  </parameter>');
            fprintf(fid,'%s\n%s\n','       </frame>');
            
            % Elements and vertices of the extruded 3d mesh
            fprintf(fid,'%s\n%s\n%s\n','       <frame key="2" label="Elements and vertices of the orginal 3d mesh">','       <section name="mesh3dprisms" />');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','          <parameter key ="numberofelements" type="',class(obj.numberofelements),'" default="',convert2str(obj.numberofelements),'">','              <help> number of elements </help>','          </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','          <parameter key ="numberofvertices" type="',class(obj.numberofvertices),'" default="',convert2str(obj.numberofvertices),'">','              <help> number of vertices </help>','          </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','          <parameter key ="elements" type="',class(obj.elements),'" default="',convert2str(obj.elements),'">','              <help> vertex indices of the mesh elements </help>','          </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="x" type="',class(obj.x),'" default="',convert2str(obj.x),'">','     <section name="mesh3dprisms" />','     <help> vertices x coordinate [m] </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="y" type="',class(obj.y),'" default="',convert2str(obj.y),'">','     <section name="mesh3dprisms" />','     <help> vertices y coordinate [m] </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="z" type="',class(obj.y),'" default="',convert2str(obj.y),'">','     <section name="mesh3dprisms" />','     <help> vertices z coordinate [m] </help>','  </parameter>');
            fprintf(fid,'%s\n%s\n','       </frame>');
     
            % properties
            fprintf(fid,'%s\n%s\n%s\n','<frame key="3" label="Properties">','<section name="mesh3dprisms" />');             
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="numberoflayers" type="',class(obj.numberoflayers),'" default="',convert2str(obj.numberoflayers),'">','     <section name="mesh3dprisms" />','     <help> number of extrusion layers </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vertexonbase" type="',class(obj.vertexonbase),'" default="',convert2str(obj.vertexonbase),'">','     <section name="mesh3dprisms" />','     <help> lower vertices flags list </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="elementontbase" type="',class(obj.elementontbase),'" default="',convert2str(obj.elementontbase),'">','     <section name="mesh3dprisms" />','     <help> lower elements flags list </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vertexonsurface" type="',class(obj.vertexonsurface),'" default="',convert2str(obj.vertexonsurface),'">','     <section name="mesh3dprisms" />','     <help> upper vertices flags list </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="elementontsurface" type="',class(obj.elementontsurface),'" default="',convert2str(obj.elementontsurface),'">','     <section name="mesh3dprisms" />','     <help> upper elements flags list </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="uppervertex" type="',class(obj.uppervertex),'" default="',convert2str(obj.uppervertex),'">','     <section name="mesh3dprisms" />','     <help> upper vertex list (NaN for vertex on the upper surface) </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="upperelements" type="',class(obj.upperelements),'" default="',convert2str(obj.upperelements),'">','     <section name="mesh3dprisms" />','     <help> upper element list (NaN for element on the upper layer) </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="lowervertex" type="',class(obj.lowervertex),'" default="',convert2str(obj.lowervertex),'">','     <section name="mesh3dprisms" />','     <help> lower vertex list (NaN for vertex on the lower surface) </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="lowerelements" type="',class(obj.lowerelements),'" default="',convert2str(obj.lowerelements),'">','     <section name="mesh3dprisms" />','     <help> element list (NaN for element on the lower layer) </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vertexonboundary" type="',class(obj.vertexonboundary),'" default="',convert2str(obj.vertexonboundary),'">','     <section name="mesh3dprisms" />','     <help> vertices on the boundary of the domain flag list </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vertexconnectivity" type="',class(obj.vertexconnectivity),'" default="',convert2str(obj.vertexconnectivity),'">','     <section name="mesh3dprisms" />','     <help> list of vertices connected to vertex_i </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="elementconnectivity" type="',class(obj.elementconnectivity),'" default="',convert2str(obj.elementconnectivity),'">','     <section name="mesh3dprisms" />','     <help> average number of vertices connected to one vertex </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="average_vertex_connectivity" type="',class(obj.average_vertex_connectivity),'" default="',convert2str(obj.average_vertex_connectivity),'">','     <section name="mesh3dprisms" />','     <help> average number of vertices connected to one vertex </help>','  </parameter>');
            fprintf(fid,'%s\n%s\n','</frame>');
            
            % Extracted model
            fprintf(fid,'%s\n%s\n%s\n','<frame key="4" label="Extracted Model">','<section name="mesh3dprisms" />'); 
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="extractedvertices" type="',class(obj.extractedvertices),'" default="',convert2str(obj.extractedvertices),'">','     <section name="mesh3dprisms" />','     <help> vertices extracted from the model </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="extractedelements" type="',class(obj.extractedelements),'" default="',convert2str(obj.extractedelements),'">','     <section name="mesh3dprisms" />','     <help> elements extracted from the model </help>','  </parameter>');
            fprintf(fid,'%s\n%s\n','</frame>');
            
            % Projection
            fprintf(fid,'%s\n%s\n%s\n','<frame key="5" label="Projection">','<section name="mesh3dprisms" />'); 
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="lat" type="',class(obj.lat),'" default="',convert2str(obj.lat),'">','     <section name="mesh3dprisms" />','     <help> vertices latitude [degrees] </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="long" type="',class(obj.long),'" default="',convert2str(obj.long),'">','     <section name="mesh3dprisms" />','     <help> verticies longitude [degrees] </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','  <parameter key ="hemisphere" type="',class(obj.hemisphere),'" default="',convert2str(obj.hemisphere),'">','     <section name="mesh3dprisms" />','     <help> Indicate hemisphere ''n'' or ''s'' </help>','  </parameter>');
            fprintf(fid,'%s\n%s\n','</frame>');
        
        end % }}}cd
		function self = mesh3dtetras(varargin) % {{{
			switch nargin
				case 0
					self=setdefaultparameters(self);
				case 1
					inputstruct=varargin{1};
					list1 = properties('mesh3dtetras');
					list2 = fieldnames(inputstruct);
					for i=1:length(list1)
						fieldname = list1{i};
						if ismember(fieldname,list2),
							self.(fieldname) = inputstruct.(fieldname);
						end
					end
				otherwise
					error('constructor not supported');
			end
		end % }}}
		function obj = setdefaultparameters(obj) % {{{

			%the connectivity is the averaged 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,'fieldname','mesh.x','NaN',1,'size',[md.mesh.numberofvertices 1]);
			md = checkfield(md,'fieldname','mesh.y','NaN',1,'size',[md.mesh.numberofvertices 1]);
			md = checkfield(md,'fieldname','mesh.z','NaN',1,'size',[md.mesh.numberofvertices 1]);
			md = checkfield(md,'fieldname','mesh.elements','NaN',1,'>',0,'values',1:md.mesh.numberofvertices);
			md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements 4]);
			if any(~ismember(1:md.mesh.numberofvertices,sort(unique(md.mesh.elements(:)))));
				md = checkmessage(md,'orphan nodes have been found. Check the mesh outline');
			end
			md = checkfield(md,'fieldname','mesh.numberoflayers','>=',0);
			md = checkfield(md,'fieldname','mesh.numberofelements','>',0);
			md = checkfield(md,'fieldname','mesh.numberofvertices','>',0);
			%no checks for numberofedges lat long and hemisphere
			md = checkfield(md,'fieldname','mesh.vertexonbase','size',[md.mesh.numberofvertices 1],'values',[0 1]);
			md = checkfield(md,'fieldname','mesh.vertexonsurface','size',[md.mesh.numberofvertices 1],'values',[0 1]);
			md = checkfield(md,'fieldname','mesh.z','>=',md.geometry.base-10^-10,'message','''mesh.z'' lower than bedrock');
			md = checkfield(md,'fieldname','mesh.z','<=',md.geometry.surface+10^-10,'message','''mesh.z'' higher than surface elevation');
			md = checkfield(md,'fieldname','mesh.average_vertex_connectivity','>=',24,'message','''mesh.average_vertex_connectivity'' should be at least 24 in 3d');
		end % }}}
		function disp(obj) % {{{
			disp(sprintf('   3D tetra Mesh:')); 

			disp(sprintf('\n      Elements and vertices of the original 2d mesh:'));
			fielddisplay(obj,'numberofelements2d','number of elements');
			fielddisplay(obj,'numberofvertices2d','number of vertices');
			fielddisplay(obj,'elements2d','vertex indices of the mesh elements');
			fielddisplay(obj,'x2d','vertices x coordinate [m]');
			fielddisplay(obj,'y2d','vertices y coordinate [m]');

			disp(sprintf('\n      Elements and vertices of the extruded 3d mesh:'));
			fielddisplay(obj,'numberofelements','number of elements');
			fielddisplay(obj,'numberofvertices','number of vertices');
			fielddisplay(obj,'elements','vertex indices of the mesh elements');
			fielddisplay(obj,'x','vertices x coordinate [m]');
			fielddisplay(obj,'y','vertices y coordinate [m]');
			fielddisplay(obj,'z','vertices z coordinate [m]');

			disp(sprintf('\n      Properties:'));
			fielddisplay(obj,'numberoflayers','number of extrusion layers');
			fielddisplay(obj,'vertexonbase','lower vertices flags list');
			fielddisplay(obj,'vertexonsurface','upper vertices flags list');
			fielddisplay(obj,'uppervertex','upper vertex list (NaN for vertex on the upper surface)');
			fielddisplay(obj,'upperelements','upper element list (NaN for element on the upper layer)');
			fielddisplay(obj,'lowervertex','lower vertex list (NaN for vertex on the lower surface)');
			fielddisplay(obj,'lowerelements','lower element list (NaN for element on the lower layer');
			fielddisplay(obj,'vertexonboundary','vertices on the boundary of the domain flag list');

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

			disp(sprintf('\n      Extracted model:'));
			fielddisplay(obj,'extractedvertices','vertices extracted from the model');
			fielddisplay(obj,'extractedelements','elements extracted from the model');

			disp(sprintf('\n      Projection:'));
			fielddisplay(obj,'lat','vertices latitude [degrees]');
			fielddisplay(obj,'long','vertices longitude [degrees]');
			fielddisplay(obj,'hemisphere','Indicate hemisphere ''n'' or ''s'' ');
		end % }}}
		function marshall(obj,md,fid) % {{{
			WriteData(fid,'enum',DomainTypeEnum(),'data',StringToEnum(['Domain' domaintype(obj)]),'format','Integer');
			WriteData(fid,'enum',DomainDimensionEnum(),'data',dimension(obj),'format','Integer');
			WriteData(fid,'enum',MeshElementtypeEnum(),'data',StringToEnum(elementtype(obj)),'format','Integer');
			WriteData(fid,'object',obj,'class','mesh','fieldname','x','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'class','mesh','fieldname','y','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'class','mesh','fieldname','z','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'class','mesh','fieldname','elements','format','DoubleMat','mattype',2);
			WriteData(fid,'object',obj,'class','mesh','fieldname','numberoflayers','format','Integer');
			WriteData(fid,'object',obj,'class','mesh','fieldname','numberofelements','format','Integer');
			WriteData(fid,'object',obj,'class','mesh','fieldname','numberofvertices','format','Integer');
			WriteData(fid,'object',obj,'class','mesh','fieldname','vertexonbase','format','BooleanMat','mattype',1);
			WriteData(fid,'object',obj,'class','mesh','fieldname','vertexonsurface','format','BooleanMat','mattype',1);
			WriteData(fid,'object',obj,'class','mesh','fieldname','lowerelements','format','DoubleMat','mattype',2);
			WriteData(fid,'object',obj,'class','mesh','fieldname','upperelements','format','DoubleMat','mattype',2);
			WriteData(fid,'object',obj,'class','mesh','fieldname','average_vertex_connectivity','format','Integer');
			WriteData(fid,'object',obj,'class','mesh','fieldname','elements2d','format','DoubleMat','mattype',3);
			WriteData(fid,'object',obj,'class','mesh','fieldname','numberofvertices2d','format','Integer');
			WriteData(fid,'object',obj,'class','mesh','fieldname','numberofelements2d','format','Integer');
		end % }}}
		function t = domaintype(obj) % {{{
			t = '3D';
		end % }}}
		function d = dimension(obj) % {{{
			d = 3;
		end % }}}
		function s = elementtype(obj) % {{{
			s = 'Tetra';
		end % }}}
	end
end
