Last active
August 20, 2024 09:17
-
-
Save jorgensd/be065cc9d333f73892da8ebaa013651f to your computer and use it in GitHub Desktop.
Submesh `ds` integral
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # Illustration of submesh facet integral with parent coefficient | |
| # Author: Jørgen S. Dokken | |
| # SPDX-License-Identifier: MIT | |
| def transfer_meshtags_to_submesh( | |
| mesh, entity_tag, submesh, sub_vertex_to_parent, sub_cell_to_parent | |
| ): | |
| """ | |
| Transfer a meshtag from a parent mesh to a sub-mesh. | |
| Function from: https://gist.github.com/jorgensd/8a5c32f491195e838f5863ca88b27bce#file-problem-py-L289-L291 | |
| SPDX-License-Identifier: MIT | |
| Author: Jø¶gen S. Dokken | |
| """ | |
| tdim = mesh.topology.dim | |
| cell_imap = mesh.topology.index_map(tdim) | |
| num_cells = cell_imap.size_local + cell_imap.num_ghosts | |
| mesh_to_submesh = np.full(num_cells, -1) | |
| mesh_to_submesh[sub_cell_to_parent] = np.arange( | |
| len(sub_cell_to_parent), dtype=np.int32 | |
| ) | |
| sub_vertex_to_parent = np.asarray(sub_vertex_to_parent) | |
| submesh.topology.create_connectivity(entity_tag.dim, 0) | |
| num_child_entities = ( | |
| submesh.topology.index_map(entity_tag.dim).size_local | |
| + submesh.topology.index_map(entity_tag.dim).num_ghosts | |
| ) | |
| submesh.topology.create_connectivity(submesh.topology.dim, entity_tag.dim) | |
| c_c_to_e = submesh.topology.connectivity(submesh.topology.dim, entity_tag.dim) | |
| c_e_to_v = submesh.topology.connectivity(entity_tag.dim, 0) | |
| child_markers = np.full(num_child_entities, 0, dtype=np.int32) | |
| mesh.topology.create_connectivity(entity_tag.dim, 0) | |
| mesh.topology.create_connectivity(entity_tag.dim, mesh.topology.dim) | |
| p_f_to_v = mesh.topology.connectivity(entity_tag.dim, 0) | |
| p_f_to_c = mesh.topology.connectivity(entity_tag.dim, mesh.topology.dim) | |
| sub_to_parent_entity_map = np.full(num_child_entities, -1, dtype=np.int32) | |
| for facet, value in zip(entity_tag.indices, entity_tag.values): | |
| facet_found = False | |
| for cell in p_f_to_c.links(facet): | |
| if facet_found: | |
| break | |
| if (child_cell := mesh_to_submesh[cell]) != -1: | |
| for child_facet in c_c_to_e.links(child_cell): | |
| child_vertices = c_e_to_v.links(child_facet) | |
| child_vertices_as_parent = sub_vertex_to_parent[child_vertices] | |
| is_facet = np.isin( | |
| child_vertices_as_parent, p_f_to_v.links(facet) | |
| ).all() | |
| if is_facet: | |
| child_markers[child_facet] = value | |
| facet_found = True | |
| sub_to_parent_entity_map[child_facet] = facet | |
| tags = dolfinx.mesh.meshtags( | |
| submesh, | |
| entity_tag.dim, | |
| np.arange(num_child_entities, dtype=np.int32), | |
| child_markers, | |
| ) | |
| tags.name = entity_tag.name | |
| return tags, sub_to_parent_entity_map | |
| from mpi4py import MPI | |
| import dolfinx | |
| import numpy as np | |
| import ufl | |
| def half_mesh(x): | |
| return x[0] < 0.5 + 1e-14 | |
| def left_boundary(x): | |
| return np.isclose(x[0], 0.0) | |
| def parallel_assemble_scalar(ufl_form): | |
| compiled_form = dolfinx.fem.form(ufl_form) | |
| comm = compiled_form.mesh.comm | |
| local_scalar = dolfinx.fem.assemble_scalar(compiled_form) | |
| return comm.allreduce(local_scalar, op=MPI.SUM) | |
| mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, dolfinx.cpp.mesh.CellType.triangle) | |
| # Create submesh | |
| tdim = mesh.topology.dim | |
| submesh_cells = dolfinx.mesh.locate_entities(mesh, tdim, half_mesh) | |
| submesh, submesh_c_map, submesh_v_map, _ = dolfinx.mesh.create_submesh(mesh, tdim, submesh_cells) | |
| # Create meshtag for facets on parent mesh | |
| facet_map = mesh.topology.index_map(tdim-1) | |
| num_facets_local = facet_map.size_local + facet_map.num_ghosts | |
| tagged_facets = dolfinx.mesh.locate_entities_boundary(mesh, tdim-1, left_boundary) | |
| facet_values = np.full(num_facets_local, 1, dtype=np.int32) | |
| facet_values[tagged_facets] = 2 | |
| ft = dolfinx.mesh.meshtags(mesh, tdim-1, np.arange(num_facets_local, dtype=np.int32), facet_values) | |
| # Transfer meshtags to submesh | |
| ft_child, facet_sub_to_parent = transfer_meshtags_to_submesh(mesh, ft, submesh, submesh_v_map, submesh_c_map) | |
| V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1)) | |
| u = dolfinx.fem.Function(V) | |
| u.interpolate(lambda x: -7*x[1]) | |
| Q = dolfinx.fem.functionspace(submesh, ("Lagrange", 1)) | |
| q = dolfinx.fem.Function(Q) | |
| q.interpolate(lambda x: 11*x[1]) | |
| # Create facet integral on submesh | |
| ds_sub = ufl.Measure("ds", domain=submesh, subdomain_data=ft_child) | |
| entity_map_sub = {mesh: submesh_c_map} | |
| J_sub = ufl.inner(u, q) * ds_sub(2) | |
| J_sub_compiled = dolfinx.fem.form(J_sub, entity_maps=entity_map_sub) | |
| val_sub = parallel_assemble_scalar(J_sub_compiled) | |
| if MPI.COMM_WORLD.rank == 0: | |
| print(val_sub) | |
| # Create facet integral on parent mesh | |
| ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft) | |
| parent_c_map = mesh.topology.index_map(tdim) | |
| parent_to_sub = np.full(parent_c_map.size_local + parent_c_map.num_ghosts, -1, dtype=np.int32) | |
| parent_to_sub[submesh_c_map] = np.arange(len(submesh_c_map), dtype=np.int32) | |
| entity_map = {submesh: parent_to_sub} | |
| J = ufl.inner(u, q) * ds(2) | |
| J_compiled = dolfinx.fem.form(J, entity_maps=entity_map) | |
| val = parallel_assemble_scalar(J_compiled) | |
| if MPI.COMM_WORLD.rank == 0: | |
| print(val) | |
| # Reference solution | |
| q_p = dolfinx.fem.Function(V) | |
| q_p.interpolate(lambda x: 11*x[1]) | |
| J_ref = ufl.inner(u, q_p) * ds(2) | |
| J_ref_compiled = dolfinx.fem.form(J_ref) | |
| val_ref = parallel_assemble_scalar(J_ref_compiled) | |
| if MPI.COMM_WORLD.rank == 0: | |
| print(val_ref) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment