usnistgov / fipy

FiPy is a Finite Volume PDE solver written in Python
http://pages.nist.gov/fipy/en/latest
Other
514 stars 150 forks source link

PermissionError preventing Gmsh2D loading *.msh files. #884

Open shinstra opened 2 years ago

shinstra commented 2 years ago

Hi,

I currently have a problem loading msh files using Gmsh2D. It looks like something is trying to access temp files created by Gmsh2D during loading - maybe it's an internal threading issue? Has anyone encountered this before or have any insights?

Here is a minimalistic example:

import gmsh
import pygmsh
from fipy import Gmsh2D, serialComm

with pygmsh.geo.Geometry() as geom:

    geom.add_circle(x0=[0, 0], radius=1, mesh_size=0.1)
    m = geom.generate_mesh(dim=2)
    gmsh.write('test_mesh.msh')
    gmsh.clear()

mesh = Gmsh2D('test_mesh.msh', communicator=serialComm)

Which results in the error:

Traceback (most recent call last):
  File "C:\Users\batemanc\AppData\Local\Continuum\Anaconda3\envs\FiPySandbox-py3.7\lib\site-packages\fipy\meshes\gmshMesh.py", line 804, in read
    facesData) = self._parseElementFile()
  File "C:\Users\batemanc\AppData\Local\Continuum\Anaconda3\envs\FiPySandbox-py3.7\lib\site-packages\fipy\meshes\gmshMesh.py", line 1134, in _parseElementFile
    currLineInts=currLineInts)
  File "C:\Users\batemanc\AppData\Local\Continuum\Anaconda3\envs\FiPySandbox-py3.7\lib\site-packages\fipy\meshes\gmshMesh.py", line 1063, in _parseTags
    numTags = currLineInts[2]
IndexError: list index out of range

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "C:/Users/batemanc/Documents/source/FiPySandbox/Poisson/gmsh_test.py", line 13, in <module>
    mesh = Gmsh2D('test_mesh.msh', communicator=serialComm)
  File "C:\Users\batemanc\AppData\Local\Continuum\Anaconda3\envs\FiPySandbox-py3.7\lib\site-packages\fipy\meshes\gmshMesh.py", line 1645, in __init__
    self._orderedCellVertexIDs_data) = self.mshFile.read()
  File "C:\Users\batemanc\AppData\Local\Continuum\Anaconda3\envs\FiPySandbox-py3.7\lib\site-packages\fipy\meshes\gmshMesh.py", line 863, in read
    os.unlink(self.elemsPath)
PermissionError: [WinError 32] The process cannot access the file because it is being used by another process: 'C:\\Users\\batemanc\\AppData\\Local\\Temp\\tmp87z94ldrElements'

Process finished with exit code 1

None of the .msh files are in protected folders, and I get the same error when loading/creating the .msh file in separate sessions.

Environment

OS: Windows 10 IDE: Pycharm Env Manager: Anaconda

python==3.7.12 (also tried 3.10) fipy==3.4.3 gmsh==4.10.5 pygmsh==7.1.17

Note: I also had to have a copy the gmsh.exe located in the same directory as the conda environment's python.exe to get around the gmsh version bug.

shinstra commented 2 years ago

A workaround I found for this in the interim is to use the *.geo_unrolled format instead of *.msh. Works for the above environment settings (however didn't work when using gmsh==3.0.6).

shinstra commented 2 years ago

A possibly more reliable workaround is to work with the fipy.meshes.mesh.Mesh class directly - which bypasses the GmshX class's need to load the mesh from file (which I have still found problematic in some cases). As I haven't worked through the 2D implemention of how to do this, here is an example for 3D volumetric meshes:

import itertools
import numpy as np

from meshio import Mesh as MeshioMesh
from fipy.meshes.mesh import Mesh as FipyMesh

def meshio_volume_to_fipy(meshio_mesh: MeshioMesh) -> FipyMesh:
    """
    Transfer volumetric meshio mesh (i.e. generated by gmsh) data to
    a fipy mesh without needing to use an intermediate save file.
    """

    # Extract cell_vertex_ids
    cell_vertex_ids = meshio_mesh.cells_dict['tetra']
    vertex_coords = meshio_mesh.points  # Note, assuming all are used across all cells. Otherwise, would need to remap cell_vertex_ids.

    # Calculate 'face_vertex_ids'
    cell_faces = [[list(vv) for vv in list(itertools.combinations(v, 3))] for v in cell_vertex_ids] # List of face vertex-tuples for each cell
    _ = [[vv.sort() for vv in v] for v in cell_faces]  # sort vertex-tuples so shared faces can be more easily identified in cell_face_ids calculation
    cell_faces = [[tuple(vv) for vv in v] for v in cell_faces]  # convert vertex-tuples from list back to tuple.
    face_vertex_ids = list(set(itertools.chain(*cell_faces))) # Create list of unique face vertex-tuples

    # Calculate 'cell_face_ids'
    cell_face_ids = [[face_vertex_ids.index(vv) for vv in v] for v in cell_faces]

    # Reshape arrays
    vertex_coords = np.array(vertex_coords, dtype=float).T
    face_vertex_ids = np.array(face_vertex_ids, dtype=int).T
    cell_face_ids = np.array(cell_face_ids, dtype=int).T

    return FipyMesh(vertex_coords, face_vertex_ids, cell_face_ids)
wd15 commented 2 years ago

@shinstra: thanks for all the reporting on your progress. I didn't have a suitable environment at hand to check this quickly, but will look through this eventually. Possibly in a few weeks time.

shinstra commented 2 years ago

Awesome, thanks @wd15. No rush on my end as the workaround I posted above got me past the blocking point for my work. In the interim I'll keep adding details of my progress that may be useful to anyone else encountering this problem.

shinstra commented 2 years ago

Here is an updated version of the above workaround that works for both 2D and 3D volumetric meshes.

import itertools
import numpy as np

from meshio import Mesh as MeshioMesh
from fipy.meshes.mesh import Mesh as FipyMesh3D
from fipy.meshes.mesh2D import Mesh2D as FipyMesh2D

from typing import Union

def meshio_to_fipy(meshio_mesh: MeshioMesh, dim: int=3) -> Union[FipyMesh2D, FipyMesh3D]:
    """
    Transfer volumetric meshio mesh (i.e. generated by gmsh) data to
    a fipy mesh without needing to use an intermediate save file.
    """

    # Collect 2D/3D specific params
    match dim:
        case 2:
            cell = 'triangle'
            FipyMesh = FipyMesh2D
        case 3:
            cell = 'tetra'
            FipyMesh = FipyMesh3D
        case _:
            raise ValueError(f'dim must have value of 2 or 3, got {dim}')

    # Extract cell_vertex_ids
    cell_vertex_ids = meshio_mesh.cells_dict[cell]
    vertex_coords = meshio_mesh.points[:, :dim]

    # Remap ids to excude orphan vertices
    keep_ids = list(set(itertools.chain(*cell_vertex_ids)))
    keep_ids.sort()
    id_map = {v: i for i, v in enumerate(keep_ids)}
    cell_vertex_ids = np.vectorize(lambda x: id_map[x])(cell_vertex_ids)
    vertex_coords = vertex_coords[keep_ids, :]

    # Calculate 'face_vertex_ids'
    cell_faces = [[list(vv) for vv in list(itertools.combinations(v, 3))] for v in cell_vertex_ids] # List of face vertex-tuples for each cell
    _ = [[vv.sort() for vv in v] for v in cell_faces]  # sort vertex-tuples so shared faces can be more easily identified in cell_face_ids calculation
    cell_faces = [[tuple(vv) for vv in v] for v in cell_faces]  # convert vertex-tuples from list back to tuple.
    face_vertex_ids = list(set(itertools.chain(*cell_faces))) # Create list of unique face vertex-tuples

    # Calculate 'cell_face_ids'
    cell_face_ids = [[face_vertex_ids.index(vv) for vv in v] for v in cell_faces]

    # Reshape arrays
    vertex_coords = np.array(vertex_coords, dtype=float).T
    face_vertex_ids = np.array(face_vertex_ids, dtype=int).T
    cell_face_ids = np.array(cell_face_ids, dtype=int).T

    return FipyMesh(vertex_coords, face_vertex_ids, cell_face_ids)