Index: /issm/trunk-jpl/src/m/classes/bamgmesh.py =================================================================== --- /issm/trunk-jpl/src/m/classes/bamgmesh.py (revision 25454) +++ /issm/trunk-jpl/src/m/classes/bamgmesh.py (revision 25455) @@ -3,9 +3,8 @@ class bamgmesh(object): - """ - BAMGMESH class definition + """BAMGMESH class definition - Usage: - bamgmesh(varargin) + Usage: + bamgmesh(varargin) """ Index: /issm/trunk-jpl/src/m/classes/basin.m =================================================================== --- /issm/trunk-jpl/src/m/classes/basin.m (revision 25454) +++ /issm/trunk-jpl/src/m/classes/basin.m (revision 25455) @@ -117,5 +117,6 @@ %recover options options=pairoptions(varargin{:}); - x=[]; y=[]; + x=[]; + y=[]; %go through boundaries, recover edges and project them in the basin epsg reference frame: @@ -127,4 +128,5 @@ y=[y;contour.y]; end + %close the contour: if x(end)~=x(1) | y(end)~=y(1), @@ -181,5 +183,4 @@ error('Basin shapefilecrop error message: epsgshapefile or projshapefile should be specified'); end - %create list of contours that have critical length > threshold (in lat,long): @@ -209,10 +210,12 @@ for i=1:length(contours), h=contours(i); - isin=inpolygon(h.x,h.y,ba.x,ba.y); + isin=inpolygon(h.x,h.y,ba.x,ba.y); if ~isempty(find(isin==0)), flags(i)=1; end end - pos=find(flags); contours(pos)=[]; + + pos=find(flags); + contours(pos)=[]; %Two options: Index: /issm/trunk-jpl/src/m/classes/basin.py =================================================================== --- /issm/trunk-jpl/src/m/classes/basin.py (revision 25454) +++ /issm/trunk-jpl/src/m/classes/basin.py (revision 25455) @@ -1,9 +1,15 @@ +from collections import OrderedDict import math + +from bamg import bamg +import matplotlib.path as path import numpy as np from boundary import boundary from epsg2proj import epsg2proj +from fielddisplay import fielddisplay from gdaltransform import gdaltransform from helpers import fileparts +from inpolygon import inpolygon from laea import laea from pairoptions import pairoptions @@ -13,10 +19,9 @@ class basin(object): #{{{ - ''' - BASIN class definition + """BASIN class definition Usage: basin = basin() - ''' + """ def __init__(self, *args): #{{{ @@ -136,15 +141,20 @@ for i in range(len(self.boundaries)): boundary = self.boundaries[i] - contour = boundary.edges() - contour.x, contour.y = gdaltransform(contour.x, contour.y, boundary.proj, self.proj) - x = [[x], [contour.x]] - y = [[y], [contour.y]] + contours = boundary.edges() # NOTE: Differs from MATLAB in that shpread returns a list + for j in range(len(contours)): + contours[j]['x'], contours[j]['y'] = gdaltransform(contours[j]['x'], contours[j]['y'], boundary.proj, self.proj) + if x == []: + x = np.array(contours[j]['x']) + y = np.array(contours[j]['y']) + else: + x = np.concatenate((x, contours[j]['x']), axis=0) + y = np.concatenate((y, contours[j]['y']), axis=0) #close the contour if x[-1] != x[0] or y[-1] != y[0]: - x.append(x[0]) - y.append(y[0]) - - out = {} + x = np.append(x, x[0]) + y = np.append(y, y[0]) + + out = OrderedDict() out['x'] = x out['y'] = y @@ -182,5 +192,5 @@ threshold = options.getfieldvalue('threshold', .65) #.65 degrees lat, long inshapefile = options.getfieldvalue('shapefile') - outshapefile = options.getfieldvalue('outshapefile', '') + outputshapefile = options.getfieldvalue('outputshapefile', '') if options.exist('epsgshapefile') and options.exist('projshapefile'): @@ -199,16 +209,21 @@ for i in range(len(contours)): contour = contours[i] - carea = polyarea(contour.x, contour.y) + carea = polyarea(contour['x'], contour['y']) clength = math.sqrt(carea) if clength < threshold: - llist = np.concatenate(llist, i) - setattr(contours, llist, []) + llist.append(i) + + contours_above_threshold = [] + for i in range(len(contours)): + if i not in llist: + contours_above_threshold.append(contours[i]) + contours = contours_above_threshold #project onto reference frame for i in range(len(contours)): h = contours[i] - h.x, h.y = gdaltransform(h.x, h.y, projshapefile, self.proj) - contours[i].x = h.x - contours[i].y = h.y + h['x'], h['y'] = gdaltransform(h['x'], h['y'], projshapefile, self.proj) + contours[i]['x'] = h['x'] + contours[i]['y'] = h['y'] #only keep the contours that are inside the basin (take centroids) @@ -217,15 +232,24 @@ for i in range(len(contours)): h = contours[i] - isin = inpolygon(h.x, h.y, ba.x, ba.y) - if len(np.where(isin == 0, isin, isin)): + isin = inpolygon(h['x'], h['y'], ba['x'], ba['y']) + if len(np.where(isin == 0)[0]): flags[i] = 1 - pos = flags.nonzero() - contours[pos] = [] + + pos = flags.nonzero()[0] + + contours_in_basin = [] + for i in range(len(contours)): + if i not in pos: + contours_in_basin.append(contours[i]) + contours = contours_in_basin #Two options + output = None if outputshapefile == '': output = contours else: shpwrite(contours, outputshapefile) + + return output #}}} #}}} Index: /issm/trunk-jpl/src/m/classes/boundary.py =================================================================== --- /issm/trunk-jpl/src/m/classes/boundary.py (revision 25454) +++ /issm/trunk-jpl/src/m/classes/boundary.py (revision 25455) @@ -4,4 +4,5 @@ from epsg2proj import epsg2proj +from fielddisplay import fielddisplay from helpers import * from pairoptions import pairoptions @@ -10,10 +11,9 @@ class boundary(object): #{{{ - ''' - BOUNDARY class definition + """BOUNDARY class definition Usage: boundary = boundary() - ''' + """ def __init__(self, *args): #{{{ @@ -32,4 +32,8 @@ self.shpfilename = options.getfieldvalue('shpfilename', '') self.orientation = options.getfieldvalue('orientation', 'normal') + + # If shppath is missing trailing slash, add it + if self.shppath[-1] != '/': + self.shppath += '/' #figure out projection string: @@ -77,13 +81,15 @@ #read domain path, name, ext = fileparts(self.shpfilename) - if ext != 'shp': - ext = 'shp' - output = shpread('{}/{}.{}'.format(self.shppath, name, ext)) + if ext != '.shp': + ext = '.shp' + output = shpread('{}{}{}'.format(self.shppath, name, ext)) #do we reverse? if self.orientation == 'reverse': for i in range(len(output)): - output[i].x = np.flipud(output[i].x) - output[i].y = np.flipud(output[i].y) + output[i]['x'] = np.flipud(output[i]['x']) + output[i]['y'] = np.flipud(output[i]['y']) + + return output #}}} @@ -107,12 +113,12 @@ #read domain path, name, ext = fileparts(self.shpfilename) - if ext != 'shp': - ext = 'shp' - domain = shpread('{}/{}.{}'.format(self.shppath, name, ext)) + if ext != '.shp': + ext = '.shp' + domain = shpread('{}{}{}'.format(self.shppath, name, ext)) #convert boundary to another reference frame #{{{ for i in range(len(domain)): try: - x, y = gdaltransform(domain[i].x, domain[i].y, self.proj, proj) + x, y = gdaltransform(domain[i]['x'], domain[i]['y'], self.proj, proj) except error as e: print(e) @@ -170,7 +176,7 @@ #read domain path, name, ext = fileparts(self.shpfilename) - if ext != 'shp': - ext = 'shp' - domain = shpread('{}/{}.{}'.format(self.shppath, name, ext)) + if ext != '.shp': + ext = '.shp' + domain = shpread('{}{}{}'.format(self.shppath, name, ext)) #TODO: Finish translating from MATLAB after test2004.py runs without plot Index: /issm/trunk-jpl/src/m/classes/giacaron.py =================================================================== --- /issm/trunk-jpl/src/m/classes/giacaron.py (revision 25455) +++ /issm/trunk-jpl/src/m/classes/giacaron.py (revision 25455) @@ -0,0 +1,53 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from project3d import project3d +from WriteData import WriteData + + +class giacaron(object): + """ + GIA class definition for Caron model (Caron et al, Geophysical Journal + International, 2017) + + Usage: + giacaron = giacaron() + """ + + def __init__(self, *args): #{{{ + #Physical constants + self.gravitational_constant = np.nan + self.surface_radius = np.nan + self.core_mantle_boudary_radius = np.nan + self.inner_core_boudary_radius = np.nan + self.stress_norm = np.nan + self.gravity_norm = np.nan + self.radius_norm = np.nan + + #Numerical parameters + self.allow_layer_deletion = np.nan + self.verbose_mode = np.nan + + #GIA problem setup + self.forcing_type = np.nan + self.isincompressible = np.nan + self.benchmark_mode = np.nan + self.calculate_sea_level = np.nan + self.calculate_rotational_feedback = np.nan + self.subtract_present_day = np.nan + self.ntime = np.nan + self.nphi = np.nan + + #Earth model + + + + + nargin = len(args) + + if nargin == 0: + self.setdefaultparameters() + else: + raise Exception('constructor not supported') + #}}} Index: /issm/trunk-jpl/src/m/classes/mesh3dsurface.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh3dsurface.py (revision 25454) +++ /issm/trunk-jpl/src/m/classes/mesh3dsurface.py (revision 25455) @@ -9,12 +9,11 @@ class mesh3dsurface(object): - ''' - MESH3DSURFACE class definition - - Usage: - mesh3dsurface = mesh3dsurface() - ''' - - def __init__(self, *args): # {{{ + """MESH3DSURFACE class definition + + Usage: + mesh3dsurface = mesh3dsurface() + """ + + def __init__(self, *args): #{{{ self.x = np.nan self.y = np.nan @@ -40,7 +39,8 @@ self.extractedelements = np.nan - if not len(args): + nargs = len(args) + if not nargs: self.setdefaultparameters() - elif len(args) == 1: + elif nargs == 1: self = mesh3dsurface() arg = args[1] @@ -52,40 +52,39 @@ else: raise RuntimeError('constructor not supported') - - # }}} - - def __repr__(self): # {{{ - string = ' 2D tria Mesh (horizontal):' - - string += '\n Elements and vertices:' - string = "%s\n%s" % (string, fielddisplay(self, 'numberofelements', 'number of elements')) - string = "%s\n%s" % (string, fielddisplay(self, 'numberofvertices', 'number of vertices')) - string = "%s\n%s" % (string, fielddisplay(self, 'elements', 'vertex indices of the mesh elements')) - string = "%s\n%s" % (string, fielddisplay(self, 'x', 'vertices x coordinate [m]')) - string = "%s\n%s" % (string, fielddisplay(self, 'y', 'vertices y coordinate [m]')) - string = "%s\n%s" % (string, fielddisplay(self, 'z', 'vertices z coordinate [m]')) - string = "%s\n%s" % (string, fielddisplay(self, 'lat', 'vertices latitude [degrees]')) - string = "%s\n%s" % (string, fielddisplay(self, 'long', 'vertices longitude [degrees]')) - string = "%s\n%s" % (string, fielddisplay(self, 'r', 'vertices radius [m]')) - - string = "%s\n%s" % (string, fielddisplay(self, 'edges', 'edges of the 2d mesh (vertex1 vertex2 element1 element2)')) - string = "%s\n%s" % (string, fielddisplay(self, 'numberofedges', 'number of edges of the 2d mesh')) - - string += '\n Properties:' - string = "%s\n%s" % (string, fielddisplay(self, 'vertexonboundary', 'vertices on the boundary of the domain flag list')) - string = "%s\n%s" % (string, fielddisplay(self, 'segments', 'edges on domain boundary (vertex1 vertex2 element)')) - string = "%s\n%s" % (string, fielddisplay(self, 'segmentmarkers', 'number associated to each segment')) - string = "%s\n%s" % (string, fielddisplay(self, 'vertexconnectivity', 'list of elements connected to vertex_i')) - string = "%s\n%s" % (string, fielddisplay(self, 'elementconnectivity', 'list of elements adjacent to element_i')) - string = "%s\n%s" % (string, fielddisplay(self, 'average_vertex_connectivity', 'average number of vertices connected to one vertex')) - - string += '\n Extracted model():' - string = "%s\n%s" % (string, fielddisplay(self, 'extractedvertices', 'vertices extracted from the model()')) - string = "%s\n%s" % (string, fielddisplay(self, 'extractedelements', 'elements extracted from the model()')) - - return string - # }}} - - def setdefaultparameters(self): # {{{ + #}}} + + def __repr__(self): #{{{ + s = ' 2D tria Mesh (surface):' + + s += '\n Elements and vertices:' + s = "%s\n%s" % (s, fielddisplay(self, 'numberofelements', 'number of elements')) + s = "%s\n%s" % (s, fielddisplay(self, 'numberofvertices', 'number of vertices')) + s = "%s\n%s" % (s, fielddisplay(self, 'elements', 'vertex indices of the mesh elements')) + s = "%s\n%s" % (s, fielddisplay(self, 'x', 'vertices x coordinate [m]')) + s = "%s\n%s" % (s, fielddisplay(self, 'y', 'vertices y coordinate [m]')) + s = "%s\n%s" % (s, fielddisplay(self, 'z', 'vertices z coordinate [m]')) + s = "%s\n%s" % (s, fielddisplay(self, 'lat', 'vertices latitude [degrees]')) + s = "%s\n%s" % (s, fielddisplay(self, 'long', 'vertices longitude [degrees]')) + s = "%s\n%s" % (s, fielddisplay(self, 'r', 'vertices radius [m]')) + + s = "%s\n%s" % (s, fielddisplay(self, 'edges', 'edges of the 2d mesh (vertex1 vertex2 element1 element2)')) + s = "%s\n%s" % (s, fielddisplay(self, 'numberofedges', 'number of edges of the 2d mesh')) + + s += '\n Properties:' + s = "%s\n%s" % (s, fielddisplay(self, 'vertexonboundary', 'vertices on the boundary of the domain flag list')) + s = "%s\n%s" % (s, fielddisplay(self, 'segments', 'edges on domain boundary (vertex1 vertex2 element)')) + s = "%s\n%s" % (s, fielddisplay(self, 'segmentmarkers', 'number associated to each segment')) + s = "%s\n%s" % (s, fielddisplay(self, 'vertexconnectivity', 'list of elements connected to vertex_i')) + s = "%s\n%s" % (s, fielddisplay(self, 'elementconnectivity', 'list of elements adjacent to element_i')) + s = "%s\n%s" % (s, fielddisplay(self, 'average_vertex_connectivity', 'average number of vertices connected to one vertex')) + + s += '\n Extracted model():' + s = "%s\n%s" % (s, fielddisplay(self, 'extractedvertices', 'vertices extracted from the model()')) + s = "%s\n%s" % (s, fielddisplay(self, 'extractedelements', 'elements extracted from the model()')) + + return s + #}}} + + def setdefaultparameters(self): #{{{ #The connectivity is the average number of nodes linked to a given node #through an edge. This connectivity is used to initially allocate @@ -94,7 +93,7 @@ #test/NightlyRun/runme.py. self.average_vertex_connectivity = 25 - # }}} - - def checkconsistency(self, md, solution, analyses): # {{{ + #}}} + + def checkconsistency(self, md, solution, analyses): #{{{ md = checkfield(md, 'fieldname', 'mesh.x', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) md = checkfield(md, 'fieldname', 'mesh.y', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) @@ -116,7 +115,7 @@ return md - # }}} - - def marshall(self, prefix, md, fid): # {{{ + #}}} + + def marshall(self, prefix, md, fid): #{{{ WriteData(fid, prefix, 'name', 'md.mesh.domain_type', 'data', 'Domain' + self.domaintype(), 'format', 'String') WriteData(fid, prefix, 'name', 'md.mesh.domain_dimension', 'data', self.dimension(), 'format', 'Integer') @@ -134,19 +133,19 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'average_vertex_connectivity', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'vertexonboundary', 'format', 'DoubleMat', 'mattype', 1) - # }}} - - def domaintype(self): # {{{ + #}}} + + def domaintype(self): #{{{ return '3Dsurface' - # }}} - - def dimension(self): # {{{ + #}}} + + def dimension(self): #{{{ return 2 - # }}} - - def elementtype(self): # {{{ + #}}} + + def elementtype(self): #{{{ return 'Tria' - # }}} - - def processmesh(self, options): # {{{ + #}}} + + def processmesh(self, options): #{{{ isplanet = 1 is2d = 0 @@ -156,8 +155,9 @@ y = self.y z = self.z - return [x, y, z, elements, is2d, isplanet] - # }}} - - def savemodeljs(self, fid, modelname): # {{{ + + return x, y, z, elements, is2d, isplanet + #}}} + + def savemodeljs(self, fid, modelname): #{{{ fid.write(' #s.mesh = new mesh3dsurface()\n' % modelname) writejs1Darray(fid, [modelname, '.mesh.x'], self.x) @@ -180,7 +180,7 @@ writejs1Darray(fid, [modelname, '.mesh.extractedvertices'], self.extractedvertices) writejs1Darray(fid, [modelname, '.mesh.extractedelements'], self.extractedelements) - # }}} - - def export(self, *args): # {{{ + #}}} + + def export(self, *args): #{{{ options = pairoptions(*args) @@ -269,3 +269,3 @@ #write style file: applyqgisstyle(filename, 'mesh') - # }}} + #}}} Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 25454) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 25455) @@ -79,10 +79,5 @@ class model(object): #properties - def __init__(self): #{{{ - ''' classtype = model.properties - for classe in dict.keys(classtype): - print classe - self.__dict__[classe] = classtype[str(classe)] - ''' + def __init__(self): #{{{ self.mesh = mesh2d() self.mask = mask() @@ -232,33 +227,32 @@ def extract(self, area): # {{{ + """EXTRACT - extract a model according to an Argus contour or flag list + + This routine extracts a submodel from a bigger model with respect to a given contour + md must be followed by the corresponding exp file or flags list + It can either be a domain file (argus type, .exp extension), or an array of element flags. + If user wants every element outside the domain to be + extract2d, add '~' to the name of the domain file (ex: '~HO.exp') + an empty string '' will be considered as an empty domain + a string 'all' will be considered as the entire domain + + Usage: + md2 = extract(md, area) + + Examples: + md2 = extract(md, 'Domain.exp') + + See also: EXTRUDE, COLLAPSE """ - extract - extract a model according to an Argus contour or flag list - - This routine extracts a submodel from a bigger model with respect to a given contour - md must be followed by the corresponding exp file or flags list - It can either be a domain file (argus type, .exp extension), or an array of element flags. - If user wants every element outside the domain to be - extract2d, add '~' to the name of the domain file (ex: '~HO.exp') - an empty string '' will be considered as an empty domain - a string 'all' will be considered as the entire domain - - Usage: - md2 = extract(md, area) - - Examples: - md2 = extract(md, 'Domain.exp') - - See also: EXTRUDE, COLLAPSE - """ - - #copy model + + #copy model md1 = copy.deepcopy(self) - #get elements that are inside area + #get elements that are inside area flag_elem = FlagElements(md1, area) if not np.any(flag_elem): raise RuntimeError("extracted model is empty") - #kick out all elements with 3 dirichlets + #kick out all elements with 3 dirichlets spc_elem = np.nonzero(np.logical_not(flag_elem))[0] spc_node = np.unique(md1.mesh.elements[spc_elem, :]) - 1 @@ -268,9 +262,9 @@ flag_elem[pos] = 0 - #extracted elements and nodes lists + #extracted elements and nodes lists pos_elem = np.nonzero(flag_elem)[0] pos_node = np.unique(md1.mesh.elements[pos_elem, :]) - 1 - #keep track of some fields + #keep track of some fields numberofvertices1 = md1.mesh.numberofvertices numberofelements1 = md1.mesh.numberofelements @@ -280,5 +274,5 @@ flag_node[pos_node] = 1 - #Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) + #Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) Pelem = np.zeros(numberofelements1, int) Pelem[pos_elem] = np.arange(1, numberofelements2 + 1) @@ -286,5 +280,5 @@ Pnode[pos_node] = np.arange(1, numberofvertices2 + 1) - #renumber the elements (some node won't exist anymore) + #renumber the elements (some node won't exist anymore) elements_1 = copy.deepcopy(md1.mesh.elements) elements_2 = elements_1[pos_elem, :] @@ -297,12 +291,12 @@ elements_2[:, 5] = Pnode[elements_2[:, 5] - 1] - #OK, now create the new model! - - #take every field from model + #OK, now create the new model! + + #take every field from model md2 = copy.deepcopy(md1) - #automatically modify fields - - #loop over model fields + #automatically modify fields + + #loop over model fields model_fields = vars(md1) for fieldi in model_fields: @@ -310,5 +304,5 @@ field = getattr(md1, fieldi) fieldsize = np.shape(field) - if hasattr(field, '__dict__') and fieldi not in ['results']: #recursive call + if hasattr(field, '__dict__') and fieldi not in ['results']: #recursive call object_fields = vars(field) for fieldj in object_fields: @@ -416,6 +410,6 @@ #recreate segments if md1.mesh.__class__.__name__ == 'mesh2d': - md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements, md2.mesh.numberofvertices)[0] - md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements, md2.mesh.vertexconnectivity)[0] + md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements, md2.mesh.numberofvertices) + md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements, md2.mesh.vertexconnectivity) md2.mesh.segments = contourenvelope(md2) md2.mesh.vertexonboundary = np.zeros(numberofvertices2, bool) @@ -423,6 +417,6 @@ else: #First do the connectivity for the contourenvelope in 2d - md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements2d, md2.mesh.numberofvertices2d)[0] - md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements2d, md2.mesh.vertexconnectivity)[0] + md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements2d, md2.mesh.numberofvertices2d) + md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements2d, md2.mesh.vertexconnectivity) segments = contourenvelope(md2) md2.mesh.vertexonboundary = np.zeros(int(numberofvertices2 / md2.mesh.numberoflayers), bool) @@ -430,6 +424,6 @@ md2.mesh.vertexonboundary = np.tile(md2.mesh.vertexonboundary, md2.mesh.numberoflayers) #Then do it for 3d as usual - md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements, md2.mesh.numberofvertices)[0] - md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements, md2.mesh.vertexconnectivity)[0] + md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements, md2.mesh.numberofvertices) + md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements, md2.mesh.vertexconnectivity) #Boundary conditions: Dirichlets on new boundary @@ -504,5 +498,5 @@ setattr(fieldr, solutionsubfield, subfield) - #Keep track of pos_node and pos_elem + #Keep track of pos_node and pos_elem md2.mesh.extractedvertices = pos_node + 1 md2.mesh.extractedelements = pos_elem + 1 @@ -512,9 +506,8 @@ def extrude(md, *args): # {{{ - """ - EXTRUDE - vertically extrude a 2d mesh - - vertically extrude a 2d mesh and create corresponding 3d mesh. - The vertical distribution can: + """EXTRUDE - vertically extrude a 2d mesh + + vertically extrude a 2d mesh and create corresponding 3d mesh. + The vertical distribution can: - follow a polynomial law - follow two polynomial laws, one for the lower part and one for the upper part of the mesh @@ -522,16 +515,17 @@ - Usage: - md = extrude(md, numlayers, extrusionexponent) - md = extrude(md, numlayers, lowerexponent, upperexponent) - md = extrude(md, listofcoefficients) - - Example: - md = extrude(md, 15, 1.3) - md = extrude(md, 15, 1.3, 1.2) - md = extrude(md, [0 0.2 0.5 0.7 0.9 0.95 1]) - - See also: MODELEXTRACT, COLLAPSE + Usage: + md = extrude(md, numlayers, extrusionexponent) + md = extrude(md, numlayers, lowerexponent, upperexponent) + md = extrude(md, listofcoefficients) + + Example: + md = extrude(md, 15, 1.3) + md = extrude(md, 15, 1.3, 1.2) + md = extrude(md, [0 0.2 0.5 0.7 0.9 0.95 1]) + + See also: MODELEXTRACT, COLLAPSE """ + #some checks on list of arguments if len(args) > 3 or len(args) < 1: @@ -698,7 +692,6 @@ # }}} - def collapse(md): #{{{ - ''' - collapses a 3d mesh into a 2d mesh + def collapse(md): #{{{ + """COLLAPSE - collapses a 3d mesh into a 2d mesh This routine collapses a 3d model into a 2d model and collapses all @@ -706,6 +699,6 @@ Usage: - md = collapse(md) - ''' + md = collapse(md) + """ #Check that the model is really a 3d model @@ -881,9 +874,8 @@ md.mesh.scale_factor = project2d(md, md.mesh.scale_factor, 1) md.mesh = mesh - md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)[0] - md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)[0] + md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices) + md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity) md.mesh.segments = contourenvelope(md) return md - #}}} Index: /issm/trunk-jpl/src/m/classes/pairoptions.py =================================================================== --- /issm/trunk-jpl/src/m/classes/pairoptions.py (revision 25454) +++ /issm/trunk-jpl/src/m/classes/pairoptions.py (revision 25455) @@ -1,11 +1,11 @@ from collections import OrderedDict + class pairoptions(object): - """ - PAIROPTIONS class definition + """PAIROPTIONS class definition - Usage: - pairoptions = pairoptions() - pairoptions = pairoptions('module', true, 'solver', false) + Usage: + pairoptions = pairoptions() + pairoptions = pairoptions('module', true, 'solver', false) """ @@ -39,5 +39,6 @@ def buildlist(self, *arg): # {{{ - """BUILDLIST - build list of objects from input""" + """BUILDLIST - build list of objects from input + """ #check length of input Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.m =================================================================== --- /issm/trunk-jpl/src/m/classes/sealevelmodel.m (revision 25454) +++ /issm/trunk-jpl/src/m/classes/sealevelmodel.m (revision 25455) @@ -21,5 +21,5 @@ mergedcaps = 0; transitions = {}; - eltransitions = {}; + eltransitions = {}; planet = ''; %}}} @@ -170,38 +170,38 @@ function intersections(self,varargin) % {{{ - options=pairoptions(varargin{:}); - force=getfieldvalue(options,'force',0); + options=pairoptions(varargin{:}); + force=getfieldvalue(options,'force',0); + + %initialize, to avoid issues of having more transitions than meshes. + self.transitions={}; self.eltransitions={}; + + %for elements: + xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3; + ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3; + ze=self.earth.mesh.z(self.earth.mesh.elements)*[1;1;1]/3; + + for i=1:length(self.icecaps), + mdi=self.icecaps{i}; + mdi=TwoDToThreeD(mdi,self.planet); - %initialize, to avoid issues of having more transitions than meshes. - self.transitions={}; self.eltransitions={}; - - %for elements: - xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3; - ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3; - ze=self.earth.mesh.z(self.earth.mesh.elements)*[1;1;1]/3; - - for i=1:length(self.icecaps), - mdi=self.icecaps{i}; - mdi=TwoDToThreeD(mdi,self.planet); - - %for elements: - xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3; - yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3; - zei=mdi.mesh.z(mdi.mesh.elements)*[1;1;1]/3; - - disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name)); + %for elements: + xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3; + yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3; + zei=mdi.mesh.z(mdi.mesh.elements)*[1;1;1]/3; - 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); - - self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force); - end + disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name)); + + 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); + + self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force); + end end % }}} function checkintersections(self) % {{{ - flags=zeros(self.earth.mesh.numberofvertices,1); - for i=1:length(self.basins), - flags(self.transitions{i})=i; - end - plotmodel(self.earth,'data',flags,'coastline','on'); + flags=zeros(self.earth.mesh.numberofvertices,1); + for i=1:length(self.basins), + flags(self.transitions{i})=i; + end + plotmodel(self.earth,'data',flags,'coastline','on'); end % }}} @@ -227,7 +227,9 @@ end continent=unique(continent); + else + %nothing to do: assume we have a list of continents end else - %nothing to do, we have a list of continents + %nothing to do: assume we have a list of continents end else @@ -416,5 +418,5 @@ eval(['self.earth.' string '=field;']); %do not forget to plug the field variable into its final location %}}} - elseif n==self.icecaps{1} .mesh.numberofelements, + elseif n==self.icecaps{1}.mesh.numberofelements, eval(['self.earth.' string '=zeros(self.earth.mesh.numberofelements,1);']); for i=1:length(self.icecaps), Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.py =================================================================== --- /issm/trunk-jpl/src/m/classes/sealevelmodel.py (revision 25454) +++ /issm/trunk-jpl/src/m/classes/sealevelmodel.py (revision 25455) @@ -7,4 +7,6 @@ from meshintersect3d import meshintersect3d from miscellaneous import miscellaneous +from model import model +from modelmerge3d import modelmerge3d from pairoptions import pairoptions from private import private @@ -13,6 +15,5 @@ class sealevelmodel(object): - ''' - SEALEVELMODEL class definition + """SEALEVELMODEL class definition Usage: @@ -27,5 +28,5 @@ 'earth', md_earth ) - ''' + """ def __init__(self, *args): #{{{ @@ -207,18 +208,20 @@ bas = options.getfieldvalue('basin', 'all') - #expand continent list #{{{ + # Expand continent list #{{{ if type(continent) == np.ndarray: if continent.shape[1] == 1: if continent[0] == 'all': - #need to transform this into a list of continents + # Need to transform this into a list of continents continent = [] for i in range(len(self.basins)): continent.append(self.basins[i].continent) continent = np.unique(continent) + else: + pass # Nothing to do: assume we have a list of continents else: - pass #nothing to do, we have a list of continents + pass # Nothing to do: assume we have a list of continents else: if continent == 'all': - #need to transform this into a list of continents + # Need to transform this into a list of continents continent = [] for i in range(len(self.basins)): @@ -226,11 +229,12 @@ continent = np.unique(continent) else: - continent = [continent] + pass # Nothing to do: assume we have a list of continents #}}} - #expand basins list using the continent list above and the extra bas discriminator #{{{ + + # Expand basins list using the continent list above and the extra bas discriminator #{{{ if type(bas) == np.ndarray: if bas.shape[1] == 1: if bas[0] == 'all': - #need to transform this into a list of basins + # Need to transform this into a list of basins baslist = [] for i in range(len(self.basins)): @@ -246,5 +250,5 @@ baslist.append(i) else: - #we have a list of basin names + # We have a list of basin names baslist = [] for i in range(len(bas)): @@ -273,2 +277,140 @@ #}}} #}}} + + def addicecap(self, md): #{{{ + if not type(md) == model: + raise Exception("addicecap method only takes a 'model' class object as input") + + self.icecaps.append(md) + #}}} + + def basinsplot3d(self, *args): #{{{ + for i in range(len(self.basins)): + self.basins[i].plot3d(*args) + #}}} + + def caticecaps(self, *args): #{{{ + # Recover options + options = pairoptions(*args) + tolerance = options.getfieldvalue('tolerance', .65) + loneedgesdetect = options.getfieldvalue('loneedgesdetect', 0) + + # Make 3D model + models = self.icecaps + for i in range(len(models)): + models[i] = TwoDToThreeD(models[i], self.planet) + + # Plug all models together + md = models[0] + for i in range(1, len(models)): + md = modelmerge3d(md, models[i], 'tolerance', tolerance) + md.private.bamg.landmask = np.vstack((md.private.bamg.landmask, models[i].private.bamg.landmask)) + + # Look for lone edges if asked for it + if loneedgesdetect: + edges = loneedges(md) + plotmodel(md, 'data', md.mask.land_levelset) + for i in range(len(edges)): + ind1 = edges(i, 1) + ind1 = edges(i, 2) + 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*-') + + # Plug into earth + self.earth = md + + # Create mesh radius + self.earth.mesh.r = planetradius('earth') + #}}} + + def viscousiterations(self): #{{{ + for i in range(len(self.icecaps)): + ic = self.icecaps[i] + mvi = ic.resutls.TransientSolution[0].StressbalanceConvergenceNumSteps + for j in range(1, len(ic.results.TransientSolution) - 1): + mvi = np.max(mvi, ic.results.TransientSolution[j].StressbalanceConvergenceNumSteps) + print("{}, {}: {}".format(i, self.icecaps[i].miscellaneous.name, mvi)) + #}}} + + def maxtimestep(self): #{{{ + for i in range(len(self.icecaps)): + ic = self.icecaps[i] + mvi = len(ic.results.TransientSolution) + timei = ic.results.TransientSolution[-1].time + print("{}, {}: {}/{}".format(i, self.icecaps[i].miscellaneous.name, mvi, timei)) + + mvi = len(self.earth.results.TransientSolution) + timei = self.earth.results.TransientSolution[-1].time + print("Earth: {}/{}", mvi, timei) + #}}} + + def transfer(self, string): #{{{ + # Recover field size in one icecap + n = np.size(getattr(self.icecaps[i], string), 0) + + if n == self.icecaps[0].mesh.numberofvertices: + setattr(self.earth, string, np.zeros((self.earth.mesh.numberofvertices, ))) + for i in range(len(self.icecaps)): + getattr(self.earth, string)[self.transitions[i]] = getattr(self.icecaps[i], string) + elif n == (self.self.icecaps[0].mesh.numberofvertices + 1): + # Dealing with transient dataset + # Check that all timetags are similar between all icecaps #{{{ + for i in range(len(self.icecaps)): + capfieldi = getattr(self.icecaps[i], string) + for j in range(1, len(self.icecaps)): + capfieldj = getattr(self.icecaps[j], string) + if capfieldi[-1, :] != capfieldj[-1, :]: + raise Exception("Time stamps for {} field is different between icecaps {} and {}".format(string, i, j)) + capfield1 = getattr(self.icecaps[0], string) + times = capfield1[-1, :] + nsteps = len(times) + #}}} + # Initialize #{{{ + field = np.zeros((self.earth.mesh.numberofvertices + 1, nsteps)) + field[-1, :] = times # Transfer the times only, not the values + #}}} + # Transfer all the time fields #{{{ + for i in range(len(self.icecaps)): + capfieldi = getattr(self.icecaps[i], string) + for j in range(nsteps): + field[self.transitions[i], j] = capfieldi[0:-2, j] # Transfer only the values, not the time + setattr(self.earth, string, field) # Do not forget to plug the field variable into its final location + #}}} + elif n == (self.icecaps[0].mesh.numberofelements): + setattr(self.earth, string, np.zeros((self.earth.mesh.numberofvertices, ))) + for i in range(len(self.icecaps)): + getattr(self.earth, string)[self.eltransitions[i]] = getattr(self.icecaps[i], string) + else: + raise Exception('not supported yet') + #}}} + + def homogenize(self, noearth=0): #{{{ + mintimestep = np.inf + + for i in range(len(self.icecaps)): + ic = self.icecaps[i] + mintimestep = np.min(mintimestep, len (ic.results.TransientSolution)) + + if not noearth: + mintimestep = np.min(mintimestep, len(self.earth.results.TransientSolution)) + + for i in range(len(self.icecaps)): + ic = self.icecaps[i] + ic.resuts.TransientSolution = ic.results.TransientSolution[:mintimestep] + self.icecaps[i] = ic + + ic = self.earth + + if not noearth: + ic.results.TransientSolution = ic.resutls.TransientSolution[:mintimestep] + + self.earth = ic + + return self + #}}} + + def initializemodels(self): #{{{ + for i in range(len(self.basins)): + md = model() + md.miscellaneous.name = self.basins[i].name + self.addicecap(md) + #}}} Index: /issm/trunk-jpl/src/m/coordsystems/epsg2proj.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/epsg2proj.m (revision 25454) +++ /issm/trunk-jpl/src/m/coordsystems/epsg2proj.m (revision 25455) @@ -1,4 +1,4 @@ function string = epsg2proj(epsg) -%FUNCTION EPSG2PROJ - uses gdalsrsinfo to provide PROJ.4 compatible string from +%EPSG2PROJ - uses gdalsrsinfo to provide PROJ.4 compatible string from %EPSG code % Index: /issm/trunk-jpl/src/m/coordsystems/epsg2proj.py =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/epsg2proj.py (revision 25454) +++ /issm/trunk-jpl/src/m/coordsystems/epsg2proj.py (revision 25455) @@ -1,9 +1,7 @@ -import shlex import subprocess def epsg2proj(epsg): #{{{ - """ - FUNCTION EPSG2PROJ - uses gdalsrsinfo to provide PROJ.4 compatible string + """EPSG2PROJ - uses gdalsrsinfo to provide PROJ.4 compatible string from EPSG code @@ -22,5 +20,5 @@ #First, get GDAL version - subproc_args = shlex.split("gdalsrsinfo --version | awk '{print $2}' | cut -d '.' -f1") + subproc_args = "gdalsrsinfo --version | awk '{print $2}' | cut -d '.' -f1" subproc = subprocess.Popen(subproc_args, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) outs, errs = subproc.communicate() @@ -30,5 +28,5 @@ version_major=int(outs) - subproc_args = shlex.split("gdalsrsinfo epsg:{} | command grep PROJ.4 | tr -d '\n' | sed 's/PROJ.4 : //'".format(epsg)) + subproc_args = "gdalsrsinfo epsg:{} | command grep PROJ.4 | tr -d '\n' | sed 's/PROJ.4 : //'".format(epsg) subproc = subprocess.Popen(subproc_args, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) outs, errs = subproc.communicate() Index: /issm/trunk-jpl/src/m/coordsystems/gdaltransform.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/gdaltransform.m (revision 25454) +++ /issm/trunk-jpl/src/m/coordsystems/gdaltransform.m (revision 25455) @@ -21,5 +21,4 @@ % % To get proj.4 string from EPSG, use gdalsrsinfo. Example: -% % gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //" @@ -32,13 +31,17 @@ fclose(fid); - [s,r]=system(['gdaltransform -s_srs "',proj_in,'" -t_srs "',proj_out,'" < ' filename_in ' > ' filename_out]); + [s,r]=system(['gdaltransform -s_srs "',proj_in,'" -t_srs "',proj_out,'" < ' filename_in ' > ' filename_out]); if s~=0 | ~isempty(deblank(r)), error(r); end + A=load(filename_out); - xout=A(:,1); xout=reshape(xout,size(x)); - yout=A(:,2); yout=reshape(yout,size(y)); %clean up delete(filename_in); delete(filename_out); + + xout=A(:,1); + xout=reshape(xout,size(x)); + yout=A(:,2); + yout=reshape(yout,size(y)); Index: /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py (revision 25454) +++ /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py (revision 25455) @@ -4,10 +4,11 @@ import tempfile +import numpy as np + from loadvars import * def gdaltransform(x, y, proj_in, proj_out): #{{{ - ''' - GDALTRANSFORM - switch from one projection system to another + """GDALTRANSFORM - switch from one projection system to another Usage: @@ -29,32 +30,42 @@ +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 - To get proj.4 string form EPSG, use gdalsrsinfo. Example: - + To get proj.4 string from EPSG, use gdalsrsinfo. Example: gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //" - - TODO: - - Follow subprocess pattern implemented in src/m/coordsystems/epsg2proj.py - ''' + """ # Give ourselves unique file names - file_in = tempfile.NamedTemporaryFile('r+b') + file_in = tempfile.NamedTemporaryFile('w', delete=False) filename_in = file_in.name - file_out = tempfile.NamedTemporaryFile('w+b') + file_out = tempfile.NamedTemporaryFile('w', delete=False) filename_out = file_out.name - file_in.write(b'%8g %8g\n' % (x.flatten(1), y.flatten(1))) + points = np.vstack((x, y)).T + np.savetxt(file_in, points, fmt='%8g %8g') + 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. + file_in = open(filename_in) # Open for reading by subprocess + + subproc_args = shlex.split("gdaltransform -s_srs '{}' -t_srs '{}'".format(proj_in, proj_out)) + subproc = subprocess.Popen(subproc_args, bufsize=-1, stdin=file_in, stdout=file_out, stderr=subprocess.PIPE, close_fds=True) + outs, errs = subproc.communicate() + if errs != '': + raise RuntimeError("gdaltransform: call to gdaltransform failed: {}".format(errs)) + + A = np.loadtxt(filename_out) + + # Clean up file_in.close() + file_out.close() + remove(filename_in) + remove(filename_out) - subproc_args = shlex.split('gdaltransform -s_srs {} -t_srs {} < {} > {}'.format(proj_in, proj_out, filename_in, filename_out)) - subprocess.check_call(subproc_args) # Will raise CalledProcessError if return code is not 0 - - A = loadvars(filename_out) - xout = A[0] - xout = xout.reshape(x.shape) - yout = A[1] - yout = yout.reshape(y.shape) - - os.remove(filename_in) - os.remove(filename_out) + if np.ndim(A) > 1: + xout = np.array(A[:,0]) + yout = np.array(A[:,1]) + # if type(x) == "np.ndarray": + # xout = xout.reshape(x.shape) # TODO: Do we need to do this? + # yout = yout.reshape(y.shape) # TODO: Do we need to do this? + else: + xout = [A[0]] + yout = [A[1]] return [xout, yout] Index: /issm/trunk-jpl/src/m/coordsystems/laea.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/laea.m (revision 25454) +++ /issm/trunk-jpl/src/m/coordsystems/laea.m (revision 25455) @@ -9,3 +9,3 @@ % return string='+proj=laea +lat_0=45 +lon_0=-90 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs' - string=sprintf('+proj=laea +lat_0=%i +lon_0=%i +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs ',lat,long); + string=sprintf('+proj=laea +lat_0=%i +lon_0=%i +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs',lat,long); Index: /issm/trunk-jpl/src/m/exp/expread.py =================================================================== --- /issm/trunk-jpl/src/m/exp/expread.py (revision 25454) +++ /issm/trunk-jpl/src/m/exp/expread.py (revision 25455) @@ -1,30 +1,43 @@ +from collections import OrderedDict import os.path + import numpy as np -from collections import OrderedDict + import MatlabFuncs as m def expread(filename): + """EXPREAD - read a exp file and build a list of OrderedDicts + + This routine reads a file .exp and builds a list of OrderedDicts containing + the fields x and y corresponding to the coordinates, one for the filename + of the exp file, for the density, for the nodes, and a field closed to + indicate if the domain is closed. The first argument is the .exp file to be + read and the second one (optional) indicate if the last point shall be read + (1 to read it, 0 not to). + + Usage: + contours = expread(filename) + + Example: + contours = expread('domainoutline.exp') + contours = expread('domainoutline.exp') + + See Also: + - EXPDOC + - EXPWRITEASVERTICES + + TODO: + - Convert returned data structure from list of OrderedDict objects to list + of OrderedStruct objects (see also src/m/shp/shpread.py). Also, modify + handling of returned data structure in, + - src/m/exp/expcoarsen.py + - src/m/exp/expdisp.py + - src/m/interp/SectionValues.py + - src/m/mesh/bamg.py + May also need to modify addressing in corresponding FetchData function, or + create new one, in src/wrappers/ContoursToNodes/ContoursToNodes.cpp. """ - EXPREAD - read a file exp and build a Structure - - This routine reads a file .exp and builds a list of dicts containing the - fields x and y corresponding to the coordinates, one for the filename of - the exp file, for the density, for the nodes, and a field closed to - indicate if the domain is closed. - The first argument is the .exp file to be read and the second one (optional) - indicate if the last point shall be read (1 to read it, 0 not to). - - Usage: - contours = expread(filename) - - Example: - contours = expread('domainoutline.exp') - contours = expread('domainoutline.exp') - - See also EXPDOC, EXPWRITEASVERTICES - - """ #some checks if not os.path.exists(filename): Index: /issm/trunk-jpl/src/m/geometry/delaunay.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/delaunay.py (revision 25455) +++ /issm/trunk-jpl/src/m/geometry/delaunay.py (revision 25455) @@ -0,0 +1,39 @@ +import numpy as np +from scipy.spatial import Delaunay # See also: https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html?highlight=delaunay#scipy.spatial.Delaunay + + +def delaunay(*args): + """DELAUNAY - Delaunay triangulation + + Replicates behaviour of MATLAB's 'delaunay'. + + Usage: + # Creates a 2-D or 3-D Delaunay triangulation from the points in a matrix P. The output DT is a three-column (for two dimensions) or four-column (for three dimensions) matrix where each row contains the row indices of the input points that make up a triangle or tetrahedron in the triangulation. + DT = delaunay(P) + + # Creates a 2-D Delaunay triangulation from the points in vectors x and y. + DT = delaunay(x, y) + + # Creates a 3-D Delaunay triangulation from the points in vectors x, y, and z. + DT = delaunay(x, y, z) + + Sources: + - https://www.mathworks.com/help/matlab/ref/delaunay.html + """ + + nargs = len(args) + + if nargs == 1: + points = args[0] + elif nargs == 2: + points = np.vstack((args[0], args[1])).T + elif nargs == 3: + # NOTE: Not tested, but it is assumed that 3-D triangulation would work. + points = np.vstack((args[0], args[1], args[1])).T + else: + print(delaunay.__doc__) + raise Exception('Wrong usage (see above)') + + simplices = Delaunay(points).simplices.astype(int) # NOTE: Need to covert to array of int so that it can be fetched properly by core + + return simplices Index: /issm/trunk-jpl/src/m/geometry/find_point.m =================================================================== --- /issm/trunk-jpl/src/m/geometry/find_point.m (revision 25454) +++ /issm/trunk-jpl/src/m/geometry/find_point.m (revision 25455) @@ -2,6 +2,5 @@ %FIND_POINT - find closest point % -% find which point of the list (tabx,taby) is -% the closest to (poinx,pointy) +% find which point of the list (tabx,taby) is the closest to (pointx,pointy) % % Usage: Index: /issm/trunk-jpl/src/m/geometry/find_point.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/find_point.py (revision 25455) +++ /issm/trunk-jpl/src/m/geometry/find_point.py (revision 25455) @@ -0,0 +1,19 @@ +import numpy as np + + +def find_point(tabx, taby, pointx, pointy): + """FIND_POINT - find closest point + + Find which point of the list (tabx, taby) is the closest to (pointx, pointy) + + Usage: + f = find_point(tabx, taby, pointx, pointy) + """ + + # Compute distance between point and cloud of points + distance = ((tabx - pointx) ** 2 + (taby - pointy) ** 2) ** 0.5 # NOTE: math.sqrt cannot be applied element-wise to a list/numpy.array + + # Find index of the minimum distance and return the first one only + f = np.where(distance == np.min(distance))[0][0] # NOTE: numpy.where returns a tuple whose first element is a list, and we want the first element of that list + + return f Index: /issm/trunk-jpl/src/m/geometry/inpolygon.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/inpolygon.py (revision 25455) +++ /issm/trunk-jpl/src/m/geometry/inpolygon.py (revision 25455) @@ -0,0 +1,27 @@ +import matplotlib.path as path +import numpy as np + +def inpolygon(xq, yq, xv, yv): #{{{ + """ + INPOLYGON - Returns points located inside polygonal region. + + Partial implementation of MATLAB's inpolygon function. Returns in + indicating if the query points specified by xq and yq are inside or on the + edge of the polygon area defined by xv and yv. + + Usage: + in_polygon = inpolygon(xq, yq, xv, yv) + + Sources: + - https://www.mathworks.com/help/matlab/ref/inpolygon.html + - https://stackoverflow.com/questions/31542843/inpolygon-for-python-examples-of-matplotlib-path-path-contains-points-method/31543337#31543337 + """ + + points = np.array((xq, yq)).T + vertices = np.array((xv, yv)).T + polygon = path.Path(vertices) + in_polygon = polygon.contains_points(points) + in_polygon = in_polygon.astype(int) # Convert from bool to int + + return in_polygon +#}}} Index: /issm/trunk-jpl/src/m/geometry/polyarea.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/polyarea.py (revision 25454) +++ /issm/trunk-jpl/src/m/geometry/polyarea.py (revision 25455) @@ -5,5 +5,5 @@ def polyarea(x, y): #{{{ - ''' + """ POLYAREA - returns the area of the 2-D polygon defined by the vertices in lists x and y @@ -24,5 +24,6 @@ - Test that output falls within some tolerance of MATLAB's polyarea function. - ''' + """ + return 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1))) #}}} Index: /issm/trunk-jpl/src/m/interp/SectionValues.py =================================================================== --- /issm/trunk-jpl/src/m/interp/SectionValues.py (revision 25454) +++ /issm/trunk-jpl/src/m/interp/SectionValues.py (revision 25455) @@ -1,14 +1,15 @@ import os + +import numpy as np + from expread import expread -import numpy as np -from project2d import project2d #from InterpFromMesh2d import InterpFromMesh2d from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d from InterpFromMeshToMesh3d import InterpFromMeshToMesh3d +from project2d import project2d def SectionValues(md, data, infile, resolution): - ''' - compute the value of a field on a section + """SECTIONVALUES - compute the value of a field on a section This routine gets the value of a given field of the model on points @@ -17,7 +18,7 @@ Usage: - [elements, x, y, z, s, data] = SectionValues(md, data, filename, resolution) - [elements, x, y, z, s, data] = SectionValues(md, data, profile_structure, resolution) - ''' + [elements, x, y, z, s, data] = SectionValues(md, data, filename, resolution) + [elements, x, y, z, s, data] = SectionValues(md, data, profile_structure, resolution) + """ if os.path.isfile(infile): @@ -84,6 +85,6 @@ #Interpolation of data on specified points - #data_interp = InterpFromMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y, data, X, Y)[0] - data_interp = InterpFromMeshToMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y, data, X, Y)[0] + #data_interp = InterpFromMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y, data, X, Y) + data_interp = InterpFromMeshToMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y, data, X, Y) #data_interp = griddata(md.mesh.x, md.mesh.y, data, X, Y) @@ -95,7 +96,7 @@ #Get base and surface for each 2d point, offset to make sure that it is inside the glacier system offset = 1.e-3 - base = InterpFromMeshToMesh2d(md.mesh.elements2d, md.mesh.x2d, md.mesh.y2d, project2d(md, md.geometry.base, 1), X, Y)[0] + offset + base = InterpFromMeshToMesh2d(md.mesh.elements2d, md.mesh.x2d, md.mesh.y2d, project2d(md, md.geometry.base, 1), X, Y) + offset base = base.reshape(-1, ) - surface = InterpFromMeshToMesh2d(md.mesh.elements2d, md.mesh.x2d, md.mesh.y2d, project2d(md, md.geometry.surface, 1), X, Y)[0] - offset + surface = InterpFromMeshToMesh2d(md.mesh.elements2d, md.mesh.x2d, md.mesh.y2d, project2d(md, md.geometry.surface, 1), X, Y) - offset surface = surface.reshape(-1, ) @@ -126,5 +127,5 @@ #Interpolation of data on specified points - data_interp = InterpFromMeshToMesh3d(md.mesh.elements, md.mesh.x, md.mesh.y, md.mesh.z, data, X3, Y3, Z3, np.nan)[0] + data_interp = InterpFromMeshToMesh3d(md.mesh.elements, md.mesh.x, md.mesh.y, md.mesh.z, data, X3, Y3, Z3, np.nan) #build outputs Index: /issm/trunk-jpl/src/m/io/loadmodel.m =================================================================== --- /issm/trunk-jpl/src/m/io/loadmodel.m (revision 25454) +++ /issm/trunk-jpl/src/m/io/loadmodel.m (revision 25455) @@ -1,4 +1,4 @@ function varargout=loadmodel(path) -%LOADMODEL - load a model using built-in load module +%LOADMODEL - load a model using built-in load module % % check that model prototype has not changed. if so, adapt to new model prototype. Index: /issm/trunk-jpl/src/m/io/loadmodel.py =================================================================== --- /issm/trunk-jpl/src/m/io/loadmodel.py (revision 25454) +++ /issm/trunk-jpl/src/m/io/loadmodel.py (revision 25455) @@ -1,21 +1,19 @@ from loadvars import loadvars -#hack to keep python 2 compatipility +# Hack to keep python 2 compatibility try: - #py3 import - from dbm import whichdb + from dbm import whichdb # Python 3 except ImportError: - #py2 import - from whichdb import whichdb + from whichdb import whichdb # Python 2 + from netCDF4 import Dataset def loadmodel(path, onlylast=False): - """ - LOADMODEL - load a model using built - in load module + """LOADMODEL - load a model - check that model prototype has not changed. if so, adapt to new model prototype. + Check that model prototype has not changed: if if has, adapt to new model prototype. - Usage: - md = loadmodel(path) + Usage: + md = loadmodel(path) """ Index: /issm/trunk-jpl/src/m/io/loadvars.py =================================================================== --- /issm/trunk-jpl/src/m/io/loadvars.py (revision 25454) +++ /issm/trunk-jpl/src/m/io/loadvars.py (revision 25455) @@ -1,21 +1,19 @@ from collections import OrderedDict +# Hack to keep python 2 compatibility +try: + from dbm import whichdb # Python 3 +except ImportError: + from whichdb import whichdb # Python 2 +from re import findall +import shelve + from netCDF4 import Dataset import numpy as np -from re import findall -import shelve from model import * -#hack to keep python 2 compatibility -try: - #py3 import - from dbm import whichdb -except ImportError: - #py2 import - from whichdb import whichdb def loadvars(*args, **kwargs): - """ - LOADVARS - function to load variables from a file. + """LOADVARS - function to load variables from a file This function loads one or more variables from a file. The names of the Index: /issm/trunk-jpl/src/m/mech/newforcing.m =================================================================== --- /issm/trunk-jpl/src/m/mech/newforcing.m (revision 25454) +++ /issm/trunk-jpl/src/m/mech/newforcing.m (revision 25455) @@ -1,16 +1,16 @@ function forcing=newforcing(t0,t1,deltaT,f0,f1,nodes) -%FUNCTION NEWFORCING Build forcing that extends temporally from t0 to t1, and in magnitude from f0 to f1. Equal time -% and magnitude spacing. +%NEWFORCING - Build forcing that extends temporally from t0 to t1, and in +%magnitude from f0 to f1. Equal time and magnitude spacing. % -% Usage: forcing=newforcing(t0,t1,deltaT,f0,f1,nodes); -% Where: -% t0:t1: time interval. -% deltaT: time step -% f0:f1: magnitude interval. -% nodes: number of vertices where we have a temporal forcing +% Usage: forcing=newforcing(t0,t1,deltaT,f0,f1,nodes); +% +% Where: +% t0:t1: time interval. +% deltaT: time step +% f0:f1: magnitude interval. +% nodes: number of vertices where we have a temporal forcing % -% Example: -% md.smb.mass_balance=newforcing(md.timestepping.start_time,md.timestepping.final_time,... -% md.timestepping.time_step,-1,+2,md.mesh.numberofvertices); +% Example: +% md.smb.mass_balance=newforcing(md.timestepping.start_time,md.timestepping.final_time,md.timestepping.time_step,-1,+2,md.mesh.numberofvertices); % Index: /issm/trunk-jpl/src/m/mech/newforcing.py =================================================================== --- /issm/trunk-jpl/src/m/mech/newforcing.py (revision 25454) +++ /issm/trunk-jpl/src/m/mech/newforcing.py (revision 25455) @@ -4,18 +4,26 @@ def newforcing(t0, t1, deltaT, f0, f1, nodes): ''' -FUNCTION NEWFORCING Build forcing that extends temporally from t0 to t1, and in magnitude from f0 to f1. Equal time + NEWFORCING - Build forcing that extends temporally from t0 to t1, and in magnitude from f0 to f1. Equal time and magnitude spacing. - Usage: forcing = newforcing(t0, t1, deltaT, f0, f1, nodes); - Where: - t0:t1: time interval. - deltaT: time step - f0:f1: magnitude interval. - nodes: number of vertices where we have a temporal forcing + Usage: forcing = newforcing(t0, t1, deltaT, f0, f1, nodes); + + Where: + t0:t1: time interval. + deltaT: time step + f0:f1: magnitude interval. + nodes: number of vertices where we have a temporal forcing - Example: - md.smb.mass_balance = newforcing(md.timestepping.start_time, md.timestepping.final_time, - md.timestepping.time_step, -1, +2, md.mesh.numberofvertices) + Example: + md.smb.mass_balance = newforcing( + md.timestepping.start_time, + md.timestepping.final_time, + md.timestepping.time_step, + -1, + +2, + md.mesh.numberofvertices + ) ''' + #Number of time steps: nsteps = (t1 - t0) / deltaT + 1 Index: /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.m (revision 25454) +++ /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.m (revision 25455) @@ -21,5 +21,5 @@ vc=md.mesh.vertexconnectivity; vb=md.mesh.vertexonboundary; - md.mesh=mesh3dsurface(); + md.mesh=mesh3dsurface(); md.mesh.lat=lat; md.mesh.long=long; Index: /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.py (revision 25454) +++ /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.py (revision 25455) @@ -1,8 +1,11 @@ import numpy as np +from epsg2proj import epsg2proj +from gdaltransform import gdaltransform from mesh3dsurface import mesh3dsurface from planetradius import planetradius -def TwoDToThreeD(md, planet): #{{{ + +def TwoDToThreeD(md, planet): # Reproject model into lat, long if necessary if md.mesh.proj != epsg2proj(4326): @@ -13,18 +16,18 @@ # We assume x and y hold the long, lat values - long = md.mesh.x - lat = md.mesh.y + longe = md.mesh.x + late = md.mesh.y # Assume spherical body - x = R * np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(long)) - y = R * np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(long)) - z = R * np.sin(np.deg2rad(lat)) + x = R * np.cos(np.deg2rad(late)) * np.cos(np.deg2rad(longe)) + y = R * np.cos(np.deg2rad(late)) * np.sin(np.deg2rad(longe)) + z = R * np.sin(np.deg2rad(late)) elements = md.mesh.elements vc = md.mesh.vertexconnectivity vb = md.mesh.vertexonboundary - md.mesh.mesh3dsurface() - md.mesh.lat = lat - md.mesh.long = long + md.mesh = mesh3dsurface() + md.mesh.lat = late + md.mesh.long = longe md.mesh.x = x md.mesh.y = y @@ -32,7 +35,8 @@ md.mesh.elements = elements md.mesh.numberofelements = len(elements) - md.mesh.numberofvertices = len(lat) - md.mesh.r = R * np.ones(md.mesh.numberofvertices) + md.mesh.numberofvertices = len(late) + md.mesh.r = R * np.ones((md.mesh.numberofvertices, )) md.mesh.vertexconnectivity = vc md.mesh.vertexonboundary = vb -#}}} + + return md Index: /issm/trunk-jpl/src/m/mesh/bamg.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.m (revision 25454) +++ /issm/trunk-jpl/src/m/mesh/bamg.m (revision 25455) @@ -80,5 +80,5 @@ domain=shpread(domainfile); else - error(['bamg error message: file ' domainfile ' format not supported (.shp or .exp)']); + error(['bamg error message: file ' domainfile ' format not supported (.exp or .shp)']); end elseif isstruct(domainfile), @@ -140,5 +140,5 @@ end - %Checks that all holes are INSIDE the principle domain outline + %Check that all holes are INSIDE the principle domain outline if i>1, flags=ContourToNodes(domain(i).x,domain(i).y,domain(1),0); @@ -149,5 +149,5 @@ %Check orientation - nods=domain(i).nods-1; %the domain are closed 1=end; + nods=domain(i).nods-1; %the domain is closed (domain[1] = domain[end]) test = sum([(domain(i).x(2:nods+1) - domain(i).x(1:nods)).*(domain(i).y(2:nods+1) + domain(i).y(1:nods))]); if (i==1 && test>0) || (i>1 && test<0), @@ -163,11 +163,13 @@ bamg_geometry.Edges =[bamg_geometry.Edges; [transpose(count+1:count+nods) transpose([count+2:count+nods count+1]) 1*ones(nods,1)]]; - %Flag how many edges we have now, that way we know which edges belong to the domain, will - %be used later for required edges when NoBoundaryRefinement is one: + % Flag how many edges we have now, that way we know which edges belong + % to the domain. Will be used later for required edges if + % NoBoundaryRefinement equals 1. new_edge_length=length(bamg_geometry.Edges); edges_required=(edge_length+1):new_edge_length; if i>1, - bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 -subdomain_ref]; subdomain_ref = subdomain_ref+1; + bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 -subdomain_ref]; + subdomain_ref = subdomain_ref+1; else bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 0]; @@ -179,15 +181,15 @@ for i=1:length(holes), - %Check that the subdomains is closed + %Check that the hole is closed if (holes(i).x(1)~=holes(i).x(end) | holes(i).y(1)~=holes(i).y(end)), error('bamg error message: all contours provided in ''holes'' should be closed'); end - %Checks that all holes are INSIDE the principle domain outline + %Check that all holes are INSIDE the principal domain (principal domain should be index 0) flags=ContourToNodes(holes(i).x,holes(i).y,domain(1),0); if any(~flags), error('bamg error message: All holes should be strictly inside the principal domain'); end %Check that hole is correctly oriented - nods=holes(i).nods-1; %the holes are closed 1=end; + nods=holes(i).nods-1; %the hole is closed (hole[1] = hole[end]) 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 disp('At least one hole was not correctly oriented and has been re-oriented'); @@ -205,20 +207,21 @@ for i=1:length(subdomains), - %Check that the subdomains is closed + %Check that the subdomain is closed if (subdomains(i).x(1)~=subdomains(i).x(end) | subdomains(i).y(1)~=subdomains(i).y(end)), error('bamg error message: all contours provided in ''subdomains'' should be closed'); end - %Checks that all holes are INSIDE the principle domain outline + %Checks that all subdomains are INSIDE the principal domain (principal domain should be index 0) flags=ContourToNodes(subdomains(i).x,subdomains(i).y,domain(1),0); if any(~flags), - error('bamg error message: All holes should be strictly inside the principal domain'); - end - - %Check that hole is correctly oriented - nods=subdomains(i).nods-1; %the subdomains are closed 1=end; + error('bamg error message: All subdomains should be strictly inside the principal domain'); + end + + %Check that subdomain is correctly oriented + nods=subdomains(i).nods-1; % the subdomains are closed (subdomains[1] = subdomains[end]) 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 disp('At least one subdomain was not correctly oriented and has been re-oriented'); - subdomains(i).x = flipud(subdomains(i).x); subdomains(i).y = flipud(subdomains(i).y); + subdomains(i).x = flipud(subdomains(i).x); + subdomains(i).y = flipud(subdomains(i).y); end @@ -227,9 +230,11 @@ bamg_geometry.Edges =[bamg_geometry.Edges; [transpose(count+1:count+nods) transpose([count+2:count+nods count+1]) 1.*ones(nods,1)]]; - bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 subdomain_ref]; subdomain_ref = subdomain_ref+1; + bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 subdomain_ref]; + subdomain_ref = subdomain_ref+1; %update counter count=count+nods; end + if getfieldvalue(options,'vertical',0), if numel(getfieldvalue(options,'Markers',[]))~=size(bamg_geometry.Edges,1), @@ -243,16 +248,6 @@ %take care of rifts if exist(options,'rifts'), - - %Check that file exists + %read rift file according to its extension: riftfile=getfieldvalue(options,'rifts'); - [pathr,namer,extr]=fileparts(riftfile); - if ~exist(riftfile,'file') - error(['bamg error message: file ' riftfile ' not found ']); - elseif strcmp(extr,'.exp'), - rift=expread(riftfile); - elseif strcmp(extr,'.shp'), - rift=shpread(riftfile); - end - %read rift file according to its extension: [path,name,ext]=fileparts(riftfile); if strcmp(ext,'.exp'), @@ -261,5 +256,5 @@ rift=shpread(riftfile); else - error(['bamg error message: file ' riftfile ' format not supported (.shp or .exp)']); + error(['bamg error message: file ' riftfile ' format not supported (.exp or .shp)']); end @@ -272,10 +267,10 @@ elseif any(~flags), - %We LOTS of work to do + %We have LOTS of work to do disp('Rift tip outside of or on the domain has been detected and is being processed...'); %check that only one point is outside (for now) if sum(~flags)~=1, - error('bamg error message: only one point outside of the domain is supported yet'); + error('bamg error message: only one point outside of the domain is supported at this time'); end @@ -290,7 +285,9 @@ end - %Get cordinate of intersection point - x1=rift(i).x(1); y1=rift(i).y(1); - x2=rift(i).x(2); y2=rift(i).y(2); + %Get coordinate of intersection point + x1=rift(i).x(1); + y1=rift(i).y(1); + x2=rift(i).x(2); + y2=rift(i).y(2); for j=1:length(domain(1).x)-1; if SegIntersect([x1 y1; x2 y2],[domain(1).x(j) domain(1).y(j); domain(1).x(j+1) domain(1).y(j+1)]), @@ -489,8 +486,6 @@ md.mesh.numberofvertices=length(md.mesh.x); md.mesh.numberofedges=size(md.mesh.edges,1); - md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; - - %Now, build the connectivity tables for this mesh. - md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; + md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); + md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; elseif getfieldvalue(options,'3dsurface',0), @@ -509,5 +504,6 @@ md.mesh.numberofvertices=length(md.mesh.x); md.mesh.numberofedges=size(md.mesh.edges,1); - md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; + md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); + md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; else @@ -524,5 +520,6 @@ md.mesh.numberofvertices=length(md.mesh.x); md.mesh.numberofedges=size(md.mesh.edges,1); - md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; + md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); + md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; end Index: /issm/trunk-jpl/src/m/mesh/bamg.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.py (revision 25454) +++ /issm/trunk-jpl/src/m/mesh/bamg.py (revision 25455) @@ -1,21 +1,23 @@ +from collections import namedtuple, OrderedDict import os.path + import numpy as np + +from bamggeom import bamggeom +from BamgMesher import BamgMesher +from bamgmesh import bamgmesh +from ContourToNodes import ContourToNodes +from expread import expread +from helpers import fileparts, OrderedStruct from mesh2d import * from mesh2dvertical import * from mesh3dsurface import * -from collections import OrderedDict from pairoptions import pairoptions -from bamggeom import bamggeom -from bamgmesh import bamgmesh -from expread import expread from SegIntersect import SegIntersect -import MatlabFuncs as m -from BamgMesher import BamgMesher -from ContourToNodes import ContourToNodes +from shpread import shpread def bamg(md, *args): - """ - BAMG - mesh generation + """BAMG - mesh generation Available options (for more details see ISSM website http://issm.jpl.nasa.gov/): @@ -43,4 +45,6 @@ 1 -> use Green formula - KeepVertices : try to keep initial vertices when adaptation is done on an existing mesh (default 1) + - NoBoundaryRefinement: do not refine boundary, only follow contour provided (default 0). Allow subdomain boundary refinement though + - NoBoundaryRefinementAllBoundaries: do not refine boundary, only follow contour provided (default 0) - maxnbv : maximum number of vertices used to allocate memory (default is 1.0e6) - maxsubdiv : maximum subdivision of exisiting elements (default is 10) @@ -63,8 +67,11 @@ will be merged - Examples: - md = bamg(md, 'domain', 'DomainOutline.exp', 'hmax', 3000) - md = bamg(md, 'field', [md.inversion.vel_obs md.geometry.thickness], 'hmax', 20000, 'hmin', 1000) - md = bamg(md, 'metric', A, 'hmin', 1000, 'hmax', 20000, 'gradation', 3, 'anisomax', 1) + Examples: + md = bamg(md, 'domain', 'DomainOutline.exp', 'hmax', 3000) + md = bamg(md, 'field', [md.inversion.vel_obs md.geometry.thickness], 'hmax', 20000, 'hmin', 1000) + md = bamg(md, 'metric', A, 'hmin', 1000, 'hmax', 20000, 'gradation', 3, 'anisomax', 1) + + TODO: + - Verify that values of third column of bamg_geometry.Vertices and bamg_geometry.Edges are valid (compare to src/m/mesh/bamg.m) """ @@ -79,4 +86,5 @@ subdomain_ref = 1 + hole_ref = 1 # Bamg Geometry parameters {{{ @@ -86,35 +94,120 @@ if type(domainfile) == str: if not os.path.exists(domainfile): - raise IOError("bamg error message: file '%s' not found" % domainfile) + raise IOError("bamg error message: file {} not found".format( domainfile)) + + # Read domain according to its extension + path, name, ext = fileparts(domainfile) + if ext == '.exp': + domain = expread(domainfile) + elif ext == '.shp': + domain = shpread(domainfile) + else: + raise Exception('bamg error message: file {} format not supported (.exp or .shp)'.format(domainfile)) domain = expread(domainfile) + elif type(domainfile) == list: + if len(domainfile): + if type(domainfile[0]) in [dict, OrderedDict]: + domain = domainfile + else: + raise Exception("bamg error message: if 'domain' is a list, its elements must be of type dict or OrderedDict") + else: + domain = domainfile + elif type(domainfile) in [dict, OrderedDict]: + domain = [domainfile] else: - domain = domainfile - - #Build geometry + raise Exception("bamg error message: 'domain' type {} not supported yet".format(type(domainfile))) + + holes = [] + if options.exist('holes'): + holesfile = options.getfieldvalue('holes') + if type(holesfile) == str: + if not os.path.exists(holesfile): + raise IOError("bamg error message: file {} not found".format(holesfile)) + + # Read holes accoridng to its extension + path, name, ext = fileparts(holesfile) + if ext == '.exp': + holes = expread(holesfile) + elif ext == '.shp': + holes = shpread(holesfile) + else: + raise Exception('bamg error message: file {} format not supported (.exp or .shp)'.format(holesfile)) + elif type(holesfile) == list: + if len(holesfile): + if type(holesfile[0]) in [dict, OrderedDict]: + holes = holesfile + else: + raise Exception("bamg error message: if 'holes' is a list, its elements must be of type dict or OrderedDict") + else: + holes = holesfile + elif type(holesfile) in [dict, OrderedDict]: + holes = [holesfile] + else: + raise Exception("'holes' type {} not supported yet".format(type(holesfile))) + + subdomains = [] + if options.exist('subdomains'): + subdomainsfile = options.getfieldvalue('subdomains') + if type(subdomainsfile) == str: + if not os.path.exists(subdomainsfile): + raise IOError("bamg error message: file {} not found".format(subdomainsfile)) + + # Read subdomains accoridng to its extension + path, name, ext = fileparts(subdomainsfile) + if ext == '.exp': + subdomains = expread(subdomainsfile) + elif ext == '.shp': + subdomains = shpread(subdomainsfile) + else: + raise Exception('bamg error message: file {} format not supported (.exp or .shp)'.format(subdomainsfile)) + elif type(subdomainsfile) == list: + if len(subdomainsfile): + if type(subdomainsfile[0]) in [dict, OrderedDict]: + subdomains = subdomainsfile + else: + raise Exception("bamg error message: if 'subdomains' is a list, its elements must be of type dict or OrderedDict") + else: + subdomains = subdomainsfile + elif type(subdomainsfile) in [dict, OrderedDict]: + subdomains = [subdomainsfile] + else: + raise Exception("'subdomains' type {} not supported yet".format(type(subdomainsfile))) + + # Build geometry count = 0 - for i, domaini in enumerate(domain): - #Check that the domain is closed - if (domaini['x'][0] != domaini['x'][-1] or domaini['y'][0] != domaini['y'][-1]): + for i in range(len(domain)): + # Check that the domain is closed + if (domain[i]['x'][0] != domain[i]['x'][-1] or domain[i]['y'][0] != domain[i]['y'][-1]): raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed") - #Checks that all holes are INSIDE the principle domain outline princial domain should be index 0 + # Check that all holes are INSIDE the principle domain outline if i: - flags = ContourToNodes(domaini['x'], domaini['y'], domainfile, 0)[0] + 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]]' if np.any(np.logical_not(flags)): raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain") - #Check orientation - nods = domaini['nods'] - 1 #the domain are closed 1 = end - test = np.sum((domaini['x'][1:nods + 1] - domaini['x'][0:nods]) * (domaini['y'][1:nods + 1] + domaini['y'][0:nods])) + # Check orientation + nods = domain[i]['nods'] - 1 # the domain is closed (domain[0] = domain[-1]) + 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])) if (i == 0 and test > 0) or (i > 0 and test < 0): print('At least one contour was not correctly oriented and has been re-oriented') - domaini['x'] = np.flipud(domaini['x']) - domaini['y'] = np.flipud(domaini['y']) - - #Add all points to bamg_geometry - nods = domaini['nods'] - 1 #the domain are closed 0 = end - bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((domaini['x'][0:nods], domaini['y'][0:nods], np.ones((nods)))).T)) + domain[i]['x'] = np.flipud(domain[i]['x']) + domain[i]['y'] = np.flipud(domain[i]['y']) + + # Flag how many edges we have so far + edge_length = len(bamg_geometry.Edges) + + # Add all points to bamg_geometry + #nods = domain[i]['nods'] - 1 #the domain are closed 0 = end + bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((domain[i]['x'][0:nods], domain[i]['y'][0:nods], np.ones((nods)))).T)) 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)) - if i: + + # Flag how many edges we have now, that way we know which edges + # belong to the subdomain. Will be used later fo required edges + # if NoBoundaryRefinement equals 1. + new_edge_length = len(bamg_geometry.Edges) + 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) + + if i: # NOTE: same as `if i > 0` (MATLAB is `if i > 1`) bamg_geometry.SubDomains = np.vstack((bamg_geometry.SubDomains, [2, count + 1, 1, -subdomain_ref])) subdomain_ref = subdomain_ref + 1 @@ -122,139 +215,123 @@ bamg_geometry.SubDomains = np.vstack((bamg_geometry.SubDomains, [2, count + 1, 1, 0])) - #update counter + # Update counter count += nods - #Deal with domain holes - if options.exist('holes'): - holesfile = options.getfieldvalue('holes') - if type(holesfile) == str: - if not os.path.exists(holesfile): - raise IOError("bamg error message: file '%s' not found" % holesfile) - holes = expread(holesfile) - else: - holes = holesfile - - #Build geometry - for i, holei in enumerate(holes): - #Check that the hole is closed - if (holei['x'][0] != holei['x'][-1] or holei['y'][0] != holei['y'][-1]): - raise RuntimeError("bamg error message: all contours provided in ''hole'' should be closed") - - #Checks that all holes are INSIDE the principle domain outline princial domain should be index 0 - flags = ContourToNodes(holei['x'], holei['y'], domainfile, 0)[0] - if np.any(np.logical_not(flags)): - raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain") - - #Check orientation - nods = holei['nods'] - 1 #the hole are closed 1 = end - test = np.sum((holei['x'][1:nods + 1] - holei['x'][0:nods]) * (holei['y'][1:nods + 1] + holei['y'][0:nods])) - if test < 0: - print('At least one hole was not correctly oriented and has been re-oriented') - holei['x'] = np.flipud(holei['x']) - holei['y'] = np.flipud(holei['y']) - - #Add all points to bamg_geometry - nods = holei['nods'] - 1 #the hole are closed 0 = end - bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((holei['x'][0:nods], holei['y'][0:nods], np.ones((nods)))).T)) - 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)) - #update counter - count += nods - - #And subdomains - if options.exist('subdomains'): - subdomainfile = options.getfieldvalue('subdomains') - if type(subdomainfile) == str: - if not os.path.exists(subdomainfile): - raise IOError("bamg error message: file '{}' not found".format(subdomainfile)) - subdomains = expread(subdomainfile) - else: - subdomains = subdomainfile - - #Build geometry - for i, subdomaini in enumerate(subdomains): - #Check that the subdomain is closed - if (subdomaini['x'][0] != subdomaini['x'][-1] or subdomaini['y'][0] != subdomaini['y'][-1]): - raise RuntimeError("bamg error message: all contours provided in ''subdomain'' should be closed") - - #Checks that all subdomains are INSIDE the principle subdomain outline princial domain should be index 0 - if i: - flags = ContourToNodes(subdomaini['x'], subdomaini['y'], domainfile, 0)[0] - if np.any(np.logical_not(flags)): - raise RuntimeError("bamg error message: All subdomains should be strictly inside the principal subdomain") - - #Check orientation - nods = subdomaini['nods'] - 1 #the subdomain are closed 1 = end - - test = np.sum((subdomaini['x'][1:nods + 1] - subdomaini['x'][0:nods]) * (subdomaini['y'][1:nods + 1] + subdomaini['y'][0:nods])) - if test > 0: - print('At least one subcontour was not correctly oriented and has been re-oriented') - subdomaini['x'] = np.flipud(subdomaini['x']) - subdomaini['y'] = np.flipud(subdomaini['y']) - - #Add all points to bamg_geometry - nods = subdomaini['nods'] - 1 #the subdomain are closed 0 = end - bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((subdomaini['x'][0:nods], subdomaini['y'][0:nods], np.ones((nods)))).T)) - 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)) - #update counter - count += nods + for i in range(len(holes)): + # heck that the hole is closed + if (holes[i]['x'][0] != holes[i]['x'][-1] or holes[i]['y'][0] != holes[i]['y'][-1]): + raise RuntimeError("bamg error message: all contours provided in ''hole'' should be closed") + + # Check that all holes are INSIDE the principal domain (principal domain should be index 0) + 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]]' + if np.any(np.logical_not(flags)): + raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain") + + # Check that hole is correctly oriented + nods = holes[i]['nods'] - 1 # the holes are closed (holes[0] = holes[-1]) + 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])) + if test < 0: + print('At least one hole was not correctly oriented and has been re-oriented') + holes[i]['x'] = np.flipud(holes[i]['x']) + holes[i]['y'] = np.flipud(holes[i]['y']) + + # Add all points to bamg_geometry + #nods = holes[i]['nods'] - 1 # the hole is closed (hole[0] = hole[-1]) + bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((holes[i]['x'][0:nods], holes[i]['y'][0:nods], np.ones((nods)))).T)) + 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)) + bamg.geometry.SubDomains = np.vstack((bamg_geometry.SubDomains, [2, count + 1, 1, -hole_ref])) + hole_ref = hole_ref + 1 + + # Update counter + count += nods + + for i in range(len(subdomains)): + # Check that the subdomain is closed + if (subdomains[i]['x'][0] != subdomains[i]['x'][-1] or subdomains[i]['y'][0] != subdomains[i]['y'][-1]): + raise RuntimeError("bamg error message: all contours provided in ''subdomains'' should be closed") + + # Check that all subdomains are INSIDE the principal domain (principal domain should be index 0) + 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]]' + if np.any(np.logical_not(flags)): + raise RuntimeError("bamg error message: All subdomains should be strictly inside the principal domain") + + # Check that subdomain is correctly oriented + nods = subdomains[i]['nods'] - 1 # the subdomains are closed (subdomains[0] = subdomains[-1]) + + 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])) + if test: + print('At least one subcontour was not correctly oriented and has been re-oriented') + subdomains[i]['x'] = np.flipud(subdomains[i]['x']) + subdomains[i]['y'] = np.flipud(subdomains[i]['y']) + + #Add all points to bamg_geometry + #nods = subdomains[i]['nods'] - 1 # the subdomains are closed (subdomains[0] = subdomains[-1]) + bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((subdomains[i]['x'][0:nods], subdomains[i]['y'][0:nods], np.ones((nods)))).T)) + 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)) + bamg_geometry.SubDomains = np.vstack((bamg_geometry.SubDomains, [2, count + 1, 1, subdomain_ref])) + subdomain_ref = subdomain_ref + 1 + + # Update counter + count += nods if options.getfieldvalue('vertical', 0): if np.size(options.getfieldvalue('Markers', [])) != np.size(bamg_geometry.Edges, 0): raise RuntimeError('for 2d vertical mesh, ''Markers'' option is required, and should be of size ' + str(np.size(bamg_geometry.Edges, 0))) - if np.size(options.getfieldvalue('Markers', [])) == np.size(bamg_geometry.Edges, 0): - bamg_geometry.Edges[:, 2] = options.getfieldvalue('Markers') - - #take care of rifts + if np.size(options.getfieldvalue('Markers', [])) == np.size(bamg_geometry.Edges, 0): + bamg_geometry.Edges[:, 2] = options.getfieldvalue('Markers') + + # Take care of rifts if options.exist('rifts'): - - #Check that file exists + # Read rift file according to its extension riftfile = options.getfieldvalue('rifts') - if not os.path.exists(riftfile): - raise IOError("bamg error message: file '%s' not found" % riftfile) - rift = expread(riftfile) - - for i, rifti in enumerate(rift): - #detect whether all points of the rift are inside the domain - flags = ContourToNodes(rifti['x'], rifti['y'], domain[0], 0)[0] + path, name, ext = fileparts(riftfile) + if ext == '.exp': + rift = expread(riftfile) + elif ext == '.shp': + rift = shpread(riftfile) + else: + raise IOError("bamg error message: file '{}' format not supported (.exp or .shp)".format(riftfile)) + + for i in range(len(rift)): + # Detect whether all points of the rift are inside the domain + 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]]' if np.all(np.logical_not(flags)): raise RuntimeError("one rift has all its points outside of the domain outline") - elif np.any(np.logical_not(flags)): - #We LOTS of work to do + # We have LOTS of work to do print("Rift tip outside of or on the domain has been detected and is being processed...") - #check that only one point is outside (for now) + + # Check that only one point is outside (for now) if np.sum(np.logical_not(flags).astype(int)) != 1: - raise RuntimeError("bamg error message: only one point outside of the domain is supported yet") - - #Move tip outside to the first position + raise RuntimeError("bamg error message: only one point outside of the domain is supported at this time") + + # Move tip outside to the first position if not flags[0]: - #OK, first point is outside (do nothing), + # OK, first point is outside (do nothing), pass elif not flags[-1]: - rifti['x'] = np.flipud(rifti['x']) - rifti['y'] = np.flipud(rifti['y']) + rift[i]['x'] = np.flipud(rift[i]['x']) + rift[i]['y'] = np.flipud(rift[i]['y']) else: raise RuntimeError("bamg error message: only a rift tip can be outside of the domain") - #Get cordinate of intersection point - x1 = rifti['x'][0] - y1 = rifti['y'][0] - x2 = rifti['x'][1] - y2 = rifti['y'][1] + # Get coordinate of intersection point + x1 = rift[i]['x'][0] + y1 = rift[i]['y'][0] + x2 = rift[i]['x'][1] + y2 = rift[i]['y'][1] for j in range(0, np.size(domain[0]['x']) - 1): 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]]])): - #Get position of the two nodes of the edge in domain + # Get position of the two nodes of the edge in domain i1 = j i2 = j + 1 - #rift is crossing edge [i1 i2] of the domain - #Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) + # Rift is crossing edge [i1, i2] of the domain + # Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) x3 = domain[0]['x'][i1] y3 = domain[0]['y'][i1] x4 = domain[0]['x'][i2] y4 = domain[0]['y'][i2] - #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]) - #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]) 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]])) 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]])) @@ -264,7 +341,7 @@ if np.min(tipdis) / segdis < options.getfieldvalue('toltip', 0): - print("moving tip - domain intersection point") - - #Get position of the closer point + print("moving tip-domain intersection point") + + # Get position of the closer point if tipdis[0] > tipdis[1]: pos = i2 @@ -272,47 +349,84 @@ pos = i1 - #This point is only in Vertices (number pos). - #OK, now we can add our own rift - nods = rifti['nods'] - 1 - bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.hstack((rifti['x'][1:].reshape(-1, ), rifti['y'][1:].reshape(-1, ), np.ones((nods, 1)))))) - bamg_geometry.Edges = np.vstack((bamg_geometry.Edges, - np.array([[pos, count + 1, (1 + i)]]), - 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)))))) + # This point is only in Vertices (number pos). + # OK, now we can add our own rift + nods = rift[i]['nods'] - 1 + 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)))))) + bamg_geometry.Edges = np.vstack(( + bamg_geometry.Edges, + np.array([[pos, count + 1, (1 + i)]]), + 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)))) + )) count += nods break - else: - #Add intersection point to Vertices - bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.array([[x, y, 1]]))) + # Add intersection point to Vertices + bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, + np.array([[x, y, 1]]) + )) count += 1 - #Decompose the crossing edge into 2 subedges + # Decompose the crossing edge into 2 subedges pos = np.nonzero(np.logical_and(bamg_geometry.Edges[:, 0] == i1, bamg_geometry.Edges[:, 1] == i2))[0] if not pos: raise RuntimeError("bamg error message: a problem occurred...") - bamg_geometry.Edges = np.vstack((bamg_geometry.Edges[0:pos - 1, :], - np.array([[bamg_geometry.Edges[pos, 0], count, bamg_geometry.Edges[pos, 2]]]), - np.array([[count, bamg_geometry.Edges[pos, 1], bamg_geometry.Edges[pos, 2]]]), - bamg_geometry.Edges[pos + 1:, :])) - - #OK, now we can add our own rift - nods = rifti['nods'] - 1 - bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.hstack((rifti['x'][1:].reshape(-1, ), rifti['y'][1:].reshape(-1, ), np.ones((nods, 1)))))) - bamg_geometry.Edges = np.vstack((bamg_geometry.Edges, - np.array([[count, count + 1, 2]]), - 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)))))) + bamg_geometry.Edges = np.vstack(( + bamg_geometry.Edges[0:pos - 1, :], + np.array([[ + bamg_geometry.Edges[pos, 0], + count, + bamg_geometry.Edges[pos, 2] + ]]), + np.array([[ + count, + bamg_geometry.Edges[pos, 1], + bamg_geometry.Edges[pos, 2] + ]]), + bamg_geometry.Edges[pos + 1:, :] + )) + + # OK, now we can add our own rift + nods = rift[i]['nods'] - 1 + 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)) + )) + )) + bamg_geometry.Edges = np.vstack(( + bamg_geometry.Edges, + np.array([[count, count + 1, 2]]), + 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)) + )) + )) count += nods break - else: - nods = rifti['nods'] - 1 - bamg_geometry.Vertices = np.vstack(bamg_geometry.Vertices, np.hstack(rifti['x'][:], rifti['y'][:], np.ones((nods + 1, 1)))) - 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)))) - count += nods + 1 - - #Deal with tracks + nods = rift[i]['nods'] - 1 + bamg_geometry.Vertices = np.vstack(( + bamg_geometry.Vertices, + np.hstack(( + rift[i]['x'][:], + rift[i]['y'][:], + np.ones((nods + 1, 1)) + )) + )) + 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)) + )) + )) + count += (nods + 1) + + # Deal with tracks if options.exist('tracks'): - - #read tracks + # Read tracks track = options.getfieldvalue('tracks') if all(isinstance(track, str)): @@ -320,40 +434,54 @@ track = np.hstack((A.x.reshape(-1, ), A.y.reshape(-1, ))) else: - track = float(track) #for some reason, it is of class "single" + track = float(track) + if np.size(track, axis=1) == 2: track = np.hstack((track, 3. * np.ones((size(track, axis=0), 1)))) - #only keep those inside - flags = ContourToNodes(track[:, 0], track[:, 1], domainfile, 0)[0] + # Only keep those inside + 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]]' track = track[np.nonzero(flags), :] - #Add all points to bamg_geometry + # Add all points to bamg_geometry nods = np.size(track, axis=0) bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, track)) - 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)))))) - - #update counter + 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)) + )) + )) + + # Update counter count += nods - #Deal with vertices that need to be kept by mesher + # Deal with vertices that need to be kept by mesher if options.exist('RequiredVertices'): - - #recover RequiredVertices - requiredvertices = options.getfieldvalue('RequiredVertices') #for some reason, it is of class "single" + # Recover RequiredVertices + requiredvertices = options.getfieldvalue('RequiredVertices') if np.size(requiredvertices, axis=1) == 2: requiredvertices = np.hstack((requiredvertices, 4. * np.ones((np.size(requiredvertices, axis=0), 1)))) - #only keep those inside - flags = ContourToNodes(requiredvertices[:, 0], requiredvertices[:, 1], domainfile, 0)[0] + # Only keep those inside + 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]]' requiredvertices = requiredvertices[np.nonzero(flags)[0], :] - #Add all points to bamg_geometry + + # Add all points to bamg_geometry nods = np.size(requiredvertices, axis=0) bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, requiredvertices)) - #update counter + # Update counter count += nods - #process geom - #bamg_geometry = processgeometry(bamg_geometry, options.getfieldvalue('tol', float(nan)), domain[0]) + # Deal with RequiredEdges + if options.getfieldvalue('NoBoundaryRefinement', 0): + bamg_geometry.RequiredEdges = edges_required + elif options.getfieldvalue('NoBoundaryRefinementAllBoundaries', 0): + bamg_geometry.RequiredEdges = np.arange(1, bamg_geometry.Edges.shape[0]).T + + # Process geom + #bamg_geometry = processgeometry(bamg_geometry, options.getfieldvalue('tol', np.nan), domain[0]) elif isinstance(md.private.bamg, dict) and 'geometry' in md.private.bamg: bamg_geometry = bamggeom(md.private.bamg['geometry'].__dict__) @@ -362,11 +490,14 @@ pass #}}} - # Bamg Mesh parameters {{{ - if not options.exist('domain') and md.mesh.numberofvertices and m.strcmp(md.mesh.elementtype(), 'Tria'): + # Bamg mesh parameters {{{ + if not options.exist('domain') and md.mesh.numberofvertices and md.mesh.elementtype() == 'Tria': if isinstance(md.private.bamg, dict) and 'mesh' in md.private.bamg: bamg_mesh = bamgmesh(md.private.bamg['mesh'].__dict__) else: - bamg_mesh.Vertices = np.vstack((md.mesh.x, md.mesh.y, np.ones((md.mesh.numberofvertices)))).T - #bamg_mesh.Vertices = np.hstack((md.mesh.x.reshape(-1, ), md.mesh.y.reshape(-1, ), np.ones((md.mesh.numberofvertices, 1)))) + bamg_mesh.Vertices = np.vstack(( + md.mesh.x, + md.mesh.y, + np.ones((md.mesh.numberofvertices)) + )).T bamg_mesh.Triangles = np.hstack((md.mesh.elements, np.ones((md.mesh.numberofelements, 1)))) @@ -374,5 +505,5 @@ raise TypeError("bamg error message: rifts not supported yet. Do meshprocessrift AFTER bamg") #}}} - # Bamg Options {{{ + # Bamg options {{{ bamg_options['Crack'] = options.getfieldvalue('Crack', 0) bamg_options['anisomax'] = options.getfieldvalue('anisomax', 10.**18) @@ -402,8 +533,8 @@ #}}} - #call Bamg + # Call Bamg bamgmesh_out, bamggeom_out = BamgMesher(bamg_mesh.__dict__, bamg_geometry.__dict__, bamg_options) - # plug results onto model + # Plug results onto model if options.getfieldvalue('vertical', 0): md.mesh = mesh2dvertical() @@ -415,13 +546,9 @@ md.mesh.segmentmarkers = bamgmesh_out['IssmSegments'][:, 3].astype(int) - #Fill in rest of fields: + # Fill in rest of fields md.mesh.numberofelements = np.size(md.mesh.elements, axis=0) md.mesh.numberofvertices = np.size(md.mesh.x) md.mesh.numberofedges = np.size(md.mesh.edges, axis=0) - md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool) - md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True - - #Now, build the connectivity tables for this mesh. Doubled in matlab for some reason - md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, ) + md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, int) md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1 @@ -437,10 +564,10 @@ md.mesh.segmentmarkers = bamgmesh_out['IssmSegments'][:, 3].astype(int) - #Fill in rest of fields: + # Fill in rest of fields md.mesh.numberofelements = np.size(md.mesh.elements, axis=0) md.mesh.numberofvertices = np.size(md.mesh.x) md.mesh.numberofedges = np.size(md.mesh.edges, axis=0) - md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool) - md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True + md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, int) + md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1 else: @@ -453,12 +580,12 @@ md.mesh.segmentmarkers = bamgmesh_out['IssmSegments'][:, 3].astype(int) - #Fill in rest of fields: + # Fill in rest of fields md.mesh.numberofelements = np.size(md.mesh.elements, axis=0) md.mesh.numberofvertices = np.size(md.mesh.x) md.mesh.numberofedges = np.size(md.mesh.edges, axis=0) - md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool) - md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True - - #Bamg private fields + md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, int) + md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1 + + # Bamg private fields md.private.bamg = OrderedDict() md.private.bamg['mesh'] = bamgmesh(bamgmesh_out) @@ -476,13 +603,14 @@ def processgeometry(geom, tol, outline): # {{{ - raise RuntimeError("bamg.py / processgeometry is not complete.") - #Deal with edges + raise RuntimeError("bamg.py::processgeometry is not complete.") + + # Deal with edges print("Checking Edge crossing...") i = 0 while (i < np.size(geom.Edges, axis=0)): - #edge counter + # Edge counter i += 1 - #Get coordinates + # Get coordinates x1 = geom.Vertices[geom.Edges[i, 0], 0] y1 = geom.Vertices[geom.Edges[i, 0], 1] @@ -491,14 +619,14 @@ color1 = geom.Edges[i, 2] - j = i #test edges located AFTER i only + j = i # Test edges located AFTER i only while (j < np.size(geom.Edges, axis=0)): - #edge counter + # Edge counter j += 1 - #Skip if the two edges already have a vertex in common + # Skip if the two edges already have a vertex in common if any(m.ismember(geom.Edges[i, 0:2], geom.Edges[j, 0:2])): continue - #Get coordinates + # Get coordinates x3 = geom.Vertices[geom.Edges[j, 0], 0] y3 = geom.Vertices[geom.Edges[j, 0], 1] @@ -507,16 +635,16 @@ color2 = geom.Edges[j, 2] - #Check if the two edges are crossing one another + # Check if the two edges are crossing one another if SegIntersect(np.array([[x1, y1], [x2, y2]]), np.array([[x3, y3], [x4, y4]])): - #Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) + # Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) 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]]))) 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]]))) - #Add vertex to the list of vertices + # Add vertex to the list of vertices geom.Vertices = np.vstack((geom.Vertices, [x, y, min(color1, color2)])) id = np.size(geom.Vertices, axis=0) - #Update edges i and j + # Update edges i and j edgei = geom.Edges[i, :].copy() edgej = geom.Edges[j, :].copy() @@ -526,28 +654,28 @@ geom.Edges = np.vstack((geom.Edges, [id, edgej(1), edgej(2)])) - #update current edge second tip + # Update current edge second tip x2 = x y2 = y - #Check point outside + # Check point outside print("Checking for points outside the domain...") i = 0 num = 0 while (i < np.size(geom.Vertices, axis=0)): - #vertex counter + # Vertex counter i += 1 - #Get coordinates + # Get coordinates x = geom.Vertices[i, 0] y = geom.Vertices[i, 1] color = geom.Vertices[i, 2] - #Check that the point is inside the domain + # Check that the point is inside the domain if color != 1 and not ContourToNodes(x, y, outline[0], 1): - #Remove points from list of Vertices + # Remove points from list of Vertices num += 1 geom.Vertices[i, :] = [] - #update edges + # Update edges posedges = np.nonzero(geom.Edges == i) geom.Edges[posedges[0], :] = [] @@ -555,5 +683,5 @@ geom.Edges[posedges] = geom.Edges[posedges] - 1 - #update counter + # Update counter i -= 1 @@ -564,5 +692,5 @@ %Check point spacing if ~isnan(tol), - disp('Checking point spacing...') + print('Checking point spacing...') i = 0 while (i < size(geom.Vertices, 1)), @@ -608,3 +736,3 @@ """ return geom - # }}} + #}}} Index: /issm/trunk-jpl/src/m/mesh/bamgflowband.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamgflowband.py (revision 25454) +++ /issm/trunk-jpl/src/m/mesh/bamgflowband.py (revision 25455) @@ -7,6 +7,5 @@ def bamgflowband(md, x, surf, base, *args): - """ - BAMGFLOWBAND - create flowband mesh with bamg + """BAMGFLOWBAND - create flowband mesh with bamg Usage: Index: /issm/trunk-jpl/src/m/mesh/findsegments.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/findsegments.m (revision 25454) +++ /issm/trunk-jpl/src/m/mesh/findsegments.m (revision 25455) @@ -2,9 +2,9 @@ %FINDSEGMENTS - build segments model field % +% Usage: +% segments=findsegments(md,varargin); +% % Optional inputs: % 'mesh.elementconnectivity' -% -% Usage: -% segments=findsegments(md,varargin); %get options @@ -14,8 +14,8 @@ mesh.elementconnectivity=getfieldvalue(options,'mesh.elementconnectivity',md.mesh.elementconnectivity); -%Now, build the connectivity tables for this mesh if not correclty done +%Now, build the connectivity tables for this mesh if not correctly done if size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements, if exist(options,'mesh.elementconnectivity'), - error(' ''mesh.elementconnectivity'' option does not have thge right size.'); + error('''mesh.elementconnectivity'' option does not have thge right size.'); else mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); @@ -39,9 +39,9 @@ els2=mesh.elementconnectivity(el1,find(mesh.elementconnectivity(el1,:))); - %el1 is connected to 2 other elements + %get nodes of 'el1' + nods1=md.mesh.elements(el1,:); + + %'el1' is connected to 2 other elements if length(els2)>1, - - %get nodes of el1 - nods1=md.mesh.elements(el1,:); %find the common vertices to the two elements connected to el1 (1 or 2) @@ -55,4 +55,5 @@ ord1=find(nods1(1)==md.mesh.elements(el1,:)); ord2=find(nods1(2)==md.mesh.elements(el1,:)); + if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ), temp=segments(count,1); @@ -60,13 +61,13 @@ segments(count,2)=temp; end + disp(segments(count,1:2)) segments(count,1:2)=fliplr(segments(count,1:2)); + disp(segments) count=count+1; - %el1 is connected to only one element + %'el1' is connected to only one element else - %get nodes of el1 - nods1=md.mesh.elements(el1,:); - %find the vertex the el1 to not share with els2 + %find the vertex the 'el1' does not share with 'els2' flag=setdiff(nods1,md.mesh.elements(els2,:)); Index: /issm/trunk-jpl/src/m/mesh/findsegments.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/findsegments.py (revision 25455) +++ /issm/trunk-jpl/src/m/mesh/findsegments.py (revision 25455) @@ -0,0 +1,80 @@ +import numpy as np + +from ElementConnectivity import ElementConnectivity +from helpers import struct +from intersect import intersect +from pairoptions import pairoptions + +def findsegments(md, *args): #{{{ + """FINDSEGMENTS - build segments model field + + Usage: + segments = findsegments(md, args) + + Optional inputs: + 'mesh.elementconnectivity' + """ + + # Get options + options = pairoptions(*args) + + # Get connectivity + mesh = struct() + mesh.elementconnectivity = options.getfieldvalue('mesh.elementconnectivity', md.mesh.elementconnectivity) + + # Now, build the connectivity tables for this mesh if not correctly done + if md.mesh.elementconnectivity.shape[0] != md.mesh.numberofelements: + if options.exist('mesh.elementconnectivity'): + raise Exception('\'mesh.elementconnectivity\' option does not have the right size.') + else: + mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity) + + # Recreate the segments + elementonboundary = np.zeros((md.mesh.numberofelements, )) + elementonboundary[np.where(mesh.elementconnectivity[:, 2] == 0)[0]] = 1 + pos = np.nonzero(elementonboundary)[0] + num_segments = len(pos) + segments = np.zeros((num_segments, 3)) + count = 0 + + # Loop over the segments + for i in range(num_segments): + # Get current element on boundary + el1 = pos[i] + + # Get elements connected to 'el1' + els2 = mesh.elementconnectivity[el1, np.nonzero(mesh.elementconnectivity[el1, :])[0]] + + # Get nodes of 'el1' + nods1 = md.mesh.elements[el1, :] + + # 'el1' is connected to 2 other elements + if len(els2) > 1: + + # Find the common vertices to the two elements connected to 'el1' (1 or 2) + flag = intersect(md.mesh.elements[els2[0] - 1, :], md.mesh.elements[els2[1] - 1, :])[0] # NOTE: Throwing away second- and third- position values returned from call + + # Get the vertices on the boundary and build segment + nods1 = np.delete(nods1, np.where(nods1 == flag)[0]) + segments[count, :] = np.append(nods1, el1 + 1) + + # Swap segment nodes if necessary + ord1 = np.where(nods1[0] == md.mesh.elements[el1, :])[0][0] + ord2 = np.where(nods1[1] == md.mesh.elements[el1, :])[0][0] + + if ((ord1 == 0 and ord2 == 1) or (ord1 == 1 and ord2 == 2) or (ord1 == 2 and ord2 == 0)): + temp = segments[count, 0] + segments[count, 0] = segments[count, 1] + segments[count, 1] = temp + + segments[count, 0:2] = np.flip(segments[count, 0:2]) # NOTE: Upper bound of index range is non-inclusive + count = count + 1 + # 'el1' is connected to only one element + else: + # Find the vertex that 'el1' does not share with 'els2' + flag = els2 + + exit() + + return segments +#}}} Index: /issm/trunk-jpl/src/m/mesh/meshconvert.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/meshconvert.py (revision 25454) +++ /issm/trunk-jpl/src/m/mesh/meshconvert.py (revision 25455) @@ -1,23 +1,25 @@ import numpy as np + +from BamgConvertMesh import BamgConvertMesh +from bamggeom import bamggeom +from bamgmesh import bamgmesh from collections import OrderedDict -from BamgConvertMesh import BamgConvertMesh from mesh2d import mesh2d -from bamgmesh import bamgmesh -from bamggeom import bamggeom def meshconvert(md, *args): - """ - CONVERTMESH - convert mesh to bamg mesh + """CONVERTMESH - convert mesh to bamg mesh - Usage: - md = meshconvert(md) - md = meshconvert(md, index, x, y) + Usage: + md = meshconvert(md) + md = meshconvert(md, index, x, y) """ - if not len(args) == 0 and not len(args) == 3: + nargs = len(args) + + if not nargs == 0 and not nargs == 3: raise TypeError("meshconvert error message: bad usage") - if not len(args): + if not nargs: index = md.mesh.elements x = md.mesh.x @@ -47,6 +49,6 @@ md.mesh.numberofvertices = np.size(md.mesh.x) md.mesh.numberofedges = np.size(md.mesh.edges, axis=0) - md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool) - md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True + md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, int) + md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1 return md Index: /issm/trunk-jpl/src/m/mesh/modelmerge3d.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/modelmerge3d.m (revision 25454) +++ /issm/trunk-jpl/src/m/mesh/modelmerge3d.m (revision 25455) @@ -8,5 +8,5 @@ options=pairoptions(varargin{:}); - tolerance=getfieldvalue(options,'tolerance',10^-5); + tolerance=getfieldvalue(options,'tolerance',1e-4); md=model(); @@ -16,5 +16,5 @@ md.private=md1.private; - %some initializatoin: + %some initialization: elements1=md1.mesh.elements; x1=md1.mesh.x; @@ -34,10 +34,15 @@ elements2=elements2+nods1; - %go into the vertices on boundary of mesh 1, and figure out which ones are common to mesh2: - verticesonboundary=find(md1.mesh.vertexonboundary); + %go into the vertices on boundary of mesh 1 and figure out which ones are common with mesh2: + verticesonboundary=find(md1.mesh.vertexonboundary); + for i=1:length(verticesonboundary), - node1=verticesonboundary(i); xnode1=x1(node1); ynode1=y1(node1); znode1=z1(node1); - %is there another node with these coordinates in mesh2? - ind=find(sqrt(((x2-xnode1).^2+(y2-ynode1).^2)+(z2-znode1).^2)1, disp('should reduce the tolerance, several vertices picked up!'); @@ -47,20 +52,21 @@ y2(ind)=NaN; z2(ind)=NaN; - pos=find(elements2==(ind+nods1)); elements2(pos)=node1; + pos=find(elements2==(ind+nods1)); + elements2(pos)=node1; end end + %go through elements2 and drop counter on each vertex that is above the x2 and y2 vertices being dropped: + indices_nan=isnan(x2); + while(~isempty(indices_nan)), + % Use the index of the first instance of 'nan' value to remove that element from 'x2', 'y2', and 'z2' + index_nan=indices_nan(1); + pos=find(elements2>(index_nan+nods1)); + elements2(pos)=elements2(pos)-1; + x2(index_nan)=[]; + y2(index_nan)=[]; + z2(index_nan)=[]; - %go through elements2 and drop counter on each vertex that is above the x2 and y2 vertices being dropped: - while( ~isempty(find(isnan(x2)))), - for i=1:length(x2), - if isnan(x2(i)), - pos=find(elements2>(i+nods1)); - elements2(pos)=elements2(pos)-1; - x2(i)=[]; - y2(i)=[]; - z2(i)=[]; - break; - end - end + % Check again in 'x2' for instances of 'nan' + indices_nan=find(isnan(x2)); end @@ -83,5 +89,7 @@ %connectivities: md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); + % disp(md.mesh.vertexconnectivity) md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); + % disp(md.mesh.elementconnectivity) %find segments: @@ -91,4 +99,6 @@ md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; + disp(md.mesh.vertexonboundary) + pause %some checks: Index: /issm/trunk-jpl/src/m/mesh/modelmerge3d.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/modelmerge3d.py (revision 25455) +++ /issm/trunk-jpl/src/m/mesh/modelmerge3d.py (revision 25455) @@ -0,0 +1,124 @@ +from copy import copy, deepcopy + +import numpy as np + +from ElementConnectivity import ElementConnectivity +from findsegments import findsegments +from model import model +from NodeConnectivity import NodeConnectivity +from pairoptions import pairoptions + + +def modelmerge3d(md1, md2, *args): + """MODELMERGE - Merge two models by merging their meshes. + + Usage: + md = modelmerge(md1, md2) + """ + + # Process options + options = pairoptions(*args) + + tolerance = options.getfieldvalue('tolerance', 1e-4) + + md = model() + + # First, copy md1 mesh into md.mesh to initialize, and additional classes + md.mesh = deepcopy(md1.mesh) + md.private = deepcopy(md1.private) + + # Some initialization + elements1 = copy(md1.mesh.elements) + x1 = copy(md1.mesh.x) + y1 = copy(md1.mesh.y) + z1 = copy(md1.mesh.z) + nods1 = copy(md1.mesh.numberofvertices) + nel1 = copy(md1.mesh.numberofelements) + + elements2 = copy(md2.mesh.elements) + x2 = copy(md2.mesh.x) + y2 = copy(md2.mesh.y) + z2 = copy(md2.mesh.z) + nods2 = copy(md2.mesh.numberofvertices) + nel2 = copy(md2.mesh.numberofelements) + + # Offset elements2 by nods1 + elements2 = elements2 + nods1 + + # Go into the vertices on boundary of mesh 1 and figure out which ones are common with mesh 2 + verticesonboundary = np.nonzero(md1.mesh.vertexonboundary)[0] + + # Do not display "RuntimeWarning: invalid value encountered in less" warning when comparing against np.nan values in 'x2', 'y2', and 'z2' + # + # TODO: Investigate if we should use some other value to represent null elements + np.warnings.filterwarnings('ignore') + + for i in range(len(verticesonboundary)): + node1 = verticesonboundary[i] + xnode1 = x1[node1] + ynode1 = y1[node1] + znode1 = z1[node1] + + # Is there another node with these coordinates in mesh 2? + ind = np.where(np.logical_and.reduce((~np.isnan(x2), ~np.isnan(y2), ~np.isnan(z2), ((x2 - xnode1) ** 2 + (y2 - ynode1) ** 2 + (z2 - znode1) ** 2) ** 0.5 < tolerance)))[0] # NOTE: math.sqrt cannot be applied element-wise to a list/numpy.array + if len(ind) > 1: + print('should reduce the tolerance, several vertices picked up!') + if len(ind): + x2[ind] = np.nan + y2[ind] = np.nan + z2[ind] = np.nan + pos = np.where(elements2 == ((ind + 1) + nods1)) + elements2[pos] = (node1 + 1) # NOTE: 'elements2' is 2D array, so 'numpy.where' returns two lists of indices + + # Go through elements2 and drop counter on each vertex that is above the x2 and y2 vertices being dropped + indices_nan = np.nonzero(np.isnan(x2).astype(int))[0] + while indices_nan.size: + # Use the index of the first instance of 'nan' value to remove that element from 'x2', 'y2', and 'z2' + index_nan = indices_nan[0] + pos = np.where(elements2 > ((index_nan + 1) + nods1)) + elements2[pos] = elements2[pos] - 1 + # TODO: Maybe set to None, then pop all None after this loop + x2 = np.delete(x2, index_nan) + y2 = np.delete(y2, index_nan) + z2 = np.delete(z2, index_nan) + + # Check again in 'x2' for instances of 'nan' + indices_nan = np.nonzero(np.isnan(x2).astype(int))[0] + + # Merge elements + elements = np.concatenate((elements1, elements2)) + + # Merge vertices + x = np.concatenate((x1, x2)) + y = np.concatenate((y1, y2)) + z = np.concatenate((z1, z2)) + + # Output + md.mesh.x = x + md.mesh.y = y + md.mesh.z = z + md.mesh.elements = elements + md.mesh.numberofvertices = len(x) + md.mesh.numberofelements = elements.shape[0] + + # Connectivities + md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices) + # print(md.mesh.vertexconnectivity) + md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity) + # print(md.mesh.elementconnectivity) + + # Find segments + md.mesh.segments = findsegments(md) + + # Vertex on boundary + md.mesh.vertexconnectivity = np.zeros((md.mesh.numberofvertices, )) + md.mesh.vertexonboundary[md.mesh.segments[:, 0:1]] = 1 + print(md.mesh.vertexonboundary) + + # Some checks + if np.max(md.mesh.elements) > md.mesh.numberofvertices: + raise Exception('issue in modelmerge, one of the element ids is > number of vertices!') + + exit() + + Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py (revision 25454) +++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py (revision 25455) @@ -1,15 +1,14 @@ import numpy as np + +from ContourToMesh import ContourToMesh from ElementsFromEdge import ElementsFromEdge import MatlabFuncs as m -from ContourToMesh import ContourToMesh def meshprocessoutsiderifts(md, domainoutline): - """ - MESHPROCESSOUTSIDERIFTS - process rifts when they touch the domain outline + """MESHPROCESSOUTSIDERIFTS - process rifts when they touch the domain outline - Usage: - md = meshprocessoutsiderifts(md, domain) - + Usage: + md = meshprocessoutsiderifts(md, domain) """ @@ -53,5 +52,5 @@ elements = np.concatenate((elements, nextelement)) #new B: - B = md.mesh.elements[nextelement - 1, np.nonzero(np.logical_not(m.ismember(md.mesh.elements[nextelement - 1, :], np.array([A, B]))))] + B = md.mesh.elements[nextelement - 1, np.nonzero(np.logical_not(m.ismember(md.mesh.elements[nextelement - 1, :], np.array([A, B]))))[0]] #take the list of elements on one side of the rift that connect to the tip, @@ -64,5 +63,5 @@ #replace tip in elements newelements = md.mesh.elements[elements - 1, :] - pos = np.nonzero(newelements == tip) + pos = np.nonzero(newelements == tip)[0] newelements[pos] = num md.mesh.elements[elements - 1, :] = newelements @@ -81,6 +80,6 @@ md.mesh.numberofelements = np.size(md.mesh.elements, axis=0) md.mesh.numberofvertices = np.size(md.mesh.x) - md.mesh.vertexonboundary = np.zeros(np.size(md.mesh.x), bool) - md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True + md.mesh.vertexonboundary = np.zeros(np.size(md.mesh.x), int) + md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1 md.rifts.numrifts = np.length(md.rifts.riftstruct) @@ -88,10 +87,9 @@ -def isconnected(elements, A, B): # {{{ - """ - ISCONNECTED: are two nodes connected by a triangulation? +def isconnected(elements, A, B): #{{{ + """ISCONNECTED: are two nodes connected by a triangulation? - Usage: flag = isconnected(elements, A, B) - + Usage: + flag = isconnected(elements, A, B) """ @@ -103,3 +101,3 @@ return flag - # }}} + #}}} Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py (revision 25454) +++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py (revision 25455) @@ -1,26 +1,25 @@ import numpy as np + +from ContourToMesh import ContourToMesh +from GetAreas import GetAreas +from meshprocessoutsiderifts import meshprocessoutsiderifts from ProcessRifts import ProcessRifts -from ContourToMesh import ContourToMesh -from meshprocessoutsiderifts import meshprocessoutsiderifts -from GetAreas import GetAreas def meshprocessrifts(md, domainoutline): - """ - MESHPROCESSRIFTS - process mesh when rifts are present + """MESHPROCESSRIFTS - process mesh when rifts are present - split rifts inside mesh (rifts are defined by presence of - segments inside the domain outline) - if domain outline is provided, check for rifts that could touch it, and open them up. + split rifts inside mesh (rifts are defined by presence of + segments inside the domain outline) + if domain outline is provided, check for rifts that could touch it, and open them up. - Usage: - md = meshprocessrifts(md, domainoutline) + Usage: + md = meshprocessrifts(md, domainoutline) - Ex: - md = meshprocessrifts(md, 'DomainOutline.exp') - + Example: + md = meshprocessrifts(md, 'DomainOutline.exp') """ - #Call MEX file + # Call Python module 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) md.mesh.elements = md.mesh.elements.astype(int) @@ -35,6 +34,6 @@ md.mesh.numberofelements = np.size(md.mesh.elements, axis=0) md.mesh.numberofvertices = np.size(md.mesh.x) - md.mesh.vertexonboundary = np.zeros(np.size(md.mesh.x), bool) - md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True + md.mesh.vertexonboundary = np.zeros(np.size(md.mesh.x), int) + md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1 #get coordinates of rift tips Index: /issm/trunk-jpl/src/m/mesh/squaremesh.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/squaremesh.py (revision 25454) +++ /issm/trunk-jpl/src/m/mesh/squaremesh.py (revision 25455) @@ -1,19 +1,19 @@ import numpy as np -from NodeConnectivity import NodeConnectivity + from ElementConnectivity import ElementConnectivity from mesh2d import mesh2d +from NodeConnectivity import NodeConnectivity def squaremesh(md, Lx, Ly, nx, ny): - """ - SQUAREMESH - create a structured square mesh + """SQUAREMESH - create a structured square mesh - This script will generate a structured square mesh - Lx and Ly are the dimension of the domain (in meters) - nx anx ny are the number of nodes in the x and y direction - The coordinates x and y returned are in meters. + This script will generate a structured square mesh + Lx and Ly are the dimension of the domain (in meters) + nx anx ny are the number of nodes in the x and y direction + The coordinates x and y returned are in meters. - Usage: - [md] = squaremesh(md, Lx, Ly, nx, ny) + Usage: + [md] = squaremesh(md, Lx, Ly, nx, ny) """ @@ -24,6 +24,6 @@ #initialization index = np.zeros((nel, 3), int) - x = np.zeros((nx * ny)) - y = np.zeros((nx * ny)) + x = np.zeros((nx * ny, )) + y = np.zeros((nx * ny, )) #create coordinates @@ -63,6 +63,6 @@ md.mesh.y = y md.mesh.numberofvertices = nods - md.mesh.vertexonboundary = np.zeros((nods), bool) - md.mesh.vertexonboundary[segments[:, 0:2] - 1] = True + md.mesh.vertexonboundary = np.zeros(nods, int) + md.mesh.vertexonboundary[segments[:, 0:2] - 1] = 1 #plug elements @@ -72,6 +72,6 @@ #Now, build the connectivity tables for this mesh. - md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)[0] - md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)[0] + md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices) + md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity) return md Index: /issm/trunk-jpl/src/m/mesh/triangle.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/triangle.py (revision 25454) +++ /issm/trunk-jpl/src/m/mesh/triangle.py (revision 25455) @@ -1,26 +1,28 @@ import os.path + import numpy as np + +from ElementConnectivity import ElementConnectivity from mesh2d import mesh2d from NodeConnectivity import NodeConnectivity -from ElementConnectivity import ElementConnectivity from Triangle_python import Triangle_python def triangle(md, domainname, *args): - """ - TRIANGLE - create model mesh using the triangle package + """TRIANGLE - create model mesh using the triangle package - This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution - where md is a @model object, domainname is the name of an Argus domain outline file, - and resolution is a characteristic length for the mesh (same unit as the domain outline - unit). Riftname is an optional argument (Argus domain outline) describing rifts. + This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution + where md is a @model object, domainname is the name of an Argus domain outline file, + and resolution is a characteristic length for the mesh (same unit as the domain outline + unit). Riftname is an optional argument (Argus domain outline) describing rifts. - Usage: - md = triangle(md, domainname, resolution) - or md = triangle(md, domainname, resolution, riftname) + Usage: + md = triangle(md, domainname, resolution) + OR + md = triangle(md, domainname, resolution, riftname) - Examples: - md = triangle(md, 'DomainOutline.exp', 1000) - md = triangle(md, 'DomainOutline.exp', 1000, 'Rifts.exp') + Examples: + md = triangle(md, 'DomainOutline.exp', 1000) + md = triangle(md, 'DomainOutline.exp', 1000, 'Rifts.exp') """ @@ -45,5 +47,5 @@ return None - area = resolution**2 + area = resolution ** 2 #Check that file exist (this is a very very common mistake) @@ -61,10 +63,10 @@ md.mesh.numberofelements = np.size(md.mesh.elements, axis=0) md.mesh.numberofvertices = np.size(md.mesh.x) - md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool) - md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True + md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, int) + md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1 #Now, build the connectivity tables for this mesh. - md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)[0] - md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)[0] + md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices) + md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity) return md Index: /issm/trunk-jpl/src/m/miscellaneous/intersect.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/intersect.py (revision 25455) +++ /issm/trunk-jpl/src/m/miscellaneous/intersect.py (revision 25455) @@ -0,0 +1,25 @@ +import numpy as np + + +def intersect(a, b): #{{{ + """INTERSECT - Python implementation of MATLAB's 'intersect' function + + Usage: + c = intersect(a, b) => Returns a list of values common to both 'a' and 'b', with no repetitions. 'c' is in sorted order. + + c, ia, ib = intersect(a, b) => Also returns index lists 'ia' and 'ib'. + + Sources: + - https://www.mathworks.com/help/matlab/ref/double.intersect.html + - https://stackoverflow.com/a/45645177 + """ + a_unique, ia = np.unique(a, return_index=True) + b_unique, ib = np.unique(b, return_index=True) + + all_unique = np.concatenate((a_unique, b_unique)) + all_unique.sort() + + c = all_unique[:-1][all_unique[1:] == all_unique[:-1]] + + return c, ia[np.isin(a_unique, c)], ib[np.isin(b_unique, c)] +#}}} Index: /issm/trunk-jpl/src/m/miscellaneous/pretty_print.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/pretty_print.m (revision 25455) +++ /issm/trunk-jpl/src/m/miscellaneous/pretty_print.m (revision 25455) @@ -0,0 +1,77 @@ +function pretty_print(data) +%PRETTY_PRINT - print longer structures as they would be under Python +% +% Utility function for debugging to print large structures. +% +% Usage: +% pretty_print(data) +% +% NOTE: +% - Currently only handles 1- and 2D structures. +% +% TODO: +% - Add an argument that allows the user to specify what constitutes a "long" +% large data structure (default is length of 6). +% - Add an argument that allows the user to specify the number of values that +% they would like to display from the head and tail of each dimension of +% 'data'. + +if ndims(data)==1 + if length(data)>6 + disp(sprintf('[%d %d %d ... %d %d %d]',data(1),data(2),data(3),data(end-2),data(end-1),data(end))); + else + disp(data); + end +elseif ndims(data)==2 + shape=size(data); + if shape(1)>6 + % if shape(2)==1 + % disp('NOTE: Single column printed as row!'); + % disp(sprintf('[%d %d %d ... %d %d %d]',data(1,1),data(2,1),data(3,1),data(end-2,1),data(end-1,1),data(end,1))); + % else + if shape(2)>6 + disp('['); + disp(sprintf('[%d %d %d ... %d %d %d]',data(1,1),data(1,2),data(1,3),data(1,end-2),data(1,end-1),data(1,end))); + disp(sprintf('[%d %d %d ... %d %d %d]',data(2,1),data(2,2),data(2,3),data(2,end-2),data(2,end-1),data(2,end))); + disp(sprintf('[%d %d %d ... %d %d %d]',data(3,1),data(3,2),data(3,3),data(3,end-2),data(3,end-1),data(3,end))); + disp('...'); + disp(sprintf('[%d %d %d ... %d %d %d]',data(end-2,1),data(end-2,2),data(end-2,3),data(end-2,end-2),data(end-2,end-1),data(end-2,end))); + disp(sprintf('[%d %d %d ... %d %d %d]',data(end-1,1),data(end-1,2),data(end-1,3),data(end-1,end-2),data(end-1,end-1),data(end-1,end))); + disp(sprintf('[%d %d %d ... %d %d %d]',data(end,1),data(end,2),data(end,3),data(end,end-2),data(end,end-1),data(end,end))); + disp(sprintf(']\n')); + else + disp('['); + disp(data(1,:)); + disp(data(2,:)); + disp(data(3,:)); + disp('...'); + disp(data(end-2,:)); + disp(data(end-1,:)); + disp(data(end,:)); + disp(sprintf(']\n')); + end + else + % if shape(2)==1 + % data=transpose(data) + % sprintf('% NOTE: Single column printed as row!'); + % % Get shape of transposed structure + % shape=size(data); + % if shape(2)>6 + % disp(sprintf('[%d %d %d ... %d %d %d]',data(1,1),data(2,1),data(3,1),data(end-2,1),data(end-1,1),data(end,1))); + % else + % disp(data); + % end + % else + if shape(2)>6 + disp('['); + for i=1:shape(1) + disp(sprintf('[%d %d %d ... %d %d %d]',data(i,1),data(i,2),data(i,3),data(i,end-2),data(i,end-1),data(i,end))); + end + disp(sprintf(']\n')); + else + disp(data); + end + end +else + disp(data); +end Index: /issm/trunk-jpl/src/m/modules/BamgMesher.py =================================================================== --- /issm/trunk-jpl/src/m/modules/BamgMesher.py (revision 25454) +++ /issm/trunk-jpl/src/m/modules/BamgMesher.py (revision 25455) @@ -3,18 +3,16 @@ def BamgMesher(bamgmesh, bamggeom, bamgoptions): - """ - BAMGMESHER - - Usage: - bamgmesh, bamggeom = BamgMesher(bamgmesh, bamggeom, bamgoptions) + """BAMGMESHER bamgmesh: input bamg mesh bamggeom: input bamg geometry for the mesh bamgoptions: options for the bamg mesh + + Usage: + bamgmesh, bamggeom = BamgMesher(bamgmesh, bamggeom, bamgoptions) """ - #Call mex module + # Call module bamgmesh, bamggeom = BamgMesher_python(bamgmesh, bamggeom, bamgoptions) - #return return bamgmesh, bamggeom Index: /issm/trunk-jpl/src/m/modules/ContourToNodes.py =================================================================== --- /issm/trunk-jpl/src/m/modules/ContourToNodes.py (revision 25454) +++ /issm/trunk-jpl/src/m/modules/ContourToNodes.py (revision 25455) @@ -3,19 +3,17 @@ def ContourToNodes(x, y, contourname, edgevalue): - """ - CONTOURTONODES - flags vertices inside contour + """CONTOURTONODES - flags vertices inside contour - Usage: - flags = ContourToNodes(x, y, contourname, edgevalue) + x, y: list of nodes + contourname: name of .exp/.shp file containing the contours, or resulting structure from call to expread/shpread + edgevalue: integer (0, 1 or 2) defining the value associated to the nodes on the edges of the polygons + flags: vector of flags (0 or 1), of size nodes - x, y: list of nodes - contourname: name of .exp file containing the contours, or resulting structure from call to expread - edgevalue: integer (0, 1 or 2) defining the value associated to the nodes on the edges of the polygons - flags: vector of flags (0 or 1), of size nodes + Usage: + flags = ContourToNodes(x, y, contourname, edgevalue) """ - #Call mex module + # Call Python module flags = ContourToNodes_python(x, y, contourname, edgevalue) - #return return flags Index: /issm/trunk-jpl/src/m/modules/ElementConnectivity.py =================================================================== --- /issm/trunk-jpl/src/m/modules/ElementConnectivity.py (revision 25454) +++ /issm/trunk-jpl/src/m/modules/ElementConnectivity.py (revision 25455) @@ -3,13 +3,12 @@ def ElementConnectivity(elements, nodeconnectivity): + """ELEMENTCONNECTIVITY - Build element connectivity using node connectivity and elements + + Usage: + elementconnectivity = ElementConnectivity(elements, nodeconnectivity) """ - ELEMENTCONNECTIVITY - Build element connectivity using node connectivity and elements - - Usage: - elementconnectivity = ElementConnectivity(elements, nodeconnectivity) - """ - #Call mex module + + # Call Python module elementconnectivity = ElementConnectivity_python(elements, nodeconnectivity) - #Return - return elementconnectivity + return elementconnectivity[0] # NOTE: Value returned from wrapper function is a tuple, the first element of which being the result we actually want Index: /issm/trunk-jpl/src/m/modules/ExpToLevelSet.py =================================================================== --- /issm/trunk-jpl/src/m/modules/ExpToLevelSet.py (revision 25454) +++ /issm/trunk-jpl/src/m/modules/ExpToLevelSet.py (revision 25455) @@ -2,24 +2,24 @@ from shpread import shpread + def ExpToLevelSet(x, y, contourname): #{{{ - ''' - EXPTOLEVELSET - Determine levelset distance between a contour and a cloud - of points + """EXPTOLEVELSET - Determine levelset distance between a contour and a + cloud of points - Usage: - distance = ExpToLevelSet(x, y, contourname) + Usage: + distance = ExpToLevelSet(x, y, contourname) - x, y: cloud point - contourname: name of .exp file containing the contours - distance: distance vector representing a levelset where the 0 - level is one of the contour segments + x, y: cloud point + contourname: name of .exp file containing the contours + distance: distance vector representing a levelset where the 0 + level is one of the contour segments - Example: - distance = ExpToLevelSet(md.mesh.x, md.mesh.y, 'Contour.exp') + Example: + distance = ExpToLevelSet(md.mesh.x, md.mesh.y, 'Contour.exp') - TODO: - - Need to compile Python version of ExpToLevelSet_matlab for this - to work as intended (see src/m/modules/ExpToLevelSet.m) - ''' + TODO: + - Need to compile Python version of ExpToLevelSet_matlab for this + to work as intended (see src/m/modules/ExpToLevelSet.m) + """ if isinstance(contourname, basestring): Index: /issm/trunk-jpl/src/m/modules/InterpFromMesh2d.m =================================================================== --- /issm/trunk-jpl/src/m/modules/InterpFromMesh2d.m (revision 25454) +++ /issm/trunk-jpl/src/m/modules/InterpFromMesh2d.m (revision 25455) @@ -4,9 +4,11 @@ % Usage: % data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime); -% or data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value); -% or data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value,contourname); +% OR +% data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value); +% OR +% data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value,contourname); % +% index: index of the mesh where data is defined % x,y: coordinates of the nodes where data is defined -% index: index of the mesh where data is defined % data: vector holding the data to be interpolated onto the points % x_prime,y_prime: coordinates of the mesh vertices onto which we interpolate @@ -32,4 +34,5 @@ data_prime=InterpFromMesh2d_matlab(varargin{1},varargin{2},varargin{3},varargin{4},varargin{5},varargin{6},varargin{7},varargin{8}); otherwise + % NOTE: Should never get here because of previous check error('InterpFromMesh2d not supported'); end Index: /issm/trunk-jpl/src/m/modules/InterpFromMesh2d.py =================================================================== --- /issm/trunk-jpl/src/m/modules/InterpFromMesh2d.py (revision 25455) +++ /issm/trunk-jpl/src/m/modules/InterpFromMesh2d.py (revision 25455) @@ -0,0 +1,41 @@ +from InterpFromMesh2d_python import InterpFromMesh2d_python + + +def InterpFromMesh2d(*args): #{{{ + """INTERPFROMMESH2D + + Usage: + data_prime = InterpFromMesh2d(index, x, y, data, x_prime, y_prime) + OR + data_prime = InterpFromMesh2d(index, x, y, data, x_prime, y_prime, default_value) + OR + data_prime = InterpFromMesh2d(index, x, y, data, x_prime, y_prime, defualt_value, contourname) + + index: index of the mesh where data is defined + x, y: coordinates of the nodes where data is defined + data: vector holding the data to be intepolated onto the points + x_prime, y_prime: coordinates of the mesh vertices onto which we interpolate + default_value: a scalar or vector of size len(x_prime) + contourname: linear interpolation will happen on all x_interp, y_interp inside the contour, default vlaue will be adopted on the rest of the mesh + data_prime: vector of prime interpolated data + """ + + # Check usage + nargs = len(args) + if nargs != 6 and nargs != 7 and nargs != 8: + print(InterpFromMesh2d.__doc__) + raise Exception('Wrong usage (see above)') + + # Call Python module + if nargs == 6: + data_prime = InterpFromMesh2d_python(args[0], args[1], args[2], args[3], args[4], args[5]) + elif nargs == 7: + data_prime = InterpFromMesh2d_python(args[0], args[1], args[2], args[3], args[4], args[5], args[6]) + elif nargs == 8: + data_prime = InterpFromMesh2d_python(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7]) + else: + # NOTE: Should never get here because of previous check + raise Exception('InterpFromMesh2d not supported') + + return data_prime[0] # NOTE: Value returned from wrapper function is a tuple, the first element of which being the result we actually want +#}}} Index: /issm/trunk-jpl/src/m/modules/InterpFromMeshToMesh2d.py =================================================================== --- /issm/trunk-jpl/src/m/modules/InterpFromMeshToMesh2d.py (revision 25454) +++ /issm/trunk-jpl/src/m/modules/InterpFromMeshToMesh2d.py (revision 25455) @@ -3,26 +3,33 @@ def InterpFromMeshToMesh2d(*args): - """ - INTERPFROMMESHTOMESH2D - Interpolation from a 2d triangular mesh onto a list of points + """INTERPFROMMESHTOMESH2D - Interpolation from a 2d triangular mesh onto a list of points - Usage: - data_interp = InterpFromMeshToMesh2d(index, x, y, data, x_interp, y_interp) - or data_interp = InterpFromMeshToMesh2d(index, x, y, data, x_interp, y_interp, OPTIONS) + Usage: + data_interp = InterpFromMeshToMesh2d(index, x, y, data, x_interp, y_interp) + OR + data_interp = InterpFromMeshToMesh2d(index, x, y, data, x_interp, y_interp, OPTIONS) - index: index of the mesh where data is defined (e.g. md.mesh.elements) - x, y: coordinates of the nodes where data is defined - data: matrix holding the data to be interpolated onto the mesh (one column per field) - x_interp, y_interp: coordinates of the points onto which we interpolate - data_interp: vector of mesh interpolated data - Available options: - default: default value if point is outsite of triangulation (instead of linear interpolation) + index: index of the mesh where data is defined (e.g. md.mesh.elements) + x, y: coordinates of the nodes where data is defined + data: matrix holding the data to be interpolated onto the mesh (one column per field) + x_interp, y_interp: coordinates of the points onto which we interpolate + data_interp: vector of mesh interpolated data + Available options: + default: default value if point is outsite of triangulation (instead of linear interpolation) Example: - load('temperature.mat') - md.initialization.temperature = InterpFromMeshToMesh2d(index, x, y, temperature, md.mesh.x, md.mesh.y) - md.initialization.temperature = InterpFromMeshToMesh2d(index, x, y, temperature, md.mesh.x, md.mesh.y, 'default', 253) + load('temperature.mat') + md.initialization.temperature = InterpFromMeshToMesh2d(index, x, y, temperature, md.mesh.x, md.mesh.y) + md.initialization.temperature = InterpFromMeshToMesh2d(index, x, y, temperature, md.mesh.x, md.mesh.y, 'default', 253) """ - # Call mex module + + # Check usage + nargs = len(args) + if nargs != 6 and nargs != 8: + print(InterpFromMeshToMesh2d.__doc__) + raise Exception('Wrong usage (see above)') + + # Call Python module data_interp = InterpFromMeshToMesh2d_python(*args) - # Return - return data_interp + + return data_interp[0] # NOTE: Value returned from wrapper function is a tuple, the first element of which being the result we actually want Index: /issm/trunk-jpl/src/m/modules/InterpFromMeshToMesh3d.py =================================================================== --- /issm/trunk-jpl/src/m/modules/InterpFromMeshToMesh3d.py (revision 25454) +++ /issm/trunk-jpl/src/m/modules/InterpFromMeshToMesh3d.py (revision 25455) @@ -3,22 +3,27 @@ def InterpFromMeshToMesh3d(index, x, y, z, data, x_prime, y_prime, z_prime, default_value): + """INTERPFROMMESHTOMESH3D - Interpolation from a 3d hexahedron mesh onto a list of points + + Usage: + index: index of the mesh where data is defined + x, y, z: coordinates of the nodes where data is defined + data: matrix holding the data to be interpolated onto the mesh + x_prime, y_prime, z_prime: coordinates of the points onto which we interpolate + default_value: default value if no data is found (holes) + data_prime: vector of mesh interpolated data + + Example: + load('temperature.mat') + md.initialization.temperature = InterpFromMeshToMesh3d(index, x, y, z, temperature, md.mesh.x, md.mesh.y, md.mesh.z, 253) """ - INTERPFROMMESHTOMESH3D - Interpolation from a 3d hexahedron mesh onto a list of points - Usage: - index: index of the mesh where data is defined - x, y, z: coordinates of the nodes where data is defined - data: matrix holding the data to be interpolated onto the mesh - x_prime, y_prime, z_prime: coordinates of the points onto which we interpolate - default_value: default value if no data is found (holes) - data_prime: vector of mesh interpolated data + # Check usage + nargs = len(args) + if nargs != 9: + print(InterpFromMeshToMesh3d.__doc__) + raise Exception('Wrong usage (see above)') - Example: - load('temperature.mat') - md.initialization.temperature = InterpFromMeshToMesh3d(index, x, y, z, temperature, md.mesh.x, md.mesh.y, md.mesh.z, 253) - """ - # Call mex module + # Call Python module data_prime = InterpFromMeshToMesh3d_python(index, x, y, z, data, x_prime, y_prime, z_prime, default_value) - # Return return data_prime Index: /issm/trunk-jpl/src/m/modules/NodeConnectivity.py =================================================================== --- /issm/trunk-jpl/src/m/modules/NodeConnectivity.py (revision 25454) +++ /issm/trunk-jpl/src/m/modules/NodeConnectivity.py (revision 25455) @@ -3,13 +3,12 @@ def NodeConnectivity(elements, numnodes): + """NODECONNECTIVITY - Build node connectivity from elements + + Usage: + connectivity = NodeConnectivity(elements, numnodes) """ - NODECONNECTIVITY - Build node connectivity from elements - Usage: - connectivity = NodeConnectivity(elements, numnodes) - """ - # Call mex module + # Call Python module connectivity = NodeConnectivity_python(elements, numnodes) - # Return - return connectivity + return connectivity[0] # NOTE: Value returned from wrapper function is a tuple, the first element of which being the result we actually want Index: /issm/trunk-jpl/src/m/parameterization/contourenvelope.py =================================================================== --- /issm/trunk-jpl/src/m/parameterization/contourenvelope.py (revision 25454) +++ /issm/trunk-jpl/src/m/parameterization/contourenvelope.py (revision 25455) @@ -9,13 +9,12 @@ def contourenvelope(md, *args): - """ - CONTOURENVELOPE - build a set of segments enveloping a contour .exp + """CONTOURENVELOPE - build a set of segments enveloping a contour .exp - Usage: - segments = contourenvelope(md, varargin) + Usage: + segments = contourenvelope(md, varargin) - Example: - segments = contourenvelope(md, 'Stream.exp') - segments = contourenvelope(md) + Example: + segments = contourenvelope(md, 'Stream.exp') + segments = contourenvelope(md) """ @@ -41,7 +40,7 @@ #Computing connectivity if np.size(md.mesh.vertexconnectivity, axis=0) != md.mesh.numberofvertices and np.size(md.mesh.vertexconnectivity, axis=0) != md.mesh.numberofvertices2d: - md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)[0] + md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices) if np.size(md.mesh.elementconnectivity, axis=0) != md.mesh.numberofelements and np.size(md.mesh.elementconnectivity, axis=0) != md.mesh.numberofelements2d: - md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)[0] + md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity) #get nodes inside profile Index: /issm/trunk-jpl/src/m/partition/adjacency.py =================================================================== --- /issm/trunk-jpl/src/m/partition/adjacency.py (revision 25454) +++ /issm/trunk-jpl/src/m/partition/adjacency.py (revision 25455) @@ -1,16 +1,18 @@ import numpy as np + +from GetAreas import * import MatlabFuncs as m from NodeConnectivity import * -from GetAreas import * def adjacency(md): - #ADJACENCY - compute adjacency matrix, list of vertices and list of weights. - # - # function to create the adjacency matrix from the connectivity table. - # - # the required output is: - # md.adj_mat (double [sparse nv x nv], vertex adjacency matrix) - # md.qmu.vertex_weight (double [nv], vertex weights) + """ADJACENCY - compute adjacency matrix, list of vertices and list of weights. + + function to create the adjacency matrix from the connectivity table. + + the required output is: + md.adj_mat (double [sparse nv x nv], vertex adjacency matrix) + md.qmu.vertex_weight (double [nv], vertex weights) + """ indi = np.array([md.mesh.elements[:, 0], md.mesh.elements[:, 1], md.mesh.elements[:, 2]]) Index: /issm/trunk-jpl/src/m/plot/processmesh.py =================================================================== --- /issm/trunk-jpl/src/m/plot/processmesh.py (revision 25454) +++ /issm/trunk-jpl/src/m/plot/processmesh.py (revision 25455) @@ -4,6 +4,5 @@ def processmesh(md, data, options): - ''' - PROCESSMESH - process mesh to be plotted + """PROCESSMESH - process mesh to be plotted Usage: @@ -16,5 +15,5 @@ - Test that output of delaunay matches output of delaunay in MATLAB (see src/m/plot/processmesh.m) - ''' + """ #some checks Index: /issm/trunk-jpl/src/m/qmu/dakota_in_data.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/dakota_in_data.py (revision 25454) +++ /issm/trunk-jpl/src/m/qmu/dakota_in_data.py (revision 25455) @@ -57,5 +57,5 @@ if strcmpi(params.analysis_driver, 'matlab') and params.analysis_components == '': - [pathstr, name, ext] = fileparts(filei) + pathstr, name, ext = fileparts(filei) params.analysis_components = fullfile(pathstr, name + '.py') Index: /issm/trunk-jpl/src/m/qmu/dakota_in_write.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/dakota_in_write.py (revision 25454) +++ /issm/trunk-jpl/src/m/qmu/dakota_in_write.py (revision 25455) @@ -56,5 +56,5 @@ filei = str(eval(input('Dakota input file to write? '))) - [pathstr, name, ext] = fileparts(filei) + pathstr, name, ext = fileparts(filei) if len(ext) == 0: # fileparts only considers '.in' to be the extension, not '.qmu.in' @@ -266,5 +266,5 @@ param_write(fidi, '\t ', 'processors_per_evaluation', '=', '\n', params) if len(params.analysis_components) != 0: - [pathstr, name, ext] = fileparts(params.analysis_components) + pathstr, name, ext = fileparts(params.analysis_components) if ext != '': ext = '.py' Index: /issm/trunk-jpl/src/m/qmu/helpers.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/helpers.py (revision 25454) +++ /issm/trunk-jpl/src/m/qmu/helpers.py (revision 25455) @@ -6,10 +6,10 @@ class struct(object): - '''An empty struct that can be assigned arbitrary attributes''' + """An empty struct that can be assigned arbitrary attributes""" pass class Lstruct(list): - ''' + """ An empty struct that can be assigned arbitrary attributes but can also be accesed as a list. Eg. x.y = 'hello', x[:] = ['w', 'o', 'r', 'l', 'd'] @@ -46,5 +46,5 @@ Sources: -https://github.com/Vectorized/Python-Attribute-List - ''' + """ def __new__(self, *args, **kwargs): @@ -64,5 +64,5 @@ class OrderedStruct(object): - ''' + """ A form of dictionary-like structure that maintains the ordering in which its fields/attributes and their corresponding values were added. @@ -113,12 +113,12 @@ Note: to access internal fields use dir(x) (input fields will be included, but are not technically internals) - ''' + """ def __init__(self, *args): - ''' + """ Provided either nothing or a series of strings, construct a structure that will, when accessed as a list, return its fields in the same order in which they were provided - ''' + """ # keys and values @@ -187,9 +187,9 @@ def __copy__(self): - ''' + """ shallow copy, hard copies of trivial attributes, references to structures like lists/OrderedDicts unless redefined as an entirely different structure - ''' + """ newInstance = type(self)() for k, v in list(self.items()): @@ -198,5 +198,5 @@ def __deepcopy__(self, memo=None): - ''' + """ hard copy of all attributes same thing but call deepcopy recursively @@ -204,5 +204,5 @@ (see https://docs.python.org/2/library/copy.html #copy.deepcopy ) but will generally work in this case - ''' + """ newInstance = type(self)() for k, v in list(self.items()): @@ -232,9 +232,9 @@ def isempty(x): - ''' + """ returns true if object is +/-infinity, NaN, None, '', has length 0, or is an array/matrix composed only of such components (includes mixtures of "empty" types) - ''' + """ if type(x) in [list, np.ndarray, tuple]: @@ -270,9 +270,9 @@ def fieldnames(x, ignore_internals=True): - ''' + """ returns a list of fields of x ignore_internals ignores all fieldnames starting with '_' and is True by default - ''' + """ result = list(vars(x).keys()) @@ -284,17 +284,17 @@ def isfield(x, y, ignore_internals=True): - ''' + """ is y is a field of x? ignore_internals ignores all fieldnames starting with '_' and is True by default - ''' + """ return str(y) in fieldnames(x, ignore_internals) def fileparts(x): - ''' + """ given: "path/path/.../file_name.ext" returns: [path, file_name, ext] (list of strings) - ''' + """ try: a = x[:x.rindex('/')] #path @@ -311,5 +311,5 @@ def fullfile(*args): - ''' + """ usage: fullfile(path, path, ... , file_name + ext) @@ -322,5 +322,5 @@ as final arguments ('file.doc') or ('file' + '.doc') will work ('final', '.doc'), and the like, will not (you'd get 'final/.doc') - ''' + """ result = str(args[0]) for i in range(len(args[1:])): @@ -334,5 +334,5 @@ def findline(fidi, s): - ''' + """ returns full first line containing s (as a string), or None @@ -340,5 +340,5 @@ use str(findline(f, s)).strip() to remove these, str() in case result is None - ''' + """ for line in fidi: if s in line: @@ -348,5 +348,5 @@ def empty_nd_list(shape, filler=0., as_numpy_ndarray=False): - ''' + """ returns a python list of the size/shape given (shape must be int or tuple) the list will be filled with the optional second argument @@ -363,5 +363,5 @@ empty_nd_list((5, 5), 0.0) # returns a 5x5 matrix of 0.0's empty_nd_list(5, None) # returns a 5 long array of NaN - ''' + """ result = np.empty(shape) result.fill(filler) Index: /issm/trunk-jpl/src/m/qmu/postqmu.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/postqmu.py (revision 25454) +++ /issm/trunk-jpl/src/m/qmu/postqmu.py (revision 25455) @@ -1,5 +1,4 @@ from os import getpid, stat from os.path import isfile -import shlex from subprocess import call @@ -62,9 +61,9 @@ # move all the individual function evalutations into zip files if not md.qmu.isdakota: - subproc_args = shlex.split('zip -mq params.in.zip params.in.[1-9]*') + subproc_args = 'zip -mq params.in.zip params.in.[1-9]*' call(subproc_args, shell=True) - subproc_args = shlex.split('zip -mq results.out.zip results.out.[1-9]*') + subproc_args = 'zip -mq results.out.zip results.out.[1-9]*' call(subproc_args, shell=True) - subproc_args = shlex.split('zip -mq matlab.out.zip matlab*.out.[1-9]*') + subproc_args = 'zip -mq matlab.out.zip matlab*.out.[1-9]*' call(subproc_args, shell=True) Index: /issm/trunk-jpl/src/m/shp/shpread.m =================================================================== --- /issm/trunk-jpl/src/m/shp/shpread.m (revision 25454) +++ /issm/trunk-jpl/src/m/shp/shpread.m (revision 25455) @@ -2,20 +2,19 @@ %SHPREAD - read a shape file and build a struct % -% This routine reads a shape file .shp and builds a structure array -% containing the fields x and y corresponding to the coordinates, one for the -% filename of the shp file, for the density, for the nodes, and a field -% closed to indicate if the domain is closed. If this initial shapefile is -% point only, the fields closed and points are omitted. -% The first argument is the .shp file to be read and the second one -% (optional) indicates if the last point shall be read (1 to read it, 0 not -% to). +% This routine reads a shape file .shp and builds a structure array +% containing the fields x and y corresponding to the coordinates, one for the +% filename of the shp file, for the density, for the nodes, and a field +% closed to indicate if the domain is closed. If this initial shapefile is +% point only, the fields closed and points are omitted. % -% Usage: -% Struct=shpread(filename) +% The first argument is the .shp file to be read and the second one +% (optional) indicates if the last point shall be read (1 to read it, 0 not +% to). % -% Example: -% Struct=shpread('domainoutline.shp') +% Usage: +% Struct=shpread(filename) % -% See also EXPDOC, EXPWRITEASVERTICES +% Example: +% Struct=shpread('domainoutline.shp') %recover options Index: /issm/trunk-jpl/src/m/shp/shpread.py =================================================================== --- /issm/trunk-jpl/src/m/shp/shpread.py (revision 25454) +++ /issm/trunk-jpl/src/m/shp/shpread.py (revision 25455) @@ -1,2 +1,3 @@ +from collections import OrderedDict from os import path try: @@ -5,25 +6,22 @@ print("could not import shapefile, PyShp has not been installed, no shapefile reading capabilities enabled") -from helpers import OrderedStruct from pairoptions import pairoptions def shpread(filename, *args): #{{{ - ''' - SHPREAD - read a shape file and build a dict + """SHPREAD - read a shapefile and build a list of shapes - This routine reads a shape file .shp and builds a dict array containing the - fields x and y corresponding to the coordinates, one for the filename - of the shp file, for the density, for the nodes, and a field closed to - indicate if the domain is closed. If this initial shapefile is point only, - the fields closed and points are ommited. - The first argument is the .shp file to be read and the second one - (optional) indicates if the last point shall be read (1 to read it, 0 not - to). + This routine reads a shapefile and builds a list of OrderedDict objects + containing the fields x and y corresponding to the coordinates, one for the + filename of the shp file, for the density, for the nodes, and a field + closed to indicate if the domain is closed. If this initial shapefile is + point only, the fields closed and points are ommitted. The first argument + is the shapefile to be read and the second argument (optional) indicates if + the last point shall be read (1 to read it, 0 not to). The current implementation of shpread depends on PyShp. Usage: - dict = shpread(filename) + list = shpread(filename) Example: @@ -31,11 +29,11 @@ a collection of three files. You specify the base filename of the shapefile or the complete filename of any of the shapefile component - files."" + files." - dict = shpread('domainoutline.shp') + list = shpread('domainoutline.shp') OR - dict = shpread('domainoutline.dbf') + list = shpread('domainoutline.dbf') OR - dict = shpread('domainoutline') + list = shpread('domainoutline') "OR any of the other 5+ formats which are potentially part of a @@ -44,13 +42,25 @@ .shp extension exists. - See also EXPDOC, EXPWRITEASVERTICES - Sources: - https://github.com/GeospatialPython/pyshp + NOTE: + - OrderedDict objects are used instead of OrderedStruct objects (although + addressing in the latter case is closer to the MATLAB struct type) in order + to remain consistent with the pattern established by src/m/exp/expread.py. + TODO: - Create class that can be used to store and pretty print shape structs - (ala OrderedStruct from src/m/qmu/helpers.py). - ''' + (ala OrderedStruct from src/m/qmu/helpers.py). + - Convert returned data structure from list of OrderedDict objects to list + of OrderedStruct objects and remove corresponding note (see also + src/m/exp/expread.py). Also, modify handling of returned data structure in, + - src/m/classes/basin.py + - src/m/classes/boundary.py + - src/m/modules/ExpToLevelSet.py + - src/m/mesh/bamg.py + May also need to modify addressing in corresponding FetchData function, or + create new one, in src/wrappers/ContoursToNodes/ContoursToNodes.cpp. + """ #recover options @@ -67,11 +77,11 @@ shapes = sf.shapes() for i in range(len(shapes)): - Struct = OrderedStruct() + Struct = OrderedDict() shape = shapes[i] if shape.shapeType == shapefile.POINT: - Struct.x = shape.points[0][0] - Struct.y = shape.points[0][1] - Struct.density = 1 - Struct.Geometry = 'Point' + Struct['x'] = shape.points[0][0] + Struct['y'] = shape.points[0][1] + Struct['density'] = 1 + Struct['Geometry'] = 'Point' elif shape.shapeType == shapefile.POLYLINE: num_points = len(shape.points) @@ -82,11 +92,11 @@ x.append(point[0]) y.append(point[1]) - Struct.x = x - Struct.y = y - Struct.nods = num_points - Struct.density = 1 - Struct.closed = 1 - Struct.BoundingBox = shape.bbox - Struct.Geometry = 'Line' + Struct['x'] = x + Struct['y'] = y + Struct['nods'] = num_points + Struct['density'] = 1 + Struct['closed'] = 1 + Struct['BoundingBox'] = shape.bbox + Struct['Geometry'] = 'Line' elif shape.shapeType == shapefile.POLYGON: num_points = len(shape.points) @@ -97,11 +107,11 @@ x.append(point[0]) y.append(point[1]) - Struct.x = x - Struct.y = y - Struct.nods = num_points - Struct.density = 1 - Struct.closed = 1 - Struct.BoundingBox = shape.bbox - Struct.Geometry = 'Polygon' + Struct['x'] = x + Struct['y'] = y + Struct['nods'] = num_points + Struct['density'] = 1 + Struct['closed'] = 1 + Struct['BoundingBox'] = shape.bbox + Struct['Geometry'] = 'Polygon' else: # NOTE: We could do this once before looping over shapes as all @@ -122,5 +132,5 @@ else: setattr(Struct, str(fieldname), sf.record(i)[j - 1]) # cast to string removes "u" from "u'fieldname'" - Struct.name = name + Struct['name'] = name Structs.append(Struct) @@ -128,6 +138,6 @@ if invert: for i in range(len(Structs)): - Structs[i].x = np.flipud(Structs[i].x) - Structs[i].y = np.flipud(Structs[i].y) + Structs[i]['x'] = np.flipud(Structs[i]['x']) + Structs[i]['y'] = np.flipud(Structs[i]['y']) return Structs Index: /issm/trunk-jpl/src/wrappers/ContourToMesh/ContourToMesh.cpp =================================================================== --- /issm/trunk-jpl/src/wrappers/ContourToMesh/ContourToMesh.cpp (revision 25454) +++ /issm/trunk-jpl/src/wrappers/ContourToMesh/ContourToMesh.cpp (revision 25455) @@ -67,5 +67,5 @@ /*Run interpolation routine: */ - ContourToMeshx( &in_nod,&in_elem,index,x,y,contours,interptype,nel,nods,edgevalue); + ContourToMeshx(&in_nod,&in_elem,index,x,y,contours,interptype,nel,nods,edgevalue); /* output: */ Index: /issm/trunk-jpl/src/wrappers/InterpFromMesh2d/InterpFromMesh2d.cpp =================================================================== --- /issm/trunk-jpl/src/wrappers/InterpFromMesh2d/InterpFromMesh2d.cpp (revision 25454) +++ /issm/trunk-jpl/src/wrappers/InterpFromMesh2d/InterpFromMesh2d.cpp (revision 25455) @@ -42,5 +42,7 @@ /*contours*/ int i; + #ifdef _HAVE_MATLAB_MODULES_ mxArray *matlabstructure = NULL; + #endif Contour **contours=NULL; int numcontours; @@ -59,8 +61,10 @@ /*checks on arguments on the matlab side: */ + #ifdef _HAVE_MATLAB_MODULES_ if(nlhs!=NLHS){ InterpFromMesh2dUsage(); _error_("InterpFromMeshToMesh2dUsage usage error"); } + #endif if((nrhs!=6) && (nrhs!=7) && (nrhs!=8)){ InterpFromMesh2dUsage(); @@ -85,5 +89,5 @@ } - if(nrhs>=8){ + if(nrhs==8){ #ifdef _HAVE_MATLAB_MODULES_ Index: /issm/trunk-jpl/src/wrappers/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp =================================================================== --- /issm/trunk-jpl/src/wrappers/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp (revision 25454) +++ /issm/trunk-jpl/src/wrappers/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp (revision 25455) @@ -31,15 +31,15 @@ /*Intermediaties*/ - int *index = NULL; - double *x_data = NULL; - double *y_data = NULL; - double *data = NULL; + int* index = NULL; + double* x_data = NULL; + double* y_data = NULL; + double* data = NULL; int nods_data,nels_data; int M_data,N_data; - double *x_interp = NULL; - double *y_interp = NULL; + double* x_interp = NULL; + double* y_interp = NULL; int N_interp; - Options *options = NULL; - double *data_interp = NULL; + Options* options = NULL; + double* data_interp = NULL; int test1,test2,test; Index: /issm/trunk-jpl/src/wrappers/python/Makefile.am =================================================================== --- /issm/trunk-jpl/src/wrappers/python/Makefile.am (revision 25454) +++ /issm/trunk-jpl/src/wrappers/python/Makefile.am (revision 25455) @@ -36,4 +36,5 @@ ElementConnectivity_python.la \ ExpToLevelSet_python.la \ + InterpFromMesh2d_python.la \ InterpFromMeshToMesh2d_python.la \ InterpFromMeshToMesh3d_python.la \ @@ -144,4 +145,8 @@ ExpToLevelSet_python_la_LIBADD = ${deps} $(PETSCLIB) $(HDF5LIB) $(BLASLAPACKLIB) $(MPILIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB) $(NEOPZLIB) +InterpFromMesh2d_python_la_SOURCES = ../InterpFromMesh2d/InterpFromMesh2d.cpp +InterpFromMesh2d_python_la_CXXFLAGS = ${AM_CXXFLAGS} +InterpFromMesh2d_python_la_LIBADD = ${deps} $(PETSCLIB) $(HDF5LIB) $(BLASLAPACKLIB) $(MPILIB) $(MULTITHREADINGLIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB) + InterpFromGridToMesh_python_la_SOURCES = ../InterpFromGridToMesh/InterpFromGridToMesh.cpp InterpFromGridToMesh_python_la_CXXFLAGS = ${AM_CXXFLAGS} Index: /issm/trunk-jpl/test/NightlyRun/test2004.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 25454) +++ /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 25455) @@ -77,7 +77,7 @@ alreadyloaded=0; %}}} -for i=sl.basinindx('basin','all'), - - bas=sl.basins{i}; +for ind=sl.basinindx('basin','all'), + + bas=sl.basins{ind}; disp(sprintf('Meshing basin %s\n',bas.name)); @@ -85,8 +85,8 @@ domain=bas.contour(); - %recover coastline inside basin, using GSHHS_c_L1. It's a lat,long file, hence epsg 4326 - coastline=bas.shapefilecrop('shapefile',gshhsshapefile,'epsgshapefile',4326,'threshold',threshold); - - %mesh: + %recover coastline inside basin, using GSHHS_c_L1. It's a lat/long file, hence epsg 4326 + coastline=bas.shapefilecrop('shapefile',gshhsshapefile,'epsgshapefile',4326,'threshold',threshold); + + %mesh: md=bamg(model,'domain',domain,'subdomains',coastline,'hmin',hmin,'hmax',hmax,defaultoptions{:}); %plotmodel(md,'data','mesh');pause(1); @@ -108,9 +108,9 @@ %Parameterization: %Parameterize ice sheets : {{{ - for ind=sl.basinindx('continent',{'antarctica'}), disp(sprintf('Parameterizing basin %s\n', sl.icecaps{ind}.miscellaneous.name)); - md=sl.icecaps{ind}; bas=sl.basins{ind}; + md=sl.icecaps{ind}; + bas=sl.basins{ind}; %masks : %{{{ %ice levelset from domain outlines: @@ -125,5 +125,5 @@ %}}} %latlong: % {{{ - [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326'); + [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326'); %}}} %geometry: {{{ @@ -137,4 +137,5 @@ if bas.iscontinentany('antarctica'), if testagainst2002, + % TODO: Check if the following works as expected: 'pos' is empty, so nothing is assigned to 'md.solidearth.surfaceload.icethicknesschange(pos)' md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); %antarctica @@ -146,11 +147,19 @@ md.solidearth.surfaceload.icethicknesschange(pos)=-100*ratio; else - delH=textread('../Data/AIS_delH_trend.txt'); - longAIS=delH(:,1); latAIS=delH(:,2); delHAIS=delH(:,3); index=delaunay(longAIS,latAIS); - lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360; + in_fileID=fopen('../Data/AIS_delH_trend.txt', 'r'); + delH=textscan(in_fileID, '%f %f %f'); + fclose(in_fileID); + longAIS=delH{:,1}; + latAIS=delH{:,2}; + delHAIS=delH{:,3}; + index=delaunay(longAIS,latAIS); + lat=md.mesh.lat; + long=md.mesh.long+360; + pos=find(long>360); + long(pos)=long(pos)-360; delHAIS=InterpFromMesh2d(index,longAIS,latAIS,delHAIS,long,lat); - northpole=find_point(md.mesh.long,md.mesh.lat,0,90); delHAIS(northpole)=0; - - md.solidearth.surfaceload.icethicknesschange=delHAIS(md.mesh.elements)*[1;1;1]/3/100; + northpole=find_point(md.mesh.long,md.mesh.lat,0,90); + delHAIS(northpole)=0; + md.solidearth.surfaceload.icethicknesschange=mean(delHAIS(md.mesh.elements),2)/100; end @@ -172,11 +181,10 @@ end %}}} + % ParameterizeContinents {{{ - -sl.basinindx('continent',{'hemisphereeast','hemispherewest'}) - for ind=sl.basinindx('continent',{'hemisphereeast','hemispherewest'}), disp(sprintf('Masks for basin %s\n', sl.icecaps{ind}.miscellaneous.name)); - md=sl.icecaps{ind}; bas=sl.basins{ind}; + md=sl.icecaps{ind}; + bas=sl.basins{ind}; %recover lat,long: @@ -184,6 +192,7 @@ %mask: %{{{ - %Figure out mask from initial mesh: deal with land and ocean masks (land include grounded ice). %{{{ - %first, transform land element mask into vertex driven one. + %Figure out mask from initial mesh: deal with land and ocean masks (land + %includes grounded ice). %{{{ + %first, transform land element mask into vertex-driven one land=md.private.bamg.landmask; land_mask=-ones(md.mesh.numberofvertices,1); @@ -192,5 +201,5 @@ land_mask(md.mesh.elements(landels,:))=1; - %gothrough edges of each land element + % Gothrough edges of each land element connectedels=md.mesh.elementconnectivity(landels,:); connectedisonocean=~land(connectedels); @@ -200,7 +209,10 @@ landelsconocean=landels(find(sumconnectedisonocean)); - ind1=[md.mesh.elements(landelsconocean,1); md.mesh.elements(landelsconocean,2); md.mesh.elements(landelsconocean,3)]; - ind2=[md.mesh.elements(landelsconocean,2); md.mesh.elements(landelsconocean,3); md.mesh.elements(landelsconocean,1)]; - + ind1=[md.mesh.elements(landelsconocean,1); + md.mesh.elements(landelsconocean,2); + md.mesh.elements(landelsconocean,3)]; + ind2=[md.mesh.elements(landelsconocean,2); + md.mesh.elements(landelsconocean,3); + md.mesh.elements(landelsconocean,1)]; %edge ind1 and ind2: @@ -223,16 +235,15 @@ % {{{ %greenland - pos=find(md.mesh.lat > 70 & md.mesh.lat < 80 & md.mesh.long>-60 & md.mesh.long<-30); + pos=find(md.mesh.lat > 70 & md.mesh.lat < 80 & md.mesh.long>-60 & md.mesh.long<-30); md.mask.ice_levelset(pos)=-1; % }}} end - % }}} %}}} %slr loading/calibration: {{{ + md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); if testagainst2002, % {{{ - md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); %greenland late=sum(md.mesh.lat(md.mesh.elements),2)/3; @@ -243,19 +254,28 @@ %correct mask: - md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; + md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % }}} else - - md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); - delH=textread('../Data/GIS_delH_trend.txt'); - longGIS=delH(:,1); latGIS=delH(:,2); delHGIS=delH(:,3); index=delaunay(longGIS,latGIS); - lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360; + longGIS=delH(:,1); + latGIS=delH(:,2); + delHGIS=delH(:,3); + index=delaunay(longGIS,latGIS); + lat=md.mesh.lat; + long=md.mesh.long+360; + pos=find(long>360); + long(pos)=long(pos)-360; delHGIS=InterpFromMeshToMesh2d(index,longGIS,latGIS,delHGIS,long,lat); delHGISe=delHGIS(md.mesh.elements)*[1;1;1]/3; - + delH=textread('../Data/GLA_delH_trend_15regions.txt'); - longGLA=delH(:,1); latGLA=delH(:,2); delHGLA=sum(delH(:,3:end),2); index=delaunay(longGLA,latGLA); - lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360; + longGLA=delH(:,1); + latGLA=delH(:,2); + delHGLA=sum(delH(:,3:end),2); + index=delaunay(longGLA,latGLA); + lat=md.mesh.lat; + long=md.mesh.long+360; + pos=find(long>360); + long(pos)=long(pos)-360; delHGLA=InterpFromMeshToMesh2d(index,longGLA,latGLA,delHGLA,long,lat); delHGLAe=delHGLA(md.mesh.elements)*[1;1;1]/3; @@ -293,4 +313,5 @@ end % }}} + %%Assemble Earth in 3D {{{ @@ -300,5 +321,5 @@ loneedgesdetect=0; -%create earth model by concatenating all the icecaps in 3d: +%create Earth model by concatenating all the icecaps in 3D: sl.caticecaps('tolerance',tolerance,'loneedgesdetect',loneedgesdetect); Index: /issm/trunk-jpl/test/NightlyRun/test423.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test423.m (revision 25454) +++ /issm/trunk-jpl/test/NightlyRun/test423.m (revision 25455) @@ -8,6 +8,6 @@ pos=find(rad==min(rad)); md.mesh.x(pos)=0.; md.mesh.y(pos)=0.; %the closest node to the center is changed to be exactly at the center -xelem=md.mesh.x(md.mesh.elements)*[1;1;1]/3.; -yelem=md.mesh.y(md.mesh.elements)*[1;1;1]/3.; +xelem=mean(md.mesh.x(md.mesh.elements),2); +yelem=mean(md.mesh.y(md.mesh.elements),2); rad=sqrt(xelem.^2+yelem.^2); flags=zeros(md.mesh.numberofelements,1); Index: /issm/trunk-jpl/test/NightlyRun/test423.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test423.py (revision 25454) +++ /issm/trunk-jpl/test/NightlyRun/test423.py (revision 25455) @@ -19,6 +19,6 @@ md.mesh.x[pos] = 0. md.mesh.y[pos] = 0. #the closest node to the center is changed to be exactly at the center -xelem = np.mean(md.mesh.x[md.mesh.elements.astype(int) - 1], axis=1) -yelem = np.mean(md.mesh.y[md.mesh.elements.astype(int) - 1], axis=1) +xelem = np.mean(md.mesh.x[md.mesh.elements - 1], axis=1) +yelem = np.mean(md.mesh.y[md.mesh.elements - 1], axis=1) rad = np.sqrt(xelem**2 + yelem**2) flags = np.zeros(md.mesh.numberofelements) Index: /issm/trunk-jpl/test/NightlyRun/test433.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test433.m (revision 25454) +++ /issm/trunk-jpl/test/NightlyRun/test433.m (revision 25455) @@ -8,6 +8,6 @@ pos=find(rad==min(rad)); md.mesh.x(pos)=0.; md.mesh.y(pos)=0.; %the closest node to the center is changed to be exactly at the center -xelem=md.mesh.x(md.mesh.elements)*[1;1;1]/3.; -yelem=md.mesh.y(md.mesh.elements)*[1;1;1]/3.; +xelem=mean(md.mesh.x(md.mesh.elements), 2); +yelem=mean(md.mesh.y(md.mesh.elements), 2); rad=sqrt(xelem.^2+yelem.^2); flags=zeros(md.mesh.numberofelements,1); Index: /issm/trunk-jpl/test/NightlyRun/test433.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test433.py (revision 25454) +++ /issm/trunk-jpl/test/NightlyRun/test433.py (revision 25455) @@ -19,6 +19,6 @@ md.mesh.x[pos] = 0. md.mesh.y[pos] = 0. #the closest node to the center is changed to be exactly at the center -xelem = np.mean(md.mesh.x[md.mesh.elements.astype(int) - 1], axis=1) -yelem = np.mean(md.mesh.y[md.mesh.elements.astype(int) - 1], axis=1) +xelem = np.mean(md.mesh.x[md.mesh.elements - 1], axis=1) +yelem = np.mean(md.mesh.y[md.mesh.elements - 1], axis=1) rad = np.sqrt(xelem**2 + yelem**2) flags = np.zeros(md.mesh.numberofelements) Index: /issm/trunk-jpl/test/Par/79North.py =================================================================== --- /issm/trunk-jpl/test/Par/79North.py (revision 25454) +++ /issm/trunk-jpl/test/Par/79North.py (revision 25455) @@ -19,8 +19,8 @@ thickness = numpy.array(archread('../Data/79North.arch', 'thickness')) -md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)[0][:, 0] -md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)[0][:, 0] -md.geometry.surface = InterpFromMeshToMesh2d(index, x, y, surface, md.mesh.x, md.mesh.y)[0][:, 0] -md.geometry.thickness = InterpFromMeshToMesh2d(index, x, y, thickness, md.mesh.x, md.mesh.y)[0][:, 0] +md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)[:, 0] +md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)[:, 0] +md.geometry.surface = InterpFromMeshToMesh2d(index, x, y, surface, md.mesh.x, md.mesh.y)[:, 0] +md.geometry.thickness = InterpFromMeshToMesh2d(index, x, y, thickness, md.mesh.x, md.mesh.y)[:, 0] md.geometry.base = md.geometry.surface - md.geometry.thickness Index: /issm/trunk-jpl/test/Par/Pig.py =================================================================== --- /issm/trunk-jpl/test/Par/Pig.py (revision 25454) +++ /issm/trunk-jpl/test/Par/Pig.py (revision 25455) @@ -20,12 +20,12 @@ bed = np.array(archread('../Data/Pig.arch', 'bed')) -md.inversion.vx_obs = InterpFromMeshToMesh2d(index, x, y, vx_obs, md.mesh.x, md.mesh.y)[0][:, 0] -md.inversion.vy_obs = InterpFromMeshToMesh2d(index, x, y, vy_obs, md.mesh.x, md.mesh.y)[0][:, 0] -md.geometry.surface = InterpFromMeshToMesh2d(index, x, y, surface, md.mesh.x, md.mesh.y)[0][:, 0] -md.geometry.thickness = InterpFromMeshToMesh2d(index, x, y, thickness, md.mesh.x, md.mesh.y)[0][:, 0] +md.inversion.vx_obs = InterpFromMeshToMesh2d(index, x, y, vx_obs, md.mesh.x, md.mesh.y)[:, 0] +md.inversion.vy_obs = InterpFromMeshToMesh2d(index, x, y, vy_obs, md.mesh.x, md.mesh.y)[:, 0] +md.geometry.surface = InterpFromMeshToMesh2d(index, x, y, surface, md.mesh.x, md.mesh.y)[:, 0] +md.geometry.thickness = InterpFromMeshToMesh2d(index, x, y, thickness, md.mesh.x, md.mesh.y)[:, 0] md.geometry.base = md.geometry.surface - md.geometry.thickness md.geometry.bed = np.array(md.geometry.base) pos = np.where(md.mask.ocean_levelset < 0.) -md.geometry.bed[pos] = InterpFromMeshToMesh2d(index, x, y, bed, md.mesh.x[pos], md.mesh.y[pos])[0][:, 0] +md.geometry.bed[pos] = InterpFromMeshToMesh2d(index, x, y, bed, md.mesh.x[pos], md.mesh.y[pos])[:, 0] md.initialization.vx = md.inversion.vx_obs md.initialization.vy = md.inversion.vy_obs Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetConstrained.py (revision 25454) +++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.py (revision 25455) @@ -28,6 +28,6 @@ index = archread('../Data/SquareSheetConstrained.arch', 'index').astype(int) -md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)[0] -md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)[0] +md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y) +md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y) md.initialization.vz = np.zeros((md.mesh.numberofvertices)) md.initialization.pressure = np.zeros((md.mesh.numberofvertices)) Index: /issm/trunk-jpl/test/Par/SquareSheetConstrainedCO2.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetConstrainedCO2.py (revision 25454) +++ /issm/trunk-jpl/test/Par/SquareSheetConstrainedCO2.py (revision 25455) @@ -44,6 +44,6 @@ index = archread('../Data/SquareSheetConstrained.arch', 'index').astype(int) -[md.initialization.vx] = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y) -[md.initialization.vy] = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y) +md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y) +md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y) md.initialization.vz = np.zeros((md.mesh.numberofvertices)) md.initialization.pressure = np.zeros((md.mesh.numberofvertices)) Index: /issm/trunk-jpl/test/Par/SquareSheetShelf.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetShelf.py (revision 25454) +++ /issm/trunk-jpl/test/Par/SquareSheetShelf.py (revision 25455) @@ -31,6 +31,6 @@ index = np.array(archread('../Data/SquareSheetShelf.arch', 'index')).astype(int) -[md.initialization.vx] = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y) -[md.initialization.vy] = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y) +md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y) +md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y) md.initialization.vz = np.zeros((md.mesh.numberofvertices)) md.initialization.pressure = np.zeros((md.mesh.numberofvertices)) Index: /issm/trunk-jpl/test/Par/SquareShelf2.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelf2.py (revision 25454) +++ /issm/trunk-jpl/test/Par/SquareShelf2.py (revision 25455) @@ -32,6 +32,6 @@ #dbg - end -[md.initialization.vx] = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y) -[md.initialization.vy] = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y) +md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y) +md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y) md.initialization.vz = np.zeros((md.mesh.numberofvertices)) md.initialization.pressure = np.zeros((md.mesh.numberofvertices)) Index: /issm/trunk-jpl/test/Par/SquareShelfConstrained.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelfConstrained.py (revision 25454) +++ /issm/trunk-jpl/test/Par/SquareShelfConstrained.py (revision 25455) @@ -28,6 +28,6 @@ index = np.array(archread('../Data/SquareShelfConstrained.arch', 'index').astype(int)) -[md.initialization.vx] = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y) -[md.initialization.vy] = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y) +md.initialization.vx = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y) +md.initialization.vy = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y) md.initialization.vz = np.zeros((md.mesh.numberofvertices)) md.initialization.pressure = np.zeros((md.mesh.numberofvertices))