Created
November 19, 2025 13:59
-
-
Save jorgensd/4e0eda5b687bea6bfda3e7ef9225ac43 to your computer and use it in GitHub Desktop.
Compute curvature at quadrature points
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
| # Author Jorgen S. Dokken | |
| # SPDX license identifier: MIT | |
| import gmsh | |
| import scifem | |
| from mpi4py import MPI | |
| import basix.ufl | |
| import dolfinx.fem.petsc | |
| import ufl | |
| import numpy as np | |
| def move_to_facet_quadrature(ufl_expr, mesh, sub_facets, scheme="default", degree=6): | |
| fdim = mesh.topology.dim - 1 | |
| # Create submesh | |
| bndry_mesh, entity_map, _, _ = dolfinx.mesh.create_submesh(mesh, fdim, sub_facets) | |
| sub_top = entity_map.sub_topology | |
| assert isinstance(sub_top, dolfinx.mesh.Topology) | |
| sub_map = sub_top.index_map(entity_map.dim) | |
| indices = np.arange(sub_map.size_local + sub_map.num_ghosts, dtype=np.int32) | |
| parent_facets = entity_map.sub_topology_to_topology(indices, inverse=False) | |
| # Create quadrature space on submesh | |
| q_el = basix.ufl.quadrature_element( | |
| bndry_mesh.basix_cell(), ufl_expr.ufl_shape, scheme, degree | |
| ) | |
| Q = dolfinx.fem.functionspace(bndry_mesh, q_el) | |
| # Compute where to evaluate expression per submesh cell | |
| integration_entities = dolfinx.fem.compute_integration_domains( | |
| dolfinx.fem.IntegralType.exterior_facet, mesh.topology, parent_facets | |
| ) | |
| compiled_expr = dolfinx.fem.Expression(ufl_expr, Q.element.interpolation_points) | |
| # Evaluate expression | |
| q = dolfinx.fem.Function(Q) | |
| q.x.array[:] = compiled_expr.eval( | |
| mesh, integration_entities.reshape(-1, 2) | |
| ).reshape(-1) | |
| return q | |
| gmsh.initialize() | |
| membrane = gmsh.model.occ.addDisk(0, 0, 0, 3, 1) | |
| gmsh.model.occ.synchronize() | |
| gdim = 2 | |
| gmsh.model.addPhysicalGroup(gdim, [membrane], 1) | |
| gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.2) | |
| gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.2) | |
| gmsh.model.mesh.generate(gdim) | |
| gmsh.model.mesh.setOrder(2) | |
| gmsh_model_rank = 0 | |
| mesh_comm = MPI.COMM_WORLD | |
| mesh_data = dolfinx.io.gmsh.model_to_mesh( | |
| gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim | |
| ) | |
| gmsh.finalize() | |
| assert mesh_data.cell_tags is not None | |
| cell_markers = mesh_data.cell_tags | |
| domain = mesh_data.mesh | |
| n = ufl.FacetNormal(domain) | |
| kappa = ufl.div(n) | |
| domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim) | |
| exterior_boundaries = dolfinx.mesh.exterior_facet_indices(domain.topology) | |
| kappa = move_to_facet_quadrature(kappa, domain, exterior_boundaries) | |
| with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "domain.xdmf", "w") as xdmf: | |
| xdmf.write_mesh(domain) | |
| scifem.xdmf.create_pointcloud("kappa.xdmf", [kappa]) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment