Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17366) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17367) @@ -31,4 +31,5 @@ parameters->AddObject(iomodel->CopyConstantObject(DamageStressThresholdEnum)); parameters->AddObject(iomodel->CopyConstantObject(DamageHealingEnum)); + parameters->AddObject(iomodel->CopyConstantObject(DamageEquivStressEnum)); } xDelete(law); @@ -395,6 +396,7 @@ IssmDouble c1,c2,c3,healing,stress_threshold; IssmDouble s_xx,s_xy,s_yy; - IssmDouble J2s,Xis,Psi,PosPsi,NegPsi; + IssmDouble J2s,Chi,Psi,PosPsi,NegPsi; IssmDouble damage,sigma_xx,sigma_xy,sigma_yy; + int equivstress; /*Fetch number of vertices and allocate output*/ @@ -418,5 +420,8 @@ Input* damage_input = element->GetInput(DamageDbarEnum); _assert_(damage_input); - /*Damage evolution z mapping: */ + /*retrieve the desired type of equivalent stress*/ + element->FindParam(&equivstress,DamageEquivStressEnum); + + /*Calculate damage evolution source term: */ Gauss* gauss=element->NewGauss(); for (int iv=0;iv=',0); md = checkfield(md,'fieldname','damage.stress_threshold','>=',0); + md = checkfield(md,'fieldname','damage.equiv_stress','numel',[1],'values',[0]); elseif strcmpi(obj.law,'undamaged'), if (solution==DamageEvolutionSolutionEnum), @@ -128,4 +131,5 @@ fielddisplay(obj,'healing','damage healing parameter 1'); fielddisplay(obj,'stress_threshold','damage stress threshold [Pa]'); + fielddisplay(obj,'equiv_stress','0: von Mises'); end @@ -151,4 +155,5 @@ WriteData(fid,'object',obj,'fieldname','stress_threshold','format','Double'); WriteData(fid,'object',obj,'fieldname','healing','format','Double'); + WriteData(fid,'object',obj,'fieldname','equiv_stress','format','Integer'); end Index: /issm/trunk-jpl/src/m/classes/damage.py =================================================================== --- /issm/trunk-jpl/src/m/classes/damage.py (revision 17366) +++ /issm/trunk-jpl/src/m/classes/damage.py (revision 17367) @@ -35,4 +35,5 @@ self.c4 = float('NaN') self.healing = float('NaN') + self.equiv_stress = float('NaN') if not len(args): @@ -62,4 +63,5 @@ s+="%s\n" % fielddisplay(self,"c4","damage parameter 4 ") s+="%s\n" % fielddisplay(self,"stress_threshold","damage stress threshold [Pa]") + s+="%s\n" % fielddisplay(self,"equiv_stresss","0: von Mises") return s @@ -95,4 +97,5 @@ self.c4=0 self.healing=0 + self.equiv_stress=0 # }}} @@ -104,5 +107,5 @@ md = checkfield(md,'fieldname','damage.spcdamage','forcing',1) - md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0,1,2]); + md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0,1,2]) md = checkfield(md,'fieldname','damage.maxiter','>=0',0); md = checkfield(md,'fieldname','damage.penalty_factor','>=0',0); @@ -111,5 +114,5 @@ if self.law == 'pralong': - md = checkfield(md,'fieldname','damage.healing','>=',0); + md = checkfield(md,'fieldname','damage.healing','>=',0) md = checkfield(md,'fieldname','damage.c1','>=',0) md = checkfield(md,'fieldname','damage.c2','>=',0) @@ -117,4 +120,5 @@ md = checkfield(md,'fieldname','damage.c4','>=',0) md = checkfield(md,'fieldname','damage.stress_threshold','>=',0) + md = checkfield(md,'fieldname','damage.equiv_stress','numel',[1],'values',[0]) elif strcmpi(self.law,'undamaged'): if (solution==DamageEvolutionSolutionEnum): @@ -142,3 +146,4 @@ WriteData(fid,'object',self,'fieldname','c4','format','Double') WriteData(fid,'object',self,'fieldname','stress_threshold','format','Double') + WriteData(fid,'object',self,'fieldname','equiv_stress','format','Integer') # }}} Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 17366) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 17367) @@ -1107,5 +1107,5 @@ disp(sprintf('%19s: %-22s -- %s','basalforcings' ,['[1x1 ' class(obj.basalforcings) ']'],'bed forcings')); disp(sprintf('%19s: %-22s -- %s','materials' ,['[1x1 ' class(obj.materials) ']'],'material properties')); - disp(sprintf('%19s: %-22s -- %s','damage' ,['[1x1 ' class(obj.damage) ']'],'damage propagation laws')); + disp(sprintf('%19s: %-22s -- %s','damage' ,['[1x1 ' class(obj.damage) ']'],'parameters for damage evolution solution')); disp(sprintf('%19s: %-22s -- %s','friction' ,['[1x1 ' class(obj.friction) ']'],'basal friction/drag properties')); disp(sprintf('%19s: %-22s -- %s','flowequation' ,['[1x1 ' class(obj.flowequation) ']'],'flow equations')); Index: /issm/trunk-jpl/src/m/enum/DamageEquivStressEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/DamageEquivStressEnum.m (revision 17367) +++ /issm/trunk-jpl/src/m/enum/DamageEquivStressEnum.m (revision 17367) @@ -0,0 +1,11 @@ +function macro=DamageEquivStressEnum() +%DAMAGEEQUIVSTRESSENUM - Enum of DamageEquivStress +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=DamageEquivStressEnum() + +macro=StringToEnum('DamageEquivStress'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17366) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17367) @@ -179,4 +179,5 @@ def DamageSpcdamageEnum(): return StringToEnum("DamageSpcdamage")[0] def DamageMaxDamageEnum(): return StringToEnum("DamageMaxDamage")[0] +def DamageEquivStressEnum(): return StringToEnum("DamageEquivStress")[0] def MaterialsRhoIceEnum(): return StringToEnum("MaterialsRhoIce")[0] def MaterialsRhoWaterEnum(): return StringToEnum("MaterialsRhoWater")[0] Index: /issm/trunk-jpl/src/m/mech/mechanicalproperties.m =================================================================== --- /issm/trunk-jpl/src/m/mech/mechanicalproperties.m (revision 17366) +++ /issm/trunk-jpl/src/m/mech/mechanicalproperties.m (revision 17367) @@ -15,5 +15,5 @@ %some checks if length(vx)~=md.mesh.numberofvertices | length(vy)~=md.mesh.numberofvertices, - error(['the input velocity should be of size ' num2str(md.mesh.numberofvertices) '!']) + %error(['the input velocity should be of size ' num2str(md.mesh.numberofvertices) '!']) end if ~(md.mesh.dimension==2) @@ -139,2 +139,6 @@ deviatoricstress.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2); md.results.deviatoricstress=deviatoricstress; + +viscosity=struct('nu',[]); +viscosity.nu=nu; +md.results.viscosity=viscosity; Index: /issm/trunk-jpl/src/m/mech/steadystateiceshelftemp.m =================================================================== --- /issm/trunk-jpl/src/m/mech/steadystateiceshelftemp.m (revision 17366) +++ /issm/trunk-jpl/src/m/mech/steadystateiceshelftemp.m (revision 17367) @@ -46,6 +46,6 @@ %calculate depth-averaged temperature (in Celsius) -%temperature(pos)=-( (Tb(pos)-Ts(pos))*ki./wi(pos) + Hi(pos).*Tb(pos) - (Hi(pos).*Ts(pos) + (Tb(pos)-Ts(pos))*ki./wi(pos)).*exp(Hi(pos).*wi(pos)/ki) )./( Hi(pos).*(exp(Hi(pos).*wi(pos)/ki)-1)); -temperature(pos)=-( ((Tb(pos)-Ts(pos))*ki./wi(pos) + Hi(pos).*Tb(pos))./exp(Hi(pos).*wi(pos)/ki) - Hi(pos).*Ts(pos) + (Tb(pos)-Ts(pos))*ki./wi(pos))./( Hi(pos).*(1-exp(-Hi(pos).*wi(pos)/ki))); +temperature(pos)=-( (Tb(pos)-Ts(pos))*ki./wi(pos) + Hi(pos).*Tb(pos) - (Hi(pos).*Ts(pos) + (Tb(pos)-Ts(pos))*ki./wi(pos)).*exp(Hi(pos).*wi(pos)/ki) )./( Hi(pos).*(exp(Hi(pos).*wi(pos)/ki)-1)); +%temperature(pos)=-( ((Tb(pos)-Ts(pos))*ki./wi(pos) + Hi(pos).*Tb(pos))./exp(Hi(pos).*wi(pos)/ki) - Hi(pos).*Ts(pos) + (Tb(pos)-Ts(pos))*ki./wi(pos))./( Hi(pos).*(1-exp(-Hi(pos).*wi(pos)/ki))); %temperature should not be less than surface temp Index: /issm/trunk-jpl/src/m/mech/thomasparams.m =================================================================== --- /issm/trunk-jpl/src/m/mech/thomasparams.m (revision 17366) +++ /issm/trunk-jpl/src/m/mech/thomasparams.m (revision 17367) @@ -1,3 +1,3 @@ -function [alpha,beta,theta,ex]=thomasparams(md,varargin) +function [alpha,beta,theta,ex,sigxx]=thomasparams(md,varargin) %THOMASPARAMS - compute Thomas' geometric parameters for an ice shelf % @@ -35,10 +35,11 @@ % the equivalent stress % -% 'exx' is the longitudinal strain rate along a coordinate system defined -% by 'coordsys' +% 'exx' is the strain rate along a coordinate system defined by 'coordsys' % -% Usage: [alpha,beta,theta,exx]=ThomasParams(md,options) +% 'sigxx' is the deviatoric stress along a coordinate system defined by 'coordsys' % -% Example: [alpha,beta,theta,exx]=ThomasParams(md,'eq','Thomas','smoothing',2,'coordsys','longitudinal') +% Usage: [alpha,beta,theta,exx,sigxx]=ThomasParams(md,options) +% +% Example: [alpha,beta,theta,exx,sigxx]=ThomasParams(md,'eq','Thomas','smoothing',2,'coordsys','longitudinal') %some checks @@ -82,4 +83,8 @@ end +% rheology +n=averaging(md,md.materials.rheology_n,0); +B=md.materials.rheology_B; + switch coordsys case 'principal' @@ -98,4 +103,5 @@ id=find(e1<0 & e2<0); a(id)=-a(id); % where both strain rates are compressive, enforce negative alpha + sigxx=(abs(ex)./((1+a+a.^2).^((n-1)/2))).^(1./n).*B; %mask=ismember(1:md.mesh.numberofvertices,id); %plotmodel(md,'data',ex,'mask',mask) @@ -114,7 +120,8 @@ a=ey./ex; b=exy./ex; - pos=find(ex<0 & ey<0); - %length(pos) - a(pos)=-a(pos); + %pos=find(ex<0 & ey<0); + %a(pos)=-a(pos); + %sigxx=(abs(ex)./((1+a+a.^2+b.^2).^((n-1)/2))).^(1./n).*B; + sigxx=abs(ex).^(1./n-1).*ex./((1+a+a.^2+b.^2).^((n-1)./(2*n))).*B; otherwise error('argument passed to "coordsys" not valid'); @@ -128,7 +135,4 @@ end a(pos)=-2+1e-3; - -% rheology -n=averaging(md,md.materials.rheology_n,0); switch eq Index: /issm/trunk-jpl/src/m/plot/colormaps/getcolormap.m =================================================================== --- /issm/trunk-jpl/src/m/plot/colormaps/getcolormap.m (revision 17366) +++ /issm/trunk-jpl/src/m/plot/colormaps/getcolormap.m (revision 17367) @@ -41,5 +41,5 @@ elseif strcmpi(map,'Rignot'), alpha=getfieldvalue(options,'alpha',1); - map = hsv; + map = hsv(128); map = rgb2hsv(map); map(:,2) = max(min( (0.1+map(:,1)).^(1/alpha) ,1),0);