From e284a28a2c7e953ffadab68cbda47e5de1a7102b Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 19 Sep 2016 23:00:29 -0500 Subject: [PATCH 1/4] Ensure unit_nodes has the right shape for the 1d case of QuadratureSimplexElementGroup --- meshmode/discretization/connection/refinement.py | 2 ++ meshmode/discretization/poly_element.py | 8 +++++++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index 4477a9bf..5ef91fe7 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 a5319bc8..54aaabef 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 -- GitLab From 55aa1ecc3825e072d89be6b379fe339639e73c68 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 19 Sep 2016 23:01:10 -0500 Subject: [PATCH 2/4] Refiner: Make 1d element refinement work --- meshmode/mesh/refinement/__init__.py | 11 ++++++----- meshmode/mesh/refinement/utils.py | 2 +- meshmode/mesh/tesselate.py | 6 ++++++ test/test_refinement.py | 19 +++++++++++++++++-- 4 files changed, 30 insertions(+), 8 deletions(-) diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index cb001e32..16639442 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 0d0476f3..d34fa14a 100644 --- a/meshmode/mesh/refinement/utils.py +++ b/meshmode/mesh/refinement/utils.py @@ -80,7 +80,7 @@ def check_nodal_adj_against_geometry(mesh, tol=1e-12): nearby_origin_vertex = mesh.vertices[ :, nearby_grp.vertex_indices[nearby_iel][0]] # noqa transformation = np.empty( - (len(mesh.vertices), nvertices_per_element-1)) + (mesh.ambient_dim, mesh.ambient_dim)) vertex_transformed = vertex - nearby_origin_vertex for inearby_vertex_index, nearby_vertex_index in enumerate( diff --git a/meshmode/mesh/tesselate.py b/meshmode/mesh/tesselate.py index 5d0e03ae..c9527802 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 9c36b975..fca3aa93 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, 20), + order=2), + 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) -- GitLab From cb97c26a29a6acb8559b6dbeb7316e8a2b145b2e Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Tue, 20 Sep 2016 14:32:11 -0500 Subject: [PATCH 3/4] check_nodal_adj_against_geometry(): Use least squares to find the barycentric coordinates. --- meshmode/mesh/refinement/utils.py | 13 ++++++++++--- test/test_refinement.py | 4 ++-- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py index d34fa14a..936f9107 100644 --- a/meshmode/mesh/refinement/utils.py +++ b/meshmode/mesh/refinement/utils.py @@ -80,7 +80,7 @@ def check_nodal_adj_against_geometry(mesh, tol=1e-12): nearby_origin_vertex = mesh.vertices[ :, nearby_grp.vertex_indices[nearby_iel][0]] # noqa transformation = np.empty( - (mesh.ambient_dim, mesh.ambient_dim)) + (len(mesh.vertices), nvertices_per_element-1)) vertex_transformed = vertex - nearby_origin_vertex for inearby_vertex_index, nearby_vertex_index in enumerate( @@ -88,10 +88,17 @@ 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 + 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/test/test_refinement.py b/test/test_refinement.py index fca3aa93..bed61eb9 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -84,8 +84,8 @@ def uniform_refine_flags(mesh): partial( make_curve_mesh, partial(ellipse, 3), - np.linspace(0, 1, 20), - order=2), + np.linspace(0, 1, 21), + order=1), uniform_refine_flags, 4), -- GitLab From 6cb47e9082491662ab33eda819162d7bf662666e Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Tue, 20 Sep 2016 14:36:13 -0500 Subject: [PATCH 4/4] residual => norm(residual) --- meshmode/mesh/refinement/utils.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py index 936f9107..751142aa 100644 --- a/meshmode/mesh/refinement/utils.py +++ b/meshmode/mesh/refinement/utils.py @@ -94,7 +94,8 @@ def check_nodal_adj_against_geometry(mesh, tol=1e-12): is_in_element_span = ( residual.size == 0 or np.linalg.norm(vertex_transformed) == 0 or - residual / np.linalg.norm(vertex_transformed) <= tol) + (np.linalg.norm(residual) / + np.linalg.norm(vertex_transformed)) <= tol) is_connected = ( is_in_element_span -- GitLab