Skip to content

Full rewrite of meshing for fractured domains#1576

Open
keileg wants to merge 325 commits intodevelopfrom
meshing_rework
Open

Full rewrite of meshing for fractured domains#1576
keileg wants to merge 325 commits intodevelopfrom
meshing_rework

Conversation

@keileg
Copy link
Contributor

@keileg keileg commented Jan 19, 2026

Proposed changes

This PR is a major update of the functionality for meshing of fractured domains. The main changes are:

  1. The fracture network classes have been almost completely rewritten. Identification of intersecting fractures is now handled by Gmsh's backport to open cascade. This also enabled us to implement truly elliptic fractures instead of the polygon approximations used previously.
  2. The algorithm for tuning of mesh size has been improved to allow for more aggressive refinement towards areas of complex geometries.
  3. The test suite for fracture meshing is to a large degree rewritten. The new (parametrized) tests take a more systematic approach to testing. Some old tests that did not fit into this framework have survived.
  4. The tutorial for fracture meshing has been rewritten to fit in the context of the multiphysics models, and to give some ideas about the new meshing framework. This will require some coordination with Update tutorials #1545; I can give an overview of this in person.

Beyond the module fracs, the changes to the code base should be limited. Some tests have been updated, partly to ensure that the tests are not too computationally expensive.

Migration to the new code

While the update should not be API breaking, the practical implications of the mesh size arguments will change. It is recommended to read the new meshing tutorial carefully before migrating to the new functionality.

Please not that fine tuning of the new algorithm may occur over the coming months.

Suggested workflow for review

  1. The natural starting point is the new tutorial.
  2. The API to the meshing, including names of mesh size controls etc. should be critically reviewed.
  3. The inner workings of the meshing algorithm is highly technical (I have tried to limit the technicalities, but the negotiations with Gmsh must be somewhere), so to some degree this is what it is. However, the interface between the components of the meshing could possibly be improved.
  4. Test coverage for fracture meshing.
  5. Check that tests outside the meshing framework have not been molested.
  6. If the functionality can be used by others as part of the review, that would be ideal.

These are my initial thoughts, they may or may not be good.

EDIT 04.02.26
This PR is downstream from #1573. The merging of that PR (when approved) into this branch will impact the following files:

  1. fracture_importer and its tests are fully rewritten in 1573, there is no need to review the changes in this branch.
  2. The tutorial will get an extra section on importing model geometries. This should be a pure add-on compared to the current text.
  3. There will be minor changes in a few other places, but my expectation is this will not be focal points of the present review.

See also todo list below.

What remains or have not been considered

  1. We need coordination with the work on Let multiphysics models import meshes (or geometries) from gmsh #1517. This includes both mentioning of that work in the new tutorial and updates of tests. I suggest we coordinate in person.
  2. There are some tests, mainly of convergence type, that fail due to changes in the underlying mesh. These should be updated, but not until we are sure the meshing algorithm is fixed.
  3. Some of the slow functional tests that are skipped in regular runs are now even slower. As with the previous point, I will update these when things are settled. Also note that some of these can be sped up by considering fewer refinement levels - I am not sure if we should do this in the same go.

EDIT: Todo list when merging in #1573

  • Fix meshing_kwargs in the new load geometry mixin
  • In set_geometry of the load geometry mixin, the spatial dimension can be inferred from the fracture network class introduced in this PR.
  • The possibility to load geometries should be mentioned in the meshing tutorial.
  • Elliptic fractures should be properly used, both in fracture_import and its test.

Types of changes

What types of changes does this PR introduce to PorePy?
Put an x in the boxes that apply.

  • Minor change (e.g., dependency bumps, broken links).
  • Bugfix (non-breaking change which fixes an issue).
  • New feature (non-breaking change which adds functionality).
  • Breaking change (fix or feature that would cause existing functionality to not work as expected).
  • Testing (contribution related to testing of existing or new functionality).
  • Documentation (contribution related to adding, improving, or fixing documentation).
  • Maintenance (e.g., improve logic and performance, remove obsolete code).
  • Other:

Checklist

Put an x in the boxes that apply or explain briefly why the box is not relevant.

  • The documentation is up-to-date.
  • Static typing is included in the update.
  • This PR does not duplicate existing functionality.
  • The update is covered by the test suite (including tests added in the PR).
  • If new skipped tests have been introduced in this PR, pytest was run with the --run-skipped flag.

keileg and others added 30 commits October 13, 2025 07:32
Merge work on gmsh object processing into main branch for new gmsh updates
New implementations, based on Gmsh OpenCASCADE to come
This is a helper module which should not be used externally
Non-convex domains are no longer supported.
The tacit support for non-convex domains is hereby dropped
Raise error messages if relevant.
Also maintenance and documentation
More systematic testing of possible scenarios. Also parametrization.
@review-notebook-app
Copy link

review-notebook-app bot commented Feb 2, 2026

View / edit / reply to this conversation on ReviewNB

IvarStefansson commented on 2026-02-02T09:04:22Z
----------------------------------------------------------------

Is the writing of all this logging default behaviour? I find it too verbose.


keileg commented on 2026-02-02T11:39:49Z
----------------------------------------------------------------

Should we make the gmsh verbosity level a user parameter?

jwboth commented on 2026-02-02T12:22:31Z
----------------------------------------------------------------

I don't know why, but I like it. So if at least there is a way of turning on info logging from gmsh, OK with me.

@review-notebook-app
Copy link

review-notebook-app bot commented Feb 2, 2026

View / edit / reply to this conversation on ReviewNB

IvarStefansson commented on 2026-02-02T09:04:22Z
----------------------------------------------------------------

Could we simplify to "It may be desirable to target the mesh resolution towards the vicinity of the fractures." ?

I think the sentence in parentheses is more than important enough to warrant a promotion out of the parentheses!


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 2, 2026

View / edit / reply to this conversation on ReviewNB

IvarStefansson commented on 2026-02-02T09:04:23Z
----------------------------------------------------------------

This seems to be a fundamental explanation that should come earlier.


jwboth commented on 2026-02-02T12:45:24Z
----------------------------------------------------------------

I agree. See my edits, that I keep for chronological understanding of my thoughts. This of course answers some of my questions.

@review-notebook-app
Copy link

review-notebook-app bot commented Feb 2, 2026

View / edit / reply to this conversation on ReviewNB

IvarStefansson commented on 2026-02-02T09:04:24Z
----------------------------------------------------------------

I think we should add a section the functionality for reading from file to manually tune grids and/or saving run time for complex geometries.


keileg commented on 2026-02-02T11:38:57Z
----------------------------------------------------------------

Agree, will do when merging the changes from that project into this PR.

Copy link
Contributor Author

keileg commented Feb 2, 2026

Agree, will do when merging the changes from that project into this PR.


View entire conversation on ReviewNB

Copy link
Contributor Author

keileg commented Feb 2, 2026

Should we make the gmsh verbosity level a user parameter?


View entire conversation on ReviewNB

Copy link
Contributor

jwboth commented Feb 2, 2026

I don't know why, but I like it. So if at least there is a way of turning on info logging from gmsh, OK with me.


View entire conversation on ReviewNB

Copy link
Contributor

jwboth commented Feb 2, 2026

I agree. See my edits, that I keep for chronological understanding of my thoughts. This of course answers some of my questions.


View entire conversation on ReviewNB

@review-notebook-app
Copy link

review-notebook-app bot commented Feb 2, 2026

View / edit / reply to this conversation on ReviewNB

jwboth commented on 2026-02-02T12:50:12Z
----------------------------------------------------------------

Use hyperlink to link to single phase flow tutorial.

It is somewhat unfortunate to communicate two ways of controlling meshing. Somewhat OK. But maybe consolidate and use the same in both tutorials?


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 2, 2026

View / edit / reply to this conversation on ReviewNB

jwboth commented on 2026-02-02T12:50:13Z
----------------------------------------------------------------

I think this is better than before! Of course, it remains a black box what the effective cell size will be, but I am OK with that. My understanding of the concrete example is. With refinement_threshold=1.0 it may have happened that some strange coarsening would take place with skewed cells with two objects just twice the target cell size apart from each other. With refinement_threshold=2.0 one manages to dictate a uniformly fine mesh with nice cell shapes in between the two objects.

Edit: It somehow does not surprise me based on my initial understanding, but maybe this is not a good example if the parameters are changed and the mesh looks the same. Better would be a visual support of the description of the target effect.


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 2, 2026

View / edit / reply to this conversation on ReviewNB

jwboth commented on 2026-02-02T12:50:14Z
----------------------------------------------------------------

maybe add some extra interpretation? Something along the lines of aiming to control the mesh within the hull of the fracture network? Just to make sure that one targets here the general mesh and delaying the transition towards the cell_size_boundary.OK also without, if this would be too much.

Edit: Looking at the subsequent plot, I am questioning what I just wrote. Which parameter was important for the control of the "hull size"? Based on the above edit, I am not sure whether it was the refinement_buffer or refinement_threshold.The latter would make more sense. Maybe flip the visual demonstration of the two keywords?


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 2, 2026

View / edit / reply to this conversation on ReviewNB

jwboth commented on 2026-02-02T12:50:15Z
----------------------------------------------------------------

As intersting for some applications as fracture propagation or fracture tip focus, maybe stress that it is neither possible to control the cell size close to fracture tips?


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 2, 2026

View / edit / reply to this conversation on ReviewNB

jwboth commented on 2026-02-02T12:50:31Z
----------------------------------------------------------------

Use hyperlink to link to single phase flow tutorial.

It is somewhat unfortunate to communicate two ways of controlling meshing. Somewhat OK. But maybe consolidate and use the same in both tutorials?


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 2, 2026

View / edit / reply to this conversation on ReviewNB

jwboth commented on 2026-02-02T12:50:32Z
----------------------------------------------------------------

I think this is better than before! Of course, it remains a black box what the effective cell size will be, but I am OK with that. My understanding of the concrete example is. With refinement_threshold=1.0 it may have happened that some strange coarsening would take place with skewed cells with two objects just twice the target cell size apart from each other. With refinement_threshold=2.0 one manages to dictate a uniformly fine mesh with nice cell shapes in between the two objects.

Edit: It somehow does not surprise me based on my initial understanding, but maybe this is not a good example if the parameters are changed and the mesh looks the same. Better would be a visual support of the description of the target effect.


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 2, 2026

View / edit / reply to this conversation on ReviewNB

jwboth commented on 2026-02-02T12:50:32Z
----------------------------------------------------------------

maybe add some extra interpretation? Something along the lines of aiming to control the mesh within the hull of the fracture network? Just to make sure that one targets here the general mesh and delaying the transition towards the cell_size_boundary.OK also without, if this would be too much.

Edit: Looking at the subsequent plot, I am questioning what I just wrote. Which parameter was important for the control of the "hull size"? Based on the above edit, I am not sure whether it was the refinement_buffer or refinement_threshold.The latter would make more sense. Maybe flip the visual demonstration of the two keywords?


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 2, 2026

View / edit / reply to this conversation on ReviewNB

jwboth commented on 2026-02-02T12:50:33Z
----------------------------------------------------------------

As intersting for some applications as fracture propagation or fracture tip focus, maybe stress that it is neither possible to control the cell size close to fracture tips?


Copy link
Contributor

@IvarStefansson IvarStefansson left a comment

Choose a reason for hiding this comment

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

Some effort! Only relatively minor comments, but quite a few of them. Let's discuss where indicated. Did not consider tests in this iteration.

gmsh.initialize()

domain_tag = network.domain_to_gmsh()
fracture_tags = network.fractures_to_gmsh()
Copy link
Contributor

Choose a reason for hiding this comment

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

Am I being unreasonable if I claim this is somewhat sparsely documented?

@@ -0,0 +1,150 @@
"""Contains classes representing two-dimensional fractures in 3D."""
Copy link
Contributor

Choose a reason for hiding this comment

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

Is there a reason the class is ellipTIC and the file ellipSE?

)

def copy(self) -> Fracture:
"""Return a copy of the fracture with the current vertices.
Copy link
Contributor

Choose a reason for hiding this comment

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

Is "current vertices" accurate?

major_axis_angle=self.major_axis_angle,
strike_angle=self.strike_angle,
dip_angle=self.dip_angle,
index=self.index,
Copy link
Contributor

Choose a reason for hiding this comment

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

What does this imply about the interpretation of the index?

order/sorting of vertices has to be implemented explicitly.
"""

def set_index(self, index: int) -> None:
Copy link
Contributor

Choose a reason for hiding this comment

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

See comment in ellipse_fracture.py. Should we specify intended use/interpretation?


# For a boundary line to be an intersection, it must be shared by at least two
# fractures.
num_lines_occ = np.bincount(np.abs(bnd_lines).astype(int))
Copy link
Contributor

Choose a reason for hiding this comment

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

Why abs?


### Get hold of lines representing fractures and boundaries.
domain_entities = gmsh.model.get_entities(nd)
# TODO: If there is more than one domain entity (the domain is split into parts
Copy link
Contributor

Choose a reason for hiding this comment

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

Please resolve.

# and assign the minimum mesh size among all occurrences of the point.
self._uniquify_mesh_size_dictionary(mesh_size_points)

# For lines that with no extra mesh size control points, assign an empty list.
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
# For lines that with no extra mesh size control points, assign an empty list.
# For lines with no extra mesh size control points, assign an empty list.

@@ -2721,6 +1147,9 @@ def to_csv(self, file_name: Path, domain: Optional[pp.Domain] = None) -> None:

Domain specification.

Raises:
Copy link
Contributor

Choose a reason for hiding this comment

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

I don't know if it is a priority, but we could create an issue about implementing it as exporting the parameters defining an elliptic fracture.


"""
pts = self.pts.T
num_pts = len(pts)
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
num_pts = len(pts)
num_pts = pts.shape[0]

Copy link
Contributor

@IvarStefansson IvarStefansson left a comment

Choose a reason for hiding this comment

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

Some effort! Only relatively minor comments, but quite a few of them. Let's discuss where indicated. Did not consider tests in this iteration.

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.

6 participants