source: issm/trunk-jpl/src/py3/mech/damagefrominversion.py@ 23670

import numpy as  np

def damagefrominversion(md):
        '''
        compute ice shelf damage from inversion results

        This routine computes damage based on the analytical formalism of Borstad et
        al. (2013, The Cryosphere).  The model must contain inversion results for
        ice rigidity.  Ice rigidity B is assumed to be parameterized by the ice
        temperature in md.materials.rheology_B. 
        
        Usage:
                damage=damagefrominversion(md)
        
        Example:
                damage=damagefrominversion(md)
        '''

        # check inputs
        if not hasattr(md.results,'strainrate'):
                raise Exception('md.results.strainrate is not present.  Calculate using md=mechanicalproperties(md,vx,vy)')
        if not '2d' in md.mesh.__doc__:
                raise Exception('only 2d (planview) model supported currently')
        if any(md.flowequation.element_equation!=2):
                raise Exception('Warning: the model has some non-SSA elements.  These will be treated like SSA elements')
        if np.ndim(md.results.StressbalanceSolution.MaterialsRheologyBbar)==2:
                Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar.reshape(-1,)
        else:
                Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar
        if np.ndim(md.materials.rheology_B)==2:
                BT=md.materials.rheology_B.reshape(-1,)
        else:
                BT=md.materials.rheology_B

        damage=np.zeros_like(Bi)

        # Damage where Bi softer than B(T)
        pos=np.nonzero(Bi<BT)[0]
        damage[pos]=1.-Bi[pos]/BT[pos]
        
        pos=np.nonzero(damage<0)
        damage[pos]=0

        return damage