Hi all,
I just upgraded dolfinx from version 0.5.2 to version 0.8.0. I tried to relaunch the code provided to me in example by @dokken (Get boundary points from custom cell tags)
from shapely import Polygon, Point
import numpy as np
import gmsh
from mpi4py import MPI
from dolfinx.io.gmshio import model_to_mesh
import dolfinx
import dolfinx.plot as plot
import pyvista
pyvista.set_jupyter_backend("static")
pyvista.start_xvfb()
rank = MPI.COMM_WORLD.rank
rect1 = np.array([[1, 1], [4, 1], [4, 4], [1, 4], [1, 1]], dtype=np.float64)
rect2 = np.array([[0, 0], [10, 0], [10, 5], [0, 5], [0, 0]], dtype=np.float64)
gmsh.initialize()
rect1_tags = []
rect2_tags = []
for point1, point2 in zip(rect1, rect2):
p1 = gmsh.model.occ.addPoint(point1[0], point1[1], 0, 1)
p2 = gmsh.model.occ.addPoint(point2[0], point2[1], 0, 1)
rect1_tags.append(p1)
rect2_tags.append(p2)
gmsh.model.occ.synchronize()
rect1_wire = []
rect2_wire = []
for ii in range(len(rect1_tags)):
if ii > 0:
l1 = gmsh.model.occ.addLine(rect1_tags[ii - 1], rect1_tags[ii])
l2 = gmsh.model.occ.addLine(rect2_tags[ii - 1], rect2_tags[ii])
rect1_wire.append(l1)
rect2_wire.append(l2)
gmsh.model.occ.synchronize()
gmsh.model.occ.addCurveLoop(rect1_wire, 1)
gmsh.model.occ.addCurveLoop(rect2_wire, 2)
gmsh.model.occ.synchronize()
gmsh.model.occ.addPlaneSurface([1], 1)
gmsh.model.occ.addPlaneSurface([1, 2], 2)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], 1)
gmsh.model.addPhysicalGroup(2, [2], 2)
gmsh.model.occ.synchronize()
# generate 2D mesh
gmsh.option.setNumber("Mesh.Algorithm", 2)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.optimize("Netgen")
gmsh.model.occ.synchronize()
gmsh.write("mesh.geo_unrolled")
model_rank = 0
mesh_comm = MPI.COMM_WORLD
cell_partitioner = dolfinx.mesh.create_cell_partitioner(dolfinx.mesh.GhostMode.shared_facet)
mesh, ct, ft = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2, partitioner=cell_partitioner)
gmsh.finalize()
def create_meshtags(
mesh: dolfinx.mesh.Mesh, poly: Polygon
) -> dolfinx.cpp.mesh.MeshTags_int32:
cells = mesh.topology.connectivity(mesh.topology.dim, 0)
id_cells_inpoly = []
for ii in range(len(cells)):
triangle = Polygon(mesh.geometry.x[cells.links(ii)])
centroid = triangle.centroid
if poly.contains(Point(centroid)):
id_cells_inpoly.append(ii)
values = np.zeros([len(cells), 1], dtype=np.int32)
values[id_cells_inpoly] = 1
meshtags = dolfinx.mesh.meshtags_from_entities(
mesh, mesh.topology.dim, cells, values
)
return meshtags
def plot_meshtags(
mesh: dolfinx.mesh.Mesh,
meshtags: dolfinx.cpp.mesh.MeshTags_int32,
filename: str = "meshtags.png",
):
# Create VTK mesh
cells, types, x = plot.create_vtk_mesh(mesh)
grid = pyvista.UnstructuredGrid(cells, types, x)
# Attach the cells tag data to the pyvita grid
grid.cell_data["Marker"] = meshtags.values
grid.set_active_scalars("Marker")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, show_scalar_bar=False)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
plotter.screenshot(filename)
# plot original cell tags
plot_meshtags(mesh, ct)
# create custom cell tags
ct = create_meshtags(mesh, Polygon(rect1))
# plot custom cell tags
plot_meshtags(mesh, ct)
cell_map = mesh.topology.index_map(mesh.topology.dim)
vertex_map = mesh.topology.index_map(0)
mesh.topology.create_entities(mesh.topology.dim-1)
facet_map = mesh.topology.index_map(mesh.topology.dim-1)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
f_to_c = mesh.topology.connectivity(mesh.topology.dim-1, mesh.topology.dim)
f_to_v = mesh.topology.connectivity(mesh.topology.dim-1, 0)
assert len(ct.indices) == cell_map.size_local + cell_map.num_ghosts
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
vertex_marker = np.zeros(vertex_map.size_local + vertex_map.num_ghosts, dtype=np.int32)
for facet in range(facet_map.size_local+ facet_map.num_ghosts):
cells = f_to_c.links(facet)
cell_markers = ct.values[cells]
# If exterior facet and connected cell is marked with 1
if 1 in cell_markers and facet in exterior_facets:
vertices = f_to_v.links(facet)
vertex_marker[vertices] = 1
# If on interface between the two domains
if 1 in cell_markers and 0 in cell_markers:
vertices = f_to_v.links(facet)
vertex_marker[vertices] = 1
vt = dolfinx.mesh.meshtags(mesh, 0, np.arange(len(vertex_marker), dtype=np.int32), vertex_marker)
with dolfinx.io.XDMFFile(mesh.comm, "vertices.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
mesh.topology.create_connectivity(0, mesh.topology.dim)
xdmf.write_meshtags(vt)
However, now I get the following error:
meshtags = dolfinx.mesh.meshtags_from_entities(
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
line 534, in meshtags_from_entities
return MeshTags(_cpp.mesh.create_meshtags(mesh.topology, dim, entities, values))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: create_meshtags(): incompatible function arguments. The following argument types are supported:
1. create_meshtags(arg0: dolfinx.cpp.mesh.Topology, arg1: int, arg2: dolfinx.cpp.graph.AdjacencyList_int32, arg3: ndarray[dtype=int8, writable=False, shape=(*), order='C'], /) -> dolfinx.cpp.mesh.MeshTags_int8
2. create_meshtags(arg0: dolfinx.cpp.mesh.Topology, arg1: int, arg2: dolfinx.cpp.graph.AdjacencyList_int32, arg3: ndarray[dtype=int32, writable=False, shape=(*), order='C'], /) -> dolfinx.cpp.mesh.MeshTags_int32
3. create_meshtags(arg0: dolfinx.cpp.mesh.Topology, arg1: int, arg2: dolfinx.cpp.graph.AdjacencyList_int32, arg3: ndarray[dtype=int64, writable=False, shape=(*), order='C'], /) -> dolfinx.cpp.mesh.MeshTags_int64
4. create_meshtags(arg0: dolfinx.cpp.mesh.Topology, arg1: int, arg2: dolfinx.cpp.graph.AdjacencyList_int32, arg3: ndarray[dtype=float64, writable=False, shape=(*), order='C'], /) -> dolfinx.cpp.mesh.MeshTags_float64
Invoked with types: dolfinx.cpp.mesh.Topology, int, dolfinx.cpp.graph.AdjacencyList_int32, ndarray
Can anyone help me? Thanks in advance.