Changeset 17763

Timestamp:

04/18/14 07:25:06 (11 years ago)

Author:

cborstad

Message:

CHG: python fixes for NR now that it is no longer necessary to specify damage

Location:

issm/trunk-jpl/src

Files:

• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (5 diffs)
• m/boundaryconditions/SetMarineIceSheetBC.py (modified) (4 diffs)
• m/classes/damage.py (modified) (7 diffs)
• m/enum/DamageIsdamageEnum.m (added)
• m/enum/EnumDefinitions.py (modified) (1 diff)

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -164 +164 @@
         MaterialsRheologyLawEnum,
         MaterialsRheologyNEnum,
+        DamageIsdamageEnum,
         DamageDEnum,
         DamageFEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -172 +172 @@
                 case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw";
                 case MaterialsRheologyNEnum : return "MaterialsRheologyN";
+                case DamageIsdamageEnum : return "DamageIsdamage";
                 case DamageDEnum : return "DamageD";
                 case DamageFEnum : return "DamageF";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -175 +175 @@
               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 +260 @@
               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 +383 @@
               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 +506 @@
               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 +629 @@
               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;

issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py

@@ -32 +32 @@
         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 +43 @@
                 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 +69 @@
         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 +91 @@
         #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
-

issm/trunk-jpl/src/m/classes/damage.py

@@ -17 +17 @@
                         
                 #damage: 
+                self.isdamage           = float('NaN')
                 self.D                                          = float('NaN')
                 self.law                                                = ''
@@ -24 +25 @@
                 #numerical
                 stabilization                           = float('NaN')
+                maxiter                                         = float('NaN')
+                elementinterp           = ''
                 penalty_threshold                       = float('NaN')
-                maxiter                                         = float('NaN')
                 penalty_lock                            = float('NaN')
                 penalty_factor                          = float('NaN')
@@ -47 +49 @@
         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 +79 @@
 
                 #damage parameters: 
+                self.isdamage=0
                 self.D=0
                 self.law='undamaged'
@@ -84 +90 @@
                 #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 +128 @@
         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 +160 @@
         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':

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -164 +164 @@
 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]