Changeset 16198 for issm/trunk-jpl/src/m/classes/damage.py

Timestamp:

09/20/13 07:39:38 (12 years ago)

Author:

cborstad

Message:

CHG: merged matlab changes to python for damage model

File:

• issm/trunk-jpl/src/m/classes/damage.py (modified) (6 diffs)

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

@@ -16 +16 @@
                         
                 #damage: 
-                isdamage = float('NaN')
-                self.D   = float('NaN')
-                self.law   = ''
+                self.D                                          = float('NaN')
+                self.law                                                = ''
+                self.spcdamage                          = float('NaN')
+                self.max_damage                 = float('NaN')
                 
-                #parameters for law 'initial': 
-                self.stress_threshold    = float('NaN')
-                self.c1                  = float('NaN')
-                self.c2                  = float('NaN')
-                self.c3                  = float('NaN')
-                self.c4                  = float('NaN')
+                #numerical
+                stabilization                           = float('NaN')
+                penalty_threshold                       = float('NaN')
+                maxiter                                         = float('NaN')
+                penalty_lock                            = float('NaN')
+                penalty_factor                          = float('NaN')
+
+                #general parameters for evolution law: 
+                self.stress_threshold   = float('NaN')
+                self.c1                 = float('NaN')
+                self.c2                 = float('NaN')
+                self.c3                 = float('NaN')
+                self.c4                 = float('NaN')
+                self.healing                            = float('NaN')
 
                 if not len(args):
@@ -38 +47 @@
                 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_lock","threshold to declare convergence of damage evolution solution (default is 0)")
+                s+="%s\n" % fielddisplay(self,"penalty_lock","scaling exponent (default is 3)")
 
                 if (self.law=='pralong'):
@@ -51 +68 @@
 
                 #damage parameters: 
-                self.isdamage=0
                 self.D=0
                 self.law='undamaged'
+
+                self.max_damage=1-1e-5 #if damage reaches 1, solve becomes singular, as viscosity becomes nil
+                
+                #Type of stabilization used
+                self.stabilization=2
+                        
+                #Maximum number of iterations
+                self.maxiter=100
+
+                #factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor
+                self.penalty_factor=3
+                        
+                #stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)
+                self.penalty_lock=0
+                        
+                #threshold to declare convergence of thermal solution (default is 0)
+                self.penalty_threshold=0
+                
+                #damage evolution parameters 
                 self.stress_threshold=0
                 self.c1=0
@@ -59 +94 @@
                 self.c3=0
                 self.c4=0
+                self.healing=0
 
         # }}}
         def checkconsistency(self,md,solution,analyses):    # {{{
 
-                md = checkfield(md,'damage.D','>=',0,'size',[md.mesh.numberofvertices])
-                if self.isdamage:
-                        md = checkfield(md,'damage.law','values',['undamaged','pralong'])
+                md = checkfield(md,'damage.D','>=',0,'<=',self.max_damage,'size',[md.mesh.numberofvertices])
+                md = checkfield(md,'damage.max_damage','<',1,'>=',0)
+                md = checkfield(md,'damage.law','values',['undamaged','pralong'])
+                md = checkfield(md,'damage.spcdamage','forcing',1)
+                        
+                md = checkfield(md,'damage.stabilization','numel',[1],'values',[0,1,2]);
+                md = checkfield(md,'damage.maxiter','>=0',0);
+                md = checkfield(md,'damage.penalty_factor','>=0',0);
+                md = checkfield(md,'damage.penalty_lock','>=0',0);
+                md = checkfield(md,'damage.penalty_threshold','>=0',0);
+
                 if self.law == 'pralong':
+                        md = checkfield(md,'damage.healing','>=',0);
                         md = checkfield(md,'damage.c1','>=',0)
                         md = checkfield(md,'damage.c2','>=',0)
@@ -72 +117 @@
                         md = checkfield(md,'damage.c4','>=',0)
                         md = checkfield(md,'damage.stress_threshold','>=',0)
+                elif strcmpi(self.law,'undamaged'):
+                        if (solution==DamageEvolutionSolutionEnum):
+                                raise RuntimeError('Invalid evolution law (md.damage.law) for a damage solution');
 
                 return md
@@ -77 +125 @@
         def marshall(self,md,fid):    # {{{
 
-                WriteData(fid,'object',self,'fieldname','isdamage','format','DoubleMat','mattype',1)
                 WriteData(fid,'object',self,'fieldname','D','format','DoubleMat','mattype',1)
-                if self.isdamage:
-                        WriteData(fid,'object',self,'fieldname','law','format','String')
+                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');
+
                 if self.law=='pralong':
                         WriteData(fid,'object',self,'fieldname','c1','format','Double')