From 0425c2244bccf514d836f61dc168e361bd64c012 Mon Sep 17 00:00:00 2001 From: jfranc Date: Wed, 20 Sep 2023 03:42:31 -0700 Subject: [PATCH 1/2] Adding adapted cartesian grid option --- geometries/README.md | 22 ++++- geometries/make_structured_mesh.py | 137 +++++++++++++++++++++++------ 2 files changed, 131 insertions(+), 28 deletions(-) diff --git a/geometries/README.md b/geometries/README.md index d3e7506..11fd8e6 100644 --- a/geometries/README.md +++ b/geometries/README.md @@ -69,7 +69,7 @@ python3 make_spe11c_geo.py --mesh-size 250 ``` -## Structured grid generation (`make_structured_mesh.py`) +## Structured grid generation (`make_structured_mesh.py -nx [nb_block_x] -ny [nb_block_y] [-nz [nb_block_z]]`) __Important__: for the variants A & B, keep in mind that the geometries are defined in the x-y plane instead of the x-z plane used in the description. See the above sections for more details. @@ -89,3 +89,23 @@ python3 make_structured_mesh.py --variant C -nx 200 -ny 200 -nz 200 Note that this script also requires the Python API of `gmsh`. Furthermore, passing the flag `--remove-cells-in-seal` creates mesh files in which the cells in the seal layers are removed. + + +### Adapted structured grid generation (`make_structured_mesh.py -rax [percent_of_domain_x nb_block_local_x ...] -ray [percent_of_domain_y nb_block_local_y ...]`) + +In addition to regular cartesian grid, an option (`--range-x or -rax`, `--range-y or -ray`) is added to generated block wise +refined grids. To be agnostic of real case dimensions, the upper bound of the block to be discretized is specified as a percent +of the total length (either _Lx_ or _Ly_ for `-rax` or `-ray`). It is followed by the number of cell along this axe for this block. + +The following command then generate a 300x200 mesh for case SPE11 variant A with 100 cells in x on the first half of _Lx_ (for variant A, _i.e._ _dx_ = 0.0125m) +and 200 cells in x on the second half of the domain (for variant A, _i.e._ _dx_ = 0.006125m). Along the y-axis, a regular spacing resulting in 200 cells is input. + +```bash +# generate a [100 200]x200 mesh for variant A in the file spe11a_structured.msh +python3 make_structured_mesh.py --variant A -rax .5 100 1. 200 -ray 1. 200 +``` + +__Note__ : `-rax` and `-ray` must be specified and it is not possible to mix and match with regular options `-nx,-ny,-nz`. +__Note__ : There is no such adaptation on the z-axis currently, _i.e._ SPE11 variant C is not accessible through this option. + + diff --git a/geometries/make_structured_mesh.py b/geometries/make_structured_mesh.py index 1cc477e..2ac4d63 100644 --- a/geometries/make_structured_mesh.py +++ b/geometries/make_structured_mesh.py @@ -16,15 +16,28 @@ import itertools import gmsh +import numpy as np from make_spe11c_geo import z_offset_at - +TOP_INDEX_SEAL = 8 PHYSICAL_INDEX_SEAL = 7 PHYSICAL_INDEX_OUTSIDE_OF_DOMAIN = 1000 PHYSICAL_NAME_OUTSIDE_OF_DOMAIN = str(PHYSICAL_INDEX_OUTSIDE_OF_DOMAIN) +def _merge_duplicates(vreal): + """ + unique vector of real normalized up to a tolerance + :param vreal: vector of non dimensional disctreization ticks between 0 and 1 + :return: vector with filtered out duplicates + """ + eps = 1e-5 + dx = np.abs(np.diff(vreal / np.max(vreal))) + i = np.where(dx < eps) + return np.delete(vreal, i) + + def _is_in_bbox(position, min, max) -> bool: return all(position[i] <= max[i] and position[i] >= min[i] for i in range(dim)) @@ -39,16 +52,39 @@ def _get_variant_geo_file(variant: str) -> str: class StructuredLattice: - def __init__(self, min: tuple, max: tuple, num_cells: tuple) -> None: + def __init__(self, min: tuple, max: tuple, num_cells: tuple, mask: list) -> None: + + if len(mask): + self._init_adapted(min, max, num_cells, mask) + else: + self._init_regular(min, max, num_cells) + + def _init_regular(self, min: tuple, max: tuple, num_cells: tuple) -> None: + self._num_cells = num_cells + self._dim = len(num_cells) + self._dx = [(max[i] - min[i]) / float(num_cells[i]) for i in range(self._dim)] + self._points = [ + ( + min[0] + i * self._dx[0], + min[1] + j * self._dx[1], + min[2] + k * (self._dx[2] if self._dim == 3 else 0.0) + ) + for k in range(num_cells[2] + 1 if self._dim == 3 else 1) + for j in range(num_cells[1] + 1) + for i in range(num_cells[0] + 1) + ] + + def _init_adapted(self, _min: tuple, _max: tuple, num_cells: tuple, mask: list) -> None: self._num_cells = num_cells self._dim = len(num_cells) - self._dx = [(max[i] - min[i])/float(num_cells[i]) for i in range(self._dim)] + self._dx = [[(_max[i] - _min[i]) * mask[i][j] for j in range(len(mask[i]))] for i in range(self._dim)] self._points = [ ( - min[0] + i*self._dx[0], - min[1] + j*self._dx[1], - min[2] + k*(self._dx[2] if self._dim == 3 else 0.0) + _min[0] + self._dx[0][i], + _min[1] + self._dx[1][j], + _min[2] + (self._dx[2][k] if self._dim == 3 else 0.0) ) + for k in range(num_cells[2] + 1 if self._dim == 3 else 1) for j in range(num_cells[1] + 1) for i in range(num_cells[0] + 1) @@ -56,17 +92,17 @@ def __init__(self, min: tuple, max: tuple, num_cells: tuple) -> None: @property def number_of_points(self) -> int: - points = (self._num_cells[0]+1)*(self._num_cells[1]+1) + points = (self._num_cells[0] + 1) * (self._num_cells[1] + 1) if self._dim == 2: return points - return points*(self._num_cells[2]+1) + return points * (self._num_cells[2] + 1) @property def number_of_cells(self) -> int: - cells = self._num_cells[0]*self._num_cells[1] + cells = self._num_cells[0] * self._num_cells[1] if self._dim == 2: return cells - return cells*self._num_cells[2] + return cells * self._num_cells[2] @property def points(self) -> list: @@ -91,13 +127,13 @@ def _get_quad_corners(p0: int) -> tuple: ) if self._dim == 2: - p0 = cell[1]*(num_cells[0] + 1) + cell[0] + p0 = cell[1] * (num_cells[0] + 1) + cell[0] return _get_quad_corners(p0) x, y, z = cell nx, ny = num_cells[0], num_cells[1] - z_layer_offset = (nx+1)*(ny+1) - p0 = z*z_layer_offset + y*(nx + 1) + x + z_layer_offset = (nx + 1) * (ny + 1) + p0 = z * z_layer_offset + y * (nx + 1) + x return _get_quad_corners(p0) + _get_quad_corners(p0 + z_layer_offset) def center(self, cell: tuple) -> tuple: @@ -108,7 +144,7 @@ def center(self, cell: tuple) -> tuple: result[i] + self._points[pidx][i] for i in range(self._dim) ]) - return tuple([(result[i]/len(corners) if i < self._dim else 0.0) for i in range(3)]) + return tuple([(result[i] / len(corners) if i < self._dim else 0.0) for i in range(3)]) class PhysicalIndexMapper: @@ -160,7 +196,6 @@ def _with_model_for_physical_index_queries(self, action): gmsh.model.setCurrent(self._model_name) return result - def _project_to_model_for_index_queries(self, position: tuple) -> tuple: if self._variant == "C": return ( @@ -246,7 +281,6 @@ def _is_included(self, cell_index: int) -> bool: return self._physical_cell_indices[cell_index] != self._exclude_physical_index - parser = argparse.ArgumentParser(description="Create a structured gmsh grid for one of the SPE11 variants") parser.add_argument( "-v", "--variant", @@ -260,9 +294,16 @@ def _is_included(self, cell_index: int) -> bool: action="store_true", help="Remove all cells within the seal layers" ) -parser.add_argument("-nx", "--number-of-cells-x", required=True, help="Desired number of cells in x-direction") -parser.add_argument("-ny", "--number-of-cells-y", required=True, help="Desired number of cells in y-direction") + +parser.add_argument("-nx", "--number-of-cells-x", required=False, help="Desired number of cells in x-direction") +parser.add_argument("-ny", "--number-of-cells-y", required=False, help="Desired number of cells in y-direction") parser.add_argument("-nz", "--number-of-cells-z", required=False, help="Desired number of cells in z-direction") + +parser.add_argument("-rax", "--range_x", required=False, nargs='+', type=float, + help="Desired mask in x for non uniform spacing in a sequence of [ percent-of-domain desired_num_cell_x ... ]") +parser.add_argument("-ray", "--range_y", required=False, nargs='+', type=float, + help="Desired mask in y for non uniform spacing in a sequence of [ percent-of-domain desired_num_cell_y ... ]") + args = vars(parser.parse_args()) variant = args["variant"] @@ -275,18 +316,60 @@ def _is_included(self, cell_index: int) -> bool: subprocess.run(["python3", "make_spe11c_geo.py", "-s", "100"], check=True) assert os.path.exists(_get_variant_geo_file(variant)) -num_cells = tuple([int(args[f"number_of_cells_{['x', 'y', 'z'][i]}"]) for i in range(dim)]) +if (args['number_of_cells_x'] is not None and args['number_of_cells_x'] is not None) and ( + args['range_x'] is None and args['range_y'] is None): + num_cells = tuple([int(args[f"number_of_cells_{['x', 'y', 'z'][i]}"]) for i in range(dim)]) +elif (args['number_of_cells_x'] is None and args['number_of_cells_x'] is None) and ( + args['range_x'] is not None and args['range_y'] is not None): + eps = .00001 + for i in range(0, len(args['range_x']), 2): + if i == 0: + _tmp = np.arange(0, args['range_x'][0] + eps, args['range_x'][0] / args['range_x'][1]) + else: + _step = (args['range_x'][i] - args['range_x'][i - 2]) / args['range_x'][i + 1] + _tmp = np.concatenate((_tmp, np.arange(_tmp[-1] + _step, args['range_x'][i] + eps, _step))) + + print(_tmp[-1], " -- ", args['range_x'][i]) + # assert (np.abs(_tmp[-1] - args['range_x'][i]) < .01 * eps) + + for i in range(0, len(args['range_y']), 2): + if i == 0: + _tmpy = np.arange(0, args['range_y'][0] + eps, args['range_y'][0] / args['range_y'][1]) + else: + _step = (args['range_y'][i] - args['range_y'][i - 2]) / args['range_y'][i + 1] + _tmpy = np.concatenate((_tmp, np.arange(_tmp[-1] + _step, args['range_y'][i] + eps, _step))) + + print(_tmpy[-1], " -- ", args['range_y'][i]) + # assert (np.abs(_tmpy[-1] - args['range_y'][i]) < .01 * eps) + + mask = [_merge_duplicates(_tmp), _merge_duplicates(_tmpy), [0.]] + num_cells = (len(mask[0]) - 1, len(mask[1]) - 1) +else: + raise NotImplementedError + +# TODO some asserts + gmsh_cell_type = ( 3 if dim == 2 # quadrangle, - else 5 # hexahedron + else 5 # hexahedron ) physical_index_mapper = PhysicalIndexMapper(variant) -lattice: StructuredLattice | FilteredLattice = StructuredLattice( - *_get_bounding_box(gmsh.model), - num_cells=num_cells -) + +if (args['range_x'] and args['range_y']): + lattice: StructuredLattice | FilteredLattice = StructuredLattice( + *_get_bounding_box(gmsh.model), + num_cells=num_cells, + mask=mask + ) +else: + lattice: StructuredLattice | FilteredLattice = StructuredLattice( + *_get_bounding_box(gmsh.model), + num_cells=num_cells, + mask=[] + ) + num_cells_total = lattice.number_of_cells print("Determining physical groups for all cells") @@ -319,7 +402,7 @@ def _is_included(self, cell_index: int) -> bool: """.format(len(physical_index_mapper.physical_groups(dim))).lstrip("\n"))) msh_file.write("{}".format( "\n".join(f'{dim} {tag} "{name}"' - for name, tag in physical_index_mapper.physical_groups(dim).items()) + for name, tag in physical_index_mapper.physical_groups(dim).items()) )) msh_file.write(textwrap.dedent(f""" $EndPhysicalNames @@ -327,7 +410,7 @@ def _is_included(self, cell_index: int) -> bool: {filtered_lattice.number_of_points} """)) for count, p in enumerate(filtered_lattice.points): - msh_file.write(f"{count+1} {' '.join(str(c) for c in p)}\n".format()) + msh_file.write(f"{count + 1} {' '.join(str(c) for c in p)}\n".format()) msh_file.write(textwrap.dedent(f""" $EndNodes $Elements @@ -335,7 +418,7 @@ def _is_included(self, cell_index: int) -> bool: """).lstrip("\n")) for cell_index, cell in enumerate(filtered_lattice.cells): phys_index = filtered_physical_cell_indices[cell_index] - msh_file.write(f"{cell_index+1} {gmsh_cell_type} 2 {phys_index} {phys_index} ") + msh_file.write(f"{cell_index + 1} {gmsh_cell_type} 2 {phys_index} {phys_index} ") msh_file.write(" ".join(str(i + 1) for i in filtered_lattice.corners(cell))) msh_file.write("\n") msh_file.write("\n$EndElements") From bfca272d0d65fd989a5f04a01d2c1289003c1856 Mon Sep 17 00:00:00 2001 From: jfranc Date: Mon, 16 Oct 2023 09:21:53 -0700 Subject: [PATCH 2/2] ready for first review --- geometries/README.md | 1 + geometries/make_structured_mesh.py | 12 ++++++------ 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/geometries/README.md b/geometries/README.md index e1fe964..7329d76 100644 --- a/geometries/README.md +++ b/geometries/README.md @@ -102,6 +102,7 @@ and 200 cells in x on the second half of the domain (for variant A, _i.e._ _dx_ ```bash # generate a [100 200]x200 mesh for variant A in the file spe11a_structured.msh python3 make_structured_mesh.py --variant A -rax .5 100 1. 200 -ray 1. 200 +python3 make_structured_mesh.py --variant B -rax 0.1 10 0.5 10 1. 15 -ray .1 15 .6 20 .8 10 1. 20 ``` __Note__ : `-rax` and `-ray` must be specified and it is not possible to mix and match with regular options `-nx,-ny,-nz`. diff --git a/geometries/make_structured_mesh.py b/geometries/make_structured_mesh.py index 2ac4d63..70fd152 100644 --- a/geometries/make_structured_mesh.py +++ b/geometries/make_structured_mesh.py @@ -29,7 +29,7 @@ def _merge_duplicates(vreal): """ unique vector of real normalized up to a tolerance - :param vreal: vector of non dimensional disctreization ticks between 0 and 1 + :param vreal: vector of non dimensional discreization ticks between 0 and 1 :return: vector with filtered out duplicates """ eps = 1e-5 @@ -324,12 +324,12 @@ def _is_included(self, cell_index: int) -> bool: eps = .00001 for i in range(0, len(args['range_x']), 2): if i == 0: - _tmp = np.arange(0, args['range_x'][0] + eps, args['range_x'][0] / args['range_x'][1]) + _tmpx = np.arange(0, args['range_x'][0] + eps, args['range_x'][0] / args['range_x'][1]) else: _step = (args['range_x'][i] - args['range_x'][i - 2]) / args['range_x'][i + 1] - _tmp = np.concatenate((_tmp, np.arange(_tmp[-1] + _step, args['range_x'][i] + eps, _step))) + _tmpx = np.concatenate((_tmpx, np.arange(_tmpx[-1] + _step, args['range_x'][i] + eps, _step))) - print(_tmp[-1], " -- ", args['range_x'][i]) + print(_tmpx[-1], " -- ", args['range_x'][i]) # assert (np.abs(_tmp[-1] - args['range_x'][i]) < .01 * eps) for i in range(0, len(args['range_y']), 2): @@ -337,12 +337,12 @@ def _is_included(self, cell_index: int) -> bool: _tmpy = np.arange(0, args['range_y'][0] + eps, args['range_y'][0] / args['range_y'][1]) else: _step = (args['range_y'][i] - args['range_y'][i - 2]) / args['range_y'][i + 1] - _tmpy = np.concatenate((_tmp, np.arange(_tmp[-1] + _step, args['range_y'][i] + eps, _step))) + _tmpy = np.concatenate((_tmpy, np.arange(_tmpy[-1] + _step, args['range_y'][i] + eps, _step))) print(_tmpy[-1], " -- ", args['range_y'][i]) # assert (np.abs(_tmpy[-1] - args['range_y'][i]) < .01 * eps) - mask = [_merge_duplicates(_tmp), _merge_duplicates(_tmpy), [0.]] + mask = [_merge_duplicates(_tmpx), _merge_duplicates(_tmpy), [0.]] num_cells = (len(mask[0]) - 1, len(mask[1]) - 1) else: raise NotImplementedError