Index: /issm/trunk-jpl/src/m/classes/mesh3dtetras.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh3dtetras.m (revision 17483) +++ /issm/trunk-jpl/src/m/classes/mesh3dtetras.m (revision 17483) @@ -0,0 +1,152 @@ +%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 + + 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 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.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,'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',MeshTypeEnum(),'data',Mesh3DtetrasEnum,'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','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/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 17482) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 17483) @@ -1103,4 +1103,49 @@ end %}}} + function md = tetras(md,varargin) % {{{ + %TETRAS - split 3d prismatic mesh into 3 tetrahedrons + % + % Usage: + % md=tetra(md) + + if ~isa(md.mesh,'mesh3dprisms') + error('mesh is not a 3d prismatic mesh'); + end + + %Initialize tetra mesh + md.mesh=mesh3dtetras(md.mesh); + + %Split in 3 tetras + subelement1 = [1 2 3 5]; + subelement2 = [4 6 5 1]; + subelement3 = [5 6 3 1]; + elements=[md.mesh.elements(:,subelement1);md.mesh.elements(:,subelement2);md.mesh.elements(:,subelement3)]; + pos_elements = repmat([1:md.mesh.numberofelements]',3,1); + + md.mesh.elements=elements; + md.mesh.numberofelements=size(elements,1); + + %p and q (same deal, except for element that are on the bedrock: ) + md.friction.p=md.friction.p(pos_elements); + md.friction.q=md.friction.q(pos_elements); + + %elementstype + if ~isnan(md.flowequation.element_equation) + oldelements_type=md.flowequation.element_equation; + md.flowequation.element_equation=zeros(number_el3d,1); + md.flowequation.element_equation=md.flowequation.element_equation(pos_elements); + end + + %connectivity + md.mesh.elementconnectivity=NaN; + + %materials + md.materials.rheology_n=md.materials.rheology_n(pos_elements); + + %increase connectivity if less than 25: + if md.mesh.average_vertex_connectivity<=25, + md.mesh.average_vertex_connectivity=100; + end + end % }}} function disp(obj) % {{{ disp(sprintf('%19s: %-22s -- %s','mesh' ,['[1x1 ' class(obj.mesh) ']'],'mesh properties'));