Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17762) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17763) @@ -164,4 +164,5 @@ MaterialsRheologyLawEnum, MaterialsRheologyNEnum, + DamageIsdamageEnum, DamageDEnum, DamageFEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17762) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17763) @@ -172,4 +172,5 @@ case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw"; case MaterialsRheologyNEnum : return "MaterialsRheologyN"; + case DamageIsdamageEnum : return "DamageIsdamage"; case DamageDEnum : return "DamageD"; case DamageFEnum : return "DamageF"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17762) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17763) @@ -175,4 +175,5 @@ else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum; else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum; + else if (strcmp(name,"DamageIsdamage")==0) return DamageIsdamageEnum; else if (strcmp(name,"DamageD")==0) return DamageDEnum; else if (strcmp(name,"DamageF")==0) return DamageFEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum; else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum; - else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum; + if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum; + else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum; else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum; else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum; else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum; - else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; + if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum; + else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; else if (strcmp(name,"FSpressure")==0) return FSpressureEnum; else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum; else if (strcmp(name,"Vel")==0) return VelEnum; - else if (strcmp(name,"Velocity")==0) return VelocityEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"VxAverage")==0) return VxAverageEnum; + if (strcmp(name,"Velocity")==0) return VelocityEnum; + else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; else if (strcmp(name,"Vx")==0) return VxEnum; else if (strcmp(name,"VxPicard")==0) return VxPicardEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; else if (strcmp(name,"SubelementMigration2")==0) return SubelementMigration2Enum; - else if (strcmp(name,"Contact")==0) return ContactEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum; + if (strcmp(name,"Contact")==0) return ContactEnum; + else if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum; else if (strcmp(name,"QmuMaskGroundediceLevelset")==0) return QmuMaskGroundediceLevelsetEnum; else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (revision 17762) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (revision 17763) @@ -32,8 +32,9 @@ else: #Guess where the ice front is - vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices,1)) - pos=numpy.nonzero(numpy.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0.,axis=1) >0.)[0] - vertexonfloatingice[md.mesh.elements[pos].astype(int)-1]=1. - vertexonicefront=numpy.logical_and(numpy.reshape(md.mesh.vertexonboundary,(-1,1)),vertexonfloatingice>0.) + vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices,)) + pos=numpy.nonzero(numpy.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0,axis=1) > 0)[0] + #vertexonfloatingice[md.mesh.elements[pos].astype(int)-1]=1. + vertexonfloatingice[md.mesh.elements[pos]-1]=1. + vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice>0.) # pos=find(md.mesh.vertexonboundary & ~vertexonicefront); @@ -42,7 +43,7 @@ print "SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually." - md.stressbalance.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) - md.stressbalance.spcvy=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) - md.stressbalance.spcvz=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + md.stressbalance.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices,)) + md.stressbalance.spcvy=float('nan')*numpy.ones((md.mesh.numberofvertices,)) + md.stressbalance.spcvz=float('nan')*numpy.ones((md.mesh.numberofvertices,)) md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6)) md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3)) @@ -68,7 +69,7 @@ else: pos=numpy.nonzero(md.mesh.vertexonboundary)[0] - md.stressbalance.spcvx[pos[:]]=0 - md.stressbalance.spcvy[pos[:]]=0 - md.stressbalance.spcvz[pos[:]]=0 + md.stressbalance.spcvx[pos]=0 + md.stressbalance.spcvy[pos]=0 + md.stressbalance.spcvz[pos]=0 #Dirichlet Values @@ -90,22 +91,21 @@ #Deal with other boundary conditions if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)): - md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1)) + md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,)) print " no balancethickness.thickening_rate specified: values set as zero" - md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) - md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) - md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,)) + md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,)) + md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,)) if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: - md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,)) if hasattr(md.mesh,'vertexonsurface'): pos=numpy.nonzero(md.mesh.vertexonsurface)[0] md.thermal.spctemperature[pos]=md.initialization.temperature[pos] #impose observed temperature on surface if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices: - md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1)) - md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)]=50.*10.**-3 #50mW/m2 + md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,)) + md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)[0]]=50.*10.**-3 #50mW/m2 else: print " no thermal boundary conditions created: no observed temperature found" return md - Index: /issm/trunk-jpl/src/m/classes/damage.py =================================================================== --- /issm/trunk-jpl/src/m/classes/damage.py (revision 17762) +++ /issm/trunk-jpl/src/m/classes/damage.py (revision 17763) @@ -17,4 +17,5 @@ #damage: + self.isdamage = float('NaN') self.D = float('NaN') self.law = '' @@ -24,6 +25,7 @@ #numerical stabilization = float('NaN') + maxiter = float('NaN') + elementinterp = '' penalty_threshold = float('NaN') - maxiter = float('NaN') penalty_lock = float('NaN') penalty_factor = float('NaN') @@ -47,15 +49,18 @@ def __repr__(self): # {{{ s =' Damage:\n' + + s+="%s\n" % fielddisplay(self,"isdamage","is damage mechanics being used? [0 (default) or 1]") + if self.isdamage: + s+="%s\n" % fielddisplay(self,"D","damage tensor (scalar for now)") + s+="%s\n" % fielddisplay(self,"law","damage law (string) from ['undamaged','pralong']") + s+="%s\n" % fielddisplay(self,"spcdamage","damage constraints (NaN means no constraint)") + s+="%s\n" % fielddisplay(self,"max_damage","maximum possible damage (0<=max_damage<1)") - s+="%s\n" % fielddisplay(self,"D","damage tensor (scalar for now)") - s+="%s\n" % fielddisplay(self,"law","damage law (string) from ['undamaged','pralong']") - s+="%s\n" % fielddisplay(self,"spcdamage","damage constraints (NaN means no constraint)") - s+="%s\n" % fielddisplay(self,"max_damage","maximum possible damage (0<=max_damage<1)") - - s+="%s\n" % fielddisplay(self,"stabilization","0: no, 1: artificial_diffusivity, 2: SUPG") - s+="%s\n" % fielddisplay(self,"maxiter","maximum number of non linear iterations") - s+="%s\n" % fielddisplay(self,"penalty_lock","stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)") - s+="%s\n" % fielddisplay(self,"penalty_threshold","threshold to declare convergence of damage evolution solution (default is 0)") - s+="%s\n" % fielddisplay(self,"penalty_factor","scaling exponent (default is 3)") + s+="%s\n" % fielddisplay(self,"stabilization","0: no, 1: artificial_diffusivity, 2: SUPG") + s+="%s\n" % fielddisplay(self,"maxiter","maximum number of non linear iterations") + s+="%s\n" % fielddisplay(self,"elementinterp","interpolation scheme for finite elements {''P1'',''P2''}") + s+="%s\n" % fielddisplay(self,"penalty_lock","stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)") + s+="%s\n" % fielddisplay(self,"penalty_threshold","threshold to declare convergence of damage evolution solution (default is 0)") + s+="%s\n" % fielddisplay(self,"penalty_factor","scaling exponent (default is 3)") if (self.law=='pralong'): @@ -74,4 +79,5 @@ #damage parameters: + self.isdamage=0 self.D=0 self.law='undamaged' @@ -84,4 +90,7 @@ #Maximum number of iterations self.maxiter=100 + + #finite element interpolation + self.elementinterp='P1' #factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor @@ -119,14 +128,17 @@ def checkconsistency(self,md,solution,analyses): # {{{ - md = checkfield(md,'fieldname','damage.D','>=',0,'<=',self.max_damage,'size',[md.mesh.numberofvertices]) - md = checkfield(md,'fieldname','damage.max_damage','<',1,'>=',0) - md = checkfield(md,'fieldname','damage.law','values',['undamaged','pralong']) - 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.maxiter','>=0',0) - md = checkfield(md,'fieldname','damage.penalty_factor','>=0',0) - md = checkfield(md,'fieldname','damage.penalty_lock','>=0',0) - md = checkfield(md,'fieldname','damage.penalty_threshold','>=0',0) + md = checkfield(md,'fieldname','damage.isdamage','numel',[1],'values',[0,1]) + if self.isdamage: + md = checkfield(md,'fieldname','damage.D','>=',0,'<=',self.max_damage,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','damage.max_damage','<',1,'>=',0) + md = checkfield(md,'fieldname','damage.law','values',['undamaged','pralong']) + 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.maxiter','>=0',0) + md = checkfield(md,'fieldname','damage.elementinterp','values',['P1','P2']) + md = checkfield(md,'fieldname','damage.penalty_factor','>=0',0) + md = checkfield(md,'fieldname','damage.penalty_lock','>=0',0) + md = checkfield(md,'fieldname','damage.penalty_threshold','>=0',0) if self.law == 'pralong': @@ -148,14 +160,17 @@ def marshall(self,md,fid): # {{{ - WriteData(fid,'object',self,'fieldname','D','format','DoubleMat','mattype',1) - WriteData(fid,'object',self,'fieldname','law','format','String') - WriteData(fid,'object',self,'fieldname','spcdamage','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) - WriteData(fid,'object',self,'fieldname','max_damage','format','Double') + WriteData(fid,'object',self,'fieldname','isdamage','format','Boolean') + if self.isdamage: + WriteData(fid,'object',self,'fieldname','D','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'fieldname','law','format','String') + WriteData(fid,'object',self,'fieldname','spcdamage','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) + WriteData(fid,'object',self,'fieldname','max_damage','format','Double') - WriteData(fid,'object',self,'fieldname','stabilization','format','Integer') - WriteData(fid,'object',self,'fieldname','penalty_threshold','format','Integer') - WriteData(fid,'object',self,'fieldname','maxiter','format','Integer') - WriteData(fid,'object',self,'fieldname','penalty_lock','format','Integer') - WriteData(fid,'object',self,'fieldname','penalty_factor','format','Double') + WriteData(fid,'object',self,'fieldname','stabilization','format','Integer') + WriteData(fid,'object',self,'fieldname','elementinterp','format','String') + WriteData(fid,'object',self,'fieldname','penalty_threshold','format','Integer') + WriteData(fid,'object',self,'fieldname','maxiter','format','Integer') + WriteData(fid,'object',self,'fieldname','penalty_lock','format','Integer') + WriteData(fid,'object',self,'fieldname','penalty_factor','format','Double') if self.law=='pralong': Index: /issm/trunk-jpl/src/m/enum/DamageIsdamageEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/DamageIsdamageEnum.m (revision 17763) +++ /issm/trunk-jpl/src/m/enum/DamageIsdamageEnum.m (revision 17763) @@ -0,0 +1,11 @@ +function macro=DamageIsdamageEnum() +%DAMAGEISDAMAGEENUM - Enum of DamageIsdamage +% +% 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=DamageIsdamageEnum() + +macro=StringToEnum('DamageIsdamage'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17762) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17763) @@ -164,4 +164,5 @@ def MaterialsRheologyLawEnum(): return StringToEnum("MaterialsRheologyLaw")[0] def MaterialsRheologyNEnum(): return StringToEnum("MaterialsRheologyN")[0] +def DamageIsdamageEnum(): return StringToEnum("DamageIsdamage")[0] def DamageDEnum(): return StringToEnum("DamageD")[0] def DamageFEnum(): return StringToEnum("DamageF")[0]