diff --git a/mesh_sphere_packing/__init__.py b/mesh_sphere_packing/__init__.py index e69de29..73f738a 100644 --- a/mesh_sphere_packing/__init__.py +++ b/mesh_sphere_packing/__init__.py @@ -0,0 +1,6 @@ +from __future__ import print_function + +from meshpy import triangle, tet +import numpy as np +from numpy import random as npr +from numpy import linalg as npl diff --git a/mesh_sphere_packing/boundarypslg.py b/mesh_sphere_packing/boundarypslg.py index 8cce4f8..047871d 100755 --- a/mesh_sphere_packing/boundarypslg.py +++ b/mesh_sphere_packing/boundarypslg.py @@ -1,12 +1,4 @@ -from __future__ import print_function - -import sys -import os -from pprint import pprint - -from meshpy import triangle, tet -import numpy as np -from numpy import linalg as npl +from mesh_sphere_packing import * WITH_PBC = True @@ -15,59 +7,6 @@ # : of a sphere surface. This is confusing... -def output_mesh_poly(points, facets, markers, holes): - with open('mesh.poly', 'w') as f: - f.write('%d 3 0 1\n' % len(points)) - for i, p in enumerate(points): - f.write('%5d %+1.15e %+1.15e %+1.15e\n' % (i, p[0], p[1], p[2])) - f.write('%d 1\n' % len(facets)) - for i, (fac, m) in enumerate(zip(facets, markers)): - f.write('1 0 %d\n%d %d %d %d\n' % (m, 3, fac[0], fac[1], fac[2])) - if len(holes): - f.write('%d\n' % len(holes)) - for i, h in enumerate(holes): - f.write('%5d %+1.15e %+1.15e %+1.15e\n' % (i, h[0], h[1], h[2])) - else: - f.write('0\n') - - -def read_geomview(fname): - with open(fname, 'r') as f: - f.readline() - counts = [int(tok) for tok in f.readline().strip().split()] - points = np.array([ - [float(tok) for tok in f.readline().strip().split()] - for i in range(counts[0]) - ]) - tris = np.array([ - [int(tok) for tok in f.readline().strip().split()[1:]] - for i in range(counts[1]) - ]) - return points, tris - - -def load_segments_geomview(prefix): - geomfiles = [ - fname for fname in os.listdir(os.getcwd()) - if fname.startswith(prefix) and os.path.splitext(fname)[1] == '.off' - ] - geomfiles.sort(key=lambda x: int(os.path.splitext(x)[0].split('_')[1])) - return [read_geomview(fname) for fname in geomfiles] - - -def get_args(argv): - """Get command line arguments - :return: sphere center coordinates, x, y, z, sphere radius, r, - domain box side lengths, Lx, Ly, Lz. - """ - try: - return float(argv[1]), float(argv[2]), float(argv[3]), argv[4] - except IndexError: - raise UserWarning('Must specify Lx Ly Lz segment_file_prefix') - except ValueError: - raise UserWarning('Invalid arguments') - - def build_boundary_PSLGs(segments, Lx, Ly, Lz): # TODO : Break up this function a bit. @@ -277,90 +216,28 @@ def triangulate_PSLGs(pslgs): return triangulated_boundaries -def build_tet_mesh(segments, boundaries, Lx, Ly, Lz): - - def duplicate_lower_boundaries(lower_boundaries, Lx, Ly, Lz): - L = np.array([Lx, Ly, Lz]) - upper_boundaries = [] - for i, (points, tris, holes) in enumerate(lower_boundaries): - translate = np.array([[L[j] if j==i else 0. for j in range(3)]]) - points_upper = points.copy() + translate - tris_upper = tris.copy() - holes_upper = holes.copy() + translate - upper_boundaries.append((points_upper, tris_upper, holes_upper)) - return lower_boundaries + upper_boundaries - - def build_point_list(segments, boundaries): - vcount = 0 - all_points = [] - for points, tris in segments: - all_points.append(points) - tris += vcount - vcount += points.shape[0] - for points, tris, _ in boundaries: - all_points.append(points) - tris += vcount - vcount += points.shape[0] - return np.vstack(all_points) - - def build_facet_list(segments, boundaries): - all_facets = [tris for _, tris, _ in boundaries] - all_markers = [ - np.full(len(all_facets[0]), 1), np.full(len(all_facets[3]), 2), - np.full(len(all_facets[1]), 3), np.full(len(all_facets[4]), 4), - np.full(len(all_facets[2]), 5), np.full(len(all_facets[5]), 6), - ] - fcount = 7 - for _, tris in segments: - all_facets.append(tris) - all_markers.append(np.full(len(tris), fcount)) - fcount += 1 - return np.vstack(all_facets), np.hstack(all_markers) - - def build_hole_list(segments): - # TODO : Ultimately each sphere segment will contain hole data - all_holes = [] - for points, _ in segments: - all_holes.append(0.5 * (points.max(axis=0) + points.min(axis=0))) - return np.vstack(all_holes) - - boundaries = duplicate_lower_boundaries(boundaries, Lx, Ly, Lz) - - points = build_point_list(segments, boundaries) - # Fix boundary points to exactly zero - for i in range(3): - points[(np.isclose(points[:,i], 0.), i)] = 0. - - facets, markers = build_facet_list(segments, boundaries) - holes = build_hole_list(segments) - - output_mesh_poly(points, facets, markers, holes) +def boundarypslg(segments, Lx, Ly, Lz): + boundary_pslgs = build_boundary_PSLGs(segments, Lx, Ly, Lz) + return triangulate_PSLGs(boundary_pslgs) - rad_edge = 1.4 - min_angle = 18. - max_volume = None # 0.00001 - # TODO : Don't mix and match between setting options with argument string - # : and option class attributes. Pick one and be consistent. - options = tet.Options('pq{}/{}Y'.format(rad_edge, min_angle)) - options.docheck = 1 - options.verbose = 1 - mesh = tet.MeshInfo() - mesh.set_points(points) - mesh.set_facets(facets.tolist(), markers=markers.tolist()) - mesh.set_holes(holes) +def get_args(argv): + """Get command line arguments + :return: sphere center coordinates, x, y, z, sphere radius, r, + domain box side lengths, Lx, Ly, Lz. + """ + try: + return float(argv[1]), float(argv[2]), float(argv[3]), argv[4] + except IndexError: + raise UserWarning('Must specify Lx Ly Lz segment_file_prefix') + except ValueError: + raise UserWarning('Invalid arguments') - return tet.build(mesh, options=options, max_volume=max_volume) +if __name__ == '__main__': + import sys, devutils -def main(): Lx, Ly, Lz, prefix = get_args(sys.argv) - segments = load_segments_geomview(prefix) - boundary_pslgs = build_boundary_PSLGs(segments, Lx, Ly, Lz) - boundaries = triangulate_PSLGs(boundary_pslgs) - mesh = build_tet_mesh(segments, boundaries, Lx, Ly, Lz) - mesh.write_vtk('mesh') - - -if __name__ == '__main__': - main() + segments = load_geomview_files(prefix) + boundaries = boundarypslg(segments, Lx, Ly, Lz) + devutils.write_boundaries_geomview(boundaries) diff --git a/mesh_sphere_packing/build.py b/mesh_sphere_packing/build.py new file mode 100644 index 0000000..e13572f --- /dev/null +++ b/mesh_sphere_packing/build.py @@ -0,0 +1,38 @@ + +import faulthandler +faulthandler.enable() + +import sys + +from mesh_sphere_packing import * +from mesh_sphere_packing.splitsphere import splitsphere +from mesh_sphere_packing.boundarypslg import boundarypslg +from mesh_sphere_packing.tetmesh import build_tetmesh + + +def build(args): + x, y, z, r, Lx, Ly, Lz = args + segments = splitsphere(args) + boundaries = boundarypslg(segments, Lx, Ly, Lz) + return build_tetmesh(segments, boundaries, Lx, Ly, Lz) + + +def get_args(argv): + # TODO : Replace with argparse. + """Get command line arguments + :return: sphere center coordinates, x, y, z, sphere radius, r, + domain box side lengths, Lx, Ly, Lz. + """ + try: + return float(argv[1]), float(argv[2]), float(argv[3]), float(argv[4]),\ + float(argv[5]), float(argv[6]), float(argv[7]) + except IndexError: + raise UserWarning('Must specify x0 y0 z0 r Lx Ly Lz') + except ValueError: + raise UserWarning('Invalid arguments') + + +if __name__ == '__main__': + args = get_args(sys.argv) + mesh = build(args) + mesh.write_vtk('mesh.vtk') diff --git a/mesh_sphere_packing/devutils.py b/mesh_sphere_packing/devutils.py index 9c17b28..a49d2ec 100644 --- a/mesh_sphere_packing/devutils.py +++ b/mesh_sphere_packing/devutils.py @@ -16,6 +16,11 @@ def write_geomview(points, tris, fname): f.write('3 %d %d %d\n' % (t[0], t[1], t[2])) +def output_segments_geomview(segments): + for i, (points, tris) in enumerate(segments): + write_geomview(points, tris, './segments_%d.off' % i) + + def output_boundaries_geomview(boundaries): for i, (points, tris, _) in enumerate(boundaries): write_geomview(points, tris, './boundary_%d.off' % i) diff --git a/mesh_sphere_packing/splitsphere.py b/mesh_sphere_packing/splitsphere.py index 42aec7a..bbbb107 100755 --- a/mesh_sphere_packing/splitsphere.py +++ b/mesh_sphere_packing/splitsphere.py @@ -1,68 +1,8 @@ -#!/usr/bin/env python -from __future__ import print_function +from mesh_sphere_packing import * -import sys - -import numpy as np -from numpy import random as npr -from numpy import linalg as npl from scipy.spatial.qhull import ConvexHull -def plot_segment_points_3d(data): - """This function is for visualising sphere points to aid debugging/development - """ - # TODO : remove this after development complete - import matplotlib.pyplot as plt - from mpl_toolkits.mplot3d import Axes3D - fig = plt.figure() - ax = fig.add_subplot(111, projection='3d') - ax.scatter(data[:,0], data[:,1], zs=data[:,2]) - plt.show() - - -def plot_polys(polys): - """This function is for visualising sphere segments to aid debugging/development - """ - # TODO : remove this after development complete - from mpl_toolkits.mplot3d.art3d import Poly3DCollection - import matplotlib.pyplot as plt - fig = plt.figure() - ax = fig.add_subplot(111, projection='3d') - fc = ["crimson" if i%2 else "gold" for i in range(len(polys))] - ax.add_collection3d(Poly3DCollection(poly3d, facecolors=fc, linewidths=1)) - plt.show() - - -def write_geomview(points, tris, fname): - with open(fname, 'w') as f: - f.write('OFF\n') - f.write('%d %d %d\n' % (len(points), len(tris), 0)) - for p in points: - f.write('%+1.15e %+1.15e %+1.15e\n' % (p[0], p[1], p[2])) - for t in tris: - f.write('3 %d %d %d\n' % (t[0], t[1], t[2])) - - -def output_segments_geomview(segments): - for i, (tris, points) in enumerate(segments): - write_geomview(points, tris, './segment_%d.off' % i) - - -def get_args(argv): - """Get command line arguments - :return: sphere center coordinates, x, y, z, sphere radius, r, - domain box side lengths, Lx, Ly, Lz. - """ - try: - return float(argv[1]), float(argv[2]), float(argv[3]), float(argv[4]),\ - float(argv[5]), float(argv[6]), float(argv[7]) - except IndexError: - raise UserWarning('Must specify x0 y0 z0 r Lx Ly Lz') - except ValueError: - raise UserWarning('Invalid arguments') - - def gen_sphere_spiral_points(x, y, z, r, num_points=200): indices = np.arange(0, num_points, dtype=float) + 0.5 @@ -222,6 +162,7 @@ def get_add_phi_points(phi_1, phi_2, r, ds, fn_x, fn_y, fn_z): return segments + def translate_segments(segments, Lx, Ly, Lz): translations = [ np.array([0., 0., 0.]), @@ -233,10 +174,11 @@ def translate_segments(segments, Lx, Ly, Lz): np.array([0., Ly, Lz]), np.array([Lx, Ly, Lz]), ] - for (tris, points), t in zip(segments, translations): + for (points, tris), t in zip(segments, translations): points += t return segments + def triangulate_segment_points(point_sets, x, y, z): """For now this function returns a set of convex hulls :return: list of scipy.spatial.qhull.ConvexHull instances @@ -275,8 +217,8 @@ def reindex_tris(tris, points): points_reindexed[count] = points[idx] count += 1 return ( + points_reindexed[:count], np.array([[reindex[idx] for idx in t] for t in tris]), - points_reindexed[:count] ) def extract_surf_from_chull(chull, x, y, z): @@ -292,34 +234,13 @@ def extract_surf_from_chull(chull, x, y, z): return [extract_surf_from_chull(h, x, y, z) for h in chulls] -def triangulate_segment_points_cgal_alt(point_sets): - from CGAL.CGAL_Kernel import Point_3 - from CGAL import CGAL_Convex_hull_3 - from CGAL.CGAL_Polyhedron_3 import Polyhedron_3 - """For now this function returns a set of convex hulls - :return: list of scipy.spatial.qhull.ConvexHull instances - """ - # TODO : this function needs to strip out only the triangles which - # : correspond to the surface of the sphere from the triangles - # : in the convex hull generated by qhull. - res = [] - for points in point_sets: - points_cgal = [Point_3(*x) for x in points] - poly = Polyhedron_3() - CGAL_Convex_hull_3.convex_hull_3(points_cgal, poly) - res.append(poly) - return res - - def smooth_segments(segments): # TODO : implement Laplacian smoothing of inner segment vertices return segments -if __name__ == '__main__': - # x, y, z, r, Lx, Ly, Lz = get_args(sys.argv) - x, y, z, r, Lx, Ly, Lz = 0., 0., 0., 0.5, 2., 2., 2. - #points = gen_sphere_hypercube_points(x, y, z, r) +def splitsphere(args): + x, y, z, r, Lx, Ly, Lz = args points = gen_sphere_spiral_points(x, y, z, r) points = filter_points(points, r) segment_point_sets = split_sphere_points_segments( @@ -327,5 +248,25 @@ def smooth_segments(segments): ) segments = triangulate_segment_points(segment_point_sets, x, y, z) segments = smooth_segments(segments) - segments = translate_segments(segments, Lx, Ly, Lz) + return translate_segments(segments, Lx, Ly, Lz) + + +def get_args(argv): + """Get command line arguments + :return: sphere center coordinates, x, y, z, sphere radius, r, + domain box side lengths, Lx, Ly, Lz. + """ + try: + return float(argv[1]), float(argv[2]), float(argv[3]), float(argv[4]),\ + float(argv[5]), float(argv[6]), float(argv[7]) + except IndexError: + raise UserWarning('Must specify x0 y0 z0 r Lx Ly Lz') + except ValueError: + raise UserWarning('Invalid arguments') + + +if __name__ == '__main__': + import sys + + segments = splitsphere(get_args(sys.argv)) output_segments_geomview(segments) diff --git a/mesh_sphere_packing/tetmesh.py b/mesh_sphere_packing/tetmesh.py new file mode 100644 index 0000000..f58e035 --- /dev/null +++ b/mesh_sphere_packing/tetmesh.py @@ -0,0 +1,98 @@ +from mesh_sphere_packing import * + + +def build_tetmesh(segments, boundaries, Lx, Ly, Lz): + + def duplicate_lower_boundaries(lower_boundaries, Lx, Ly, Lz): + L = np.array([Lx, Ly, Lz]) + upper_boundaries = [] + for i, (points, tris, holes) in enumerate(lower_boundaries): + translate = np.array([[L[j] if j==i else 0. for j in range(3)]]) + points_upper = points.copy() + translate + tris_upper = tris.copy() + holes_upper = holes.copy() + translate + upper_boundaries.append((points_upper, tris_upper, holes_upper)) + return lower_boundaries + upper_boundaries + + def build_point_list(segments, boundaries): + vcount = 0 + all_points = [] + for points, tris in segments: + all_points.append(points) + tris += vcount + vcount += points.shape[0] + for points, tris, _ in boundaries: + all_points.append(points) + tris += vcount + vcount += points.shape[0] + return np.vstack(all_points) + + def build_facet_list(segments, boundaries): + all_facets = [tris for _, tris, _ in boundaries] + all_markers = [ + np.full(len(all_facets[0]), 1), np.full(len(all_facets[3]), 2), + np.full(len(all_facets[1]), 3), np.full(len(all_facets[4]), 4), + np.full(len(all_facets[2]), 5), np.full(len(all_facets[5]), 6), + ] + fcount = 7 + for _, tris in segments: + all_facets.append(tris) + all_markers.append(np.full(len(tris), fcount)) + fcount += 1 + return np.vstack(all_facets), np.hstack(all_markers) + + def build_hole_list(segments): + # TODO : Ultimately each sphere segment will contain hole data + all_holes = [] + for points, _ in segments: + all_holes.append(0.5 * (points.max(axis=0) + points.min(axis=0))) + return np.vstack(all_holes) + + boundaries = duplicate_lower_boundaries(boundaries, Lx, Ly, Lz) + + points = build_point_list(segments, boundaries) + # Fix boundary points to exactly zero + for i in range(3): + points[(np.isclose(points[:,i], 0.), i)] = 0. + + facets, markers = build_facet_list(segments, boundaries) + holes = build_hole_list(segments) + + rad_edge = 1.4 + min_angle = 18. + max_volume = None # 0.00001 + # TODO : Don't mix and match between setting options with argument string + # : and option class attributes. Pick one and be consistent. + options = tet.Options('pq{}/{}Y'.format(rad_edge, min_angle)) + options.docheck = 1 + options.verbose = 1 + + mesh = tet.MeshInfo() + mesh.set_points(points) + mesh.set_facets(facets.tolist(), markers=markers.tolist()) + mesh.set_holes(holes) + + return tet.build(mesh, options=options, max_volume=max_volume) + + +def get_args(argv): + """Get command line arguments + :return: sphere center coordinates, x, y, z, sphere radius, r, + domain box side lengths, Lx, Ly, Lz. + """ + try: + return float(argv[1]), float(argv[2]), float(argv[3]), argv[4], argv[5] + except IndexError: + raise UserWarning('Must specify Lx Ly Lz segment_file_prefix') + except ValueError: + raise UserWarning('Invalid arguments') + + +if __name__ == '__main__': + import sys, devutils + + Lx, Ly, Lz, segment_prefix, boundary_prefix = get_args(sys.argv) + segments = devutils.load_geomview_files(segment_prefix) + boundaries = devutils.load_geomview_files(boundaries_prefix) + mesh = build_tetmesh(segments, boundaries, Lx, Ly, Lz) + mesh.write_vtk('mesh')