Changeset 25455


Ignore:
Timestamp:
08/25/20 00:32:13 (5 years ago)
Author:
jdquinn
Message:

CHG: Saving chnages so that Basile has access to potential fix to solidearthmodel class

Location:
issm/trunk-jpl
Files:
9 added
68 edited

Legend:

Unmodified
Added
Removed

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

    r24213 r25455  
    33
    44class bamgmesh(object):
    5     """
    6     BAMGMESH class definition
     5    """BAMGMESH class definition
    76
    8        Usage:
    9           bamgmesh(varargin)
     7    Usage:
     8        bamgmesh(varargin)
    109    """
    1110

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

    r25065 r25455  
    117117                        %recover options
    118118                        options=pairoptions(varargin{:});
    119                         x=[]; y=[];
     119                        x=[];
     120                        y=[];
    120121
    121122                        %go through boundaries, recover edges and project them in the basin epsg reference frame:
    … …  
    127128                                y=[y;contour.y];
    128129                        end
     130
    129131                        %close the contour:
    130132                        if x(end)~=x(1) | y(end)~=y(1),
    … …  
    181183                                error('Basin shapefilecrop error message: epsgshapefile or projshapefile should be specified');
    182184                        end
    183 
    184185
    185186                        %create list of contours that have critical length > threshold (in lat,long):
    … …  
    209210                        for i=1:length(contours),
    210211                                h=contours(i);
    211                                 isin=inpolygon(h.x,h.y,ba.x,ba.y); 
     212                                isin=inpolygon(h.x,h.y,ba.x,ba.y);
    212213                                if ~isempty(find(isin==0)),
    213214                                        flags(i)=1;
    214215                                end
    215216                        end
    216                         pos=find(flags);  contours(pos)=[];
     217
     218                        pos=find(flags);
     219                        contours(pos)=[];
    217220
    218221                        %Two options:

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

    r25125 r25455  
     1from collections import OrderedDict
    12import math
     3
     4from bamg import bamg
     5import matplotlib.path as path
    26import numpy as np
    37
    48from boundary import boundary
    59from epsg2proj import epsg2proj
     10from fielddisplay import fielddisplay
    611from gdaltransform import gdaltransform
    712from helpers import fileparts
     13from inpolygon import inpolygon
    814from laea import laea
    915from pairoptions import pairoptions
    … …  
    1319
    1420class basin(object): #{{{
    15     '''
    16     BASIN class definition
     21    """BASIN class definition
    1722
    1823    Usage:
    1924        basin = basin()
    20     '''
     25    """
    2126
    2227    def __init__(self, *args): #{{{
    … …  
    136141        for i in range(len(self.boundaries)):
    137142            boundary = self.boundaries[i]
    138             contour = boundary.edges()
    139             contour.x, contour.y = gdaltransform(contour.x, contour.y, boundary.proj, self.proj)
    140             x = [[x], [contour.x]]
    141             y = [[y], [contour.y]]
     143            contours = boundary.edges() # NOTE: Differs from MATLAB in that shpread returns a list
     144            for j in range(len(contours)):
     145                contours[j]['x'], contours[j]['y'] = gdaltransform(contours[j]['x'], contours[j]['y'], boundary.proj, self.proj)
     146                if x == []:
     147                    x = np.array(contours[j]['x'])
     148                    y = np.array(contours[j]['y'])
     149                else:
     150                    x = np.concatenate((x, contours[j]['x']), axis=0)
     151                    y = np.concatenate((y, contours[j]['y']), axis=0)
    142152
    143153        #close the contour
    144154        if x[-1] != x[0] or y[-1] != y[0]:
    145             x.append(x[0])
    146             y.append(y[0])
    147 
    148         out = {}
     155            x = np.append(x, x[0])
     156            y = np.append(y, y[0])
     157
     158        out = OrderedDict()
    149159        out['x'] = x
    150160        out['y'] = y
    … …  
    182192        threshold = options.getfieldvalue('threshold', .65) #.65 degrees lat, long
    183193        inshapefile = options.getfieldvalue('shapefile')
    184         outshapefile = options.getfieldvalue('outshapefile', '')
     194        outputshapefile = options.getfieldvalue('outputshapefile', '')
    185195
    186196        if options.exist('epsgshapefile') and options.exist('projshapefile'):
    … …  
    199209        for i in range(len(contours)):
    200210            contour = contours[i]
    201             carea = polyarea(contour.x, contour.y)
     211            carea = polyarea(contour['x'], contour['y'])
    202212            clength = math.sqrt(carea)
    203213            if clength < threshold:
    204                 llist = np.concatenate(llist, i)
    205         setattr(contours, llist, [])
     214                llist.append(i)
     215
     216        contours_above_threshold = []
     217        for i in range(len(contours)):
     218            if i not in llist:
     219                contours_above_threshold.append(contours[i])
     220        contours = contours_above_threshold
    206221
    207222        #project onto reference frame
    208223        for i in range(len(contours)):
    209224            h = contours[i]
    210             h.x, h.y = gdaltransform(h.x, h.y, projshapefile, self.proj)
    211             contours[i].x = h.x
    212             contours[i].y = h.y
     225            h['x'], h['y'] = gdaltransform(h['x'], h['y'], projshapefile, self.proj)
     226            contours[i]['x'] = h['x']
     227            contours[i]['y'] = h['y']
    213228
    214229        #only keep the contours that are inside the basin (take centroids)
    … …  
    217232        for i in range(len(contours)):
    218233            h = contours[i]
    219             isin = inpolygon(h.x, h.y, ba.x, ba.y)
    220             if len(np.where(isin == 0, isin, isin)):
     234            isin = inpolygon(h['x'], h['y'], ba['x'], ba['y'])
     235            if len(np.where(isin == 0)[0]):
    221236                flags[i] = 1
    222         pos = flags.nonzero()
    223         contours[pos] = []
     237
     238        pos = flags.nonzero()[0]
     239
     240        contours_in_basin = []
     241        for i in range(len(contours)):
     242            if i not in pos:
     243                contours_in_basin.append(contours[i])
     244        contours = contours_in_basin
    224245
    225246        #Two options
     247        output = None
    226248        if outputshapefile == '':
    227249            output = contours
    228250        else:
    229251            shpwrite(contours, outputshapefile)
     252
     253        return output
    230254    #}}}
    231255#}}}

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

    r25125 r25455  
    44
    55from epsg2proj import epsg2proj
     6from fielddisplay import fielddisplay
    67from helpers import *
    78from pairoptions import pairoptions
    … …  
    1011
    1112class boundary(object): #{{{
    12     '''
    13     BOUNDARY class definition
     13    """BOUNDARY class definition
    1414
    1515    Usage:
    1616        boundary = boundary()
    17     '''
     17    """
    1818
    1919    def __init__(self, *args): #{{{
    … …  
    3232            self.shpfilename = options.getfieldvalue('shpfilename', '')
    3333            self.orientation = options.getfieldvalue('orientation', 'normal')
     34
     35            # If shppath is missing trailing slash, add it
     36            if self.shppath[-1] != '/':
     37                self.shppath += '/'
    3438
    3539            #figure out projection string:
    … …  
    7781        #read domain
    7882        path, name, ext = fileparts(self.shpfilename)
    79         if ext != 'shp':
    80             ext = 'shp'
    81         output = shpread('{}/{}.{}'.format(self.shppath, name, ext))
     83        if ext != '.shp':
     84            ext = '.shp'
     85        output = shpread('{}{}{}'.format(self.shppath, name, ext))
    8286
    8387        #do we reverse?
    8488        if self.orientation == 'reverse':
    8589            for i in range(len(output)):
    86                 output[i].x = np.flipud(output[i].x)
    87                 output[i].y = np.flipud(output[i].y)
     90                output[i]['x'] = np.flipud(output[i]['x'])
     91                output[i]['y'] = np.flipud(output[i]['y'])
     92
     93        return output
    8894    #}}}
    8995
    … …  
    107113        #read domain
    108114        path, name, ext = fileparts(self.shpfilename)
    109         if ext != 'shp':
    110             ext = 'shp'
    111         domain = shpread('{}/{}.{}'.format(self.shppath, name, ext))
     115        if ext != '.shp':
     116            ext = '.shp'
     117        domain = shpread('{}{}{}'.format(self.shppath, name, ext))
    112118
    113119        #convert boundary to another reference frame #{{{
    114120        for i in range(len(domain)):
    115121            try:
    116                 x, y = gdaltransform(domain[i].x, domain[i].y, self.proj, proj)
     122                x, y = gdaltransform(domain[i]['x'], domain[i]['y'], self.proj, proj)
    117123            except error as e:
    118124                print(e)
    … …  
    170176        #read domain
    171177        path, name, ext = fileparts(self.shpfilename)
    172         if ext != 'shp':
    173             ext = 'shp'
    174         domain = shpread('{}/{}.{}'.format(self.shppath, name, ext))
     178        if ext != '.shp':
     179            ext = '.shp'
     180        domain = shpread('{}{}{}'.format(self.shppath, name, ext))
    175181
    176182        #TODO: Finish translating from MATLAB after test2004.py runs without plot

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

    r25158 r25455  
    99
    1010class mesh3dsurface(object):
    11     '''
    12     MESH3DSURFACE class definition
    13 
    14         Usage:
    15             mesh3dsurface = mesh3dsurface()
    16     '''
    17 
    18     def __init__(self, *args): # {{{
     11    """MESH3DSURFACE class definition
     12
     13    Usage:
     14        mesh3dsurface = mesh3dsurface()
     15    """
     16
     17    def __init__(self, *args): #{{{
    1918        self.x = np.nan
    2019        self.y = np.nan
    … …  
    4039        self.extractedelements = np.nan
    4140
    42         if not len(args):
     41        nargs = len(args)
     42        if not nargs:
    4343            self.setdefaultparameters()
    44         elif len(args) == 1:
     44        elif nargs == 1:
    4545            self = mesh3dsurface()
    4646            arg = args[1]
    … …  
    5252        else:
    5353            raise RuntimeError('constructor not supported')
    54 
    55     # }}}
    56 
    57     def __repr__(self): # {{{
    58         string = '   2D tria Mesh (horizontal):'
    59 
    60         string += '\n      Elements and vertices:'
    61         string = "%s\n%s" % (string, fielddisplay(self, 'numberofelements', 'number of elements'))
    62         string = "%s\n%s" % (string, fielddisplay(self, 'numberofvertices', 'number of vertices'))
    63         string = "%s\n%s" % (string, fielddisplay(self, 'elements', 'vertex indices of the mesh elements'))
    64         string = "%s\n%s" % (string, fielddisplay(self, 'x', 'vertices x coordinate [m]'))
    65         string = "%s\n%s" % (string, fielddisplay(self, 'y', 'vertices y coordinate [m]'))
    66         string = "%s\n%s" % (string, fielddisplay(self, 'z', 'vertices z coordinate [m]'))
    67         string = "%s\n%s" % (string, fielddisplay(self, 'lat', 'vertices latitude [degrees]'))
    68         string = "%s\n%s" % (string, fielddisplay(self, 'long', 'vertices longitude [degrees]'))
    69         string = "%s\n%s" % (string, fielddisplay(self, 'r', 'vertices radius [m]'))
    70 
    71         string = "%s\n%s" % (string, fielddisplay(self, 'edges', 'edges of the 2d mesh (vertex1 vertex2 element1 element2)'))
    72         string = "%s\n%s" % (string, fielddisplay(self, 'numberofedges', 'number of edges of the 2d mesh'))
    73 
    74         string += '\n      Properties:'
    75         string = "%s\n%s" % (string, fielddisplay(self, 'vertexonboundary', 'vertices on the boundary of the domain flag list'))
    76         string = "%s\n%s" % (string, fielddisplay(self, 'segments', 'edges on domain boundary (vertex1 vertex2 element)'))
    77         string = "%s\n%s" % (string, fielddisplay(self, 'segmentmarkers', 'number associated to each segment'))
    78         string = "%s\n%s" % (string, fielddisplay(self, 'vertexconnectivity', 'list of elements connected to vertex_i'))
    79         string = "%s\n%s" % (string, fielddisplay(self, 'elementconnectivity', 'list of elements adjacent to element_i'))
    80         string = "%s\n%s" % (string, fielddisplay(self, 'average_vertex_connectivity', 'average number of vertices connected to one vertex'))
    81 
    82         string += '\n      Extracted model():'
    83         string = "%s\n%s" % (string, fielddisplay(self, 'extractedvertices', 'vertices extracted from the model()'))
    84         string = "%s\n%s" % (string, fielddisplay(self, 'extractedelements', 'elements extracted from the model()'))
    85 
    86         return string
    87     # }}}
    88 
    89     def setdefaultparameters(self): # {{{
     54    #}}}
     55
     56    def __repr__(self): #{{{
     57        s = '   2D tria Mesh (surface):'
     58
     59        s += '\n      Elements and vertices:'
     60        s = "%s\n%s" % (s, fielddisplay(self, 'numberofelements', 'number of elements'))
     61        s = "%s\n%s" % (s, fielddisplay(self, 'numberofvertices', 'number of vertices'))
     62        s = "%s\n%s" % (s, fielddisplay(self, 'elements', 'vertex indices of the mesh elements'))
     63        s = "%s\n%s" % (s, fielddisplay(self, 'x', 'vertices x coordinate [m]'))
     64        s = "%s\n%s" % (s, fielddisplay(self, 'y', 'vertices y coordinate [m]'))
     65        s = "%s\n%s" % (s, fielddisplay(self, 'z', 'vertices z coordinate [m]'))
     66        s = "%s\n%s" % (s, fielddisplay(self, 'lat', 'vertices latitude [degrees]'))
     67        s = "%s\n%s" % (s, fielddisplay(self, 'long', 'vertices longitude [degrees]'))
     68        s = "%s\n%s" % (s, fielddisplay(self, 'r', 'vertices radius [m]'))
     69
     70        s = "%s\n%s" % (s, fielddisplay(self, 'edges', 'edges of the 2d mesh (vertex1 vertex2 element1 element2)'))
     71        s = "%s\n%s" % (s, fielddisplay(self, 'numberofedges', 'number of edges of the 2d mesh'))
     72
     73        s += '\n      Properties:'
     74        s = "%s\n%s" % (s, fielddisplay(self, 'vertexonboundary', 'vertices on the boundary of the domain flag list'))
     75        s = "%s\n%s" % (s, fielddisplay(self, 'segments', 'edges on domain boundary (vertex1 vertex2 element)'))
     76        s = "%s\n%s" % (s, fielddisplay(self, 'segmentmarkers', 'number associated to each segment'))
     77        s = "%s\n%s" % (s, fielddisplay(self, 'vertexconnectivity', 'list of elements connected to vertex_i'))
     78        s = "%s\n%s" % (s, fielddisplay(self, 'elementconnectivity', 'list of elements adjacent to element_i'))
     79        s = "%s\n%s" % (s, fielddisplay(self, 'average_vertex_connectivity', 'average number of vertices connected to one vertex'))
     80
     81        s += '\n      Extracted model():'
     82        s = "%s\n%s" % (s, fielddisplay(self, 'extractedvertices', 'vertices extracted from the model()'))
     83        s = "%s\n%s" % (s, fielddisplay(self, 'extractedelements', 'elements extracted from the model()'))
     84
     85        return s
     86    #}}}
     87
     88    def setdefaultparameters(self): #{{{
    9089        #The connectivity is the average number of nodes linked to a given node
    9190        #through an edge. This connectivity is used to initially allocate
    … …  
    9493        #test/NightlyRun/runme.py.
    9594        self.average_vertex_connectivity = 25
    96     # }}}
    97 
    98     def checkconsistency(self, md, solution, analyses): # {{{
     95    #}}}
     96
     97    def checkconsistency(self, md, solution, analyses): #{{{
    9998        md = checkfield(md, 'fieldname', 'mesh.x', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    10099        md = checkfield(md, 'fieldname', 'mesh.y', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    … …  
    116115
    117116        return md
    118     # }}}
    119 
    120     def marshall(self, prefix, md, fid): # {{{
     117    #}}}
     118
     119    def marshall(self, prefix, md, fid): #{{{
    121120        WriteData(fid, prefix, 'name', 'md.mesh.domain_type', 'data', 'Domain' + self.domaintype(), 'format', 'String')
    122121        WriteData(fid, prefix, 'name', 'md.mesh.domain_dimension', 'data', self.dimension(), 'format', 'Integer')
    … …  
    134133        WriteData(fid, prefix, 'object', self, 'fieldname', 'average_vertex_connectivity', 'format', 'Integer')
    135134        WriteData(fid, prefix, 'object', self, 'fieldname', 'vertexonboundary', 'format', 'DoubleMat', 'mattype', 1)
    136     # }}}
    137 
    138     def domaintype(self): # {{{
     135    #}}}
     136
     137    def domaintype(self): #{{{
    139138        return '3Dsurface'
    140     # }}}
    141 
    142     def dimension(self): # {{{
     139    #}}}
     140
     141    def dimension(self): #{{{
    143142        return 2
    144     # }}}
    145 
    146     def elementtype(self): # {{{
     143    #}}}
     144
     145    def elementtype(self): #{{{
    147146        return 'Tria'
    148     # }}}
    149 
    150     def processmesh(self, options): # {{{
     147    #}}}
     148
     149    def processmesh(self, options): #{{{
    151150        isplanet = 1
    152151        is2d = 0
    … …  
    156155        y = self.y
    157156        z = self.z
    158         return [x, y, z, elements, is2d, isplanet]
    159     # }}}
    160 
    161     def savemodeljs(self, fid, modelname): # {{{
     157
     158        return x, y, z, elements, is2d, isplanet
     159    #}}}
     160
     161    def savemodeljs(self, fid, modelname): #{{{
    162162        fid.write('  #s.mesh = new mesh3dsurface()\n' % modelname)
    163163        writejs1Darray(fid, [modelname, '.mesh.x'], self.x)
    … …  
    180180        writejs1Darray(fid, [modelname, '.mesh.extractedvertices'], self.extractedvertices)
    181181        writejs1Darray(fid, [modelname, '.mesh.extractedelements'], self.extractedelements)
    182     # }}}
    183 
    184     def export(self, *args): # {{{
     182    #}}}
     183
     184    def export(self, *args): #{{{
    185185        options = pairoptions(*args)
    186186
    … …  
    269269        #write style file:
    270270        applyqgisstyle(filename, 'mesh')
    271     # }}}
     271    #}}}

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

    r25172 r25455  
    7979class model(object):
    8080    #properties
    81     def __init__(self):  #{{{
    82         ''' classtype = model.properties
    83          for classe in dict.keys(classtype):
    84                print classe
    85                self.__dict__[classe] = classtype[str(classe)]
    86         '''
     81    def __init__(self): #{{{
    8782        self.mesh = mesh2d()
    8883        self.mask = mask()
    … …  
    232227
    233228    def extract(self, area):  # {{{
     229        """EXTRACT - extract a model according to an Argus contour or flag list
     230
     231        This routine extracts a submodel from a bigger model with respect to a given contour
     232        md must be followed by the corresponding exp file or flags list
     233        It can either be a domain file (argus type, .exp extension), or an array of element flags.
     234        If user wants every element outside the domain to be
     235        extract2d, add '~' to the name of the domain file (ex: '~HO.exp')
     236        an empty string '' will be considered as an empty domain
     237        a string 'all' will be considered as the entire domain
     238
     239        Usage:
     240            md2 = extract(md, area)
     241
     242        Examples:
     243            md2 = extract(md, 'Domain.exp')
     244
     245        See also: EXTRUDE, COLLAPSE
    234246        """
    235         extract - extract a model according to an Argus contour or flag list
    236 
    237            This routine extracts a submodel from a bigger model with respect to a given contour
    238            md must be followed by the corresponding exp file or flags list
    239            It can either be a domain file (argus type, .exp extension), or an array of element flags.
    240            If user wants every element outside the domain to be
    241            extract2d, add '~' to the name of the domain file (ex: '~HO.exp')
    242            an empty string '' will be considered as an empty domain
    243            a string 'all' will be considered as the entire domain
    244 
    245            Usage:
    246               md2 = extract(md, area)
    247 
    248            Examples:
    249               md2 = extract(md, 'Domain.exp')
    250 
    251            See also: EXTRUDE, COLLAPSE
    252         """
    253 
    254     #copy model
     247
     248        #copy model
    255249        md1 = copy.deepcopy(self)
    256250
    257     #get elements that are inside area
     251        #get elements that are inside area
    258252        flag_elem = FlagElements(md1, area)
    259253        if not np.any(flag_elem):
    260254            raise RuntimeError("extracted model is empty")
    261255
    262     #kick out all elements with 3 dirichlets
     256        #kick out all elements with 3 dirichlets
    263257        spc_elem = np.nonzero(np.logical_not(flag_elem))[0]
    264258        spc_node = np.unique(md1.mesh.elements[spc_elem, :]) - 1
    … …  
    268262        flag_elem[pos] = 0
    269263
    270     #extracted elements and nodes lists
     264        #extracted elements and nodes lists
    271265        pos_elem = np.nonzero(flag_elem)[0]
    272266        pos_node = np.unique(md1.mesh.elements[pos_elem, :]) - 1
    273267
    274     #keep track of some fields
     268        #keep track of some fields
    275269        numberofvertices1 = md1.mesh.numberofvertices
    276270        numberofelements1 = md1.mesh.numberofelements
    … …  
    280274        flag_node[pos_node] = 1
    281275
    282     #Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
     276        #Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
    283277        Pelem = np.zeros(numberofelements1, int)
    284278        Pelem[pos_elem] = np.arange(1, numberofelements2 + 1)
    … …  
    286280        Pnode[pos_node] = np.arange(1, numberofvertices2 + 1)
    287281
    288     #renumber the elements (some node won't exist anymore)
     282        #renumber the elements (some node won't exist anymore)
    289283        elements_1 = copy.deepcopy(md1.mesh.elements)
    290284        elements_2 = elements_1[pos_elem, :]
    … …  
    297291            elements_2[:, 5] = Pnode[elements_2[:, 5] - 1]
    298292
    299     #OK, now create the new model!
    300 
    301     #take every field from model
     293        #OK, now create the new model!
     294
     295        #take every field from model
    302296        md2 = copy.deepcopy(md1)
    303297
    304     #automatically modify fields
    305 
    306     #loop over model fields
     298        #automatically modify fields
     299
     300        #loop over model fields
    307301        model_fields = vars(md1)
    308302        for fieldi in model_fields:
    … …  
    310304            field = getattr(md1, fieldi)
    311305            fieldsize = np.shape(field)
    312             if hasattr(field, '__dict__') and fieldi not in ['results']:  #recursive call
     306            if hasattr(field, '__dict__') and fieldi not in ['results']: #recursive call
    313307                object_fields = vars(field)
    314308                for fieldj in object_fields:
    … …  
    416410        #recreate segments
    417411        if md1.mesh.__class__.__name__ == 'mesh2d':
    418             md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements, md2.mesh.numberofvertices)[0]
    419             md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements, md2.mesh.vertexconnectivity)[0]
     412            md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements, md2.mesh.numberofvertices)
     413            md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements, md2.mesh.vertexconnectivity)
    420414            md2.mesh.segments = contourenvelope(md2)
    421415            md2.mesh.vertexonboundary = np.zeros(numberofvertices2, bool)
    … …  
    423417        else:
    424418            #First do the connectivity for the contourenvelope in 2d
    425             md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements2d, md2.mesh.numberofvertices2d)[0]
    426             md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements2d, md2.mesh.vertexconnectivity)[0]
     419            md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements2d, md2.mesh.numberofvertices2d)
     420            md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements2d, md2.mesh.vertexconnectivity)
    427421            segments = contourenvelope(md2)
    428422            md2.mesh.vertexonboundary = np.zeros(int(numberofvertices2 / md2.mesh.numberoflayers), bool)
    … …  
    430424            md2.mesh.vertexonboundary = np.tile(md2.mesh.vertexonboundary, md2.mesh.numberoflayers)
    431425            #Then do it for 3d as usual
    432             md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements, md2.mesh.numberofvertices)[0]
    433             md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements, md2.mesh.vertexconnectivity)[0]
     426            md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements, md2.mesh.numberofvertices)
     427            md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements, md2.mesh.vertexconnectivity)
    434428
    435429        #Boundary conditions: Dirichlets on new boundary
    … …  
    504498                                    setattr(fieldr, solutionsubfield, subfield)
    505499
    506     #Keep track of pos_node and pos_elem
     500        #Keep track of pos_node and pos_elem
    507501        md2.mesh.extractedvertices = pos_node + 1
    508502        md2.mesh.extractedelements = pos_elem + 1
    … …  
    512506
    513507    def extrude(md, *args):  # {{{
    514         """
    515         EXTRUDE - vertically extrude a 2d mesh
    516 
    517            vertically extrude a 2d mesh and create corresponding 3d mesh.
    518            The vertical distribution can:
     508        """EXTRUDE - vertically extrude a 2d mesh
     509
     510        vertically extrude a 2d mesh and create corresponding 3d mesh.
     511        The vertical distribution can:
    519512        - follow a polynomial law
    520513        - follow two polynomial laws, one for the lower part and one for the upper part of the mesh
    … …  
    522515
    523516
    524            Usage:
    525               md = extrude(md, numlayers, extrusionexponent)
    526               md = extrude(md, numlayers, lowerexponent, upperexponent)
    527               md = extrude(md, listofcoefficients)
    528 
    529            Example:
    530                         md = extrude(md, 15, 1.3)
    531                         md = extrude(md, 15, 1.3, 1.2)
    532                         md = extrude(md, [0 0.2 0.5 0.7 0.9 0.95 1])
    533 
    534            See also: MODELEXTRACT, COLLAPSE
     517        Usage:
     518            md = extrude(md, numlayers, extrusionexponent)
     519            md = extrude(md, numlayers, lowerexponent, upperexponent)
     520            md = extrude(md, listofcoefficients)
     521
     522        Example:
     523            md = extrude(md, 15, 1.3)
     524            md = extrude(md, 15, 1.3, 1.2)
     525            md = extrude(md, [0 0.2 0.5 0.7 0.9 0.95 1])
     526
     527        See also: MODELEXTRACT, COLLAPSE
    535528        """
     529
    536530        #some checks on list of arguments
    537531        if len(args) > 3 or len(args) < 1:
    … …  
    698692    # }}}
    699693
    700     def collapse(md):  #{{{
    701         '''
    702         collapses a 3d mesh into a 2d mesh
     694    def collapse(md): #{{{
     695        """COLLAPSE - collapses a 3d mesh into a 2d mesh
    703696
    704697        This routine collapses a 3d model into a 2d model and collapses all
    … …  
    706699
    707700        Usage:
    708                 md = collapse(md)
    709         '''
     701            md = collapse(md)
     702        """
    710703
    711704        #Check that the model is really a 3d model
    … …  
    881874                md.mesh.scale_factor = project2d(md, md.mesh.scale_factor, 1)
    882875        md.mesh = mesh
    883         md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)[0]
    884         md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)[0]
     876        md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
     877        md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
    885878        md.mesh.segments = contourenvelope(md)
    886879
    887880        return md
    888 
    889881    #}}}

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

    r25125 r25455  
    11from collections import OrderedDict
    22
     3
    34class pairoptions(object):
    4     """
    5     PAIROPTIONS class definition
     5    """PAIROPTIONS class definition
    66
    7        Usage:
    8           pairoptions = pairoptions()
    9           pairoptions = pairoptions('module', true, 'solver', false)
     7    Usage:
     8        pairoptions = pairoptions()
     9        pairoptions = pairoptions('module', true, 'solver', false)
    1010    """
    1111
    … …  
    3939
    4040    def buildlist(self, *arg):  # {{{
    41         """BUILDLIST - build list of objects from input"""
     41        """BUILDLIST - build list of objects from input
     42        """
    4243
    4344        #check length of input

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

    r25136 r25455  
    2121                mergedcaps       = 0;
    2222                transitions      = {};
    23                 eltransitions      = {};
     23                eltransitions    = {};
    2424                planet           = '';
    2525                %}}}
    … …  
    170170                function intersections(self,varargin) % {{{
    171171
    172                 options=pairoptions(varargin{:});
    173                 force=getfieldvalue(options,'force',0);
     172                        options=pairoptions(varargin{:});
     173                        force=getfieldvalue(options,'force',0);
     174                       
     175                        %initialize, to avoid issues of having more transitions than meshes.
     176                        self.transitions={}; self.eltransitions={};
     177
     178                        %for elements:
     179                        xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3;
     180                        ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
     181                        ze=self.earth.mesh.z(self.earth.mesh.elements)*[1;1;1]/3;
     182
     183                        for i=1:length(self.icecaps),
     184                                mdi=self.icecaps{i};
     185                                mdi=TwoDToThreeD(mdi,self.planet);
    174186               
    175                 %initialize, to avoid issues of having more transitions than meshes.
    176                 self.transitions={}; self.eltransitions={};
    177 
    178                 %for elements:
    179                 xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3;
    180                 ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
    181                 ze=self.earth.mesh.z(self.earth.mesh.elements)*[1;1;1]/3;
    182 
    183                 for i=1:length(self.icecaps),
    184                         mdi=self.icecaps{i};
    185                         mdi=TwoDToThreeD(mdi,self.planet);
    186        
    187                         %for elements:
    188                         xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3;
    189                         yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3;
    190                         zei=mdi.mesh.z(mdi.mesh.elements)*[1;1;1]/3;
    191        
    192                         disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
     187                                %for elements:
     188                                xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3;
     189                                yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3;
     190                                zei=mdi.mesh.z(mdi.mesh.elements)*[1;1;1]/3;
    193191               
    194                         self.transitions{end+1}=meshintersect3d(self.earth.mesh.x,self.earth.mesh.y,self.earth.mesh.z,mdi.mesh.x,mdi.mesh.y,mdi.mesh.z,'force',force);
    195 
    196                         self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force);
    197                 end
     192                                disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
     193                       
     194                                self.transitions{end+1}=meshintersect3d(self.earth.mesh.x,self.earth.mesh.y,self.earth.mesh.z,mdi.mesh.x,mdi.mesh.y,mdi.mesh.z,'force',force);
     195
     196                                self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force);
     197                        end
    198198
    199199                end % }}}
    200200                function checkintersections(self) % {{{
    201                 flags=zeros(self.earth.mesh.numberofvertices,1);
    202                 for i=1:length(self.basins),
    203                         flags(self.transitions{i})=i;
    204                 end
    205                 plotmodel(self.earth,'data',flags,'coastline','on');
     201                        flags=zeros(self.earth.mesh.numberofvertices,1);
     202                        for i=1:length(self.basins),
     203                                flags(self.transitions{i})=i;
     204                        end
     205                        plotmodel(self.earth,'data',flags,'coastline','on');
    206206
    207207                end % }}}
    … …  
    227227                                                end
    228228                                                continent=unique(continent);
     229                                        else
     230                                                %nothing to do: assume we have a list of continents
    229231                                        end
    230232                                else
    231                                         %nothing to do, we have a list of continents
     233                                        %nothing to do: assume we have a list of continents
    232234                                end
    233235                        else
    … …  
    416418                                eval(['self.earth.' string '=field;']); %do not forget to plug the field variable into its final location
    417419                                %}}}
    418                         elseif n==self.icecaps{1} .mesh.numberofelements,
     420                        elseif n==self.icecaps{1}.mesh.numberofelements,
    419421                                eval(['self.earth.' string '=zeros(self.earth.mesh.numberofelements,1);']);
    420422                                for i=1:length(self.icecaps),

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

    r25158 r25455  
    77from meshintersect3d import meshintersect3d
    88from miscellaneous import miscellaneous
     9from model import model
     10from modelmerge3d import modelmerge3d
    911from pairoptions import pairoptions
    1012from private import private
    … …  
    1315
    1416class sealevelmodel(object):
    15     '''
    16     SEALEVELMODEL class definition
     17    """SEALEVELMODEL class definition
    1718
    1819    Usage:
    … …  
    2728            'earth', md_earth
    2829        )
    29     '''
     30    """
    3031
    3132    def __init__(self, *args): #{{{
    … …  
    207208        bas = options.getfieldvalue('basin', 'all')
    208209
    209         #expand continent list #{{{
     210        # Expand continent list #{{{
    210211        if type(continent) == np.ndarray:
    211212            if continent.shape[1] == 1:
    212213                if continent[0] == 'all':
    213                     #need to transform this into a list of continents
     214                    # Need to transform this into a list of continents
    214215                    continent = []
    215216                    for i in range(len(self.basins)):
    216217                        continent.append(self.basins[i].continent)
    217218                    continent = np.unique(continent)
     219                else:
     220                    pass # Nothing to do: assume we have a list of continents
    218221            else:
    219                 pass #nothing to do, we have a list of continents
     222                pass # Nothing to do: assume we have a list of continents
    220223        else:
    221224            if continent == 'all':
    222                 #need to transform this into a list of continents
     225                # Need to transform this into a list of continents
    223226                continent = []
    224227                for i in range(len(self.basins)):
    … …  
    226229                continent = np.unique(continent)
    227230            else:
    228                 continent = [continent]
     231                pass # Nothing to do: assume we have a list of continents
    229232        #}}}
    230         #expand basins list using the continent list above and the extra bas discriminator #{{{
     233
     234        # Expand basins list using the continent list above and the extra bas discriminator #{{{
    231235        if type(bas) == np.ndarray:
    232236            if bas.shape[1] == 1:
    233237                if bas[0] == 'all':
    234                     #need to transform this into a list of basins
     238                    # Need to transform this into a list of basins
    235239                    baslist = []
    236240                    for i in range(len(self.basins)):
    … …  
    246250                                baslist.append(i)
    247251            else:
    248                 #we have a list of basin names
     252                # We have a list of basin names
    249253                baslist = []
    250254                for i in range(len(bas)):
    … …  
    273277        #}}}
    274278    #}}}
     279
     280    def addicecap(self, md): #{{{
     281        if not type(md) ==  model:
     282            raise Exception("addicecap method only takes a 'model' class object as input")
     283
     284        self.icecaps.append(md)
     285    #}}}
     286
     287    def basinsplot3d(self, *args): #{{{
     288        for i in range(len(self.basins)):
     289            self.basins[i].plot3d(*args)
     290    #}}}
     291
     292    def caticecaps(self, *args): #{{{
     293        # Recover options
     294        options = pairoptions(*args)
     295        tolerance = options.getfieldvalue('tolerance', .65)
     296        loneedgesdetect = options.getfieldvalue('loneedgesdetect', 0)
     297
     298        # Make 3D model
     299        models = self.icecaps
     300        for i in range(len(models)):
     301            models[i] = TwoDToThreeD(models[i], self.planet)
     302
     303        # Plug all models together
     304        md = models[0]
     305        for i in range(1, len(models)):
     306            md = modelmerge3d(md, models[i], 'tolerance', tolerance)
     307            md.private.bamg.landmask = np.vstack((md.private.bamg.landmask, models[i].private.bamg.landmask))
     308
     309        # Look for lone edges if asked for it
     310        if loneedgesdetect:
     311            edges = loneedges(md)
     312            plotmodel(md, 'data', md.mask.land_levelset)
     313            for i in range(len(edges)):
     314                ind1 = edges(i, 1)
     315                ind1 = edges(i, 2)
     316                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*-')
     317
     318        # Plug into earth
     319        self.earth = md
     320
     321        # Create mesh radius
     322        self.earth.mesh.r = planetradius('earth')
     323    #}}}
     324
     325    def viscousiterations(self): #{{{
     326        for i in range(len(self.icecaps)):
     327            ic = self.icecaps[i]
     328            mvi = ic.resutls.TransientSolution[0].StressbalanceConvergenceNumSteps
     329            for j in range(1, len(ic.results.TransientSolution) - 1):
     330                mvi = np.max(mvi, ic.results.TransientSolution[j].StressbalanceConvergenceNumSteps)
     331            print("{}, {}: {}".format(i, self.icecaps[i].miscellaneous.name, mvi))
     332    #}}}
     333
     334    def maxtimestep(self): #{{{
     335        for i in range(len(self.icecaps)):
     336            ic = self.icecaps[i]
     337            mvi = len(ic.results.TransientSolution)
     338            timei = ic.results.TransientSolution[-1].time
     339            print("{}, {}: {}/{}".format(i, self.icecaps[i].miscellaneous.name, mvi, timei))
     340
     341        mvi = len(self.earth.results.TransientSolution)
     342        timei = self.earth.results.TransientSolution[-1].time
     343        print("Earth: {}/{}", mvi, timei)
     344    #}}}
     345   
     346    def transfer(self, string): #{{{
     347        # Recover field size in one icecap
     348        n = np.size(getattr(self.icecaps[i], string), 0)
     349
     350        if n == self.icecaps[0].mesh.numberofvertices:
     351            setattr(self.earth, string, np.zeros((self.earth.mesh.numberofvertices, )))
     352            for i in range(len(self.icecaps)):
     353                getattr(self.earth, string)[self.transitions[i]] = getattr(self.icecaps[i], string)
     354        elif n == (self.self.icecaps[0].mesh.numberofvertices + 1):
     355            # Dealing with transient dataset
     356            # Check that all timetags are similar between all icecaps #{{{
     357            for i in range(len(self.icecaps)):
     358                capfieldi = getattr(self.icecaps[i], string)
     359                for j in range(1, len(self.icecaps)):
     360                    capfieldj = getattr(self.icecaps[j], string)
     361                    if capfieldi[-1, :] != capfieldj[-1, :]:
     362                        raise Exception("Time stamps for {} field is different between icecaps {} and {}".format(string, i, j))
     363            capfield1 = getattr(self.icecaps[0], string)
     364            times = capfield1[-1, :]
     365            nsteps = len(times)
     366            #}}}
     367            # Initialize #{{{
     368            field = np.zeros((self.earth.mesh.numberofvertices + 1, nsteps))
     369            field[-1, :] = times # Transfer the times only, not the values
     370            #}}}
     371            # Transfer all the time fields #{{{
     372            for i in range(len(self.icecaps)):
     373                capfieldi = getattr(self.icecaps[i], string)
     374                for j in range(nsteps):
     375                    field[self.transitions[i], j] = capfieldi[0:-2, j] # Transfer only the values, not the time
     376            setattr(self.earth, string, field) # Do not forget to plug the field variable into its final location
     377            #}}}
     378        elif n == (self.icecaps[0].mesh.numberofelements):
     379            setattr(self.earth, string, np.zeros((self.earth.mesh.numberofvertices, )))
     380            for i in range(len(self.icecaps)):
     381                getattr(self.earth, string)[self.eltransitions[i]] = getattr(self.icecaps[i], string)
     382        else:
     383            raise Exception('not supported yet')
     384    #}}}
     385
     386    def homogenize(self, noearth=0): #{{{
     387        mintimestep = np.inf
     388
     389        for i in range(len(self.icecaps)):
     390            ic = self.icecaps[i]
     391            mintimestep = np.min(mintimestep, len (ic.results.TransientSolution))
     392
     393        if not noearth:
     394            mintimestep = np.min(mintimestep, len(self.earth.results.TransientSolution))
     395
     396        for i in range(len(self.icecaps)):
     397            ic = self.icecaps[i]
     398            ic.resuts.TransientSolution = ic.results.TransientSolution[:mintimestep]
     399            self.icecaps[i] = ic
     400
     401        ic = self.earth
     402
     403        if not noearth:
     404            ic.results.TransientSolution = ic.resutls.TransientSolution[:mintimestep]
     405
     406        self.earth = ic
     407
     408        return self
     409    #}}}
     410
     411    def initializemodels(self): #{{{
     412        for i in range(len(self.basins)):
     413            md = model()
     414            md.miscellaneous.name = self.basins[i].name
     415            self.addicecap(md)
     416    #}}}

  • issm/trunk-jpl/src/m/coordsystems/epsg2proj.m

    r25189 r25455  
    11function string = epsg2proj(epsg)
    2 %FUNCTION EPSG2PROJ - uses gdalsrsinfo to provide PROJ.4 compatible string from
     2%EPSG2PROJ - uses gdalsrsinfo to provide PROJ.4 compatible string from
    33%EPSG code
    44%

  • issm/trunk-jpl/src/m/coordsystems/epsg2proj.py

    r25346 r25455  
    1 import shlex
    21import subprocess
    32
    43
    54def epsg2proj(epsg): #{{{
    6     """
    7     FUNCTION EPSG2PROJ - uses gdalsrsinfo to provide PROJ.4 compatible string
     5    """EPSG2PROJ - uses gdalsrsinfo to provide PROJ.4 compatible string
    86    from EPSG code
    97
    … …  
    2220
    2321    #First, get GDAL version
    24     subproc_args = shlex.split("gdalsrsinfo --version | awk '{print $2}' | cut -d '.' -f1")
     22    subproc_args = "gdalsrsinfo --version | awk '{print $2}' | cut -d '.' -f1"
    2523    subproc = subprocess.Popen(subproc_args, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    2624    outs, errs = subproc.communicate()
    … …  
    3028    version_major=int(outs)
    3129
    32     subproc_args = shlex.split("gdalsrsinfo epsg:{} | command grep PROJ.4 | tr -d '\n' | sed 's/PROJ.4 : //'".format(epsg))
     30    subproc_args = "gdalsrsinfo epsg:{} | command grep PROJ.4 | tr -d '\n' | sed 's/PROJ.4 : //'".format(epsg)
    3331    subproc = subprocess.Popen(subproc_args, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    3432    outs, errs = subproc.communicate()

  • issm/trunk-jpl/src/m/coordsystems/gdaltransform.m

    r25065 r25455  
    2121%
    2222%       To get proj.4 string from EPSG, use gdalsrsinfo. Example:
    23 %
    2423%               gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //"
    2524
    … …  
    3231        fclose(fid);
    3332
    34         [s,r]=system(['gdaltransform -s_srs "',proj_in,'" -t_srs "',proj_out,'"  < ' filename_in ' > ' filename_out]);
     33        [s,r]=system(['gdaltransform -s_srs "',proj_in,'" -t_srs "',proj_out,'" < ' filename_in ' > ' filename_out]);
    3534        if s~=0 | ~isempty(deblank(r)),
    3635                error(r);
    3736        end
     37
    3838        A=load(filename_out);
    39         xout=A(:,1); xout=reshape(xout,size(x));
    40         yout=A(:,2); yout=reshape(yout,size(y));
    4139
    4240        %clean up
    4341        delete(filename_in);
    4442        delete(filename_out);
     43
     44        xout=A(:,1);
     45        xout=reshape(xout,size(x));
     46        yout=A(:,2);
     47        yout=reshape(yout,size(y));

  • issm/trunk-jpl/src/m/coordsystems/gdaltransform.py

    r25163 r25455  
    44import tempfile
    55
     6import numpy as np
     7
    68from loadvars import *
    79
    810
    911def gdaltransform(x, y, proj_in, proj_out): #{{{
    10     '''
    11     GDALTRANSFORM - switch from one projection system to another
     12    """GDALTRANSFORM - switch from one projection system to another
    1213
    1314    Usage:
    … …  
    2930        +proj = stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
    3031
    31     To get proj.4 string form EPSG, use gdalsrsinfo. Example:
    32 
     32    To get proj.4 string from EPSG, use gdalsrsinfo. Example:
    3333        gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //"
    34 
    35     TODO:
    36     - Follow subprocess pattern implemented in src/m/coordsystems/epsg2proj.py
    37     '''
     34    """
    3835
    3936    # Give ourselves unique file names
    40     file_in = tempfile.NamedTemporaryFile('r+b')
     37    file_in = tempfile.NamedTemporaryFile('w', delete=False)
    4138    filename_in = file_in.name
    42     file_out = tempfile.NamedTemporaryFile('w+b')
     39    file_out = tempfile.NamedTemporaryFile('w', delete=False)
    4340    filename_out = file_out.name
    4441
    45     file_in.write(b'%8g %8g\n' % (x.flatten(1), y.flatten(1)))
     42    points = np.vstack((x, y)).T
     43    np.savetxt(file_in, points, fmt='%8g %8g')
     44    file_in.close() # NOTE: Opening file in 'r+' or 'w+' mode does not allow subsequent reading by subprocess. We therefore need to close it and reopen it.
     45    file_in = open(filename_in) # Open for reading by subprocess
     46
     47    subproc_args = shlex.split("gdaltransform -s_srs '{}' -t_srs '{}'".format(proj_in, proj_out))
     48    subproc = subprocess.Popen(subproc_args, bufsize=-1, stdin=file_in, stdout=file_out, stderr=subprocess.PIPE, close_fds=True)
     49    outs, errs = subproc.communicate()
     50    if errs != '':
     51        raise RuntimeError("gdaltransform: call to gdaltransform failed: {}".format(errs))
     52
     53    A = np.loadtxt(filename_out)
     54
     55    # Clean up
    4656    file_in.close()
     57    file_out.close()
     58    remove(filename_in)
     59    remove(filename_out)
    4760
    48     subproc_args = shlex.split('gdaltransform -s_srs {} -t_srs {} < {} > {}'.format(proj_in, proj_out, filename_in, filename_out))
    49     subprocess.check_call(subproc_args) # Will raise CalledProcessError if return code is not 0
    50 
    51     A = loadvars(filename_out)
    52     xout = A[0]
    53     xout = xout.reshape(x.shape)
    54     yout = A[1]
    55     yout = yout.reshape(y.shape)
    56 
    57     os.remove(filename_in)
    58     os.remove(filename_out)
     61    if np.ndim(A) > 1:
     62        xout = np.array(A[:,0])
     63        yout = np.array(A[:,1])
     64        # if type(x) == "np.ndarray":
     65        #     xout = xout.reshape(x.shape) # TODO: Do we need to do this?
     66        #     yout = yout.reshape(y.shape) # TODO: Do we need to do this?
     67    else:
     68        xout = [A[0]]
     69        yout = [A[1]]
    5970
    6071    return [xout, yout]

  • issm/trunk-jpl/src/m/coordsystems/laea.m

    r24777 r25455  
    99%      return string='+proj=laea +lat_0=45 +lon_0=-90 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs'
    1010
    11         string=sprintf('+proj=laea +lat_0=%i +lon_0=%i +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs ',lat,long);
     11        string=sprintf('+proj=laea +lat_0=%i +lon_0=%i +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs',lat,long);

  • issm/trunk-jpl/src/m/exp/expread.py

    r24254 r25455  
     1from collections import OrderedDict
    12import os.path
     3
    24import numpy as np
    3 from collections import OrderedDict
     5
    46import MatlabFuncs as m
    57
    68
    79def expread(filename):
     10    """EXPREAD - read a exp file and build a list of OrderedDicts
     11
     12    This routine reads a file .exp and builds a list of OrderedDicts containing
     13    the fields x and y corresponding to the coordinates, one for the filename
     14    of the exp file, for the density, for the nodes, and a field closed to
     15    indicate if the domain is closed. The first argument is the .exp file to be
     16    read and the second one (optional) indicate if the last point shall be read
     17    (1 to read it, 0 not to).
     18
     19    Usage:
     20        contours = expread(filename)
     21
     22    Example:
     23        contours = expread('domainoutline.exp')
     24        contours = expread('domainoutline.exp')
     25
     26    See Also:
     27    - EXPDOC
     28    - EXPWRITEASVERTICES
     29
     30    TODO:
     31    - Convert returned data structure from list of OrderedDict objects to list
     32    of OrderedStruct objects (see also src/m/shp/shpread.py). Also, modify
     33    handling of returned data structure in,
     34        - src/m/exp/expcoarsen.py
     35        - src/m/exp/expdisp.py
     36        - src/m/interp/SectionValues.py
     37        - src/m/mesh/bamg.py
     38    May also need to modify addressing in corresponding FetchData function, or
     39    create new one, in src/wrappers/ContoursToNodes/ContoursToNodes.cpp.
    840    """
    941
    10     EXPREAD - read a file exp and build a Structure
    11 
    12        This routine reads a file .exp and builds a list of dicts containing the
    13        fields x and y corresponding to the coordinates, one for the filename of
    14        the exp file, for the density, for the nodes, and a field closed to
    15        indicate if the domain is closed.
    16        The first argument is the .exp file to be read and the second one (optional)
    17        indicate if the last point shall be read (1 to read it, 0 not to).
    18 
    19        Usage:
    20           contours = expread(filename)
    21 
    22        Example:
    23           contours = expread('domainoutline.exp')
    24           contours = expread('domainoutline.exp')
    25 
    26        See also EXPDOC, EXPWRITEASVERTICES
    27 
    28     """
    2942    #some checks
    3043    if not os.path.exists(filename):

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

    r21067 r25455  
    22%FIND_POINT - find closest point
    33%
    4 %   find which point of the list (tabx,taby) is
    5 %   the closest to (poinx,pointy)
     4%   find which point of the list (tabx,taby) is the closest to (pointx,pointy)
    65%
    76%   Usage:

  • issm/trunk-jpl/src/m/geometry/polyarea.py

    r25125 r25455  
    55
    66def polyarea(x, y): #{{{
    7     '''
     7    """
    88    POLYAREA - returns the area of the 2-D polygon defined by the vertices in
    99    lists x and y
    … …  
    2424    - Test that output falls within some tolerance of MATLAB's polyarea
    2525    function.
    26     '''
     26    """
     27
    2728    return 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1)))
    2829#}}}

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

    r24256 r25455  
    11import os
     2
     3import numpy as np
     4
    25from expread import expread
    3 import numpy as np
    4 from project2d import project2d
    56#from InterpFromMesh2d import InterpFromMesh2d
    67from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
    78from InterpFromMeshToMesh3d import InterpFromMeshToMesh3d
     9from project2d import project2d
    810
    911
    1012def SectionValues(md, data, infile, resolution):
    11     '''
    12     compute the value of a field on a section
     13    """SECTIONVALUES - compute the value of a field on a section
    1314
    1415    This routine gets the value of a given field of the model on points
    … …  
    1718
    1819    Usage:
    19     [elements, x, y, z, s, data] = SectionValues(md, data, filename, resolution)
    20     [elements, x, y, z, s, data] = SectionValues(md, data, profile_structure, resolution)
    21     '''
     20        [elements, x, y, z, s, data] = SectionValues(md, data, filename, resolution)
     21        [elements, x, y, z, s, data] = SectionValues(md, data, profile_structure, resolution)
     22    """
    2223
    2324    if os.path.isfile(infile):
    … …  
    8485
    8586        #Interpolation of data on specified points
    86         #data_interp = InterpFromMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y, data, X, Y)[0]
    87         data_interp = InterpFromMeshToMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y, data, X, Y)[0]
     87        #data_interp = InterpFromMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y, data, X, Y)
     88        data_interp = InterpFromMeshToMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y, data, X, Y)
    8889    #data_interp = griddata(md.mesh.x, md.mesh.y, data, X, Y)
    8990
    … …  
    9596        #Get base and surface for each 2d point, offset to make sure that it is inside the glacier system
    9697        offset = 1.e-3
    97         base = InterpFromMeshToMesh2d(md.mesh.elements2d, md.mesh.x2d, md.mesh.y2d, project2d(md, md.geometry.base, 1), X, Y)[0] + offset
     98        base = InterpFromMeshToMesh2d(md.mesh.elements2d, md.mesh.x2d, md.mesh.y2d, project2d(md, md.geometry.base, 1), X, Y) + offset
    9899        base = base.reshape(-1, )
    99         surface = InterpFromMeshToMesh2d(md.mesh.elements2d, md.mesh.x2d, md.mesh.y2d, project2d(md, md.geometry.surface, 1), X, Y)[0] - offset
     100        surface = InterpFromMeshToMesh2d(md.mesh.elements2d, md.mesh.x2d, md.mesh.y2d, project2d(md, md.geometry.surface, 1), X, Y) - offset
    100101        surface = surface.reshape(-1, )
    101102
    … …  
    126127
    127128    #Interpolation of data on specified points
    128         data_interp = InterpFromMeshToMesh3d(md.mesh.elements, md.mesh.x, md.mesh.y, md.mesh.z, data, X3, Y3, Z3, np.nan)[0]
     129        data_interp = InterpFromMeshToMesh3d(md.mesh.elements, md.mesh.x, md.mesh.y, md.mesh.z, data, X3, Y3, Z3, np.nan)
    129130
    130131    #build outputs

  • issm/trunk-jpl/src/m/io/loadmodel.m

    r13007 r25455  
    11function varargout=loadmodel(path)
    2 %LOADMODEL - load a model using built-in load module
     2%LOADMODEL - load a model using built-in load module 
    33%
    44%   check that model prototype has not changed. if so, adapt to new model prototype.

  • issm/trunk-jpl/src/m/io/loadmodel.py

    r25346 r25455  
    11from loadvars import loadvars
    2 #hack to keep python 2 compatipility
     2# Hack to keep python 2 compatibility
    33try:
    4     #py3 import
    5     from dbm import whichdb
     4    from dbm import whichdb # Python 3
    65except ImportError:
    7     #py2 import
    8     from whichdb import whichdb
     6    from whichdb import whichdb # Python 2
     7
    98from netCDF4 import Dataset
    109
    1110
    1211def loadmodel(path, onlylast=False):
    13     """
    14     LOADMODEL - load a model using built - in load module
     12    """LOADMODEL - load a model
    1513
    16        check that model prototype has not changed. if so, adapt to new model prototype.
     14    Check that model prototype has not changed: if if has, adapt to new model prototype.
    1715
    18        Usage:
    19           md = loadmodel(path)
     16    Usage:
     17        md = loadmodel(path)
    2018    """
    2119

  • issm/trunk-jpl/src/m/io/loadvars.py

    r25346 r25455  
    11from collections import OrderedDict
     2# Hack to keep python 2 compatibility
     3try:
     4    from dbm import whichdb # Python 3
     5except ImportError:
     6    from whichdb import whichdb # Python 2
     7from re import findall
     8import shelve
     9
    210from netCDF4 import Dataset
    311import numpy as np
    4 from re import findall
    5 import shelve
    612
    713from model import *
    8 #hack to keep python 2 compatibility
    9 try:
    10     #py3 import
    11     from dbm import whichdb
    12 except ImportError:
    13     #py2 import
    14     from whichdb import whichdb
    1514
    1615
    1716def loadvars(*args, **kwargs):
    18     """
    19     LOADVARS - function to load variables from a file.
     17    """LOADVARS - function to load variables from a file
    2018
    2119    This function loads one or more variables from a file. The names of the

  • issm/trunk-jpl/src/m/mech/newforcing.m

    r20186 r25455  
    11function forcing=newforcing(t0,t1,deltaT,f0,f1,nodes)
    2 %FUNCTION NEWFORCING Build forcing that extends temporally from t0 to t1, and in magnitude from f0 to f1. Equal time
    3 %                    and magnitude spacing.
     2%NEWFORCING - Build forcing that extends temporally from t0 to t1, and in
     3%magnitude from f0 to f1. Equal time and magnitude spacing.
    44%
    5 %       Usage: forcing=newforcing(t0,t1,deltaT,f0,f1,nodes); 
    6 %       Where:
    7 %          t0:t1: time interval.
    8 %          deltaT: time step
    9 %          f0:f1: magnitude interval.
    10 %          nodes: number of vertices where we have a temporal forcing
     5%   Usage: forcing=newforcing(t0,t1,deltaT,f0,f1,nodes); 
     6%   
     7%   Where:
     8%      t0:t1: time interval.
     9%      deltaT: time step
     10%      f0:f1: magnitude interval.
     11%      nodes: number of vertices where we have a temporal forcing
    1112%
    12 %       Example:
    13 %           md.smb.mass_balance=newforcing(md.timestepping.start_time,md.timestepping.final_time,...
    14 %                                          md.timestepping.time_step,-1,+2,md.mesh.numberofvertices);
     13%   Example:
     14%      md.smb.mass_balance=newforcing(md.timestepping.start_time,md.timestepping.final_time,md.timestepping.time_step,-1,+2,md.mesh.numberofvertices);
    1515%
    1616

  • issm/trunk-jpl/src/m/mech/newforcing.py

    r24256 r25455  
    44def newforcing(t0, t1, deltaT, f0, f1, nodes):
    55    '''
    6 FUNCTION NEWFORCING Build forcing that extends temporally from t0 to t1, and in magnitude from f0 to f1. Equal time
     6    NEWFORCING - Build forcing that extends temporally from t0 to t1, and in magnitude from f0 to f1. Equal time
    77                    and magnitude spacing.
    88
    9        Usage: forcing = newforcing(t0, t1, deltaT, f0, f1, nodes);
    10        Where:
    11           t0:t1: time interval.
    12           deltaT: time step
    13           f0:f1: magnitude interval.
    14           nodes: number of vertices where we have a temporal forcing
     9    Usage: forcing = newforcing(t0, t1, deltaT, f0, f1, nodes);
     10   
     11    Where:
     12        t0:t1: time interval.
     13        deltaT: time step
     14        f0:f1: magnitude interval.
     15        nodes: number of vertices where we have a temporal forcing
    1516
    16        Example:
    17            md.smb.mass_balance = newforcing(md.timestepping.start_time, md.timestepping.final_time,
    18                                           md.timestepping.time_step, -1, +2, md.mesh.numberofvertices)
     17    Example:
     18        md.smb.mass_balance = newforcing(
     19            md.timestepping.start_time,
     20            md.timestepping.final_time,
     21            md.timestepping.time_step,
     22            -1,
     23            +2,
     24            md.mesh.numberofvertices
     25            )
    1926    '''
     27   
    2028    #Number of time steps:
    2129    nsteps = (t1 - t0) / deltaT + 1

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

    r24999 r25455  
    2121        vc=md.mesh.vertexconnectivity;
    2222        vb=md.mesh.vertexonboundary;
    23         md.mesh=mesh3dsurface(); 
     23        md.mesh=mesh3dsurface();
    2424        md.mesh.lat=lat;
    2525        md.mesh.long=long;

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

    r25035 r25455  
    11import numpy as np
    22
     3from epsg2proj import epsg2proj
     4from gdaltransform import gdaltransform
    35from mesh3dsurface import mesh3dsurface
    46from planetradius import planetradius
    57
    6 def TwoDToThreeD(md, planet): #{{{
     8
     9def TwoDToThreeD(md, planet):
    710    # Reproject model into lat, long if necessary
    811    if md.mesh.proj != epsg2proj(4326):
    … …  
    1316
    1417    # We assume x and y hold the long, lat values
    15     long = md.mesh.x
    16     lat = md.mesh.y
     18    longe = md.mesh.x
     19    late = md.mesh.y
    1720
    1821    # Assume spherical body
    19     x = R * np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(long))
    20     y = R * np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(long))
    21     z = R * np.sin(np.deg2rad(lat))
     22    x = R * np.cos(np.deg2rad(late)) * np.cos(np.deg2rad(longe))
     23    y = R * np.cos(np.deg2rad(late)) * np.sin(np.deg2rad(longe))
     24    z = R * np.sin(np.deg2rad(late))
    2225
    2326    elements = md.mesh.elements
    2427    vc = md.mesh.vertexconnectivity
    2528    vb = md.mesh.vertexonboundary
    26     md.mesh.mesh3dsurface()
    27     md.mesh.lat = lat
    28     md.mesh.long = long
     29    md.mesh = mesh3dsurface()
     30    md.mesh.lat = late
     31    md.mesh.long = longe
    2932    md.mesh.x = x
    3033    md.mesh.y = y
    … …  
    3235    md.mesh.elements = elements
    3336    md.mesh.numberofelements = len(elements)
    34     md.mesh.numberofvertices = len(lat)
    35     md.mesh.r = R * np.ones(md.mesh.numberofvertices)
     37    md.mesh.numberofvertices = len(late)
     38    md.mesh.r = R * np.ones((md.mesh.numberofvertices, ))
    3639    md.mesh.vertexconnectivity = vc
    3740    md.mesh.vertexonboundary = vb
    38 #}}}
     41
     42    return md

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

    r24865 r25455  
    8080                        domain=shpread(domainfile);
    8181                else
    82                         error(['bamg error message: file ' domainfile ' format not supported (.shp or .exp)']);
     82                        error(['bamg error message: file ' domainfile ' format not supported (.exp or .shp)']);
    8383                end
    8484        elseif isstruct(domainfile),
    … …  
    140140                end
    141141
    142                 %Checks that all holes are INSIDE the principle domain outline
     142                %Check that all holes are INSIDE the principle domain outline
    143143                if i>1,
    144144                        flags=ContourToNodes(domain(i).x,domain(i).y,domain(1),0);
    … …  
    149149
    150150                %Check orientation
    151                 nods=domain(i).nods-1; %the domain are closed 1=end;
     151                nods=domain(i).nods-1; %the domain is closed (domain[1] = domain[end])
    152152                test = sum([(domain(i).x(2:nods+1) - domain(i).x(1:nods)).*(domain(i).y(2:nods+1) + domain(i).y(1:nods))]);
    153153                if (i==1 && test>0) || (i>1 && test<0),
    … …  
    163163                bamg_geometry.Edges   =[bamg_geometry.Edges;    [transpose(count+1:count+nods) transpose([count+2:count+nods count+1])  1*ones(nods,1)]];
    164164
    165                 %Flag how many edges we have now, that way we know which edges belong to the domain, will
    166                 %be used later for required edges when NoBoundaryRefinement is one:
     165                % Flag how many edges we have now, that way we know which edges belong
     166                % to the domain. Will be used later for required edges if
     167                % NoBoundaryRefinement equals 1.
    167168                new_edge_length=length(bamg_geometry.Edges);
    168169                edges_required=(edge_length+1):new_edge_length;
    169170
    170171                if i>1,
    171                         bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 -subdomain_ref]; subdomain_ref = subdomain_ref+1;
     172                        bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 -subdomain_ref];
     173                        subdomain_ref = subdomain_ref+1;
    172174                else
    173175                        bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 0];
    … …  
    179181        for i=1:length(holes),
    180182
    181                 %Check that the subdomains is closed
     183                %Check that the hole is closed
    182184                if (holes(i).x(1)~=holes(i).x(end) | holes(i).y(1)~=holes(i).y(end)),
    183185                        error('bamg error message: all contours provided in ''holes'' should be closed');
    184186                end
    185187
    186                 %Checks that all holes are INSIDE the principle domain outline
     188                %Check that all holes are INSIDE the principal domain (principal domain should be index 0)
    187189                flags=ContourToNodes(holes(i).x,holes(i).y,domain(1),0);
    188190                if any(~flags), error('bamg error message: All holes should be strictly inside the principal domain'); end
    189191
    190192                %Check that hole is correctly oriented
    191                 nods=holes(i).nods-1; %the holes are closed 1=end;
     193                nods=holes(i).nods-1; %the hole is closed (hole[1] = hole[end])
    192194                if(sum([(holes(i).x(2:nods+1) - holes(i).x(1:nods)).*(holes(i).y(2:nods+1) + holes(i).y(1:nods))]))<0
    193195                        disp('At least one hole was not correctly oriented and has been re-oriented');
    … …  
    205207        for i=1:length(subdomains),
    206208
    207                 %Check that the subdomains is closed
     209                %Check that the subdomain is closed
    208210                if (subdomains(i).x(1)~=subdomains(i).x(end) | subdomains(i).y(1)~=subdomains(i).y(end)),
    209211                        error('bamg error message: all contours provided in ''subdomains'' should be closed');
    210212                end
    211213
    212                 %Checks that all holes are INSIDE the principle domain outline
     214                %Checks that all subdomains are INSIDE the principal domain (principal domain should be index 0)
    213215                flags=ContourToNodes(subdomains(i).x,subdomains(i).y,domain(1),0);
    214216                if any(~flags),
    215                         error('bamg error message: All holes should be strictly inside the principal domain');
    216                 end
    217 
    218                 %Check that hole is correctly oriented
    219                 nods=subdomains(i).nods-1; %the subdomains are closed 1=end;
     217                        error('bamg error message: All subdomains should be strictly inside the principal domain');
     218                end
     219
     220                %Check that subdomain is correctly oriented
     221                nods=subdomains(i).nods-1; % the subdomains are closed (subdomains[1] = subdomains[end])
    220222                if(sum([(subdomains(i).x(2:nods+1) - subdomains(i).x(1:nods)).*(subdomains(i).y(2:nods+1) + subdomains(i).y(1:nods))]))>0
    221223                        disp('At least one subdomain was not correctly oriented and has been re-oriented');
    222                         subdomains(i).x = flipud(subdomains(i).x); subdomains(i).y = flipud(subdomains(i).y);
     224                        subdomains(i).x = flipud(subdomains(i).x);
     225                        subdomains(i).y = flipud(subdomains(i).y);
    223226                end
    224227
    … …  
    227230                bamg_geometry.Edges   =[bamg_geometry.Edges;    [transpose(count+1:count+nods) transpose([count+2:count+nods count+1])  1.*ones(nods,1)]];
    228231               
    229                 bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 subdomain_ref]; subdomain_ref = subdomain_ref+1;
     232                bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 subdomain_ref];
     233                subdomain_ref = subdomain_ref+1;
    230234               
    231235                %update counter
    232236                count=count+nods;
    233237        end
     238
    234239        if getfieldvalue(options,'vertical',0),
    235240                if numel(getfieldvalue(options,'Markers',[]))~=size(bamg_geometry.Edges,1),
    … …  
    243248        %take care of rifts
    244249        if exist(options,'rifts'),
    245 
    246                 %Check that file exists
     250                %read rift file according to its extension:
    247251                riftfile=getfieldvalue(options,'rifts');
    248                 [pathr,namer,extr]=fileparts(riftfile);
    249                 if ~exist(riftfile,'file')
    250                         error(['bamg error message: file ' riftfile ' not found ']);
    251                 elseif strcmp(extr,'.exp'),
    252                         rift=expread(riftfile);
    253                 elseif strcmp(extr,'.shp'),
    254                         rift=shpread(riftfile);
    255                 end
    256                 %read rift file according to its extension:
    257252                [path,name,ext]=fileparts(riftfile);
    258253                if strcmp(ext,'.exp'),
    … …  
    261256                        rift=shpread(riftfile);
    262257                else
    263                         error(['bamg error message: file ' riftfile ' format not supported (.shp or .exp)']);
     258                        error(['bamg error message: file ' riftfile ' format not supported (.exp or .shp)']);
    264259                end
    265260
    … …  
    272267
    273268                        elseif any(~flags),
    274                                 %We LOTS of work to do
     269                                %We have LOTS of work to do
    275270                                disp('Rift tip outside of or on the domain has been detected and is being processed...');
    276271
    277272                                %check that only one point is outside (for now)
    278273                                if sum(~flags)~=1,
    279                                         error('bamg error message: only one point outside of the domain is supported yet');
     274                                        error('bamg error message: only one point outside of the domain is supported at this time');
    280275                                end
    281276
    … …  
    290285                                end
    291286
    292                                 %Get cordinate of intersection point
    293                                 x1=rift(i).x(1); y1=rift(i).y(1);
    294                                 x2=rift(i).x(2); y2=rift(i).y(2);
     287                                %Get coordinate of intersection point
     288                                x1=rift(i).x(1);
     289                                y1=rift(i).y(1);
     290                                x2=rift(i).x(2);
     291                                y2=rift(i).y(2);
    295292                                for j=1:length(domain(1).x)-1;
    296293                                        if SegIntersect([x1 y1; x2 y2],[domain(1).x(j) domain(1).y(j); domain(1).x(j+1) domain(1).y(j+1)]),
    … …  
    489486        md.mesh.numberofvertices=length(md.mesh.x);
    490487        md.mesh.numberofedges=size(md.mesh.edges,1);
    491         md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
    492 
    493         %Now, build the connectivity tables for this mesh.
    494         md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
     488        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
     489        md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
    495490
    496491elseif getfieldvalue(options,'3dsurface',0),
    … …  
    509504        md.mesh.numberofvertices=length(md.mesh.x);
    510505        md.mesh.numberofedges=size(md.mesh.edges,1);
    511         md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
     506        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
     507        md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
    512508
    513509else
    … …  
    524520        md.mesh.numberofvertices=length(md.mesh.x);
    525521        md.mesh.numberofedges=size(md.mesh.edges,1);
    526         md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
     522        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
     523        md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
    527524end
    528525

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

    r24262 r25455  
     1from collections import namedtuple, OrderedDict
    12import os.path
     3
    24import numpy as np
     5
     6from bamggeom import bamggeom
     7from BamgMesher import BamgMesher
     8from bamgmesh import bamgmesh
     9from ContourToNodes import ContourToNodes
     10from expread import expread
     11from helpers import fileparts, OrderedStruct
    312from mesh2d import *
    413from mesh2dvertical import *
    514from mesh3dsurface import *
    6 from collections import OrderedDict
    715from pairoptions import pairoptions
    8 from bamggeom import bamggeom
    9 from bamgmesh import bamgmesh
    10 from expread import expread
    1116from SegIntersect import SegIntersect
    12 import MatlabFuncs as m
    13 from BamgMesher import BamgMesher
    14 from ContourToNodes import ContourToNodes
     17from shpread import shpread
    1518
    1619
    1720def bamg(md, *args):
    18     """
    19     BAMG - mesh generation
     21    """BAMG - mesh generation
    2022
    2123    Available options (for more details see ISSM website http://issm.jpl.nasa.gov/):
    … …  
    4345                            1 -> use Green formula
    4446    - KeepVertices :      try to keep initial vertices when adaptation is done on an existing mesh (default 1)
     47    - NoBoundaryRefinement: do not refine boundary, only follow contour provided (default 0). Allow subdomain boundary refinement though
     48    - NoBoundaryRefinementAllBoundaries: do not refine boundary, only follow contour provided (default 0)
    4549    - maxnbv :            maximum number of vertices used to allocate memory (default is 1.0e6)
    4650    - maxsubdiv :         maximum subdivision of exisiting elements (default is 10)
    … …  
    6367                             will be merged
    6468
    65        Examples:
    66           md = bamg(md, 'domain', 'DomainOutline.exp', 'hmax', 3000)
    67           md = bamg(md, 'field', [md.inversion.vel_obs md.geometry.thickness], 'hmax', 20000, 'hmin', 1000)
    68           md = bamg(md, 'metric', A, 'hmin', 1000, 'hmax', 20000, 'gradation', 3, 'anisomax', 1)
     69    Examples:
     70        md = bamg(md, 'domain', 'DomainOutline.exp', 'hmax', 3000)
     71        md = bamg(md, 'field', [md.inversion.vel_obs md.geometry.thickness], 'hmax', 20000, 'hmin', 1000)
     72        md = bamg(md, 'metric', A, 'hmin', 1000, 'hmax', 20000, 'gradation', 3, 'anisomax', 1)
     73
     74    TODO:
     75    - Verify that values of third column of bamg_geometry.Vertices and bamg_geometry.Edges are valid (compare to src/m/mesh/bamg.m)
    6976    """
    7077
    … …  
    7986
    8087    subdomain_ref = 1
     88    hole_ref = 1
    8189
    8290    # Bamg Geometry parameters {{{
    … …  
    8694        if type(domainfile) == str:
    8795            if not os.path.exists(domainfile):
    88                 raise IOError("bamg error message: file '%s' not found" % domainfile)
     96                raise IOError("bamg error message: file {} not found".format( domainfile))
     97
     98            # Read domain according to its extension
     99            path, name, ext = fileparts(domainfile)
     100            if ext == '.exp':
     101                domain = expread(domainfile)
     102            elif ext == '.shp':
     103                domain = shpread(domainfile)
     104            else:
     105                raise Exception('bamg error message: file {} format not supported (.exp or .shp)'.format(domainfile))
    89106            domain = expread(domainfile)
     107        elif type(domainfile) == list:
     108            if len(domainfile):
     109                if type(domainfile[0]) in [dict, OrderedDict]:
     110                    domain = domainfile
     111                else:
     112                    raise Exception("bamg error message: if 'domain' is a list, its elements must be of type dict or OrderedDict")
     113            else:
     114                domain = domainfile
     115        elif type(domainfile) in [dict, OrderedDict]:
     116            domain = [domainfile]
    90117        else:
    91             domain = domainfile
    92 
    93         #Build geometry
     118            raise Exception("bamg error message: 'domain' type {} not supported yet".format(type(domainfile)))
     119
     120        holes = []
     121        if options.exist('holes'):
     122            holesfile = options.getfieldvalue('holes')
     123            if type(holesfile) == str:
     124                if not os.path.exists(holesfile):
     125                    raise IOError("bamg error message: file {} not found".format(holesfile))
     126
     127                # Read holes accoridng to its extension
     128                path, name, ext = fileparts(holesfile)
     129                if ext == '.exp':
     130                    holes = expread(holesfile)
     131                elif ext == '.shp':
     132                    holes = shpread(holesfile)
     133                else:
     134                    raise Exception('bamg error message: file {} format not supported (.exp or .shp)'.format(holesfile))
     135            elif type(holesfile) == list:
     136                if len(holesfile):
     137                    if type(holesfile[0]) in [dict, OrderedDict]:
     138                        holes = holesfile
     139                    else:
     140                        raise Exception("bamg error message: if 'holes' is a list, its elements must be of type dict or OrderedDict")
     141                else:
     142                    holes = holesfile
     143            elif type(holesfile) in [dict, OrderedDict]:
     144                holes = [holesfile]
     145            else:
     146                raise Exception("'holes' type {} not supported yet".format(type(holesfile)))
     147
     148        subdomains = []
     149        if options.exist('subdomains'):
     150            subdomainsfile = options.getfieldvalue('subdomains')
     151            if type(subdomainsfile) == str:
     152                if not os.path.exists(subdomainsfile):
     153                    raise IOError("bamg error message: file {} not found".format(subdomainsfile))
     154
     155                # Read subdomains accoridng to its extension
     156                path, name, ext = fileparts(subdomainsfile)
     157                if ext == '.exp':
     158                    subdomains = expread(subdomainsfile)
     159                elif ext == '.shp':
     160                    subdomains = shpread(subdomainsfile)
     161                else:
     162                    raise Exception('bamg error message: file {} format not supported (.exp or .shp)'.format(subdomainsfile))
     163            elif type(subdomainsfile) == list:
     164                if len(subdomainsfile):
     165                    if type(subdomainsfile[0]) in [dict, OrderedDict]:
     166                        subdomains = subdomainsfile
     167                    else:
     168                        raise Exception("bamg error message: if 'subdomains' is a list, its elements must be of type dict or OrderedDict")
     169                else:
     170                    subdomains = subdomainsfile
     171            elif type(subdomainsfile) in [dict, OrderedDict]:
     172                subdomains = [subdomainsfile]
     173            else:
     174                raise Exception("'subdomains' type {} not supported yet".format(type(subdomainsfile)))
     175
     176        # Build geometry
    94177        count = 0
    95         for i, domaini in enumerate(domain):
    96             #Check that the domain is closed
    97             if (domaini['x'][0] != domaini['x'][-1] or domaini['y'][0] != domaini['y'][-1]):
     178        for i in range(len(domain)):
     179            # Check that the domain is closed
     180            if (domain[i]['x'][0] != domain[i]['x'][-1] or domain[i]['y'][0] != domain[i]['y'][-1]):
    98181                raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed")
    99182
    100             #Checks that all holes are INSIDE the principle domain outline princial domain should be index 0
     183            # Check that all holes are INSIDE the principle domain outline
    101184            if i:
    102                 flags = ContourToNodes(domaini['x'], domaini['y'], domainfile, 0)[0]
     185                flags = ContourToNodes(domain[i]['x'], domain[i]['y'], [domain[0]], 0)[0] # NOTE: Down stack call to FetchPythonData requires contour to be a list of struct if we are not passing in a path to a file, hence the odd addressing: '[domain[0]]'
    103186                if np.any(np.logical_not(flags)):
    104187                    raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
    105188
    106             #Check orientation
    107             nods = domaini['nods'] - 1  #the domain are closed 1 = end
    108             test = np.sum((domaini['x'][1:nods + 1] - domaini['x'][0:nods]) * (domaini['y'][1:nods + 1] + domaini['y'][0:nods]))
     189            # Check orientation
     190            nods = domain[i]['nods'] - 1  # the domain is closed (domain[0] = domain[-1])
     191            test = np.sum((domain[i]['x'][1:nods + 1] - domain[i]['x'][0:nods]) * (domain[i]['y'][1:nods + 1] + domain[i]['y'][0:nods]))
    109192            if (i == 0 and test > 0) or (i > 0 and test < 0):
    110193                print('At least one contour was not correctly oriented and has been re-oriented')
    111                 domaini['x'] = np.flipud(domaini['x'])
    112                 domaini['y'] = np.flipud(domaini['y'])
    113 
    114             #Add all points to bamg_geometry
    115             nods = domaini['nods'] - 1  #the domain are closed 0 = end
    116             bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((domaini['x'][0:nods], domaini['y'][0:nods], np.ones((nods)))).T))
     194                domain[i]['x'] = np.flipud(domain[i]['x'])
     195                domain[i]['y'] = np.flipud(domain[i]['y'])
     196
     197            # Flag how many edges we have so far
     198            edge_length = len(bamg_geometry.Edges)
     199
     200            # Add all points to bamg_geometry
     201            #nods = domain[i]['nods'] - 1  #the domain are closed 0 = end
     202            bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((domain[i]['x'][0:nods], domain[i]['y'][0:nods], np.ones((nods)))).T))
    117203            bamg_geometry.Edges = np.vstack((bamg_geometry.Edges, np.vstack((np.arange(count + 1, count + nods + 1), np.hstack((np.arange(count + 2, count + nods + 1), count + 1)), 1. * np.ones((nods)))).T))
    118             if i:
     204
     205            # Flag how many edges we have now, that way we know which edges
     206            # belong to the subdomain. Will be used later fo required edges
     207            # if NoBoundaryRefinement equals 1.
     208            new_edge_length = len(bamg_geometry.Edges)
     209            edges_required = range((edge_length + 1), (new_edge_length + 1)) # NOTE: Upper bound of range is non-inclusive (compare to src/m/mesh/bamg.m)
     210
     211            if i: # NOTE: same as `if i > 0` (MATLAB is `if i > 1`)
    119212                bamg_geometry.SubDomains = np.vstack((bamg_geometry.SubDomains, [2, count + 1, 1, -subdomain_ref]))
    120213                subdomain_ref = subdomain_ref + 1
    … …  
    122215                bamg_geometry.SubDomains = np.vstack((bamg_geometry.SubDomains, [2, count + 1, 1, 0]))
    123216
    124             #update counter
     217            # Update counter
    125218            count += nods
    126219
    127         #Deal with domain holes
    128         if options.exist('holes'):
    129             holesfile = options.getfieldvalue('holes')
    130             if type(holesfile) == str:
    131                 if not os.path.exists(holesfile):
    132                     raise IOError("bamg error message: file '%s' not found" % holesfile)
    133                 holes = expread(holesfile)
    134             else:
    135                 holes = holesfile
    136 
    137             #Build geometry
    138             for i, holei in enumerate(holes):
    139                 #Check that the hole is closed
    140                 if (holei['x'][0] != holei['x'][-1] or holei['y'][0] != holei['y'][-1]):
    141                     raise RuntimeError("bamg error message: all contours provided in ''hole'' should be closed")
    142 
    143                 #Checks that all holes are INSIDE the principle domain outline princial domain should be index 0
    144                 flags = ContourToNodes(holei['x'], holei['y'], domainfile, 0)[0]
    145                 if np.any(np.logical_not(flags)):
    146                     raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
    147 
    148                 #Check orientation
    149                 nods = holei['nods'] - 1  #the hole are closed 1 = end
    150                 test = np.sum((holei['x'][1:nods + 1] - holei['x'][0:nods]) * (holei['y'][1:nods + 1] + holei['y'][0:nods]))
    151                 if test < 0:
    152                     print('At least one hole was not correctly oriented and has been re-oriented')
    153                     holei['x'] = np.flipud(holei['x'])
    154                     holei['y'] = np.flipud(holei['y'])
    155 
    156                 #Add all points to bamg_geometry
    157                 nods = holei['nods'] - 1  #the hole are closed 0 = end
    158                 bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((holei['x'][0:nods], holei['y'][0:nods], np.ones((nods)))).T))
    159                 bamg_geometry.Edges = np.vstack((bamg_geometry.Edges, np.vstack((np.arange(count + 1, count + nods + 1), np.hstack((np.arange(count + 2, count + nods + 1), count + 1)), 1. * np.ones((nods)))).T))
    160                 #update counter
    161                 count += nods
    162 
    163         #And subdomains
    164         if options.exist('subdomains'):
    165             subdomainfile = options.getfieldvalue('subdomains')
    166             if type(subdomainfile) == str:
    167                 if not os.path.exists(subdomainfile):
    168                     raise IOError("bamg error message: file '{}' not found".format(subdomainfile))
    169                 subdomains = expread(subdomainfile)
    170             else:
    171                 subdomains = subdomainfile
    172 
    173             #Build geometry
    174             for i, subdomaini in enumerate(subdomains):
    175                 #Check that the subdomain is closed
    176                 if (subdomaini['x'][0] != subdomaini['x'][-1] or subdomaini['y'][0] != subdomaini['y'][-1]):
    177                     raise RuntimeError("bamg error message: all contours provided in ''subdomain'' should be closed")
    178 
    179                 #Checks that all subdomains are INSIDE the principle subdomain outline princial domain should be index 0
    180                 if i:
    181                     flags = ContourToNodes(subdomaini['x'], subdomaini['y'], domainfile, 0)[0]
    182                     if np.any(np.logical_not(flags)):
    183                         raise RuntimeError("bamg error message: All subdomains should be strictly inside the principal subdomain")
    184 
    185                 #Check orientation
    186                 nods = subdomaini['nods'] - 1  #the subdomain are closed 1 = end
    187 
    188                 test = np.sum((subdomaini['x'][1:nods + 1] - subdomaini['x'][0:nods]) * (subdomaini['y'][1:nods + 1] + subdomaini['y'][0:nods]))
    189                 if test > 0:
    190                     print('At least one subcontour was not correctly oriented and has been re-oriented')
    191                     subdomaini['x'] = np.flipud(subdomaini['x'])
    192                     subdomaini['y'] = np.flipud(subdomaini['y'])
    193 
    194                 #Add all points to bamg_geometry
    195                 nods = subdomaini['nods'] - 1  #the subdomain are closed 0 = end
    196                 bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((subdomaini['x'][0:nods], subdomaini['y'][0:nods], np.ones((nods)))).T))
    197                 bamg_geometry.Edges = np.vstack((bamg_geometry.Edges, np.vstack((np.arange(count + 1, count + nods + 1), np.hstack((np.arange(count + 2, count + nods + 1), count + 1)), 1. * np.ones((nods)))).T))
    198                 #update counter
    199                 count += nods
     220        for i in range(len(holes)):
     221            # heck that the hole is closed
     222            if (holes[i]['x'][0] != holes[i]['x'][-1] or holes[i]['y'][0] != holes[i]['y'][-1]):
     223                raise RuntimeError("bamg error message: all contours provided in ''hole'' should be closed")
     224
     225            # Check that all holes are INSIDE the principal domain (principal domain should be index 0)
     226            flags = ContourToNodes(holes[i]['x'], holes[i]['y'], [domain[0]], 0)[0] # NOTE: Down stack call to FetchPythonData requires contour to be a list of struct if we are not passing in a path to a file, hence the odd addressing: '[domain[0]]'
     227            if np.any(np.logical_not(flags)):
     228                raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
     229
     230            # Check that hole is correctly oriented
     231            nods = holes[i]['nods'] - 1  # the holes are closed (holes[0] = holes[-1])
     232            test = np.sum((holes[i]['x'][1:nods + 1] - holes[i]['x'][0:nods]) * (holes[i]['y'][1:nods + 1] + holes[i]['y'][0:nods]))
     233            if test < 0:
     234                print('At least one hole was not correctly oriented and has been re-oriented')
     235                holes[i]['x'] = np.flipud(holes[i]['x'])
     236                holes[i]['y'] = np.flipud(holes[i]['y'])
     237
     238            # Add all points to bamg_geometry
     239            #nods = holes[i]['nods'] - 1  # the hole is closed (hole[0] = hole[-1])
     240            bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((holes[i]['x'][0:nods], holes[i]['y'][0:nods], np.ones((nods)))).T))
     241            bamg_geometry.Edges = np.vstack((bamg_geometry.Edges, np.vstack((np.arange(count + 1, count + nods + 1), np.hstack((np.arange(count + 2, count + nods + 1), count + 1)), 1. * np.ones((nods)))).T))
     242            bamg.geometry.SubDomains = np.vstack((bamg_geometry.SubDomains, [2, count + 1, 1, -hole_ref]))
     243            hole_ref = hole_ref + 1
     244
     245            # Update counter
     246            count += nods
     247
     248        for i in range(len(subdomains)):
     249            # Check that the subdomain is closed
     250            if (subdomains[i]['x'][0] != subdomains[i]['x'][-1] or subdomains[i]['y'][0] != subdomains[i]['y'][-1]):
     251                raise RuntimeError("bamg error message: all contours provided in ''subdomains'' should be closed")
     252
     253            # Check that all subdomains are INSIDE the principal domain (principal domain should be index 0)
     254            flags = ContourToNodes(subdomains[i]['x'], subdomains[i]['y'], [domain[0]], 0)[0] # NOTE: Down stack call to FetchPythonData requires contour to be a list of struct if we are not passing in a path to a file, hence the odd addressing: '[domain[0]]'
     255            if np.any(np.logical_not(flags)):
     256                raise RuntimeError("bamg error message: All subdomains should be strictly inside the principal domain")
     257
     258            # Check that subdomain is correctly oriented
     259            nods = subdomains[i]['nods'] - 1  # the subdomains are closed (subdomains[0] = subdomains[-1])
     260
     261            test = np.sum((subdomains[i]['x'][1:nods + 1] - subdomains[i]['x'][0:nods]) * (subdomains[i]['y'][1:nods + 1] + subdomains[i]['y'][0:nods]))
     262            if test:
     263                print('At least one subcontour was not correctly oriented and has been re-oriented')
     264                subdomains[i]['x'] = np.flipud(subdomains[i]['x'])
     265                subdomains[i]['y'] = np.flipud(subdomains[i]['y'])
     266
     267            #Add all points to bamg_geometry
     268            #nods = subdomains[i]['nods'] - 1  # the subdomains are closed (subdomains[0] = subdomains[-1])
     269            bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((subdomains[i]['x'][0:nods], subdomains[i]['y'][0:nods], np.ones((nods)))).T))
     270            bamg_geometry.Edges = np.vstack((bamg_geometry.Edges, np.vstack((np.arange(count + 1, count + nods + 1), np.hstack((np.arange(count + 2, count + nods + 1), count + 1)), 1. * np.ones((nods)))).T))
     271            bamg_geometry.SubDomains = np.vstack((bamg_geometry.SubDomains, [2, count + 1, 1, subdomain_ref]))
     272            subdomain_ref = subdomain_ref + 1
     273
     274            # Update counter
     275            count += nods
    200276
    201277        if options.getfieldvalue('vertical', 0):
    202278            if np.size(options.getfieldvalue('Markers', [])) != np.size(bamg_geometry.Edges, 0):
    203279                raise RuntimeError('for 2d vertical mesh, ''Markers'' option is required, and should be of size ' + str(np.size(bamg_geometry.Edges, 0)))
    204             if np.size(options.getfieldvalue('Markers', [])) == np.size(bamg_geometry.Edges, 0):
    205                 bamg_geometry.Edges[:, 2] = options.getfieldvalue('Markers')
    206 
    207         #take care of rifts
     280        if np.size(options.getfieldvalue('Markers', [])) == np.size(bamg_geometry.Edges, 0):
     281            bamg_geometry.Edges[:, 2] = options.getfieldvalue('Markers')
     282
     283        # Take care of rifts
    208284        if options.exist('rifts'):
    209 
    210             #Check that file exists
     285            # Read rift file according to its extension
    211286            riftfile = options.getfieldvalue('rifts')
    212             if not os.path.exists(riftfile):
    213                 raise IOError("bamg error message: file '%s' not found" % riftfile)
    214             rift = expread(riftfile)
    215 
    216             for i, rifti in enumerate(rift):
    217                 #detect whether all points of the rift are inside the domain
    218                 flags = ContourToNodes(rifti['x'], rifti['y'], domain[0], 0)[0]
     287            path, name, ext = fileparts(riftfile)
     288            if ext == '.exp':
     289                rift = expread(riftfile)
     290            elif ext == '.shp':
     291                rift = shpread(riftfile)
     292            else:
     293                raise IOError("bamg error message: file '{}' format not supported (.exp or .shp)".format(riftfile))
     294
     295            for i in range(len(rift)):
     296                # Detect whether all points of the rift are inside the domain
     297                flags = ContourToNodes(rift[i]['x'], rift[i]['y'], [domain[0]], 0)[0] # NOTE: Down stack call to FetchPythonData requires contour to be a list of struct if we are not passing in a path to a file, hence the odd addressing: '[domain[0]]'
    219298                if np.all(np.logical_not(flags)):
    220299                    raise RuntimeError("one rift has all its points outside of the domain outline")
    221 
    222300                elif np.any(np.logical_not(flags)):
    223                     #We LOTS of work to do
     301                    # We have LOTS of work to do
    224302                    print("Rift tip outside of or on the domain has been detected and is being processed...")
    225                     #check that only one point is outside (for now)
     303
     304                    # Check that only one point is outside (for now)
    226305                    if np.sum(np.logical_not(flags).astype(int)) != 1:
    227                         raise RuntimeError("bamg error message: only one point outside of the domain is supported yet")
    228 
    229                     #Move tip outside to the first position
     306                        raise RuntimeError("bamg error message: only one point outside of the domain is supported at this time")
     307
     308                    # Move tip outside to the first position
    230309                    if not flags[0]:
    231                         #OK, first point is outside (do nothing),
     310                        # OK, first point is outside (do nothing),
    232311                        pass
    233312                    elif not flags[-1]:
    234                         rifti['x'] = np.flipud(rifti['x'])
    235                         rifti['y'] = np.flipud(rifti['y'])
     313                        rift[i]['x'] = np.flipud(rift[i]['x'])
     314                        rift[i]['y'] = np.flipud(rift[i]['y'])
    236315                    else:
    237316                        raise RuntimeError("bamg error message: only a rift tip can be outside of the domain")
    238317
    239                     #Get cordinate of intersection point
    240                     x1 = rifti['x'][0]
    241                     y1 = rifti['y'][0]
    242                     x2 = rifti['x'][1]
    243                     y2 = rifti['y'][1]
     318                    # Get coordinate of intersection point
     319                    x1 = rift[i]['x'][0]
     320                    y1 = rift[i]['y'][0]
     321                    x2 = rift[i]['x'][1]
     322                    y2 = rift[i]['y'][1]
    244323                    for j in range(0, np.size(domain[0]['x']) - 1):
    245324                        if SegIntersect(np.array([[x1, y1], [x2, y2]]), np.array([[domain[0]['x'][j], domain[0]['y'][j]], [domain[0]['x'][j + 1], domain[0]['y'][j + 1]]])):
    246325
    247                             #Get position of the two nodes of the edge in domain
     326                            # Get position of the two nodes of the edge in domain
    248327                            i1 = j
    249328                            i2 = j + 1
    250329
    251                             #rift is crossing edge [i1 i2] of the domain
    252                             #Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
     330                            # Rift is crossing edge [i1, i2] of the domain
     331                            # Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
    253332                            x3 = domain[0]['x'][i1]
    254333                            y3 = domain[0]['y'][i1]
    255334                            x4 = domain[0]['x'][i2]
    256335                            y4 = domain[0]['y'][i2]
    257                             #x = det([det([x1 y1; x2 y2])  x1 - x2;det([x3 y3; x4 y4])  x3 - x4]) / det([x1 - x2 y1 - y2;x3 - x4 y3 - y4])
    258                             #y = det([det([x1 y1; x2 y2])  y1 - y2;det([x3 y3; x4 y4])  y3 - y4]) / det([x1 - x2 y1 - y2;x3 - x4 y3 - y4])
    259336                            x = np.linalg.det(np.array([[np.linalg.det(np.array([[x1, y1], [x2, y2]])), x1 - x2], [np.linalg.det(np.array([[x3, y3], [x4, y4]])), x3 - x4]])) / np.linalg.det(np.array([[x1 - x2, y1 - y2], [x3 - x4, y3 - y4]]))
    260337                            y = np.linalg.det(np.array([[np.linalg.det(np.array([[x1, y1], [x2, y2]])), y1 - y2], [np.linalg.det(np.array([[x3, y3], [x4, y4]])), y3 - y4]])) / np.linalg.det(np.array([[x1 - x2, y1 - y2], [x3 - x4, y3 - y4]]))
    … …  
    264341
    265342                            if np.min(tipdis) / segdis < options.getfieldvalue('toltip', 0):
    266                                 print("moving tip - domain intersection point")
    267 
    268                                 #Get position of the closer point
     343                                print("moving tip-domain intersection point")
     344
     345                                # Get position of the closer point
    269346                                if tipdis[0] > tipdis[1]:
    270347                                    pos = i2
    … …  
    272349                                    pos = i1
    273350
    274                                 #This point is only in Vertices (number pos).
    275                                 #OK, now we can add our own rift
    276                                 nods = rifti['nods'] - 1
    277                                 bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.hstack((rifti['x'][1:].reshape(-1, ), rifti['y'][1:].reshape(-1, ), np.ones((nods, 1))))))
    278                                 bamg_geometry.Edges = np.vstack((bamg_geometry.Edges,
    279                                                                  np.array([[pos, count + 1, (1 + i)]]),
    280                                                                  np.hstack((np.arange(count + 1, count + nods).reshape(-1, ), np.arange(count + 2, count + nods + 1).reshape(-1, ), (1 + i) * np.ones((nods - 1, 1))))))
     351                                # This point is only in Vertices (number pos).
     352                                # OK, now we can add our own rift
     353                                nods = rift[i]['nods'] - 1
     354                                bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.hstack((rift[i]['x'][1:].reshape(-1, ), rift[i]['y'][1:].reshape(-1, ), np.ones((nods, 1))))))
     355                                bamg_geometry.Edges = np.vstack((
     356                                    bamg_geometry.Edges,
     357                                    np.array([[pos, count + 1, (1 + i)]]),
     358                                    np.hstack((np.arange(count + 1, count + nods).reshape(-1, ), np.arange(count + 2, count + nods + 1).reshape(-1, ), (1 + i) * np.ones((nods - 1, 1))))
     359                                ))
    281360                                count += nods
    282361                                break
    283 
    284362                            else:
    285                                 #Add intersection point to Vertices
    286                                 bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.array([[x, y, 1]])))
     363                                # Add intersection point to Vertices
     364                                bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices,
     365                                    np.array([[x, y, 1]])
     366                                ))
    287367                                count += 1
    288368
    289                                 #Decompose the crossing edge into 2 subedges
     369                                # Decompose the crossing edge into 2 subedges
    290370                                pos = np.nonzero(np.logical_and(bamg_geometry.Edges[:, 0] == i1, bamg_geometry.Edges[:, 1] == i2))[0]
    291371                                if not pos:
    292372                                    raise RuntimeError("bamg error message: a problem occurred...")
    293                                 bamg_geometry.Edges = np.vstack((bamg_geometry.Edges[0:pos - 1, :],
    294                                                                  np.array([[bamg_geometry.Edges[pos, 0], count, bamg_geometry.Edges[pos, 2]]]),
    295                                                                  np.array([[count, bamg_geometry.Edges[pos, 1], bamg_geometry.Edges[pos, 2]]]),
    296                                                                  bamg_geometry.Edges[pos + 1:, :]))
    297 
    298                                 #OK, now we can add our own rift
    299                                 nods = rifti['nods'] - 1
    300                                 bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.hstack((rifti['x'][1:].reshape(-1, ), rifti['y'][1:].reshape(-1, ), np.ones((nods, 1))))))
    301                                 bamg_geometry.Edges = np.vstack((bamg_geometry.Edges,
    302                                                                  np.array([[count, count + 1, 2]]),
    303                                                                  np.hstack((np.arange(count + 1, count + nods).reshape(-1, ), np.arange(count + 2, count + nods + 1).reshape(-1, ), (1 + i) * np.ones((nods - 1, 1))))))
     373                                bamg_geometry.Edges = np.vstack((
     374                                    bamg_geometry.Edges[0:pos - 1, :],
     375                                    np.array([[
     376                                        bamg_geometry.Edges[pos, 0],
     377                                        count,
     378                                        bamg_geometry.Edges[pos, 2]
     379                                    ]]),
     380                                    np.array([[
     381                                        count,
     382                                        bamg_geometry.Edges[pos, 1],
     383                                        bamg_geometry.Edges[pos, 2]
     384                                    ]]),
     385                                    bamg_geometry.Edges[pos + 1:, :]
     386                                ))
     387
     388                                # OK, now we can add our own rift
     389                                nods = rift[i]['nods'] - 1
     390                                bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices,
     391                                    np.hstack((
     392                                        rift[i]['x'][1:].reshape(-1, ),
     393                                        rift[i]['y'][1:].reshape(-1, ),
     394                                        np.ones((nods, 1))
     395                                    ))
     396                                ))
     397                                bamg_geometry.Edges = np.vstack((
     398                                    bamg_geometry.Edges,
     399                                    np.array([[count, count + 1, 2]]),
     400                                    np.hstack((
     401                                        np.arange(count + 1, count + nods).reshape(-1, ),
     402                                        np.arange(count + 2, count + nods + 1).reshape(-1, ),
     403                                        (1 + i) * np.ones((nods - 1, 1))
     404                                    ))
     405                                ))
    304406                                count += nods
    305407                                break
    306 
    307408                else:
    308                     nods = rifti['nods'] - 1
    309                     bamg_geometry.Vertices = np.vstack(bamg_geometry.Vertices, np.hstack(rifti['x'][:], rifti['y'][:], np.ones((nods + 1, 1))))
    310                     bamg_geometry.Edges = np.vstack(bamg_geometry.Edges, np.hstack(np.arange(count + 1, count + nods).reshape(-1, ), np.arange(count + 2, count + nods + 1).reshape(-1, ), i * np.ones((nods, 1))))
    311                     count += nods + 1
    312 
    313         #Deal with tracks
     409                    nods = rift[i]['nods'] - 1
     410                    bamg_geometry.Vertices = np.vstack((
     411                        bamg_geometry.Vertices,
     412                        np.hstack((
     413                            rift[i]['x'][:],
     414                            rift[i]['y'][:],
     415                            np.ones((nods + 1, 1))
     416                        ))
     417                    ))
     418                    bamg_geometry.Edges = np.vstack((
     419                        bamg_geometry.Edges,
     420                        np.hstack((
     421                            np.arange(count + 1, count + nods).reshape(-1, ),
     422                            np.arange(count + 2, count + nods + 1).reshape(-1, ),
     423                            i * np.ones((nods, 1))
     424                        ))
     425                    ))
     426                    count += (nods + 1)
     427
     428        # Deal with tracks
    314429        if options.exist('tracks'):
    315 
    316             #read tracks
     430            # Read tracks
    317431            track = options.getfieldvalue('tracks')
    318432            if all(isinstance(track, str)):
    … …  
    320434                track = np.hstack((A.x.reshape(-1, ), A.y.reshape(-1, )))
    321435            else:
    322                 track = float(track)  #for some reason, it is of class "single"
     436                track = float(track)
     437
    323438            if np.size(track, axis=1) == 2:
    324439                track = np.hstack((track, 3. * np.ones((size(track, axis=0), 1))))
    325440
    326             #only keep those inside
    327             flags = ContourToNodes(track[:, 0], track[:, 1], domainfile, 0)[0]
     441            # Only keep those inside
     442            flags = ContourToNodes(track[:, 0], track[:, 1], [domain[0]], 0)[0] # NOTE: Down stack call to FetchPythonData requires contour to be a list of struct if we are not passing in a path to a file, hence the odd addressing: '[domain[0]]'
    328443            track = track[np.nonzero(flags), :]
    329444
    330             #Add all points to bamg_geometry
     445            # Add all points to bamg_geometry
    331446            nods = np.size(track, axis=0)
    332447            bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, track))
    333             bamg_geometry.Edges = np.vstack((bamg_geometry.Edges, np.hstack((np.arange(count + 1, count + nods).reshape(-1, ), np.arange(count + 2, count + nods + 1).reshape(-1, ), 3. * np.ones((nods - 1, 1))))))
    334 
    335             #update counter
     448            bamg_geometry.Edges = np.vstack((
     449                bamg_geometry.Edges,
     450                np.hstack((
     451                    np.arange(count + 1, count + nods).reshape(-1, ),
     452                    np.arange(count + 2, count + nods + 1).reshape(-1, ),
     453                    3. * np.ones((nods - 1, 1))
     454                ))
     455            ))
     456
     457            # Update counter
    336458            count += nods
    337459
    338         #Deal with vertices that need to be kept by mesher
     460        # Deal with vertices that need to be kept by mesher
    339461        if options.exist('RequiredVertices'):
    340 
    341             #recover RequiredVertices
    342             requiredvertices = options.getfieldvalue('RequiredVertices')  #for some reason, it is of class "single"
     462            # Recover RequiredVertices
     463            requiredvertices = options.getfieldvalue('RequiredVertices')
    343464            if np.size(requiredvertices, axis=1) == 2:
    344465                requiredvertices = np.hstack((requiredvertices, 4. * np.ones((np.size(requiredvertices, axis=0), 1))))
    345466
    346             #only keep those inside
    347             flags = ContourToNodes(requiredvertices[:, 0], requiredvertices[:, 1], domainfile, 0)[0]
     467            # Only keep those inside
     468            flags = ContourToNodes(requiredvertices[:, 0], requiredvertices[:, 1], [domain[0]], 0)[0] # NOTE: Down stack call to FetchPythonData requires contour to be a list of struct if we are not passing in a path to a file, hence the odd addressing: '[domain[0]]'
    348469            requiredvertices = requiredvertices[np.nonzero(flags)[0], :]
    349             #Add all points to bamg_geometry
     470
     471            # Add all points to bamg_geometry
    350472            nods = np.size(requiredvertices, axis=0)
    351473            bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, requiredvertices))
    352474
    353             #update counter
     475            # Update counter
    354476            count += nods
    355477
    356     #process geom
    357     #bamg_geometry = processgeometry(bamg_geometry, options.getfieldvalue('tol', float(nan)), domain[0])
     478        # Deal with RequiredEdges
     479        if options.getfieldvalue('NoBoundaryRefinement', 0):
     480            bamg_geometry.RequiredEdges = edges_required
     481        elif options.getfieldvalue('NoBoundaryRefinementAllBoundaries', 0):
     482            bamg_geometry.RequiredEdges = np.arange(1, bamg_geometry.Edges.shape[0]).T
     483
     484        # Process geom
     485        #bamg_geometry = processgeometry(bamg_geometry, options.getfieldvalue('tol', np.nan), domain[0])
    358486    elif isinstance(md.private.bamg, dict) and 'geometry' in md.private.bamg:
    359487        bamg_geometry = bamggeom(md.private.bamg['geometry'].__dict__)
    … …  
    362490        pass
    363491    #}}}
    364     # Bamg Mesh parameters {{{
    365     if not options.exist('domain') and md.mesh.numberofvertices and m.strcmp(md.mesh.elementtype(), 'Tria'):
     492    # Bamg mesh parameters {{{
     493    if not options.exist('domain') and md.mesh.numberofvertices and md.mesh.elementtype() == 'Tria':
    366494        if isinstance(md.private.bamg, dict) and 'mesh' in md.private.bamg:
    367495            bamg_mesh = bamgmesh(md.private.bamg['mesh'].__dict__)
    368496        else:
    369             bamg_mesh.Vertices = np.vstack((md.mesh.x, md.mesh.y, np.ones((md.mesh.numberofvertices)))).T
    370             #bamg_mesh.Vertices = np.hstack((md.mesh.x.reshape(-1, ), md.mesh.y.reshape(-1, ), np.ones((md.mesh.numberofvertices, 1))))
     497            bamg_mesh.Vertices = np.vstack((
     498                md.mesh.x,
     499                md.mesh.y,
     500                np.ones((md.mesh.numberofvertices))
     501            )).T
    371502            bamg_mesh.Triangles = np.hstack((md.mesh.elements, np.ones((md.mesh.numberofelements, 1))))
    372503
    … …  
    374505            raise TypeError("bamg error message: rifts not supported yet. Do meshprocessrift AFTER bamg")
    375506    #}}}
    376     # Bamg Options {{{
     507    # Bamg options {{{
    377508    bamg_options['Crack'] = options.getfieldvalue('Crack', 0)
    378509    bamg_options['anisomax'] = options.getfieldvalue('anisomax', 10.**18)
    … …  
    402533    #}}}
    403534
    404     #call Bamg
     535    # Call Bamg
    405536    bamgmesh_out, bamggeom_out = BamgMesher(bamg_mesh.__dict__, bamg_geometry.__dict__, bamg_options)
    406537
    407     # plug results onto model
     538    # Plug results onto model
    408539    if options.getfieldvalue('vertical', 0):
    409540        md.mesh = mesh2dvertical()
    … …  
    415546        md.mesh.segmentmarkers = bamgmesh_out['IssmSegments'][:, 3].astype(int)
    416547
    417         #Fill in rest of fields:
     548        # Fill in rest of fields
    418549        md.mesh.numberofelements = np.size(md.mesh.elements, axis=0)
    419550        md.mesh.numberofvertices = np.size(md.mesh.x)
    420551        md.mesh.numberofedges = np.size(md.mesh.edges, axis=0)
    421         md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool)
    422         md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True
    423 
    424         #Now, build the connectivity tables for this mesh. Doubled in matlab for some reason
    425         md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, )
     552        md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, int)
    426553        md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1
    427554
    … …  
    437564        md.mesh.segmentmarkers = bamgmesh_out['IssmSegments'][:, 3].astype(int)
    438565
    439         #Fill in rest of fields:
     566        # Fill in rest of fields
    440567        md.mesh.numberofelements = np.size(md.mesh.elements, axis=0)
    441568        md.mesh.numberofvertices = np.size(md.mesh.x)
    442569        md.mesh.numberofedges = np.size(md.mesh.edges, axis=0)
    443         md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool)
    444         md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True
     570        md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, int)
     571        md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1
    445572
    446573    else:
    … …  
    453580        md.mesh.segmentmarkers = bamgmesh_out['IssmSegments'][:, 3].astype(int)
    454581
    455         #Fill in rest of fields:
     582        # Fill in rest of fields
    456583        md.mesh.numberofelements = np.size(md.mesh.elements, axis=0)
    457584        md.mesh.numberofvertices = np.size(md.mesh.x)
    458585        md.mesh.numberofedges = np.size(md.mesh.edges, axis=0)
    459         md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool)
    460         md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True
    461 
    462     #Bamg private fields
     586        md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, int)
     587        md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1
     588
     589    # Bamg private fields
    463590    md.private.bamg = OrderedDict()
    464591    md.private.bamg['mesh'] = bamgmesh(bamgmesh_out)
    … …  
    476603
    477604def processgeometry(geom, tol, outline):  # {{{
    478     raise RuntimeError("bamg.py / processgeometry is not complete.")
    479     #Deal with edges
     605    raise RuntimeError("bamg.py::processgeometry is not complete.")
     606
     607    # Deal with edges
    480608    print("Checking Edge crossing...")
    481609    i = 0
    482610    while (i < np.size(geom.Edges, axis=0)):
    483         #edge counter
     611        # Edge counter
    484612        i += 1
    485613
    486         #Get coordinates
     614        # Get coordinates
    487615        x1 = geom.Vertices[geom.Edges[i, 0], 0]
    488616        y1 = geom.Vertices[geom.Edges[i, 0], 1]
    … …  
    491619        color1 = geom.Edges[i, 2]
    492620
    493         j = i  #test edges located AFTER i only
     621        j = i # Test edges located AFTER i only
    494622        while (j < np.size(geom.Edges, axis=0)):
    495             #edge counter
     623            # Edge counter
    496624            j += 1
    497625
    498             #Skip if the two edges already have a vertex in common
     626            # Skip if the two edges already have a vertex in common
    499627            if any(m.ismember(geom.Edges[i, 0:2], geom.Edges[j, 0:2])):
    500628                continue
    501629
    502             #Get coordinates
     630            # Get coordinates
    503631            x3 = geom.Vertices[geom.Edges[j, 0], 0]
    504632            y3 = geom.Vertices[geom.Edges[j, 0], 1]
    … …  
    507635            color2 = geom.Edges[j, 2]
    508636
    509             #Check if the two edges are crossing one another
     637            # Check if the two edges are crossing one another
    510638            if SegIntersect(np.array([[x1, y1], [x2, y2]]), np.array([[x3, y3], [x4, y4]])):
    511639
    512                 #Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
     640                # Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
    513641                x = np.linalg.det(np.array([np.linalg.det(np.array([[x1, y1], [x2, y2]])), x1 - x2], [np.linalg.det(np.array([[x3, y3], [x4, y4]])), x3 - x4]) / np.linalg.det(np.array([[x1 - x2, y1 - y2], [x3 - x4, y3 - y4]])))
    514642                y = np.linalg.det(np.array([np.linalg.det(np.array([[x1, y1], [x2, y2]])), y1 - y2], [np.linalg.det(np.array([[x3, y3], [x4, y4]])), y3 - y4]) / np.linalg.det(np.array([[x1 - x2, y1 - y2], [x3 - x4, y3 - y4]])))
    515643
    516                 #Add vertex to the list of vertices
     644                # Add vertex to the list of vertices
    517645                geom.Vertices = np.vstack((geom.Vertices, [x, y, min(color1, color2)]))
    518646                id = np.size(geom.Vertices, axis=0)
    519647
    520                 #Update edges i and j
     648                # Update edges i and j
    521649                edgei = geom.Edges[i, :].copy()
    522650                edgej = geom.Edges[j, :].copy()
    … …  
    526654                geom.Edges = np.vstack((geom.Edges, [id, edgej(1), edgej(2)]))
    527655
    528                 #update current edge second tip
     656                # Update current edge second tip
    529657                x2 = x
    530658                y2 = y
    531659
    532     #Check point outside
     660    # Check point outside
    533661    print("Checking for points outside the domain...")
    534662    i = 0
    535663    num = 0
    536664    while (i < np.size(geom.Vertices, axis=0)):
    537         #vertex counter
     665        # Vertex counter
    538666        i += 1
    539667
    540         #Get coordinates
     668        # Get coordinates
    541669        x = geom.Vertices[i, 0]
    542670        y = geom.Vertices[i, 1]
    543671        color = geom.Vertices[i, 2]
    544672
    545         #Check that the point is inside the domain
     673        # Check that the point is inside the domain
    546674        if color != 1 and not ContourToNodes(x, y, outline[0], 1):
    547             #Remove points from list of Vertices
     675            # Remove points from list of Vertices
    548676            num += 1
    549677            geom.Vertices[i, :] = []
    550678
    551             #update edges
     679            # Update edges
    552680            posedges = np.nonzero(geom.Edges == i)
    553681            geom.Edges[posedges[0], :] = []
    … …  
    555683            geom.Edges[posedges] = geom.Edges[posedges] - 1
    556684
    557             #update counter
     685            # Update counter
    558686            i -= 1
    559687
    … …  
    564692    %Check point spacing
    565693    if ~isnan(tol),
    566             disp('Checking point spacing...')
     694            print('Checking point spacing...')
    567695            i = 0
    568696            while (i < size(geom.Vertices, 1)),
    … …  
    608736    """
    609737    return geom
    610     # }}}
     738    #}}}

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

    r24260 r25455  
    77
    88def bamgflowband(md, x, surf, base, *args):
    9     """
    10     BAMGFLOWBAND - create flowband mesh with bamg
     9    """BAMGFLOWBAND - create flowband mesh with bamg
    1110
    1211    Usage:

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

    r13009 r25455  
    22%FINDSEGMENTS - build segments model field
    33%
     4%   Usage:
     5%      segments=findsegments(md,varargin);
     6%
    47%   Optional inputs:
    58%      'mesh.elementconnectivity'
    6 %
    7 %   Usage:
    8 %      segments=findsegments(md,varargin);
    99
    1010%get options
    … …  
    1414mesh.elementconnectivity=getfieldvalue(options,'mesh.elementconnectivity',md.mesh.elementconnectivity);
    1515
    16 %Now, build the connectivity tables for this mesh if not correclty done
     16%Now, build the connectivity tables for this mesh if not correctly done
    1717if size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements,
    1818        if exist(options,'mesh.elementconnectivity'),
    19                 error(' ''mesh.elementconnectivity'' option does not have thge right size.');
     19                error('''mesh.elementconnectivity'' option does not have thge right size.');
    2020        else
    2121                mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
    … …  
    3939        els2=mesh.elementconnectivity(el1,find(mesh.elementconnectivity(el1,:)));
    4040
    41         %el1 is connected to 2 other elements
     41        %get nodes of 'el1'
     42        nods1=md.mesh.elements(el1,:);
     43
     44        %'el1' is connected to 2 other elements
    4245        if length(els2)>1,
    43 
    44                 %get nodes of el1
    45                 nods1=md.mesh.elements(el1,:);
    4646
    4747                %find the common vertices to the two elements connected to el1 (1 or 2)
    … …  
    5555                ord1=find(nods1(1)==md.mesh.elements(el1,:));
    5656                ord2=find(nods1(2)==md.mesh.elements(el1,:));
     57
    5758                if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
    5859                        temp=segments(count,1);
    … …  
    6061                        segments(count,2)=temp;
    6162                end
     63                disp(segments(count,1:2))
    6264                segments(count,1:2)=fliplr(segments(count,1:2));
     65                disp(segments)
    6366                count=count+1;
    6467
    65         %el1 is connected to only one element
     68        %'el1' is connected to only one element
    6669        else
    67                 %get nodes of el1
    68                 nods1=md.mesh.elements(el1,:);
    6970
    70                 %find the vertex  the el1 to not share with els2
     71                %find the vertex the 'el1' does not share with 'els2'
    7172                flag=setdiff(nods1,md.mesh.elements(els2,:));
    7273

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

    r24213 r25455  
    11import numpy as np
     2
     3from BamgConvertMesh import BamgConvertMesh
     4from bamggeom import bamggeom
     5from bamgmesh import bamgmesh
    26from collections import OrderedDict
    3 from BamgConvertMesh import BamgConvertMesh
    47from mesh2d import mesh2d
    5 from bamgmesh import bamgmesh
    6 from bamggeom import bamggeom
    78
    89
    910def meshconvert(md, *args):
    10     """
    11     CONVERTMESH - convert mesh to bamg mesh
     11    """CONVERTMESH - convert mesh to bamg mesh
    1212
    13        Usage:
    14           md = meshconvert(md)
    15           md = meshconvert(md, index, x, y)
     13    Usage:
     14        md = meshconvert(md)
     15        md = meshconvert(md, index, x, y)
    1616    """
    1717
    18     if not len(args) == 0 and not len(args) == 3:
     18    nargs = len(args)
     19
     20    if not nargs == 0 and not nargs == 3:
    1921        raise TypeError("meshconvert error message: bad usage")
    2022
    21     if not len(args):
     23    if not nargs:
    2224        index = md.mesh.elements
    2325        x = md.mesh.x
    … …  
    4749    md.mesh.numberofvertices = np.size(md.mesh.x)
    4850    md.mesh.numberofedges = np.size(md.mesh.edges, axis=0)
    49     md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool)
    50     md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True
     51    md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, int)
     52    md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1
    5153
    5254    return md

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

    r25000 r25455  
    88        options=pairoptions(varargin{:});
    99       
    10         tolerance=getfieldvalue(options,'tolerance',10^-5);
     10        tolerance=getfieldvalue(options,'tolerance',1e-4);
    1111       
    1212        md=model();
    … …  
    1616        md.private=md1.private;
    1717
    18         %some initializatoin:
     18        %some initialization:
    1919        elements1=md1.mesh.elements;
    2020        x1=md1.mesh.x;
    … …  
    3434        elements2=elements2+nods1;
    3535
    36         %go into the vertices on boundary of mesh 1, and figure out which ones are common to mesh2:
    37         verticesonboundary=find(md1.mesh.vertexonboundary);
     36        %go into the vertices on boundary of mesh 1 and figure out which ones are common with mesh2:
     37        verticesonboundary=find(md1.mesh.vertexonboundary);
     38
    3839        for i=1:length(verticesonboundary),
    39                 node1=verticesonboundary(i); xnode1=x1(node1); ynode1=y1(node1); znode1=z1(node1);
    40                 %is there another node with these coordinates in mesh2?
    41                 ind=find(sqrt(((x2-xnode1).^2+(y2-ynode1).^2)+(z2-znode1).^2)<tolerance);
     40                node1=verticesonboundary(i);
     41                xnode1=x1(node1);
     42                ynode1=y1(node1);
     43                znode1=z1(node1);
     44
     45                %is there another node with these coordinates in mesh 2?
     46                ind=find(sqrt((x2-xnode1).^2+(y2-ynode1).^2+(z2-znode1).^2)<tolerance);
    4247                if length(ind)>1,
    4348                        disp('should reduce the tolerance, several vertices picked up!');
    … …  
    4752                        y2(ind)=NaN;
    4853                        z2(ind)=NaN;
    49                         pos=find(elements2==(ind+nods1)); elements2(pos)=node1;
     54                        pos=find(elements2==(ind+nods1));
     55                        elements2(pos)=node1;
    5056                end
    5157        end
     58        %go through elements2 and drop counter on each vertex that is above the x2 and y2 vertices being dropped:
     59        indices_nan=isnan(x2);
     60        while(~isempty(indices_nan)),
     61                % Use the index of the first instance of 'nan' value to remove that element from 'x2', 'y2', and 'z2'
     62                index_nan=indices_nan(1);
     63                pos=find(elements2>(index_nan+nods1));
     64                elements2(pos)=elements2(pos)-1;
     65                x2(index_nan)=[];
     66                y2(index_nan)=[];
     67                z2(index_nan)=[];
    5268
    53         %go through elements2 and drop counter on each vertex that is above the x2 and y2 vertices being dropped:
    54         while( ~isempty(find(isnan(x2)))),
    55                 for i=1:length(x2),
    56                         if isnan(x2(i)),
    57                                 pos=find(elements2>(i+nods1));
    58                                 elements2(pos)=elements2(pos)-1;
    59                                 x2(i)=[];
    60                                 y2(i)=[];
    61                                 z2(i)=[];
    62                                 break;
    63                         end
    64                 end
     69                % Check again in 'x2' for instances of 'nan'
     70                indices_nan=find(isnan(x2));
    6571        end
    6672
    … …  
    8389        %connectivities:
    8490        md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
     91        % disp(md.mesh.vertexconnectivity)
    8592        md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
     93        % disp(md.mesh.elementconnectivity)
    8694
    8795        %find segments:
    … …  
    9199        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
    92100        md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
     101        disp(md.mesh.vertexonboundary)
     102        pause
    93103
    94104        %some checks:

  • issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py

    r24213 r25455  
    11import numpy as np
     2
     3from ContourToMesh import ContourToMesh
    24from ElementsFromEdge import ElementsFromEdge
    35import MatlabFuncs as m
    4 from ContourToMesh import ContourToMesh
    56
    67
    78def meshprocessoutsiderifts(md, domainoutline):
    8     """
    9     MESHPROCESSOUTSIDERIFTS - process rifts when they touch the domain outline
     9    """MESHPROCESSOUTSIDERIFTS - process rifts when they touch the domain outline
    1010
    11        Usage:
    12           md = meshprocessoutsiderifts(md, domain)
    13 
     11    Usage:
     12        md = meshprocessoutsiderifts(md, domain)
    1413    """
    1514
    … …  
    5352                elements = np.concatenate((elements, nextelement))
    5453                #new B:
    55                 B = md.mesh.elements[nextelement - 1, np.nonzero(np.logical_not(m.ismember(md.mesh.elements[nextelement - 1, :], np.array([A, B]))))]
     54                B = md.mesh.elements[nextelement - 1, np.nonzero(np.logical_not(m.ismember(md.mesh.elements[nextelement - 1, :], np.array([A, B]))))[0]]
    5655
    5756            #take the list of elements on one side of the rift that connect to the tip,
    … …  
    6463            #replace tip in elements
    6564            newelements = md.mesh.elements[elements - 1, :]
    66             pos = np.nonzero(newelements == tip)
     65            pos = np.nonzero(newelements == tip)[0]
    6766            newelements[pos] = num
    6867            md.mesh.elements[elements - 1, :] = newelements
    … …  
    8180    md.mesh.numberofelements = np.size(md.mesh.elements, axis=0)
    8281    md.mesh.numberofvertices = np.size(md.mesh.x)
    83     md.mesh.vertexonboundary = np.zeros(np.size(md.mesh.x), bool)
    84     md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True
     82    md.mesh.vertexonboundary = np.zeros(np.size(md.mesh.x), int)
     83    md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1
    8584    md.rifts.numrifts = np.length(md.rifts.riftstruct)
    8685
    … …  
    8887
    8988
    90 def isconnected(elements, A, B):  # {{{
    91     """
    92     ISCONNECTED: are two nodes connected by a triangulation?
     89def isconnected(elements, A, B):  #{{{
     90    """ISCONNECTED: are two nodes connected by a triangulation?
    9391
    94        Usage: flag = isconnected(elements, A, B)
    95 
     92    Usage:
     93        flag = isconnected(elements, A, B)
    9694    """
    9795
    … …  
    103101
    104102    return flag
    105     # }}}
     103    #}}}

  • issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py

    r24256 r25455  
    11import numpy as np
     2
     3from ContourToMesh import ContourToMesh
     4from GetAreas import GetAreas
     5from meshprocessoutsiderifts import meshprocessoutsiderifts
    26from ProcessRifts import ProcessRifts
    3 from ContourToMesh import ContourToMesh
    4 from meshprocessoutsiderifts import meshprocessoutsiderifts
    5 from GetAreas import GetAreas
    67
    78
    89def meshprocessrifts(md, domainoutline):
    9     """
    10     MESHPROCESSRIFTS - process mesh when rifts are present
     10    """MESHPROCESSRIFTS - process mesh when rifts are present
    1111
    12        split rifts inside mesh (rifts are defined by presence of
    13        segments inside the domain outline)
    14        if domain outline is provided, check for rifts that could touch it, and open them up.
     12    split rifts inside mesh (rifts are defined by presence of
     13    segments inside the domain outline)
     14    if domain outline is provided, check for rifts that could touch it, and open them up.
    1515
    16        Usage:
    17           md = meshprocessrifts(md, domainoutline)
     16    Usage:
     17        md = meshprocessrifts(md, domainoutline)
    1818
    19        Ex:
    20           md = meshprocessrifts(md, 'DomainOutline.exp')
    21 
     19    Example:
     20        md = meshprocessrifts(md, 'DomainOutline.exp')
    2221    """
    2322
    24     #Call MEX file
     23    # Call Python module
    2524    md.mesh.elements, md.mesh.x, md.mesh.y, md.mesh.segments, md.mesh.segmentmarkers, md.rifts.riftstruct = ProcessRifts(md.mesh.elements, md.mesh.x, md.mesh.y, md.mesh.segments, md.mesh.segmentmarkers)
    2625    md.mesh.elements = md.mesh.elements.astype(int)
    … …  
    3534    md.mesh.numberofelements = np.size(md.mesh.elements, axis=0)
    3635    md.mesh.numberofvertices = np.size(md.mesh.x)
    37     md.mesh.vertexonboundary = np.zeros(np.size(md.mesh.x), bool)
    38     md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True
     36    md.mesh.vertexonboundary = np.zeros(np.size(md.mesh.x), int)
     37    md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1
    3938
    4039    #get coordinates of rift tips

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

    r24213 r25455  
    11import numpy as np
    2 from NodeConnectivity import NodeConnectivity
     2
    33from ElementConnectivity import ElementConnectivity
    44from mesh2d import mesh2d
     5from NodeConnectivity import NodeConnectivity
    56
    67
    78def squaremesh(md, Lx, Ly, nx, ny):
    8     """
    9     SQUAREMESH - create a structured square mesh
     9    """SQUAREMESH - create a structured square mesh
    1010
    11        This script will generate a structured square mesh
    12        Lx and Ly are the dimension of the domain (in meters)
    13        nx anx ny are the number of nodes in the x and y direction
    14        The coordinates x and y returned are in meters.
     11    This script will generate a structured square mesh
     12    Lx and Ly are the dimension of the domain (in meters)
     13    nx anx ny are the number of nodes in the x and y direction
     14    The coordinates x and y returned are in meters.
    1515
    16        Usage:
    17           [md] = squaremesh(md, Lx, Ly, nx, ny)
     16    Usage:
     17        [md] = squaremesh(md, Lx, Ly, nx, ny)
    1818    """
    1919
    … …  
    2424    #initialization
    2525    index = np.zeros((nel, 3), int)
    26     x = np.zeros((nx * ny))
    27     y = np.zeros((nx * ny))
     26    x = np.zeros((nx * ny, ))
     27    y = np.zeros((nx * ny, ))
    2828
    2929    #create coordinates
    … …  
    6363    md.mesh.y = y
    6464    md.mesh.numberofvertices = nods
    65     md.mesh.vertexonboundary = np.zeros((nods), bool)
    66     md.mesh.vertexonboundary[segments[:, 0:2] - 1] = True
     65    md.mesh.vertexonboundary = np.zeros(nods, int)
     66    md.mesh.vertexonboundary[segments[:, 0:2] - 1] = 1
    6767
    6868    #plug elements
    … …  
    7272
    7373    #Now, build the connectivity tables for this mesh.
    74     md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)[0]
    75     md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)[0]
     74    md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
     75    md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
    7676
    7777    return md

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

    r24213 r25455  
    11import os.path
     2
    23import numpy as np
     4
     5from ElementConnectivity import ElementConnectivity
    36from mesh2d import mesh2d
    47from NodeConnectivity import NodeConnectivity
    5 from ElementConnectivity import ElementConnectivity
    68from Triangle_python import Triangle_python
    79
    810
    911def triangle(md, domainname, *args):
    10     """
    11     TRIANGLE - create model mesh using the triangle package
     12    """TRIANGLE - create model mesh using the triangle package
    1213
    13        This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
    14        where md is a @model object, domainname is the name of an Argus domain outline file,
    15        and resolution is a characteristic length for the mesh (same unit as the domain outline
    16        unit). Riftname is an optional argument (Argus domain outline) describing rifts.
     14    This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
     15    where md is a @model object, domainname is the name of an Argus domain outline file,
     16    and resolution is a characteristic length for the mesh (same unit as the domain outline
     17    unit). Riftname is an optional argument (Argus domain outline) describing rifts.
    1718
    18        Usage:
    19           md = triangle(md, domainname, resolution)
    20        or md = triangle(md, domainname, resolution, riftname)
     19    Usage:
     20        md = triangle(md, domainname, resolution)
     21        OR
     22        md = triangle(md, domainname, resolution, riftname)
    2123
    22        Examples:
    23           md = triangle(md, 'DomainOutline.exp', 1000)
    24           md = triangle(md, 'DomainOutline.exp', 1000, 'Rifts.exp')
     24    Examples:
     25        md = triangle(md, 'DomainOutline.exp', 1000)
     26        md = triangle(md, 'DomainOutline.exp', 1000, 'Rifts.exp')
    2527    """
    2628
    … …  
    4547            return None
    4648
    47     area = resolution**2
     49    area = resolution ** 2
    4850
    4951    #Check that file exist (this is a very very common mistake)
    … …  
    6163    md.mesh.numberofelements = np.size(md.mesh.elements, axis=0)
    6264    md.mesh.numberofvertices = np.size(md.mesh.x)
    63     md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool)
    64     md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True
     65    md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, int)
     66    md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1
    6567
    6668    #Now, build the connectivity tables for this mesh.
    67     md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)[0]
    68     md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)[0]
     69    md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
     70    md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
    6971
    7072    return md

  • issm/trunk-jpl/src/m/modules/BamgMesher.py

    r24213 r25455  
    33
    44def BamgMesher(bamgmesh, bamggeom, bamgoptions):
    5     """
    6     BAMGMESHER
    7 
    8     Usage:
    9             bamgmesh, bamggeom = BamgMesher(bamgmesh, bamggeom, bamgoptions)
     5    """BAMGMESHER
    106
    117    bamgmesh: input bamg mesh
    128    bamggeom: input bamg geometry for the mesh
    139    bamgoptions: options for the bamg mesh
     10
     11    Usage:
     12        bamgmesh, bamggeom = BamgMesher(bamgmesh, bamggeom, bamgoptions)
    1413    """
    1514
    16     #Call mex module
     15    # Call module
    1716    bamgmesh, bamggeom = BamgMesher_python(bamgmesh, bamggeom, bamgoptions)
    1817
    19     #return
    2018    return bamgmesh, bamggeom

  • issm/trunk-jpl/src/m/modules/ContourToNodes.py

    r24213 r25455  
    33
    44def ContourToNodes(x, y, contourname, edgevalue):
    5     """
    6     CONTOURTONODES - flags vertices inside contour
     5    """CONTOURTONODES - flags vertices inside contour
    76
    8         Usage:
    9             flags = ContourToNodes(x, y, contourname, edgevalue)
     7    x, y:           list of nodes
     8    contourname:    name of .exp/.shp file containing the contours, or resulting structure from call to expread/shpread
     9    edgevalue:      integer (0, 1 or 2) defining the value associated to the nodes on the edges of the polygons
     10    flags:          vector of flags (0 or 1), of size nodes
    1011
    11         x, y: list of nodes
    12         contourname: name of .exp file containing the contours, or resulting structure from call to expread
    13         edgevalue: integer (0, 1 or 2) defining the value associated to the nodes on the edges of the polygons
    14         flags: vector of flags (0 or 1), of size nodes
     12    Usage:
     13        flags = ContourToNodes(x, y, contourname, edgevalue)
    1514    """
    1615
    17     #Call mex module
     16    # Call Python module
    1817    flags = ContourToNodes_python(x, y, contourname, edgevalue)
    1918
    20     #return
    2119    return flags

  • issm/trunk-jpl/src/m/modules/ElementConnectivity.py

    r24213 r25455  
    33
    44def ElementConnectivity(elements, nodeconnectivity):
     5    """ELEMENTCONNECTIVITY - Build element connectivity using node connectivity and elements
     6
     7    Usage:
     8        elementconnectivity = ElementConnectivity(elements, nodeconnectivity)
    59    """
    6     ELEMENTCONNECTIVITY - Build element connectivity using node connectivity and elements
    7 
    8         Usage:
    9             elementconnectivity = ElementConnectivity(elements, nodeconnectivity)
    10     """
    11     #Call mex module
     10   
     11    # Call Python module
    1212    elementconnectivity = ElementConnectivity_python(elements, nodeconnectivity)
    1313
    14     #Return
    15     return elementconnectivity
     14    return elementconnectivity[0] # NOTE: Value returned from wrapper function is a tuple, the first element of which being the result we actually want

  • issm/trunk-jpl/src/m/modules/ExpToLevelSet.py

    r25163 r25455  
    22from shpread import shpread
    33
     4
    45def ExpToLevelSet(x, y, contourname): #{{{
    5     '''
    6     EXPTOLEVELSET - Determine levelset distance between a contour and a cloud
    7     of points
     6    """EXPTOLEVELSET - Determine levelset distance between a contour and a
     7    cloud of points
    88
    9         Usage:
    10             distance = ExpToLevelSet(x, y, contourname)
     9    Usage:
     10        distance = ExpToLevelSet(x, y, contourname)
    1111
    12         x, y:           cloud point
    13         contourname:    name of .exp file containing the contours
    14         distance:       distance vector representing a levelset where the 0
    15                         level is one of the contour segments
     12    x, y:           cloud point
     13    contourname:    name of .exp file containing the contours
     14    distance:       distance vector representing a levelset where the 0
     15                    level is one of the contour segments
    1616
    17         Example:
    18             distance = ExpToLevelSet(md.mesh.x, md.mesh.y, 'Contour.exp')
     17    Example:
     18        distance = ExpToLevelSet(md.mesh.x, md.mesh.y, 'Contour.exp')
    1919
    20         TODO:
    21         - Need to compile Python version of ExpToLevelSet_matlab for this
    22         to work as intended (see src/m/modules/ExpToLevelSet.m)
    23     '''
     20    TODO:
     21    - Need to compile Python version of ExpToLevelSet_matlab for this
     22    to work as intended (see src/m/modules/ExpToLevelSet.m)
     23    """
    2424
    2525    if isinstance(contourname, basestring):

  • issm/trunk-jpl/src/m/modules/InterpFromMesh2d.m

    r20875 r25455  
    44%   Usage:
    55%      data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime);
    6 %      or data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value);
    7 %      or data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value,contourname);
     6%      OR
     7%      data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value);
     8%      OR
     9%      data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value,contourname);
    810%
     11%   index:      index of the mesh where data is defined
    912%   x,y:        coordinates of the nodes where data is defined
    10 %   index:      index of the mesh where data is defined
    1113%   data:       vector holding the data to be interpolated onto the points
    1214%   x_prime,y_prime:    coordinates of the mesh vertices onto which we interpolate
    … …  
    3234                data_prime=InterpFromMesh2d_matlab(varargin{1},varargin{2},varargin{3},varargin{4},varargin{5},varargin{6},varargin{7},varargin{8});
    3335        otherwise
     36                % NOTE: Should never get here because of previous check
    3437                error('InterpFromMesh2d not supported');
    3538end

  • issm/trunk-jpl/src/m/modules/InterpFromMeshToMesh2d.py

    r24213 r25455  
    33
    44def InterpFromMeshToMesh2d(*args):
    5     """
    6     INTERPFROMMESHTOMESH2D - Interpolation from a 2d triangular mesh onto a list of points
     5    """INTERPFROMMESHTOMESH2D - Interpolation from a 2d triangular mesh onto a list of points
    76
    8             Usage:
    9                     data_interp = InterpFromMeshToMesh2d(index, x, y, data, x_interp, y_interp)
    10                     or data_interp = InterpFromMeshToMesh2d(index, x, y, data, x_interp, y_interp, OPTIONS)
     7    Usage:
     8        data_interp = InterpFromMeshToMesh2d(index, x, y, data, x_interp, y_interp)
     9        OR
     10        data_interp = InterpFromMeshToMesh2d(index, x, y, data, x_interp, y_interp, OPTIONS)
    1111
    12             index:  index of the mesh where data is defined (e.g. md.mesh.elements)
    13             x, y:    coordinates of the nodes where data is defined
    14             data:   matrix holding the data to be interpolated onto the mesh (one column per field)
    15             x_interp, y_interp:      coordinates of the points onto which we interpolate
    16             data_interp:    vector of mesh interpolated data
    17             Available options:
    18                     default:        default value if point is outsite of triangulation (instead of linear interpolation)
     12    index:              index of the mesh where data is defined (e.g. md.mesh.elements)
     13    x, y:               coordinates of the nodes where data is defined
     14    data:               matrix holding the data to be interpolated onto the mesh (one column per field)
     15    x_interp, y_interp: coordinates of the points onto which we interpolate
     16    data_interp:        vector of mesh interpolated data
     17    Available options:
     18            default:    default value if point is outsite of triangulation (instead of linear interpolation)
    1919
    2020    Example:
    21             load('temperature.mat')
    22             md.initialization.temperature = InterpFromMeshToMesh2d(index, x, y, temperature, md.mesh.x, md.mesh.y)
    23             md.initialization.temperature = InterpFromMeshToMesh2d(index, x, y, temperature, md.mesh.x, md.mesh.y, 'default', 253)
     21        load('temperature.mat')
     22        md.initialization.temperature = InterpFromMeshToMesh2d(index, x, y, temperature, md.mesh.x, md.mesh.y)
     23        md.initialization.temperature = InterpFromMeshToMesh2d(index, x, y, temperature, md.mesh.x, md.mesh.y, 'default', 253)
    2424    """
    25     # Call mex module
     25
     26    # Check usage
     27    nargs = len(args)
     28    if nargs != 6 and nargs != 8:
     29        print(InterpFromMeshToMesh2d.__doc__)
     30        raise Exception('Wrong usage (see above)')
     31
     32    # Call Python module
    2633    data_interp = InterpFromMeshToMesh2d_python(*args)
    27     # Return
    28     return data_interp
     34
     35    return data_interp[0] # NOTE: Value returned from wrapper function is a tuple, the first element of which being the result we actually want

  • issm/trunk-jpl/src/m/modules/InterpFromMeshToMesh3d.py

    r24213 r25455  
    33
    44def InterpFromMeshToMesh3d(index, x, y, z, data, x_prime, y_prime, z_prime, default_value):
     5    """INTERPFROMMESHTOMESH3D - Interpolation from a 3d hexahedron mesh onto a list of points
     6
     7    Usage:
     8        index:                      index of the mesh where data is defined
     9        x, y, z:                    coordinates of the nodes where data is defined
     10        data:                       matrix holding the data to be interpolated onto the mesh
     11        x_prime, y_prime, z_prime:  coordinates of the points onto which we interpolate
     12        default_value:              default value if no data is found (holes)
     13        data_prime:                 vector of mesh interpolated data
     14
     15    Example:
     16        load('temperature.mat')
     17        md.initialization.temperature = InterpFromMeshToMesh3d(index, x, y, z, temperature, md.mesh.x, md.mesh.y, md.mesh.z, 253)
    518    """
    6     INTERPFROMMESHTOMESH3D - Interpolation from a 3d hexahedron mesh onto a list of points
    719
    8    Usage:
    9       index:    index of the mesh where data is defined
    10       x, y, z:    coordinates of the nodes where data is defined
    11       data:    matrix holding the data to be interpolated onto the mesh
    12       x_prime, y_prime, z_prime:    coordinates of the points onto which we interpolate
    13       default_value:    default value if no data is found (holes)
    14       data_prime:    vector of mesh interpolated data
     20    # Check usage
     21    nargs = len(args)
     22    if nargs != 9:
     23        print(InterpFromMeshToMesh3d.__doc__)
     24        raise Exception('Wrong usage (see above)')
    1525
    16    Example:
    17       load('temperature.mat')
    18       md.initialization.temperature = InterpFromMeshToMesh3d(index, x, y, z, temperature, md.mesh.x, md.mesh.y, md.mesh.z, 253)
    19     """
    20     # Call mex module
     26    # Call Python module
    2127    data_prime = InterpFromMeshToMesh3d_python(index, x, y, z, data, x_prime, y_prime, z_prime, default_value)
    2228
    23     # Return
    2429    return data_prime

  • issm/trunk-jpl/src/m/modules/NodeConnectivity.py

    r24213 r25455  
    33
    44def NodeConnectivity(elements, numnodes):
     5    """NODECONNECTIVITY - Build node connectivity from elements
     6
     7    Usage:
     8        connectivity = NodeConnectivity(elements, numnodes)
    59    """
    6     NODECONNECTIVITY - Build node connectivity from elements
    710
    8         Usage:
    9             connectivity = NodeConnectivity(elements, numnodes)
    10     """
    11     # Call mex module
     11    # Call Python module
    1212    connectivity = NodeConnectivity_python(elements, numnodes)
    1313
    14     # Return
    15     return connectivity
     14    return connectivity[0] # NOTE: Value returned from wrapper function is a tuple, the first element of which being the result we actually want

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

    r24213 r25455  
    99
    1010def contourenvelope(md, *args):
    11     """
    12     CONTOURENVELOPE - build a set of segments enveloping a contour .exp
     11    """CONTOURENVELOPE - build a set of segments enveloping a contour .exp
    1312
    14        Usage:
    15           segments = contourenvelope(md, varargin)
     13    Usage:
     14        segments = contourenvelope(md, varargin)
    1615
    17        Example:
    18           segments = contourenvelope(md, 'Stream.exp')
    19           segments = contourenvelope(md)
     16    Example:
     17        segments = contourenvelope(md, 'Stream.exp')
     18        segments = contourenvelope(md)
    2019    """
    2120
    … …  
    4140    #Computing connectivity
    4241    if np.size(md.mesh.vertexconnectivity, axis=0) != md.mesh.numberofvertices and np.size(md.mesh.vertexconnectivity, axis=0) != md.mesh.numberofvertices2d:
    43         md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)[0]
     42        md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
    4443    if np.size(md.mesh.elementconnectivity, axis=0) != md.mesh.numberofelements and np.size(md.mesh.elementconnectivity, axis=0) != md.mesh.numberofelements2d:
    45         md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)[0]
     44        md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
    4645
    4746    #get nodes inside profile

  • issm/trunk-jpl/src/m/partition/adjacency.py

    r24254 r25455  
    11import numpy as np
     2
     3from GetAreas import *
    24import MatlabFuncs as m
    35from NodeConnectivity import *
    4 from GetAreas import *
    56
    67
    78def adjacency(md):
    8     #ADJACENCY - compute adjacency matrix, list of vertices and list of weights.
    9     #
    10     #  function to create the adjacency matrix from the connectivity table.
    11     #
    12     #  the required output is:
    13     #    md.adj_mat     (double [sparse nv x nv], vertex adjacency matrix)
    14     #    md.qmu.vertex_weight        (double [nv], vertex weights)
     9    """ADJACENCY - compute adjacency matrix, list of vertices and list of weights.
     10
     11    function to create the adjacency matrix from the connectivity table.
     12   
     13    the required output is:
     14        md.adj_mat     (double [sparse nv x nv], vertex adjacency matrix)
     15        md.qmu.vertex_weight        (double [nv], vertex weights)
     16    """
    1517
    1618    indi = np.array([md.mesh.elements[:, 0], md.mesh.elements[:, 1], md.mesh.elements[:, 2]])

  • issm/trunk-jpl/src/m/plot/processmesh.py

    r25163 r25455  
    44
    55def processmesh(md, data, options):
    6     '''
    7     PROCESSMESH - process mesh to be plotted
     6    """PROCESSMESH - process mesh to be plotted
    87
    98    Usage:
    … …  
    1615    - Test that output of delaunay matches output of delaunay in MATLAB (see
    1716    src/m/plot/processmesh.m)
    18     '''
     17    """
    1918
    2019    #some checks

  • issm/trunk-jpl/src/m/qmu/dakota_in_data.py

    r25163 r25455  
    5757
    5858    if strcmpi(params.analysis_driver, 'matlab') and params.analysis_components == '':
    59         [pathstr, name, ext] = fileparts(filei)
     59        pathstr, name, ext = fileparts(filei)
    6060        params.analysis_components = fullfile(pathstr, name + '.py')
    6161

  • issm/trunk-jpl/src/m/qmu/dakota_in_write.py

    r24870 r25455  
    5656        filei = str(eval(input('Dakota input file to write?  ')))
    5757
    58     [pathstr, name, ext] = fileparts(filei)
     58    pathstr, name, ext = fileparts(filei)
    5959    if len(ext) == 0:
    6060        # fileparts only considers '.in' to be the extension, not '.qmu.in'
    … …  
    266266            param_write(fidi, '\t  ', 'processors_per_evaluation', '=', '\n', params)
    267267        if len(params.analysis_components) != 0:
    268             [pathstr, name, ext] = fileparts(params.analysis_components)
     268            pathstr, name, ext = fileparts(params.analysis_components)
    269269            if ext != '':
    270270                ext = '.py'

  • issm/trunk-jpl/src/m/qmu/helpers.py

    r25065 r25455  
    66
    77class struct(object):
    8     '''An empty struct that can be assigned arbitrary attributes'''
     8    """An empty struct that can be assigned arbitrary attributes"""
    99    pass
    1010
    1111
    1212class Lstruct(list):
    13     '''
     13    """
    1414    An empty struct that can be assigned arbitrary attributes but can also be
    1515    accesed as a list. Eg. x.y = 'hello', x[:] = ['w', 'o', 'r', 'l', 'd']
    … …  
    4646    Sources:
    4747    -https://github.com/Vectorized/Python-Attribute-List
    48     '''
     48    """
    4949
    5050    def __new__(self, *args, **kwargs):
    … …  
    6464
    6565class OrderedStruct(object):
    66     '''
     66    """
    6767    A form of dictionary-like structure that maintains the ordering in which
    6868    its fields/attributes and their corresponding values were added.
    … …  
    113113    Note: to access internal fields use dir(x) (input fields will be included,
    114114    but are not technically internals)
    115     '''
     115    """
    116116
    117117    def __init__(self, *args):
    118         '''
     118        """
    119119        Provided either nothing or a series of strings, construct a structure
    120120        that will, when accessed as a list, return its fields in the same order
    121121        in which they were provided
    122         '''
     122        """
    123123
    124124        # keys and values
    … …  
    187187
    188188    def __copy__(self):
    189         '''
     189        """
    190190        shallow copy, hard copies of trivial attributes,
    191191        references to structures like lists/OrderedDicts
    192192        unless redefined as an entirely different structure
    193         '''
     193        """
    194194        newInstance = type(self)()
    195195        for k, v in list(self.items()):
    … …  
    198198
    199199    def __deepcopy__(self, memo=None):
    200         '''
     200        """
    201201        hard copy of all attributes
    202202        same thing but call deepcopy recursively
    … …  
    204204        (see https://docs.python.org/2/library/copy.html  #copy.deepcopy )
    205205        but will generally work in this case
    206         '''
     206        """
    207207        newInstance = type(self)()
    208208        for k, v in list(self.items()):
    … …  
    232232
    233233def isempty(x):
    234     '''
     234    """
    235235    returns true if object is +/-infinity, NaN, None, '', has length 0, or is
    236236    an array/matrix composed only of such components (includes mixtures of
    237237    "empty" types)
    238     '''
     238    """
    239239
    240240    if type(x) in [list, np.ndarray, tuple]:
    … …  
    270270
    271271def fieldnames(x, ignore_internals=True):
    272     '''
     272    """
    273273    returns a list of fields of x
    274274    ignore_internals ignores all fieldnames starting with '_' and is True by
    275275    default
    276     '''
     276    """
    277277    result = list(vars(x).keys())
    278278
    … …  
    284284
    285285def isfield(x, y, ignore_internals=True):
    286     '''
     286    """
    287287    is y is a field of x?
    288288    ignore_internals ignores all fieldnames starting with '_' and is True by
    289289    default
    290     '''
     290    """
    291291    return str(y) in fieldnames(x, ignore_internals)
    292292
    293293
    294294def fileparts(x):
    295     '''
     295    """
    296296    given:   "path/path/.../file_name.ext"
    297297    returns: [path, file_name, ext] (list of strings)
    298     '''
     298    """
    299299    try:
    300300        a = x[:x.rindex('/')]  #path
    … …  
    311311
    312312def fullfile(*args):
    313     '''
     313    """
    314314    usage:
    315315        fullfile(path, path, ... , file_name + ext)
    … …  
    322322        as final arguments ('file.doc') or ('file' + '.doc') will work
    323323        ('final', '.doc'), and the like, will not (you'd get 'final/.doc')
    324     '''
     324    """
    325325    result = str(args[0])
    326326    for i in range(len(args[1:])):
    … …  
    334334
    335335def findline(fidi, s):
    336     '''
     336    """
    337337    returns full first line containing s (as a string), or None
    338338
    … …  
    340340    use str(findline(f, s)).strip() to remove these, str() in case result is
    341341    None
    342     '''
     342    """
    343343    for line in fidi:
    344344        if s in line:
    … …  
    348348
    349349def empty_nd_list(shape, filler=0., as_numpy_ndarray=False):
    350     '''
     350    """
    351351    returns a python list of the size/shape given (shape must be int or tuple)
    352352    the list will be filled with the optional second argument
    … …  
    363363        empty_nd_list((5, 5), 0.0)  # returns a 5x5 matrix of 0.0's
    364364        empty_nd_list(5, None)  # returns a 5 long array of NaN
    365     '''
     365    """
    366366    result = np.empty(shape)
    367367    result.fill(filler)

  • issm/trunk-jpl/src/m/qmu/postqmu.py

    r25163 r25455  
    11from os import getpid, stat
    22from os.path import isfile
    3 import shlex
    43from subprocess import call
    54
    … …  
    6261    # move all the individual function evalutations into zip files
    6362    if not md.qmu.isdakota:
    64         subproc_args = shlex.split('zip -mq params.in.zip params.in.[1-9]*')
     63        subproc_args = 'zip -mq params.in.zip params.in.[1-9]*'
    6564        call(subproc_args, shell=True)
    66         subproc_args = shlex.split('zip -mq results.out.zip results.out.[1-9]*')
     65        subproc_args = 'zip -mq results.out.zip results.out.[1-9]*'
    6766        call(subproc_args, shell=True)
    68         subproc_args = shlex.split('zip -mq matlab.out.zip matlab*.out.[1-9]*')
     67        subproc_args = 'zip -mq matlab.out.zip matlab*.out.[1-9]*'
    6968        call(subproc_args, shell=True)
    7069

  • issm/trunk-jpl/src/m/shp/shpread.m

    r25125 r25455  
    22%SHPREAD - read a shape file and build a struct
    33%
    4 %       This routine reads a shape file .shp and builds a structure array
    5 %       containing the fields x and y corresponding to the coordinates, one for the
    6 %       filename of the shp file, for the density, for the nodes, and a field
    7 %       closed to indicate if the domain is closed. If this initial shapefile is
    8 %       point only, the fields closed and points are omitted.
    9 %       The first argument is the .shp file to be read and the second one
    10 %       (optional) indicates if the last point shall be read (1 to read it, 0 not
    11 %       to).
     4%   This routine reads a shape file .shp and builds a structure array
     5%   containing the fields x and y corresponding to the coordinates, one for the
     6%   filename of the shp file, for the density, for the nodes, and a field
     7%   closed to indicate if the domain is closed. If this initial shapefile is
     8%   point only, the fields closed and points are omitted.
    129%
    13 %       Usage:
    14 %               Struct=shpread(filename)
     10%   The first argument is the .shp file to be read and the second one
     11%   (optional) indicates if the last point shall be read (1 to read it, 0 not
     12%   to).
    1513%
    16 %       Example:
    17 %               Struct=shpread('domainoutline.shp')
     14%   Usage:
     15%      Struct=shpread(filename)
    1816%
    19 %       See also EXPDOC, EXPWRITEASVERTICES
     17%   Example:
     18%      Struct=shpread('domainoutline.shp')
    2019
    2120%recover options

  • issm/trunk-jpl/src/m/shp/shpread.py

    r25125 r25455  
     1from collections import OrderedDict
    12from os import path
    23try:
    … …  
    56    print("could not import shapefile, PyShp has not been installed, no shapefile reading capabilities enabled")
    67
    7 from helpers import OrderedStruct
    88from pairoptions import pairoptions
    99
    1010
    1111def shpread(filename, *args): #{{{
    12     '''
    13     SHPREAD - read a shape file and build a dict
     12    """SHPREAD - read a shapefile and build a list of shapes
    1413
    15     This routine reads a shape file .shp and builds a dict array containing the
    16     fields x and y corresponding to the coordinates, one for the filename
    17     of the shp file, for the density, for the nodes, and a field closed to
    18     indicate if the domain is closed. If this initial shapefile is point only,
    19     the fields closed and points are ommited.
    20     The first argument is the .shp file to be read and the second one
    21     (optional) indicates if the last point shall be read (1 to read it, 0 not
    22     to).
     14    This routine reads a shapefile and builds a list of OrderedDict objects
     15    containing the fields x and y corresponding to the coordinates, one for the
     16    filename of the shp file, for the density, for the nodes, and a field
     17    closed to indicate if the domain is closed. If this initial shapefile is
     18    point only, the fields closed and points are ommitted. The first argument
     19    is the shapefile to be read and the second argument (optional) indicates if
     20    the last point shall be read (1 to read it, 0 not to).
    2321
    2422    The current implementation of shpread depends on PyShp.
    2523
    2624    Usage:
    27         dict = shpread(filename)
     25        list = shpread(filename)
    2826
    2927    Example:
    … …  
    3129        a collection of three files. You specify the base filename of the
    3230        shapefile or the complete filename of any of the shapefile component
    33         files.""
     31        files."
    3432
    35         dict = shpread('domainoutline.shp')
     33        list = shpread('domainoutline.shp')
    3634        OR
    37         dict = shpread('domainoutline.dbf')
     35        list = shpread('domainoutline.dbf')
    3836        OR
    39         dict = shpread('domainoutline')
     37        list = shpread('domainoutline')
    4038
    4139        "OR any of the other 5+ formats which are potentially part of a
    … …  
    4442        .shp extension exists.
    4543
    46     See also EXPDOC, EXPWRITEASVERTICES
    47 
    4844    Sources:
    4945    - https://github.com/GeospatialPython/pyshp
    5046
     47    NOTE:
     48    - OrderedDict objects are used instead of OrderedStruct objects (although
     49    addressing in the latter case is closer to the MATLAB struct type) in order
     50    to remain consistent with the pattern established by src/m/exp/expread.py.
     51
    5152    TODO:
    5253    - Create class that can be used to store and pretty print shape structs
    53       (ala OrderedStruct from src/m/qmu/helpers.py).
    54     '''
     54    (ala OrderedStruct from src/m/qmu/helpers.py).
     55    - Convert returned data structure from list of OrderedDict objects to list
     56    of OrderedStruct objects and remove corresponding note (see also
     57    src/m/exp/expread.py). Also, modify handling of returned data structure in,
     58        - src/m/classes/basin.py
     59        - src/m/classes/boundary.py
     60        - src/m/modules/ExpToLevelSet.py
     61        - src/m/mesh/bamg.py
     62    May also need to modify addressing in corresponding FetchData function, or
     63    create new one, in src/wrappers/ContoursToNodes/ContoursToNodes.cpp.
     64    """
    5565
    5666    #recover options
    … …  
    6777    shapes = sf.shapes()
    6878    for i in range(len(shapes)):
    69         Struct = OrderedStruct()
     79        Struct = OrderedDict()
    7080        shape = shapes[i]
    7181        if shape.shapeType == shapefile.POINT:
    72             Struct.x = shape.points[0][0]
    73             Struct.y = shape.points[0][1]
    74             Struct.density = 1
    75             Struct.Geometry = 'Point'
     82            Struct['x'] = shape.points[0][0]
     83            Struct['y'] = shape.points[0][1]
     84            Struct['density'] = 1
     85            Struct['Geometry'] = 'Point'
    7686        elif shape.shapeType == shapefile.POLYLINE:
    7787            num_points = len(shape.points)
    … …  
    8292                x.append(point[0])
    8393                y.append(point[1])
    84             Struct.x = x
    85             Struct.y = y
    86             Struct.nods = num_points
    87             Struct.density = 1
    88             Struct.closed = 1
    89             Struct.BoundingBox = shape.bbox
    90             Struct.Geometry = 'Line'
     94            Struct['x'] = x
     95            Struct['y'] = y
     96            Struct['nods'] = num_points
     97            Struct['density'] = 1
     98            Struct['closed'] = 1
     99            Struct['BoundingBox'] = shape.bbox
     100            Struct['Geometry'] = 'Line'
    91101        elif shape.shapeType == shapefile.POLYGON:
    92102            num_points = len(shape.points)
    … …  
    97107                x.append(point[0])
    98108                y.append(point[1])
    99             Struct.x = x
    100             Struct.y = y
    101             Struct.nods = num_points
    102             Struct.density = 1
    103             Struct.closed = 1
    104             Struct.BoundingBox = shape.bbox
    105             Struct.Geometry = 'Polygon'
     109            Struct['x'] = x
     110            Struct['y'] = y
     111            Struct['nods'] = num_points
     112            Struct['density'] = 1
     113            Struct['closed'] = 1
     114            Struct['BoundingBox'] = shape.bbox
     115            Struct['Geometry'] = 'Polygon'
    106116        else:
    107117            # NOTE: We could do this once before looping over shapes as all
    … …  
    122132            else:
    123133                setattr(Struct, str(fieldname), sf.record(i)[j - 1]) # cast to string removes "u" from "u'fieldname'"
    124         Struct.name = name
     134        Struct['name'] = name
    125135        Structs.append(Struct)
    126136
    … …  
    128138    if invert:
    129139        for i in range(len(Structs)):
    130             Structs[i].x = np.flipud(Structs[i].x)
    131             Structs[i].y = np.flipud(Structs[i].y)
     140            Structs[i]['x'] = np.flipud(Structs[i]['x'])
     141            Structs[i]['y'] = np.flipud(Structs[i]['y'])
    132142
    133143    return Structs

  • issm/trunk-jpl/src/wrappers/ContourToMesh/ContourToMesh.cpp

    r20855 r25455  
    6767
    6868        /*Run interpolation routine: */
    69         ContourToMeshx( &in_nod,&in_elem,index,x,y,contours,interptype,nel,nods,edgevalue);
     69        ContourToMeshx(&in_nod,&in_elem,index,x,y,contours,interptype,nel,nods,edgevalue);
    7070
    7171        /* output: */

  • issm/trunk-jpl/src/wrappers/InterpFromMesh2d/InterpFromMesh2d.cpp

    r22731 r25455  
    4242        /*contours*/
    4343        int i;
     44        #ifdef _HAVE_MATLAB_MODULES_
    4445        mxArray *matlabstructure = NULL;
     46        #endif
    4547        Contour<double> **contours=NULL;
    4648        int numcontours;
    … …  
    5961
    6062        /*checks on arguments on the matlab side: */
     63        #ifdef _HAVE_MATLAB_MODULES_
    6164        if(nlhs!=NLHS){
    6265                InterpFromMesh2dUsage();
    6366                _error_("InterpFromMeshToMesh2dUsage usage error");
    6467        }
     68        #endif
    6569        if((nrhs!=6) && (nrhs!=7) && (nrhs!=8)){
    6670                InterpFromMesh2dUsage();
    … …  
    8589        }
    8690
    87         if(nrhs>=8){
     91        if(nrhs==8){
    8892
    8993                #ifdef _HAVE_MATLAB_MODULES_

  • issm/trunk-jpl/src/wrappers/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp

    r23168 r25455  
    3131
    3232        /*Intermediaties*/
    33         int     *index              = NULL;
    34         double  *x_data             = NULL;
    35         double  *y_data             = NULL;
    36         double  *data               = NULL;
     33        int*     index              = NULL;
     34        double*  x_data             = NULL;
     35        double*  y_data             = NULL;
     36        double*  data               = NULL;
    3737        int      nods_data,nels_data;
    3838        int      M_data,N_data;
    39         double  *x_interp           = NULL;
    40         double  *y_interp           = NULL;
     39        double*  x_interp           = NULL;
     40        double*  y_interp           = NULL;
    4141        int      N_interp;
    42         Options *options   = NULL;
    43         double  *data_interp = NULL;
     42        Options* options   = NULL;
     43        double*  data_interp = NULL;
    4444        int      test1,test2,test;
    4545

  • issm/trunk-jpl/src/wrappers/python/Makefile.am

    r24919 r25455  
    3636        ElementConnectivity_python.la \
    3737        ExpToLevelSet_python.la \
     38        InterpFromMesh2d_python.la \
    3839        InterpFromMeshToMesh2d_python.la \
    3940        InterpFromMeshToMesh3d_python.la \
    … …  
    144145ExpToLevelSet_python_la_LIBADD = ${deps} $(PETSCLIB) $(HDF5LIB) $(BLASLAPACKLIB) $(MPILIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB) $(NEOPZLIB)
    145146
     147InterpFromMesh2d_python_la_SOURCES = ../InterpFromMesh2d/InterpFromMesh2d.cpp
     148InterpFromMesh2d_python_la_CXXFLAGS = ${AM_CXXFLAGS}
     149InterpFromMesh2d_python_la_LIBADD = ${deps} $(PETSCLIB) $(HDF5LIB) $(BLASLAPACKLIB) $(MPILIB) $(MULTITHREADINGLIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
     150
    146151InterpFromGridToMesh_python_la_SOURCES = ../InterpFromGridToMesh/InterpFromGridToMesh.cpp
    147152InterpFromGridToMesh_python_la_CXXFLAGS = ${AM_CXXFLAGS}

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

    r25147 r25455  
    7777alreadyloaded=0;
    7878%}}}
    79 for i=sl.basinindx('basin','all'),
    80 
    81         bas=sl.basins{i};
     79for ind=sl.basinindx('basin','all'),
     80
     81        bas=sl.basins{ind};
    8282        disp(sprintf('Meshing basin %s\n',bas.name));
    8383
    … …  
    8585        domain=bas.contour();
    8686
    87         %recover coastline inside basin, using GSHHS_c_L1. It's a lat,long file, hence epsg 4326
    88         coastline=bas.shapefilecrop('shapefile',gshhsshapefile,'epsgshapefile',4326,'threshold',threshold); 
    89 
    90         %mesh: 
     87        %recover coastline inside basin, using GSHHS_c_L1. It's a lat/long file, hence epsg 4326
     88        coastline=bas.shapefilecrop('shapefile',gshhsshapefile,'epsgshapefile',4326,'threshold',threshold);
     89
     90        %mesh:
    9191        md=bamg(model,'domain',domain,'subdomains',coastline,'hmin',hmin,'hmax',hmax,defaultoptions{:});
    9292        %plotmodel(md,'data','mesh');pause(1);
    … …  
    108108%Parameterization:
    109109%Parameterize ice sheets : {{{
    110 
    111110for ind=sl.basinindx('continent',{'antarctica'}),
    112111        disp(sprintf('Parameterizing basin %s\n', sl.icecaps{ind}.miscellaneous.name));
    113112
    114         md=sl.icecaps{ind}; bas=sl.basins{ind};
     113        md=sl.icecaps{ind};
     114        bas=sl.basins{ind};
    115115        %masks :  %{{{
    116116        %ice levelset from domain outlines:
    … …  
    125125        %}}}
    126126        %latlong:  % {{{
    127         [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326'); 
     127        [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326');
    128128        %}}}
    129129        %geometry: {{{
    … …  
    137137        if bas.iscontinentany('antarctica'),
    138138                if testagainst2002,
     139                        % TODO: Check if the following works as expected: 'pos' is empty, so nothing is assigned to 'md.solidearth.surfaceload.icethicknesschange(pos)'
    139140                        md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
    140141                        %antarctica
    … …  
    146147                        md.solidearth.surfaceload.icethicknesschange(pos)=-100*ratio;
    147148                else
    148                         delH=textread('../Data/AIS_delH_trend.txt');
    149                         longAIS=delH(:,1); latAIS=delH(:,2); delHAIS=delH(:,3); index=delaunay(longAIS,latAIS);
    150                         lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360;
     149                        in_fileID=fopen('../Data/AIS_delH_trend.txt', 'r');
     150                        delH=textscan(in_fileID, '%f %f %f');
     151                        fclose(in_fileID);
     152                        longAIS=delH{:,1};
     153                        latAIS=delH{:,2};
     154                        delHAIS=delH{:,3};
     155                        index=delaunay(longAIS,latAIS);
     156                        lat=md.mesh.lat;
     157                        long=md.mesh.long+360;
     158                        pos=find(long>360);
     159                        long(pos)=long(pos)-360;
    151160                        delHAIS=InterpFromMesh2d(index,longAIS,latAIS,delHAIS,long,lat);
    152                         northpole=find_point(md.mesh.long,md.mesh.lat,0,90); delHAIS(northpole)=0;
    153                        
    154                         md.solidearth.surfaceload.icethicknesschange=delHAIS(md.mesh.elements)*[1;1;1]/3/100;
     161                        northpole=find_point(md.mesh.long,md.mesh.lat,0,90);
     162                        delHAIS(northpole)=0;
     163                        md.solidearth.surfaceload.icethicknesschange=mean(delHAIS(md.mesh.elements),2)/100;
    155164                end
    156165
    … …  
    172181end
    173182%}}}
     183
    174184% ParameterizeContinents {{{
    175 
    176 sl.basinindx('continent',{'hemisphereeast','hemispherewest'})
    177 
    178185for ind=sl.basinindx('continent',{'hemisphereeast','hemispherewest'}),
    179186        disp(sprintf('Masks for basin %s\n', sl.icecaps{ind}.miscellaneous.name));
    180         md=sl.icecaps{ind}; bas=sl.basins{ind};
     187        md=sl.icecaps{ind};
     188        bas=sl.basins{ind};
    181189
    182190        %recover lat,long:
    … …  
    184192
    185193        %mask:  %{{{
    186         %Figure out mask from initial mesh: deal with land and ocean masks (land include grounded ice).  %{{{
    187         %first, transform land element  mask into vertex driven one. 
     194        %Figure out mask from initial mesh: deal with land and ocean masks (land
     195        %includes grounded ice).  %{{{
     196        %first, transform land element mask into vertex-driven one
    188197        land=md.private.bamg.landmask;
    189198        land_mask=-ones(md.mesh.numberofvertices,1);
    … …  
    192201        land_mask(md.mesh.elements(landels,:))=1;
    193202
    194         %gothrough edges of each land element
     203        % Gothrough edges of each land element
    195204        connectedels=md.mesh.elementconnectivity(landels,:);
    196205        connectedisonocean=~land(connectedels);
    … …  
    200209        landelsconocean=landels(find(sumconnectedisonocean));
    201210
    202         ind1=[md.mesh.elements(landelsconocean,1); md.mesh.elements(landelsconocean,2); md.mesh.elements(landelsconocean,3)];
    203         ind2=[md.mesh.elements(landelsconocean,2); md.mesh.elements(landelsconocean,3); md.mesh.elements(landelsconocean,1)];
    204 
     211        ind1=[md.mesh.elements(landelsconocean,1);
     212        md.mesh.elements(landelsconocean,2);
     213        md.mesh.elements(landelsconocean,3)];
     214        ind2=[md.mesh.elements(landelsconocean,2);
     215        md.mesh.elements(landelsconocean,3);
     216        md.mesh.elements(landelsconocean,1)];
    205217
    206218        %edge ind1 and ind2:
    … …  
    223235                % {{{
    224236                %greenland
    225                 pos=find(md.mesh.lat > 70 &  md.mesh.lat < 80 & md.mesh.long>-60 & md.mesh.long<-30);
     237                pos=find(md.mesh.lat > 70 & md.mesh.lat < 80 & md.mesh.long>-60 & md.mesh.long<-30);
    226238                md.mask.ice_levelset(pos)=-1;
    227239                % }}}
    228240        end
    229 
    230241        % }}}
    231242        %}}}
    232243        %slr loading/calibration:  {{{
     244        md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
    233245
    234246        if testagainst2002,
    235247                % {{{
    236                 md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
    237248                %greenland
    238249                late=sum(md.mesh.lat(md.mesh.elements),2)/3;
    … …  
    243254
    244255                %correct mask:
    245                 md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; 
     256                md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
    246257                % }}}
    247258        else
    248 
    249                 md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
    250 
    251259                delH=textread('../Data/GIS_delH_trend.txt');
    252                 longGIS=delH(:,1); latGIS=delH(:,2); delHGIS=delH(:,3); index=delaunay(longGIS,latGIS);
    253                 lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360;
     260                longGIS=delH(:,1);
     261                latGIS=delH(:,2);
     262                delHGIS=delH(:,3);
     263                index=delaunay(longGIS,latGIS);
     264                lat=md.mesh.lat;
     265                long=md.mesh.long+360;
     266                pos=find(long>360);
     267                long(pos)=long(pos)-360;
    254268                delHGIS=InterpFromMeshToMesh2d(index,longGIS,latGIS,delHGIS,long,lat);
    255269                delHGISe=delHGIS(md.mesh.elements)*[1;1;1]/3;
    256        
     270
    257271                delH=textread('../Data/GLA_delH_trend_15regions.txt');
    258                 longGLA=delH(:,1); latGLA=delH(:,2); delHGLA=sum(delH(:,3:end),2); index=delaunay(longGLA,latGLA);
    259                 lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360;
     272                longGLA=delH(:,1);
     273                latGLA=delH(:,2);
     274                delHGLA=sum(delH(:,3:end),2);
     275                index=delaunay(longGLA,latGLA);
     276                lat=md.mesh.lat;
     277                long=md.mesh.long+360;
     278                pos=find(long>360);
     279                long(pos)=long(pos)-360;
    260280                delHGLA=InterpFromMeshToMesh2d(index,longGLA,latGLA,delHGLA,long,lat);
    261281                delHGLAe=delHGLA(md.mesh.elements)*[1;1;1]/3;
    … …  
    293313end
    294314% }}}
     315
    295316%%Assemble Earth in 3D {{{
    296317
    … …  
    300321loneedgesdetect=0;
    301322
    302 %create earth model by concatenating all the icecaps in 3d:
     323%create Earth model by concatenating all the icecaps in 3D:
    303324sl.caticecaps('tolerance',tolerance,'loneedgesdetect',loneedgesdetect);
    304325

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

    r24862 r25455  
    88pos=find(rad==min(rad));
    99md.mesh.x(pos)=0.; md.mesh.y(pos)=0.; %the closest node to the center is changed to be exactly at the center
    10 xelem=md.mesh.x(md.mesh.elements)*[1;1;1]/3.;
    11 yelem=md.mesh.y(md.mesh.elements)*[1;1;1]/3.;
     10xelem=mean(md.mesh.x(md.mesh.elements),2);
     11yelem=mean(md.mesh.y(md.mesh.elements),2);
    1212rad=sqrt(xelem.^2+yelem.^2);
    1313flags=zeros(md.mesh.numberofelements,1);

  • issm/trunk-jpl/test/NightlyRun/test423.py

    r24862 r25455  
    1919md.mesh.x[pos] = 0.
    2020md.mesh.y[pos] = 0.  #the closest node to the center is changed to be exactly at the center
    21 xelem = np.mean(md.mesh.x[md.mesh.elements.astype(int) - 1], axis=1)
    22 yelem = np.mean(md.mesh.y[md.mesh.elements.astype(int) - 1], axis=1)
     21xelem = np.mean(md.mesh.x[md.mesh.elements - 1], axis=1)
     22yelem = np.mean(md.mesh.y[md.mesh.elements - 1], axis=1)
    2323rad = np.sqrt(xelem**2 + yelem**2)
    2424flags = np.zeros(md.mesh.numberofelements)

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

    r24862 r25455  
    88pos=find(rad==min(rad));
    99md.mesh.x(pos)=0.; md.mesh.y(pos)=0.; %the closest node to the center is changed to be exactly at the center
    10 xelem=md.mesh.x(md.mesh.elements)*[1;1;1]/3.;
    11 yelem=md.mesh.y(md.mesh.elements)*[1;1;1]/3.;
     10xelem=mean(md.mesh.x(md.mesh.elements), 2);
     11yelem=mean(md.mesh.y(md.mesh.elements), 2);
    1212rad=sqrt(xelem.^2+yelem.^2);
    1313flags=zeros(md.mesh.numberofelements,1);

  • issm/trunk-jpl/test/NightlyRun/test433.py

    r24862 r25455  
    1919md.mesh.x[pos] = 0.
    2020md.mesh.y[pos] = 0.  #the closest node to the center is changed to be exactly at the center
    21 xelem = np.mean(md.mesh.x[md.mesh.elements.astype(int) - 1], axis=1)
    22 yelem = np.mean(md.mesh.y[md.mesh.elements.astype(int) - 1], axis=1)
     21xelem = np.mean(md.mesh.x[md.mesh.elements - 1], axis=1)
     22yelem = np.mean(md.mesh.y[md.mesh.elements - 1], axis=1)
    2323rad = np.sqrt(xelem**2 + yelem**2)
    2424flags = np.zeros(md.mesh.numberofelements)

  • issm/trunk-jpl/test/Par/79North.py

    r24862 r25455  
    1919thickness = numpy.array(archread('../Data/79North.arch', 'thickness'))
    2020
    21 md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)[0][:, 0]
    22 md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)[0][:, 0]
    23 md.geometry.surface = InterpFromMeshToMesh2d(index, x, y, surface, md.mesh.x, md.mesh.y)[0][:, 0]
    24 md.geometry.thickness = InterpFromMeshToMesh2d(index, x, y, thickness, md.mesh.x, md.mesh.y)[0][:, 0]
     21md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)[:, 0]
     22md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)[:, 0]
     23md.geometry.surface = InterpFromMeshToMesh2d(index, x, y, surface, md.mesh.x, md.mesh.y)[:, 0]
     24md.geometry.thickness = InterpFromMeshToMesh2d(index, x, y, thickness, md.mesh.x, md.mesh.y)[:, 0]
    2525md.geometry.base = md.geometry.surface - md.geometry.thickness
    2626

  • issm/trunk-jpl/test/Par/Pig.py

    r24862 r25455  
    2020bed = np.array(archread('../Data/Pig.arch', 'bed'))
    2121
    22 md.inversion.vx_obs = InterpFromMeshToMesh2d(index, x, y, vx_obs, md.mesh.x, md.mesh.y)[0][:, 0]
    23 md.inversion.vy_obs = InterpFromMeshToMesh2d(index, x, y, vy_obs, md.mesh.x, md.mesh.y)[0][:, 0]
    24 md.geometry.surface = InterpFromMeshToMesh2d(index, x, y, surface, md.mesh.x, md.mesh.y)[0][:, 0]
    25 md.geometry.thickness = InterpFromMeshToMesh2d(index, x, y, thickness, md.mesh.x, md.mesh.y)[0][:, 0]
     22md.inversion.vx_obs = InterpFromMeshToMesh2d(index, x, y, vx_obs, md.mesh.x, md.mesh.y)[:, 0]
     23md.inversion.vy_obs = InterpFromMeshToMesh2d(index, x, y, vy_obs, md.mesh.x, md.mesh.y)[:, 0]
     24md.geometry.surface = InterpFromMeshToMesh2d(index, x, y, surface, md.mesh.x, md.mesh.y)[:, 0]
     25md.geometry.thickness = InterpFromMeshToMesh2d(index, x, y, thickness, md.mesh.x, md.mesh.y)[:, 0]
    2626md.geometry.base = md.geometry.surface - md.geometry.thickness
    2727md.geometry.bed = np.array(md.geometry.base)
    2828pos = np.where(md.mask.ocean_levelset < 0.)
    29 md.geometry.bed[pos] = InterpFromMeshToMesh2d(index, x, y, bed, md.mesh.x[pos], md.mesh.y[pos])[0][:, 0]
     29md.geometry.bed[pos] = InterpFromMeshToMesh2d(index, x, y, bed, md.mesh.x[pos], md.mesh.y[pos])[:, 0]
    3030md.initialization.vx = md.inversion.vx_obs
    3131md.initialization.vy = md.inversion.vy_obs

  • issm/trunk-jpl/test/Par/SquareSheetConstrained.py

    r25129 r25455  
    2828index = archread('../Data/SquareSheetConstrained.arch', 'index').astype(int)
    2929
    30 md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)[0]
    31 md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)[0]
     30md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)
     31md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)
    3232md.initialization.vz = np.zeros((md.mesh.numberofvertices))
    3333md.initialization.pressure = np.zeros((md.mesh.numberofvertices))

  • issm/trunk-jpl/test/Par/SquareSheetConstrainedCO2.py

    r24862 r25455  
    4444index = archread('../Data/SquareSheetConstrained.arch', 'index').astype(int)
    4545
    46 [md.initialization.vx] = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)
    47 [md.initialization.vy] = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)
     46md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)
     47md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)
    4848md.initialization.vz = np.zeros((md.mesh.numberofvertices))
    4949md.initialization.pressure = np.zeros((md.mesh.numberofvertices))

  • issm/trunk-jpl/test/Par/SquareSheetShelf.py

    r24862 r25455  
    3131index = np.array(archread('../Data/SquareSheetShelf.arch', 'index')).astype(int)
    3232
    33 [md.initialization.vx] = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)
    34 [md.initialization.vy] = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)
     33md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)
     34md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)
    3535md.initialization.vz = np.zeros((md.mesh.numberofvertices))
    3636md.initialization.pressure = np.zeros((md.mesh.numberofvertices))

  • issm/trunk-jpl/test/Par/SquareShelf2.py

    r24862 r25455  
    3232#dbg - end
    3333
    34 [md.initialization.vx] = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)
    35 [md.initialization.vy] = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)
     34md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)
     35md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)
    3636md.initialization.vz = np.zeros((md.mesh.numberofvertices))
    3737md.initialization.pressure = np.zeros((md.mesh.numberofvertices))

  • issm/trunk-jpl/test/Par/SquareShelfConstrained.py

    r24862 r25455  
    2828index = np.array(archread('../Data/SquareShelfConstrained.arch', 'index').astype(int))
    2929
    30 [md.initialization.vx] = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)
    31 [md.initialization.vy] = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)
     30md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)
     31md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)
    3232md.initialization.vz = np.zeros((md.mesh.numberofvertices))
    3333md.initialization.pressure = np.zeros((md.mesh.numberofvertices))
Note: See TracChangeset for help on using the changeset viewer.