.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/mesh_generation/4_tree_mesh.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_mesh_generation_4_tree_mesh.py: Tree Meshes =========== Compared to tensor meshes, tree meshes are able to provide higher levels of discretization in certain regions while reducing the total number of cells. Tree meshes belong to the class (:class:`~discretize.TreeMesh`). Tree meshes can be defined in 2 or 3 dimensions. Here we demonstrate: - How to create basic tree meshes in 2D and 3D - Strategies for local mesh refinement - How to plot tree meshes - How to extract properties from tree meshes To create a tree mesh, we first define the base tensor mesh (a mesh comprised entirely of the smallest cells). Next we choose the level of discretization around certain points or within certain regions. When creating tree meshes, we must remember certain rules: - The number of base mesh cells in x, y and z must all be powers of 2 - We cannot refine the mesh to create cells smaller than those defining the base mesh - The range of cell sizes in the tree mesh depends on the number of base mesh cells in x, y and z .. GENERATED FROM PYTHON SOURCE LINES 28-33 Import Packages --------------- Here we import the packages required for this tutorial. .. GENERATED FROM PYTHON SOURCE LINES 34-42 .. code-block:: Python from discretize import TreeMesh from discretize.utils import mkvc import matplotlib.pyplot as plt import numpy as np # sphinx_gallery_thumbnail_number = 4 .. GENERATED FROM PYTHON SOURCE LINES 43-51 Basic Example ------------- Here we demonstrate the basic two step process for creating a 2D tree mesh (QuadTree mesh). The region of highest discretization if defined within a rectangular box. We use the keyword argument *octree_levels* to define the rate of cell width increase outside the box. .. GENERATED FROM PYTHON SOURCE LINES 51-71 .. code-block:: Python dh = 5 # minimum cell width (base mesh cell width) nbc = 64 # number of base mesh cells in x # Define base mesh (domain and finest discretization) h = dh * np.ones(nbc) mesh = TreeMesh([h, h]) # Define corner points for rectangular box xp, yp = np.meshgrid([120.0, 240.0], [80.0, 160.0]) xy = np.c_[mkvc(xp), mkvc(yp)] # mkvc creates vectors # Discretize to finest cell size within rectangular box, with padding in the z direction # at the finest and second finest levels. padding = [[0, 2], [0, 2]] mesh.refine_bounding_box(xy, level=-1, padding_cells_by_level=padding) mesh.plot_grid(show_it=True) .. image-sg:: /tutorials/mesh_generation/images/sphx_glr_4_tree_mesh_001.png :alt: 4 tree mesh :srcset: /tutorials/mesh_generation/images/sphx_glr_4_tree_mesh_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 72-84 Intermediate Example and Plotting --------------------------------- The widths of the base mesh cells do not need to be the same in x and y. However the number of base mesh cells in x and y each needs to be a power of 2. Here we show topography-based mesh refinement and refinement about a set of points. We also show some aspect of customizing plots. We use the keyword argument *octree_levels* to define the rate of cell width increase relative to our surface and the set of discrete points about which we are refining. .. GENERATED FROM PYTHON SOURCE LINES 84-123 .. code-block:: Python dx = 5 # minimum cell width (base mesh cell width) in x dy = 5 # minimum cell width (base mesh cell width) in y x_length = 300.0 # domain width in x y_length = 300.0 # domain width in y # Compute number of base mesh cells required in x and y nbcx = 2 ** int(np.round(np.log(x_length / dx) / np.log(2.0))) nbcy = 2 ** int(np.round(np.log(y_length / dy) / np.log(2.0))) # Define the base mesh hx = [(dx, nbcx)] hy = [(dy, nbcy)] mesh = TreeMesh([hx, hy], x0="CC") # Refine surface topography xx = mesh.nodes_x yy = -3 * np.exp((xx**2) / 100**2) + 50.0 pts = np.c_[mkvc(xx), mkvc(yy)] padding = [[0, 2], [0, 2]] mesh.refine_surface(pts, padding_cells_by_level=padding, finalize=False) # Refine mesh near points xx = np.array([0.0, 10.0, 0.0, -10.0]) yy = np.array([-20.0, -10.0, 0.0, -10]) pts = np.c_[mkvc(xx), mkvc(yy)] mesh.refine_points(pts, padding_cells_by_level=[2, 2], finalize=False) mesh.finalize() # We can apply the plot_grid method and output to a specified axes object fig = plt.figure(figsize=(6, 6)) ax = fig.add_subplot(111) mesh.plot_grid(ax=ax) ax.set_xbound(mesh.x0[0], mesh.x0[0] + np.sum(mesh.h[0])) ax.set_ybound(mesh.x0[1], mesh.x0[1] + np.sum(mesh.h[1])) ax.set_title("QuadTree Mesh") .. image-sg:: /tutorials/mesh_generation/images/sphx_glr_4_tree_mesh_002.png :alt: QuadTree Mesh :srcset: /tutorials/mesh_generation/images/sphx_glr_4_tree_mesh_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'QuadTree Mesh') .. GENERATED FROM PYTHON SOURCE LINES 124-130 Extracting Mesh Properties -------------------------- Once the mesh is created, you may want to extract certain properties. Here, we show some properties that can be extracted from a QuadTree mesh. .. GENERATED FROM PYTHON SOURCE LINES 130-183 .. code-block:: Python dx = 5 # minimum cell width (base mesh cell width) in x dy = 5 # minimum cell width (base mesh cell width) in y x_length = 300.0 # domain width in x y_length = 300.0 # domain width in y # Compute number of base mesh cells required in x and y nbcx = 2 ** int(np.round(np.log(x_length / dx) / np.log(2.0))) nbcy = 2 ** int(np.round(np.log(y_length / dy) / np.log(2.0))) # Define the base mesh hx = [(dx, nbcx)] hy = [(dy, nbcy)] mesh = TreeMesh([hx, hy], x0="CC") # Refine surface topography xx = mesh.nodes_x yy = -3 * np.exp((xx**2) / 100**2) + 50.0 pts = np.c_[mkvc(xx), mkvc(yy)] padding = [[0, 2], [0, 2]] mesh.refine_surface(pts, padding_cells_by_level=padding, finalize=False) # Refine near points xx = np.array([0.0, 10.0, 0.0, -10.0]) yy = np.array([-20.0, -10.0, 0.0, -10]) pts = np.c_[mkvc(xx), mkvc(yy)] mesh.refine_points(pts, padding_cells_by_level=[2, 2], finalize=False) mesh.finalize() # The bottom west corner x0 = mesh.x0 # The total number of cells nC = mesh.nC # An (nC, 2) array containing the cell-center locations cc = mesh.gridCC # A boolean array specifying which cells lie on the boundary bInd = mesh.cell_boundary_indices # The cell areas (2D "volume") s = mesh.cell_volumes fig = plt.figure(figsize=(6, 6)) ax = fig.add_subplot(111) mesh.plot_image(np.log10(s), grid=True, ax=ax) ax.set_xbound(mesh.x0[0], mesh.x0[0] + np.sum(mesh.h[0])) ax.set_ybound(mesh.x0[1], mesh.x0[1] + np.sum(mesh.h[1])) ax.set_title("Log of Cell Areas") .. image-sg:: /tutorials/mesh_generation/images/sphx_glr_4_tree_mesh_003.png :alt: Log of Cell Areas :srcset: /tutorials/mesh_generation/images/sphx_glr_4_tree_mesh_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Log of Cell Areas') .. GENERATED FROM PYTHON SOURCE LINES 184-190 3D Example ---------- Here we show how the same approach can be used to create and extract properties from a 3D tree mesh. .. GENERATED FROM PYTHON SOURCE LINES 190-242 .. code-block:: Python dx = 5 # minimum cell width (base mesh cell width) in x dy = 5 # minimum cell width (base mesh cell width) in y dz = 5 # minimum cell width (base mesh cell width) in z x_length = 300.0 # domain width in x y_length = 300.0 # domain width in y z_length = 300.0 # domain width in y # Compute number of base mesh cells required in x and y nbcx = 2 ** int(np.round(np.log(x_length / dx) / np.log(2.0))) nbcy = 2 ** int(np.round(np.log(y_length / dy) / np.log(2.0))) nbcz = 2 ** int(np.round(np.log(z_length / dz) / np.log(2.0))) # Define the base mesh hx = [(dx, nbcx)] hy = [(dy, nbcy)] hz = [(dz, nbcz)] mesh = TreeMesh([hx, hy, hz], x0="CCC") # Refine surface topography [xx, yy] = np.meshgrid(mesh.nodes_x, mesh.nodes_y) zz = -3 * np.exp((xx**2 + yy**2) / 100**2) + 50.0 pts = np.c_[mkvc(xx), mkvc(yy), mkvc(zz)] padding = [[0, 0, 2], [0, 0, 2]] mesh.refine_surface(pts, padding_cells_by_level=padding, finalize=False) # Refine box xp, yp, zp = np.meshgrid([-40.0, 40.0], [-40.0, 40.0], [-60.0, 0.0]) xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)] mesh.refine_bounding_box(xyz, padding_cells_by_level=padding, finalize=False) mesh.finalize() # The bottom west corner x0 = mesh.x0 # The total number of cells nC = mesh.nC # An (nC, 3) array containing the cell-center locations cc = mesh.gridCC # A boolean array specifying which cells lie on the boundary bInd = mesh.cell_boundary_indices # Cell volumes v = mesh.cell_volumes fig = plt.figure(figsize=(6, 6)) ax = fig.add_subplot(111) mesh.plot_slice(np.log10(v), normal="Y", ax=ax, ind=int(mesh.h[1].size / 2), grid=True) ax.set_title("Cell Log-Volumes at Y = 0 m") .. image-sg:: /tutorials/mesh_generation/images/sphx_glr_4_tree_mesh_004.png :alt: Cell Log-Volumes at Y = 0 m :srcset: /tutorials/mesh_generation/images/sphx_glr_4_tree_mesh_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Cell Log-Volumes at Y = 0 m') .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.810 seconds) .. _sphx_glr_download_tutorials_mesh_generation_4_tree_mesh.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 4_tree_mesh.ipynb <4_tree_mesh.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 4_tree_mesh.py <4_tree_mesh.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 4_tree_mesh.zip <4_tree_mesh.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_