Index: /issm/trunk-jpl/src/m/classes/friction.py =================================================================== --- /issm/trunk-jpl/src/m/classes/friction.py (revision 23869) +++ /issm/trunk-jpl/src/m/classes/friction.py (revision 23870) @@ -1,3 +1,2 @@ -import numpy as np from fielddisplay import fielddisplay from project3d import project3d @@ -5,71 +4,76 @@ from WriteData import WriteData + class friction(object): - """ - FRICTION class definition + """ + FRICTION class definition - Usage: - friction=friction() - """ + Usage: + friction=friction() + """ - def __init__(self): # {{{ - self.coefficient = float('NaN') - self.p = float('NaN') - self.q = float('NaN') - self.coupling = 0 - self.effective_pressure = float('NaN') - #set defaults - self.setdefaultparameters() + def __init__(self): # {{{ + self.coefficient = float('NaN') + self.p = float('NaN') + self.q = float('NaN') + self.coupling = 0 + self.effective_pressure = float('NaN') + #set defaults + self.setdefaultparameters() + #}}} - #}}} - def __repr__(self): # {{{ - string="Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b,\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*base, r=q/p and s=1/p)" + def __repr__(self): # {{{ + string = "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b,\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*base, r=q/p and s=1/p)" - string="%s\n%s"%(string,fielddisplay(self,"coefficient","friction coefficient [SI]")) - string="%s\n%s"%(string,fielddisplay(self,"p","p exponent")) - string="%s\n%s"%(string,fielddisplay(self,"q","q exponent")) - string="%s\n%s"%(string,fielddisplay(self,'coupling','Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)')) - string="%s\n%s"%(string,fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]')) - return string - #}}} - def extrude(self,md): # {{{ - self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1) - self.p=project3d(md,'vector',self.p,'type','element') - self.q=project3d(md,'vector',self.q,'type','element') - #if self.coupling==0: #doesnt work with empty loop, so just skip it? - if self.coupling in[3,4]: - self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1) - elif self.coupling > 4: - raise ValueError('md.friction.coupling larger than 4, not supported yet') - return self - #}}} - def setdefaultparameters(self): # {{{ - return self - #}}} - def checkconsistency(self,md,solution,analyses): # {{{ + string = "%s\n%s" % (string, fielddisplay(self, "coefficient", "friction coefficient [SI]")) + string = "%s\n%s" % (string, fielddisplay(self, "p", "p exponent")) + string = "%s\n%s" % (string, fielddisplay(self, "q", "q exponent")) + string = "%s\n%s" % (string, fielddisplay(self, 'coupling', 'Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)')) + string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]')) + return string + #}}} - #Early return - if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses: - return md + def extrude(self, md): # {{{ + self.coefficient = project3d(md, 'vector', self.coefficient, 'type', 'node', 'layer', 1) + self.p = project3d(md, 'vector', self.p, 'type', 'element') + self.q = project3d(md, 'vector', self.q, 'type', 'element') + #if self.coupling==0: #doesnt work with empty loop, so just skip it? + if self.coupling in[3, 4]: + self.effective_pressure = project3d(md, 'vector', self.effective_pressure, 'type', 'node', 'layer', 1) + elif self.coupling > 4: + raise ValueError('md.friction.coupling larger than 4, not supported yet') + return self + #}}} - md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1) - md = checkfield(md,'fieldname','friction.q','NaN',1,'Inf',1,'size',[md.mesh.numberofelements]) - md = checkfield(md,'fieldname','friction.p','NaN',1,'Inf',1,'size',[md.mesh.numberofelements]) - md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0,1,2,3,4]) - if self.coupling==3: - md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1) - elif self.coupling > 4: - raise ValueError('md.friction.coupling larger than 4, not supported yet') - return md - # }}} - def marshall(self,prefix,md,fid): # {{{ - WriteData(fid,prefix,'name','md.friction.law','data',1,'format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1) - WriteData(fid,prefix,'object',self,'fieldname','p','format','DoubleMat','mattype',2) - WriteData(fid,prefix,'object',self,'fieldname','q','format','DoubleMat','mattype',2) - WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer') - if self.coupling in[3,4]: - WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - elif self.coupling > 4: - raise ValueError('md.friction.coupling larger than 4, not supported yet') - # }}} + def setdefaultparameters(self): # {{{ + return self + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + + #Early return + if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses: + return md + + md = checkfield(md, 'fieldname', 'friction.coefficient', 'timeseries', 1, 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'friction.q', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofelements]) + md = checkfield(md, 'fieldname', 'friction.p', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofelements]) + md = checkfield(md, 'fieldname', 'friction.coupling', 'numel', [1], 'values', [0, 1, 2, 3, 4]) + if self.coupling == 3: + md = checkfield(md, 'fieldname', 'friction.effective_pressure', 'NaN', 1, 'Inf', 1, 'timeseries', 1) + elif self.coupling > 4: + raise ValueError('md.friction.coupling larger than 4, not supported yet') + return md + # }}} + + def marshall(self, prefix, md, fid): # {{{ + WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 1, 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'p', 'format', 'DoubleMat', 'mattype', 2) + WriteData(fid, prefix, 'object', self, 'fieldname', 'q', 'format', 'DoubleMat', 'mattype', 2) + WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'coupling', 'format', 'Integer') + if self.coupling in[3, 4]: + WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + elif self.coupling > 4: + raise ValueError('md.friction.coupling larger than 4, not supported yet') + # }}} Index: /issm/trunk-jpl/src/m/classes/hydrologydc.py =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologydc.py (revision 23869) +++ /issm/trunk-jpl/src/m/classes/hydrologydc.py (revision 23870) @@ -5,254 +5,259 @@ from WriteData import WriteData + class hydrologydc(object): - """ - Hydrologydc class definition - - Usage: - hydrologydc=hydrologydc(); - """ - - def __init__(self): # {{{ - self.water_compressibility = 0 - self.isefficientlayer = 0 - self.penalty_factor = 0 - self.penalty_lock = 0 - self.rel_tol = 0 - self.max_iter = 0 - self.steps_per_step = 0 - self.sedimentlimit_flag = 0 - self.sedimentlimit = 0 - self.transfer_flag = 0 - self.unconfined_flag = 0 - self.leakage_factor = 0 - self.basal_moulin_input = np.nan - self.requested_outputs = [] - - self.spcsediment_head = np.nan - self.mask_thawed_node = np.nan - self.sediment_transmitivity = np.nan - self.sediment_compressibility = 0 - self.sediment_porosity = 0 - self.sediment_thickness = 0 - - self.spcepl_head = np.nan - self.mask_eplactive_node = np.nan - self.epl_compressibility = 0 - self.epl_porosity = 0 - self.epl_initial_thickness = 0 - self.epl_colapse_thickness = 0 - self.epl_thick_comp = 0 - self.epl_max_thickness = 0 - self.epl_conductivity = 0 - self.eplflip_lock = 0 - - #set defaults - self.setdefaultparameters() - #}}} - def __repr__(self): # {{{ - string=' hydrology Dual Porous Continuum Equivalent parameters:' - string=' - general parameters' - string="%s\n%s"%(string,fielddisplay(self,'water_compressibility','compressibility of water [Pa^-1]')) - string="%s\n%s"%(string,fielddisplay(self,'isefficientlayer','do we use an efficient drainage system [1: true 0: false]')) - string="%s\n%s"%(string,fielddisplay(self,'penalty_factor','exponent of the value used in the penalisation method [dimensionless]')) - string="%s\n%s"%(string,fielddisplay(self,'penalty_lock','stabilize unstable constraints that keep zigzagging after n iteration (default is 0, no stabilization)')) - string="%s\n%s"%(string,fielddisplay(self,'rel_tol','tolerance of the nonlinear iteration for the transfer between layers [dimensionless]')) - string="%s\n%s"%(string,fielddisplay(self,'max_iter','maximum number of nonlinear iteration')) - string="%s\n%s"%(string,fielddisplay(self,'steps_per_step','number of hydrology steps per time step')) - string="%s\n%s"%(string,fielddisplay(self,'basal_moulin_input','water flux at a given point [m3 s-1]')) - string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) - string="%s\n%s"%(string,fielddisplay(self,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer')) - string="%s\n\t\t%s"%(string,'0: no limit') - string="%s\n\t\t%s"%(string,'1: user defined sedimentlimit') - string="%s\n\t\t%s"%(string,'2: hydrostatic pressure') - string="%s\n\t\t%s"%(string,'3: normal stress') - - if self.sedimentlimit_flag==1: - string="%s\n%s"%(string,fielddisplay(self,'sedimentlimit','user defined upper limit for the inefficient layer [m]')) - - string="%s\n%s"%(string,fielddisplay(self,'transfer_flag','what kind of transfer method is applied between the layers')) - string="%s\n\t\t%s"%(string,'0: no transfer') - string="%s\n\t\t%s"%(string,'1: constant leakage factor: leakage_factor') - - if self.transfer_flag is 1: - string="%s\n%s"%(string,fielddisplay(self,'leakage_factor','user defined leakage factor [m]')) - - string="%s\n%s"%(string,fielddisplay(self,'unconfined_flag','using an unconfined scheme or not (transitory)')) - string="%s\n\t\t%s"%(string,'0: Confined only') - string="%s\n\t\t%s"%(string,'1: Confined-Unconfined') - - string="%s\n%s"%(string,' - for the sediment layer') - string="%s\n%s"%(string,fielddisplay(self,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]')) - string="%s\n%s"%(string,fielddisplay(self,'sediment_compressibility','sediment compressibility [Pa^-1]')) - string="%s\n%s"%(string,fielddisplay(self,'sediment_porosity','sediment [dimensionless]')) - string="%s\n%s"%(string,fielddisplay(self,'sediment_thickness','sediment thickness [m]')) - string="%s\n%s"%(string,fielddisplay(self,'sediment_transmitivity','sediment transmitivity [m^2/s]')) - string="%s\n%s"%(string,fielddisplay(self,'mask_thawed_node','IDS is deactivaed (0) on frozen nodes')) - - if self.isefficientlayer==1: - string="%s\n%s"%(string,' - for the epl layer') - string="%s\n%s"%(string,fielddisplay(self,'spcepl_head','epl water head constraints (NaN means no constraint) [m above MSL]')) - string="%s\n%s"%(string,fielddisplay(self,'mask_eplactive_node','active (1) or not (0) EPL')) - string="%s\n%s"%(string,fielddisplay(self,'epl_compressibility','epl compressibility [Pa^-1]')) - string="%s\n%s"%(string,fielddisplay(self,'epl_porosity','epl [dimensionless]')) - string="%s\n%s"%(string,fielddisplay(self,'epl_max_thickness','epl maximal thickness [m]')) - string="%s\n%s"%(string,fielddisplay(self,'epl_initial_thickness','epl initial thickness [m]')) - string="%s\n%s"%(string,fielddisplay(self,'epl_colapse_thickness','epl colapsing thickness [m]')) - string="%s\n%s"%(string,fielddisplay(self,'epl_thick_comp','epl thickness computation flag')) - string="%s\n%s"%(string,fielddisplay(self,'epl_conductivity','epl conductivity [m^2/s]')) - string="%s\n%s"%(string,fielddisplay(self,'eplflip_lock','lock epl activity to avoid flip-floping (default is 0, no stabilization)')) - return string -#}}} - def extrude(self,md): # {{{ - self.spcsediment_head=project3d(md,'vector',self.spcsediment_head,'type','node','layer',1) - self.sediment_transmitivity=project3d(md,'vector',self.sediment_transmitivity,'type','node','layer',1) - self.basal_moulin_input=project3d(md,'vector',self.basal_moulin_input,'type','node','layer',1) - self.mask_thawed_node=project3d(md,'vector',self.mask_thawed_node,'type','node','layer',1) - if self.isefficientlayer==1 : - self.spcepl_head=project3d(md,'vector',self.spcepl_head,'type','node','layer',1) - self.mask_eplactive_node=project3d(md,'vector',self.mask_eplactive_node,'type','node','layer',1) - return self - #}}} - def setdefaultparameters(self): #{{{ - #Parameters from de Fleurian 2014 - self.water_compressibility = 5.04e-10 - self.isefficientlayer = 1 - self.penalty_factor = 3 - self.penalty_lock = 0 - self.rel_tol = 1.0e-06 - self.max_iter = 100 - self.steps_per_step = 1 - self.sedimentlimit_flag = 0 - self.sedimentlimit = 0 - self.transfer_flag = 0 - self.unconfined_flag = 0 - self.leakage_factor = 10.0 - self.requested_outputs = ['default'] - - self.sediment_compressibility = 1.0e-08 - self.sediment_porosity = 0.4 - self.sediment_thickness = 20.0 - self.sediment_transmitivity = 8.0e-04 - - self.epl_compressibility = 1.0e-08 - self.epl_conductivity = 8.0e-02 - self.epl_porosity = 0.4 - self.epl_initial_thickness = 1.0 - self.epl_colapse_thickness = self.sediment_transmitivity/self.epl_conductivity - self.epl_thick_comp = 1 - self.epl_max_thickness = 5.0 - self.eplflip_lock = 0 - - return self - # }}} - - def defaultoutputs(self,md): # {{{ - list = ['SedimentHeadHydrostep','SedimentHeadResidual','EffectivePressureHydrostep'] - if self.isefficientlayer==1: - list.extend(['EplHeadHydrostep','HydrologydcMaskEplactiveNode','HydrologydcMaskEplactiveElt','EplHeadSlopeX','EplHeadSlopeY','HydrologydcEplThicknessHydrostep']) - if self.steps_per_step>1: - list.extend(['EffectivePressure','SedimentHead']) - if self.isefficientlayer==1: - list.extend(['EplHead','HydrologydcEplThickness']) - return list - #}}} - - def initialize(self,md): # {{{ - self.epl_colapse_thickness = self.sediment_transmitivity/self.epl_conductivity - if np.all(np.isnan(self.basal_moulin_input)): - self.basal_moulin_input=np.zeros((md.mesh.numberofvertices)) - print(" no hydrology.basal_moulin_input specified: values set as zero") - - return self - # }}} - def checkconsistency(self,md,solution,analyses): #{{{ - - #Early return - if 'HydrologyDCInefficientAnalysis' not in analyses and 'HydrologyDCEfficientAnalysis' not in analyses: - return md - - md = checkfield(md,'fieldname','hydrology.water_compressibility','numel',[1],'>',0.) - md = checkfield(md,'fieldname','hydrology.isefficientlayer','numel',[1],'values',[0,1]) - md = checkfield(md,'fieldname','hydrology.penalty_factor','>',0.,'numel',[1]) - md = checkfield(md,'fieldname','hydrology.penalty_lock','>=',0.,'numel',[1]) - md = checkfield(md,'fieldname','hydrology.rel_tol','>',0.,'numel',[1]) - md = checkfield(md,'fieldname','hydrology.max_iter','>',0.,'numel',[1]) - md = checkfield(md,'fieldname','hydrology.steps_per_step','>',0.,'numel',[1]) - md = checkfield(md,'fieldname','hydrology.sedimentlimit_flag','numel',[1],'values',[0,1,2,3]) - md = checkfield(md,'fieldname','hydrology.transfer_flag','numel',[1],'values',[0,1]) - md = checkfield(md,'fieldname','hydrology.unconfined_flag','numel',[1],'values',[0,1]) - md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1) - - if self.sedimentlimit_flag==1: - md = checkfield(md,'fieldname','hydrology.sedimentlimit','>',0.,'numel',[1]) - - if self.transfer_flag==1: - md = checkfield(md,'fieldname','hydrology.leakage_factor','>',0.,'numel',[1]) - - md = checkfield(md,'fieldname','hydrology.basal_moulin_input','NaN',1,'Inf',1,'timeseries',1) - md = checkfield(md,'fieldname','hydrology.spcsediment_head','Inf',1,'timeseries',1) - md = checkfield(md,'fieldname','hydrology.sediment_compressibility','>',0.,'numel',[1]) - md = checkfield(md,'fieldname','hydrology.sediment_porosity','>',0.,'numel',[1]) - md = checkfield(md,'fieldname','hydrology.sediment_thickness','>',0.,'numel',[1]) - md = checkfield(md,'fieldname','hydrology.sediment_transmitivity','>=',0,'size',[md.mesh.numberofvertices]) - md = checkfield(md,'fieldname','hydrology.mask_thawed_node','size',[md.mesh.numberofvertices],'values',[0,1]) - if self.isefficientlayer==1: - md = checkfield(md,'fieldname','hydrology.spcepl_head','Inf',1,'timeseries',1) - md = checkfield(md,'fieldname','hydrology.mask_eplactive_node','size',[md.mesh.numberofvertices],'values',[0,1]) - md = checkfield(md,'fieldname','hydrology.epl_compressibility','>',0.,'numel',[1]) - md = checkfield(md,'fieldname','hydrology.epl_porosity','>',0.,'numel',[1]) - md = checkfield(md,'fieldname','hydrology.epl_max_thickness','numel',[1],'>',0.) - md = checkfield(md,'fieldname','hydrology.epl_initial_thickness','numel',[1],'>',0.) - md = checkfield(md,'fieldname','hydrology.epl_colapse_thickness','numel',[1],'>',0.) - md = checkfield(md,'fieldname','hydrology.epl_thick_comp','numel',[1],'values',[0,1]) - md = checkfield(md,'fieldname','hydrology.eplflip_lock','>=',0.,'numel',[1]) - if self.epl_colapse_thickness > self.epl_initial_thickness: - md.checkmessage('Colapsing thickness for EPL larger than initial thickness') - md = checkfield(md,'fieldname','hydrology.epl_conductivity','numel',[1],'>',0.) - # }}} - def marshall(self,prefix,md,fid): #{{{ - WriteData(fid,prefix,'name','md.hydrology.model','data',1,'format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','water_compressibility','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','isefficientlayer','format','Boolean') - WriteData(fid,prefix,'object',self,'fieldname','penalty_factor','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','penalty_lock','format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','rel_tol','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','max_iter','format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','sedimentlimit_flag','format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','transfer_flag','format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','unconfined_flag','format','Integer') - if self.sedimentlimit_flag==1: - WriteData(fid,prefix,'object',self,'fieldname','sedimentlimit','format','Double') - - if self.transfer_flag==1: - WriteData(fid,prefix,'object',self,'fieldname','leakage_factor','format','Double') - - WriteData(fid,prefix,'object',self,'fieldname','basal_moulin_input','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - WriteData(fid,prefix,'object',self,'fieldname','spcsediment_head','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - WriteData(fid,prefix,'object',self,'fieldname','sediment_compressibility','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','sediment_porosity','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','sediment_thickness','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','sediment_transmitivity','format','DoubleMat','mattype',1) - WriteData(fid,prefix,'object',self,'fieldname','mask_thawed_node','format','DoubleMat','mattype',1) - - if self.isefficientlayer==1: - WriteData(fid,prefix,'object',self,'fieldname','spcepl_head','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - WriteData(fid,prefix,'object',self,'fieldname','mask_eplactive_node','format','DoubleMat','mattype',1) - WriteData(fid,prefix,'object',self,'fieldname','epl_compressibility','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','epl_porosity','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','epl_max_thickness','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','epl_initial_thickness','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','epl_colapse_thickness','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','epl_thick_comp','format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','epl_conductivity','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','eplflip_lock','format','Integer') - - #process requested outputs - outputs = self.requested_outputs - indices = [i for i, x in enumerate(outputs) if x == 'default'] - if len(indices) > 0: - outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:] - outputs =outputscopy - WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray') - # }}} + """ + Hydrologydc class definition + + Usage: + hydrologydc=hydrologydc(); + """ + + def __init__(self): # {{{ + self.water_compressibility = 0 + self.isefficientlayer = 0 + self.penalty_factor = 0 + self.penalty_lock = 0 + self.rel_tol = 0 + self.max_iter = 0 + self.steps_per_step = 0 + self.sedimentlimit_flag = 0 + self.sedimentlimit = 0 + self.transfer_flag = 0 + self.unconfined_flag = 0 + self.leakage_factor = 0 + self.basal_moulin_input = np.nan + self.requested_outputs = [] + + self.spcsediment_head = np.nan + self.mask_thawed_node = np.nan + self.sediment_transmitivity = np.nan + self.sediment_compressibility = 0 + self.sediment_porosity = 0 + self.sediment_thickness = 0 + + self.spcepl_head = np.nan + self.mask_eplactive_node = np.nan + self.epl_compressibility = 0 + self.epl_porosity = 0 + self.epl_initial_thickness = 0 + self.epl_colapse_thickness = 0 + self.epl_thick_comp = 0 + self.epl_max_thickness = 0 + self.epl_conductivity = 0 + self.eplflip_lock = 0 + + #set defaults + self.setdefaultparameters() + #}}} + + def __repr__(self): # {{{ + string = ' hydrology Dual Porous Continuum Equivalent parameters:' + string = ' - general parameters' + string = "%s\n%s" % (string, fielddisplay(self, 'water_compressibility', 'compressibility of water [Pa^-1]')) + string = "%s\n%s" % (string, fielddisplay(self, 'isefficientlayer', 'do we use an efficient drainage system [1: true 0: false]')) + string = "%s\n%s" % (string, fielddisplay(self, 'penalty_factor', 'exponent of the value used in the penalisation method [dimensionless]')) + string = "%s\n%s" % (string, fielddisplay(self, 'penalty_lock', 'stabilize unstable constraints that keep zigzagging after n iteration (default is 0, no stabilization)')) + string = "%s\n%s" % (string, fielddisplay(self, 'rel_tol', 'tolerance of the nonlinear iteration for the transfer between layers [dimensionless]')) + string = "%s\n%s" % (string, fielddisplay(self, 'max_iter', 'maximum number of nonlinear iteration')) + string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of hydrology steps per time step')) + string = "%s\n%s" % (string, fielddisplay(self, 'basal_moulin_input', 'water flux at a given point [m3 s-1]')) + string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested')) + string = "%s\n%s" % (string, fielddisplay(self, 'sedimentlimit_flag', 'what kind of upper limit is applied for the inefficient layer')) + string = "%s\n\t\t%s" % (string, '0: no limit') + string = "%s\n\t\t%s" % (string, '1: user defined sedimentlimit') + string = "%s\n\t\t%s" % (string, '2: hydrostatic pressure') + string = "%s\n\t\t%s" % (string, '3: normal stress') + + if self.sedimentlimit_flag == 1: + string = "%s\n%s" % (string, fielddisplay(self, 'sedimentlimit', 'user defined upper limit for the inefficient layer [m]')) + + string = "%s\n%s" % (string, fielddisplay(self, 'transfer_flag', 'what kind of transfer method is applied between the layers')) + string = "%s\n\t\t%s" % (string, '0: no transfer') + string = "%s\n\t\t%s" % (string, '1: constant leakage factor: leakage_factor') + + if self.transfer_flag == 1: + string = "%s\n%s" % (string, fielddisplay(self, 'leakage_factor', 'user defined leakage factor [m]')) + + string = "%s\n%s" % (string, fielddisplay(self, 'unconfined_flag', 'using an unconfined scheme or not (transitory)')) + string = "%s\n\t\t%s" % (string, '0: Confined only') + string = "%s\n\t\t%s" % (string, '1: Confined-Unconfined') + + string = "%s\n%s" % (string, ' - for the sediment layer') + string = "%s\n%s" % (string, fielddisplay(self, 'spcsediment_head', 'sediment water head constraints (NaN means no constraint) [m above MSL]')) + string = "%s\n%s" % (string, fielddisplay(self, 'sediment_compressibility', 'sediment compressibility [Pa^-1]')) + string = "%s\n%s" % (string, fielddisplay(self, 'sediment_porosity', 'sediment [dimensionless]')) + string = "%s\n%s" % (string, fielddisplay(self, 'sediment_thickness', 'sediment thickness [m]')) + string = "%s\n%s" % (string, fielddisplay(self, 'sediment_transmitivity', 'sediment transmitivity [m^2/s]')) + string = "%s\n%s" % (string, fielddisplay(self, 'mask_thawed_node', 'IDS is deactivaed (0) on frozen nodes')) + + if self.isefficientlayer == 1: + string = "%s\n%s" % (string, ' - for the epl layer') + string = "%s\n%s" % (string, fielddisplay(self, 'spcepl_head', 'epl water head constraints (NaN means no constraint) [m above MSL]')) + string = "%s\n%s" % (string, fielddisplay(self, 'mask_eplactive_node', 'active (1) or not (0) EPL')) + string = "%s\n%s" % (string, fielddisplay(self, 'epl_compressibility', 'epl compressibility [Pa^-1]')) + string = "%s\n%s" % (string, fielddisplay(self, 'epl_porosity', 'epl [dimensionless]')) + string = "%s\n%s" % (string, fielddisplay(self, 'epl_max_thickness', 'epl maximal thickness [m]')) + string = "%s\n%s" % (string, fielddisplay(self, 'epl_initial_thickness', 'epl initial thickness [m]')) + string = "%s\n%s" % (string, fielddisplay(self, 'epl_colapse_thickness', 'epl colapsing thickness [m]')) + string = "%s\n%s" % (string, fielddisplay(self, 'epl_thick_comp', 'epl thickness computation flag')) + string = "%s\n%s" % (string, fielddisplay(self, 'epl_conductivity', 'epl conductivity [m^2/s]')) + string = "%s\n%s" % (string, fielddisplay(self, 'eplflip_lock', 'lock epl activity to avoid flip-floping (default is 0, no stabilization)')) + return string + + def extrude(self, md): # {{{ + self.spcsediment_head = project3d(md, 'vector', self.spcsediment_head, 'type', 'node', 'layer', 1) + self.sediment_transmitivity = project3d(md, 'vector', self.sediment_transmitivity, 'type', 'node', 'layer', 1) + self.basal_moulin_input = project3d(md, 'vector', self.basal_moulin_input, 'type', 'node', 'layer', 1) + self.mask_thawed_node = project3d(md, 'vector', self.mask_thawed_node, 'type', 'node', 'layer', 1) + if self.isefficientlayer == 1: + self.spcepl_head = project3d(md, 'vector', self.spcepl_head, 'type', 'node', 'layer', 1) + self.mask_eplactive_node = project3d(md, 'vector', self.mask_eplactive_node, 'type', 'node', 'layer', 1) + return self + #}}} + + def setdefaultparameters(self): #{{{ + #Parameters from de Fleurian 2014 + self.water_compressibility = 5.04e-10 + self.isefficientlayer = 1 + self.penalty_factor = 3 + self.penalty_lock = 0 + self.rel_tol = 1.0e-06 + self.max_iter = 100 + self.steps_per_step = 1 + self.sedimentlimit_flag = 0 + self.sedimentlimit = 0 + self.transfer_flag = 0 + self.unconfined_flag = 0 + self.leakage_factor = 10.0 + self.requested_outputs = ['default'] + + self.sediment_compressibility = 1.0e-08 + self.sediment_porosity = 0.4 + self.sediment_thickness = 20.0 + self.sediment_transmitivity = 8.0e-04 + + self.epl_compressibility = 1.0e-08 + self.epl_conductivity = 8.0e-02 + self.epl_porosity = 0.4 + self.epl_initial_thickness = 1.0 + self.epl_colapse_thickness = self.sediment_transmitivity / self.epl_conductivity + self.epl_thick_comp = 1 + self.epl_max_thickness = 5.0 + self.eplflip_lock = 0 + + return self + # }}} + + def defaultoutputs(self, md): # {{{ + list = ['SedimentHeadHydrostep', 'SedimentHeadResidual', 'EffectivePressureHydrostep'] + if self.isefficientlayer == 1: + list.extend(['EplHeadHydrostep', 'HydrologydcMaskEplactiveNode', 'HydrologydcMaskEplactiveElt', 'EplHeadSlopeX', 'EplHeadSlopeY', 'HydrologydcEplThicknessHydrostep']) + if self.steps_per_step > 1: + list.extend(['EffectivePressure', 'SedimentHead']) + if self.isefficientlayer == 1: + list.extend(['EplHead', 'HydrologydcEplThickness']) + return list + #}}} + + def initialize(self, md): # {{{ + self.epl_colapse_thickness = self.sediment_transmitivity / self.epl_conductivity + if np.all(np.isnan(self.basal_moulin_input)): + self.basal_moulin_input = np.zeros((md.mesh.numberofvertices)) + print(" no hydrology.basal_moulin_input specified: values set as zero") + + return self + # }}} + + def checkconsistency(self, md, solution, analyses): #{{{ + + #Early return + if 'HydrologyDCInefficientAnalysis' not in analyses and 'HydrologyDCEfficientAnalysis' not in analyses: + return md + + md = checkfield(md, 'fieldname', 'hydrology.water_compressibility', 'numel', [1], '>', 0.) + md = checkfield(md, 'fieldname', 'hydrology.isefficientlayer', 'numel', [1], 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'hydrology.penalty_factor', '>', 0., 'numel', [1]) + md = checkfield(md, 'fieldname', 'hydrology.penalty_lock', '>=', 0., 'numel', [1]) + md = checkfield(md, 'fieldname', 'hydrology.rel_tol', '>', 0., 'numel', [1]) + md = checkfield(md, 'fieldname', 'hydrology.max_iter', '>', 0., 'numel', [1]) + md = checkfield(md, 'fieldname', 'hydrology.steps_per_step', '>', 0., 'numel', [1]) + md = checkfield(md, 'fieldname', 'hydrology.sedimentlimit_flag', 'numel', [1], 'values', [0, 1, 2, 3]) + md = checkfield(md, 'fieldname', 'hydrology.transfer_flag', 'numel', [1], 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'hydrology.unconfined_flag', 'numel', [1], 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'hydrology.requested_outputs', 'stringrow', 1) + + if self.sedimentlimit_flag == 1: + md = checkfield(md, 'fieldname', 'hydrology.sedimentlimit', '>', 0., 'numel', [1]) + + if self.transfer_flag == 1: + md = checkfield(md, 'fieldname', 'hydrology.leakage_factor', '>', 0., 'numel', [1]) + + md = checkfield(md, 'fieldname', 'hydrology.basal_moulin_input', 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'hydrology.spcsediment_head', 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'hydrology.sediment_compressibility', '>', 0., 'numel', [1]) + md = checkfield(md, 'fieldname', 'hydrology.sediment_porosity', '>', 0., 'numel', [1]) + md = checkfield(md, 'fieldname', 'hydrology.sediment_thickness', '>', 0., 'numel', [1]) + md = checkfield(md, 'fieldname', 'hydrology.sediment_transmitivity', '>=', 0, 'size', [md.mesh.numberofvertices]) + md = checkfield(md, 'fieldname', 'hydrology.mask_thawed_node', 'size', [md.mesh.numberofvertices], 'values', [0, 1]) + if self.isefficientlayer == 1: + md = checkfield(md, 'fieldname', 'hydrology.spcepl_head', 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'hydrology.mask_eplactive_node', 'size', [md.mesh.numberofvertices], 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'hydrology.epl_compressibility', '>', 0., 'numel', [1]) + md = checkfield(md, 'fieldname', 'hydrology.epl_porosity', '>', 0., 'numel', [1]) + md = checkfield(md, 'fieldname', 'hydrology.epl_max_thickness', 'numel', [1], '>', 0.) + md = checkfield(md, 'fieldname', 'hydrology.epl_initial_thickness', 'numel', [1], '>', 0.) + md = checkfield(md, 'fieldname', 'hydrology.epl_colapse_thickness', 'numel', [1], '>', 0.) + md = checkfield(md, 'fieldname', 'hydrology.epl_thick_comp', 'numel', [1], 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'hydrology.eplflip_lock', '>=', 0., 'numel', [1]) + if self.epl_colapse_thickness > self.epl_initial_thickness: + md.checkmessage('Colapsing thickness for EPL larger than initial thickness') + md = checkfield(md, 'fieldname', 'hydrology.epl_conductivity', 'numel', [1], '>', 0.) + # }}} + + def marshall(self, prefix, md, fid): #{{{ + WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 1, 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'water_compressibility', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'isefficientlayer', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'penalty_factor', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'penalty_lock', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'rel_tol', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'max_iter', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'sedimentlimit_flag', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'transfer_flag', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'unconfined_flag', 'format', 'Integer') + if self.sedimentlimit_flag == 1: + WriteData(fid, prefix, 'object', self, 'fieldname', 'sedimentlimit', 'format', 'Double') + + if self.transfer_flag == 1: + WriteData(fid, prefix, 'object', self, 'fieldname', 'leakage_factor', 'format', 'Double') + + WriteData(fid, prefix, 'object', self, 'fieldname', 'basal_moulin_input', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'spcsediment_head', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'sediment_compressibility', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'sediment_porosity', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'sediment_thickness', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'sediment_transmitivity', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'mask_thawed_node', 'format', 'DoubleMat', 'mattype', 1) + + if self.isefficientlayer == 1: + WriteData(fid, prefix, 'object', self, 'fieldname', 'spcepl_head', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'mask_eplactive_node', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'epl_compressibility', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'epl_porosity', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'epl_max_thickness', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'epl_initial_thickness', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'epl_colapse_thickness', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'epl_thick_comp', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'epl_conductivity', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'eplflip_lock', 'format', 'Integer') + + #process requested outputs + outputs = self.requested_outputs + indices = [i for i, x in enumerate(outputs) if x == 'default'] + if len(indices) > 0: + outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:] + outputs = outputscopy + WriteData(fid, prefix, 'data', outputs, 'name', 'md.hydrology.requested_outputs', 'format', 'StringArray') + # }}} Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 23869) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 23870) @@ -753,6 +753,7 @@ if md.inversion.max_parameters.size > 1: md.inversion.max_parameters = project2d(md, md.inversion.max_parameters, md.mesh.numberoflayers) - if not np.isnan(md.smb.mass_balance).all(): - md.smb.mass_balance = project2d(md, md.smb.mass_balance, md.mesh.numberoflayers) + if hasattr(md.smb, 'mass_balance'): + if not np.isnan(md.smb.mass_balance).all(): + md.smb.mass_balance = project2d(md, md.smb.mass_balance, md.mesh.numberoflayers) #results Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 23869) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 23870) @@ -1,75 +1,76 @@ -def AnalysisConfiguration(solutiontype): #{{{ - """ - ANALYSISCONFIGURATION - return type of analyses, number of analyses +def AnalysisConfiguration(solutiontype): #{{{ + """ + ANALYSISCONFIGURATION - return type of analyses, number of analyses - Usage: - [analyses]=AnalysisConfiguration(solutiontype); - """ + Usage: + [analyses]=AnalysisConfiguration(solutiontype); + """ - if solutiontype == 'StressbalanceSolution': - analyses=['StressbalanceAnalysis','StressbalanceVerticalAnalysis','StressbalanceSIAAnalysis','L2ProjectionBaseAnalysis'] - elif solutiontype == 'SteadystateSolution': - analyses=['StressbalanceAnalysis','StressbalanceVerticalAnalysis','StressbalanceSIAAnalysis','L2ProjectionBaseAnalysis','ThermalAnalysis','MeltingAnalysis','EnthalpyAnalysis'] - elif solutiontype == 'ThermalSolution': - analyses=['EnthalpyAnalysis','ThermalAnalysis','MeltingAnalysis'] - elif solutiontype == 'MasstransportSolution': - analyses=['MasstransportAnalysis'] - elif solutiontype == 'BalancethicknessSolution': - analyses=['BalancethicknessAnalysis'] - elif solutiontype == 'SurfaceSlopeSolution': - analyses=['L2ProjectionBaseAnalysis'] - elif solutiontype == 'BalancevelocitySolution': - analyses=['BalancevelocityAnalysis'] - elif solutiontype == 'BedSlopeSolution': - analyses=['L2ProjectionBaseAnalysis'] - elif solutiontype == 'GiaSolution': - analyses=['GiaIvinsAnalysis'] - elif solutiontype == 'LoveSolution': - analyses=['LoveAnalysis'] - elif solutiontype == 'TransientSolution': - analyses=['StressbalanceAnalysis','StressbalanceVerticalAnalysis','StressbalanceSIAAnalysis','L2ProjectionBaseAnalysis','ThermalAnalysis','MeltingAnalysis','EnthalpyAnalysis','MasstransportAnalysis'] - elif solutiontype == 'HydrologySolution': - analyses=['L2ProjectionBaseAnalysis','HydrologyShreveAnalysis','HydrologyDCInefficientAnalysis','HydrologyDCEfficientAnalysis'] - elif 'DamageEvolutionSolution': - analyses=['DamageEvolutionAnalysis'] + if solutiontype == 'StressbalanceSolution': + analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis'] + elif solutiontype == 'SteadystateSolution': + analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis'] + elif solutiontype == 'ThermalSolution': + analyses = ['EnthalpyAnalysis', 'ThermalAnalysis', 'MeltingAnalysis'] + elif solutiontype == 'MasstransportSolution': + analyses = ['MasstransportAnalysis'] + elif solutiontype == 'BalancethicknessSolution': + analyses = ['BalancethicknessAnalysis'] + elif solutiontype == 'SurfaceSlopeSolution': + analyses = ['L2ProjectionBaseAnalysis'] + elif solutiontype == 'BalancevelocitySolution': + analyses = ['BalancevelocityAnalysis'] + elif solutiontype == 'BedSlopeSolution': + analyses = ['L2ProjectionBaseAnalysis'] + elif solutiontype == 'GiaSolution': + analyses = ['GiaIvinsAnalysis'] + elif solutiontype == 'LoveSolution': + analyses = ['LoveAnalysis'] + elif solutiontype == 'TransientSolution': + analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis'] + elif solutiontype == 'HydrologySolution': + analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis'] + elif 'DamageEvolutionSolution': + analyses = ['DamageEvolutionAnalysis'] - else: - raise TypeError("solution type: '%s' not supported yet!" % solutiontype) + else: + raise TypeError("solution type: '%s' not supported yet!" % solutiontype) - return analyses + return analyses #}}} + def ismodelselfconsistent(md): - """ - ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem. + """ + ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem. - Usage: - ismodelselfconsistent(md), - """ + Usage: + ismodelselfconsistent(md), + """ - #initialize consistency as true - md.private.isconsistent=True + #initialize consistency as true + md.private.isconsistent = True - #Get solution and associated analyses - solution=md.private.solution - analyses=AnalysisConfiguration(solution) + #Get solution and associated analyses + solution = md.private.solution + analyses = AnalysisConfiguration(solution) - #Go through a model fields, check that it is a class, and call checkconsistency - fields=vars(md) -# for field in fields.iterkeys(): - for field in md.properties(): + #Go through a model fields, check that it is a class, and call checkconsistency + #fields = vars(md) + #for field in fields.iterkeys(): + for field in md.properties(): - #Some properties do not need to be checked - if field in ['results','debug','radaroverlay']: - continue + #Some properties do not need to be checked + if field in ['results', 'debug', 'radaroverlay']: + continue - #Check that current field is an object - if not hasattr(getattr(md,field),'checkconsistency'): - md.checkmessage("field '%s' is not an object." % field) + #Check that current field is an object + if not hasattr(getattr(md, field), 'checkconsistency'): + md.checkmessage("field '%s' is not an object." % field) - #Check consistency of the object - exec("md.{}.checkconsistency(md,solution,analyses)".format(field)) + #Check consistency of the object + exec("md.{}.checkconsistency(md, solution, analyses)".format(field)) - #error message if mode is not consistent - if not md.private.isconsistent: - raise RuntimeError('Model not consistent, see messages above.') + #error message if mode is not consistent + if not md.private.isconsistent: + raise RuntimeError('Model not consistent, see messages above.') Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 23869) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 23870) @@ -9,5 +9,5 @@ if path.exists(filename): print('File {} allready exist'.format(filename)) - newname = eval(input('Give a new name or "delete" to replace: ')) + newname = input('Give a new name or "delete" to replace: ') if newname == 'delete': remove(filename) @@ -27,8 +27,9 @@ dimindex = 1 dimlist = [2, md.mesh.numberofelements, md.mesh.numberofvertices, np.shape(md.mesh.elements)[1]] + print('===Creating dimensions===') for i in range(0, 4): if dimlist[i] not in list(DimDict.keys()): dimindex += 1 - NewDim = NCData.createDimension('DimNum'+str(dimindex), dimlist[i]) + NewDim = NCData.createDimension('DimNum' + str(dimindex), dimlist[i]) DimDict[len(NewDim)] = 'DimNum' + str(dimindex) typelist = [bool, str, str, int, float, complex, @@ -37,4 +38,5 @@ groups = dict.keys(md.__dict__) # get all model classes and create respective groups + print('===Creating and populating groups===') for group in groups: NCgroup = NCData.createGroup(str(group)) @@ -65,17 +67,17 @@ Listgroup = Subgroup.createGroup(str(md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__[naming])) except AttributeError: - Listgroup=Subgroup.createGroup(str(md.__dict__[group].__dict__[field].__class__.__name__)+str(listindex)) - Listgroup.__setattr__('classtype',md.__dict__[group].__dict__[field].__getitem__(listindex).__class__.__name__) + Listgroup = Subgroup.createGroup(str(md.__dict__[group].__dict__[field].__class__.__name__) + str(listindex)) + Listgroup.__setattr__('classtype', md.__dict__[group].__dict__[field].__getitem__(listindex).__class__.__name__) try: - subfields=dict.keys(md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__) + subfields = dict.keys(md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__) except AttributeError: - subfields=dict.keys(md.__dict__[group].__dict__[field].__getitem__(listindex)) + subfields = dict.keys(md.__dict__[group].__dict__[field].__getitem__(listindex)) for subfield in subfields: - if subfield!='outlog': + if subfield != 'outlog': try: - Var=md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__[subfield] + Var = md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__[subfield] except AttributeError: Var = md.__dict__[group].__dict__[field].__getitem__(listindex)[subfield] - DimDict = CreateVar(NCData,Var,subfield,Listgroup,DimDict,md.__dict__[group],field,listindex) + DimDict = CreateVar(NCData, Var, subfield, Listgroup, DimDict, md.__dict__[group], field, listindex) # No subgroup, we directly treat the variable elif type(md.__dict__[group].__dict__[field]) in typelist or field == 'bamg': @@ -90,5 +92,5 @@ NCgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__) Var = md.__dict__[group].__dict__[field].data - DimDict = CreateVar(NCData,Var,field,NCgroup,DimDict) + DimDict = CreateVar(NCData, Var, field, NCgroup, DimDict) else: NCgroup.__setattr__('classtype', str(group)) @@ -186,5 +188,5 @@ except KeyError: index = len(DimDict) + 1 # if the dimension does not exist, increment naming - NewDim = NCData.createDimension('DimNum'+str(index), val_shape) # create dimension + NewDim = NCData.createDimension('DimNum' + str(index), val_shape) # create dimension DimDict[len(NewDim)] = 'DimNum' + str(index) # and update the dimension dictionary output = [str(DimDict[val_shape])] + [DimDict[2]] # now proceed with the shape of the value Index: /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (revision 23869) +++ /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (revision 23870) @@ -1,204 +1,264 @@ import numpy as np -import os -import model -import glob -def exportVTK(filename,model,*args): - ''' - vtk export - function exportVTK(filename,model) - creates a directory with the vtk files for displays in paraview - (only work for triangle and wedges based on their number of nodes) - - Give only the results for nw but could be extended to geometry, mask... - - input: filename destination - (string) - ------------------------------------------------------------------ -model this is md - ------------------------------------------------------------------ - By default only the results are exported, you can add whichever - field you need as a string: - add 'geometry' to export md.geometry - - Basile de Fleurian: - ''' - Dir=os.path.basename(filename) - Path=filename[:-len(Dir)] - - if os.path.exists(filename): - print(('File {} allready exist'.format(filename))) - newname=eval(input('Give a new name or "delete" to replace: ')) - if newname=='delete': - filelist = glob.glob(filename+'/*') - for oldfile in filelist: - os.remove(oldfile) - else: - print(('New file name is {}'.format(newname))) - filename=newname - os.mkdir(filename) - else: - os.mkdir(filename) - - # {{{ get the element related variables - if 'z' in dict.keys(model.mesh.__dict__): - points=np.column_stack((model.mesh.x,model.mesh.y,model.mesh.z)) - dim=3 - else: - points=np.column_stack((model.mesh.x,model.mesh.y,np.zeros(np.shape(model.mesh.x)))) - dim=2 - - num_of_points=np.size(model.mesh.x) - num_of_elt=np.shape(model.mesh.elements)[0] - point_per_elt=np.shape(model.mesh.elements)[1] - # }}} - # {{{ Select the type of element function of the number of nodes per elements - if point_per_elt==3: - celltype=5 #triangles - elif point_per_elt==6: - celltype=13 #wedges - else: - error('Your Element definition is not taken into account \n') - # }}} - # {{{ this is the result structure - res_struct=model.results - if (len(res_struct.__dict__)>0): - #Getting all the solutions of the model - solnames=(dict.keys(res_struct.__dict__)) - num_of_sols=len(solnames) - num_of_timesteps=1 - #%building solutionstructure - for solution in solnames: - #looking for multiple time steps - if (np.size(res_struct.__dict__[solution])>num_of_timesteps): - num_of_timesteps=np.size(res_struct.__dict__[solution]) - num_of_timesteps=int(num_of_timesteps) - else: - num_of_timesteps=1 - # }}} - # {{{ write header and mesh - for step in range(0,num_of_timesteps): - timestep=step - fid=open((filename +'/Timestep.vtk'+str(timestep)+'.vtk'),'w+') - fid.write('# vtk DataFile Version 2.0 \n') - fid.write('Data for run %s \n' % model.miscellaneous.name) - fid.write('ASCII \n') - fid.write('DATASET UNSTRUCTURED_GRID \n') - fid.write('POINTS %d float\n' % num_of_points) - if(dim==3): - for point in points: - fid.write('%f %f %f \n'%(point[0], point[1], point[2])) - elif(dim==2): - for point in points: - fid.write('%f %f %f \n'%(point[0], point[1], point[2])) - - fid.write('CELLS %d %d\n' %(num_of_elt, num_of_elt*(point_per_elt+1))) - - if point_per_elt==3: - for elt in range(0, num_of_elt): - fid.write('3 %d %d %d\n' %(model.mesh.elements[elt,0]-1,model.mesh.elements[elt,1]-1,model.mesh.elements[elt,2]-1)) - elif point_per_elt==6: - for elt in range(0, num_of_elt): - fid.write('6 %d %d %d %d %d %d\n' %(model.mesh.elements[elt,0]-1,model.mesh.elements[elt,1]-1,model.mesh.elements[elt,2]-1,model.mesh.elements[elt,3]-1,model.mesh.elements[elt,4]-1,model.mesh.elements[elt,5]-1)) - else: - print('Number of nodes per element not supported') - - fid.write('CELL_TYPES %d\n' %num_of_elt) - for elt in range(0, num_of_elt): - fid.write('%d\n' %celltype) - - fid.write('POINT_DATA %s \n' %str(num_of_points)) - # }}} - # {{{ loop over the different solution structures - if 'solnames' in locals(): - for sol in solnames: - #dealing with results on different timesteps - if(np.size(res_struct.__dict__[sol])>timestep): - timestep = step - else: - timestep = np.size(res_struct.__dict__[sol]) - - #getting the fields in the solution - if(np.size(res_struct.__dict__[sol])>1): - fieldnames=dict.keys(res_struct.__dict__[sol].__getitem__(timestep).__dict__) - else: - fieldnames=dict.keys(res_struct.__dict__[sol].__dict__) - #check which field is a real result and print - for field in fieldnames: - if(np.size(res_struct.__dict__[sol])>1): - fieldstruct=res_struct.__dict__[sol].__getitem__(timestep).__dict__[field] - else: - fieldstruct=res_struct.__dict__[sol].__dict__[field] - - if ((np.size(fieldstruct))==num_of_points): - fid.write('SCALARS %s float 1 \n' % field) - fid.write('LOOKUP_TABLE default\n') - for node in range(0,num_of_points): - #paraview does not like NaN, replacing - if np.isnan(fieldstruct[node]): - fid.write('%e\n' % -9999.9999) - #also checking for verry small value that mess up - elif (abs(fieldstruct[node])<1.0e-20): - fid.write('%e\n' % 0.0) - else: - fid.write('%e\n' % fieldstruct[node]) - # }}} - # {{{ loop on arguments, if something other than result is asked, do it now - for other in args: - other_struct=model.__dict__[other] - othernames=(dict.keys(other_struct.__dict__)) - for field in othernames: - if np.ndim(other_struct.__dict__[field])==1: - if np.size(other_struct.__dict__[field])==num_of_points: - fid.write('SCALARS %s float 1 \n' % field) - fid.write('LOOKUP_TABLE default\n') - for node in range(0,num_of_points): - #paraview does not like NaN, replacing - if np.isnan(other_struct.__dict__[field][node]): - fid.write('%e\n' % -9999.9999) - #also checking for verry small value that mess up - elif (abs(other_struct.__dict__[field][node])<1.0e-20): - fid.write('%e\n' % 0.0) - else: - fid.write('%e\n' % other_struct.__dict__[field][node]) - elif np.ndim(other_struct.__dict__[field])==2: - #deal with forcings - if np.shape(other_struct.__dict__[field])[0]==num_of_points+1: - current_time=res_struct.__dict__[sol].__getitem__(timestep).__dict__['time']/model.__dict__['constants'].__dict__['yts'] - times=other_struct.__dict__[field][-1,:] - if np.any(times==current_time): - time_loc=np.where(times==current_time) - current_force=other_struct.__dict__[field][:-1,time_loc] - else: - precede_time_loc=np.where(timescurrent_time)[0][0] - time_scaling=(current_time-times[precede_time_loc])/(times[follow_time_loc]-times[precede_time_loc]) - current_force=other_struct.__dict__[field][:-1,precede_time_loc]+(other_struct.__dict__[field][:-1,follow_time_loc]-other_struct.__dict__[field][:-1,precede_time_loc])*time_scaling - fid.write('SCALARS %s float 1 \n' % field) - fid.write('LOOKUP_TABLE default\n') - for node in range(0,num_of_points): - #paraview does not like NaN, replacing - if np.isnan(current_force[node]): - fid.write('%e\n' % -9999.9999) - #also checking for verry small value that mess up - elif (abs(current_force[node])<1.0e-20): - fid.write('%e\n' % 0.0) - else: - fid.write('%e\n' % current_force[node]) - # reloaded variable are generally of dim 2 - elif np.shape(other_struct.__dict__[field])[0]==num_of_points: - # we want only vector - if np.shape(other_struct.__dict__[field])[1]==1: - fid.write('SCALARS %s float 1 \n' % field) - fid.write('LOOKUP_TABLE default\n') - for node in range(0,num_of_points): - #paraview does not like NaN, replacing - print((other_struct.__dict__[field][node])) - if np.isnan(other_struct.__dict__[field][node]): - fid.write('%e\n' % -9999.9999) - #also checking for verry small value that mess up - elif (abs(other_struct.__dict__[field][node])<1.0e-20): - fid.write('%e\n' % 0.0) - else: - fid.write('%e\n' % other_struct.__dict__[field][node]) - # }}} - fid.close(); +from os import path, remove, mkdir +from glob import glob + + +def exportVTK(filename, md, *args, enveloppe=False): + ''' + vtk export + function exportVTK(filename,md) + creates a directory with the vtk files for displays in paraview + (only work for triangle and wedges based on their number of nodes) + + Usage: + exportVTK('DirName',md) + exportVTK('DirName',md,'geometry','mesh') + exportVTK('DirName',md,'geometry','mesh',enveloppe=True) + + DirName is the name of the output directory, each timestep then has it + own file ('Timestep.vtkX.vtk') with X the number of the output step + enveloppe is an option keeping only the enveloppe of the md (it is False by default) + + TODO: - make time easily accessible + - make evolving geometry + + Basile de Fleurian: + ''' + + # File checking and creation {{{ + Dir = path.basename(filename) + Path = filename[:-len(Dir)] + if path.exists(filename): + print(('File {} allready exist'.format(filename))) + newname = input('Give a new name or "delete" to replace: ') + if newname == 'delete': + filelist = glob(filename + '/*') + for oldfile in filelist: + remove(oldfile) + else: + print(('New file name is {}'.format(newname))) + filename = newname + mkdir(filename) + else: + mkdir(filename) + # }}} + + # this is the result structure {{{ + print('Getting accessorie variables') + res_struct = md.results + moving_mesh = False + if(type(res_struct) != list): + #Getting all the solutions of the md + solnames = dict.keys(res_struct.__dict__) + num_of_sols = len(solnames) + num_of_timesteps = 1 + #%building solutionstructure + for solution in solnames: + #looking for multiple time steps + if (np.size(res_struct.__dict__[solution]) > num_of_timesteps): + num_of_timesteps = np.size(res_struct.__dict__[solution]) + num_of_timesteps = int(num_of_timesteps) + if 'Surface' in dict.keys(res_struct.__dict__[solution][0].__dict__): + moving_mesh = True + else: + num_of_timesteps = 1 + # }}} + + # get the element related variables {{{ + print('Now treating the mesh') + #first get the general things + dim = int(md.mesh.domaintype()[0]) + every_nodes = md.mesh.numberofvertices + every_cells = md.mesh.numberofelements + + if np.shape(md.mesh.elements)[1] == 3 or enveloppe: + point_per_elt = 3 + celltype = 5 #triangles + elif np.shape(md.mesh.elements)[1] == 6: + point_per_elt = 6 + celltype = 13 #wedges + else: + raise BadDimension('exportVTK does not support your element type') + + if enveloppe: + if dim == 3: + is_enveloppe = np.logical_or(md.mesh.vertexonbase, md.mesh.vertexonsurface) + enveloppe_index = np.where(is_enveloppe)[0] + convert_index = np.nan * np.ones(np.shape(md.mesh.x)) + convert_index = np.asarray([[i, np.where(enveloppe_index == i)[0][0]] for i, val in enumerate(convert_index) if any(enveloppe_index == i)]) + points = np.column_stack((md.mesh.x[enveloppe_index], + md.mesh.y[enveloppe_index], + md.mesh.z[enveloppe_index])) + + num_of_elt = np.size(np.where(np.isnan(md.mesh.lowerelements))) + np.size(np.where(np.isnan(md.mesh.upperelements))) + connect = md.mesh.elements[np.where(is_enveloppe[md.mesh.elements - 1])].reshape(int(num_of_elt), 3) - 1 + for elt in range(0, num_of_elt): + connect[elt, 0] = convert_index[np.where(convert_index == connect[elt, 0])[0], 1][0] + connect[elt, 1] = convert_index[np.where(convert_index == connect[elt, 1])[0], 1][0] + connect[elt, 2] = convert_index[np.where(convert_index == connect[elt, 2])[0], 1][0] + + num_of_points = np.size(enveloppe_index) + else: + raise BadDimension("exportVTK can't get an enveloppe for dimension {}".format(dim)) + + else: + #we get all the mesh, mainly defining dummies + num_of_elt = every_cells + connect = md.mesh.elements - 1 + enveloppe_index = np.arange(0, np.size(md.mesh.x)) + num_of_points = every_nodes + if dim == 2: + points = np.column_stack((md.mesh.x, md.mesh.y, md.geometry.surface)) + elif dim == 3: + points = np.column_stack((md.mesh.x, md.mesh.y, md.mesh.z)) + else: + raise BadDimension('exportVTK does not support dimension {}'.format(dim)) + # }}} + # write header and mesh {{{ + print('Now starting to write stuff') + for step in range(0, num_of_timesteps): + print('Writing for step {}'.format(step)) + saved_cells = {} + timestep = step + fid = open((filename + '/Timestep.vtk' + str(timestep) + '.vtk'), 'w+') + fid.write('# vtk DataFile Version 3.0 \n') + fid.write('Data for run {} \n'.format(md.miscellaneous.name)) + fid.write('ASCII \n') + fid.write('DATASET UNSTRUCTURED_GRID \n') + fid.write('POINTS {:d} float\n'.format(num_of_points)) + #updating z for mesh evolution + if moving_mesh: + base = np.squeeze(res_struct.__dict__['TransientSolution'][step].__dict__['Base'][enveloppe_index]) + thick_change_ratio = (np.squeeze(res_struct.__dict__['TransientSolution'][step].__dict__['Thickness'][enveloppe_index]) / md.geometry.thickness[enveloppe_index]) + above_bed = points[:, 2] - md.geometry.base[enveloppe_index] + altitude = base + thick_change_ratio * above_bed + else: + altitude = points[:, 2] + for index, point in enumerate(points): + fid.write('{:f} {:f} {:f} \n'.format(point[0], point[1], altitude[index])) + + fid.write('CELLS {:d} {:d}\n'.format(num_of_elt, num_of_elt * (point_per_elt + 1))) + + for elt in range(0, num_of_elt): + if celltype == 5: + fid.write('3 {:d} {:d} {:d}\n'.format(connect[elt, 0], + connect[elt, 1], + connect[elt, 2])) + elif celltype == 13: + fid.write('6 {:d} {:d} {:d} {:d} {:d} {:d}\n'.format(connect[elt, 0], + connect[elt, 1], + connect[elt, 2], + connect[elt, 3], + connect[elt, 4], + connect[elt, 5])) + + fid.write('CELL_TYPES {:d}\n'.format(num_of_elt)) + for elt in range(0, num_of_elt): + fid.write('{:d}\n'.format(celltype)) + + fid.write('POINT_DATA {:s} \n'.format(str(num_of_points))) + # }}} + # loop over the different solution structures{{{ + # first check if there are solutions to grab + if 'solnames' in locals(): + for sol in solnames: + treated_res = [] + #dealing with results on different timesteps + if(np.size(res_struct.__dict__[sol]) > timestep): + timestep = step + else: + timestep = np.size(res_struct.__dict__[sol]) + + #getting the fields in the solution + if(type(res_struct.__dict__[sol]) == list): + spe_res_struct = res_struct.__dict__[sol].__getitem__(timestep) + fieldnames = dict.keys(spe_res_struct.__dict__) + else: + spe_res_struct = res_struct.__dict__[sol] + fieldnames = dict.keys(spe_res_struct.__dict__) + + #Sorting scalars, vectors and tensors + tensors = [field for field in fieldnames if field[-2:] in ['xx', 'yy', 'xy', 'zz', 'xz', 'yz']] + non_tensor = [field for field in fieldnames if field not in tensors] + vectors = [field for field in non_tensor if field[-1] in ['x', 'y', 'z']] + + #check which field is a real result and print + for field in fieldnames: + if field in treated_res: + continue + elif field in vectors: + try: + Vxstruct = np.squeeze(spe_res_struct.__dict__[field[:-1] + 'x']) + Vystruct = np.squeeze(spe_res_struct.__dict__[field[:-1] + 'y']) + if dim == 3: + Vzstruct = np.squeeze(spe_res_struct.__dict__[field[:-1] + 'z']) + treated_res += [field[:-1] + 'x', field[:-1] + 'y', field[:-1] + 'z'] + elif dim == 2: + treated_res += [field[:-1] + 'x', field[:-1] + 'y'] + except KeyError: + fieldnames += field + vectors.remove(field) + + fid.write('VECTORS {} float \n'.format(field[:-1])) + for node in range(0, num_of_points): + Vx = cleanOutliers(Vxstruct[enveloppe_index[node]]) + Vy = cleanOutliers(Vystruct[enveloppe_index[node]]) + if dim == 3: + Vz = cleanOutliers(Vzstruct[enveloppe_index[node]]) + fid.write('{:f} {:f} {:f}\n'.format(Vx, Vy, Vz)) + elif dim == 2: + fid.write('{:f} {:f} {:f}\n'.format(Vx, Vy, 0)) + + elif field in tensors: + print("nothing") + else: + if ((np.size(spe_res_struct.__dict__[field])) == every_nodes): + fid.write('SCALARS {} float 1 \n'.format(field)) + fid.write('LOOKUP_TABLE default\n') + for node in range(0, num_of_points): + outval = cleanOutliers(np.squeeze(spe_res_struct.__dict__[field][enveloppe_index[node]])) + fid.write('{:f}\n'.format(outval)) + elif ((np.size(spe_res_struct.__dict__[field])) == every_cells): + saved_cells[field] = np.squeeze(spe_res_struct.__dict__[field]) + # }}} + # loop on arguments, if something other than result is asked, do it now {{{ + for other in args: + other_struct = md.__dict__[other] + othernames = (dict.keys(other_struct.__dict__)) + for field in othernames: + if (np.size(other_struct.__dict__[field]) == every_nodes): + fid.write('SCALARS {} float 1 \n'.format(field)) + fid.write('LOOKUP_TABLE default\n') + for node in range(0, num_of_points): + outval = cleanOutliers(other_struct.__dict__[field][enveloppe_index[node]]) + fid.write('{:f}\n'.format(outval)) + elif (np.size(other_struct.__dict__[field]) == every_cells): + saved_cells[field] = other_struct.__dict__[field] + # }}} + # Now writting cell variables {{{ + if np.size(list(saved_cells.keys())) > 0: + fid.write('CELL_DATA {:d} \n'.format(num_of_elt)) + for key in list(saved_cells.keys()): + fid.write('SCALARS {} float 1 \n'.format(key)) + fid.write('LOOKUP_TABLE default\n') + for cell in range(0, num_of_elt): + outval = cleanOutliers(saved_cells[key][cell]) + fid.write('{:f}\n'.format(outval)) + # }}} + fid.close() + + +def cleanOutliers(Val): + #paraview does not like NaN, replacing + if np.isnan(Val): + CleanVal = -9999.999 + #also checking for very small value that mess up + elif (abs(Val) < 1.0e-20): + CleanVal = 0.0 + else: + CleanVal = Val + return CleanVal + + +class BadDimension(Exception): + """The required dimension is not supported yet.""" Index: /issm/trunk-jpl/src/m/extrusion/project3d.py =================================================================== --- /issm/trunk-jpl/src/m/extrusion/project3d.py (revision 23869) +++ /issm/trunk-jpl/src/m/extrusion/project3d.py (revision 23870) @@ -3,119 +3,116 @@ def project3d(md,*args): - """ - PROJECT3D - vertically project a vector from 2d mesh + """ + PROJECT3D - vertically project a vector from 2d mesh - vertically project a vector from 2d mesh (split in noncoll and coll areas) into a 3d mesh. - This vector can be a node vector of size (md.mesh.numberofvertices2d,N/A) or an - element vector of size (md.mesh.numberofelements2d,N/A). - arguments: - 'vector': 2d vector - 'type': 'element' or 'node'. - options: - 'layer' a layer number where vector should keep its values. If not specified, all layers adopt the - value of the 2d vector. - 'padding': default to 0 (value adopted by other 3d layers not being projected + vertically project a vector from 2d mesh (split in noncoll and coll areas) into a 3d mesh. + This vector can be a node vector of size (md.mesh.numberofvertices2d,N/A) or an + element vector of size (md.mesh.numberofelements2d,N/A). + arguments: + 'vector': 2d vector + 'type': 'element' or 'node'. + options: + 'layer' a layer number where vector should keep its values. If not specified, all layers adopt the + value of the 2d vector. + 'padding': default to 0 (value adopted by other 3d layers not being projected - Examples: - extruded_vector=project3d(md,'vector',vector2d,'type','node','layer',1,'padding',NaN) - extruded_vector=project3d(md,'vector',vector2d,'type','element','padding',0) - extruded_vector=project3d(md,'vector',vector2d,'type','node') - """ + Examples: + extruded_vector=project3d(md,'vector',vector2d,'type','node','layer',1,'padding',NaN) + extruded_vector=project3d(md,'vector',vector2d,'type','element','padding',0) + extruded_vector=project3d(md,'vector',vector2d,'type','node') + """ - #some regular checks - if not md: - raise TypeError("bad usage") - if md.mesh.domaintype().lower() != '3d': - raise TypeError("input model is not 3d") + #some regular checks + if not md: + raise TypeError("bad usage") + if md.mesh.domaintype().lower() != '3d': + raise TypeError("input model is not 3d") - #retrieve parameters from options. - options = pairoptions(*args) - vector2d = options.getfieldvalue('vector') #mandatory - vectype = options.getfieldvalue('type') #mandatory - layer = options.getfieldvalue('layer',0) #optional (do all layers otherwise) - paddingvalue = options.getfieldvalue('padding',0) #0 by default + #retrieve parameters from options. + options = pairoptions(*args) + vector2d = options.getfieldvalue('vector') #mandatory + vectype = options.getfieldvalue('type') #mandatory + layer = options.getfieldvalue('layer', 0) #optional (do all layers otherwise) + paddingvalue = options.getfieldvalue('padding', 0) #0 by default - vector1d=False - if isinstance(vector2d,np.ndarray) and np.ndim(vector2d)==1: - vector1d=True - vector2d=vector2d.reshape(-1,) + vector1d = False + if isinstance(vector2d, np.ndarray) and np.ndim(vector2d) == 1: + vector1d = True + vector2d = vector2d.reshape(-1,) - if isinstance(vector2d,(bool,int,float)) or np.size(vector2d)==1: - projected_vector=vector2d + if isinstance(vector2d, (bool, int, float)) or np.size(vector2d) == 1: + projected_vector = vector2d - elif vectype.lower()=='node': + elif vectype.lower() == 'node': + #Initialize 3d vector + if np.ndim(vector2d) == 1: + if vector2d.shape[0] == md.mesh.numberofvertices2d: + projected_vector = (paddingvalue * np.ones((md.mesh.numberofvertices))).astype(vector2d.dtype) + elif vector2d.shape[0] == md.mesh.numberofvertices2d + 1: + projected_vector = (paddingvalue * np.ones((md.mesh.numberofvertices + 1))).astype(vector2d.dtype) + projected_vector[-1] = vector2d[-1] + vector2d = vector2d[:-1] + else: + raise TypeError("vector length not supported") + #Fill in + if layer == 0: + for i in range(md.mesh.numberoflayers): + projected_vector[(i * md.mesh.numberofvertices2d):((i + 1) * md.mesh.numberofvertices2d)] = vector2d + else: + projected_vector[((layer - 1) * md.mesh.numberofvertices2d):(layer * md.mesh.numberofvertices2d)] = vector2d + else: + if vector2d.shape[0] == md.mesh.numberofvertices2d: + projected_vector = (paddingvalue * np.ones((md.mesh.numberofvertices, np.size(vector2d, axis=1)))).astype(vector2d.dtype) + elif vector2d.shape[0] == md.mesh.numberofvertices2d + 1: + projected_vector = (paddingvalue * np.ones((md.mesh.numberofvertices + 1, np.size(vector2d, axis=1)))).astype(vector2d.dtype) + projected_vector[-1, :] = vector2d[-1, :] + vector2d = vector2d[:-1, :] + else: + raise TypeError("vector length not supported") + #Fill in + if layer == 0: + for i in range(md.mesh.numberoflayers): + projected_vector[(i * md.mesh.numberofvertices2d):((i + 1) * md.mesh.numberofvertices2d), :] = vector2d + else: + projected_vector[((layer - 1) * md.mesh.numberofvertices2d):(layer * md.mesh.numberofvertices2d), :] = vector2d - #Initialize 3d vector - if np.ndim(vector2d)==1: - if vector2d.shape[0]==md.mesh.numberofvertices2d: - projected_vector=(paddingvalue*np.ones((md.mesh.numberofvertices))).astype(vector2d.dtype) - elif vector2d.shape[0]==md.mesh.numberofvertices2d+1: - projected_vector=(paddingvalue*np.ones((md.mesh.numberofvertices+1))).astype(vector2d.dtype) - projected_vector[-1]=vector2d[-1] - vector2d=vector2d[:-1] - else: - raise TypeError("vector length not supported") - #Fill in - if layer==0: - for i in range(md.mesh.numberoflayers): - projected_vector[(i*md.mesh.numberofvertices2d):((i+1)*md.mesh.numberofvertices2d)]=vector2d - else: - projected_vector[((layer-1)*md.mesh.numberofvertices2d):(layer*md.mesh.numberofvertices2d)]=vector2d - else: - if vector2d.shape[0]==md.mesh.numberofvertices2d: - projected_vector=(paddingvalue*np.ones((md.mesh.numberofvertices,np.size(vector2d,axis=1)))).astype(vector2d.dtype) - elif vector2d.shape[0]==md.mesh.numberofvertices2d+1: - projected_vector=(paddingvalue*np.ones((md.mesh.numberofvertices+1,np.size(vector2d,axis=1)))).astype(vector2d.dtype) - projected_vector[-1,:]=vector2d[-1,:] - vector2d=vector2d[:-1,:] - else: - raise TypeError("vector length not supported") - #Fill in - if layer==0: - for i in range(md.mesh.numberoflayers): - projected_vector[(i*md.mesh.numberofvertices2d):((i+1)*md.mesh.numberofvertices2d),:]=vector2d - else: - projected_vector[((layer-1)*md.mesh.numberofvertices2d):(layer*md.mesh.numberofvertices2d),:]=vector2d + elif vectype.lower() == 'element': + #Initialize 3d vector + if np.ndim(vector2d) == 1: + if vector2d.shape[0] == md.mesh.numberofelements2d: + projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements))).astype(vector2d.dtype) + elif vector2d.shape[0] == md.mesh.numberofelements2d + 1: + projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements + 1))).astype(vector2d.dtype) + projected_vector[-1] = vector2d[-1] + vector2d = vector2d[:-1] + else: + raise TypeError("vector length not supported") + #Fill in + if layer == 0: + for i in range(md.mesh.numberoflayers - 1): + projected_vector[(i * md.mesh.numberofelements2d):((i + 1) * md.mesh.numberofelements2d)] = vector2d + else: + projected_vector[((layer - 1) * md.mesh.numberofelements2d):(layer * md.mesh.numberofelements2d)] = vector2d + else: + if vector2d.shape[0] == md.mesh.numberofelements2d: + projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements, np.size(vector2d, axis=1)))).astype(vector2d.dtype) + elif vector2d.shape[0] == md.mesh.numberofelements2d + 1: + projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements + 1, np.size(vector2d, axis=1)))).astype(vector2d.dtype) + projected_vector[-1, :] = vector2d[-1, :] + vector2d = vector2d[:-1, :] + else: + raise TypeError("vector length not supported") + #Fill in + if layer == 0: + for i in range(md.mesh.numberoflayers - 1): + projected_vector[(i * md.mesh.numberofelements2d):((i + 1) * md.mesh.numberofelements2d), :] = vector2d + else: + projected_vector[((layer - 1) * md.mesh.numberofelements2d):(layer * md.mesh.numberofelements2d), :] = vector2d + else: + raise TypeError("project3d error message: unknown projection type") - elif vectype.lower()=='element': + if vector1d: + projected_vector = projected_vector.reshape(-1,) - #Initialize 3d vector - if np.ndim(vector2d)==1: - if vector2d.shape[0]==md.mesh.numberofelements2d: - projected_vector=(paddingvalue*np.ones((md.mesh.numberofelements))).astype(vector2d.dtype) - elif vector2d.shape[0]==md.mesh.numberofelements2d+1: - projected_vector=(paddingvalue*np.ones((md.mesh.numberofelements+1))).astype(vector2d.dtype) - projected_vector[-1]=vector2d[-1] - vector2d=vector2d[:-1] - else: - raise TypeError("vector length not supported") - #Fill in - if layer==0: - for i in range(md.mesh.numberoflayers-1): - projected_vector[(i*md.mesh.numberofelements2d):((i+1)*md.mesh.numberofelements2d)]=vector2d - else: - projected_vector[((layer-1)*md.mesh.numberofelements2d):(layer*md.mesh.numberofelements2d)]=vector2d - else: - if vector2d.shape[0]==md.mesh.numberofelements2d: - projected_vector=(paddingvalue*np.ones((md.mesh.numberofelements, np.size(vector2d,axis=1)))).astype(vector2d.dtype) - elif vector2d.shape[0]==md.mesh.numberofelements2d+1: - projected_vector=(paddingvalue*np.ones((md.mesh.numberofelements+1,np.size(vector2d,axis=1)))).astype(vector2d.dtype) - projected_vector[-1,:]=vector2d[-1,:] - vector2d=vector2d[:-1,:] - else: - raise TypeError("vector length not supported") - #Fill in - if layer==0: - for i in range(md.mesh.numberoflayers-1): - projected_vector[(i*md.mesh.numberofelements2d):((i+1)*md.mesh.numberofelements2d),:]=vector2d - else: - projected_vector[((layer-1)*md.mesh.numberofelements2d):(layer*md.mesh.numberofelements2d),:]=vector2d - - else: - raise TypeError("project3d error message: unknown projection type") - - if vector1d: - projected_vector=projected_vector.reshape(-1,) - - return projected_vector + return projected_vector Index: /issm/trunk-jpl/src/m/plot/plotmodel.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plotmodel.py (revision 23869) +++ /issm/trunk-jpl/src/m/plot/plotmodel.py (revision 23870) @@ -1,108 +1,105 @@ import numpy as np from plotoptions import plotoptions -from plotdoc import plotdoc from plot_manager import plot_manager from math import ceil, sqrt try: - import pylab as p - import matplotlib.pyplot as plt - from mpl_toolkits.axes_grid1 import ImageGrid, AxesGrid - from mpl_toolkits.mplot3d import Axes3D + import matplotlib.pyplot as plt + from mpl_toolkits.axes_grid1 import ImageGrid except ImportError: - print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled") - -def plotmodel(md,*args): - ''' at command prompt, type 'plotdoc()' for additional documentation - ''' - - #First process options - options=plotoptions(*args) - - #get number of subplots - subplotwidth=ceil(sqrt(options.numberofplots)) - #Get figure number and number of plots - figurenumber=options.figurenumber - numberofplots=options.numberofplots - - #get hold - hold=options.list[0].getfieldvalue('hold',False) - - #if nrows and ncols specified, then bypass - if options.list[0].exist('nrows'): - nrows=options.list[0].getfieldvalue('nrows') - nr=True - else: - nrows=np.ceil(numberofplots/subplotwidth) - nr=False - - if options.list[0].exist('ncols'): - ncols=options.list[0].getfieldvalue('ncols') - nc=True - else: - ncols=int(subplotwidth) - nc=False - ncols=int(ncols) - nrows=int(nrows) - - #check that nrows and ncols were given at the same time! - if not nr==nc: - raise Exception('error: nrows and ncols need to be specified together, or not at all') - - #Go through plots - if numberofplots: - #if plt.fignum_exists(figurenumber): - # plt.cla() - - #if figsize specified - if options.list[0].exist('figsize'): - figsize=options.list[0].getfieldvalue('figsize') - fig=plt.figure(figurenumber,figsize=(figsize[0],figsize[1]))#,tight_layout=True) - else: - fig=plt.figure(figurenumber)#,tight_layout=True) - fig.clf() - - backgroundcolor=options.list[0].getfieldvalue('backgroundcolor',(0.7,0.7,0.7)) - fig.set_facecolor(backgroundcolor) + print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled") - translator={'on':'each', - 'off':'None', - 'one':'single'} - # options needed to define plot grid - plotnum=options.numberofplots - if plotnum==1: - plotnum=None - direction=options.list[0].getfieldvalue('direction','row') # row,column - axes_pad=options.list[0].getfieldvalue('axes_pad',0.25) - add_all=options.list[0].getfieldvalue('add_all',True) # True,False - share_all=options.list[0].getfieldvalue('share_all',True) # True,False - label_mode=options.list[0].getfieldvalue('label_mode','L') # 1,L,all - colorbar=options.list[0].getfieldvalue('colorbar','on') # on, off (single) - cbar_mode=translator[colorbar] - cbar_location=options.list[0].getfieldvalue('colorbarpos','right') # right,top - cbar_size=options.list[0].getfieldvalue('colorbarsize','5%') - cbar_pad=options.list[0].getfieldvalue('colorbarpad',0.025) # None or % +def plotmodel(md, *args): + ''' at command prompt, type 'plotdoc()' for additional documentation + ''' - axgrid=ImageGrid(fig,111, - nrows_ncols=(nrows,ncols), - ngrids=plotnum, - direction=direction, - axes_pad=axes_pad, - add_all=add_all, - share_all=share_all, - label_mode=label_mode, - cbar_mode=cbar_mode, - cbar_location=cbar_location, - cbar_size=cbar_size, - cbar_pad=cbar_pad) + #First process options + options = plotoptions(*args) - if cbar_mode=='None': - for ax in axgrid.cbar_axes: - fig._axstack.remove(ax) + #get number of subplots + subplotwidth = ceil(sqrt(options.numberofplots)) + #Get figure number and number of plots + figurenumber = options.figurenumber + numberofplots = options.numberofplots - for i,ax in enumerate(axgrid.axes_all): - plot_manager(options.list[i].getfieldvalue('model',md),options.list[i],fig,axgrid,i) - fig.show() - else: - raise Exception('plotmodel error message: no output data found.') + #get hold + hold = options.list[0].getfieldvalue('hold', False) + + #if nrows and ncols specified, then bypass + if options.list[0].exist('nrows'): + nrows = options.list[0].getfieldvalue('nrows') + nr = True + else: + nrows = np.ceil(numberofplots / subplotwidth) + nr = False + + if options.list[0].exist('ncols'): + ncols = options.list[0].getfieldvalue('ncols') + nc = True + else: + ncols = int(subplotwidth) + nc = False + ncols = int(ncols) + nrows = int(nrows) + + #check that nrows and ncols were given at the same time! + if not nr == nc: + raise Exception('error: nrows and ncols need to be specified together, or not at all') + + #Go through plots + if numberofplots: + #if plt.fignum_exists(figurenumber): + # plt.cla() + + #if figsize specified + if options.list[0].exist('figsize'): + figsize = options.list[0].getfieldvalue('figsize') + fig = plt.figure(figurenumber, figsize=(figsize[0], figsize[1])) + else: + fig = plt.figure(figurenumber) + fig.clf() + + backgroundcolor = options.list[0].getfieldvalue('backgroundcolor', (0.7, 0.7, 0.7)) + fig.set_facecolor(backgroundcolor) + + translator = {'on': 'each', + 'off': 'None', + 'one': 'single'} + # options needed to define plot grid + plotnum = options.numberofplots + if plotnum == 1: + plotnum = None + direction = options.list[0].getfieldvalue('direction', 'row') # row, column + axes_pad = options.list[0].getfieldvalue('axes_pad', 0.25) + add_all = options.list[0].getfieldvalue('add_all', True) # True, False + share_all = options.list[0].getfieldvalue('share_all', True) # True, False + label_mode = options.list[0].getfieldvalue('label_mode', 'L') # 1, L, all + colorbar = options.list[0].getfieldvalue('colorbar', 'on') # on, off (single) + cbar_mode = translator[colorbar] + cbar_location = options.list[0].getfieldvalue('colorbarpos', 'right') # right, top + cbar_size = options.list[0].getfieldvalue('colorbarsize', '5%') + cbar_pad = options.list[0].getfieldvalue('colorbarpad', 0.025) # None or % + + axgrid = ImageGrid(fig, 111, + nrows_ncols=(nrows, ncols), + #ngrids=plotnum, + direction=direction, + axes_pad=axes_pad, + add_all=add_all, + share_all=share_all, + label_mode=label_mode, + cbar_mode=cbar_mode, + cbar_location=cbar_location, + cbar_size=cbar_size, + cbar_pad=cbar_pad) + + if cbar_mode == 'None': + for ax in axgrid.cbar_axes: + fig._axstack.remove(ax) + + for i, ax in enumerate(axgrid.axes_all): + plot_manager(options.list[i].getfieldvalue('model', md), options.list[i], fig, axgrid, i) + fig.show() + else: + raise Exception('plotmodel error message: no output data found.')