diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index 4477a9bfd0fdbf4e24659c059997e1e6156fb894..5ef91fe7e41dc464061e73438796afa174756560 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -51,6 +51,8 @@ def _map_unit_nodes_to_children(unit_nodes, tesselation): """ ref_vertices = np.array(tesselation.ref_vertices, dtype=np.float) + assert len(unit_nodes.shape) == 2 + for child_element in tesselation.children: center = np.vstack(ref_vertices[child_element[0]]) # Scale by 1/2 since sides in the tesselation have length 2. diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index a5319bc84f47ab6689f502629ff6009893033832..54aaabef97ac5f9bfc97b642ebe1a476608e0603 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -144,7 +144,13 @@ class QuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase): @property @memoize_method def unit_nodes(self): - return self._quadrature_rule().nodes + result = self._quadrature_rule().nodes + if len(result.shape) == 1: + result = np.array([result]) + + dim2, nunit_nodes = result.shape + assert dim2 == self.mesh_el_group.dim + return result @property @memoize_method diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index cb001e3266192a9280e313ae95ef6b93ff988d7d..16639442bdbaf10d6dd2def5c87bd8ff3254c25a 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -73,17 +73,19 @@ class Refiner(object): # {{{ constructor def __init__(self, mesh): - from meshmode.mesh.tesselate import tesselatetet, tesselatetri + from meshmode.mesh.tesselate import tesselateseg, tesselatetet, tesselatetri self.lazy = False self.seen_tuple = {} self.group_refinement_records = [] + seg_node_tuples, seg_result = tesselateseg() tri_node_tuples, tri_result = tesselatetri() tet_node_tuples, tet_result = tesselatetet() #quadrilateral_node_tuples = [ #print tri_result, tet_result - self.simplex_node_tuples = [None, None, tri_node_tuples, tet_node_tuples] + self.simplex_node_tuples = [ + None, seg_node_tuples, tri_node_tuples, tet_node_tuples] # Dimension-parameterized tesselations for refinement - self.simplex_result = [None, None, tri_result, tet_result] + self.simplex_result = [None, seg_result, tri_result, tet_result] #print tri_node_tuples, tri_result #self.simplex_node_tuples, self.simplex_result = tesselatetet() self.last_mesh = mesh @@ -594,8 +596,7 @@ class Refiner(object): midpoint_vertices = [] vertex_indices = grp.vertex_indices[iel_grp] #if simplex - if len(grp.vertex_indices[iel_grp]) == len(self.last_mesh.vertices)+1: - + if len(vertex_indices) == grp.dim + 1: # {{{ Get midpoints for all pairs of vertices for i in range(len(vertex_indices)): diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py index 0d0476f3454a7079ebd4df3b1ace1b604355edde..751142aa6d67ff7717b017ebe83f13a18e216ef3 100644 --- a/meshmode/mesh/refinement/utils.py +++ b/meshmode/mesh/refinement/utils.py @@ -88,10 +88,18 @@ def check_nodal_adj_against_geometry(mesh, tol=1e-12): nearby_vertex = mesh.vertices[:, nearby_vertex_index] transformation[:, inearby_vertex_index] = \ nearby_vertex - nearby_origin_vertex - bary_coord = np.linalg.solve(transformation, vertex_transformed) + bary_coord, residual = \ + np.linalg.lstsq(transformation, vertex_transformed)[0:2] + + is_in_element_span = ( + residual.size == 0 or + np.linalg.norm(vertex_transformed) == 0 or + (np.linalg.norm(residual) / + np.linalg.norm(vertex_transformed)) <= tol) is_connected = ( - np.sum(bary_coord) <= 1+tol + is_in_element_span + and np.sum(bary_coord) <= 1+tol and (bary_coord >= -tol).all()) el1 = group_and_iel_to_global_iel(nearby_igrp, nearby_iel) el2 = group_and_iel_to_global_iel(igrp, iel_grp) diff --git a/meshmode/mesh/tesselate.py b/meshmode/mesh/tesselate.py index 5d0e03aef4b9b2a8ac9ead1571dbdea1290ea2a7..c9527802526b4a5e6f34660a2d3200229b889ec5 100644 --- a/meshmode/mesh/tesselate.py +++ b/meshmode/mesh/tesselate.py @@ -6,6 +6,12 @@ def add_tuples(a, b): return tuple(ac+bc for ac, bc in zip(a, b)) +def tesselateseg(): + node_tuples = [(0,), (1,), (2,)] + result = [(0, 1), (1, 2)] + return [node_tuples, result] + + def tesselatetri(): result = [] diff --git a/test/test_refinement.py b/test/test_refinement.py index 9c36b9750ebe87a899a8dd5c5f17f4191ed1b0b4..bed61eb94c900db0f78b526cfe43ffababe4d825 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -31,7 +31,7 @@ from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) from meshmode.mesh.generation import ( # noqa - generate_icosahedron, generate_box_mesh) + generate_icosahedron, generate_box_mesh, make_curve_mesh, ellipse) from meshmode.mesh.refinement.utils import check_nodal_adj_against_geometry from meshmode.mesh.refinement import Refiner @@ -80,6 +80,15 @@ def uniform_refine_flags(mesh): # partial(random_refine_flags, 0.4), # 3), + ("3_to_1_ellipse_unif", + partial( + make_curve_mesh, + partial(ellipse, 3), + np.linspace(0, 1, 21), + order=1), + uniform_refine_flags, + 4), + ("rect2d_rand", partial(generate_box_mesh, ( np.linspace(0, 1, 3), @@ -138,6 +147,7 @@ def test_refinement(case_name, mesh_gen, flag_gen, num_generations): PolynomialEquidistantGroupFactory ]) @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ + ("circle", 1, [20, 30, 40]), ("blob", 2, [1e-1, 8e-2, 5e-2]), ("warp", 2, [4, 5, 6]), ("warp", 3, [4, 5, 6]), @@ -171,7 +181,12 @@ def test_refinement_connection( for mesh_par in mesh_pars: # {{{ get mesh - if mesh_name == "blob": + if mesh_name == "circle": + assert dim == 1 + h = 1 / mesh_par + mesh = make_curve_mesh( + partial(ellipse, 1), np.linspace(0, 1, mesh_par + 1), order=1) + elif mesh_name == "blob": assert dim == 2 h = mesh_par mesh = gen_blob_mesh(h)