part.py 14.2 KB
Newer Older
Baptiste Durand's avatar
Baptiste Durand committed
1
# coding: utf8
2 3 4 5 6 7 8 9
"""
Created on 17/11/2018
@author: baptiste

"""

import logging
import matplotlib.pyplot as plt
10
import numpy as np
11
import ho_homog.toolbox_FEniCS as fetools
12
import dolfin as fe
13
import ho_homog.materials as mat
14
from subprocess import run
15
from pathlib import Path
16
from ho_homog.toolbox_gmsh import msh_conversion
17 18 19

plt.ioff()

20
logger = logging.getLogger(__name__)  # http://sametmax.com/ecrire-des-logs-en-python/
21

22

23
class FenicsPart(object):
24 25 26
    """
    Contrat : Créer un couple maillage + matériaux pour des géométries 2D, planes.
    """
27 28 29 30

    def __init__(
        self, mesh, materials, subdomains, global_dimensions=None, facet_regions=None
    ):
31 32 33 34 35 36 37
        """
        Parameters
        ----------
        subdomains : dolfin.MeshFunction
            Indicates the subdomains that have been defined in the mesh.
        """
        self.mesh = mesh
38
        self.dim = mesh.topology().dim()  # dimension d'espace de depart
39
        self.materials = materials
40
        self.subdomains = subdomains
41
        self.global_dimensions = global_dimensions
42
        self.facet_regions = facet_regions
43 44 45 46
        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(
47 48
                self.subdomains, self.materials, self.dim
            )
49
        else:
50 51 52
            raise TypeError(
                "materials parameter must be an instance of Material or a dictionnary that contains Material instances."
            )
53 54 55 56 57

    def mat_area(self):
        try:
            return self._mat_area
        except AttributeError:
58
            self._mat_area = fe.assemble(fe.Constant(1) * fe.dx(self.mesh))
59 60 61 62 63 64
            return self._mat_area

    def global_area(self):
        if not self.global_dimensions is None:
            return np.linalg.det(self.global_dimensions)
        else:
65 66 67
            raise AttributeError(
                f"global_size information is lacking for FenicsPart {self}."
            )
68 69

    @staticmethod
70 71 72 73 74 75 76 77
    def file_2_FenicsPart(
        mesh_path,
        materials,
        global_dimensions=None,
        subdomains_import=False,
        plots=True,
        explicit_subdo_val=0,
    ):
78
        """Generate an instance of FenicsPart from a .xml or .msh file that contains the mesh.
79 80 81 82

        Parameters
        ----------
        mesh_path : string or Path
83
            Relative or absolute path to the mesh file (format Dolfin XML, XDMF or MSH version 2)
84 85
        global_dimensions : 2D-array
            shape : 2×2 if 2D problem
86
            dimensions of the RVE
87
        materials : dict or Material
88 89 90 91 92 93 94 95
            [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
        -------
96
        FenicsPart instance
97 98 99 100 101

        """
        if not isinstance(mesh_path, Path):
            mesh_path = Path(mesh_path)
        name = mesh_path.stem
102 103 104 105 106 107 108 109 110 111 112
        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)
113
            suffix = mesh_path.suffix
114 115

        # Each supported mesh format -> one if structure
Baptiste Durand's avatar
Baptiste Durand committed
116
        subdomains, facets = None, None
117
        if suffix == ".xml":
118
            mesh = fe.Mesh(str(mesh_path))
Baptiste Durand's avatar
Baptiste Durand committed
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
            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."
                    )
141 142 143 144 145 146 147 148 149 150 151 152 153

        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"
                )

154 155
        if plots:
            plt.figure()
156
            fe.plot(mesh, title=f"Mesh of {name}")
157 158
            if subdomains_import:
                plt.figure()
159 160 161
                subdo_plt = fe.plot(
                    subdomains, mesh, title=f"Subdomains imported for {name}"
                )
162 163
                plt.colorbar(subdo_plt)
                plt.figure()
164 165
                facets_val_range = max(facets_val[1]) - min(facets_val[1])
                cmap = plt.cm.get_cmap("viridis", facets_val_range)
166 167
                facets_plt = fetools.facet_plot2d(facets, mesh, cmap=cmap)
                plt.colorbar(facets_plt[0])
168
            plt.show()
169
        logger.info(f"Import of the mesh : DONE")
170

171
        return FenicsPart(mesh, materials, subdomains, global_dimensions, facets)
172

173

174 175 176 177
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.
    """
178 179 180 181

    def __init__(
        self, mesh, generating_vectors, material_dict, subdomains, facet_regions
    ):
182 183 184 185 186 187 188
        """
        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
        """
189
        self.mesh = mesh
190 191
        self.gen_vect = generating_vectors
        self.rve_area = np.linalg.det(self.gen_vect)
192 193
        self.mat_area = fe.assemble(fe.Constant(1) * fe.dx(mesh))
        self.mesh_dim = mesh.topology().dim()  # dimension d'espace de depart
194 195 196
        self.materials = material_dict
        self.subdomains = subdomains
        self.facet_regions = facet_regions
197

198 199 200 201
        if isinstance(self.materials, mat.Material):
            self.C_per = fe.as_matrix(self.materials.get_C())
        else:
            self.C_per = mat.mat_per_subdomains(
202 203
                self.subdomains, self.materials, self.mesh_dim
            )
204

205
    def epsilon(self, u):
206
        return mat.epsilon(u)
207

208
    def sigma(self, eps):
209
        return mat.sigma(self.C_per, eps)
210

211 212 213
    def StrainCrossEnergy(self, sig, eps):
        return mat.strain_cross_energy(sig, eps, self.mesh, self.rve_area)

214
    @staticmethod
215
    def gmsh_2_Fenics_2DRVE(gmsh_2D_RVE, material_dict, plots=True):
216 217
        """
        Generate an instance of Fenics2DRVE from a instance of the Gmsh2DRVE class.
218

219
        """
220 221 222 223 224 225 226
        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()
        )
227 228 229 230
        run(cmd, shell=True, check=True)
        mesh = fe.Mesh(xml_path.as_posix())
        name = xml_path.stem
        subdomains = fe.MeshFunction(
231 232
            "size_t", mesh, xml_path.with_name(name + "_physical_region.xml").as_posix()
        )
233
        facets = fe.MeshFunction(
234 235
            "size_t", mesh, xml_path.with_name(name + "_facet_region.xml").as_posix()
        )
236 237
        subdo_val = fetools.get_MeshFunction_val(subdomains)
        facets_val = fetools.get_MeshFunction_val(facets)
238 239 240 241 242 243
        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]}"
        )
244
        if plots:
245 246 247 248
            plt.figure()
            subdo_plt = fe.plot(subdomains)
            plt.colorbar(subdo_plt)
            plt.figure()
249
            cmap = plt.cm.get_cmap("viridis", max(facets_val[1]) - min(facets_val[1]))
250
            facets_plt = fetools.facet_plot2d(facets, mesh, cmap=cmap)
251
            plt.colorbar(facets_plt[0])
252
            # clrbar.set_ticks(facets_val)
253
            plt.draw()
254
        logger.info(f"Import of the mesh : DONE")
255

256 257 258
        generating_vectors = gmsh_2D_RVE.gen_vect
        return Fenics2DRVE(mesh, generating_vectors, material_dict, subdomains, facets)

259 260
    @staticmethod
    def file_2_Fenics_2DRVE(mesh_path, generating_vectors, material_dict, plots=True):
261 262
        """Generate an instance of Fenics2DRVE from a .xml, .msh or .xdmf
        file that contains the mesh.
263

264 265 266
        Parameters
        ----------
        mesh_path : string or Path
267 268
            Relative or absolute path to the mesh file
            (format Dolfin XML, XDMF or MSH version 2)
269 270 271 272 273 274 275 276 277 278
        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
279

280 281 282 283 284
        """
        if not isinstance(mesh_path, Path):
            mesh_path = Path(mesh_path)
        name = mesh_path.stem

285
        if mesh_path.suffix == ".xml":
286
            mesh = fe.Mesh(str(mesh_path))
287
        elif mesh_path.suffix == ".xdmf":
288 289 290 291
            mesh = fe.Mesh()
            with fe.XDMFFile(str(mesh_path)) as file_in:
                file_in.read(mesh)
        else:
292
            cmd = f"dolfin-convert {mesh_path} {mesh_path.with_suffix('.xml')}"
293
            run(cmd, shell=True, check=True)
294
            mesh_path = mesh_path.with_suffix(".xml")
Baptiste Durand's avatar
Baptiste Durand committed
295
            mesh = fe.Mesh(str(mesh_path))
296

297 298
        subdo_path = mesh_path.with_name(name + "_physical_region.xml")
        facet_path = mesh_path.with_name(name + "_facet_region.xml")
299
        if subdo_path.exists():
300
            subdomains = fe.MeshFunction("size_t", mesh, subdo_path.as_posix())
301
            subdo_val = fetools.get_MeshFunction_val(subdomains)
302 303 304
            logger.info(
                f"{subdo_val[0]} physical regions imported. Tags : {subdo_val[1]}"
            )
305
        else:
306 307 308
            logger.info(
                f"For mesh file {mesh_path.name}, _physical_region.xml file is missing."
            )
309 310
            subdomains = None
        if facet_path.exists():
311
            facets = fe.MeshFunction("size_t", mesh, facet_path.as_posix())
312
            facets_val = fetools.get_MeshFunction_val(facets)
313 314 315
            logger.info(
                f"{facets_val[0]} facet regions imported. Tags : {facets_val[1]}"
            )
316
        else:
317 318 319
            logger.info(
                f"For mesh file {mesh_path.name}, _facet_region.xml file is missing."
            )
320
            facets = None
321

322 323
        if plots:
            plt.figure()
324 325 326 327 328 329 330
            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()
331 332 333
                cmap = plt.cm.get_cmap(
                    "viridis", max(facets_val[1]) - min(facets_val[1])
                )
334 335
                facets_plt = fetools.facet_plot2d(facets, mesh, cmap=cmap)
                plt.colorbar(facets_plt[0])
336
            plt.draw()
337
        logger.info(f"Import of the mesh : DONE")
338

339
        return Fenics2DRVE(mesh, generating_vectors, material_dict, subdomains, facets)
340

341

342
if __name__ == "__main__":
343 344
    pass
    # geo.init_geo_tools()
345

346 347
    # a = 1
    # b, k = a, a/3
348
    # panto_test = Gmsh2DRVE.pantograph(a, b, k, 0.1, nb_cells=(2, 3), soft_mat=False, name='panto_test')
349 350 351 352
    # 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 &")

353 354 355 356 357 358
    # 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 &")
359

360 361 362 363 364 365 366
    # 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 &")
367

368 369 370 371 372 373 374 375 376
    # 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 &")
377

378 379
    # gmsh.option.setNumber('Mesh.SurfaceFaces',1) #Display faces of surface mesh?
    # gmsh.fltk.run()
380 381 382 383 384 385 386 387 388 389

# 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 &")