diff --git a/doc/discretization.rst b/doc/discretization.rst index b9c6f6539185664877cdb0b7f07f322ab1ce69c6..719a6a713ca2f977853190de71125896f6fcf56c 100644 --- a/doc/discretization.rst +++ b/doc/discretization.rst @@ -12,6 +12,11 @@ Composite polynomial discretization .. automodule:: meshmode.discretization.poly_element Resampling ------------------------------------ +---------- .. automodule:: meshmode.discretization.resampling + +Visualization +------------- + +.. automodule:: meshmode.discretization.visualization diff --git a/doc/index.rst b/doc/index.rst index badd1c497f30dc690bae2260c5c4f421667c66c7..f652d1ba9f72dffe21584f07737ebf1a7b40c4a9 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -5,6 +5,7 @@ Contents: .. toctree:: :maxdepth: 2 + mesh discretization diff --git a/doc/mesh.rst b/doc/mesh.rst index 9a3d4722666bdb8341181d990fc4f389801a7969..a6b9a5bce75ee17f5b7a8dad96250e1ff5bdf9d6 100644 --- a/doc/mesh.rst +++ b/doc/mesh.rst @@ -16,4 +16,14 @@ Mesh generation .. automodule:: pytential.mesh.generation +Mesh input/output +----------------- + +.. automodule:: pytential.mesh.io + +Mesh connections (for interpolation) +------------------------------------ + +.. automodule:: pytential.mesh.connection + .. vim: sw=4 diff --git a/meshmode/discretization/connection.py b/meshmode/discretization/connection.py new file mode 100644 index 0000000000000000000000000000000000000000..ff5fb0212faf6c37e760314e9673ad65552c8ae3 --- /dev/null +++ b/meshmode/discretization/connection.py @@ -0,0 +1,144 @@ +from __future__ import division + +__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np +import pyopencl as cl +import pyopencl.array # noqa + + +__doc__ = """ +.. autoclass:: DiscretizationConnection + +.. autofunction:: make_same_mesh_connection + +Implementation details +^^^^^^^^^^^^^^^^^^^^^^ + +.. autoclass:: InterpolationGroup + +.. autoclass:: DiscretizationConnectionElementGroup +""" + + +class InterpolationGroup(object): + """ + .. attribute:: source_element_indices + + A :class:`numpy.ndarray` of length ``nelements``, containing the + element index from which this "*to*" element's data will be + interpolated. + + .. attribute:: target_element_indices + + A :class:`numpy.ndarray` of length ``nelements``, containing the + element index to which this "*to*" element's data will be + interpolated. + + .. attribute:: source_interpolation_nodes + + A :class:`numpy.ndarray` of shape + ``(from_group.dim,to_group.nelements,to_group.nunit_nodes)`` + storing the coordinates of the nodes (in unit coordinates + of the *from* reference element) from which the node + locations of this element should be interpolated. + """ + def __init__(self, source_element_indices, + target_element_indices, source_interpolation_nodes): + self.source_element_indices = source_element_indices + self.target_element_indices = target_element_indices + self.source_interpolation_nodes = source_interpolation_nodes + + @property + def nelements(self): + return len(self.source_element_indices) + + +class DiscretizationConnectionElementGroup(object): + """ + .. attribute:: interpolation_groups + + A list of :class:`InterpolationGroup` instances. + """ + def __init__(self, interpolation_groups): + self.interpolation_groups = interpolation_groups + + +class DiscretizationConnection(object): + """ + .. attribute:: from_discr + + .. attribute:: to_discr + + .. attribute:: groups + + a list of :class:`MeshConnectionGroup` instances, with + a one-to-one correspondence to the groups in + :attr:`from_discr` and :attr:`to_discr`. + """ + + def __init__(self, from_discr, to_discr, groups): + if from_discr.cl_context != to_discr.cl_context: + raise ValueError("from_discr and to_discr must live in the " + "same OpenCL context") + + self.cl_context = from_discr.cl_context + + self.from_discr = from_discr + self.to_discr = to_discr + self.groups = groups + + def __call__(self, queue, field): + pass + + +# {{{ constructor functions + +def make_same_mesh_connection(queue, from_discr, to_discr): + if from_discr.mesh is not to_discr.mesh: + raise ValueError("from_discr and to_discr must be based on " + "the same mesh") + + assert queue.context == from_discr.cl_context + assert queue.context == to_discr.cl_context + + groups = [] + for fgrp, tgrp in zip(from_discr.groups, to_discr.groups): + all_elements = cl.array.arange(queue, + fgrp.nelements, + dtype=np.intp) + igroup = InterpolationGroup( + source_element_indices=all_elements, + target_element_indices=all_elements, + source_interpolation_nodes=tgrp.unit_nodes) + + groups.append( + DiscretizationConnectionElementGroup( + igroup)) + + return DiscretizationConnection( + from_discr, to_discr, groups) + +# }}} + +# vim: foldmethod=marker diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index c8016f8755dad6a613807a0158a5c612ab1db58e..114d6ace4482aa0c3ea92cedeb94d0ddd739c4c1 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -107,16 +107,16 @@ class PolynomialElementGroupBase(object): # }}} -# {{{ element group +# {{{ concrete element groups -class PolynomialElementGroup(PolynomialElementGroupBase): +class PolynomialQuadratureElementGroup(PolynomialElementGroupBase): @memoize_method def _quadrature_rule(self): dims = self.mesh_el_group.dim if dims == 1: return mp.LegendreGaussQuadrature(self.order) else: - return mp.VioreanuRokhlinSimplexQuad(self.order, dims) + return mp.VioreanuRokhlinSimplexQuadrature(self.order, dims) @property @memoize_method @@ -128,12 +128,25 @@ class PolynomialElementGroup(PolynomialElementGroupBase): def weights(self): return self._quadrature_rule().weights + +class PolynomialWarpAndBlendElementGroup(PolynomialElementGroupBase): + @property + @memoize_method + def unit_nodes(self): + dims = self.mesh_el_group.dim + return mp.warp_and_blend_nodes(dims, self.order) + + @property + @memoize_method + def weights(self): + raise NotImplementedError() + # }}} # {{{ discretization -class PolynomialElementDiscretization(Discretization): +class PolynomialElementDiscretizationBase(Discretization): """An (unstructured) composite polynomial discretization without any specific opinion on how to evaluate layer potentials. @@ -280,9 +293,16 @@ class PolynomialElementDiscretization(Discretization): return result - group_class = PolynomialElementGroup - # }}} +class PolynomialQuadratureElementDiscretization( + PolynomialElementDiscretizationBase): + group_class = PolynomialQuadratureElementGroup + + +class PolynomialWarpAndBlendElementDiscretization( + PolynomialElementDiscretizationBase): + group_class = PolynomialWarpAndBlendElementGroup + # vim: fdm=marker diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py new file mode 100644 index 0000000000000000000000000000000000000000..fd928a4a992d3d996e61a4954f875f2454b069fe --- /dev/null +++ b/meshmode/discretization/visualization.py @@ -0,0 +1,204 @@ +from __future__ import division + +__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np +from pytools import memoize_method +import pyopencl as cl + +__doc__ = """ + +.. autofunction:: make_visualizer + +.. autoclass:: Visualizer +""" + + +def separate_by_real_and_imag(data, real_only): + for name, field in data: + from pytools.obj_array import log_shape + ls = log_shape(field) + + if ls != () and ls[0] > 1: + assert len(ls) == 1 + from pytools.obj_array import ( + oarray_real_copy, oarray_imag_copy, + with_object_array_or_scalar) + + if field[0].dtype.kind == "c": + if real_only: + yield (name, + with_object_array_or_scalar(oarray_real_copy, field)) + else: + yield (name+"_r", + with_object_array_or_scalar(oarray_real_copy, field)) + yield (name+"_i", + with_object_array_or_scalar(oarray_imag_copy, field)) + else: + yield (name, field) + else: + # ls == () + if field.dtype.kind == "c": + yield (name+"_r", field.real.copy()) + yield (name+"_i", field.imag.copy()) + else: + yield (name, field) + + +class Visualizer(object): + def __init__(self, discr, vis_discr, connection): + self.discr = discr + self.vis_discr = vis_discr + self.connection = connection + + @memoize_method + def _vis_connectivity(self): + """ + :return: an array of shape + ``(vis_discr.nelements,nsubelements,primitive_element_size)`` + """ + # Assume that we're using modepy's default node ordering. + + from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ + as gnitstam, single_valued + vis_order = single_valued( + group.order for group in self.vis_discr.groups) + node_tuples = list(gnitstam(vis_order, self.vis_discr.dim)) + + from modepy.tools import submesh + el_connectivity = np.array( + submesh(node_tuples), + dtype=np.intp) + + nelements = sum(group.nelements for group in self.vis_discr.groups) + vis_connectivity = np.empty( + (nelements,) + el_connectivity.shape, dtype=np.intp) + + el_nr_base = 0 + for group in self.vis_discr.groups: + assert len(node_tuples) == group.nunit_nodes + vis_connectivity[el_nr_base:el_nr_base+group.nelements] = ( + np.arange( + el_nr_base*group.nunit_nodes, + (el_nr_base+group.nelements)*group.nunit_nodes, + group.nunit_nodes + )[:, np.newaxis, np.newaxis] + + el_connectivity) + el_nr_base += group.nelements + + return vis_connectivity + + def show_scalar_in_mayavi(self, field, **kwargs): + import mayavi.mlab as mlab + + do_show = kwargs.pop("do_show", True) + + with cl.CommandQueue(self.vis_discr.cl_context) as queue: + nodes = self.vis_discr.nodes().with_queue(queue).get() + + assert nodes.shape[0] == self.vis_discr.ambient_dim + #mlab.points3d(nodes[0], nodes[1], 0*nodes[0]) + + vis_connectivity = self._vis_connectivity() + + if self.vis_discr.dim == 2: + nodes = list(nodes) + # pad to 3D with zeros + while len(nodes) < 3: + nodes.append(0*nodes[0]) + + args = tuple(nodes) + (vis_connectivity.reshape(-1, 3),) + + mlab.triangular_mesh(*args, **kwargs) + + else: + raise RuntimeError("meshes of bulk dimension %d are currently " + "unsupported" % self.vis_discr.dim) + + if do_show: + mlab.show() + + def write_vtk_file(self, file_name, fields_and_names, compressor=None, + real_only=False): + + from pyvisfile.vtk import ( + UnstructuredGrid, DataArray, + AppendedDataXMLGenerator, + VTK_LINE, VTK_TRIANGLE, VTK_TETRA, + VF_LIST_OF_COMPONENTS) + el_types = { + 1: VTK_LINE, + 2: VTK_TRIANGLE, + 3: VTK_TETRA, + } + + el_type = el_types[self.vis_discr.dim] + + with cl.CommandQueue(self.vis_discr.cl_context) as queue: + nodes = self.vis_discr.nodes().with_queue(queue).get() + + connectivity = self._vis_connectivity() + + nprimitive_elements = ( + connectivity.shape[0] + * connectivity.shape[1]) + + grid = UnstructuredGrid( + (self.vis_discr.nnodes, + DataArray("points", + nodes.reshape(self.vis_discr.ambient_dim, -1), + vector_format=VF_LIST_OF_COMPONENTS)), + cells=connectivity.reshape(-1), + cell_types=np.asarray([el_type] * nprimitive_elements, + dtype=np.uint8)) + + # for name, field in separate_by_real_and_imag(cell_data, real_only): + # grid.add_celldata(DataArray(name, field, + # vector_format=VF_LIST_OF_COMPONENTS)) + + for name, field in separate_by_real_and_imag(fields_and_names, real_only): + grid.add_pointdata(DataArray(name, field, + vector_format=VF_LIST_OF_COMPONENTS)) + + from os.path import exists + if exists(file_name): + raise RuntimeError("output file '%s' already exists" + % file_name) + + with open(file_name, "w") as outf: + AppendedDataXMLGenerator(compressor)(grid).write(outf) + + +def make_visualizer(queue, discr, vis_order): + from meshmode.discretization.poly_element import \ + PolynomialWarpAndBlendElementDiscretization + vis_discr = PolynomialWarpAndBlendElementDiscretization( + discr.cl_context, discr.mesh, order=vis_order, + real_dtype=discr.real_dtype) + from meshmode.discretization.connection import \ + make_same_mesh_connection + cnx = make_same_mesh_connection(queue, discr, vis_discr) + + return Visualizer(discr, vis_discr, cnx) + +# vim: foldmethod=marker diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 3434c4424685d0de0fd3015a58a342ad7e8e1d98..2ca445bc33618247010bd22b93908752f41be4cc 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -82,6 +82,10 @@ class MeshElementGroup(Record): """ if unit_nodes is None: + if dim is None: + raise TypeError("'dim' must be passed " + "if 'unit_nodes' is not passed") + unit_nodes = mp.warp_and_blend_nodes(dim, order) dims = unit_nodes.shape[0] diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py new file mode 100644 index 0000000000000000000000000000000000000000..c7758b61051b95c92924867073dd479d07ff782b --- /dev/null +++ b/meshmode/mesh/io.py @@ -0,0 +1,198 @@ +from __future__ import division + +__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np + +from meshpy.gmsh_reader import ( # noqa + GmshMeshReceiverBase, FileSource, LiteralSource) + + +__doc__ = """ + +.. autoclass:: FileSource +.. autoclass:: LiteralSource + +.. autofunction:: read_gmsh +.. autofunction:: generate_gmsh + +""" + + +# {{{ gmsh receiver + +class GmshMeshReceiver(GmshMeshReceiverBase): + def __init__(self): + # Use data fields similar to meshpy.triangle.MeshInfo and + # meshpy.tet.MeshInfo + self.points = None + self.elements = None + self.element_types = None + self.element_markers = None + self.tags = None + + def set_up_nodes(self, count): + # Preallocate array of nodes within list; treat None as sentinel value. + # Preallocation not done for performance, but to assign values at indices + # in random order. + self.points = [None] * count + + def add_node(self, node_nr, point): + self.points[node_nr] = point + + def finalize_nodes(self): + self.points = np.array(self.points, dtype=np.float64) + + def set_up_elements(self, count): + # Preallocation of arrays for assignment elements in random order. + self.element_vertices = [None] * count + self.element_nodes = [None] * count + self.element_types = [None] * count + self.element_markers = [None] * count + self.tags = [] + + def add_element(self, element_nr, element_type, vertex_nrs, + lexicographic_nodes, tag_numbers): + self.element_vertices[element_nr] = vertex_nrs + self.element_nodes[element_nr] = lexicographic_nodes + self.element_types[element_nr] = element_type + self.element_markers[element_nr] = tag_numbers + + def finalize_elements(self): + pass + + def add_tag(self, name, index, dimension): + pass + + def finalize_tags(self): + pass + + def get_mesh(self): + el_type_hist = {} + for el_type in self.element_types: + el_type_hist[el_type] = el_type_hist.get(el_type, 0) + 1 + + groups = self.groups = [] + + ambient_dim = self.points.shape[-1] + + mesh_bulk_dim = max( + el_type.dimensions for el_type in el_type_hist.iterkeys()) + + # {{{ build vertex numbering + + vertex_index_gmsh_to_mine = {} + for el_vertices, el_type in zip( + self.element_vertices, self.element_types): + for gmsh_vertex_nr in el_vertices: + if gmsh_vertex_nr not in vertex_index_gmsh_to_mine: + vertex_index_gmsh_to_mine[gmsh_vertex_nr] = \ + len(vertex_index_gmsh_to_mine) + + # }}} + + # {{{ build vertex array + + gmsh_vertex_indices, my_vertex_indices = \ + zip(*vertex_index_gmsh_to_mine.iteritems()) + vertices = np.empty( + (ambient_dim, len(vertex_index_gmsh_to_mine)), dtype=np.float64) + vertices[:, np.array(my_vertex_indices, np.intp)] = \ + self.points[np.array(gmsh_vertex_indices, np.intp)].T + + # }}} + + from meshmode.mesh import MeshElementGroup, Mesh + + for group_el_type, ngroup_elements in el_type_hist.iteritems(): + if group_el_type.dimensions != mesh_bulk_dim: + continue + + nodes = np.empty((ambient_dim, ngroup_elements, el_type.node_count()), + np.float64) + el_vertex_count = group_el_type.vertex_count() + vertex_indices = np.empty( + (ngroup_elements, el_vertex_count), + np.float64) + i = 0 + + for el_vertices, el_nodes, el_type in zip( + self.element_vertices, self.element_nodes, self.element_types): + if el_type is not group_el_type: + continue + + nodes[:, i] = self.points[el_nodes].T + vertex_indices[i] = [vertex_index_gmsh_to_mine[v_nr] + for v_nr in el_vertices] + + i += 1 + + unit_nodes = (np.array(group_el_type.lexicographic_node_tuples(), + dtype=np.float64).T/group_el_type.order)*2 - 1 + + groups.append(MeshElementGroup( + group_el_type.order, + vertex_indices, + nodes, + unit_nodes=unit_nodes + )) + + return Mesh(vertices, groups) + +# }}} + + +def read_gmsh(filename, force_ambient_dimension=None): + """Read a gmsh mesh file from *filename* and return a + :class:`meshmode.mesh.Mesh`. + + :arg force_dimension: if not None, truncate point coordinates to + this many dimensions. + """ + from meshpy.gmsh_reader import read_gmsh + recv = GmshMeshReceiver() + read_gmsh(recv, filename, force_dimension=force_ambient_dimension) + return recv.get_mesh() + + +def generate_gmsh(source, dimensions, order=None, other_options=[], + extension="geo", gmsh_executable="gmsh", force_ambient_dimension=None): + """Run :cmd:`gmsh` on the input given by *source*, and return a + :class:`meshmode.mesh.Mesh` based on the result. + + :arg source: an instance of either :class:`FileSource` or + :class:`LiteralSource` + :arg force_ambient_dimension: if not None, truncate point coordinates to + this many dimensions. + """ + recv = GmshMeshReceiver() + + from meshpy.gmsh import GmshRunner + from meshpy.gmsh_reader import parse_gmsh + with GmshRunner(source, dimensions, order=order, + other_options=other_options, extension=extension, + gmsh_executable=gmsh_executable) as runner: + parse_gmsh(recv, runner.output_file, + force_dimension=force_ambient_dimension) + + return recv.get_mesh()