Commit 528c5c3a authored by Baptiste Durand's avatar Baptiste Durand

Demo script : Mesh import from xdmf

Including subdomains and with plots.
+ minor corrections linked to the msh_conversion function.
parent 0e4920a9
# 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()
...@@ -86,7 +86,7 @@ gmsh.logger.stop() ...@@ -86,7 +86,7 @@ gmsh.logger.stop()
gmsh.model.mesh.renumberNodes() gmsh.model.mesh.renumberNodes()
gmsh.model.mesh.renumberElements() gmsh.model.mesh.renumberElements()
gmsh.write(str(rve_geo.mesh_abs_path)) gmsh.write(str(rve_geo.mesh_abs_path))
rve_path, *_ = 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 # * Step 3 : Build the mesh of the part from the mesh of the RVE
...@@ -94,7 +94,7 @@ gmsh.logger.start() ...@@ -94,7 +94,7 @@ gmsh.logger.start()
part_geo = mesh_generate_2D.Gmsh2DPartFromRVE(rve_geo, (75, 1)) part_geo = mesh_generate_2D.Gmsh2DPartFromRVE(rve_geo, (75, 1))
process_gmsh_log(gmsh.logger.get()) process_gmsh_log(gmsh.logger.get())
gmsh.logger.stop() gmsh.logger.stop()
part_path, *_ = msh_conversion(part_geo.mesh_abs_path, ".xdmf") part_path, = msh_conversion(part_geo.mesh_abs_path, ".xdmf")
# * Step 4 : Defining the material properties # * Step 4 : Defining the material properties
E, nu = 1.0, 0.3 E, nu = 1.0, 0.3
......
...@@ -292,7 +292,7 @@ def xdmf_mesh(mesh_file, import_subdomains=False, facet_file="", physical_file=" ...@@ -292,7 +292,7 @@ def xdmf_mesh(mesh_file, import_subdomains=False, facet_file="", physical_file="
facet_vc = fe.MeshValueCollection("size_t", mesh, dim - 1) facet_vc = fe.MeshValueCollection("size_t", mesh, dim - 1)
with fe.XDMFFile(str(facet_file)) as f_in: with fe.XDMFFile(str(facet_file)) as f_in:
f_in.read(facet_vc, "facet_data") f_in.read(facet_vc, "facet_data")
facet_function = fe.cpp.mesh.MeshFunctionSizet(mesh, facet_vc) facet_regions = fe.cpp.mesh.MeshFunctionSizet(mesh, facet_vc)
if not physical_file.exists(): if not physical_file.exists():
subdomains = None subdomains = None
else: else:
......
...@@ -44,7 +44,7 @@ def process_gmsh_log(gmsh_log: list, detect_error=True): ...@@ -44,7 +44,7 @@ def process_gmsh_log(gmsh_log: list, detect_error=True):
def msh_conversion( def msh_conversion(
mesh, format_: str = ".xdmf", output_dir=None, subdomain: bool = False, dim: int = 2 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. Convert a ".msh" mesh generated with Gmsh to a format suitable for FEniCS.
...@@ -62,7 +62,7 @@ def msh_conversion( ...@@ -62,7 +62,7 @@ def msh_conversion(
Path of the directory where the converted mesh file must be written. Path of the directory where the converted mesh file must be written.
If None, the converted file is written in the same directory If None, the converted file is written in the same directory
as the input file. (default: None) as the input file. (default: None)
subdomain : bool, optional subdomains : bool, optional
If True, extra files are created to store information about subdomains. If True, extra files are created to store information about subdomains.
(default: False) (default: False)
dim: int, optional dim: int, optional
...@@ -118,7 +118,7 @@ def msh_conversion( ...@@ -118,7 +118,7 @@ def msh_conversion(
else: else:
ValueError ValueError
meshio.write(str(mesh_path), geo_only) meshio.write(str(mesh_path), geo_only)
if subdomain: if subdomains:
cell_function = meshio.Mesh( cell_function = meshio.Mesh(
points=mesh.points, points=mesh.points,
cells={cell: mesh.cells[cell]}, cells={cell: mesh.cells[cell]},
...@@ -131,11 +131,11 @@ def msh_conversion( ...@@ -131,11 +131,11 @@ def msh_conversion(
cell_data={face: {"facet_data": mesh.cell_data[face]["gmsh:physical"]}}, cell_data={face: {"facet_data": mesh.cell_data[face]["gmsh:physical"]}},
) )
meshio.write(str(facet_region), facet_function) meshio.write(str(facet_region), facet_function)
if subdomain: if subdomains:
return ( return (
mesh_path, mesh_path,
physical_region if physical_region.exist() else None, physical_region if physical_region.exists() else None,
facet_region if facet_region.exist() else None, facet_region if facet_region.exists() else None,
) )
else: else:
return mesh_path return mesh_path
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment