From 13456597bfcd03294c80c3de7e18a4e4cfaa1e2e Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sat, 11 Apr 2015 19:12:53 -0500
Subject: [PATCH] Add generate_box_mesh
---
meshmode/mesh/generation.py | 131 +++++++++++++++++++++++++++---------
test/test_meshmode.py | 5 ++
2 files changed, 105 insertions(+), 31 deletions(-)
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 5b925dcb..9a72c537 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -53,6 +53,11 @@ Surfaces
.. autofunction:: generate_icosahedron
.. autofunction:: generate_icosphere
.. autofunction:: generate_torus
+
+Volumes
+-------
+
+.. autofunction:: generate_box_mesh
.. autofunction:: generate_regular_rect_mesh
"""
@@ -373,51 +378,95 @@ def generate_torus(r_outer, r_inner, n_outer=20, n_inner=10, order=1):
return mesh
-# {{{ generate_regular_rect_mesh
+# {{{ generate_box_mesh
-def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1):
- """Create a semi-structured rectangular mesh.
+def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64):
+ """Create a semi-structured mesh.
- :param a: the lower left hand point of the rectangle
- :param b: the upper right hand point of the rectangle
- :param n: a tuple of integers indicating the total number of points
- on [a,b].
+ :param axis_coords: a tuple with a number of entries corresponding
+ to the number of dimensions, with each entry a numpy array
+ specifying the coordinates to be used along that axis.
"""
- if min(n) < 2:
- raise ValueError("need at least two points in each direction")
- node_dict = {}
- vertices = []
- points_1d = [np.linspace(a_i, b_i, n_i)
- for a_i, b_i, n_i in zip(a, b, n)]
+ for iaxis, axc in enumerate(axis_coords):
+ if len(axc) < 2:
+ raise ValueError("need at least two points along axis %d"
+ % (iaxis+1))
- for j in range(n[1]):
- for i in range(n[0]):
- node_dict[i, j] = len(vertices)
- vertices.append(np.array([points_1d[0][i], points_1d[1][j]]))
+ dim = len(axis_coords)
- vertices = np.array(vertices).T.copy()
+ shape = tuple(len(axc) for axc in axis_coords)
- vertex_indices = []
+ from pytools import product
+ nvertices = product(shape)
- for i in range(n[0]-1):
- for j in range(n[1]-1):
+ vertex_indices = np.arange(nvertices).reshape(*shape, order="F")
- # c--d
- # | |
+ vertices = np.empty((dim,)+shape, dtype=coord_dtype)
+ for idim in range(dim):
+ vshape = (shape[idim],) + (1,)*idim
+ vertices[idim] = axis_coords[idim].reshape(*vshape)
+
+ vertices = vertices.reshape(dim, -1)
+
+ el_vertices = []
+
+ if dim == 1:
+ for i in range(shape[0]-1):
# a--b
- a = node_dict[i, j]
- b = node_dict[i+1, j]
- c = node_dict[i, j+1]
- d = node_dict[i+1, j+1]
+ a = vertex_indices[i]
+ b = vertex_indices[i+1]
- vertex_indices.append((a, b, c))
- vertex_indices.append((d, c, b))
+ el_vertices.append((a, b,))
- vertex_indices = np.array(vertex_indices, dtype=np.int32)
+ elif dim == 2:
+ for i in range(shape[0]-1):
+ for j in range(shape[1]-1):
- grp = make_group_from_vertices(vertices, vertex_indices, order)
+ # c--d
+ # | |
+ # a--b
+
+ a = vertex_indices[i, j]
+ b = vertex_indices[i+1, j]
+ c = vertex_indices[i, j+1]
+ d = vertex_indices[i+1, j+1]
+
+ el_vertices.append((a, b, c))
+ el_vertices.append((d, c, b))
+
+ elif dim == 3:
+ for i in range(shape[0]-1):
+ for j in range(shape[1]-1):
+ for k in range(shape[2]-1):
+
+ a000 = vertex_indices[i, j, k]
+ a001 = vertex_indices[i, j, k+1]
+ a010 = vertex_indices[i, j+1, k]
+ a011 = vertex_indices[i, j+1, k+1]
+
+ a100 = vertex_indices[i+1, j, k]
+ a101 = vertex_indices[i+1, j, k+1]
+ a110 = vertex_indices[i+1, j+1, k]
+ a111 = vertex_indices[i+1, j+1, k+1]
+
+ el_vertices.append((a000, a100, a010, a001))
+ el_vertices.append((a101, a100, a001, a010))
+ el_vertices.append((a101, a011, a010, a001))
+
+ el_vertices.append((a100, a010, a101, a110))
+ el_vertices.append((a011, a010, a110, a101))
+ el_vertices.append((a011, a111, a101, a110))
+
+ else:
+ raise NotImplementedError("box meshes of dimension %d"
+ % dim)
+
+ el_vertices = np.array(el_vertices, dtype=np.int32)
+
+ grp = make_group_from_vertices(
+ vertices.reshape(dim, -1), el_vertices, order)
from meshmode.mesh import Mesh
return Mesh(vertices, [grp],
@@ -426,4 +475,24 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1):
# }}}
+# {{{ generate_regular_rect_mesh
+
+def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1):
+ """Create a semi-structured rectangular mesh.
+
+ :param a: the lower left hand point of the rectangle
+ :param b: the upper right hand point of the rectangle
+ :param n: a tuple of integers indicating the total number of points
+ on [a,b].
+ """
+ if min(n) < 2:
+ raise ValueError("need at least two points in each direction")
+
+ axis_coords = [np.linspace(a_i, b_i, n_i)
+ for a_i, b_i, n_i in zip(a, b, n)]
+
+ return generate_box_mesh(axis_coords, order=order)
+
+# }}}
+
# vim: fdm=marker
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 25b26dab..e7d3e941 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -354,6 +354,11 @@ def test_rect_mesh(do_plot=False):
pt.show()
+def test_box_mesh():
+ from meshmode.mesh.generation import generate_box_mesh
+ generate_box_mesh(3*(np.linspace(0, 1, 5),))
+
+
def test_as_python():
from meshmode.mesh.generation import make_curve_mesh, cloverleaf
mesh = make_curve_mesh(cloverleaf, np.linspace(0, 1, 100), order=3)
--
GitLab