From cd697b376d7875fc174f2a978de9ea7cf7280c98 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sat, 21 Jun 2014 13:40:50 -0500
Subject: [PATCH] Boundary finding appears to work
---
meshmode/discretization/connection.py | 119 ++++++++++++++++++++++-
meshmode/discretization/visualization.py | 17 +++-
meshmode/mesh/__init__.py | 46 ++++++++-
meshmode/mesh/io.py | 2 +-
4 files changed, 179 insertions(+), 5 deletions(-)
diff --git a/meshmode/discretization/connection.py b/meshmode/discretization/connection.py
index b20b267..4c251aa 100644
--- a/meshmode/discretization/connection.py
+++ b/meshmode/discretization/connection.py
@@ -23,6 +23,7 @@ THE SOFTWARE.
"""
import numpy as np
+import modepy as mp
import pyopencl as cl
import pyopencl.array # noqa
from pytools import memoize_method, memoize_method_nested
@@ -188,11 +189,127 @@ def make_same_mesh_connection(queue, to_discr, from_discr):
from_discr, to_discr, groups)
-def make_boundary_extractor(queue, discr):
+def make_boundary_extractor(queue, discr, group_factory):
"""
:return: a tuple ``(bdry_mesh, bdry_discr, connection)``
"""
+ # {{{ build face_map
+
+ # maps (igrp, el_grp, face_id) to a frozenset of vertex IDs
+ face_map = {}
+
+ for igrp, mgrp in enumerate(discr.mesh.groups):
+ grp_face_vertex_indices = mgrp.face_vertex_indices()
+ for iel_grp in xrange(mgrp.nelements):
+ for fid, loc_face_vertices in enumerate(grp_face_vertex_indices):
+ face_vertices = frozenset(
+ mgrp.vertex_indices[iel_grp, fvi]
+ for fvi in loc_face_vertices
+ )
+ face_map.setdefault(face_vertices, []).append(
+ (igrp, iel_grp, fid))
+
+ del face_vertices
+
+ # }}}
+
+ boundary_faces = [
+ face_ids[0]
+ for face_vertices, face_ids in face_map.iteritems()
+ if len(face_ids) == 1]
+
+ from pytools import flatten
+ bdry_vertex_vol_nrs = sorted(set(flatten(face_map.iterkeys())))
+
+ vol_to_bdry_vertices = np.empty(
+ discr.mesh.vertices.shape[-1],
+ discr.mesh.vertices.dtype)
+ vol_to_bdry_vertices.fill(-1)
+ vol_to_bdry_vertices[bdry_vertex_vol_nrs] = np.arange(
+ len(bdry_vertex_vol_nrs))
+
+ bdry_vertices = discr.mesh.vertices[:, bdry_vertex_vol_nrs]
+
+ from meshmode.mesh import Mesh, SimplexElementGroup
+ bdry_mesh_groups = []
+ for igrp, grp in enumerate(discr.groups):
+ mgrp = grp.mesh_el_group
+ group_boundary_faces = [
+ (ibface_group, ibface_el, ibface_face)
+ for ibface_group, ibface_el, ibface_face in boundary_faces
+ if ibface_group == igrp]
+
+ if not isinstance(mgrp, SimplexElementGroup):
+ raise NotImplementedError("can only take boundary of "
+ "SimplexElementGroup-based meshes")
+
+ # {{{ Preallocate arrays for mesh group
+
+ ngroup_bdry_elements = len(group_boundary_faces)
+ vertex_indices = np.empty(
+ (ngroup_bdry_elements, mgrp.dim+1-1),
+ mgrp.vertex_indices.dtype)
+
+ bdry_unit_nodes = mp.warp_and_blend_nodes(mgrp.dim-1, mgrp.order)
+ bdry_unit_nodes_01 = (bdry_unit_nodes + 1)*0.5
+
+ vol_basis = mp.simplex_onb(mgrp.dim, mgrp.order)
+ nbdry_unit_nodes = bdry_unit_nodes_01.shape[-1]
+ nodes = np.empty(
+ (discr.ambient_dim, ngroup_bdry_elements, nbdry_unit_nodes),
+ dtype=np.float64)
+
+ # }}}
+
+ grp_face_vertex_indices = mgrp.face_vertex_indices()
+ grp_vertex_unit_coordinates = mgrp.vertex_unit_coordinates()
+
+ for ibdry_el, (ibface_group, ibface_el, ibface_face) in enumerate(
+ group_boundary_faces):
+
+ # Find boundary vertex indices
+ loc_face_vertices = list(grp_face_vertex_indices[ibface_face])
+ glob_face_vertices = mgrp.vertex_indices[ibface_el, loc_face_vertices]
+ vertex_indices[ibdry_el] = vol_to_bdry_vertices[glob_face_vertices]
+
+ # Find unit nodes for boundary element
+ face_vertex_unit_coordinates = \
+ grp_vertex_unit_coordinates[loc_face_vertices]
+
+ # Find A, b such that A [e_1 e_2] + b = [r_1 r_2]
+ # (Notation assumes that the volume is 3D and the face is 2D.)
+ b = face_vertex_unit_coordinates[0]
+ A = (
+ face_vertex_unit_coordinates[1:]
+ - face_vertex_unit_coordinates[0]).T
+
+ face_unit_nodes = (np.dot(A, bdry_unit_nodes_01).T + b).T
+
+ resampling_mat = mp.resampling_matrix(
+ vol_basis,
+ face_unit_nodes, mgrp.unit_nodes)
+
+ nodes[:, ibdry_el, :] = np.einsum(
+ "ij,dj->di",
+ resampling_mat,
+ mgrp.nodes[:, ibface_el, :])
+
+ bdry_mesh_group = SimplexElementGroup(
+ mgrp.order, vertex_indices, nodes, unit_nodes=bdry_unit_nodes)
+ bdry_mesh_groups.append(bdry_mesh_group)
+
+ bdry_mesh = Mesh(bdry_vertices, bdry_mesh_groups)
+
+ from meshmode.discretization import Discretization
+ bdry_discr = Discretization(
+ discr.cl_context, bdry_mesh, group_factory)
+
+ # FIXME
+ connection = None
+
+ return bdry_mesh, bdry_discr, connection
+
# }}}
# vim: foldmethod=marker
diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py
index 9bdf0e0..815df95 100644
--- a/meshmode/discretization/visualization.py
+++ b/meshmode/discretization/visualization.py
@@ -131,7 +131,22 @@ class Visualizer(object):
vis_connectivity = self._vis_connectivity()
- if self.vis_discr.dim == 2:
+ if self.vis_discr.dim == 1:
+ nodes = list(nodes)
+ # pad to 3D with zeros
+ while len(nodes) < 3:
+ nodes.append(0*nodes[0])
+ assert len(nodes) == 3
+
+ args = tuple(nodes) + (field,)
+
+ # http://docs.enthought.com/mayavi/mayavi/auto/example_plotting_many_lines.html # noqa
+ src = mlab.pipeline.scalar_scatter(*args)
+ src.mlab_source.dataset.lines = vis_connectivity.reshape(-1, 2)
+ lines = mlab.pipeline.stripper(src)
+ mlab.pipeline.surface(lines, **kwargs)
+
+ elif self.vis_discr.dim == 2:
nodes = list(nodes)
# pad to 3D with zeros
while len(nodes) < 3:
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 3ff935a..91f3a7b 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -22,8 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
-
-#import numpy as np
+import numpy as np
import modepy as mp
#import numpy.linalg as la
from pytools import Record
@@ -130,6 +129,9 @@ class SimplexElementGroup(MeshElementGroup):
automatically assigned.
"""
+ if not issubclass(vertex_indices.dtype.type, np.integer):
+ raise TypeError("vertex_indices must be integral")
+
if unit_nodes is None:
if dim is None:
raise TypeError("'dim' must be passed "
@@ -147,6 +149,46 @@ class SimplexElementGroup(MeshElementGroup):
MeshElementGroup.__init__(self, order, vertex_indices, nodes,
element_nr_base, node_nr_base, unit_nodes, dim)
+ def face_vertex_indices(self):
+ if self.dim == 1:
+ return ((0, 1))
+ elif self.dim == 2:
+ return (
+ (0, 1),
+ (1, 2),
+ (0, 2),
+ )
+ elif self.dim == 3:
+ return (
+ (0, 1, 2),
+ (0, 1, 3),
+ (0, 2, 3),
+ (1, 2, 3)
+ )
+ else:
+ raise NotImplementedError("dim=%d" % self.dim)
+
+ def vertex_unit_coordinates(self):
+ if self.dim == 1:
+ return np.array([
+ [-1], [1]
+ ], dtype=np.float64)
+ elif self.dim == 2:
+ return np.array([
+ [-1, -1],
+ [1, -1],
+ [-1, 1],
+ ], dtype=np.float64)
+ elif self.dim == 3:
+ return np.array([
+ [-1, -1, -1],
+ [1, -1, -1],
+ [-1, 1, -1],
+ [-1, -1, 1],
+ ], dtype=np.float64)
+ else:
+ raise NotImplementedError("dim=%d" % self.dim)
+
# }}}
diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py
index 160153c..af83e49 100644
--- a/meshmode/mesh/io.py
+++ b/meshmode/mesh/io.py
@@ -133,7 +133,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase):
el_vertex_count = group_el_type.vertex_count()
vertex_indices = np.empty(
(ngroup_elements, el_vertex_count),
- np.float64)
+ np.int32)
i = 0
for el_vertices, el_nodes, el_type in zip(
--
GitLab