Context Navigation

mesh3dsurface.m@ 25129

Last change on this file since 25129 was 25129, checked in by jdquinn, 5 years ago
CHG: Saving progress on MATLAB to Python translations (still need to track down marshalling bug on Python side)
File size: 11.5 KB

Rev	Line
[19076]	1	%MESH3DSURFACE class definition
	2	%
	3	% Usage:
	4	% mesh3dsurface=mesh3dsurface();
	5
	6	classdef mesh3dsurface
	7	properties (SetAccess=public)
	8	x = NaN;
	9	y = NaN;
	10	z = NaN;
	11	elements = NaN;
	12	numberofelements = 0;
	13	numberofvertices = 0;
	14	numberofedges = 0;
	15
	16	lat = NaN;
	17	long = NaN;
	18	r = NaN;
[24156]	19	area = NaN;
[19076]	20
	21	vertexonboundary = NaN;
	22	edges = NaN;
	23	segments = NaN;
	24	segmentmarkers = NaN;
	25	vertexconnectivity = NaN;
	26	elementconnectivity = NaN;
	27	average_vertex_connectivity = 0;
	28
	29	extractedvertices = NaN;
	30	extractedelements = NaN;
	31	end
	32	methods (Static)
	33	function self = loadobj(self) % {{{
	34	% This function is directly called by matlab when a model selfect is
	35	% loaded. Update old properties here
	36
	37	%2014 Oct. 1st
	38	if isstruct(self),
	39	oldself=self;
	40	%Assign property values from struct
	41	self=structtoobj(mesh3dsurface(),oldself);
[24778]	42	%if isfield(oldself,'hemisphere'), %NOTE: not ever needed, epsg is not a field anymore.
	43	% disp('md.mesh.hemisphere has been automatically converted to EPSG code');
	44	% if strcmpi(oldself.hemisphere,'n'),
	45	% self.epsg=3413;
	46	% else
	47	% self.epsg=3031;
	48	% end
	49	%end
[19076]	50	end
	51
	52	end% }}}
	53	end
	54	methods
	55	function self = mesh3dsurface(varargin) % {{{
	56	switch nargin
	57	case 0
	58	self=setdefaultparameters(self);
	59	case 1
	60	self=mesh3dsurface();
	61	object=varargin{1};
	62	fields=fieldnames(object);
	63	for i=1:length(fields)
	64	field=fields{i};
	65	if ismember(field,properties('mesh3dsurface')),
	66	self.(field)=object.(field);
	67	end
	68	end
	69	otherwise
	70	error('constructor not supported');
	71	end
	72	end % }}}
	73	function obj = setdefaultparameters(obj) % {{{
	74
[25129]	75	%The connectivity is the average number of nodes linked to a given
	76	%node through an edge. This connectivity is used to initially
	77	%allocate memory to the stiffness matrix. A value of 16 seems to
	78	%give a good memory/time ratio. This value can be checked in
	79	%test/NightlyRun/runme.py.
[19076]	80	obj.average_vertex_connectivity=25;
	81	end % }}}
[20706]	82	function md = checkconsistency(obj,md,solution,analyses) % {{{
[19076]	83
[22004]	84	if strcmpi(solution,'LoveSolution'), return; end
	85
[19897]	86	md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
	87	md = checkfield(md,'fieldname','mesh.y','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
	88	md = checkfield(md,'fieldname','mesh.z','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
	89	md = checkfield(md,'fieldname','mesh.lat','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
	90	md = checkfield(md,'fieldname','mesh.long','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
	91	md = checkfield(md,'fieldname','mesh.r','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
	92	md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',1:md.mesh.numberofvertices);
[19076]	93	md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements 3]);
	94	if any(~ismember(1:md.mesh.numberofvertices,sort(unique(md.mesh.elements(:)))));
[24919]	95	md = checkmessage(md,'orphan nodes have been found; check the mesh outline');
[19076]	96	end
	97	md = checkfield(md,'fieldname','mesh.numberofelements','>',0);
	98	md = checkfield(md,'fieldname','mesh.numberofvertices','>',0);
	99	md = checkfield(md,'fieldname','mesh.average_vertex_connectivity','>=',9,'message','''mesh.average_vertex_connectivity'' should be at least 9 in 2d');
	100
[21049]	101	if strcmp(solution,'ThermalSolution')
	102	md = checkmessage(md,'thermal not supported for 2d mesh');
[19076]	103	end
	104	end % }}}
	105	function disp(obj) % {{{
[22955]	106	disp(sprintf(' 3D tria Mesh (surface):'));
[19076]	107
	108	disp(sprintf('\n Elements and vertices:'));
	109	fielddisplay(obj,'numberofelements','number of elements');
	110	fielddisplay(obj,'numberofvertices','number of vertices');
	111	fielddisplay(obj,'elements','vertex indices of the mesh elements');
	112	fielddisplay(obj,'x','vertices x coordinate [m]');
	113	fielddisplay(obj,'y','vertices y coordinate [m]');
	114	fielddisplay(obj,'z','vertices z coordinate [m]');
	115	fielddisplay(obj,'lat','vertices latitude [degrees]');
	116	fielddisplay(obj,'long','vertices longitude [degrees]');
	117	fielddisplay(obj,'r','vertices radius [m]');
[24156]	118	fielddisplay(obj,'area','elemental areas [m^2]');
[19076]	119
	120	fielddisplay(obj,'edges','edges of the 2d mesh (vertex1 vertex2 element1 element2)');
	121	fielddisplay(obj,'numberofedges','number of edges of the 2d mesh');
	122
	123	disp(sprintf('\n Properties:'));
	124	fielddisplay(obj,'vertexonboundary','vertices on the boundary of the domain flag list');
	125	fielddisplay(obj,'segments','edges on domain boundary (vertex1 vertex2 element)');
	126	fielddisplay(obj,'segmentmarkers','number associated to each segment');
[22879]	127	fielddisplay(obj,'vertexconnectivity','list of elements connected to vertex_i');
	128	fielddisplay(obj,'elementconnectivity','list of elements adjacent to element_i');
[19076]	129	fielddisplay(obj,'average_vertex_connectivity','average number of vertices connected to one vertex');
	130
	131	disp(sprintf('\n Extracted model:'));
	132	fielddisplay(obj,'extractedvertices','vertices extracted from the model');
	133	fielddisplay(obj,'extractedelements','elements extracted from the model');
	134	end % }}}
[20923]	135	function marshall(obj,prefix,md,fid) % {{{
[20970]	136	WriteData(fid,prefix,'name','md.mesh.domain_type','data',['Domain' domaintype(obj)],'format','String');
[20690]	137	WriteData(fid,prefix,'name','md.mesh.domain_dimension','data',dimension(obj),'format','Integer');
[20970]	138	WriteData(fid,prefix,'name','md.mesh.elementtype','data',elementtype(obj),'format','String');
[20690]	139	WriteData(fid,prefix,'object',obj,'fieldname','x','format','DoubleMat','mattype',1);
	140	WriteData(fid,prefix,'object',obj,'fieldname','y','format','DoubleMat','mattype',1);
	141	WriteData(fid,prefix,'object',obj,'fieldname','z','format','DoubleMat','mattype',1);
	142	WriteData(fid,prefix,'object',obj,'fieldname','lat','data',obj.lat,'format','DoubleMat','mattype',1);
	143	WriteData(fid,prefix,'object',obj,'fieldname','long','data',obj.long,'format','DoubleMat','mattype',1);
	144	WriteData(fid,prefix,'object',obj,'fieldname','r','format','DoubleMat','mattype',1);
	145	WriteData(fid,prefix,'name','md.mesh.z','data',zeros(obj.numberofvertices,1),'format','DoubleMat','mattype',1);
	146	WriteData(fid,prefix,'object',obj,'fieldname','elements','format','DoubleMat','mattype',2);
	147	WriteData(fid,prefix,'object',obj,'fieldname','numberofelements','format','Integer');
	148	WriteData(fid,prefix,'object',obj,'fieldname','numberofvertices','format','Integer');
	149	WriteData(fid,prefix,'object',obj,'fieldname','average_vertex_connectivity','format','Integer');
	150	WriteData(fid,prefix,'object',obj,'fieldname','vertexonboundary','format','DoubleMat','mattype',1);
[19076]	151	end % }}}
	152	function t = domaintype(obj) % {{{
	153	t = '3Dsurface';
	154	end % }}}
	155	function d = dimension(obj) % {{{
	156	d = 2;
	157	end % }}}
	158	function s = elementtype(obj) % {{{
	159	s = 'Tria';
	160	end % }}}
	161	function [x y z elements is2d isplanet] = processmesh(self,options) % {{{
	162
	163	isplanet = 1;
	164	is2d = 0;
	165
	166	elements = self.elements;
	167	x = self.x;
	168	y = self.y;
	169	z = self.z;
	170	end % }}}
[20181]	171	function savemodeljs(self,fid,modelname) % {{{
	172
[20196]	173	fprintf(fid,'%s.mesh=new mesh3dsurface();\n',modelname);
[20181]	174	writejs1Darray(fid,[modelname '.mesh.x'],self.x);
	175	writejs1Darray(fid,[modelname '.mesh.y'],self.y);
	176	writejs1Darray(fid,[modelname '.mesh.z'],self.z);
	177	writejs2Darray(fid,[modelname '.mesh.elements'],self.elements);
	178	writejsdouble(fid,[modelname '.mesh.numberofelements'],self.numberofelements);
	179	writejsdouble(fid,[modelname '.mesh.numberofvertices'],self.numberofvertices);
	180	writejsdouble(fid,[modelname '.mesh.numberofedges'],self.numberofedges);
	181	writejs1Darray(fid,[modelname '.mesh.lat'],self.lat);
	182	writejs1Darray(fid,[modelname '.mesh.long'],self.long);
	183	writejs1Darray(fid,[modelname '.mesh.r'],self.r);
[24156]	184	writejs1Darray(fid,[modelname '.mesh.area'],self.area);
[20181]	185	writejs1Darray(fid,[modelname '.mesh.vertexonboundary'],self.vertexonboundary);
	186	writejs2Darray(fid,[modelname '.mesh.edges'],self.edges);
	187	writejs2Darray(fid,[modelname '.mesh.segments'],self.segments);
	188	writejs2Darray(fid,[modelname '.mesh.segmentmarkers'],self.segmentmarkers);
	189	writejs2Darray(fid,[modelname '.mesh.vertexconnectivity'],self.vertexconnectivity);
	190	writejs2Darray(fid,[modelname '.mesh.elementconnectivity'],self.elementconnectivity);
	191	writejsdouble(fid,[modelname '.mesh.average_vertex_connectivity'],self.average_vertex_connectivity);
	192	writejs1Darray(fid,[modelname '.mesh.extractedvertices'],self.extractedvertices);
	193	writejs1Darray(fid,[modelname '.mesh.extractedelements'],self.extractedelements);
	194
	195	end % }}}
[24852]	196	function export(self,varargin) % {{{
	197
	198	options=pairoptions(varargin{:});
	199	filename=getfieldvalue(options,'filename');
	200	format=getfieldvalue(options,'format','shp');
	201	geometry=getfieldvalue(options,'geometry','line');
[24885]	202	index=getfieldvalue(options,'index',[]);
	203	proj=getfieldvalue(options,'projection','');
[24852]	204
	205	%prepare contours:
	206	contours= struct([]);
	207	if strcmpi(geometry,'point'),
	208	for i=1:self.numberofvertices,
	209	contours(i).x = self.long(i);
	210	contours(i).y = self.lat(i);
	211	contours(i).id = i;
	212	contours(i).Geometry = 'Point';
	213	end
	214	elseif strcmpi(geometry,'line'),
[24897]	215	count=1;
[24852]	216	for i=1:self.numberofelements,
	217	el=self.elements(i,:);
	218	%first line:
[24897]	219	contours(count).x = [self.long(el(1)) self.long(el(2))];
	220	contours(count).y = [self.lat(el(1)) self.lat(el(2))];
	221	contours(count).Geometry = 'Line';
[24852]	222
	223	%second line:
[24897]	224	contours(count+1).x = [self.long(el(2)) self.long(el(3))];
	225	contours(count+1).y = [self.lat(el(2)) self.lat(el(3))];
	226	contours(count+1).Geometry = 'Line';
[24852]	227
	228	%third line:
[24897]	229	contours(count+2).x = [self.long(el(3)) self.long(el(1))];
	230	contours(count+2).y = [self.lat(el(3)) self.lat(el(1))];
	231	contours(count+2).Geometry = 'Line';
[24852]	232
[24897]	233	%increase count:
[24900]	234	count = count+3;
[24852]	235	end
	236	elseif strcmpi(geometry,'polygon'),
[24894]	237	% TODO: Refactor the following to reduce repeated code, or
	238	% leave as is because it is more readable?
[24885]	239	if isempty(index),
	240	for i=1:self.numberofelements,
	241	el=self.elements(i,:);
	242	contours(i).x=[self.long(el(1)) self.long(el(2)) self.long(el(3)) self.long(el(1))];
	243	contours(i).y=[self.lat(el(1)) self.lat(el(2)) self.lat(el(3)) self.lat(el(1))];
[24894]	244	contours(i).Id = i;
[24885]	245	contours(i).Geometry = 'Polygon';
	246	end
	247	else
	248	for i=1:length(index),
	249	el=self.elements(index(i),:);
	250	contours(i).x=[self.long(el(1)) self.long(el(2)) self.long(el(3)) self.long(el(1))];
	251	contours(i).y=[self.lat(el(1)) self.lat(el(2)) self.lat(el(3)) self.lat(el(1))];
	252	contours(i).id = index(i);
	253	contours(i).Geometry = 'Polygon';
	254	end
[24852]	255	end
	256	else
[24894]	257	error(sprintf('mesh3dsurface ''export'' error message: geometry %s not supported yet (should be ''point'' or ''line'' or ''polygon'')',geometry));
[24852]	258	end
	259
	260	%write file:
	261	if strcmpi(format,'shp'),
	262	shpwrite(contours,filename);
	263	elseif strcmpi(format,'exp'),
	264	expwrite(contours,filename);
	265	else
	266	error(sprintf('mesh3dsurface ''export'' error message: file format %s not supported yet',format));
	267	end
	268
[24885]	269	%write projection file:
	270	if ~isempty(proj),
	271	proj2shpprj(filename,proj);
	272	end
[24852]	273
[24885]	274	%write style file:
	275	applyqgisstyle(filename,'mesh');
[24852]	276
[24885]	277
	278
	279
[24852]	280	end % }}}
[19076]	281	end
	282	end

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: issm/trunk-jpl/src/m/classes/mesh3dsurface.m@ 25129

Download in other formats: