Changeset 25499


Ignore:
Timestamp:
08/31/20 14:54:17 (5 years ago)
Author:
jdquinn
Message:

CHG: Committing changes in support of test2004.py (no longer failing, but still need to track down sources of large errors)

Location:
issm/trunk-jpl
Files:
3 added
27 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/classes/friction.m

    r25403 r25499  
    7979                end % }}}
    8080                function marshall(self,prefix,md,fid) % {{{
    81                         yts=md.constants.yts;
    82 
    8381                        WriteData(fid,prefix,'name','md.friction.law','data',1,'format','Integer');
    8482                        if(size(self.coefficient,1)==md.mesh.numberofvertices | size(self.coefficient,1)==md.mesh.numberofvertices+1),
    85                                 mattype=1; tsl = md.mesh.numberofvertices;
     83                                mattype=1;
     84                                tsl = md.mesh.numberofvertices;
    8685                        else
    87                                 mattype=2; tsl = md.mesh.numberofelements;
     86                                mattype=2;
     87                                tsl = md.mesh.numberofelements;
    8888                        end
    8989                        WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',mattype,'timeserieslength',tsl+1,'yts',md.constants.yts);

  • issm/trunk-jpl/src/m/classes/friction.py

    r25161 r25499  
     1import numpy as np
     2
    13from checkfield import checkfield
    24from fielddisplay import fielddisplay
    … …  
    68
    79class friction(object):
    8     '''
    9     FRICTION class definition
     10    """FRICTION class definition
    1011
    11         Usage:
    12             friction = friction()
    13     '''
     12    Usage:
     13        friction = friction()
     14    """
    1415
    1516    def __init__(self):  # {{{
    16         self.coefficient = float('NaN')
    17         self.p = float('NaN')
    18         self.q = float('NaN')
     17        self.coefficient = np.nan
     18        self.p = np.nan
     19        self.q = np.nan
    1920        self.coupling = 0
    20         self.effective_pressure = float('NaN')
     21        self.effective_pressure = np.nan
    2122        self.effective_pressure_limit = 0
     23
    2224        #set defaults
    2325        self.setdefaultparameters()
    24         #self.requested_outputs = []
    2526    #}}}
    2627
    2728    def __repr__(self):  # {{{
    28         string = "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b, \n(effective stress Neff = rho_ice * g * thickness + rho_water * g * base, r = q / p and s = 1 / p)"
     29        s = "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b, \n(effective stress Neff = rho_ice * g * thickness + rho_water * g * base, r = q / p and s = 1 / p)\n"
     30        s = "{}\n".format(fielddisplay(self, "coefficient", "friction coefficient [SI]"))
     31        s = "{}\n".format(fielddisplay(self, "p", "p exponent"))
     32        s = "{}\n".format(fielddisplay(self, "q", "q exponent"))
     33        s = "{}\n".format(fielddisplay(self, 'coupling', 'Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)'))
     34        s = "{}\n".format(fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
     35        s = "{}\n".format(fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
     36       
     37        return s
     38    #}}}
    2939
    30         string = "%s\n%s" % (string, fielddisplay(self, "coefficient", "friction coefficient [SI]"))
    31         string = "%s\n%s" % (string, fielddisplay(self, "p", "p exponent"))
    32         string = "%s\n%s" % (string, fielddisplay(self, "q", "q exponent"))
    33         string = "%s\n%s" % (string, fielddisplay(self, 'coupling', 'Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)'))
    34         string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
    35         string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
    36         #string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    37         return string
     40    def setdefaultparameters(self):  # {{{
     41        self.effective_pressure_limit = 0
     42
     43        return self
    3844    #}}}
    3945
    … …  
    5056    #}}}
    5157
    52     def setdefaultparameters(self):  # {{{
    53         #self.requested_outputs = ['default']
    54         self.effective_pressure_limit = 0
    55         return self
    56     #}}}
    57 
    5858    def defaultoutputs(self, md):  # {{{
    5959        list = []
    … …  
    6464        #Early return
    6565        if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
     66            return md
     67        if solution == 'TransientSoluition' and not md.transient.isstressbalance and not md.transient.isthermal:
    6668            return md
    6769
    … …  
    7476            md = checkfield(md, 'fieldname', 'friction.effective_pressure', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    7577        elif self.coupling > 4:
    76             raise ValueError('md.friction.coupling larger than 4, not supported yet')
    77         #md = checkfield(md, 'fieldname', 'friction.requested_outputs', 'stringrow', 1)
     78            raise ValueError('not supported yet')
     79
    7880        return md
    7981    # }}}
    … …  
    8183    def marshall(self, prefix, md, fid):  # {{{
    8284        WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 1, 'format', 'Integer')
    83         WriteData(fid, prefix, 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', 1)
     85
     86        if type(self.coefficient) in [np.ndarray] and (self.coefficient.shape[0] == md.mesh.numberofvertices or self.coefficient.shape[0] == (md.mesh.numberofvertices + 1)):
     87            mattype = 1
     88            tsl = md.mesh.numberofvertices
     89        else:
     90            mattype = 2
     91            tsl = md.mesh.numberofelements
     92
     93        WriteData(fid, prefix, 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', mattype, 'timeserieslength', tsl + 1, 'yts', md.constants.yts)
    8494        WriteData(fid, prefix, 'object', self, 'fieldname', 'p', 'format', 'DoubleMat', 'mattype', 2)
    8595        WriteData(fid, prefix, 'object', self, 'fieldname', 'q', 'format', 'DoubleMat', 'mattype', 2)
    … …  
    8898        if self.coupling == 3:
    8999            WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    90         if self.coupling == 4:
    91             WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1)
     100        elif self.coupling == 4:
     101            WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    92102        elif self.coupling > 4:
    93             raise ValueError('md.friction.coupling larger than 4, not supported yet')
    94 
    95     # #process requested outputs
    96     #     outputs = self.requested_outputs
    97     #     indices = [i for i, x in enumerate(outputs) if x == 'default']
    98     #     if len(indices) > 0:
    99     #         outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
    100     #         outputs = outputscopy
    101     #     WriteData(fid, prefix, 'data', outputs, 'name', 'md.friction.requested_outputs', 'format', 'StringArray')
     103            raise ValueError('not supported yet')
    102104    # }}}

  • issm/trunk-jpl/src/m/classes/frictioncoulomb.py

    r24670 r25499  
     1import numpy as np
     2
     3from checkfield import checkfield
    14from fielddisplay import fielddisplay
    25from project3d import project3d
    3 from checkfield import checkfield
    46from WriteData import WriteData
    57
    68
    79class frictioncoulomb(object):
    8     """
    9     FRICTIONCOULOMB class definition
     10    """FRICTIONCOULOMB class definition
    1011
    1112    Usage:
    12     frictioncoulomb = frictioncoulomb()
     13        frictioncoulomb = frictioncoulomb()
    1314    """
    1415
    1516    def __init__(self):  # {{{
    16         self.coefficient = float('NaN')
    17         self.coefficientcoulomb = float('NaN')
    18         self.p = float('NaN')
    19         self.q = float('NaN')
     17        self.coefficient = np.nan
     18        self.coefficientcoulomb = np.nan
     19        self.p = np.nan
     20        self.q = np.nan
    2021        self.coupling = 0
    21         self.effective_pressure = float('NaN')
     22        self.effective_pressure = np.nan
    2223        self.effective_pressure_limit = 0
    23     #set defaults
     24
     25        #set defaults
    2426        self.setdefaultparameters()
    2527    #}}}
    2628
    2729    def __repr__(self):  # {{{
    28         string = "Basal shear stress parameters: Sigma_b = min(coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b, \n coefficientcoulomb^2 * rho_i * g * (h - h_f)), (effective stress Neff = rho_ice * g * thickness + rho_water * g * bed, r = q / p and s = 1 / p)."
    29         string = "%s\n%s" % (string, fielddisplay(self, "coefficient", "power law (Weertman) friction coefficient [SI]"))
    30         string = "%s\n%s" % (string, fielddisplay(self, "coefficientcoulomb", "Coulomb friction coefficient [SI]"))
    31         string = "%s\n%s" % (string, fielddisplay(self, "p", "p exponent"))
    32         string = "%s\n%s" % (string, fielddisplay(self, "q", "q exponent"))
    33         string = "%s\n%s" % (string, fielddisplay(self, 'coupling', 'Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure)  and 2 for coupled(not implemented yet)'))
    34         string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
    35         string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
    36         return string
     30        s = "Basal shear stress parameters: Sigma_b = min(coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b, \n coefficientcoulomb^2 * rho_i * g * (h - h_f)), (effective stress Neff = rho_ice * g * thickness + rho_water * g * bed, r = q / p and s = 1 / p).\n"
     31        s = "{}\n".format(fielddisplay(self, "coefficient", "power law (Weertman) friction coefficient [SI]"))
     32        s = "{}\n".format(fielddisplay(self, "coefficientcoulomb", "Coulomb friction coefficient [SI]"))
     33        s = "{}\n".format(fielddisplay(self, "p", "p exponent"))
     34        s = "{}\n".format(fielddisplay(self, "q", "q exponent"))
     35        s = "{}\n".format(fielddisplay(self, 'coupling', 'Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure)  and 2 for coupled(not implemented yet)'))
     36        s = "{}\n".format(fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
     37        s = "{}\n".format(fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
     38       
     39        return s
     40    #}}}
     41
     42    def setdefaultparameters(self):  # {{{
     43        self.effective_pressure_limit = 0
     44
     45        return self
    3746    #}}}
    3847
    … …  
    4554            self.effective_pressure = project3d(md, 'vector', self.effective_pressure, 'type', 'node', 'layer', 1)
    4655        elif self.coupling == 2:
    47             raise ValueError('coupling not supported yet')
     56            raise ValueError('not implemented yet')
    4857        elif self.coupling > 2:
    49             raise ValueError('md.friction.coupling larger than 2, not supported yet')
     58            raise ValueError('not supported yet')
     59
    5060        return self
    5161    #}}}
    52 
    53     def setdefaultparameters(self):  # {{{
    54         self.effective_pressure_limit = 0
    55         return self
    56     #}}}
    57 
    5862    def checkconsistency(self, md, solution, analyses):  # {{{
    5963        #Early return
    … …  
    6569        md = checkfield(md, 'fieldname', 'friction.q', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofelements])
    6670        md = checkfield(md, 'fieldname', 'friction.p', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofelements])
     71        md = checkfield(md, 'fieldname', 'friction.coupling', 'numel', [1], 'values', [0, 1, 2])
    6772        md = checkfield(md, 'fieldname', 'friction.effective_pressure_limit', 'numel', [1], '>=', 0)
    6873        if self.coupling == 1:
    6974            md = checkfield(md, 'fieldname', 'friction.effective_pressure', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    7075        elif self.coupling == 2:
    71             raise ValueError('coupling not supported yet')
     76            raise ValueError('not implemented yet')
    7277        elif self.coupling > 2:
    73             raise ValueError('md.friction.coupling larger than 2, not supported yet')
     78            raise ValueError('not supported yet')
     79
    7480        return md
    7581    # }}}
    … …  
    8692            WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    8793        elif self.coupling == 2:
    88             raise ValueError('coupling not supported yet')
     94            raise ValueError('not implemented yet')
    8995        elif self.coupling > 2:
    90             raise ValueError('md.friction.coupling larger than 2, not supported yet')
     96            raise ValueError('not supported yet')
    9197    # }}}

  • issm/trunk-jpl/src/m/classes/frictionshakti.m

    r23020 r25499  
    3636                end % }}}
    3737                function marshall(self,prefix,md,fid) % {{{
    38                         yts=md.constants.yts;
    39 
    4038                        WriteData(fid,prefix,'name','md.friction.law','data',8,'format','Integer');
    4139                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);

  • issm/trunk-jpl/src/m/classes/frictionshakti.py

    r24213 r25499  
     1import numpy as np
     2
     3from checkfield import checkfield
    14from fielddisplay import fielddisplay
    25from project3d import project3d
    3 from checkfield import checkfield
    46from WriteData import WriteData
    57
    68
    79class frictionshakti(object):
    8     """
    9     FRICTIONSHAKTI class definition
     10    """FRICTIONSHAKTI class definition
    1011
    1112    Usage:
    12     friction = frictionshakti()
     13        friction = frictionshakti()
    1314    """
    1415
    1516    def __init__(self, md):  # {{{
    16         self.coefficient = md.friction.coefficient
    17     #set defaults
     17        self.coefficient = np.nan
     18        #set defaults
    1819        self.setdefaultparameters()
    1920    #}}}
    2021
    2122    def __repr__(self):  # {{{
    22         string = "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff = rho_ice * g * thickness + rho_water * g * (head - b))"
    23         string = "%s\n%s" % (string, fielddisplay(self, "coefficient", "friction coefficient [SI]"))
    24         return string
     23        s = "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff = rho_ice * g * thickness + rho_water * g * (head - b))\n"
     24        s += "{}\n".format(fielddisplay(self, "coefficient", "friction coefficient [SI]"))
     25
     26        return s
     27    #}}}
     28
     29    def setdefaultparameters(self):  # {{{
     30        return self
    2531    #}}}
    2632
    … …  
    3036    #}}}
    3137
    32     def setdefaultparameters(self):  # {{{
    33         return self
    34     #}}}
    35 
    3638    def checkconsistency(self, md, solution, analyses):  # {{{
    37         #Early return
     39        # Early return
    3840        if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
    3941            return md

  • issm/trunk-jpl/src/m/classes/frictionwaterlayer.py

    r24213 r25499  
     1import numpy as np
     2
     3from checkfield import checkfield
     4from fielddisplay import fielddisplay
    15from project3d import project3d
    2 from fielddisplay import fielddisplay
    3 from checkfield import checkfield
    46from WriteData import WriteData
    57
    68
    79class frictionwaterlayer(object):
    8     """
    9     frictionwaterlayer class definition
     10    """FRICTIONWATERLAYER class definition
    1011
    1112    Usage:
    1213        friction = frictionwaterlayer(md)
    1314    """
     15
    1416    def __init__(self, md):  # {{{
    15         self.coefficient = md.friction.coefficient
    16         self.f = float('NaN')
    17         self.p = md.friction.p
    18         self.q = md.friction.q
    19         self.water_layer = float('NaN')
     17        self.coefficient = np.nan
     18        self.f = np.nan
     19        self.p = np.nan
     20        self.q = np.nan
     21        self.water_layer = np.nan
     22    #}}}
     23
     24    def __repr__(self):  # {{{
     25        s = 'Basal shear stress parameters: tau_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b * 1 / f(T)\n(effective stress Neff = rho_ice * g * thickness + rho_water * g * (bed + water_layer), r = q / p and s = 1 / p)\n'
     26        s = "{}\n".format(fielddisplay(self, 'coefficient', 'frictiontemp coefficient [SI]'))
     27        s = "{}\n".format(fielddisplay(self, 'f', 'f variable for effective pressure'))
     28        s = "{}\n".format(fielddisplay(self, 'p', 'p exponent'))
     29        s = "{}\n".format(fielddisplay(self, 'q', 'q exponent'))
     30        s = "{}\n".format(fielddisplay(self, 'water_layer', 'water thickness at the base of the ice (m)'))
     31
     32        return s
     33    #}}}
     34
     35    def setdefaultparameters(self):  # {{{
     36        return self
    2037    #}}}
    2138
    … …  
    3754        self.q = project3d(md, 'vector', self.q, 'type', 'element')
    3855        self.water_layer = project3d(md, 'vector', self.water_layer, 'type', 'node', 'layer', 1)
     56
    3957        return self
    4058    # }}}
    4159
    42     def __repr__(self):  # {{{
    43         string = 'Basal shear stress parameters: tau_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b * 1 / f(T)\n(effective stress Neff = rho_ice * g * thickness + rho_water * g * (bed + water_layer), r = q / p and s = 1 / p)'
    44         string = "%s\n%s" % (string, fielddisplay(self, 'coefficient', 'frictiontemp coefficient [SI]'))
    45         string = "%s\n%s" % (string, fielddisplay(self, 'f', 'f variable for effective pressure'))
    46         string = "%s\n%s" % (string, fielddisplay(self, 'p', 'p exponent'))
    47         string = "%s\n%s" % (string, fielddisplay(self, 'q', 'q exponent'))
    48         string = "%s\n%s" % (string, fielddisplay(self, 'water_layer', 'water thickness at the base of the ice (m)'))
    49 
    50         return string
    51     #}}}
    52 
    5360    def marshall(self, prefix, md, fid):  #{{{
    54         yts = md.constants.yts
    5561        WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 5, 'format', 'Integer')
    56         WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
     62        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    5763        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'f', 'format', 'Double')
    5864        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'p', 'format', 'DoubleMat', 'mattype', 2)
    5965        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'q', 'format', 'DoubleMat', 'mattype', 2)
    60         WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'water_layer', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
     66        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'water_layer', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    6167    #}}}

  • issm/trunk-jpl/src/m/classes/mesh3dsurface.py

    r25455 r25499  
    207207
    208208                #first line:
    209                 contour                 = OrderedStruct()
    210                 contour.x               = [self.long[el[0]], self.long[el[1]]]
    211                 contour.y               = [self.lat[el[0]], self.lat[el[1]]]
    212                 contour.Geometry        = 'Line'
    213                 contours[count]       = contour
     209                contour             = OrderedStruct()
     210                contour.x           = [self.long[el[0]], self.long[el[1]]]
     211                contour.y           = [self.lat[el[0]], self.lat[el[1]]]
     212                contour.Geometry    = 'Line'
     213                contours[count]     = contour
    214214
    215215                #second line:
    216                 contour                 = OrderedStruct()
    217                 contour.x               = [self.long[el[1]], self.long[el[2]]]
    218                 contour.y               = [self.lat[el[1]], self.lat[el[2]]]
    219                 contour.Geometry        = 'Line'
    220                 contours[count + 1]   = contour
     216                contour             = OrderedStruct()
     217                contour.x           = [self.long[el[1]], self.long[el[2]]]
     218                contour.y           = [self.lat[el[1]], self.lat[el[2]]]
     219                contour.Geometry    = 'Line'
     220                contours[count + 1] = contour
    221221
    222222                #second line:
    223                 contour                 = OrderedStruct()
    224                 contour.x               = [self.long[el[2]], self.long[el[0]]]
    225                 contour.y               = [self.lat[el[2]], self.lat[el[0]]]
    226                 contour.Geometry        = 'Line'
    227                 contours[count + 2]   = contour
     223                contour             = OrderedStruct()
     224                contour.x           = [self.long[el[2]], self.long[el[0]]]
     225                contour.y           = [self.lat[el[2]], self.lat[el[0]]]
     226                contour.Geometry    = 'Line'
     227                contours[count + 2] = contour
    228228
    229229                #increase count:

  • issm/trunk-jpl/src/m/classes/model.m

    r25172 r25499  
    214214                        %
    215215                        %   This routine collapses a 3d model into a 2d model
    216                         %   and collapses all the fileds of the 3d model by
     216                        %   and collapses all the fields of the 3d model by
    217217                        %   taking their depth-averaged values
    218218                        %
    … …  
    261261
    262262                        %observations
    263                         if ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); end;
    264                         if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end;
    265                         if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end;
    266                         if ~isnan(md.inversion.thickness_obs), md.inversion.thickness_obs=project2d(md,md.inversion.thickness_obs,md.mesh.numberoflayers); end;
    267                         if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end;
    268                         if numel(md.inversion.min_parameters)>1, md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end;
    269                         if numel(md.inversion.max_parameters)>1, md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end;
     263                        if ~isnan(md.inversion.vx_obs),
     264                                md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers);
     265                        end
     266                        if ~isnan(md.inversion.vy_obs),
     267                                md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers);
     268                        end
     269                        if ~isnan(md.inversion.vel_obs),
     270                                md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers);
     271                        end
     272                        if ~isnan(md.inversion.thickness_obs),
     273                                md.inversion.thickness_obs=project2d(md,md.inversion.thickness_obs,md.mesh.numberoflayers);
     274                        end
     275                        if ~isnan(md.inversion.cost_functions_coefficients),
     276                                md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers);
     277                        end
     278                        if numel(md.inversion.min_parameters)>1,
     279                                md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers);
     280                        end
     281                        if numel(md.inversion.max_parameters)>1,
     282                                md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers);
     283                        end
    270284                        if isa(md.smb,'SMBforcing') & ~isnan(md.smb.mass_balance),
    271285                                md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers);
    272286                        elseif isa(md.smb,'SMBhenning') & ~isnan(md.smb.smbref),
    273287                                md.smb.smbref=project2d(md,md.smb.smbref,md.mesh.numberoflayers);
    274                         end;
     288                        end
    275289
    276290                        %results
    277                         if ~isnan(md.initialization.vx),md.initialization.vx=DepthAverage(md,md.initialization.vx);end;
    278                         if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end;
    279                         if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end;
    280                         if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end;
    281                         if ~isnan(md.initialization.temperature),md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end;
    282                         if ~isnan(md.initialization.pressure),md.initialization.pressure=project2d(md,md.initialization.pressure,1);end;
    283                         if ~isnan(md.initialization.sediment_head),md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1);end;
    284                         if ~isnan(md.initialization.epl_head),md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);end;
    285                         if ~isnan(md.initialization.epl_thickness),md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);end;
    286                         if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project2d(md,md.initialization.waterfraction,1);end;
    287                         if ~isnan(md.initialization.watercolumn),md.initialization.watercolumn=project2d(md,md.initialization.watercolumn,1);end;
     291                        if ~isnan(md.initialization.vx),
     292                                md.initialization.vx=DepthAverage(md,md.initialization.vx);
     293                        end
     294                        if ~isnan(md.initialization.vy),
     295                                md.initialization.vy=DepthAverage(md,md.initialization.vy);
     296                        end
     297                        if ~isnan(md.initialization.vz),
     298                                md.initialization.vz=DepthAverage(md,md.initialization.vz);
     299                        end
     300                        if ~isnan(md.initialization.vel),
     301                                md.initialization.vel=DepthAverage(md,md.initialization.vel);
     302                        end
     303                        if ~isnan(md.initialization.temperature),
     304                                md.initialization.temperature=DepthAverage(md,md.initialization.temperature);
     305                        end
     306                        if ~isnan(md.initialization.pressure),
     307                                md.initialization.pressure=project2d(md,md.initialization.pressure,1);
     308                        end
     309                        if ~isnan(md.initialization.sediment_head),
     310                                md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1);
     311                        end
     312                        if ~isnan(md.initialization.epl_head),
     313                                md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);
     314                        end
     315                        if ~isnan(md.initialization.epl_thickness),
     316                                md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);
     317                        end
     318                        if ~isnan(md.initialization.waterfraction),
     319                                md.initialization.waterfraction=project2d(md,md.initialization.waterfraction,1);
     320                        end
     321                        if ~isnan(md.initialization.watercolumn),
     322                                md.initialization.watercolumn=project2d(md,md.initialization.watercolumn,1);
     323                        end
     324
    288325                        %giaivins
    289326                        if isa(md.gia,'giaivins'),
    … …  
    299336                                md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1);
    300337                                md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1);
    301                         end     
     338                        end
    302339
    303340                        %boundary conditions
    … …  
    307344                        md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers);
    308345                        md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers);
    309                         if numel(md.masstransport.spcthickness)>1 md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers); end
    310                         if numel(md.damage.spcdamage)>1, md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers); end
    311                         if numel(md.levelset.spclevelset)>1, md.levelset.spclevelset=project2d(md,md.levelset.spclevelset,md.mesh.numberoflayers); end
     346                        if numel(md.masstransport.spcthickness)>1,
     347                                md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers);
     348                        end
     349                        if numel(md.damage.spcdamage)>1,
     350                                md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers);
     351                        end
     352                        if numel(md.levelset.spclevelset)>1,
     353                                md.levelset.spclevelset=project2d(md,md.levelset.spclevelset,md.mesh.numberoflayers);
     354                        end
    312355                        md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
    313356
    … …  
    362405                                md.geometry.bed=project2d(md,md.geometry.bed,1);
    363406                        end
    364 
    365407                        if ~isnan(md.mask.ocean_levelset),
    366408                                md.mask.ocean_levelset=project2d(md,md.mask.ocean_levelset,1);
    … …  
    371413
    372414                        %lat long
    373                         if numel(md.mesh.lat) ==md.mesh.numberofvertices,  md.mesh.lat=project2d(md,md.mesh.lat,1); end
    374                         if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end
     415                        if numel(md.mesh.lat)==md.mesh.numberofvertices,
     416                                md.mesh.lat=project2d(md,md.mesh.lat,1);
     417                        end
     418                        if numel(md.mesh.long)==md.mesh.numberofvertices,
     419                                md.mesh.long=project2d(md,md.mesh.long,1);
     420                        end
    375421
    376422                        %outputdefinitions
    … …  
    388434                        end
    389435
    390                         %Initialize with the 2d mesh
     436                        %Initialize the 2d mesh
    391437                        mesh=mesh2d();
    392438                        mesh.x=md.mesh.x2d;
    … …  
    395441                        mesh.numberofelements=md.mesh.numberofelements2d;
    396442                        mesh.elements=md.mesh.elements2d;
    397                         if numel(md.mesh.lat) ==md.mesh.numberofvertices,  mesh.lat=project2d(md,md.mesh.lat,1); end
    398                         if numel(md.mesh.long)==md.mesh.numberofvertices,       mesh.long=project2d(md,md.mesh.long,1); end
     443                        if numel(md.mesh.lat)==md.mesh.numberofvertices,
     444                                mesh.lat=project2d(md,md.mesh.lat,1);
     445                        end
     446                        if numel(md.mesh.long)==md.mesh.numberofvertices,
     447                                mesh.long=project2d(md,md.mesh.long,1);
     448                        end
    399449                        mesh.epsg=md.mesh.epsg;
    400                         if numel(md.mesh.scale_factor)==md.mesh.numberofvertices,       mesh.scale_factor=project2d(md,md.mesh.scale_factor,1); end
    401                         if ~isnan(md.mesh.vertexonboundary), mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); end
    402                         if ~isnan(md.mesh.elementconnectivity), mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1); end
     450                        if numel(md.mesh.scale_factor)==md.mesh.numberofvertices,
     451                                mesh.scale_factor=project2d(md,md.mesh.scale_factor,1);
     452                        end
     453                        if ~isnan(md.mesh.vertexonboundary),
     454                                mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1);
     455                        end
     456                        if ~isnan(md.mesh.elementconnectivity),
     457                                mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1);
     458                        end
    403459                        md.mesh=mesh;
    404460                        md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
    … …  
    622678                                md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
    623679                                md2.mesh.segments=contourenvelope(md2.mesh);
    624                                 md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
     680                                md2.mesh.vertexonboundary=zeros(numberofvertices2,1);
     681                                md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
    625682                        else
    626683                                %First do the connectivity for the contourenvelope in 2d
    … …  
    628685                                md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity);
    629686                                segments=contourenvelope(md2.mesh);
    630                                 md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(segments(:,1:2))=1;
     687                                md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1);
     688                                md2.mesh.vertexonboundary(segments(:,1:2))=1;
    631689                                md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1);
    632690                                %Then do it for 3d as usual

  • issm/trunk-jpl/src/m/classes/model.py

    r25455 r25499  
    412412            md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements, md2.mesh.numberofvertices)
    413413            md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements, md2.mesh.vertexconnectivity)
    414             md2.mesh.segments = contourenvelope(md2)
     414            md2.mesh.segments = contourenvelope(md2.mesh)
    415415            md2.mesh.vertexonboundary = np.zeros(numberofvertices2, bool)
    416416            md2.mesh.vertexonboundary[md2.mesh.segments[:, 0:2] - 1] = True
    … …  
    419419            md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements2d, md2.mesh.numberofvertices2d)
    420420            md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements2d, md2.mesh.vertexconnectivity)
    421             segments = contourenvelope(md2)
     421            msegments = contourenvelope(md2.mesh)
    422422            md2.mesh.vertexonboundary = np.zeros(int(numberofvertices2 / md2.mesh.numberoflayers), bool)
    423423            md2.mesh.vertexonboundary[segments[:, 0:2] - 1] = True
    … …  
    696696
    697697        This routine collapses a 3d model into a 2d model and collapses all
    698         the fileds of the 3d model by taking their depth - averaged values
     698        the fields of the 3d model by taking their depth-averaged values
    699699
    700700        Usage:
    701701            md = collapse(md)
     702
     703        See also: EXTRUDE, MODELEXTRACT
    702704        """
    703705
    704         #Check that the model is really a 3d model
    705         if md.mesh.domaintype().lower() != '3d':
    706             raise Exception("only a 3D model can be collapsed")
    707 
    708         #dealing with the friction law
    709         #drag is limited to nodes that are on the bedrock.
    710         if hasattr(md.friction, 'coefficient'):
     706        # Check that the model is really a 3d model
     707        if md.mesh.elementtype() != 'Penta':
     708            raise Exception("collapse error message: only a 3d mesh can be collapsed")
     709
     710        # Start with changing all the fields from the 3d mesh
     711
     712        # Dealing with the friction law
     713        # Drag is limited to nodes that are on the bedrock.
     714        if md.friction.__class__.__name__ == 'friction':
    711715            md.friction.coefficient = project2d(md, md.friction.coefficient, 1)
    712 
    713         #p and q (same deal, except for element that are on the bedrock: )
    714         if hasattr(md.friction, 'p'):
    715716            md.friction.p = project2d(md, md.friction.p, 1)
    716         if hasattr(md.friction, 'q'):
    717717            md.friction.q = project2d(md, md.friction.q, 1)
    718 
    719         if hasattr(md.friction, 'coefficientcoulomb'):
     718        elif md.friction.__class__.__name__ == 'frictioncoulomb':
     719            md.friction.coefficient = project2d(md, md.friction.coefficient, 1)
    720720            md.friction.coefficientcoulomb = project2d(md, md.friction.coefficientcoulomb, 1)
    721         if hasattr(md.friction, 'C'):
     721            md.friction.p = project2d(md, md.friction.p, 1)
     722            md.friction.q = project2d(md, md.friction.q, 1)
     723        elif md.friction.__class__.__name__ == 'frictionhydro':
     724            md.friction.q = project2d(md, md.friction.q, 1)
    722725            md.friction.C = project2d(md, md.friction.C, 1)
    723         if hasattr(md.friction, 'As'):
    724726            md.friction.As = project2d(md, md.friction.As, 1)
    725         if hasattr(md.friction, 'effective_pressure') and not np.isnan(md.friction.effective_pressure).all():
    726727            md.friction.effective_pressure = project2d(md, md.friction.effective_pressure, 1)
    727         if hasattr(md.friction, 'water_layer'):
     728        elif md.friction.__class__.__name__ == 'frictionwaterlayer':
     729            md.friction.coefficient = project2d(md, md.friction.coefficient, 1)
     730            md.friction.p = project2d(md, md.friction.p, 1)
     731            md.friction.q = project2d(md, md.friction.q, 1)
    728732            md.friction.water_layer = project2d(md, md.friction.water_layer, 1)
    729         if hasattr(md.friction, 'm'):
     733        elif md.friction.__class__.__name__ == 'frictionweertman':
     734            md.friction.C = project2d(md, md.friction.C, 1)
    730735            md.friction.m = project2d(md, md.friction.m, 1)
    731 
    732         #observations
     736        elif md.friction.__class__.__name__ == 'frictionweertmantemp':
     737            md.friction.C = project2d(md, md.friction.C, 1)
     738            md.friction.m = project2d(md, md.friction.m, 1)
     739        else:
     740            print('friction type not supported')
     741
     742        # Observations
    733743        if not np.isnan(md.inversion.vx_obs).all():
    734744            md.inversion.vx_obs = project2d(md, md.inversion.vx_obs, md.mesh.numberoflayers)
    … …  
    741751        if not np.isnan(md.inversion.cost_functions_coefficients).all():
    742752            md.inversion.cost_functions_coefficients = project2d(md, md.inversion.cost_functions_coefficients, md.mesh.numberoflayers)
    743         if isinstance(md.inversion.min_parameters, np.ndarray):
    744             if md.inversion.min_parameters.size > 1:
    745                 md.inversion.min_parameters = project2d(md, md.inversion.min_parameters, md.mesh.numberoflayers)
    746         if isinstance(md.inversion.max_parameters, np.ndarray):
    747             if md.inversion.max_parameters.size > 1:
    748                 md.inversion.max_parameters = project2d(md, md.inversion.max_parameters, md.mesh.numberoflayers)
    749         if hasattr(md.smb, 'mass_balance'):
    750             if not np.isnan(md.smb.mass_balance).all():
     753        if isinstance(md.inversion.min_parameters, np.ndarray) and md.inversion.min_parameters.size > 1:
     754            md.inversion.min_parameters = project2d(md, md.inversion.min_parameters, md.mesh.numberoflayers)
     755        if isinstance(md.inversion.max_parameters, np.ndarray) and md.inversion.max_parameters.size > 1:
     756            md.inversion.max_parameters = project2d(md, md.inversion.max_parameters, md.mesh.numberoflayers)
     757        if md.smb.__class__.__name__ == 'SMBforcing' and not np.isnan(md.smb.mass_balance).all():
    751758                md.smb.mass_balance = project2d(md, md.smb.mass_balance, md.mesh.numberoflayers)
    752 
    753         #results
     759        elif md.smb.__class__.__name__ == 'SMBhenning' and not np.isnan(md.smb.smbref).all():
     760                md.smb.smbref = project2d(md, md.smb.smbref, md.mesh.numberoflayers)
     761
     762        # Results
    754763        if not np.isnan(md.initialization.vx).all():
    755764            md.initialization.vx = DepthAverage(md, md.initialization.vx)
    … …  
    770779        if not np.isnan(md.initialization.epl_thickness).all():
    771780            md.initialization.epl_thickness = project2d(md, md.initialization.epl_thickness, 1)
    772 
    773         #giaivins
     781        if not np.isnan(md.initialization.waterfraction).all():
     782            md.initialization.waterfraction = project2d(md, md.initialization.waterfraction, 1)
     783        if not np.isnan(md.initialization.watercolumn).all():
     784            md.initialization.watercolumn = project2d(md, md.initialization.watercolumn, 1)
     785
     786        # giaivins
    774787        if md.gia.__class__.__name__ == 'giaivins':
    775788            if not np.isnan(md.gia.mantle_viscosity).all():
    … …  
    778791                md.gia.lithosphere_thickness = project2d(md, md.gia.lithosphere_thickness, 1)
    779792
    780         #elementstype
     793        # elementstype
    781794        if not np.isnan(md.flowequation.element_equation).all():
    782795            md.flowequation.element_equation = project2d(md, md.flowequation.element_equation, 1)
    … …  
    786799            md.flowequation.borderFS = project2d(md, md.flowequation.borderFS, 1)
    787800
    788         # Hydrologydc variables
    789         hydrofields = md.hydrology.__dict__.keys()
    790         for field in hydrofields:
    791             try:
    792                 isvector = np.logical_or(np.shape(md.hydrology.__dict__[field])[0] == md.mesh.numberofelements,
    793                                          np.shape(md.hydrology.__dict__[field])[0] == md.mesh.numberofvertices)
    794             except IndexError:
    795                 isvector = False
    796             #we colpase only fields that are vertices or element based
    797             if isvector:
    798                 md.hydrology.__dict__[field] = project2d(md, md.hydrology.__dict__[field], 1)
    799 
    800         #boundary conditions
     801        # Boundary conditions
    801802        md.stressbalance.spcvx = project2d(md, md.stressbalance.spcvx, md.mesh.numberoflayers)
    802803        md.stressbalance.spcvy = project2d(md, md.stressbalance.spcvy, md.mesh.numberoflayers)
    … …  
    804805        md.stressbalance.referential = project2d(md, md.stressbalance.referential, md.mesh.numberoflayers)
    805806        md.stressbalance.loadingforce = project2d(md, md.stressbalance.loadingforce, md.mesh.numberoflayers)
     807        # TODO:
     808        # - Check if md.mesh.numberoflayershould really be offset by 1.
     809        # - Find out why md.masstransport.spcthickness is not offset, but the
     810        #   other fields are.
     811        # - If offset is required, figure out if it can be abstarcted away to
     812        #   another part of the API.
    806813        if np.size(md.masstransport.spcthickness) > 1:
    807814            md.masstransport.spcthickness = project2d(md, md.masstransport.spcthickness, md.mesh.numberoflayers)
    … …  
    812819        md.thermal.spctemperature = project2d(md, md.thermal.spctemperature, md.mesh.numberoflayers - 1)
    813820
    814         #materials
     821        # Hydrologydc variables
     822        if md.hydrology.__class__.__name__ == 'hydrologydc':
     823            md.hydrology.spcsediment_head = project2d(md, md.hydrology.spcsediment_head, 1)
     824            md.hydrology.mask_eplactive_node = project2d(md, md.hydrology.mask_eplactive_node, 1)
     825            md.hydrology.sediment_transmitivity = project2d(md, md.hydrology.sediment_transmitivity, 1)
     826            md.hydrology.basal_moulin_input = project2d(md, md.hydrology.basal_moulin_input, 1)
     827            if md.hydrology.isefficientlayer == 1:
     828                md.hydrology.spcepl_head = project2d(md, md.hydrology.spcepl_head, 1)
     829            # hydrofields = md.hydrology.__dict__.keys()
     830            # for field in hydrofields:
     831            #     try:
     832            #         isvector = np.logical_or(np.shape(md.hydrology.__dict__[field])[0] == md.mesh.numberofelements,
     833            #                                  np.shape(md.hydrology.__dict__[field])[0] == md.mesh.numberofvertices)
     834            #     except IndexError:
     835            #         isvector = False
     836            #     #we collapse only fields that are vertices or element based
     837            #     if isvector:
     838            #         md.hydrology.__dict__[field] = project2d(md, md.hydrology.__dict__[field], 1)
     839
     840        # Materials
    815841        md.materials.rheology_B = DepthAverage(md, md.materials.rheology_B)
    816842        md.materials.rheology_n = project2d(md, md.materials.rheology_n, 1)
    817843
    818         #damage:
     844        # Damage
    819845        if md.damage.isdamage:
    820846            md.damage.D = DepthAverage(md, md.damage.D)
    821847
    822         #special for thermal modeling:
    823         md.basalforcings.groundedice_melting_rate = project2d(md, md.basalforcings.groundedice_melting_rate, 1)
    824         md.basalforcings.floatingice_melting_rate = project2d(md, md.basalforcings.floatingice_melting_rate, 1)
     848        # Special for thermal modeling
     849        if not np.isnan(md.basalforcings.groundedice_melting_rate).all():
     850            md.basalforcings.groundedice_melting_rate = project2d(md, md.basalforcings.groundedice_melting_rate, 1)
     851        if hasattr(md.basalforcings, 'floatingice_melting_rate') and not np.isnan(md.basalforcings.floatingice_melting_rate).all():
     852            md.basalforcings.floatingice_melting_rate = project2d(md, md.basalforcings.floatingice_melting_rate, 1)
    825853        md.basalforcings.geothermalflux = project2d(md, md.basalforcings.geothermalflux, 1)  #bedrock only gets geothermal flux
    826854
    827         #update of connectivity matrix
     855        if hasattr(md.calving, 'coeff') and not np.isnan(md.calving.coeff).all():
     856            md.calving.coeff = project2d(md, md.calving.coeff, 1)
     857        if hasattr(md.frontalforcings, 'meltingrate') and not np.isnan(md.frontalforcings.meltingrate).all():
     858            md.frontalforcings.meltingrate = project2d(md, md.frontalforcings.meltingrate, 1)
     859
     860        # Update of connectivity matrix
    828861        md.mesh.average_vertex_connectivity = 25
    829862
    830         #parameters
     863        # Collapse the mesh
     864        nodes2d = md.mesh.numberofvertices2d
     865        elements2d = md.mesh.numberofelements2d
     866
     867        # Parameters
    831868        md.geometry.surface = project2d(md, md.geometry.surface, 1)
    832869        md.geometry.thickness = project2d(md, md.geometry.thickness, 1)
    833870        md.geometry.base = project2d(md, md.geometry.base, 1)
    834         if isinstance(md.geometry.bed, np.ndarray):
     871        if not np.isnan(md.geometry.bed).all():
    835872            md.geometry.bed = project2d(md, md.geometry.bed, 1)
    836         md.mask.ocean_levelset = project2d(md, md.mask.ocean_levelset, 1)
    837         md.mask.ice_levelset = project2d(md, md.mask.ice_levelset, 1)
    838 
    839         #OutputDefinitions
     873        if not np.isnan(md.mask.ocean_levelset).all():
     874            md.mask.ocean_levelset = project2d(md, md.mask.ocean_levelset, 1)
     875        if not np.isnan(md.mask.ice_levelset).all():
     876            md.mask.ice_levelset = project2d(md, md.mask.ice_levelset, 1)
     877
     878        # Lat/long
     879        if np.size(md.mesh.lat) == md.mesh.numberofvertices:
     880            md.mesh.lat = project2d(md, md.mesh.lat, 1)
     881        if np.size(md.mesh.long) == md.mesh.numberofvertices:
     882            md.mesh.long = project2d(md, md.mesh.long, 1)
     883
     884        # OutputDefinitions
    840885        if md.outputdefinition.definitions:
    841886            for solutionfield, field in list(md.outputdefinition.__dict__.items()):
    842887                if isinstance(field, list):
    843                     #get each definition
     888                    # Get each definition
    844889                    for i, fieldi in enumerate(field):
    845890                        if fieldi:
    846891                            fieldr = getattr(md.outputdefinition, solutionfield)[i]
    847                             #get subfields
     892                            # Get subfields
    848893                            for solutionsubfield, subfield in list(fieldi.__dict__.items()):
    849894                                if np.size(subfield) == md.mesh.numberofvertices:
    … …  
    852897                                    setattr(fieldr, solutionsubfield, project2d(md, subfield, 1))
    853898
    854         #Initialize with the 2d mesh
     899        # Initialize he 2d mesh
    855900        mesh = mesh2d()
    856901        mesh.x = md.mesh.x2d
    … …  
    859904        mesh.numberofelements = md.mesh.numberofelements2d
    860905        mesh.elements = md.mesh.elements2d
     906        # if not np.isnan(md.mesh.vertexonboundary).all():
     907        #     mesh.vertexonboundary = project2d(md, md.mesh.vertexonboundary, 1)
     908        # if not np.isnan(md.mesh.elementconnectivity).all():
     909        #     mesh.elementconnectivity = project2d(md, md.mesh.elementconnectivity, 1)
     910        if np.size(md.mesh.lat) == md.mesh.numberofvertices:
     911            mesh.lat = project2d(md, md.mesh.lat, 1)
     912        if np.size(md.mesh.long) == md.mesh.numberofvertices:
     913            mesh.long = project2d(md, md.mesh.long, 1)
     914        mesh.epsg = md.mesh.epsg
     915        if np.size(md.mesh.scale_factor) == md.mesh.numberofvertices:
     916            mesh.scale_factor = project2d(md, md.mesh.scale_factor, 1)
    861917        if not np.isnan(md.mesh.vertexonboundary).all():
    862             mesh.vertexonboundary = project2d(md, md.mesh.vertexonboundary, 1)
    863         if not np.isnan(md.mesh.elementconnectivity).all():
    864             mesh.elementconnectivity = project2d(md, md.mesh.elementconnectivity, 1)
    865         if isinstance(md.mesh.lat, np.ndarray):
    866             if md.mesh.lat.size == md.mesh.numberofvertices:
    867                 mesh.lat = project2d(md, md.mesh.lat, 1)
    868         if isinstance(md.mesh.long, np.ndarray):
    869             if md.mesh.long.size == md.mesh.numberofvertices:
    870                 md.mesh.long = project2d(md, md.mesh.long, 1)
    871         mesh.epsg = md.mesh.epsg
    872         if isinstance(md.mesh.scale_factor, np.ndarray):
    873             if md.mesh.scale_factor.size == md.mesh.numberofvertices:
    874                 md.mesh.scale_factor = project2d(md, md.mesh.scale_factor, 1)
     918            mesh.vertexonboundary= project2d(md, md.mesh.vertexonboundary, 1)
     919        if not np.isnan(md.mesh.elementonboundary).all():
     920            mesh.elementonboundary = project2d(md, md.mesh.elementonboundary, 1)
    875921        md.mesh = mesh
    876922        md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
    877923        md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
    878         md.mesh.segments = contourenvelope(md)
     924        md.mesh.segments = contourenvelope(md.mesh)
    879925
    880926        return md

  • issm/trunk-jpl/src/m/classes/sealevelmodel.m

    r25455 r25499  
    174174                       
    175175                        %initialize, to avoid issues of having more transitions than meshes.
    176                         self.transitions={}; self.eltransitions={};
     176                        self.transitions={};
     177                        self.eltransitions={};
    177178
    178179                        %for elements:
    … …  
    180181                        ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
    181182                        ze=self.earth.mesh.z(self.earth.mesh.elements)*[1;1;1]/3;
    182 
     183                       
    183184                        for i=1:length(self.icecaps),
    184185                                mdi=self.icecaps{i};
    … …  
    196197                                self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force);
    197198                        end
    198 
    199199                end % }}}
    200200                function checkintersections(self) % {{{
    … …  
    397397                                                eval(['capfieldj= self.icecaps{' num2str(j) '}.' string ';']);
    398398                                                if ~isequal(capfieldi(end,:),capfieldj(end,:)),
    399                                                         error(['Time stamps for ' string ' field is different between icecaps ' num2str(i) ' and ' num2str(j)]);
     399                                                        error(['Time stamps for ' string ' field are different between icecaps ' num2str(i) ' and ' num2str(j)]);
    400400                                                end
    401401                                        end

  • issm/trunk-jpl/src/m/classes/sealevelmodel.py

    r25477 r25499  
     1from copy import deepcopy
    12import warnings
    23import numpy as np
    34from fielddisplay import fielddisplay
    4 from fieldnames import fieldnames
    55from generic import generic
     6from getsubattr import getsubattr
    67from issmsettings import issmsettings
    78from meshintersect3d import meshintersect3d
    … …  
    1011from modelmerge3d import modelmerge3d
    1112from pairoptions import pairoptions
     13from planetradius import planetradius
     14from plotmodel import plotmodel
    1215from private import private
     16from setsubattr import setsubattr
    1317from TwoDToThreeD import TwoDToThreeD
    14 from plotmodel import plotmodel
    15 from plot3 import plot3
    16 from loneedges import loneedges
    17 from planetradius import planetradius
     18
    1819
    1920class sealevelmodel(object):
    … …  
    174175
    175176        # For elements
    176         onesmatrix = np.array([[1], [1], [1]])
    177         xe = self.earth.mesh.x[self.earth.mesh.elements] * onesmatrix / 3
    178         ye = self.earth.mesh.y[self.earth.mesh.elements] * onesmatrix / 3
    179         ze = self.earth.mesh.z[self.earth.mesh.elements] * onesmatrix / 3
     177        xe = np.mean(self.earth.mesh.x[self.earth.mesh.elements - 1], axis=1)
     178        ye = np.mean(self.earth.mesh.y[self.earth.mesh.elements - 1], axis=1)
     179        ze = np.mean(self.earth.mesh.z[self.earth.mesh.elements - 1], axis=1)
    180180
    181181        for i in range(len(self.icecaps)):
    … …  
    184184
    185185            # For elements
    186             xei = mdi.mesh.x[mdi.mesh.elements] * onesmatrix / 3
    187             yei = mdi.mesh.y[mdi.mesh.elements] * onesmatrix / 3
    188             zei = mdi.mesh.z[mdi.mesh.elements] * onesmatrix / 3
     186            xei = np.mean(mdi.mesh.x[mdi.mesh.elements - 1], axis=1)
     187            yei = np.mean(mdi.mesh.y[mdi.mesh.elements - 1], axis=1)
     188            zei = np.mean(mdi.mesh.z[mdi.mesh.elements - 1], axis=1)
    189189
    190190            print('Computing vertex intersections for basin {}'.format(self.basins[i].name))
    … …  
    300300
    301301        # Make 3D model
    302         models = self.icecaps
     302        models = deepcopy(self.icecaps)
    303303        for i in range(len(models)):
    304304            models[i] = TwoDToThreeD(models[i], self.planet)
    … …  
    308308        for i in range(1, len(models)):
    309309            md = modelmerge3d(md, models[i], 'tolerance', tolerance)
    310             md.private.bamg.landmask = np.vstack((md.private.bamg.landmask, models[i].private.bamg.landmask))
     310            md.private.bamg.landmask = np.hstack((md.private.bamg.landmask, models[i].private.bamg.landmask))
    311311
    312312        # Look for lone edges if asked for it
    … …  
    317317                ind1 = edges(i, 1)
    318318                ind2 = edges(i, 2)
     319                # TODO: Reconfigure the following in the process of bringing plotting online
    319320                plot3([md.mesh.x[ind1], md.mesh.x[ind2]], [md.mesh.y[ind1], md.mesh.y[ind2]], [md.mesh.z[ind1], md.mesh.z[ind2]], 'g*-')
    320321
    … …  
    323324
    324325        # Create mesh radius
    325         self.earth.mesh.r = planetradius('earth')
     326        self.earth.mesh.r = planetradius('earth') * np.ones((md.mesh.numberofvertices, ))
    326327    #}}}
    327328
    … …  
    349350    def transfer(self, string):  #{{{
    350351        # Recover field size in one icecap
    351         n = np.size(getattr(self.icecaps[0], string), 0)
     352        n = getsubattr(self.icecaps[0], string).shape[0]
    352353
    353354        if n == self.icecaps[0].mesh.numberofvertices:
    354             setattr(self.earth, string, np.zeros((self.earth.mesh.numberofvertices, )))
     355            setsubattr(self.earth, string, np.zeros((self.earth.mesh.numberofvertices, ))) # Assign array of zeros to target attribute
     356            earth_attr = getsubattr(self.earth, string) # Retrieve reference to target attribute
    355357            for i in range(len(self.icecaps)):
    356                 getattr(self.earth, string)[self.transitions[i]] = getattr(self.icecaps[i], string)
    357         elif n == (self.self.icecaps[0].mesh.numberofvertices + 1):
     358                earth_attr[self.transitions[i]] = getsubattr(self.icecaps[i], string)
     359        elif n == (self.icecaps[0].mesh.numberofvertices + 1):
    358360            # Dealing with transient dataset
    359361            # Check that all timetags are similar between all icecaps #{{{
    360362            for i in range(len(self.icecaps)):
    361                 capfieldi = getattr(self.icecaps[i], string)
    362                 for j in range(1, len(self.icecaps)):
    363                     capfieldj = getattr(self.icecaps[j], string)
     363                capfieldi = getsubattr(self.icecaps[i], string)
     364                for j in range(i + 1, len(self.icecaps)):
     365                    capfieldj = getsubattr(self.icecaps[j], string)
    364366                    if capfieldi[-1, :] != capfieldj[-1, :]:
    365                         raise Exception("Time stamps for {} field is different between icecaps {} and {}".format(string, i, j))
    366             capfield1 = getattr(self.icecaps[0], string)
     367                        raise Exception("Time stamps for {} field are different between icecaps {} and {}".format(string, i, j))
     368            capfield1 = getsubattr(self.icecaps[0], string)
    367369            times = capfield1[-1, :]
    368370            nsteps = len(times)
    … …  
    374376            # Transfer all the time fields #{{{
    375377            for i in range(len(self.icecaps)):
    376                 capfieldi = getattr(self.icecaps[i], string)
     378                capfieldi = getsubattr(self.icecaps[i], string)
    377379                for j in range(nsteps):
    378                     field[self.transitions[i], j] = capfieldi[0:-2, j]  # Transfer only the values, not the time
    379             setattr(self.earth, string, field)  # Do not forget to plug the field variable into its final location
     380                    field[self.transitions[i], j] = capfieldi[0:-1, j]  # Transfer only the values, not the time
     381            setsubattr(self.earth, string, field)  # Do not forget to plug the field variable into its final location
    380382            #}}}
    381383        elif n == (self.icecaps[0].mesh.numberofelements):
    382             setattr(self.earth, string, np.zeros((self.earth.mesh.numberofvertices, )))
     384            setsubattr(self.earth, string, np.zeros((self.earth.mesh.numberofelements, ))) # Assign array of zeros to target attribute
     385            earth_attr = getsubattr(self.earth, string) # Retrieve reference to target attribute
    383386            for i in range(len(self.icecaps)):
    384                 getattr(self.earth, string)[self.eltransitions[i]] = getattr(self.icecaps[i], string)
     387                earth_attr[self.eltransitions[i]] = getsubattr(self.icecaps[i], string)
    385388        else:
    386389            raise Exception('not supported yet')

  • issm/trunk-jpl/src/m/classes/surfaceload.py

    r25171 r25499  
    77
    88class surfaceload(object):
    9     '''
    10     SURFACELOAD class definition
     9    """SURFACELOAD class definition
    1110
    12         Usage:
    13             surfaceload = surfaceload()
    14     '''
     11    Usage:
     12        surfaceload = surfaceload()
     13    """
    1514
    1615    def __init__(self, *args): #{{{
    … …  
    5857    def marshall(self, prefix, md, fid): #{{{
    5958        if len(self.icethicknesschange) == 0:
    60             self.icethicknesschange = np.zeros(md.mesh.numberofelements + 1)
     59            self.icethicknesschange = np.zeros((md.mesh.numberofelements + 1, ))
    6160
    6261        if len(self.waterheightchange) == 0:
    63             self.waterheightchange = np.zeros(md.mesh.numberofelements + 1)
     62            self.waterheightchange = np.zeros((md.mesh.numberofelements + 1, ))
    6463
    6564        if len(self.other) == 0:
    66             self.other = np.zeros(md.mesh.numberofelements + 1)
     65            self.other = np.zeros((md.mesh.numberofelements + 1, ))
    6766
    6867        WriteData(fid, prefix, 'object', self, 'fieldname', 'icethicknesschange', 'name', 'md.solidearth.surfaceload.icethicknesschange', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', md.constants.yts)

  • issm/trunk-jpl/src/m/consistency/checkfield.py

    r25385 r25499  
    88
    99def checkfield(md, *args):
    10     '''
    11     CHECKFIELD - check field consistency
    12 
    13         Used to check model consistency
     10    """CHECKFIELD - check field consistency
     11
     12    Used to check model consistency
    1413       
    15         Requires:
    16             'field' or 'fieldname' option. If 'fieldname' is provided, it will retrieve it from the model md. (md.(fieldname))
    17             If 'field' is provided, it will assume the argument following
    18             'field' is a numeric array.
    19 
    20         Available options:
     14    Requires:
     15        'field' or 'fieldname' option. If 'fieldname' is provided, it will
     16        retrieve it from the model md. (md.(fieldname))
     17        If 'field' is provided, it will assume the argument following 'field'
     18        is a numeric array.
     19
     20    Available options:
    2121        - NaN: 1 if check that there is no NaN
    2222        - size: [lines cols], NaN for non checked dimensions, or 'universal'
    … …  
    3434        - message: overloaded error message
    3535
    36         Usage:
    37             md = checkfield(md, fieldname, options)
    38     '''
     36    Usage:
     37        md = checkfield(md, fieldname, options)
     38    """
    3939
    4040    #get options
    … …  
    7676                #Check that vector size will not be confusing for ModelProcessorx
    7777                if (md.mesh.numberofvertices == md.mesh.numberofelements):
    78                     raise RuntimeError('number of vertices is the same as number of elements')
     78                    raise Exception('number of vertices is the same as number of elements')
    7979                elif (md.mesh.numberofvertices + 1 == md.mesh.numberofelements):
    80                     raise RuntimeError('number of vertices + 1 is the same as number of elements')
     80                    raise Exception('number of vertices + 1 is the same as number of elements')
    8181                elif (md.mesh.numberofvertices == md.mesh.numberofelements + 1):
    82                     raise RuntimeError('number of vertices is the same as number of elements + 1')
     82                    raise Exception('number of vertices is the same as number of elements + 1')
    8383
    8484                #Uniform field
    85                 if (np.size(field, 0) == 1):
     85                if (field.shape[0] == 1):
    8686                    if (np.ndim(field) > 1 and np.shape(field)[1] != 1):
    8787                        md = md.checkmessage(options.getfieldvalue('message', "field '{}' is not supported".format(fieldname)))
    … …  
    112112
    113113            else:
    114                 raise RuntimeError("fieldsize '{}' not supported yet".format(fieldsize))
     114                raise Exception("fieldsize '{}' not supported yet".format(fieldsize))
    115115
    116116        else:
    … …  
    155155                md = md.checkmessage(options.getfieldvalue('message', "field '{}' value should be '{}'".format(fieldname, fieldvalues[0])))
    156156            elif len(fieldvalues) == 2:
    157                 md = md.checkmessage(options.getfieldvalue('message', "field '{}' values should be '{}' '{}' or '%s'".format(fieldname, fieldvalues[0], fieldvalues[1])))
     157                md = md.checkmessage(options.getfieldvalue('message', "field '{}' values should be '{}' '{}' or '{}'".format(fieldname, fieldvalues[0], fieldvalues[1])))
    158158            else:
    159159                md = md.checkmessage(options.getfieldvalue('message', "field '{}' should have values in {}".format(fieldname, fieldvalues)))
    … …  
    166166        if np.size(lowerbound) > 1:  #checking elementwise
    167167            if any(field < lowerbound):
    168                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values above %d" % (fieldname, lowerbound)))
     168                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
    169169        else:
    170170            minval = np.nanmin(field)
    … …  
    178178
    179179            if minval < lowerbound:
    180                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values above %d" % (fieldname, lowerbound)))
     180                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
    181181
    182182    if options.exist('>'):
    … …  
    186186        if np.size(lowerbound) > 1:  #checking elementwise
    187187            if any(field <= lowerbound):
    188                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values above %d" % (fieldname, lowerbound)))
     188                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
    189189        else:
    190190            minval = np.nanmin(field)
    … …  
    198198
    199199            if minval <= lowerbound:
    200                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values above %d" % (fieldname, lowerbound)))
     200                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
    201201
    202202    #check smaller
    … …  
    207207        if np.size(upperbound) > 1:  #checking elementwise
    208208            if any(field > upperbound):
    209                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values below %d" % (fieldname, upperbound)))
     209                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
    210210        else:
    211211            maxval = np.nanmax(field)
    … …  
    220220                maxval = field.fov_forward_indices[0]
    221221            if maxval > upperbound:
    222                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values below %d" % (fieldname, upperbound)))
     222                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
    223223
    224224    if options.exist('<'):
    … …  
    228228        if np.size(upperbound) > 1:  #checking elementwise
    229229            if any(field >= upperbound):
    230                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values below %d" % (fieldname, upperbound)))
     230                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
    231231
    232232        else:
    … …  
    241241
    242242                if maxval >= upperbound:
    243                     md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values below %d" % (fieldname, upperbound)))
     243                    md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
    244244
    245245    #check file
    246246    if options.getfieldvalue('file', 0):
    247247        if not os.path.exists(field):
    248             md = md.checkmessage("file provided in '%s': '%s' does not exist" % (fieldname, field))
     248            md = md.checkmessage("file provided in {}: {} does not exist".format(fieldname, field))
    249249
    250250    #Check row of strings
    251251    if options.exist('stringrow'):
    252252        if not isinstance(field, list):
    253             md = md.checkmessage(options.getfieldvalue('message', "field '%s' should be a list" % fieldname))
     253            md = md.checkmessage(options.getfieldvalue('message', "field {} should be a list".format(fieldname)))
    254254
    255255    #Check forcings (size and times)
    256256    if options.getfieldvalue('timeseries', 0):
    257         if np.size(field, 0) == md.mesh.numberofvertices or np.size(field, 0) == md.mesh.numberofelements:
     257        if field.shape[0] == md.mesh.numberofvertices or field.shape[0] == md.mesh.numberofelements:
    258258            if np.ndim(field) > 1 and not np.size(field, 1) == 1:
    259                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname))
    260         elif np.size(field, 0) == md.mesh.numberofvertices + 1 or np.size(field, 0) == md.mesh.numberofelements + 1:
     259                md = md.checkmessage(options.getfieldvalue('message', "field {} should have only one column as there are md.mesh.numberofvertices lines".format(fieldname)))
     260        elif field.shape[0] == md.mesh.numberofvertices + 1 or field.shape[0] == md.mesh.numberofelements + 1:
    261261            if np.ndim(field) > 1 and not all(field[-1, :] == np.sort(field[-1, :])):
    262                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' columns should be sorted chronologically" % fieldname))
     262                md = md.checkmessage(options.getfieldvalue('message', "field {} columns should be sorted chronologically".format(fieldname)))
    263263            if np.ndim(field) > 1 and any(field[-1, 0:-1] == field[-1, 1:]):
    264                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' columns must not contain duplicate timesteps" % fieldname))
    265         else:
    266             md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have md.mesh.numberofvertices or md.mesh.numberofvertices + 1 lines" % fieldname))
     264                md = md.checkmessage(options.getfieldvalue('message', "field {} columns must not contain duplicate timesteps".format(fieldname)))
     265        else:
     266            md = md.checkmessage(options.getfieldvalue('message', "field {} should have md.mesh.numberofvertices or md.mesh.numberofvertices + 1 lines".format(fieldname)))
    267267
    268268    #Check single value forcings (size and times)
    269269    if options.getfieldvalue('singletimeseries', 0):
    270         if np.size(field, 0) == 2:
     270        if field.shape[0] == 2:
    271271            if not all(field[-1, :] == np.sort(field[-1, :])):
    272                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' columns should be sorted chronologically" % fieldname))
     272                md = md.checkmessage(options.getfieldvalue('message', "field {} columns should be sorted chronologically".format(fieldname)))
    273273            if any(field[-1, 0:-1] == field[-1, 1:]):
    274                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' columns must not contain duplicate timesteps" % fieldname))
    275         elif np.size(field, 0) == 1:
     274                md = md.checkmessage(options.getfieldvalue('message', "field {} columns must not contain duplicate timesteps".format(fieldname)))
     275        elif field.shape[0] == 1:
    276276            if np.ndim(field) > 1 and not np.size(field, 1) == 1:
    277                 md = md.checkmessage(options.getfieldvalue('message', "field '%s' should be either a scalar or have 2 lines" % fieldname))
    278         else:
    279             md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have 2 lines or be a scalar" % fieldname))
     277                md = md.checkmessage(options.getfieldvalue('message', "field {} should be either a scalar or have 2 lines".format(fieldname)))
     278        else:
     279            md = md.checkmessage(options.getfieldvalue('message', "field {} should have 2 lines or be a scalar".format(fieldname)))
    280280
    281281    return md

  • issm/trunk-jpl/src/m/geometry/GetAreas3DTria.m

    r19311 r25499  
    22%GETAREAS3DTRIA - compute areas of triangles with 3D coordinates
    33%
    4 %   compute areas of trianguls with 3D coordinates
     4%   Compute areas of triangles with 3D coordinates
    55%
    66%   Usage:
    7 %      areas  =GetAreas3DTria(index,x,y,z);
     7%      areas=GetAreas3DTria(index,x,y,z);
    88%
    99%   Examples:
    10 %      areas  =GetAreas3DTria(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z);
     10%      areas=GetAreas3DTria(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z);
    1111
    1212%get number of elements and number of nodes
    13 nels=size(index,1); 
    14 nods=length(x); 
     13nels=size(index,1);
     14nods=length(x);
    1515
    1616%some checks
    1717if nargout~=1 | (nargin~=3 & nargin~=4),
    1818        help GetAreas3DTria
    19         error('GetAreas error message: bad usage')
     19        error('GetAreas3DTria error message: bad usage')
    2020end
    2121if ((length(y)~=nods) | (nargin==4 & length(z)~=nods)),
    22         error('GetAreas3DTria error message: x,y and z do not have the same length')
     22        error('GetAreas3DTria error message: x, y, and z do not have the same length')
    2323end
     24
    2425if max(index(:))>nods,
    2526        error(['GetAreas3DTria error message: index should not have values above ' num2str(nods) ])
    2627end
    2728if (nargin==4 & size(index,2)~=3),
    28         error('GetAreas3DTria error message: index should have 3 columns for 2d meshes.')
     29        error('GetAreas3DTria error message: index should have 3 columns for 2d meshes')
    2930end
    3031
    … …  
    3738%compute the volume of each element
    3839if nargin==4,
    39    % area of triangles with 3D coordinats
    40    for j=1:nels
    41                 m1=[x1(j) x2(j) x3(j); y1(j) y2(j) y3(j); 1 1 1];
    42            m2=[y1(j) y2(j) y3(j); z1(j) z2(j) z3(j); 1 1 1];
    43            m3=[z1(j) z2(j) z3(j); x1(j) x2(j) x3(j); 1 1 1];
    44       areas(j)=sqrt(det(m1)^2 + det(m2)^2 + det(m3)^2)/2;
    45    end
    46 end
    47 
    48 
     40        % area of triangles with 3D coordinates
     41        for i=1:nels
     42                m1=[x1(i) x2(i) x3(i); y1(i) y2(i) y3(i); 1 1 1];
     43                m2=[y1(i) y2(i) y3(i); z1(i) z2(i) z3(i); 1 1 1];
     44                m3=[z1(i) z2(i) z3(i); x1(i) x2(i) x3(i); 1 1 1];
     45                areas(i)=sqrt(det(m1)^2 + det(m2)^2 + det(m3)^2)/2;
     46        end
     47end

  • issm/trunk-jpl/src/m/interp/averaging.m

    r17686 r25499  
    4343end
    4444
    45 %load some variables (it is much faster if the variabes are loaded from md once for all)
     45%load some variables (it is much faster if the variables are loaded from md once for all)
    4646if layer==0,
    4747        index=md.mesh.elements;
    … …  
    8787end
    8888
    89 %return output as a full matrix (C code do not like sparse matrices)
     89%return output as a full matrix (C code does not like sparse matrices)
    9090average=full(average_node);

  • issm/trunk-jpl/src/m/interp/averaging.py

    r24256 r25499  
    11import numpy as np
    2 from GetAreas import GetAreas
    32try:
    43    from scipy.sparse import csc_matrix
    … …  
    65    print("could not import scipy, no averaging capabilities enabled")
    76
     7from GetAreas import GetAreas
     8
    89
    910def averaging(md, data, iterations, layer=0):
    10     '''
    11     AVERAGING - smooths the input over the mesh
     11    """AVERAGING - smooths the input over the mesh
    1212
    13        This routine takes a list over the elements or the nodes in input
    14        and return a list over the nodes.
    15        For each iterations it computes the average over each element (average
    16        of the vertices values) and then computes the average over each node
    17        by taking the average of the element around a node weighted by the
    18        elements volume
    19        For 3d mesh, a last argument can be added to specify the layer to be averaged on.
     13    This routine takes a list of the elements or the nodes in input and return
     14    a list over the nodes.
     15    For each iterations it computes the average over each element (average of
     16    the vertices values) and then computes the average over each node by taking
     17    the average of the element around a node weighted by the elements volume.
     18    For 3d mesh, a last argument can be added to specify the layer to be
     19    averaged on.
    2020
    21        Usage:
    22           smoothdata = averaging(md, data, iterations)
    23           smoothdata = averaging(md, data, iterations, layer)
     21    Usage:
     22        smoothdata = averaging(md, data, iterations)
     23        smoothdata = averaging(md, data, iterations, layer)
    2424
    25        Examples:
    26           velsmoothed = averaging(md, md.initialization.vel, 4)
    27           pressure = averaging(md, md.initialization.pressure, 0)
    28           temperature = averaging(md, md.initialization.temperature, 1, 1)
    29     '''
     25    Examples:
     26        velsmoothed = averaging(md, md.initialization.vel, 4)
     27        pressure = averaging(md, md.initialization.pressure, 0)
     28        temperature = averaging(md, md.initialization.temperature, 1, 1)
     29    """
    3030
    3131    if len(data) != md.mesh.numberofelements and len(data) != md.mesh.numberofvertices:
    … …  
    3737        layer = 0
    3838
    39     #initialization
     39    # Initialization
    4040    if layer == 0:
    4141        weights = np.zeros(md.mesh.numberofvertices, )
    … …  
    4545        data = data[(layer - 1) * md.mesh.numberofvertices2d + 1:layer * md.mesh.numberofvertices2d, :]
    4646
    47     #load some variables (it is much faster if the variabes are loaded from md once for all)
     47    # Load some variables (it is much faster if the variables are loaded from md once for all)
    4848    if layer == 0:
    4949        index = md.mesh.elements
    … …  
    5555        numberofelements = md.mesh.numberofelements2d
    5656
    57     #build some variables
     57    # Build some variables
    5858    if md.mesh.dimension() == 3 and layer == 0:
    5959        rep = 6
    … …  
    6666        areas = GetAreas(index, md.mesh.x2d, md.mesh.y2d)
    6767
    68     index = index - 1  # since python indexes starting from zero
     68    index = index - 1  # Python indexes from zero
    6969    line = index.flatten(1)
    7070    areas = np.vstack(areas).reshape(-1, )
    … …  
    7272    linesize = rep * numberofelements
    7373
    74     #update weights that holds the volume of all the element holding the node i
     74    # Update weights that holds the volume of all the element holding the node i
    7575    weights = csc_matrix((np.tile(areas, (rep, 1)).reshape(-1, ), (line, np.zeros(linesize, ))), shape=(numberofnodes, 1))
    7676
    77     #initialization
     77    # Initialization
    7878    if len(data) == numberofelements:
    7979        average_node = csc_matrix((np.tile(areas * data, (rep, 1)).reshape(-1, ), (line, np.zeros(linesize, ))), shape=(numberofnodes, 1))
    … …  
    8383        average_node = csc_matrix(data.reshape(-1, 1))
    8484
    85     #loop over iteration
     85    # Loop over iteration
    8686    for i in np.arange(1, iterations + 1):
    8787        average_el = np.asarray(np.dot(average_node.todense()[index].reshape(numberofelements, rep), np.vstack(summation))).reshape(-1, )
    … …  
    9090        average_node = csc_matrix(average_node)
    9191
    92     #return output as a full matrix (C code do not like sparse matrices)
     92    # Return output as a full matrix (C code does not like sparse matrices)
    9393    average = np.asarray(average_node.todense()).reshape(-1, )
    9494

  • issm/trunk-jpl/src/m/mesh/TwoDToThreeD.m

    r25455 r25499  
    11function md=TwoDToThreeD(md,planet)
    2 
    32        %reproject model into lat,long if necessary:
    43        if ~strcmpi(md.mesh.proj,epsg2proj(4326)),

  • issm/trunk-jpl/src/m/mesh/findsegments.m

    r25457 r25499  
    6666        %'el1' is connected to only one element
    6767        else
    68 
    69                 %find the vertex the 'el1' does not share with 'els2'
     68                %find the vertex that 'el1' does not share with 'els2'
    7069                flag=setdiff(nods1,md.mesh.elements(els2,:));
    7170
    7271                for j=1:3,
    73                         nods=nods1; nods(j)=[];
     72                        nods=nods1;
     73                        nods(j)=[];
    7474                        if any(ismember(flag,nods)),
    7575

  • issm/trunk-jpl/src/m/mesh/findsegments.py

    r25455 r25499  
    3535    pos = np.nonzero(elementonboundary)[0]
    3636    num_segments = len(pos)
    37     segments = np.zeros((num_segments, 3))
     37    segments = np.zeros((num_segments, 3)).astype(int)
    3838    count = 0
    3939
    … …  
    5656
    5757            # Get the vertices on the boundary and build segment
    58             nods1 = np.delete(nods1, np.where(nods1 == flag)[0])
     58            nods1 = np.delete(nods1, np.where(np.in1d(nods1, flag, assume_unique=True)))
    5959            segments[count, :] = np.append(nods1, el1 + 1)
    6060
    … …  
    7272        # 'el1' is connected to only one element
    7373        else:
     74            # NOTE: This block is untested as it does not get touched by test2004 (remove this note once it has been tested)
     75
    7476            # Find the vertex that 'el1' does not share with 'els2'
    75             flag = els2
    76        
    77         exit()
     77            flag = np.setdiff1d(nods, md.mesh.elements[els2, :])
     78
     79            for j in range(3):
     80                nods = nods1
     81                nods = np.delete(nods, j)
     82                if np.any(np.in1d(flag, nods)):
     83                    segments[count, :] = np.append(nods, el1 + 1)
     84
     85                    # Swap segment nodes if necessary
     86                    ord1 = np.where(nods1[0] == md.mesh.elements[el1, :])[0][0]
     87                    ord2 = np.where(nods1[1] == md.mesh.elements[el1, :])[0][0]
     88
     89                    if ((ord1 == 0 and ord2 == 1) or (ord1 == 1 and ord2 == 2) or (ord1 == 2 and ord2 == 0)):
     90                        temp = segments[count, 0]
     91                        segments[count, 0] = segments[count, 1]
     92                        segments[count, 1] = temp
     93
     94                    segments[count, 0:2] = np.flip(segments[count, 0:2]) # NOTE: Upper bound of index range is non-inclusive
     95                    count = count + 1
    7896
    7997    return segments

  • issm/trunk-jpl/src/m/mesh/meshintersect3d.py

    r25065 r25499  
    1 import math
    2 
    31import matplotlib.pyplot as plt
    42from mpl_toolkits.mplot3d import Axes3D
    … …  
    97
    108def meshintersect3d(x, y, z, xs, ys, zs, *args): #{{{
    11     '''
    12     MESHINTERSECT - returns indices (into x, y, and z) of common values between
    13     (x, y, z) and (xs, ys, zs) (i.e. x(index) = xs; y(index) = ys).
    14     '''
     9    """MESHINTERSECT - returns indices (into x, y, and z) of common values
     10    between (x, y, z) and (xs, ys, zs) (i.e. x(index) = xs; y(index) = ys).
     11    """
    1512
    1613    # Process options
    … …  
    1916    # Retrieve tolerance
    2017    maxtol = options.getfieldvalue('maxtol', 100000) # 100 km
    21     tolincrement = options.getfieldvalue('tolincrement', [])
     18    tolincrement = options.getfieldvalue('tolincrement', 10)
    2219    force = options.getfieldvalue('force', 0)
    2320
    24     # Go through lats, longs and find within tolerance, the index of the corresponding value in lat, long
    25     indices = np.zeros(len(xs))
     21    # Go through lats, longs and find, within tolerance, the index of the corresponding value in lat, long
     22    indices = np.zeros((len(xs), ))
    2623
    2724    for i in range(len(xs)):
    2825        tolerance = 0
    29         distance = math.sqrt((x - xs[i]) ** 2 + (y - ys[i]) ** 2 + (z - zs[i]) ** 2)
     26        distance = ((x - xs[i]) ** 2 + (y - ys[i]) ** 2 + (z - zs[i]) ** 2) ** 0.5 # NOTE: math.sqrt cannot be applied element-wise to a list/numpy.array
     27        s = np.where(distance == 0)[0]
    3028
    31         s = np.where(distance == 0)[0]
    32         if s.size():
    33             if s.size() > 1:
     29        if len(s):
     30            if len(s) > 1:
    3431                # We have two vertices that are coincident! Not good
     32                #
     33                # TODO: Reconfigure the following in the process of bringing plotting online
    3534                for j in range(len(s)):
    3635                    plot(x[s[j]], y[s[j]], z[s[j]], c='cyan', s=40)
    … …  
    4645            # We could not find a 0 distance, find the lowest tolerance that generates a find
    4746            count = 1
    48             while not s.size():
     47            while not len(s):
    4948                if count > 1000:
    50                     print('could not find a vertex matching vertex %i of input mesh!' % i)
    51                     print('Might think anbout changing tolerance increment')
    52                     raise RuntimeError('')
     49                    raise Exception('could not find a vertex matching vertex {} of input mesh!\nMight think anbout changing tolerance increment'.format(i))
    5350                tolerance = tolerance + tolincrement
    54                 s = np.where(distance < tolerance)
     51                s = np.where(distance < tolerance)[0]
    5552                count = count + 1
    5653            if tolerance > maxtol:
    57                 print('found matching vertices %i in output mesh for input mesh vertex %i' % (s, i))
    58                 print('however, these vertices are farther than the max tolerance allowed!')
    59                 raise RuntimeError('')
     54                raise Exception('found matching vertices {} in output mesh for input mesh vertex {}\nhowever, these vertices are farther than the max tolerance allowed!'.format(s, i))
    6055
    6156            # Recover minimum distance
    6257            sf = distance[s]
    63             pos = np.where(sf == sf.min())
     58            pos = np.where(sf == np.min(sf))[0]
    6459            s = s[pos]
    6560            indices[i] = s
    6661
    67     if np.where(indices == 0).size():
     62    if len(np.where(indices == 0)[0]) > 1: # NOTE: This check is different than the corresponding one under MATLAB as one index may indeed be '0'
    6863        raise RuntimeError('issue with transition vector having one empty slot')
     64
     65    # Convert results to type 'int' to avoid modifying structures to which they are assigned
     66    indices = indices.astype(int)
     67
     68    return indices
    6969#}}}

  • issm/trunk-jpl/src/m/mesh/modelmerge3d.py

    r25455 r25499  
    112112
    113113    # Vertex on boundary
    114     md.mesh.vertexconnectivity = np.zeros((md.mesh.numberofvertices, ))
    115     md.mesh.vertexonboundary[md.mesh.segments[:, 0:1]] = 1
    116     print(md.mesh.vertexonboundary)
     114    md.mesh.vertexonboundary = np.zeros((md.mesh.numberofvertices, ))
     115    md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1
    117116
    118117    # Some checks
    … …  
    120119        raise Exception('issue in modelmerge, one of the element ids is > number of vertices!')
    121120
    122     exit()
     121    return md
    123122
    124123

  • issm/trunk-jpl/src/m/parameterization/contourenvelope.m

    r19955 r25499  
    2727                isfile=0;
    2828        else
    29                 error('contourenvelope error message:  second argument should be a file or an elements flag');
     29                error('contourenvelope error message: second argument should be a file or an elements flag');
    3030        end
    3131end

  • issm/trunk-jpl/src/m/parameterization/contourenvelope.py

    r25455 r25499  
     1import copy
    12import os.path
    23import numpy as np
    3 import copy
     4
     5from ContourToMesh import ContourToMesh
     6from ElementConnectivity import ElementConnectivity
     7import MatlabFuncs as m
    48from NodeConnectivity import NodeConnectivity
    5 from ElementConnectivity import ElementConnectivity
    6 from ContourToMesh import ContourToMesh
    7 import MatlabFuncs as m
    89
    910
    10 def contourenvelope(md, *args):
     11def contourenvelope(mh, *args):
    1112    """CONTOURENVELOPE - build a set of segments enveloping a contour .exp
    1213
    1314    Usage:
    14         segments = contourenvelope(md, varargin)
     15        segments = contourenvelope(mh, *args)
    1516
    1617    Example:
    17         segments = contourenvelope(md, 'Stream.exp')
    18         segments = contourenvelope(md)
     18        segments = contourenvelope(mh, 'Stream.exp')
     19        segments = contourenvelope(mh)
    1920    """
    2021
    21     #some checks
    22     if len(args) > 1:
    23         raise RuntimeError("contourenvelope error message: bad usage")
     22    # Some checks
     23    nargs = len(args)
    2424
    25     if len(args) == 1:
     25    if nargs > 1:
     26        print(contourenvelope.__doc__)
     27        raise Exception("contourenvelope error message: bad usage")
     28
     29    if nargs == 1:
    2630        flags = args[0]
    2731
    … …  
    3539            isfile = 0
    3640        else:
    37             raise TypeError("contourenvelope error message:  second argument should be a file or an elements flag")
     41            raise TypeError("contourenvelope error message: second argument should be a file or an elements flag")
    3842
    39     #Now, build the connectivity tables for this mesh.
    40     #Computing connectivity
    41     if np.size(md.mesh.vertexconnectivity, axis=0) != md.mesh.numberofvertices and np.size(md.mesh.vertexconnectivity, axis=0) != md.mesh.numberofvertices2d:
    42         md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
    43     if np.size(md.mesh.elementconnectivity, axis=0) != md.mesh.numberofelements and np.size(md.mesh.elementconnectivity, axis=0) != md.mesh.numberofelements2d:
    44         md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
     43    # Now, build the connectivity tables for this mesh
     44    # Computing connectivity
     45    if np.size(mh.vertexconnectivity, axis=0) != mh.numberofvertices and np.size(mh.vertexconnectivity, axis=0) != mh.numberofvertices2d:
     46        mh.vertexconnectivity = NodeConnectivity(mh.elements, mh.numberofvertices)
     47    if np.size(mh.elementconnectivity, axis=0) != mh.numberofelements and np.size(mh.elementconnectivity, axis=0) != mh.numberofelements2d:
     48        mh.elementconnectivity = ElementConnectivity(mh.elements, mh.vertexconnectivity)
    4549
    4650    #get nodes inside profile
    47     elementconnectivity = copy.deepcopy(md.mesh.elementconnectivity)
    48     if md.mesh.dimension() == 2:
    49         elements = copy.deepcopy(md.mesh.elements)
    50         x = copy.deepcopy(md.mesh.x)
    51         y = copy.deepcopy(md.mesh.y)
    52         numberofvertices = copy.deepcopy(md.mesh.numberofvertices)
    53         numberofelements = copy.deepcopy(md.mesh.numberofelements)
     51    elementconnectivity = copy.deepcopy(mh.elementconnectivity)
     52    if mh.dimension() == 2:
     53        elements = copy.deepcopy(mh.elements)
     54        x = copy.deepcopy(mh.x)
     55        y = copy.deepcopy(mh.y)
     56        numberofvertices = copy.deepcopy(mh.numberofvertices)
     57        numberofelements = copy.deepcopy(mh.numberofelements)
    5458    else:
    55         elements = copy.deepcopy(md.mesh.elements2d)
    56         x = copy.deepcopy(md.mesh.x2d)
    57         y = copy.deepcopy(md.mesh.y2d)
    58         numberofvertices = copy.deepcopy(md.mesh.numberofvertices2d)
    59         numberofelements = copy.deepcopy(md.mesh.numberofelements2d)
     59        elements = copy.deepcopy(mh.elements2d)
     60        x = copy.deepcopy(mh.x2d)
     61        y = copy.deepcopy(mh.y2d)
     62        numberofvertices = copy.deepcopy(mh.numberofvertices2d)
     63        numberofelements = copy.deepcopy(mh.numberofelements2d)
    6064
    6165    if len(args) == 1:
    62 
    6366        if isfile:
    64             #get flag list of elements and nodes inside the contour
     67            # Get flag list of elements and nodes inside the contour
    6568            nodein = ContourToMesh(elements, x, y, file, 'node', 1)
    6669            elemin = (np.sum(nodein(elements), axis=1) == np.size(elements, axis=1))
    67     #modify element connectivity
     70            # Modify element connectivity
    6871            elemout = np.nonzero(np.logical_not(elemin))[0]
    6972            elementconnectivity[elemout, :] = 0
    7073            elementconnectivity[np.nonzero(m.ismember(elementconnectivity, elemout + 1))] = 0
    7174        else:
    72             #get flag list of elements and nodes inside the contour
     75            # Get flag list of elements and nodes inside the contour
    7376            nodein = np.zeros(numberofvertices)
    7477            elemin = np.zeros(numberofelements)
    … …  
    7881            nodein[elements[pos, :] - 1] = 1
    7982
    80     #modify element connectivity
     83            # Modify element connectivity
    8184            elemout = np.nonzero(np.logical_not(elemin))[0]
    8285            elementconnectivity[elemout, :] = 0
    8386            elementconnectivity[np.nonzero(m.ismember(elementconnectivity, elemout + 1))] = 0
    8487
    85     #Find element on boundary
    86     #First: find elements on the boundary of the domain
     88    # Find element on boundary
     89    # First: find elements on the boundary of the domain
    8790    flag = copy.deepcopy(elementconnectivity)
    8891    if len(args) == 1:
    … …  
    9093    elementonboundary = np.logical_and(np.prod(flag, axis=1) == 0, np.sum(flag, axis=1) > 0)
    9194
    92     #Find segments on boundary
     95    # Find segments on boundary
    9396    pos = np.nonzero(elementonboundary)[0]
    9497    num_segments = np.size(pos)

  • issm/trunk-jpl/src/m/solve/marshall.m

    r25315 r25499  
    4545% % Uncomment the following to make a copy of the binary input file for debugging
    4646% % purposes (can be fed into scripts/ReadBin.py).
    47 % copyfile([md.miscellaneous.name '.bin'], [md.miscellaneous.name '.m.bin'])
     47copyfile([md.miscellaneous.name '.bin'], [md.miscellaneous.name '.m.bin'])
    4848if st==-1,
    4949        error(['marshall error message: could not close file ' [md.miscellaneous.name '.bin']]);

  • issm/trunk-jpl/src/m/solve/marshall.py

    r25315 r25499  
    4444        # # Uncomment the following to make a copy of the binary input file for
    4545        # # debugging purposes (can be fed into scripts/ReadBin.py).
    46         # copy_cmd = "cp {}.bin {}.py.bin".format(md.miscellaneous.name, md.miscellaneous.name)
    47         # subprocess.call(copy_cmd, shell=True)
     46        copy_cmd = "cp {}.bin {}.py.bin".format(md.miscellaneous.name, md.miscellaneous.name)
     47        subprocess.call(copy_cmd, shell=True)
    4848    except IOError as e:
    4949        print("marshall error message: could not close file '{}.bin' due to:".format(md.miscellaneous.name), e)

  • issm/trunk-jpl/src/m/solve/solve.py

    r25300 r25499  
    1212
    1313def solve(md, solutionstring, *args):
    14     '''
    15     SOLVE - apply solution sequence for this model
     14    """SOLVE - apply solution sequence for this model
    1615
    17         Usage:
    18             md = solve(md, solutionstring, varargin)
    19             where varargin is a list of paired arguments of string OR enums
     16    Usage:
     17        md = solve(md, solutionstring, varargin)
     18        where varargin is a list of paired arguments of string OR enums
    2019
    2120        solution types available comprise:
    … …  
    4342                                  directory) where the restart file is located
    4443
    45         Examples:
    46             md = solve(md, 'Stressbalance')
    47             md = solve(md, 'sb')
    48     '''
     44    Examples:
     45        md = solve(md, 'Stressbalance')
     46        md = solve(md, 'sb')
     47    """
    4948
    5049    #recover and process solve options

  • issm/trunk-jpl/test/NightlyRun/test2004.m

    r25455 r25499  
    9393
    9494        %miscellaneous:
    95         md.mesh.proj=bas.proj; md.miscellaneous.name=bas.name;
     95        md.mesh.proj=bas.proj;
     96        md.miscellaneous.name=bas.name;
    9697
    9798        %recover mask where we have land:
    … …  
    324325sl.caticecaps('tolerance',tolerance,'loneedgesdetect',loneedgesdetect);
    325326
    326 %figure out how each icecap's mesh connects to the larger earth mesh:
     327%figure out how each icecap's mesh connects to the larger Earth mesh:
    327328sl.intersections('force',1);
    328329
Note: See TracChangeset for help on using the changeset viewer.