Index: /issm/trunk-jpl/src/m/classes/friction.m =================================================================== --- /issm/trunk-jpl/src/m/classes/friction.m (revision 21738) +++ /issm/trunk-jpl/src/m/classes/friction.m (revision 21739) @@ -9,4 +9,6 @@ p = NaN; q = NaN; + coupling = 0; %Silje + effective_pressure = NaN; %Silje end methods @@ -15,4 +17,13 @@ self.p=project3d(md,'vector',self.p,'type','element'); self.q=project3d(md,'vector',self.q,'type','element'); + switch self.coupling + case 0 + case 1 + self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1); + case 2 + error('not implemented yet'); + otherwise + error('not supported yet'); + end end % }}} function self = friction(varargin) % {{{ @@ -36,4 +47,14 @@ md = checkfield(md,'fieldname','friction.q','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); md = checkfield(md,'fieldname','friction.p','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); + md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0 1 2]);%Silje + switch self.coupling + case 0 + case 1 + md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1); + case 2 + error('not implemented yet'); + otherwise + error('not supported yet'); + end end % }}} function disp(self) % {{{ @@ -42,4 +63,6 @@ fielddisplay(self,'p','p exponent'); fielddisplay(self,'q','q exponent'); + fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]');%Silje + fielddisplay(self,'coupling','Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure) and 2 for coupled(not implemented yet)');%Silje end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -50,6 +73,14 @@ 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'); %Silje + switch self.coupling + case 0 + case 1 + WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + case 2 + error('not implemented yet'); + otherwise + error('not supported yet'); + end end % }}} function savemodeljs(self,fid,modelname) % {{{ @@ -58,4 +89,6 @@ writejs1Darray(fid,[modelname '.friction.p'],self.p); writejs1Darray(fid,[modelname '.friction.q'],self.q); + writejs1Darray(fid,[modelname '.friction.coupling'],self.coupling);%Silje + writejs1Darray(fid,[modelname '.friction.effective_pressure'],self.effective_pressure);%Silje end % }}} Index: /issm/trunk-jpl/src/m/classes/friction.py =================================================================== --- /issm/trunk-jpl/src/m/classes/friction.py (revision 21738) +++ /issm/trunk-jpl/src/m/classes/friction.py (revision 21739) @@ -1,2 +1,3 @@ +import numpy as np from fielddisplay import fielddisplay from project3d import project3d @@ -16,5 +17,6 @@ self.p = float('NaN') self.q = float('NaN') - + self.coupling = 0 + self.effective_pressure = float('NaN') #set defaults self.setdefaultparameters() @@ -22,9 +24,11 @@ #}}} 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*bed, r=q/p and s=1/p)" + 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 for default, 1 for forcing(provide md.friction.effective_pressure) and 2 for coupled(not implemented yet)')) + string="%s\n%s"%(string,fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]')) return string #}}} @@ -33,4 +37,11 @@ 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==1: + self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1) + elif self.coupling==2: + raise ValueError('coupling not supported yet') + elif self.coupling > 2: + raise ValueError('md.friction.coupling larger than 2, not supported yet') return self #}}} @@ -47,5 +58,11 @@ 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]) + if self.coupling==1: + md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1) + elif self.coupling==2: + raise ValueError('coupling not supported yet') + elif self.coupling > 2: + raise ValueError('md.friction.coupling larger than 2, not supported yet') return md # }}} @@ -55,3 +72,10 @@ 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==1: + 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==2: + raise ValueError('coupling not supported yet') + elif self.coupling > 2: + raise ValueError('md.friction.coupling larger than 2, not supported yet') # }}} Index: /issm/trunk-jpl/src/m/classes/frictionhydro.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionhydro.m (revision 21738) +++ /issm/trunk-jpl/src/m/classes/frictionhydro.m (revision 21739) @@ -6,5 +6,5 @@ classdef frictionhydro properties (SetAccess=public) - Coupling = 0; + coupling = 0; q = NaN; C = NaN; @@ -28,10 +28,17 @@ %Early return if ~ismember('StressbalanceAnalysis',analyses) & ~ismember('ThermalAnalysis',analyses), return; end - md = checkfield(md,'fieldname','friction.Coupling','numel',[1],'values',[0 1]); + md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0 1 2]); md = checkfield(md,'fieldname','friction.q','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); md = checkfield(md,'fieldname','friction.C','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); md = checkfield(md,'fieldname','friction.As','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); - if self.Coupling==0, - md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1); + switch self.coupling + case 0 + case 1 + md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1); + case 2 + error('not implemented yet'); + otherwise + error('not supported yet'); + end end end % }}} @@ -40,11 +47,17 @@ self.C=project3d(md,'vector',self.C,'type','element'); self.As=project3d(md,'vector',self.As,'type','element'); - if self.Coupling==0, - self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1); + switch self.coupling + case 0 + case 1 + self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1); + case 2 + error('not implemented yet'); + otherwise + error('not supported yet'); end end % }}} function disp(self) % {{{ disp(sprintf('Effective Pressure based friction law described in Gagliardini 2007')); - fielddisplay(self,'Coupling','Coupling flag, 1 for coupling and 0 for forcing'); + fielddisplay(self,'coupling','Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure) and 2 for coupled(not implemented yet)'); fielddisplay(self,'q','friction law exponent q>=1'); fielddisplay(self,'C','friction law max value [SI]'); @@ -54,10 +67,16 @@ function marshall(self,prefix,md,fid) % {{{ WriteData(fid,prefix,'name','md.friction.law','data',3,'format','Integer'); - WriteData(fid,prefix,'class','friction','object',self,'fieldname','Coupling','format','Integer'); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer'); WriteData(fid,prefix,'class','friction','object',self,'fieldname','q','format','DoubleMat','mattype',2); WriteData(fid,prefix,'class','friction','object',self,'fieldname','C','format','DoubleMat','mattype',2); WriteData(fid,prefix,'class','friction','object',self,'fieldname','As','format','DoubleMat','mattype',2); - if self.Coupling==0, - WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + switch self.coupling + case 0 + case 1 + WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + case 2 + error('not implemented yet'); + otherwise + error('not supported yet'); end end % }}} Index: /issm/trunk-jpl/src/m/classes/frictionhydro.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionhydro.py (revision 21739) +++ /issm/trunk-jpl/src/m/classes/frictionhydro.py (revision 21739) @@ -0,0 +1,93 @@ +import numpy as np +from project3d import project3d +from fielddisplay import fielddisplay +from checkfield import checkfield +from WriteData import WriteData + + +class frictionhydro(object): + """ + FRICTIONWEERTMAN class definition + + Usage: + friction=frictionhydro(); + """ + def __init__(self): # {{{ + self.coupling = 0 + self.q = float('NaN') + self.C = float('NaN') + self.As = float('NaN') + self.effective_pressure = float('NaN') + #set defaults + self.setdefaultparameters() + #}}} + + #what is this: + #function self = frictionhydro(varargin) % {{{ + # switch nargin + # case 0 + # self=setdefaultparameters(self); + # otherwise + # error('constructor not supported'); + # end + #end % }}} + + def checkconsistency(self,md,solution,analyses): #{{{ + + #Early return + if 'StressbalanceAnalysis' in analyses and 'ThermalAnalysis' in analyses: + return md + + md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0,1,2]) + md = checkfield(md,'fieldname','friction.q','NaN',1,'Inf',1,'size',[md.mesh.numberofelements,1]) + md = checkfield(md,'fieldname','friction.C','NaN',1,'Inf',1,'size',[md.mesh.numberofelements,1]) + md = checkfield(md,'fieldname','friction.As','NaN',1,'Inf',1,'size',[md.mesh.numberofelements,1]) + if self.coupling==1: + md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1) + elif self.coupling==2: + raise ValueError('coupling not supported yet') + elif self.coupling > 2: + raise ValueError('md.friction.coupling larger than 2, not supported yet') + # }}} + + def extrude(self,md): # {{{ + self.q=project3d(md,'vector',self.q,'type','element') + self.C=project3d(md,'vector',self.C,'type','element') + self.As=project3d(md,'vector',self.As,'type','element') + if self.coupling==1: + self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1) + elif self.coupling==2: + raise ValueError('coupling not supported yet') + elif self.coupling > 2: + raise ValueError('md.friction.coupling larger than 2, not supported yet') + return self + # }}} + + + def __repr__(self): # {{{ + string='Effective Pressure based friction law described in Gagliardini 2007' + string="%s\n%s"%(string,fielddisplay(self,'coupling','Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure) and 2 for coupled(not implemented yet)')) + string="%s\n%s"%(string,fielddisplay(self,'q','friction law exponent q>=1')) + string="%s\n%s"%(string,fielddisplay(self,'C','friction law max value [SI]')) + string="%s\n%s"%(string,fielddisplay(self,'As','Sliding Parameter without cavitation [m Pa^-n s^-1]')) + string="%s\n%s"%(string,fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]')) + + return string + #}}} + + + def marshall(self,prefix,md,fid): #{{{ + WriteData(fid,prefix,'name','md.friction.law','data',3,'format','Integer') + WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer') + WriteData(fid,prefix,'class','friction','object',self,'fieldname','q','format','DoubleMat','mattype',2) + WriteData(fid,prefix,'class','friction','object',self,'fieldname','C','format','DoubleMat','mattype',2) + WriteData(fid,prefix,'class','friction','object',self,'fieldname','As','format','DoubleMat','mattype',2) + if self.coupling==1: + 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==2: + raise ValueError('coupling not supported yet') + elif self.coupling > 2: + raise ValueError('md.friction.coupling larger than 2, not supported yet') + #}}} + +