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

Conversation

pierricmora
Copy link
Contributor

Relates to issue #3332

Add checks to the gmshio interface to ensure that all cells are tagged once, to avoid undesired behaviour or crashes.

Compared to the discussion in issue #3332, the impact on the code is slightly more as not only model_to_mesh but also extract_topology_and_markers have been modified. This is because storing entity_tags in extract_topology_and_markers enables to call the relatively expensive np.unique on a smaller array. Then, the overcost due to these checks is below 5%, compared to 15% - 10% reported before.

@@ -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()
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

@pierricmora
Copy link
Contributor Author

I corrected two formatting errors but I still get Lint failed. This time the error is not very verbose... Any idea? When I do ruff check and ruff format in my computer I get no error.

@mscroggs mscroggs added this to the 0.9.0 milestone Sep 3, 2024
@@ -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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants