diff --git a/flopy/discretization/vertexgrid.py b/flopy/discretization/vertexgrid.py index 1dd4514374..9d11f589b2 100644 --- a/flopy/discretization/vertexgrid.py +++ b/flopy/discretization/vertexgrid.py @@ -32,13 +32,14 @@ class for a vertex model grid returns vertices for a single cell at cellid. """ - def __init__(self, vertices=None, cell2d=None, top=None, botm=None, idomain=None, - lenuni=None, epsg=None, proj4=None, prj=None, xoff=0.0, - yoff=0.0, angrot=0.0, grid_type='vertex', - nlay=None, ncpl=None): + def __init__(self, vertices=None, cell2d=None, top=None, + botm=None, idomain=None, lenuni=None, epsg=None, proj4=None, + prj=None, xoff=0.0, yoff=0.0, angrot=0.0, grid_type='vertex', + nlay=None, ncpl=None, cell1d=None): super(VertexGrid, self).__init__(grid_type, top, botm, idomain, lenuni, epsg, proj4, prj, xoff, yoff, angrot) self._vertices = vertices + self._cell1d = cell1d self._cell2d = cell2d self._top = top self._botm = botm @@ -52,26 +53,32 @@ def __init__(self, vertices=None, cell2d=None, top=None, botm=None, idomain=None @property def is_valid(self): - if self._vertices is not None and self._cell2d is not None: + if self._vertices is not None and (self._cell2d is not None or + self._cell1d is not None): return True return False @property def is_complete(self): - if self._vertices is not None and self._cell2d is not None and \ + if self._vertices is not None and (self._cell2d is not None or + self._cell1d is not None) and \ super(VertexGrid, self).is_complete: return True return False @property def nlay(self): - if self._botm is not None: + if self._cell1d is not None: + return 1 + elif self._botm is not None: return len(self._botm) else: return self._nlay @property def ncpl(self): + if self._cell1d is not None: + return len(self._cell1d) if self._botm is not None: return len(self._botm[0]) else: @@ -233,34 +240,62 @@ def _build_grid_geometry_info(self): cache_index_cc = 'cellcenters' cache_index_vert = 'xyzgrid' - vertexdict = {v[0]: [v[1], v[2]] - for v in self._vertices} xcenters = [] ycenters = [] xvertices = [] yvertices = [] - # build xy vertex and cell center info - for cell2d in self._cell2d: - cell2d = tuple(cell2d) - xcenters.append(cell2d[1]) - ycenters.append(cell2d[2]) - - vert_number = [] - for i in cell2d[4:]: - if i is not None: - vert_number.append(int(i)) + if self._cell1d is not None: + zcenters = [] + zvertices = [] + vertexdict = {v[0]: [v[1], v[2], v[3]] + for v in self._vertices} + for cell1d in self._cell1d: + cell1d = tuple(cell1d) + xcenters.append(cell1d[1]) + ycenters.append(cell1d[2]) + zcenters.append(cell1d[3]) + + vert_number = [] + for i in cell1d[3:]: + if i is not None: + vert_number.append(int(i)) + + xcellvert = [] + ycellvert = [] + zcellvert = [] + for ix in vert_number: + xcellvert.append(vertexdict[ix][0]) + ycellvert.append(vertexdict[ix][1]) + zcellvert.append(vertexdict[ix][2]) + xvertices.append(xcellvert) + yvertices.append(ycellvert) + zvertices.append(zcellvert) - xcellvert = [] - ycellvert = [] - for ix in vert_number: - xcellvert.append(vertexdict[ix][0]) - ycellvert.append(vertexdict[ix][1]) - xvertices.append(xcellvert) - yvertices.append(ycellvert) - - # build z cell centers - zvertices, zcenters = self._zcoords() + else: + vertexdict = {v[0]: [v[1], v[2]] + for v in self._vertices} + # build xy vertex and cell center info + for cell2d in self._cell2d: + cell2d = tuple(cell2d) + xcenters.append(cell2d[1]) + ycenters.append(cell2d[2]) + + vert_number = [] + for i in cell2d[4:]: + if i is not None: + vert_number.append(int(i)) + + xcellvert = [] + ycellvert = [] + for ix in vert_number: + xcellvert.append(vertexdict[ix][0]) + ycellvert.append(vertexdict[ix][1]) + xvertices.append(xcellvert) + yvertices.append(ycellvert) + + # build z cell centers + zvertices, zcenters = self._zcoords() if self._has_ref_coordinates: # transform x and y @@ -270,7 +305,8 @@ def _build_grid_geometry_info(self): # vertices are a list within a list for xcellvertices, ycellvertices in zip(xvertices, yvertices): xcellvertices, \ - ycellvertices = self.get_coords(xcellvertices, ycellvertices) + ycellvertices = self.get_coords(xcellvertices, + ycellvertices) xvertxform.append(xcellvertices) yvertxform.append(ycellvertices) xvertices = xvertxform diff --git a/flopy/export/shapefile_utils.py b/flopy/export/shapefile_utils.py index df3c15cd1f..0e8c348c73 100755 --- a/flopy/export/shapefile_utils.py +++ b/flopy/export/shapefile_utils.py @@ -51,9 +51,9 @@ def write_gridlines_shapefile(filename, mg): " instead.", category=DeprecationWarning) else: - grid_lines = mg.grid_lines() + grid_lines = mg.grid_lines for i, line in enumerate(grid_lines): - wr.poly([line]) + wr.line([line]) wr.record(i) wr.close() diff --git a/flopy/mf6/coordinates/modeldimensions.py b/flopy/mf6/coordinates/modeldimensions.py index fb919b4e98..c2411b342a 100644 --- a/flopy/mf6/coordinates/modeldimensions.py +++ b/flopy/mf6/coordinates/modeldimensions.py @@ -336,6 +336,10 @@ def _create_model_grid(self, grid_type): elif grid_type == DiscretizationType.DISU: self._model_grid = UnstructuredModelGrid(self.model_name, self.simulation_data) + elif grid_type == DiscretizationType.DISL: + self._model_grid = ModelGrid(self.model_name, + self.simulation_data, + DiscretizationType.DISL) else: self._model_grid = ModelGrid(self.model_name, self.simulation_data, diff --git a/flopy/mf6/coordinates/modelgrid.py b/flopy/mf6/coordinates/modelgrid.py index 76840f2530..d44b0c050e 100644 --- a/flopy/mf6/coordinates/modelgrid.py +++ b/flopy/mf6/coordinates/modelgrid.py @@ -401,6 +401,10 @@ def get_grid_type(simulation_data, model_name): 'disu{}'.format(structure.get_version_string()), 0) is not None: return DiscretizationType.DISU + elif package_recarray.search_data( + 'disl{}'.format(structure.get_version_string()), + 0) is not None: + return DiscretizationType.DISL return DiscretizationType.UNDEFINED @@ -411,6 +415,9 @@ def get_idomain(self): elif self._grid_type == DiscretizationType.DISV: return self._simulation_data.mfdata[ (self._model_name, 'disv', 'griddata', 'idomain')].get_data() + elif self._grid_type == DiscretizationType.DISL: + return self._simulation_data.mfdata[ + (self._model_name, 'disl', 'griddata', 'idomain')].get_data() elif self._grid_type == DiscretizationType.DISU: except_str = 'ERROR: Can not return idomain for model {}. This ' \ 'model uses a DISU grid that does not ' \ @@ -447,10 +454,11 @@ def get_horizontal_cross_section_dim_arrays(self): np.arange(1, self.num_columns() + 1, 1, np.int32)] elif self.grid_type() == DiscretizationType.DISV: return [np.arange(1, self.num_cells_per_layer() + 1, 1, np.int32)] - elif self.grid_type() == DiscretizationType.DISU: + elif self.grid_type() == DiscretizationType.DISU or \ + self.grid_type() == DiscretizationType.DISL: except_str = 'ERROR: Can not get horizontal plane arrays for ' \ - 'model "{}" DISU grid. DISU grids do not support ' \ - 'individual layers.'.format(self._model_name) + 'model "{}" grid. DISU and DISL grids do not ' \ + 'support individual layers.'.format(self._model_name) print(except_str) raise MFGridException(except_str) @@ -459,7 +467,8 @@ def get_model_dim(self): return [self.num_layers(), self.num_rows(), self.num_columns()] elif self.grid_type() == DiscretizationType.DISV: return [self.num_layers(), self.num_cells_per_layer()] - elif self.grid_type() == DiscretizationType.DISU: + elif self.grid_type() == DiscretizationType.DISU or \ + self.grid_type() == DiscretizationType.DISL: return [self.num_cells()] def get_model_dim_arrays(self): @@ -470,7 +479,8 @@ def get_model_dim_arrays(self): elif self.grid_type() == DiscretizationType.DISV: return [np.arange(1, self.num_layers() + 1, 1, np.int32), np.arange(1, self.num_cells_per_layer() + 1, 1, np.int32)] - elif self.grid_type() == DiscretizationType.DISU: + elif self.grid_type() == DiscretizationType.DISU or \ + self.grid_type() == DiscretizationType.DISL: return [np.arange(1, self.num_cells() + 1, 1, np.int32)] def get_row_array(self): @@ -487,7 +497,8 @@ def get_horizontal_cross_section_dim_names(self): return ['row', 'column'] elif self.grid_type() == DiscretizationType.DISV: return ['layer_cell_num'] - elif self.grid_type() == DiscretizationType.DISU: + elif self.grid_type() == DiscretizationType.DISU or \ + self.grid_type() == DiscretizationType.DISL: except_str = 'ERROR: Can not get layer dimension name for model ' \ '"{}" DISU grid. DISU grids do not support ' \ 'layers.'.format(self._model_name) @@ -499,7 +510,8 @@ def get_model_dim_names(self): return ['layer', 'row', 'column'] elif self.grid_type() == DiscretizationType.DISV: return ['layer', 'layer_cell_num'] - elif self.grid_type() == DiscretizationType.DISU: + elif self.grid_type() == DiscretizationType.DISU or \ + self.grid_type() == DiscretizationType.DISL: return ['node'] def get_num_spatial_coordinates(self): @@ -507,15 +519,10 @@ def get_num_spatial_coordinates(self): return 3 elif self.grid_type() == DiscretizationType.DISV: return 2 - elif self.grid_type() == DiscretizationType.DISU: + elif self.grid_type() == DiscretizationType.DISU or \ + self.grid_type() == DiscretizationType.DISL: return 1 - def change_grid_spacing(self, spacing_factor): - self.test = 1 - - def change_discretization_type(self, new_dis_type): - self.test = 1 - def num_rows(self): if self.grid_type() != DiscretizationType.DIS: except_str = 'ERROR: Model "{}" does not have rows. Can not ' \ @@ -567,7 +574,8 @@ def num_layers(self): elif self.grid_type() == DiscretizationType.DISV: return self._simulation_data.mfdata[ (self._model_name, 'disv', 'dimensions', 'nlay')].get_data() - elif self.grid_type() == DiscretizationType.DISU: + elif self.grid_type() == DiscretizationType.DISU or \ + self.grid_type() == DiscretizationType.DISL: return None def num_cells(self): @@ -578,6 +586,9 @@ def num_cells(self): elif self.grid_type() == DiscretizationType.DISU: return self._simulation_data.mfdata[ (self._model_name, 'disu', 'dimensions', 'nodes')].get_data() + elif self.grid_type() == DiscretizationType.DISL: + return self._simulation_data.mfdata[ + (self._model_name, 'disl', 'dimensions', 'nodes')].get_data() def get_all_model_cells(self): model_cells = [] @@ -592,7 +603,8 @@ def get_all_model_cells(self): for layer_cellid in range(0, self.num_rows()): model_cells.append((layer + 1, layer_cellid + 1)) return model_cells - elif self.grid_type() == DiscretizationType.DISU: + elif self.grid_type() == DiscretizationType.DISU or \ + self.grid_type() == DiscretizationType.DISL: for node in range(0, self.num_cells()): model_cells.append(node + 1) return model_cells diff --git a/flopy/mf6/data/mfdatastorage.py b/flopy/mf6/data/mfdatastorage.py index ff0e7ea4e1..b193e745b5 100644 --- a/flopy/mf6/data/mfdatastorage.py +++ b/flopy/mf6/data/mfdatastorage.py @@ -1535,7 +1535,8 @@ def _verify_list(self, data): model_grid = self.data_dimensions.get_model_grid() cellid_size = model_grid.\ get_num_spatial_coordinates() - if len(data_line[index]) != cellid_size: + if cellid_size != 1 and \ + len(data_line[index]) != cellid_size: message = 'Cellid "{}" contains {} integer(s). ' \ 'Expected a cellid containing {} ' \ 'integer(s) for grid type' \ diff --git a/flopy/mf6/data/mfstructure.py b/flopy/mf6/data/mfstructure.py index 46cc66f898..e9dac9f527 100644 --- a/flopy/mf6/data/mfstructure.py +++ b/flopy/mf6/data/mfstructure.py @@ -65,8 +65,9 @@ def __init__(self): # FIX: Transport - multi packages are hard coded self.multi_package = {'gwfmvr': 0, 'exggwfgwf': 0, 'gwfchd': 0, 'gwfrch': 0, 'gwfwel': 0, - 'gwfdrn': 0, 'gwfriv': 0, 'utlobs': 0, - 'utlts': 0, 'utltas': 0} + 'gwfdrn': 0, 'gwfriv': 0, 'lnfcgeo': 0, + 'lnfrgeo': 0, 'lnfngeo':0, + 'utlobs': 0, 'utlts': 0, 'utltas': 0} def get_file_list(self): file_order = ['sim-nam', # dfn completed tex updated @@ -77,6 +78,7 @@ def get_file_list(self): 'gwf-dis', # dfn completed tex updated 'gwf-disv', # dfn completed tex updated 'gwf-disu', # dfn completed tex updated + 'lnf-disl', # dfn completed tex updated 'gwf-ic', # dfn completed tex updated 'gwf-npf', # dfn completed tex updated 'gwf-sto', # dfn completed tex updated diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py index cf5d857c66..6b32581917 100644 --- a/flopy/mf6/mfmodel.py +++ b/flopy/mf6/mfmodel.py @@ -368,6 +368,33 @@ def modelgrid(self): epsg=self._modelgrid.epsg, xoff=self._modelgrid.xoffset, yoff=self._modelgrid.yoffset, angrot=self._modelgrid.angrot, nodes=dis.nodes.get_data()) + elif self.get_grid_type() == DiscretizationType.DISL: + dis = self.get_package('disl') + if not hasattr(dis, '_init_complete'): + if not hasattr(dis, 'cell1d'): + # disv package has not yet been initialized + return self._modelgrid + else: + # disv package has been partially initialized + self._modelgrid = VertexGrid(vertices=dis.vertices.array, + cell1d=dis.cell1d.array, + top=None, + botm=None, + idomain=None, + lenuni=None, + proj4=self._modelgrid.proj4, + epsg=self._modelgrid.epsg, + xoff=self._modelgrid.xoffset, + yoff=self._modelgrid.yoffset, + angrot=self._modelgrid.angrot) + else: + self._modelgrid = VertexGrid( + vertices=dis.vertices.array, cell1d=dis.cell1d.array, + top=dis.top.array, botm=dis.botm.array, + idomain=dis.idomain.array, lenuni=dis.length_units.array, + proj4=self._modelgrid.proj4, epsg=self._modelgrid.epsg, + xoff=self._modelgrid.xoffset, yoff=self._modelgrid.yoffset, + angrot=self._modelgrid.angrot) else: return self._modelgrid @@ -671,6 +698,10 @@ def get_grid_type(self): 'disu{}'.format(structure.get_version_string()), 0) is not None: return DiscretizationType.DISU + elif package_recarray.search_data( + 'disl{}'.format(structure.get_version_string()), + 0) is not None: + return DiscretizationType.DISL return DiscretizationType.UNDEFINED diff --git a/flopy/mf6/utils/mfenums.py b/flopy/mf6/utils/mfenums.py index dee5ea7c34..3b0b0a164c 100644 --- a/flopy/mf6/utils/mfenums.py +++ b/flopy/mf6/utils/mfenums.py @@ -9,3 +9,4 @@ class DiscretizationType(Enum): DIS = 1 DISV = 2 DISU = 3 + DISL = 4