Merge pull request #14 from BAMresearch/7-fix-meshing-with-gmsh

7 fix meshing with gmsh

.github/workflows/wall.yml

@@ -0,0 +1,52 @@
+# This workflow will install the environment, run the example/Wall
+name: wall
+
+on:
+  push:
+    branches-ignore:
+
+  # Allows you to run this workflow manually from the Actions tab
+  workflow_dispatch:
+
+env:
+  CACHE_NUMBER: 1 # increase to reset cache manually
+
+jobs:
+  tests:
+    runs-on: ubuntu-latest
+
+    steps:
+      - name: checkout repo content
+        uses: actions/checkout@v2
+      - name: Setup Mambaforge
+        uses: conda-incubator/setup-miniconda@v2
+        with:
+            miniforge-variant: Mambaforge
+            miniforge-version: latest
+            activate-environment: amworkflow
+            use-mamba: true
+
+      - name: Set cache date
+        run: echo "DATE=$(date +'%Y%m%d')" >> $GITHUB_ENV
+
+      - uses: actions/cache@v2
+        with:
+          path: "/usr/share/miniconda3/envs/amworkflow"
+          key: conda-${{ hashFiles('environment.yml') }}-${{ env.DATE }}-${{ env.CACHE_NUMBER }}
+        id: cache
+
+      - name: Update environment
+        run: mamba env update -n amworkflow -f environment.yml
+        if: steps.cache.outputs.cache-hit != 'true'
+
+      - name: Install package
+        run: git clone https://github.com/tpaviot/pythonocc-utils.git && pip install ./pythonocc-utils
+
+      - name: Install amworkflow
+        run: python -m pip install .
+
+      - name: run dodo_wall
+        shell: bash -l {0} #new shell
+        run: |
+          doit -f examples/Wall/dodo_wall.py
+

amworkflow/meshing/meshing.py

@@ -2,6 +2,8 @@
 import typing
 from pathlib import Path
 
+import meshio
+
 import gmsh
 import multiprocessing
 from OCC.Core.TopoDS import TopoDS_Solid
@@ -69,7 +71,8 @@ def create(
         assert step_file.is_file(), f"Step file {step_file} does not exist."
 
         shape = read_step_file(filename=str(step_file))
-        solid = None  # TODO: convert shape to solid
+        solid = occ_helpers.solid_maker(shape)
+
         assert isinstance(solid, TopoDS_Solid), "Must be TopoDS_Shape object to mesh."
 
         gmsh.initialize()
@@ -86,9 +89,8 @@ def create(
 
         model = gmsh.model()
         threads_count = multiprocessing.cpu_count()
-        gmsh.option.setNumber("General.NumThreads", threads_count)
-        # model.add(model_name) # TODO: required?
-
+        # gmsh.option.setNumber("General.NumThreads", threads_count) # FIX: Conflict with doit. Will looking for solutions.
+        # model.add("model name") # TODO: required? Not necessarily but perhaps for output .msh
         layers = model.occ.importShapesNativePointer(int(geo.this), highestDimOnly=True)
         model.occ.synchronize()
         for layer in layers:
@@ -101,18 +103,46 @@ def create(
         phy_gp = model.getPhysicalGroups()
         model_name = model.get_current()
 
-        # save
-        msh, cell_markers, facet_markers = gmshio.model_to_mesh(model, MPI.COMM_SELF, 0)
-        msh.name = model_name
-        cell_markers.name = f"{msh.name}_cells"
-        facet_markers.name = f"{msh.name}_facets"
-        with XDMFFile(msh.comm, out_xdmf, "w") as file:
-            file.write_mesh(msh)
-            file.write_meshtags(cell_markers)
-            msh.topology.create_connectivity(msh.topology.dim - 1, msh.topology.dim)
-            file.write_meshtags(facet_markers)
+        # # save
+        # msh, cell_markers, facet_markers = gmshio.model_to_mesh(model, MPI.COMM_SELF, 0)
+        # msh.name = model_name
+        # cell_markers.name = f"{msh.name}_cells"
+        # facet_markers.name = f"{msh.name}_facets"
+        # with XDMFFile(msh.comm, out_xdmf, "w") as file:
+        #     file.write_mesh(msh)
+        #     file.write_meshtags(cell_markers)
+        #     msh.topology.create_connectivity(msh.topology.dim - 1, msh.topology.dim)
+        #     file.write_meshtags(facet_markers)
+
+        #workaround for since gmshio.model_to_mesh is not working
+        out_msh = out_xdmf.with_suffix(".msh")
+        gmsh.write(str(out_msh))
+        msh = meshio.read(out_msh)
+        mesh = self.create_mesh(msh, "tetra")
+        meshio.write(out_xdmf, mesh)
 
         if out_vtk:
-            gmsh.write(out_vtk)
+            gmsh.write(str(out_vtk))
+
 
         return
+
+    def create_mesh(self, mesh, cell_type: str, prune_z: bool = False) -> meshio.Mesh:
+        """Convert meshio mesh to fenics compatible mesh.
+            based on https://jsdokken.com/dolfinx-tutorial/chapter3/subdomains.html?highlight=read_mesh
+
+        Args:
+            mesh: Mesh read by meshio from msh file (meshio.read(file_msh)).
+            cell_type: Type of cell to be meshed (e.g. 'tetra','triangle' ...).
+            prune_z: True for 2D meshes - removes z coordinate.
+
+        Returns:
+            out_mesh: Mesh which can be saved as xdmf file. (meshio.write(file_xdmf, out_mesh))
+        """
+        cells = mesh.get_cells_type(cell_type)
+        cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
+        points = mesh.points[:, :2] if prune_z else mesh.points
+        out_mesh = meshio.Mesh(
+            points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]}
+        )
+        return out_mesh

environment.yml

@@ -4,12 +4,13 @@ channels:
   - defaults
 dependencies:
   - pythonocc-core=7.7.0
-  - python=3.10
+  - python >= 3.10
   - doit
   - pandas
   - sqlite
   - numpy
-  - fenics-dolfinx
+  - fenics-dolfinx=0.6.0
+  - meshio
   - mpich
   - pyvista
   - pip
@@ -19,4 +20,4 @@ dependencies:
     - pyqt5
     - sphinx
     - sqlalchemy
-    - ruamel-yaml
+    - ruamel-yaml