Skip to content

Meshing fails for T-type intersection #1694

Description

@keileg

When trying to do meshing with the following example (slightly modified from the single phase flow tutorial):

import porepy as pp
import numpy as np
from porepy.models.fluid_mass_balance import SinglePhaseFlow
from porepy.applications.md_grids.domains import nd_cube_domain

class ModifiedGeometry:
    def set_domain(self) -> None:
        size = self.units.convert_units(2, "m")
        self._domain = nd_cube_domain(2, size)

    def set_fractures(self) -> None:
        fractures = [
            pp.LineFracture(np.array([[0.2, 1.8], [0.2, 1.8]])),
            pp.LineFracture(np.array([[1.0, 0.3], [1.0, 0.3]])),
        ]
        self._fractures = fractures

    def grid_type(self) -> str:
        return self.params.get("grid_type", "simplex")

    def meshing_arguments(self) -> dict:
        cell_size = self.units.convert_units(0.25, "m")
        mesh_args: dict[str, float] = {"cell_size": cell_size}
        return mesh_args

class SinglePhaseFlowGeometry(ModifiedGeometry, SinglePhaseFlow):
    ...

model_params = {}
model = SinglePhaseFlowGeometry(model_params)
pp.ModelRunner(model).run()
print(model.mdg)

the output shows an MDG with a single fracture grid (representing the first fracture) instead of two. Also, there are 2 0d grids instead of the expected one grid , with the extra 0d grid representing the endpoint of the skipped frature.

I have traced the behavior down to the fragmentation of fractures by Gmsh's OpenCASCADE kernel (called in the method _fragmen_fractures in class FractureNetwork ): This falsely identifies the line of the second fracture as being part of the first fracture. With this misidentification, the rest of the code stand no chance - the second fracture will get lost. The extra 0d grid is strange, but should disappear if the missing fracture is somehow recovered.

Right now, it is tempting to blame OpenCASCADE or the Gmsh interface to it, but that is too early to tell. I am not sure what is the next step; possibly to try to dig inside Gmsh to understand what happens.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions