# coding: utf8 """ Created on 17/11/2018 @author: baptiste """ import logging import matplotlib.pyplot as plt import numpy as np import ho_homog.toolbox_FEniCS as fetools import dolfin as fe import ho_homog.materials as mat from subprocess import run from pathlib import Path from ho_homog.toolbox_gmsh import msh_conversion plt.ioff() logger = logging.getLogger(__name__) # http://sametmax.com/ecrire-des-logs-en-python/ class FenicsPart(object): """ Contrat : Créer un couple maillage + matériaux pour des géométries 2D, planes. """ def __init__( self, mesh, materials, subdomains, global_dimensions=None, facet_regions=None ): """ Parameters ---------- subdomains : dolfin.MeshFunction Indicates the subdomains that have been defined in the mesh. """ self.mesh = mesh self.dim = mesh.topology().dim() # dimension d'espace de depart self.materials = materials self.subdomains = subdomains self.global_dimensions = global_dimensions self.facet_regions = facet_regions if isinstance(materials, mat.Material): self.elasticity_tensor = fe.as_matrix(materials.get_C()) elif isinstance(materials, dict): self.elasticity_tensor = mat.mat_per_subdomains( self.subdomains, self.materials, self.dim ) else: raise TypeError( "materials parameter must be an instance of Material or a dictionnary that contains Material instances." ) def mat_area(self): try: return self._mat_area except AttributeError: self._mat_area = fe.assemble(fe.Constant(1) * fe.dx(self.mesh)) return self._mat_area def global_area(self): if not self.global_dimensions is None: return np.linalg.det(self.global_dimensions) else: raise AttributeError( f"global_size information is lacking for FenicsPart {self}." ) @staticmethod def file_2_FenicsPart( mesh_path, materials, global_dimensions=None, subdomains_import=False, plots=True, explicit_subdo_val=0, ): """Generate an instance of FenicsPart from a .xml or .msh file that contains the mesh. Parameters ---------- mesh_path : string or Path Relative or absolute path to the mesh file (format Dolfin XML, XDMF or MSH version 2) global_dimensions : 2D-array shape : 2×2 if 2D problem dimensions of the RVE materials : dict or Material [description] #TODO : à compléter subdomains_import : bool Import subdomains data ? plots : bool, optional If True (default) the physical regions and the facet regions are plotted at the end of the import. Returns ------- FenicsPart instance """ if not isinstance(mesh_path, Path): mesh_path = Path(mesh_path) name = mesh_path.stem suffix = mesh_path.suffix if suffix not in fetools.SUPPORTED_MESH_SUFFIX: mesh_file_paths = msh_conversion( mesh_path, format_=".xml", subdomains=subdomains_import ) try: mesh_path = mesh_file_paths[0] except IndexError as error: mesh_path = mesh_file_paths logger.warning(error) suffix = mesh_path.suffix # Each supported mesh format -> one if structure subdomains, facets = None, None if suffix == ".xml": mesh = fe.Mesh(str(mesh_path)) if subdomains_import: subdo_path = mesh_path.with_name(f"{name}_physical_region.xml") facet_path = mesh_path.with_name(f"{name}_facet_region.xml") if subdo_path.exists(): subdomains = fe.MeshFunction("size_t", mesh, str(subdo_path)) subdo_val = fetools.get_MeshFunction_val(subdomains) logger.info(f"{subdo_val[0]} physical regions imported.") logger.info(f"The values of their tags are : {subdo_val[1]}") else: logger.info( f"For mesh file {mesh_path.name}, _physical_region.xml file is missing." ) if facet_path.exists(): facets = fe.MeshFunction("size_t", mesh, str(facet_path)) facets_val = fetools.get_MeshFunction_val(facets) logger.info(f"{facets_val[0]} facet regions imported.") logger.info(f"The values of their tags are : {facets_val[1]}") else: logger.info( f"For mesh file {mesh_path.name}, _facet_region.xml file is missing." ) if suffix == ".xdmf": if subdomains_import: mesh, subdomains, facets = fetools.xdmf_mesh(mesh_path, True) logger.info( f"Import of a mesh from {mesh_path} file, with subdomains data" ) else: mesh = fetools.xdmf_mesh(mesh_path, False) logger.info( f"Import of a mesh from {mesh_path} file, without subdomains data" ) if plots: plt.figure() fe.plot(mesh, title=f"Mesh of {name}") if subdomains_import: plt.figure() subdo_plt = fe.plot( subdomains, mesh, title=f"Subdomains imported for {name}" ) plt.colorbar(subdo_plt) plt.figure() facets_val_range = max(facets_val[1]) - min(facets_val[1]) cmap = plt.cm.get_cmap("viridis", facets_val_range) facets_plt = fetools.facet_plot2d(facets, mesh, cmap=cmap) plt.colorbar(facets_plt[0]) plt.show() logger.info(f"Import of the mesh : DONE") return FenicsPart(mesh, materials, subdomains, global_dimensions, facets) class Fenics2DRVE(FenicsPart): """ Contrat : Créer un couple maillage + matériaux pour des RVE 2D, plans, comportant au plus 2 matériaux constitutifs et pouvant contenir plusieurs cellules. """ def __init__( self, mesh, generating_vectors, material_dict, subdomains, facet_regions ): """ Parameters ---------- subdomains : dolfin.MeshFunction Indicates the subdomains that have been defined in the mesh. #! L'ordre des facet function à probablement de l'importance pour la suite des opérations """ self.mesh = mesh self.gen_vect = generating_vectors self.rve_area = np.linalg.det(self.gen_vect) self.mat_area = fe.assemble(fe.Constant(1) * fe.dx(mesh)) self.mesh_dim = mesh.topology().dim() # dimension d'espace de depart self.materials = material_dict self.subdomains = subdomains self.facet_regions = facet_regions if isinstance(self.materials, mat.Material): self.C_per = fe.as_matrix(self.materials.get_C()) else: self.C_per = mat.mat_per_subdomains( self.subdomains, self.materials, self.mesh_dim ) def epsilon(self, u): return mat.epsilon(u) def sigma(self, eps): return mat.sigma(self.C_per, eps) def StrainCrossEnergy(self, sig, eps): return mat.strain_cross_energy(sig, eps, self.mesh, self.rve_area) @staticmethod def gmsh_2_Fenics_2DRVE(gmsh_2D_RVE, material_dict, plots=True): """ Generate an instance of Fenics2DRVE from a instance of the Gmsh2DRVE class. """ xml_path = gmsh_2D_RVE.mesh_abs_path.with_suffix(".xml") cmd = ( "dolfin-convert " + gmsh_2D_RVE.mesh_abs_path.as_posix() + " " + xml_path.as_posix() ) run(cmd, shell=True, check=True) mesh = fe.Mesh(xml_path.as_posix()) name = xml_path.stem subdomains = fe.MeshFunction( "size_t", mesh, xml_path.with_name(name + "_physical_region.xml").as_posix() ) facets = fe.MeshFunction( "size_t", mesh, xml_path.with_name(name + "_facet_region.xml").as_posix() ) subdo_val = fetools.get_MeshFunction_val(subdomains) facets_val = fetools.get_MeshFunction_val(facets) logger.info( f"{subdo_val[0]} physical regions imported. The values of their tags are : {subdo_val[1]}" ) logger.info( f"{facets_val[0]} facet regions imported. The values of their tags are : {facets_val[1]}" ) if plots: plt.figure() subdo_plt = fe.plot(subdomains) plt.colorbar(subdo_plt) plt.figure() cmap = plt.cm.get_cmap("viridis", max(facets_val[1]) - min(facets_val[1])) facets_plt = fetools.facet_plot2d(facets, mesh, cmap=cmap) plt.colorbar(facets_plt[0]) # clrbar.set_ticks(facets_val) plt.draw() logger.info(f"Import of the mesh : DONE") generating_vectors = gmsh_2D_RVE.gen_vect return Fenics2DRVE(mesh, generating_vectors, material_dict, subdomains, facets) @staticmethod def file_2_Fenics_2DRVE(mesh_path, generating_vectors, material_dict, plots=True): """Generate an instance of Fenics2DRVE from a .xml, .msh or .xdmf file that contains the mesh. Parameters ---------- mesh_path : string or Path Relative or absolute path to the mesh file (format Dolfin XML, XDMF or MSH version 2) generating_vectors : 2D-array dimensions of the RVE material_dict : dictionnary [description] #TODO : à compléter plots : bool, optional If True (default) the physical regions and the facet regions are plotted at the end of the import. Returns ------- Fenics2DRVE instance """ if not isinstance(mesh_path, Path): mesh_path = Path(mesh_path) name = mesh_path.stem if mesh_path.suffix == ".xml": mesh = fe.Mesh(str(mesh_path)) elif mesh_path.suffix == ".xdmf": mesh = fe.Mesh() with fe.XDMFFile(str(mesh_path)) as file_in: file_in.read(mesh) else: cmd = f"dolfin-convert {mesh_path} {mesh_path.with_suffix('.xml')}" run(cmd, shell=True, check=True) mesh_path = mesh_path.with_suffix(".xml") mesh = fe.Mesh(str(mesh_path)) subdo_path = mesh_path.with_name(name + "_physical_region.xml") facet_path = mesh_path.with_name(name + "_facet_region.xml") if subdo_path.exists(): subdomains = fe.MeshFunction("size_t", mesh, subdo_path.as_posix()) subdo_val = fetools.get_MeshFunction_val(subdomains) logger.info( f"{subdo_val[0]} physical regions imported. Tags : {subdo_val[1]}" ) else: logger.info( f"For mesh file {mesh_path.name}, _physical_region.xml file is missing." ) subdomains = None if facet_path.exists(): facets = fe.MeshFunction("size_t", mesh, facet_path.as_posix()) facets_val = fetools.get_MeshFunction_val(facets) logger.info( f"{facets_val[0]} facet regions imported. Tags : {facets_val[1]}" ) else: logger.info( f"For mesh file {mesh_path.name}, _facet_region.xml file is missing." ) facets = None if plots: plt.figure() fe.plot(mesh) if subdomains is not None: plt.figure() subdo_plt = fe.plot(subdomains) plt.colorbar(subdo_plt) if facets is not None: plt.figure() cmap = plt.cm.get_cmap( "viridis", max(facets_val[1]) - min(facets_val[1]) ) facets_plt = fetools.facet_plot2d(facets, mesh, cmap=cmap) plt.colorbar(facets_plt[0]) plt.draw() logger.info(f"Import of the mesh : DONE") return Fenics2DRVE(mesh, generating_vectors, material_dict, subdomains, facets) if __name__ == "__main__": pass # geo.init_geo_tools() # a = 1 # b, k = a, a/3 # panto_test = Gmsh2DRVE.pantograph(a, b, k, 0.1, nb_cells=(2, 3), soft_mat=False, name='panto_test') # panto_test.main_mesh_refinement((0.1,0.5),(0.03,0.3),False) # panto_test.mesh_generate() # os.system(f"gmsh {panto_test.name}.msh &") # a = 1 # b, k = a, a/3 # panto_test_offset = Gmsh2DRVE.pantograph(a, b, k, 0.1, nb_cells=(2,3), offset=(0.25,0.25), soft_mat=False, name='panto_test_offset') # panto_test_offset.main_mesh_refinement((0.1,0.5),(0.03,0.3),False) # panto_test_offset.mesh_generate() # os.system(f"gmsh {panto_test_offset.name}.msh &") # L, t = 1, 0.05 # a = L-3*t # aux_sqr_test = Gmsh2DRVE.auxetic_square(a, L, t, nb_cells=(4,3), soft_mat=False, name='aux_square_test') # os.system(f"gmsh {aux_sqr_test.name}.brep &") # aux_sqr_test.main_mesh_refinement((0.1,0.3), (0.01,0.05), False) # aux_sqr_test.mesh_generate() # os.system(f"gmsh {aux_sqr_test.name}.msh &") # a = 1 # b = a # w = a/50 # r = 4*w # beam_panto_test = Gmsh2DRVE.beam_pantograph(a, b, w, r, nb_cells=(1, 1), offset=(0., 0.), soft_mat=False, name='beam_panto_test') # os.system(f"gmsh {beam_panto_test.name}.brep &") # beam_panto_test.main_mesh_refinement((5*w, a/2),(w/5, w), True) # beam_panto_test.mesh_generate() # os.system(f"gmsh {beam_panto_test.name}.msh &") # gmsh.option.setNumber('Mesh.SurfaceFaces',1) #Display faces of surface mesh? # gmsh.fltk.run() # msh.set_background_mesh(field) # gmsh.option.setNumber('Mesh.CharacteristicLengthExtendFromBoundary',0) # geo.PhysicalGroup.set_group_mesh(True) # model.mesh.generate(1) # model.mesh.generate(2) # gmsh.write(f"{self.name}.msh") # os.system(f"gmsh {self.name}.msh &")