Index: /issm/trunk-jpl/src/m/classes/mesh.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh.m (revision 12942) +++ /issm/trunk-jpl/src/m/classes/mesh.m (revision 12943) @@ -57,5 +57,5 @@ function obj = setdefaultparameters(obj) % {{{ - %the connectivity is the avergaded number of nodes linked to a + %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 Index: /issm/trunk-jpl/src/m/classes/mesh.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh.py (revision 12942) +++ /issm/trunk-jpl/src/m/classes/mesh.py (revision 12943) @@ -1,4 +1,8 @@ #module imports +import numpy from fielddisplay import fielddisplay +from checkfield import * +from EnumDefinitions import * +from MatlabFuncs import * class mesh: @@ -65,10 +69,5 @@ else: string="\n%s"%(" Elements and vertices:") - string="%s\n%s"%(string,fielddisplay(obj,"numberofelements","number of elements")) - - - - string="%s\n%s"%(string,fielddisplay(obj,"numberofvertices","number of vertices")) string="%s\n%s"%(string,fielddisplay(obj,"elements","index into (x,y,z), coordinates of the vertices")) @@ -80,5 +79,4 @@ string="%s%s"%(string,"\n Properties:") - string="%s\n%s"%(string,fielddisplay(obj,"dimension","mesh dimension (2d or 3d)")) string="%s\n%s"%(string,fielddisplay(obj,"numberoflayers","number of extrusion layers")) @@ -92,5 +90,4 @@ string="%s\n%s"%(string,fielddisplay(obj,"lowerelements","lower element list (NaN for element on the lower layer")) string="%s\n%s"%(string,fielddisplay(obj,"vertexonboundary","vertices on the boundary of the domain flag list")) - string="%s\n%s"%(string,fielddisplay(obj,"segments","edges on domain boundary (vertex1 vertex2 element)")) string="%s\n%s"%(string,fielddisplay(obj,"segmentmarkers","number associated to each segment")) @@ -100,5 +97,4 @@ string="%s%s"%(string,"\n Extracted model:") - string="%s\n%s"%(string,fielddisplay(obj,"extractedvertices","vertices extracted from the model")) string="%s\n%s"%(string,fielddisplay(obj,"extractedelements","elements extracted from the model")) @@ -114,5 +110,5 @@ # {{{setdefaultparameters - #the connectivity is the avergaded number of nodes linked to a + #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 @@ -124,2 +120,75 @@ #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + + md = checkfield(md,'mesh.x','NaN',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'mesh.y','NaN',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'mesh.z','NaN',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'mesh.elements','NaN',1,'>',0,'values',range(1,md.mesh.numberofvertices+1)) + if md.mesh.dimension==2: + md = checkfield(md,'mesh.elements','size',[md.mesh.numberofelements,3]) + else: + md = checkfield(md,'mesh.elements','size',[md.mesh.numberofelements,6]) + if any(numpy.logical_not(ismember(range(1,md.mesh.numberofvertices+1),md.mesh.elements))): + md = checkmessage(md,"orphan nodes have been found. Check the mesh outline") + md = checkfield(md,'mesh.dimension','values',[2,3]) + md = checkfield(md,'mesh.numberoflayers','>=',0) + md = checkfield(md,'mesh.numberofelements','>',0) + md = checkfield(md,'mesh.numberofvertices','>',0) + #no checks for numberofedges lat long and hemisphere + md = checkfield(md,'mesh.elementonbed','size',[md.mesh.numberofelements],'values',[0,1]) + md = checkfield(md,'mesh.elementonsurface','size',[md.mesh.numberofelements],'values',[0,1]) + md = checkfield(md,'mesh.vertexonbed','size',[md.mesh.numberofvertices],'values',[0,1]) + md = checkfield(md,'mesh.vertexonsurface','size',[md.mesh.numberofvertices],'values',[0,1]) + if md.mesh.dimension==2: + md = checkfield(md,'mesh.average_vertex_connectivity','>=',9,'message',"'mesh.average_vertex_connectivity' should be at least 9 in 2d") + else: + md = checkfield(md,'mesh.average_vertex_connectivity','>=',24,'message',"'mesh.average_vertex_connectivity' should be at least 24 in 3d") + md = checkfield(md,'mesh.elementconnectivity','size',[md.mesh.numberofelements,3],'NaN',1) + + #Solution specific checks + if solution==PrognosticSolutionEnum: + if md.prognostic.stabilization==3: + md = checkfield(md,'mesh.dimension','values',2,'message',"Discontinuous Galerkin only supported for 2d meshes") + md = checkfield(md,'mesh.edges','size',[float('NaN'),4]) + md = checkfield(md,'mesh.edges[:,1:3]','>',0) + elif solution==BalancethicknessSolutionEnum: + if md.balancethickness.stabilization==3: + md = checkfield(md,'mesh.dimension','values',2,'message',"Discontinuous Galerkin only supported for 2d meshes") + md = checkfield(md,'mesh.edges','size',[float('NaN'),4]) + md = checkfield(md,'mesh.edges[:,1:3]','>',0) + elif solution==TransientSolutionEnum: + if md.transient.isprognostic and md.prognostic.stabilization==3: + md = checkfield(md,'mesh.dimension','values',2,'message',"Discontinuous Galerkin only supported for 2d meshes") + md = checkfield(md,'mesh.edges','size',[float('NaN'),4]) + md = checkfield(md,'mesh.edges[:,1:3]','>',0) + elif solution==ThermalSolutionEnum: + md = checkfield(md,'mesh.dimension','values',3,'message','thermal solution only supported on 3d meshes') + + return md + # }}} + + def marshall(self,fid): # {{{ + WriteData(fid,'object',self,'fieldname','x','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'fieldname','y','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'fieldname','z','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'fieldname','elements','format','DoubleMat','mattype',2) + WriteData(fid,'object',self,'fieldname','dimension','format','Integer') + WriteData(fid,'object',self,'fieldname','numberoflayers','format','Integer') + WriteData(fid,'object',self,'fieldname','numberofelements','format','Integer') + WriteData(fid,'object',self,'fieldname','numberofvertices','format','Integer') + WriteData(fid,'object',self,'fieldname','numberofedges','format','Integer') + WriteData(fid,'object',self,'fieldname','elementonbed','format','BooleanMat','mattype',2) + WriteData(fid,'object',self,'fieldname','elementonsurface','format','BooleanMat','mattype',2) + WriteData(fid,'object',self,'fieldname','vertexonbed','format','BooleanMat','mattype',1) + WriteData(fid,'object',self,'fieldname','vertexonsurface','format','BooleanMat','mattype',1) + WriteData(fid,'object',self,'fieldname','lowerelements','format','DoubleMat','mattype',2) + WriteData(fid,'object',self,'fieldname','upperelements','format','DoubleMat','mattype',2) + WriteData(fid,'object',self,'fieldname','edges','format','DoubleMat','mattype',3) + WriteData(fid,'object',self,'fieldname','elementconnectivity','format','DoubleMat','mattype',3) + WriteData(fid,'object',self,'fieldname','average_vertex_connectivity','format','Integer') + WriteData(fid,'object',self,'fieldname','elements2d','format','DoubleMat','mattype',3) + WriteData(fid,'object',self,'fieldname','numberofvertices2d','format','Integer') + WriteData(fid,'object',self,'fieldname','numberofelements2d','format','Integer') + # }}} + Index: /issm/trunk-jpl/src/m/classes/private.m =================================================================== --- /issm/trunk-jpl/src/m/classes/private.m (revision 12942) +++ /issm/trunk-jpl/src/m/classes/private.m (revision 12943) @@ -31,5 +31,5 @@ fielddisplay(obj,'isconsistent','is model self consistent'); fielddisplay(obj,'runtimename','name of the run launched'); - fielddisplay(obj,'bamg','structure with mesh properties construced if bamg is used to mesh the domain'); + fielddisplay(obj,'bamg','structure with mesh properties constructed if bamg is used to mesh the domain'); fielddisplay(obj,'solution','type of solution launched'); Index: /issm/trunk-jpl/src/m/classes/private.py =================================================================== --- /issm/trunk-jpl/src/m/classes/private.py (revision 12942) +++ /issm/trunk-jpl/src/m/classes/private.py (revision 12943) @@ -6,7 +6,8 @@ def __init__(self): # {{{ Properties - self.runtimename = '' - self.bamg = {} - self.solution = ''; + self.isconsistent = True; + self.runtimename = '' + self.bamg = {} + self.solution = ''; #set defaults @@ -17,6 +18,8 @@ # {{{ Display string=' private parameters: do not change' + + string="%s\n%s"%(string,fielddisplay(obj,'isconsistent','is model self consistent')) string="%s\n%s"%(string,fielddisplay(obj,'runtimename','name of the run launched')) - string="%s\n%s"%(string,fielddisplay(obj,'bamg','structure with mesh properties construced if bamg is used to mesh the domain')) + string="%s\n%s"%(string,fielddisplay(obj,'bamg','structure with mesh properties constructed if bamg is used to mesh the domain')) string="%s\n%s"%(string,fielddisplay(obj,'solution','type of solution launched')) return string Index: /issm/trunk-jpl/src/m/model/WriteData.py =================================================================== --- /issm/trunk-jpl/src/m/model/WriteData.py (revision 12942) +++ /issm/trunk-jpl/src/m/model/WriteData.py (revision 12943) @@ -9,9 +9,9 @@ Usage: - WriteData(fid,*args) + WriteData(fid,varargin) """ #process options - options=pairoptions(args) + options=pairoptions(*args) #Get data properties Index: /issm/trunk-jpl/src/m/model/ismodelselfconsistent.py =================================================================== --- /issm/trunk-jpl/src/m/model/ismodelselfconsistent.py (revision 12942) +++ /issm/trunk-jpl/src/m/model/ismodelselfconsistent.py (revision 12943) @@ -17,7 +17,6 @@ #Go through a model fields, check that it is a class, and call checkconsistency -# fields=vars('model') - fields=dir(md) - for field in fields: + fields=vars(md) + for field in fields.iterkeys(): #Some properties do not need to be checked Index: /issm/trunk-jpl/src/m/model/marshall.py =================================================================== --- /issm/trunk-jpl/src/m/model/marshall.py (revision 12942) +++ /issm/trunk-jpl/src/m/model/marshall.py (revision 12943) @@ -26,5 +26,5 @@ fields=vars(md) - for field in fields.interkeys(): + for field in fields.iterkeys(): #Some properties do not need to be marshalled Index: /issm/trunk-jpl/src/m/utils/consistency/checkfield.py =================================================================== --- /issm/trunk-jpl/src/m/utils/consistency/checkfield.py (revision 12942) +++ /issm/trunk-jpl/src/m/utils/consistency/checkfield.py (revision 12943) @@ -8,27 +8,27 @@ CHECKFIELD - check field consistency - Used to check model consistency. - Available options: - - NaN: 1 if check that there is no NaN - - size: [lines cols], NaN for non checked dimensions - - >: greater than provided value - - >=: greater or equal to provided value - - <: smallerthan provided value - - <=: smaller or equal to provided value - - < vec: smallerthan provided values on each vertex - - forcing: 1 if check forcing consistency (size and time) - - values: cell of strings or vector of acceptable values - - numel: list of acceptable number of elements - - cell: 1 if check that is cell - - empty: 1 if check that non empty - - message: overloaded error message + Used to check model consistency. + Available options: + - NaN: 1 if check that there is no NaN + - size: [lines cols], NaN for non checked dimensions + - >: greater than provided value + - >=: greater or equal to provided value + - <: smallerthan provided value + - <=: smaller or equal to provided value + - < vec: smallerthan provided values on each vertex + - forcing: 1 if check forcing consistency (size and time) + - values: cell of strings or vector of acceptable values + - numel: list of acceptable number of elements + - cell: 1 if check that is cell + - empty: 1 if check that non empty + - message: overloaded error message - Usage: - md = checkfield(md,fieldname,options); + Usage: + md = checkfield(md,fieldname,options); - Example: - md = checkfield(md,'mesh.elementonbed','size',[md.mesh.numberofelements 1],'values',[0 1]); - md = checkfield(md,'diagnostic.icefront','size',[NaN 4],'NaN',1); - md = checkfield(md,'diagnostic.icefront(:,end)','values',[0 1 2]); + Example: + md = checkfield(md,'mesh.elementonbed','size',[md.mesh.numberofelements 1],'values',[0 1]); + md = checkfield(md,'diagnostic.icefront','size',[NaN 4],'NaN',1); + md = checkfield(md,'diagnostic.icefront(:,end)','values',[0 1 2]); """ @@ -37,8 +37,9 @@ #get field from model - field=getattr(md,fieldname) +# field=getattr(md,fieldname) + exec("field=md.%s" % fieldname) #check empty - if 'empty' in options: + if options.exist('empty'): if not field: md = md.checkmessage(options.getfieldvalue('message',\ @@ -46,21 +47,26 @@ #Check size - if 'size' in options: + if options.exist('size'): fieldsize=options.getfieldvalue('size') - if numpy.isnan(fieldsize[0]): - if not numpy.size(field,1)==fieldsize[1]: + if len(fieldsize) == 1: + if (not numpy.size(field,0)==fieldsize[0]): md = md.checkmessage(options.getfieldvalue('message',\ - "field '%s' should have %d columns" % (fieldname,fieldsize[1]))) - elif numpy.isnan(fieldsize[1]): - if not numpy.size(field,0)==fieldsize[0]: - md = md.checkmessage(options.getfieldvalue('message',\ - "field '%s' should have %d lines" % (fieldname,fieldsize[0]))) - else: - if (not numpy.size(field,0)==fieldsize[0]) or (not numpy.size(field,1)==fieldsize[1]): - md = md.checkmessage(options.getfieldvalue('message',\ - "field '%s' size should be %d x %d" % (fieldname,fieldsize[0],fieldsize[1]))) + "field '%s' size should be %d" % (fieldname,fieldsize[0]))) + elif len(fieldsize) == 2: + if numpy.isnan(fieldsize[0]): + if not numpy.size(field,1)==fieldsize[1]: + md = md.checkmessage(options.getfieldvalue('message',\ + "field '%s' should have %d columns" % (fieldname,fieldsize[1]))) + elif numpy.isnan(fieldsize[1]): + if not numpy.size(field,0)==fieldsize[0]: + md = md.checkmessage(options.getfieldvalue('message',\ + "field '%s' should have %d lines" % (fieldname,fieldsize[0]))) + else: + if (not numpy.size(field,0)==fieldsize[0]) or (not numpy.size(field,1)==fieldsize[1]): + md = md.checkmessage(options.getfieldvalue('message',\ + "field '%s' size should be %d x %d" % (fieldname,fieldsize[0],fieldsize[1]))) #Check numel - if 'numel' in options: + if options.exist('numel'): fieldnumel=options.getfieldvalue('numel') if not numpy.size(field) in fieldnumel: @@ -88,5 +94,5 @@ #check values - if 'values' in options: + if options.exist('values'): fieldvalues=options.getfieldvalue('values') if False in ismember(field,fieldvalues): @@ -102,10 +108,10 @@ #check greater - if '>=' in options: + if options.exist('>='): lowerbound=options.getfieldvalue('>=') if numpy.any(field' in options: + if options.exist('>'): lowerbound=options.getfieldvalue('>') if numpy.any(field<=lowerbound): @@ -114,10 +120,10 @@ #check smaller - if '<=' in options: + if options.exist('<='): upperbound=options.getfieldvalue('<=') if numpy.any(field>upperbound): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have values below %d" % (fieldname,upperbound))) - if '<' in options: + if options.exist('<'): upperbound=options.getfieldvalue('<') if numpy.any(field>=upperbound):