diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 3f52ec51..2e650e68 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -933,6 +933,7 @@ def make_mesh( node_vertex_consistency_tolerance: Optional[float] = None, skip_element_orientation_test: bool = False, force_positive_orientation: bool = False, + face_vertex_indices_to_tags=None, ) -> "Mesh": """Construct a new mesh from a given list of *groups*. @@ -1021,6 +1022,15 @@ def make_mesh( nodal_adjacency = ( NodalAdjacency(neighbors_starts=nb_starts, neighbors=nbs)) + face_vert_ind_to_tags_local = None + if face_vertex_indices_to_tags is not None: + face_vert_ind_to_tags_local = face_vertex_indices_to_tags.copy() + + if (facial_adjacency_groups is False or facial_adjacency_groups is None): + if face_vertex_indices_to_tags is not None: + facial_adjacency_groups = _compute_facial_adjacency_from_vertices( + groups, np.int32, np.int8, face_vertex_indices_to_tags) + if ( facial_adjacency_groups is not False and facial_adjacency_groups is not None): @@ -1047,8 +1057,13 @@ def make_mesh( if force_positive_orientation: if mesh.dim == mesh.ambient_dim: import meshmode.mesh.processing as mproc + mesh_making_kwargs = { + "face_vertex_indices_to_tags": face_vert_ind_to_tags_local + } mesh = mproc.perform_flips( - mesh, mproc.find_volume_mesh_element_orientations(mesh) < 0) + mesh=mesh, + flip_flags=mproc.find_volume_mesh_element_orientations(mesh) < 0, + skip_tests=False, mesh_making_kwargs=mesh_making_kwargs) else: raise ValueError("cannot enforce positive element orientation " "on non-volume meshes") diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 84c0c863..22148f9a 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -248,6 +248,7 @@ def get_mesh(self, return_tag_to_elements_map=False): # compute facial adjacency for Mesh if there is tag information facial_adjacency_groups = None + face_vert_ind_to_tags_local = face_vertex_indices_to_tags.copy() if is_conforming and self.tags: from meshmode.mesh import _compute_facial_adjacency_from_vertices facial_adjacency_groups = _compute_facial_adjacency_from_vertices( @@ -257,6 +258,7 @@ def get_mesh(self, return_tag_to_elements_map=False): vertices, groups, is_conforming=is_conforming, facial_adjacency_groups=facial_adjacency_groups, + face_vertex_indices_to_tags=face_vert_ind_to_tags_local, **self.mesh_construction_kwargs) return (mesh, tag_to_elements) if return_tag_to_elements_map else mesh diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index a61349da..74a17726 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -799,7 +799,8 @@ def flip_element_group( def perform_flips( mesh: Mesh, flip_flags: np.ndarray, - skip_tests: bool = False) -> Mesh: + skip_tests: bool = False, + mesh_making_kwargs = None) -> Mesh: """ :arg flip_flags: A :class:`numpy.ndarray` with :attr:`meshmode.mesh.Mesh.nelements` entries @@ -808,6 +809,8 @@ def perform_flips( """ if mesh.vertices is None: raise ValueError("Mesh must have vertices to perform flips") + if mesh_making_kwargs is None: + mesh_making_kwargs = {} flip_flags = flip_flags.astype(bool) @@ -823,9 +826,8 @@ def perform_flips( new_groups.append(new_grp) return make_mesh( - mesh.vertices, new_groups, skip_tests=skip_tests, - is_conforming=mesh.is_conforming, - ) + mesh.vertices, groups=new_groups, skip_tests=skip_tests, + is_conforming=mesh.is_conforming, **mesh_making_kwargs) # }}} diff --git a/test/3x3.msh b/test/3x3.msh new file mode 100644 index 00000000..faa3b135 --- /dev/null +++ b/test/3x3.msh @@ -0,0 +1,38 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +1 +2 10 "fluid" +$EndPhysicalNames +$Nodes +16 +1 0 0 0 +2 1 0 0 +3 1 0.5 0 +4 0 0.5 0 +5 0.3333333333333333 0 0 +6 0.6666666666666666 0 0 +7 1 0.1666666666666667 0 +8 1 0.3333333333333333 0 +9 0.6666666666666667 0.5 0 +10 0.3333333333333334 0.5 0 +11 0 0.3333333333333334 0 +12 0 0.1666666666666667 0 +13 0.3333333333333334 0.1666666666666667 0 +14 0.3333333333333334 0.3333333333333333 0 +15 0.6666666666666667 0.1666666666666667 0 +16 0.6666666666666667 0.3333333333333333 0 +$EndNodes +$Elements +9 +1 3 2 10 1 1 5 13 12 +2 3 2 10 1 12 13 14 11 +3 3 2 10 1 11 14 10 4 +4 3 2 10 1 5 6 15 13 +5 3 2 10 1 13 15 16 14 +6 3 2 10 1 14 16 9 10 +7 3 2 10 1 6 2 7 15 +8 3 2 10 1 15 7 8 16 +9 3 2 10 1 16 8 3 9 +$EndElements diff --git a/test/3x3_bound.msh b/test/3x3_bound.msh new file mode 100644 index 00000000..de74e581 --- /dev/null +++ b/test/3x3_bound.msh @@ -0,0 +1,54 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +5 +1 2 "Left" +1 3 "Right" +1 4 "Bottom" +1 5 "Top" +2 1 "fluid" +$EndPhysicalNames +$Nodes +16 +1 0 0 0 +2 1 0 0 +3 1 0.5 0 +4 0 0.5 0 +5 0.3333333333333333 0 0 +6 0.6666666666666666 0 0 +7 1 0.1666666666666667 0 +8 1 0.3333333333333333 0 +9 0.6666666666666667 0.5 0 +10 0.3333333333333334 0.5 0 +11 0 0.3333333333333334 0 +12 0 0.1666666666666667 0 +13 0.3333333333333334 0.1666666666666667 0 +14 0.3333333333333334 0.3333333333333333 0 +15 0.6666666666666667 0.1666666666666667 0 +16 0.6666666666666667 0.3333333333333333 0 +$EndNodes +$Elements +21 +1 1 2 2 1 1 12 +2 1 2 2 1 12 11 +3 1 2 2 1 11 4 +4 1 2 3 1 7 2 +5 1 2 3 1 8 7 +6 1 2 3 1 3 8 +7 1 2 4 1 5 1 +8 1 2 4 1 6 5 +9 1 2 4 1 2 6 +10 1 2 5 1 3 9 +11 1 2 5 1 9 10 +12 1 2 5 1 10 4 +13 3 2 1 1 12 13 5 1 +14 3 2 1 1 11 14 13 12 +15 3 2 1 1 4 10 14 11 +16 3 2 1 1 13 15 6 5 +17 3 2 1 1 14 16 15 13 +18 3 2 1 1 10 9 16 14 +19 3 2 1 1 15 7 2 6 +20 3 2 1 1 16 8 7 15 +21 3 2 1 1 9 3 8 16 +$EndElements diff --git a/test/3x3_minus.msh b/test/3x3_minus.msh new file mode 100644 index 00000000..b7192593 --- /dev/null +++ b/test/3x3_minus.msh @@ -0,0 +1,38 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +1 +2 10 "fluid" +$EndPhysicalNames +$Nodes +16 +1 0 0 0 +2 1 0 0 +3 1 0.5 0 +4 0 0.5 0 +5 0.3333333333333333 0 0 +6 0.6666666666666666 0 0 +7 1 0.1666666666666667 0 +8 1 0.3333333333333333 0 +9 0.6666666666666667 0.5 0 +10 0.3333333333333334 0.5 0 +11 0 0.3333333333333334 0 +12 0 0.1666666666666667 0 +13 0.3333333333333334 0.1666666666666667 0 +14 0.3333333333333334 0.3333333333333333 0 +15 0.6666666666666667 0.1666666666666667 0 +16 0.6666666666666667 0.3333333333333333 0 +$EndNodes +$Elements +9 +1 3 2 10 1 12 13 5 1 +2 3 2 10 1 14 13 12 11 +3 3 2 10 1 10 14 11 4 +4 3 2 10 1 13 15 6 5 +5 3 2 10 1 14 16 15 13 +6 3 2 10 1 16 14 10 9 +7 3 2 10 1 6 15 7 2 +8 3 2 10 1 15 16 8 7 +9 3 2 10 1 8 16 9 3 +$EndElements diff --git a/test/3x3_twisted.msh b/test/3x3_twisted.msh new file mode 100644 index 00000000..e0607bef --- /dev/null +++ b/test/3x3_twisted.msh @@ -0,0 +1,38 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +1 +2 10 "fluid" +$EndPhysicalNames +$Nodes +16 +1 0 0 0 +2 1 0 0 +3 1 0.5 0 +4 0 0.5 0 +5 0.3333333333333333 0 0 +6 0.6666666666666666 0 0 +7 1 0.1666666666666667 0 +8 1 0.3333333333333333 0 +9 0.6666666666666667 0.5 0 +10 0.3333333333333334 0.5 0 +11 0 0.3333333333333334 0 +12 0 0.1666666666666667 0 +13 0.3333333333333334 0.1666666666666667 0 +14 0.3333333333333334 0.3333333333333333 0 +15 0.6666666666666667 0.1666666666666667 0 +16 0.6666666666666667 0.3333333333333333 0 +$EndNodes +$Elements +9 +1 3 2 10 1 1 5 13 12 +2 3 2 10 1 11 12 13 14 +3 3 2 10 1 4 11 14 10 +4 3 2 10 1 5 6 15 13 +5 3 2 10 1 13 15 16 14 +6 3 2 10 1 9 10 14 16 +7 3 2 10 1 2 7 15 6 +8 3 2 10 1 7 8 16 15 +9 3 2 10 1 3 9 16 8 +$EndElements diff --git a/test/test_mesh.py b/test/test_mesh.py index 22cd615e..8cc0e3ee 100644 --- a/test/test_mesh.py +++ b/test/test_mesh.py @@ -571,12 +571,14 @@ def test_merge_and_map(actx_factory, group_cls, visualize=False): # {{{ element orientation -@pytest.mark.parametrize("case", ["blob", "gh-394"]) +@pytest.mark.parametrize("case", ["blob", "gh-394", "3x3", "3x3_twisted", + "3x3_minus", "3x3_bound"]) def test_element_orientation_via_flipping(case): from meshmode.mesh.io import FileSource, generate_gmsh mesh_order = 3 + meshfile = f"{thisdir}/{case}.msh" if case == "blob": mesh = generate_gmsh( FileSource(str(thisdir / "blob-2d.step")), 2, order=mesh_order, @@ -586,13 +588,62 @@ def test_element_orientation_via_flipping(case): ) elif case == "gh-394": mesh = mio.read_gmsh( - str(thisdir / "gh-394.msh"), + meshfile, + force_ambient_dim=2, + mesh_construction_kwargs={"skip_tests": True}) + elif case == "3x3": # regular ole rectangular 3x3 tensor product els (TPE) + mesh = mio.read_gmsh( + meshfile, + force_ambient_dim=2, + mesh_construction_kwargs={"skip_tests": True}) + elif case == "3x3_twisted": # TPEs, rotated connectivities, all positive + mesh = mio.read_gmsh( + meshfile, + force_ambient_dim=2, + mesh_construction_kwargs={"skip_tests": True}) + elif case == "3x3_minus": # TPEs with negative orientation (clockwise conn) + mesh = mio.read_gmsh( + meshfile, + force_ambient_dim=2, + mesh_construction_kwargs={"skip_tests": True}) + elif case == "3x3_bound": # TPEs (clockwise conn, w/boundaries) + mesh = mio.read_gmsh( + meshfile, force_ambient_dim=2, mesh_construction_kwargs={"skip_tests": True}) else: raise ValueError(f"unknown case: {case}") + boundary_tags = set() + for igrp in range(len(mesh.groups)): + bdry_fagrps = [ + fagrp for fagrp in mesh.facial_adjacency_groups[igrp] + if isinstance(fagrp, BoundaryAdjacencyGroup)] + for bdry_fagrp in bdry_fagrps: + print(f"Boundary tag: {bdry_fagrp.boundary_tag}") + boundary_tags.add(bdry_fagrp.boundary_tag) + mesh_orient = mproc.find_volume_mesh_element_orientations(mesh) + if not (mesh_orient > 0).all(): + logger.info(f"Mesh({meshfile}) is negative, trying to reorient.") + mesh = mio.read_gmsh( + meshfile, + force_ambient_dim=2, + mesh_construction_kwargs={ + "skip_tests": True, + "force_positive_orientation": True}) + + mesh_orient = mproc.find_volume_mesh_element_orientations(mesh) + boundary_tags_reoriented = set() + for igrp in range(len(mesh.groups)): + bdry_fagrps = [ + fagrp for fagrp in mesh.facial_adjacency_groups[igrp] + if isinstance(fagrp, BoundaryAdjacencyGroup)] + for bdry_fagrp in bdry_fagrps: + boundary_tags_reoriented.add(bdry_fagrp.boundary_tag) + + # Make sure rotation doesn't lose boundaries + assert boundary_tags == boundary_tags_reoriented assert (mesh_orient > 0).all()