Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add checks of cell tag consistency in gmshio #3342

Open
wants to merge 15 commits into
base: main
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 32 additions & 1 deletion python/dolfinx/io/gmshio.py
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,9 @@ def extract_topology_and_markers(model, name: typing.Optional[str] = None):
# Create marker array of length of number of tagged cells
marker = np.full_like(entity_tags[0], tag)

# Keep track of entity tags to check tag consistency
entity_tags = entity_tags[0]

# Group element topology and markers of the same entity type
entity_type = entity_types[0]
if entity_type in topologies.keys():
Expand All @@ -161,8 +164,15 @@ def extract_topology_and_markers(model, name: typing.Optional[str] = None):
topologies[entity_type]["cell_data"] = np.hstack(
[topologies[entity_type]["cell_data"], marker]
)
topologies[entity_type]["entity_tags"] = np.hstack(
[topologies[entity_type]["entity_tags"], entity_tags]
)
else:
topologies[entity_type] = {"topology": topology, "cell_data": marker}
topologies[entity_type] = {
"topology": topology,
"cell_data": marker,
"entity_tags": entity_tags,
}

return topologies

Expand Down Expand Up @@ -257,6 +267,27 @@ def model_to_mesh(
# Sort elements by ascending dimension
perm_sort = np.argsort(cell_dimensions)

# Check that all cells are tagged once
_d = model.getDimension()
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could be used to remove gdim from parameters. What about this simplification? @jorgensd

if comm.rank == rank:
    gdim = model.getDimension()
    gdim = comm.bcast(gdim, root=rank)

# ...

else:
    gdim = comm.bcast(None, root=rank)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure that is safe, as you might want: A flat manifold in some plane, i.e.

In [1]: import gmsh

In [2]: gmsh.initialize()

In [3]: idx = gmsh.model.occ.add_rectangle(0,0,0.1,1,1)

In [4]: gmsh.model.occ.synchronize()

In [5]: gmsh.model.add_physical_group(2,[idx])
Out[5]: 1

In [6]: gmsh.model.mesh.generate(2)
Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 30%] Meshing curve 2 (Line)
Info    : [ 60%] Meshing curve 3 (Line)
Info    : [ 80%] Meshing curve 4 (Line)
Info    : Done meshing 1D (Wall 0.000247945s, CPU 0.000433s)
Info    : Meshing 2D...
Info    : Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.00361243s, CPU 0.003612s)
Info    : 98 nodes 198 elements

In [7]: gmsh.model.getDimension()
Out[7]: 2

In [8]: gmsh.model.mesh.generate(3)
Info    : Meshing 3D...
Info    : Done meshing 3D (Wall 2.2143e-05s, CPU 2.1e-05s)
Info    : 98 nodes 198 elements

In [9]: gmsh.model.getDimension()
Out[9]: 2

ref:

In [10]: gmsh.model.occ.add_rectangle?
Signature: gmsh.model.occ.add_rectangle(x, y, z, dx, dy, tag=-1, roundedRadius=0.0)
Docstring:
gmsh.model.occ.addRectangle(x, y, z, dx, dy, tag=-1, roundedRadius=0.)

Add a rectangle in the OpenCASCADE CAD representation, with lower left
corner at (`x', `y', `z') and upper right corner at (`x' + `dx', `y' +
`dy', `z'). If `tag' is positive, set the tag explicitly; otherwise a new
tag is selected automatically. Round the corners if `roundedRadius' is
nonzero. Return the tag of the rectangle.

Return an integer.

Types:
- `x': double
- `y': double
- `z': double
- `dx': double
- `dy': double
- `tag': integer
- `roundedRadius': double

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, good point. I'm not sure to understand in details what should be the desired behaviour here -- I guess you would expect 3 as output, and another method that does not exist but which could be named gmsh.model.getTopologicalDimension() to return 2, am I right?
I think I am done with my code. Could you have a review?

_elementTypes, _elementTags, _nodeTags = model.mesh.getElements(dim=_d, tag=-1)
# assert only one type of elements
# assert len(_elementTypes) == 1 # NOTE: already checked in extract_topology_and_markers
_elementType_dim = _elementTypes[0]
if _elementType_dim not in topologies.keys():
raise RuntimeError("All cells are expected to be tagged once; none found")
nbcells = len(_elementTags[0])
nbcells_tagged = len(topologies[_elementType_dim]["entity_tags"])
if nbcells != nbcells_tagged:
e = (
"All cells are expected to be tagged once;"
f"found: {nbcells_tagged}, expected: {nbcells}"
)
raise RuntimeError(e)
nbcells_tagged_once = len(np.unique(topologies[_elementType_dim]["entity_tags"]))
if nbcells_tagged != nbcells_tagged_once:
e = "All cells are expected to be tagged once; found duplicates"
raise RuntimeError(e)

# Broadcast cell type data and geometric dimension
cell_id = cell_information[perm_sort[-1]]["id"]
tdim = cell_information[perm_sort[-1]]["dim"]
Expand Down
Loading