Merge branch 'xdmf_management' into 'master'

Xdmf management

See merge request !3

demo/demo_import_mesh_xdmf.py

+# coding: utf-8
+"""
+Created on 01/07/2019
+@author: baptiste
+
+1- Conversion of a mesh file MSH format (version 4) -> XDMF format,
+2- Import of the mesh from the XDMF file,
+3- Show facet regions and subdomains with matplotlib.
+"""
+
+import logging
+
+import gmsh
+import matplotlib.pyplot as plt
+import dolfin as fe
+from ho_homog import geometry, mesh_generate_2D, toolbox_FEniCS, toolbox_gmsh
+
+# * Logging
+logger = logging.getLogger(__name__)  # http://sametmax.com/ecrire-des-logs-en-python/
+logger.setLevel(logging.INFO)
+
+formatter = logging.Formatter("%(asctime)s :: %(levelname)s :: %(message)s", "%H:%M")
+stream_handler = logging.StreamHandler()
+stream_handler.setLevel(logging.INFO)
+stream_handler.setFormatter(formatter)
+logger.addHandler(stream_handler)
+
+geometry.init_geo_tools()
+geometry.set_gmsh_option("Mesh.MshFileVersion", 4.1)
+
+logger.info("Generating the geometry model")
+a = 1
+b, k = a, a / 3
+r = 0.01 * a
+rve_geo = mesh_generate_2D.Gmsh2DRVE.pantograph(a, b, k, r, soft_mat=True, name="panto")
+
+logger.info("Generating the mesh")
+lc_ratio = 1 / 3.0
+d_min_max = (2 * r, a)
+lc_min_max = (lc_ratio * r, lc_ratio * a)
+rve_geo.main_mesh_refinement(d_min_max, lc_min_max, False)
+rve_geo.soft_mesh_refinement(d_min_max, lc_min_max, False)
+rve_geo.mesh_generate()
+
+logger.info("Saving mesh with MSH 4 format")
+gmsh.model.mesh.renumberNodes()
+gmsh.model.mesh.renumberElements()
+gmsh.write(str(rve_geo.mesh_abs_path))
+
+logger.info("Mesh conversion MSH -> XDMF")
+mesh_path, *_ = toolbox_gmsh.msh_conversion(
+    rve_geo.mesh_abs_path, ".xdmf", subdomains=True
+)
+
+logger.info("Import of mesh and partitions as MeshFunction instances")
+mesh, subdomains, facet_regions = toolbox_FEniCS.xdmf_mesh(mesh_path, True)
+
+
+plt.figure()
+fe.plot(mesh, title="Mesh only")
+plt.figure()
+subdo_plt = fe.plot(subdomains, title="Subdomains")
+plt.colorbar(subdo_plt)
+
+nb_val, facets_val = toolbox_FEniCS.get_MeshFunction_val(facet_regions)
+facets_val = facets_val[facets_val != 18446744073709551615]
+facets_val = facets_val[facets_val != 0]
+print(facets_val)
+delta = max(facets_val) - max(facets_val)
+cmap = plt.cm.get_cmap("viridis", max(facets_val) - min(facets_val)+1)
+fig, ax = plt.subplots()
+facets_plt = toolbox_FEniCS.facet_plot2d(
+    facet_regions, mesh, cmap=cmap, exclude_val=(0, 18446744073709551615)
+)
+ax.set_title("Facet regions")
+cbar = fig.colorbar(facets_plt[0], ticks=list(facets_val))
+cbar.ax.set_xticklabels(list(facets_val))
+
+plt.show()

demo/demo_pantograph_compare_same_mesh.py

@@ -29,7 +29,7 @@ from ho_homog import (
     periodicity,
 )
 from ho_homog.toolbox_FEniCS import function_errornorm
-from ho_homog.toolbox_gmsh import process_gmsh_log
+from ho_homog.toolbox_gmsh import process_gmsh_log, msh_conversion
 
 logger = logging.getLogger("demo_full_compare")
 logger_root = logging.getLogger()
@@ -86,7 +86,7 @@ gmsh.logger.stop()
 gmsh.model.mesh.renumberNodes()
 gmsh.model.mesh.renumberElements()
 gmsh.write(str(rve_geo.mesh_abs_path))
-rve_path, *_ = mesh_tools.msh_conversion(rve_geo.mesh_abs_path, ".xdmf")
+rve_path, = msh_conversion(rve_geo.mesh_abs_path, ".xdmf")
 
 
 # * Step 3 : Build the mesh of the part from the mesh of the RVE
@@ -94,7 +94,7 @@ gmsh.logger.start()
 part_geo = mesh_generate_2D.Gmsh2DPartFromRVE(rve_geo, (75, 1))
 process_gmsh_log(gmsh.logger.get())
 gmsh.logger.stop()
-part_path, *_ = mesh_tools.msh_conversion(part_geo.mesh_abs_path, ".xdmf")
+part_path, = msh_conversion(part_geo.mesh_abs_path, ".xdmf")
 
 # * Step 4 : Defining the material properties
 E, nu = 1.0, 0.3

ho_homog/mesh_tools.py

@@ -17,9 +17,6 @@ import numpy as np
 
 import ho_homog.geometry as geo
 import gmsh
-from pathlib import Path
-import meshio
-from subprocess import run
 
 # nice shortcuts
 model = gmsh.model
@@ -440,94 +437,4 @@ def order_curves(curves, dir_v, orientation=False):
     return None
 
 
-def msh_conversion(
-    mesh, format_: str = ".xdmf", output_dir=None, subdomain: bool = False, dim: int = 2
-):
-    """
-    Convert a ".msh" mesh generated with Gmsh to a format suitable for FEniCS.
-
-    Parameters
-    ----------
-    mesh : Path or str
-        Path that points to the existing mesh file.
-    format : str
-        Suffix desired for the mesh file. (default: ".xdmf")
-        Supported suffixes :
-            - ".xdmf"
-            - ".xml" (DOLFIN xml format)
-    output_dir : Path, optional
-        Path of the directory where the converted mesh file must be written.
-        If None, the converted file is written in the same directory
-        as the input file. (default: None)
-    subdomain : bool, optional
-        If True, extra files are created to store information about subdomains.
-        (default: False)
-    dim: int, optional
-        Geometrical dimension of the mesh (2D or 3D). (default: 2)
 
-    Returns
-    -------
-    tuple
-        First element : Path to the main mesh file
-        Then paths to the extra files if subdomain conversion is requested.
-
-    Warning
-    -------
-    A specific version of the MSH format should be use in accordance with the
-    desired output format :
-        - ".xml" output -> MSH file format version 2;
-        - ".xdmf" output -> MSH file format version 4.
-    """
-
-    input_path = Path(mesh)
-    name = input_path.stem
-    if format_ not in (".xml", ".xdmf"):
-        raise TypeError
-    mesh_path = input_path.with_suffix(format_)
-    if output_dir:
-        mesh_path = output_dir.joinpath(mesh_path.name)
-    physical_region = mesh_path.with_name(name + "_physical_region" + format_)
-    facet_region = mesh_path.with_name(name + "_facet_region" + format_)
-    if physical_region.exists():
-        physical_region.unlink()
-    if facet_region.exists():
-        facet_region.unlink()
-    if format_ == ".xml":
-        cmd = f"dolfin-convert {input_path} {mesh_path}"
-        run(cmd, shell=True, check=True)
-    elif format_ == ".xdmf":
-        mesh = meshio.read(str(input_path))
-
-        if dim == 2:
-            mesh.points = mesh.points[:, :2]
-            geo_only = meshio.Mesh(
-                points=mesh.points, cells={"triangle": mesh.cells["triangle"]}
-            )
-            cell = "triangle"
-            face = "line"
-        elif dim == 3:
-            raise NotImplementedError("3D meshes are not supported yet.")
-        else:
-            ValueError
-        meshio.write(str(mesh_path), geo_only)
-        if subdomain:
-            cell_function = meshio.Mesh(
-                points=mesh.points,
-                cells={cell: mesh.cells[cell]},
-                cell_data={cell: {"cell_data": mesh.cell_data[cell]["gmsh:physical"]}},
-            )
-            meshio.write(str(physical_region), cell_function)
-            facet_function = meshio.Mesh(
-                points=mesh.points,
-                cells={face: mesh.cells[face]},
-                cell_data={face: {"facet_data": mesh.cell_data[face]["gmsh:physical"]}},
-            )
-            meshio.write(str(facet_region), facet_function)
-    if subdomain:
-        return (
-            mesh_path,
-            physical_region if physical_region.exist() else None,
-            facet_region if facet_region.exist() else None,
-        )
-    else:
-        return (mesh_path,)

ho_homog/toolbox_FEniCS.py

@@ -11,6 +11,7 @@ import logging
 import dolfin as fe
 import matplotlib.pyplot as plt
 import numpy as np
+from pathlib import Path, PurePath
 
 logger = logging.getLogger(__name__)
 
@@ -221,3 +222,83 @@ def _wrap_in_list(obj, name, types=type):
                 "expected a (list of) %s as '%s' argument" % (str(types), name)
             )
     return lst
+
+
+def xdmf_mesh(mesh_file, import_subdomains=False, facet_file="", physical_file=""):
+    """Create a FeniCS mesh from a mesh file with xdmf format.
+
+    Parameters
+    ----------
+    mesh_file : str or Path
+        Path of the file mesh (xdmf file)
+    import_subdomains : bool, optional
+        True if information about subdomains have to be imported.
+        The paths of the auxiliary files that contains information about subdomains
+        can be indicated with facet_file and physical_file.
+        The paths used by default are :
+            - "<mesh path>_facet_region.xdmf" and
+            - "<mesh path>_physical_region.xdmf" (for subdomains)
+    facet_file : str or Path, optional
+        Path to the mesh auxiliary file that contains subdomains data.
+        Defaults to "" i.e. the default path will be used.
+    physical_file : str or Path, optional
+        Path to the mesh auxiliary file that contains facet regions data.
+        Defaults to "" i.e. the default path will be used.
+
+    Returns
+    -------
+    Mesh / tuple
+        If subdomains are not requested :
+            - The Mesh instance
+        If subdomains are requested :
+            - The Mesh instance,
+            - The MeshFunction for subdomains if it exists else None,
+            - The MeshFunction for facets if it exists else None,
+
+    Source
+    ------
+    Gist meshtagging_mvc.py, June 2018, Michal Habera
+    https://gist.github.com/michalhabera/bbe8a17f788192e53fd758a67cbf3bed
+    """
+    if not isinstance(mesh_file, PurePath):
+        mesh_file = Path(mesh_file)
+    if not mesh_file.suffix == ".xdmf":
+        raise TypeError("Wrong suffix for the path to the mesh.")
+    mesh = fe.Mesh()
+    with fe.XDMFFile(str(mesh_file)) as f_in:
+        f_in.read(mesh)
+
+    if not import_subdomains:
+        return mesh
+
+    dim = mesh.geometric_dimension()
+
+    if not facet_file:
+        facet_file = mesh_file.with_name(f"{mesh_file.stem}_facet_region.xdmf")
+    else:
+        facet_file = Path(facet_file)
+        if not facet_file.suffix == ".xdmf":
+            raise TypeError("Wrong suffix for the path to facet regions.")
+    if not physical_file:
+        physical_file = mesh_file.with_name(f"{mesh_file.stem}_physical_region.xdmf")
+    else:
+        physical_file = Path(physical_file)
+        if not physical_file.suffix == ".xdmf":
+            raise TypeError("Wrong suffix for the path to subdomains.")
+
+    if not facet_file.exists():
+        facet_regions = None
+    else:
+        facet_vc = fe.MeshValueCollection("size_t", mesh, dim - 1)
+        with fe.XDMFFile(str(facet_file)) as f_in:
+            f_in.read(facet_vc, "facet_data")
+        facet_regions = fe.cpp.mesh.MeshFunctionSizet(mesh, facet_vc)
+    if not physical_file.exists():
+        subdomains = None
+    else:
+        cell_vc = fe.MeshValueCollection("size_t", mesh, dim - 1)
+        with fe.XDMFFile(str(physical_file)) as f_in:
+            f_in.read(cell_vc, "cell_data")
+        subdomains = fe.cpp.mesh.MeshFunctionSizet(mesh, cell_vc)
+
+    return mesh, subdomains, facet_regions

ho_homog/toolbox_gmsh.py

@@ -8,6 +8,9 @@ Collection of tools designed to help users working with gmsh python API.
 """
 
 import logging
+from pathlib import Path
+import meshio
+from subprocess import run
 
 gmsh_logger = logging.getLogger("ho_homog.gmsh")
 
@@ -38,3 +41,101 @@ def process_gmsh_log(gmsh_log: list, detect_error=True):
     gmsh_logger.info("**********")
     if detect_error and err_msg:
         raise AssertionError("Gmsh logging messages signal errors.")
+
+
+def msh_conversion(
+    mesh, format_: str = ".xdmf", output_dir=None, subdomains: bool = False, dim: int = 2
+):
+    """
+    Convert a ".msh" mesh generated with Gmsh to a format suitable for FEniCS.
+
+    Parameters
+    ----------
+    mesh : Path or str
+        Path that points to the existing mesh file.
+    format : str
+        Suffix desired for the mesh file. (default: ".xdmf")
+        Supported suffixes :
+            - ".xdmf"
+            - ".xml" (DOLFIN xml format)
+    output_dir : Path, optional
+        Path of the directory where the converted mesh file must be written.
+        If None, the converted file is written in the same directory
+        as the input file. (default: None)
+    subdomains : bool, optional
+        If True, extra files are created to store information about subdomains.
+        (default: False)
+    dim: int, optional
+        Geometrical dimension of the mesh (2D or 3D). (default: 2)
+
+    Returns
+    -------
+    Path / tuple
+    If subdomain conversion is not requested :
+        - Path to the mesh file
+    If subdomain conversion is requested :
+        - Path to the mesh file,
+        - Path to the extra files for subdomains if it exists else None,
+        - Path to the extra files for facet regions if it exists else None,
+
+    Warning
+    -------
+    A specific version of the MSH format should be use in accordance with the
+    desired output format :
+        - ".xml" output -> MSH file format version 2;
+        - ".xdmf" output -> MSH file format version 4.
+    """
+    input_path = Path(mesh)
+    name = input_path.stem
+
+    if format_ not in (".xml", ".xdmf"):
+        raise TypeError
+    mesh_path = input_path.with_suffix(format_)
+    if output_dir:
+        mesh_path = output_dir.joinpath(mesh_path.name)
+    physical_region = mesh_path.with_name(name + "_physical_region" + format_)
+    facet_region = mesh_path.with_name(name + "_facet_region" + format_)
+    if physical_region.exists():
+        physical_region.unlink()
+    if facet_region.exists():
+        facet_region.unlink()
+
+    if format_ == ".xml":
+        cmd = f"dolfin-convert {input_path} {mesh_path}"
+        run(cmd, shell=True, check=True)
+    elif format_ == ".xdmf":
+        mesh = meshio.read(str(input_path))
+
+        if dim == 2:
+            mesh.points = mesh.points[:, :2]
+            geo_only = meshio.Mesh(
+                points=mesh.points, cells={"triangle": mesh.cells["triangle"]}
+            )
+            cell = "triangle"
+            face = "line"
+        elif dim == 3:
+            raise NotImplementedError("3D meshes are not supported yet.")
+        else:
+            ValueError
+        meshio.write(str(mesh_path), geo_only)
+        if subdomains:
+            cell_function = meshio.Mesh(
+                points=mesh.points,
+                cells={cell: mesh.cells[cell]},
+                cell_data={cell: {"cell_data": mesh.cell_data[cell]["gmsh:physical"]}},
+            )
+            meshio.write(str(physical_region), cell_function)
+            facet_function = meshio.Mesh(
+                points=mesh.points,
+                cells={face: mesh.cells[face]},
+                cell_data={face: {"facet_data": mesh.cell_data[face]["gmsh:physical"]}},
+            )
+            meshio.write(str(facet_region), facet_function)
+    if subdomains:
+        return (
+            mesh_path,
+            physical_region if physical_region.exists() else None,
+            facet_region if facet_region.exists() else None,
+        )
+    else:
+        return mesh_path

test/test_mesh_generate_2D.py

@@ -5,8 +5,8 @@ Created on 12/06/2019
 
 """
 from ho_homog import geometry as geo
-from ho_homog import mesh_tools
 from ho_homog import mesh_generate_2D as mesh2D
+from ho_homog.toolbox_gmsh import msh_conversion
 import gmsh
 
 geo.init_geo_tools()
@@ -30,5 +30,5 @@ def test_rve_2_part():
     gmsh.model.mesh.renumberElements()
     gmsh.write(str(panto_rve.mesh_abs_path))
     panto_part = mesh2D.Gmsh2DPartFromRVE(panto_rve, (2, 3))
-    mesh_tools.msh_conversion(panto_rve.mesh_abs_path, ".xdmf")
-    mesh_tools.msh_conversion(panto_part.mesh_abs_path, ".xdmf")
+    msh_conversion(panto_rve.mesh_abs_path, ".xdmf")
+    msh_conversion(panto_part.mesh_abs_path, ".xdmf")