From 94f5c62bf64ec42abbe0a5ac8465a6fffe049a47 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Thu, 27 Nov 2014 23:53:07 -0600
Subject: [PATCH] Add mesh merge and affine map
---
meshmode/discretization/__init__.py | 8 ++-
meshmode/mesh/processing.py | 80 +++++++++++++++++++++++++++++
setup.cfg | 2 +-
test/test_meshmode.py | 34 +++++++++++-
4 files changed, 120 insertions(+), 4 deletions(-)
diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py
index e7a09c2..6088753 100644
--- a/meshmode/discretization/__init__.py
+++ b/meshmode/discretization/__init__.py
@@ -241,10 +241,14 @@ class Discretization(object):
result[d, k, i] = \
sum(j, resampling_mat[i, j] * nodes[d, k, j])
""",
- name="nodes")
+ name="nodes",
+ default_offset=lp.auto)
knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
- return lp.tag_inames(knl, dict(k="g.0"))
+ knl = lp.tag_inames(knl, dict(k="g.0"))
+ knl = lp.tag_data_axes(knl, "result",
+ "stride:auto,stride:auto,stride:auto")
+ return knl
result = self.empty(self.real_dtype, extra_dims=(self.ambient_dim,))
diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index dfa4607..94cd8b0 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -34,6 +34,8 @@ __doc__ = """
.. autofunction:: find_volume_mesh_element_orientations
.. autofunction:: perform_flips
.. autofunction:: find_bounding_box
+.. autofunction:: merge_dijsoint_meshes
+.. autofunction:: affine_map
"""
@@ -215,4 +217,82 @@ def find_bounding_box(mesh):
# }}}
+
+# {{{ merging
+
+def merge_dijsoint_meshes(meshes, skip_tests=False):
+ if not meshes:
+ raise ValueError("must pass at least one mesh")
+
+ from pytools import is_single_valued
+ if not is_single_valued(mesh.ambient_dim for mesh in meshes):
+ raise ValueError("all meshes must share the same ambient dimension")
+
+ # {{{ assemble combined vertex array
+
+ ambient_dim = meshes[0].ambient_dim
+ nvertices = sum(
+ mesh.vertices.shape[-1]
+ for mesh in meshes)
+
+ vert_dtype = np.find_common_type(
+ [mesh.vertices.dtype for mesh in meshes],
+ [])
+ vertices = np.empty(
+ (ambient_dim, nvertices), vert_dtype)
+
+ current_vert_base = 0
+ vert_bases = []
+ for mesh in meshes:
+ mesh_nvert = mesh.vertices.shape[-1]
+ vertices[:, current_vert_base:current_vert_base+mesh_nvert] = \
+ mesh.vertices
+
+ vert_bases.append(current_vert_base)
+ current_vert_base += mesh_nvert
+
+ # }}}
+
+ # {{{ assemble new groups list
+
+ new_groups = []
+
+ for mesh, vert_base in zip(meshes, vert_bases):
+ for group in mesh.groups:
+ new_vertex_indices = group.vertex_indices + vert_base
+ new_group = group.copy(vertex_indices=new_vertex_indices)
+ new_groups.append(new_group)
+
+ # }}}
+
+ from meshmode.mesh import Mesh
+ return Mesh(vertices, new_groups, skip_tests=skip_tests)
+
+# }}}
+
+
+# {{{ affine map
+
+def affine_map(mesh, A=None, b=None):
+ """Apply the affine map *f(x)=Ax+b* to the geometry of *mesh*."""
+
+ vertices = np.einsum("ij,jv->iv", A, mesh.vertices) + b[:, np.newaxis]
+
+ # {{{ assemble new groups list
+
+ new_groups = []
+
+ for group in mesh.groups:
+ nodes = (
+ np.einsum("ij,jen->ien", A, group.nodes)
+ + b[:, np.newaxis, np.newaxis])
+ new_groups.append(group.copy(nodes=nodes))
+
+ # }}}
+
+ from meshmode.mesh import Mesh
+ return Mesh(vertices, new_groups, skip_tests=True)
+
+# }}}
+
# vim: foldmethod=marker
diff --git a/setup.cfg b/setup.cfg
index eca7b43..6faef2e 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,3 +1,3 @@
[flake8]
-ignore = E126,E127,E128,E123,E226,E241,E242
+ignore = E126,E127,E128,E123,E226,E241,E242,E265
max-line-length=85
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 3c8f20d..0a24243 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -117,6 +117,37 @@ def test_element_orientation():
assert ((mesh_orient < 0) == (flippy > 0)).all()
+def test_merge_and_map(ctx_getter, visualize=False):
+ from meshmode.mesh.io import generate_gmsh, FileSource
+
+ mesh_order = 3
+
+ mesh = generate_gmsh(
+ FileSource("blob-2d.step"), 2, order=mesh_order,
+ force_ambient_dim=2,
+ other_options=["-string", "Mesh.CharacteristicLengthMax = 0.02;"]
+ )
+
+ from meshmode.mesh.processing import merge_dijsoint_meshes, affine_map
+ mesh2 = affine_map(mesh, A=np.eye(2), b=np.array([5, 0]))
+
+ mesh3 = merge_dijsoint_meshes((mesh2, mesh))
+
+ if visualize:
+ from meshmode.discretization import Discretization
+ from meshmode.discretization.poly_element import \
+ PolynomialWarpAndBlendGroupFactory
+ cl_ctx = ctx_getter()
+ queue = cl.CommandQueue(cl_ctx)
+
+ discr = Discretization(cl_ctx, mesh3,
+ PolynomialWarpAndBlendGroupFactory(3))
+
+ from meshmode.discretization.visualization import make_visualizer
+ vis = make_visualizer(queue, discr, 1)
+ vis.write_vtk_file("merged.vtu", [])
+
+
@pytest.mark.parametrize("dim", [2, 3])
@pytest.mark.parametrize("order", [1, 3])
def test_sanity_single_element(ctx_getter, dim, order, visualize=False):
@@ -242,7 +273,8 @@ def test_sanity_balls(ctx_getter, src_file, dim, mesh_order,
from meshmode.discretization.connection import make_boundary_restriction
bdry_mesh, bdry_discr, bdry_connection = make_boundary_restriction(
- queue, vol_discr, InterpolatoryQuadratureSimplexGroupFactory(quad_order))
+ queue, vol_discr,
+ InterpolatoryQuadratureSimplexGroupFactory(quad_order))
# }}}
--
GitLab