From 07f8e10d5f2005493863dd26b9b3eef37a84d90a Mon Sep 17 00:00:00 2001
From: Shivam Gupta <sgupta72@illinois.edu>
Date: Sat, 23 Apr 2016 23:23:57 -0500
Subject: [PATCH] Done
---
examples/plot-connectivity.py | 311 +++++++++++++++++++++++++++++++---
1 file changed, 291 insertions(+), 20 deletions(-)
diff --git a/examples/plot-connectivity.py b/examples/plot-connectivity.py
index c5a6439b..1de09830 100644
--- a/examples/plot-connectivity.py
+++ b/examples/plot-connectivity.py
@@ -4,7 +4,7 @@ import numpy as np # noqa
import pyopencl as cl
import random
import os
-order = 4
+order = 1
def main():
@@ -13,22 +13,20 @@ def main():
from meshmode.mesh.generation import ( # noqa
generate_icosphere, generate_icosahedron,
- generate_torus)
+ generate_torus, generate_box_mesh)
#mesh = generate_icosphere(1, order=order)
#mesh = generate_icosahedron(1, order=order)
- mesh = generate_torus(3, 1, order=order)
+ #mesh = generate_torus(3, 1, order=order)
+ mesh = generate_box_mesh(3*(np.linspace(0, 1, 5),))
from meshmode.mesh.refinement import Refiner
r = Refiner(mesh)
#mesh = r.refine(0)
- from meshmode.discretization import Discretization
from meshmode.discretization.poly_element import \
PolynomialWarpAndBlendGroupFactory
- discr = Discretization(
- cl_ctx, mesh, PolynomialWarpAndBlendGroupFactory(order))
-
- from meshmode.discretization.visualization import make_visualizer
- vis = make_visualizer(queue, discr, order=1)
+ discr = Discretization(cl_ctx, mesh, PolynomialWarpAndBlendGroupFactory(order))
+ from meshmode.discretization.visualization import make_visualizer
+ vis = make_visualizer(queue, discr, order)
os.remove("geometry.vtu")
os.remove("connectivity.vtu")
vis.write_vtk_file("geometry.vtu", [
@@ -41,28 +39,223 @@ def main():
write_mesh_connectivity_vtk_file("connectivity.vtu",
mesh)
+#construct vertex vertex_index
+def get_vertex(mesh, vertex_index):
+ vertex = np.empty([len(mesh.vertices)])
+ for i_cur_dim, cur_dim_coords in enumerate(mesh.vertices):
+ vertex[i_cur_dim] = cur_dim_coords[vertex_index]
+ return vertex
+
+def test_refiner_connectivity(mesh):
+ def group_and_iel_to_global_iel(igrp, iel):
+ return mesh.groups[igrp].element_nr_base + iel
+
+ from meshmode.mesh.tools import make_element_lookup_tree
+ tree = make_element_lookup_tree(mesh)
+
+ from meshmode.mesh.processing import find_bounding_box
+ bbox_min, bbox_max = find_bounding_box(mesh)
+
+ extent = bbox_max-bbox_min
+ cnx = mesh.element_connectivity
+ nvertices_per_element = len(mesh.groups[0].vertex_indices[0])
+ #stores connectivity obtained from geometry using barycentric coordinates
+ connected_to_element_geometry = np.zeros([mesh.nelements, mesh.nelements], dtype=bool)
+ #store connectivity obtained from mesh's connectivity
+ connected_to_element_connectivity = np.zeros([mesh.nelements, mesh.nelements], dtype=bool)
+
+ import six
+ for igrp, grp in enumerate(mesh.groups):
+ iel_base = grp.element_nr_base
+ for iel_grp in six.moves.range(grp.nelements):
+
+ iel_g = group_and_iel_to_global_iel(igrp, iel_grp)
+ nb_starts = cnx.neighbors_starts
+ for nb_iel_g in cnx.neighbors[nb_starts[iel_g]:nb_starts[iel_g+1]]:
+ connected_to_element_connectivity[nb_iel_g, iel_g] = True
+ connected_to_element_connectivity[iel_g, nb_iel_g] = True
+
+ for vertex_index in grp.vertex_indices[iel_grp]:
+ vertex = get_vertex(mesh, vertex_index)
+ #check which elements touch this vertex
+ for bounding_igrp, bounding_iel in tree.generate_matches(vertex):
+ if bounding_igrp == igrp and bounding_iel == iel_grp:
+ continue
+ bounding_grp = mesh.groups[bounding_igrp]
+
+ last_bounding_vertex = get_vertex(mesh, \
+ bounding_grp.vertex_indices[bounding_iel][nvertices_per_element-1])
+ transformation = np.empty([len(mesh.vertices), nvertices_per_element-1])
+ vertex_transformed = vertex - last_bounding_vertex
+ for ibounding_vertex_index, bounding_vertex_index in \
+ enumerate(bounding_grp.vertex_indices[bounding_iel][:nvertices_per_element-1]):
+ bounding_vertex = get_vertex(mesh, bounding_vertex_index)
+ transformation[:,ibounding_vertex_index] = bounding_vertex - last_bounding_vertex
+ barycentric_coordinates = np.linalg.solve(transformation, vertex_transformed)
+ is_connected = True
+ sum_of_coords = 0.0
+ for coord in barycentric_coordinates:
+ if coord < 0:
+ is_connected = False
+ sum_of_coords += coord
+ if sum_of_coords > 1:
+ is_connected = False
+ if is_connected:
+ connected_to_element_geometry[group_and_iel_to_global_iel(bounding_igrp, bounding_iel),\
+ group_and_iel_to_global_iel(igrp, iel_grp)] = True
+ connected_to_element_geometry[group_and_iel_to_global_iel(igrp, iel_grp),\
+ group_and_iel_to_global_iel(bounding_igrp, bounding_iel)] = True
+ print "GEOMETRY: "
+ print connected_to_element_geometry
+ print "CONNECTIVITY: "
+ print connected_to_element_connectivity
+ cmpmatrix = (connected_to_element_geometry == connected_to_element_connectivity)
+ print cmpmatrix
+ print np.where(cmpmatrix == False)
+ if not (connected_to_element_geometry == connected_to_element_connectivity).all():
+ cmpmatrix = connected_to_element_connectivity == connected_to_element_geometry
+ for ii, i in enumerate(cmpmatrix):
+ for ij, j in enumerate(i):
+ if j == False:
+ print ii, ij, "GEOMETRY: ", connected_to_element_geometry[ii][ij], "CONNECTIVITY: ", connected_to_element_connectivity[ii][ij]
+ #assert (connected_to_element_geometry == connected_to_element_connectivity).all()
+
+def test_refiner_connectivity_efficient(mesh):
+ def group_and_iel_to_global_iel(igrp, iel):
+ return mesh.groups[igrp].element_nr_base + iel
+
+ from meshmode.mesh.tools import make_element_lookup_tree
+ tree = make_element_lookup_tree(mesh)
+
+ from meshmode.mesh.processing import find_bounding_box
+ bbox_min, bbox_max = find_bounding_box(mesh)
+
+ extent = bbox_max-bbox_min
+ cnx = mesh.element_connectivity
+ nvertices_per_element = len(mesh.groups[0].vertex_indices[0])
+ #stores connectivity obtained from geometry using barycentric coordinates
+ #connected_to_element_geometry = np.zeros([mesh.nelements, mesh.nelements], dtype=bool)
+ #store connectivity obtained from mesh's connectivity
+ #connected_to_element_connectivity = np.zeros([mesh.nelements, mesh.nelements], dtype=bool)
+ connected_to_element_geometry = []
+ connected_to_element_connectivity = []
+ import six
+ for i in six.moves.range(mesh.nelements):
+ connected_to_element_geometry.append([])
+ connected_to_element_connectivity.append([])
+
+ for igrp, grp in enumerate(mesh.groups):
+ iel_base = grp.element_nr_base
+ for iel_grp in six.moves.range(grp.nelements):
+
+ iel_g = group_and_iel_to_global_iel(igrp, iel_grp)
+ nb_starts = cnx.neighbors_starts
+ for nb_iel_g in cnx.neighbors[nb_starts[iel_g]:nb_starts[iel_g+1]]:
+ if iel_g not in connected_to_element_connectivity[nb_iel_g]:
+ connected_to_element_connectivity[nb_iel_g].append(iel_g)
+ if nb_iel_g not in connected_to_element_connectivity[iel_g]:
+ connected_to_element_connectivity[iel_g].append(nb_iel_g)
+ #connected_to_element_connectivity[nb_iel_g, iel_g] = True
+ #connected_to_element_connectivity[iel_g, nb_iel_g] = True
+
+ for vertex_index in grp.vertex_indices[iel_grp]:
+ vertex = get_vertex(mesh, vertex_index)
+ #check which elements touch this vertex
+ for bounding_igrp, bounding_iel in tree.generate_matches(vertex):
+ if bounding_igrp == igrp and bounding_iel == iel_grp:
+ continue
+ bounding_grp = mesh.groups[bounding_igrp]
+
+ last_bounding_vertex = get_vertex(mesh, \
+ bounding_grp.vertex_indices[bounding_iel][nvertices_per_element-1])
+ transformation = np.empty([len(mesh.vertices), nvertices_per_element-1])
+ vertex_transformed = vertex - last_bounding_vertex
+ for ibounding_vertex_index, bounding_vertex_index in \
+ enumerate(bounding_grp.vertex_indices[bounding_iel][:nvertices_per_element-1]):
+ bounding_vertex = get_vertex(mesh, bounding_vertex_index)
+ transformation[:,ibounding_vertex_index] = bounding_vertex - last_bounding_vertex
+ barycentric_coordinates = np.linalg.solve(transformation, vertex_transformed)
+ is_connected = True
+ sum_of_coords = 0.0
+ for coord in barycentric_coordinates:
+ if coord < 0:
+ is_connected = False
+ sum_of_coords += coord
+ if sum_of_coords > 1:
+ is_connected = False
+ if is_connected:
+ el1 = group_and_iel_to_global_iel(bounding_igrp, bounding_iel)
+ el2 = group_and_iel_to_global_iel(igrp, iel_grp)
+ if el2 not in connected_to_element_geometry[el1]:
+ connected_to_element_geometry[el1].append(el2)
+ if el1 not in connected_to_element_geometry[el2]:
+ connected_to_element_geometry[el2].append(el1)
+ from collections import Counter
+ for i in six.moves.range(mesh.nelements):
+ assert(Counter(connected_to_element_connectivity[i]) == Counter(connected_to_element_geometry[i]))
+
def main2():
cl_ctx = cl.create_some_context()
queue = cl.CommandQueue(cl_ctx)
-
from meshmode.mesh.generation import ( # noqa
generate_icosphere, generate_icosahedron,
- generate_torus, generate_regular_rect_mesh)
- #mesh = generate_icosphere(1, order=order)
+ generate_torus, generate_regular_rect_mesh,
+ generate_box_mesh)
+# mesh = generate_icosphere(1, order=order)
#mesh = generate_icosahedron(1, order=order)
-# mesh = generate_torus(3, 1, order=order)
- mesh = generate_regular_rect_mesh()
+ #mesh = generate_torus(3, 1, order=order)
+ #mesh = generate_regular_rect_mesh()
+ mesh = generate_box_mesh(3*(np.linspace(0, 1, 3),))
+ print "NELEMENTS: ", mesh.nelements
+ #print mesh
+ test_refiner_connectivity(mesh);
from meshmode.mesh.refinement import Refiner
r = Refiner(mesh)
- times = random.randint(1, 1)
+ #random.seed(0)
+ times = random.randint(4, 4)
+ count = 0
+ poss_flags = np.ones(len(mesh.groups[0].vertex_indices))
for time in xrange(times):
+ print "NELS:", mesh.nelements
+ count = 0
flags = np.zeros(len(mesh.groups[0].vertex_indices))
-# for i in xrange(0, len(flags)):
-# flags[i] = random.randint(0, 1)
+ for i in xrange(0, len(flags)):
+ flags[i] = random.randint(0, 1)
+ if flags[i] == 1:
+ #print "REFINING:", i
+ count += 1
mesh = r.refine(flags)
+ #from meshmode.mesh.visualization import draw_2d_mesh
+ #draw_2d_mesh(mesh, True, True, True, fill=None)
+ import matplotlib.pyplot as pt
+ pt.show()
+
+ poss_flags = np.zeros(len(mesh.groups[0].vertex_indices))
+ for i in xrange(0, len(flags)):
+ poss_flags[i] = flags[i]
+ for i in xrange(len(flags), len(poss_flags)):
+ poss_flags[i] = 1
+ '''
+ flags = np.zeros(len(mesh.groups[0].vertex_indices))
+ flags[0] = 1
+ flags[1] = 1
+ mesh = r.refine(flags)
+ flags = np.zeros(len(mesh.groups[0].vertex_indices))
+ flags[0] = 1
+ flags[1] = 1
+ flags[2] = 1
+ mesh = r.refine(flags)
+ '''
+ test_refiner_connectivity(mesh)
+ #r.print_rays(70)
+ #r.print_rays(117)
+ #r.print_hanging_elements(10)
+ #r.print_hanging_elements(117)
+ #r.print_hanging_elements(757)
+ '''
from meshmode.mesh.visualization import draw_2d_mesh
- draw_2d_mesh(mesh, True, True, True);
+ draw_2d_mesh(mesh, True, True, True, fill=None)
import matplotlib.pyplot as pt
pt.show()
'''
@@ -76,6 +269,7 @@ def main2():
from meshmode.discretization.visualization import make_visualizer
vis = make_visualizer(queue, discr, order)
os.remove("connectivity2.vtu")
+ os.remove("geometry2.vtu")
vis.write_vtk_file("geometry2.vtu", [
("f", discr.nodes()[0]),
])
@@ -85,7 +279,84 @@ def main2():
write_mesh_connectivity_vtk_file("connectivity2.vtu",
mesh)
- '''
+
+def all_refine(num_mesh, depth, fname):
+ import six
+ from meshmode.mesh.generation import ( # noqa
+ generate_icosphere, generate_icosahedron,
+ generate_torus, generate_regular_rect_mesh,
+ generate_box_mesh)
+ import timeit
+ nelements = []
+ runtimes = []
+ for el_fact in six.moves.range(2, num_mesh+2):
+ mesh = generate_box_mesh(3*(np.linspace(0, 1, el_fact),))
+ from meshmode.mesh.refinement import Refiner
+ r = Refiner(mesh)
+ for time in xrange(depth):
+ flags = np.ones(len(mesh.groups[0].vertex_indices))
+ if time < depth-1:
+ mesh = r.refine(flags)
+ else:
+ start = timeit.default_timer()
+ mesh = r.refine(flags)
+ stop = timeit.default_timer()
+ nelements.append(mesh.nelements)
+ runtimes.append(stop-start)
+ #test_refiner_connectivity_efficient(mesh)
+ import matplotlib.pyplot as pt
+ pt.plot(nelements, runtimes, "o")
+ pt.savefig(fname)
+ pt.clf()
+ #pt.show()
+
+def uniform_refine(num_mesh, fract, depth, fname):
+ import six
+ from meshmode.mesh.generation import ( # noqa
+ generate_icosphere, generate_icosahedron,
+ generate_torus, generate_regular_rect_mesh,
+ generate_box_mesh)
+ import timeit
+ nelements = []
+ runtimes = []
+ for el_fact in six.moves.range(2, num_mesh+2):
+ mesh = generate_box_mesh(3*(np.linspace(0, 1, el_fact),))
+ from meshmode.mesh.refinement import Refiner
+ r = Refiner(mesh)
+ all_els = range(mesh.nelements)
+ for time in xrange(depth):
+ flags = np.zeros(mesh.nelements)
+ from random import shuffle
+ shuffle(all_els)
+ nels_this_round = 0
+ for i in six.moves.range(len(all_els)):
+ if i / len(flags) > fract:
+ break
+ flags[all_els[i]] = 1
+ nels_this_round += 1
+
+ if time < depth-1:
+ mesh = r.refine(flags)
+ else:
+ start = timeit.default_timer()
+ mesh = r.refine(flags)
+ stop = timeit.default_timer()
+ nelements.append(mesh.nelements)
+ runtimes.append(stop-start)
+ all_els = []
+ for i in six.moves.range(len(flags)):
+ if flags[i]:
+ all_els.append(i)
+ for i in six.moves.range(len(flags), mesh.nelements):
+ all_els.append(i)
+ test_refiner_connectivity(mesh)
+ import matplotlib.pyplot as pt
+ pt.plot(nelements, runtimes, "o")
+ pt.savefig(fname)
+ pt.clf()
if __name__ == "__main__":
-# main()
+ print "HEREERERE"
+ #all_refine(3, 2, 'all_a.pdf')
+ #all_refine(3, 3, 'all_b.pdf')
+ #uniform_refine(3, 0.2, 3, 'uniform_a.pdf')
main2()
--
GitLab