Index: /issm/trunk-jpl/src/m/classes/mesh3dprisms.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh3dprisms.m (revision 17556) +++ /issm/trunk-jpl/src/m/classes/mesh3dprisms.m (revision 17556) @@ -0,0 +1,150 @@ +%MESH3DPRISMS class definition +% +% Usage: +% mesh=mesh3dprisms(); + +classdef mesh3dprisms + properties (SetAccess=public) + x = NaN; + y = NaN; + z = NaN + elements = NaN + numberoflayers = 0; + numberofelements = 0; + numberofvertices = 0; + + lat = NaN + long = NaN + hemisphere = NaN + + elementonbed = NaN + elementonsurface = NaN + vertexonbed = 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 obj = mesh(varargin) % {{{ + switch nargin + case 0 + obj=setdefaultparameters(obj); + 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 6]); + 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.elementonbed','size',[md.mesh.numberofelements 1],'values',[0 1]); + md = checkfield(md,'fieldname','mesh.elementonsurface','size',[md.mesh.numberofelements 1],'values',[0 1]); + md = checkfield(md,'fieldname','mesh.vertexonbed','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.bed-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(' 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,'vertexonbed','lower vertices flags list'); + fielddisplay(obj,'elementonbed','lower elements flags list'); + fielddisplay(obj,'vertexonsurface','upper vertices flags list'); + fielddisplay(obj,'elementonsurface','upper elements 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',MeshTypeEnum(),'data',StringToEnum(['Mesh' meshtype(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','elementonbed','format','BooleanMat','mattype',2); + WriteData(fid,'object',obj,'class','mesh','fieldname','elementonsurface','format','BooleanMat','mattype',2); + WriteData(fid,'object',obj,'class','mesh','fieldname','vertexonbed','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 type = meshtype(obj) % {{{ + type = '3D'; + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/mesh3dprisms.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh3dprisms.py (revision 17556) +++ /issm/trunk-jpl/src/m/classes/mesh3dprisms.py (revision 17556) @@ -0,0 +1,153 @@ +import numpy +from fielddisplay import fielddisplay +from EnumDefinitions import * +from checkfield import * +from MatlabFuncs import * + +class mesh3dprisms(object): + """ + MESH3DPRISMS class definition + + Usage: + mesh3d=mesh3dprisms(); + """ + + def __init__(self): # {{{ + self.x = float('NaN'); + self.y = float('NaN'); + self.z = float('NaN'); + self.elements = float('NaN'); + self.dimension = 0; + self.numberoflayers = 0; + self.numberofelements = 0; + self.numberofvertices = 0; + + self.lat = float('NaN'); + self.long = float('NaN'); + self.hemisphere = float('NaN'); + + self.elementonbed = float('NaN'); + self.elementonsurface = float('NaN'); + self.vertexonbed = float('NaN'); + self.vertexonsurface = float('NaN'); + self.lowerelements = float('NaN'); + self.lowervertex = float('NaN'); + self.upperelements = float('NaN'); + self.uppervertex = float('NaN'); + self.vertexonboundary = float('NaN'); + + self.vertexconnectivity = float('NaN'); + self.elementconnectivity = float('NaN'); + self.average_vertex_connectivity = 0; + + self.x2d = float('NaN'); + self.y2d = float('NaN'); + self.elements2d = float('NaN'); + self.numberofvertices2d = 0; + self.numberofelements2d = 0; + + self.extractedvertices = float('NaN'); + self.extractedelements = float('NaN'); + + #set defaults + self.setdefaultparameters() + #}}} + def __repr__(self): # {{{ + string=" Mesh 3D prisms:" + + string="%s\n%s"%(string,"\n Elements and vertices of the original 2d mesh3dprisms:") + + string="%s\n%s"%(string,fielddisplay(self,"numberofelements2d","number of elements")) + string="%s\n%s"%(string,fielddisplay(self,"numberofvertices2d","number of vertices")) + string="%s\n%s"%(string,fielddisplay(self,"elements2d","vertex indices of the mesh3dprisms elements")) + string="%s\n%s"%(string,fielddisplay(self,"x2d","vertices x coordinate [m]")) + string="%s\n%s"%(string,fielddisplay(self,"y2d","vertices y coordinate [m]")) + + string="%s\n%s"%(string,"\n\n Elements and vertices of the extruded 3d mesh3dprisms:") + string="%s\n%s"%(string,fielddisplay(self,"numberofelements","number of elements")) + string="%s\n%s"%(string,fielddisplay(self,"numberofvertices","number of vertices")) + string="%s\n%s"%(string,fielddisplay(self,"elements","vertex indices of the mesh3dprisms elements")) + string="%s\n%s"%(string,fielddisplay(self,"x","vertices x coordinate [m]")) + string="%s\n%s"%(string,fielddisplay(self,"y","vertices y coordinate [m]")) + string="%s\n%s"%(string,fielddisplay(self,"z","vertices z coordinate [m]")) + + string="%s%s"%(string,"\n\n Properties:") + string="%s\n%s"%(string,fielddisplay(self,"numberoflayers","number of extrusion layers")) + string="%s\n%s"%(string,fielddisplay(self,"vertexonbed","lower vertices flags list")) + string="%s\n%s"%(string,fielddisplay(self,"elementonbed","lower elements flags list")) + string="%s\n%s"%(string,fielddisplay(self,"vertexonsurface","upper vertices flags list")) + string="%s\n%s"%(string,fielddisplay(self,"elementonsurface","upper elements flags list")) + string="%s\n%s"%(string,fielddisplay(self,"uppervertex","upper vertex list (-1 for vertex on the upper surface)")) + string="%s\n%s"%(string,fielddisplay(self,"upperelements","upper element list (-1 for element on the upper layer)")) + string="%s\n%s"%(string,fielddisplay(self,"lowervertex","lower vertex list (-1 for vertex on the lower surface)")) + string="%s\n%s"%(string,fielddisplay(self,"lowerelements","lower element list (-1 for element on the lower layer)")) + string="%s\n%s"%(string,fielddisplay(self,"vertexonboundary","vertices on the boundary of the domain flag list")) + string="%s\n%s"%(string,fielddisplay(self,"vertexconnectivity","list of vertices connected to vertex_i")) + string="%s\n%s"%(string,fielddisplay(self,"elementconnectivity","list of vertices connected to element_i")) + string="%s\n%s"%(string,fielddisplay(self,"average_vertex_connectivity","average number of vertices connected to one vertex")) + + string="%s%s"%(string,"\n\n Extracted model:") + string="%s\n%s"%(string,fielddisplay(self,"extractedvertices","vertices extracted from the model")) + string="%s\n%s"%(string,fielddisplay(self,"extractedelements","elements extracted from the model")) + + string="%s%s"%(string,"\n\n Projection:") + string="%s\n%s"%(string,fielddisplay(self,"lat","vertices latitude [degrees]")) + string="%s\n%s"%(string,fielddisplay(self,"long","vertices longitude [degrees]")) + string="%s\n%s"%(string,fielddisplay(self,"hemisphere","Indicate hemisphere 'n' or 's'")) + return string + #}}} + def setdefaultparameters(self): # {{{ + + #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 + self.average_vertex_connectivity=25 + + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + + md = checkfield(md,'fieldname','mesh.x','NaN',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','mesh.y','NaN',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','mesh.z','NaN',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','mesh.elements','NaN',1,'>',0,'values',numpy.arange(1,md.mesh.numberofvertices+1)) + md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements,6]) + if numpy.any(numpy.logical_not(ismember(numpy.arange(1,md.mesh.numberofvertices+1),md.mesh.elements))): + md.checkmessage("orphan nodes have been found. Check the mesh3dprisms outline") + 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.elementonbed','size',[md.mesh.numberofelements],'values',[0,1]) + md = checkfield(md,'fieldname','mesh.elementonsurface','size',[md.mesh.numberofelements],'values',[0,1]) + md = checkfield(md,'fieldname','mesh.vertexonbed','size',[md.mesh.numberofvertices],'values',[0,1]) + md = checkfield(md,'fieldname','mesh.vertexonsurface','size',[md.mesh.numberofvertices],'values',[0,1]) + md = checkfield(md,'fieldname','mesh.average_vertex_connectivity','>=',24,'message',"'mesh.average_vertex_connectivity' should be at least 24 in 3d") + + return md + # }}} + def meshtype(self): # {{{ + return "3D" + #}}} + def marshall(self,md,fid): # {{{ + WriteData(fid,'enum',MeshTypeEnum(),'data',StringToEnum("Mesh"+self.meshtype())[0],'format','Integer'); + WriteData(fid,'object',self,'class','mesh','fieldname','x','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','mesh','fieldname','y','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','mesh','fieldname','z','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','mesh','fieldname','elements','format','DoubleMat','mattype',2) + WriteData(fid,'object',self,'class','mesh','fieldname','numberoflayers','format','Integer') + WriteData(fid,'object',self,'class','mesh','fieldname','numberofelements','format','Integer') + WriteData(fid,'object',self,'class','mesh','fieldname','numberofvertices','format','Integer') + WriteData(fid,'object',self,'class','mesh','fieldname','elementonbed','format','BooleanMat','mattype',2) + WriteData(fid,'object',self,'class','mesh','fieldname','elementonsurface','format','BooleanMat','mattype',2) + WriteData(fid,'object',self,'class','mesh','fieldname','vertexonbed','format','BooleanMat','mattype',1) + WriteData(fid,'object',self,'class','mesh','fieldname','vertexonsurface','format','BooleanMat','mattype',1) + WriteData(fid,'object',self,'class','mesh','fieldname','lowerelements','format','DoubleMat','mattype',2) + WriteData(fid,'object',self,'class','mesh','fieldname','upperelements','format','DoubleMat','mattype',2) + WriteData(fid,'object',self,'class','mesh','fieldname','average_vertex_connectivity','format','Integer') + WriteData(fid,'object',self,'class','mesh','fieldname','elements2d','format','DoubleMat','mattype',3) + WriteData(fid,'object',self,'class','mesh','fieldname','numberofvertices2d','format','Integer') + WriteData(fid,'object',self,'class','mesh','fieldname','numberofelements2d','format','Integer') + # }}}