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.

Sign up for GitHub

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

Jump to bottom

add checks of cell tag consistency in gmshio #3342

Open
wants to merge 15 commits into
base: main
Choose a base branch
from
Open

add checks of cell tag consistency in gmshio #3342

Changes from 1 commit
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
c1e8103
add checks of cell tag consistency in gmshio
pierricmora Aug 8, 2024
35531da
RuntimeError rather than assert
pierricmora Aug 8, 2024
8390cd4
Merge branch 'main' into PR_check_celltags_in_gmshio
jorgensd Aug 14, 2024
de4b5f6
Merge branch 'main' into PR_check_celltags_in_gmshio
jorgensd Aug 27, 2024
6775c8a
Merge branch 'FEniCS:main' into PR_check_celltags_in_gmshio
pierricmora Aug 29, 2024
97a2023
E501
pierricmora Aug 29, 2024
37436e2
E713
pierricmora Aug 29, 2024
24d087f
Merge branch 'main' into PR_check_celltags_in_gmshio
garth-wells Sep 10, 2024
56f371c
Merge branch 'main' into PR_check_celltags_in_gmshio
jorgensd Sep 15, 2024
ffefa65
Merge branch 'main' into PR_check_celltags_in_gmshio
garth-wells Sep 19, 2024
9b5fb55
Merge branch 'main' into PR_check_celltags_in_gmshio
jorgensd Sep 24, 2024
3dafc47
ruff
pierricmora Sep 24, 2024
4919588
Merge branch 'main' into PR_check_celltags_in_gmshio
pierricmora Sep 24, 2024
4e8ad33
bug fix: read the appropriate key of topologies
pierricmora Sep 24, 2024
815ee7f
typo
pierricmora Sep 24, 2024
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
25 changes: 24 additions & 1 deletion python/dolfinx/io/gmshio.py
View file
Open in desktop
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,14 @@ 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 +266,20 @@ 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?

assert _d in topologies.keys(), 'All cells are expected to be tagged once; none found'
Copy link
Member

Choose a reason for hiding this comment

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

As these are user-facing errors they should raise Exceptions of the appropriate type directly, probably RuntimeError, rather than using assert.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Done.

Suggestion: could the first RuntimeError (no cell is tagged) be converted into a RuntimeWarning by defaulting a uniform tag for all cells?

Copy link
Member

Choose a reason for hiding this comment

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

Currently most warnings get swallowed by FFCx (FEniCS/ffcx#712) so I am not sure this is a good idea

_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
nbcells = len(_elementTags[0])
nbcells_tagged = len(topologies[_d]["entity_tags"])
assert nbcells == nbcells_tagged, \
f'All cells are expected to be tagged once; found: {nbcells_tagged}, expected: {nbcells}'
nbcells_tagged_once = len(np.unique(topologies[_d]['entity_tags']))
assert nbcells_tagged == nbcells_tagged_once, \
'All cells are expected to be tagged once; found duplicates'

# 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