__init__.py 19.5 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13
# coding: utf-8
"""
Created on 17/01/2019
@author: Baptiste Durand, baptiste.durand@enpc.fr

Tools and classes designed to create 2D gmsh model of a periodic metamaterial represented on a certain domain.

The Gmsh2DRVE class has been specially designed for the generation of meshes that represent RVE with a given microstructure.
"""

import copy
import logging
import math
14
from itertools import chain, product
15
from pathlib import Path, PurePath
16

Baptiste Durand's avatar
Baptiste Durand committed
17
import gmsh
18
import numpy as np
19
from more_itertools import flatten, one, padded
20

21 22
import ho_homog.geometry as geo
import ho_homog.mesh_tools as msh
23 24 25 26 27 28

# nice shortcuts
model = gmsh.model
factory = model.occ


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


32 33
__all__ = [
    "pantograph",
34
    "kagome",
35 36 37 38 39 40 41 42
    "duplicate_pattern",
    "offset_pattern",
    "Gmsh2DRVE",
    "Gmsh2DPart",
    "Gmsh2DPartFromRVE",
]


43 44 45 46 47 48 49 50 51 52 53 54
def duplicate_pattern(cell_ll, nb_cells, gen_vect):
    """
    Propager la cellule de base dans selon les directions de l'espace.

    Parameters
    ----------
    cell_ll : list
        The LineLoops that will be replicated. They define the pattern for a unit cell.
    np_cells : tuple, dimension 2 or 3
        Number of cells in each direction.
    gen_vect : array
        The generating vectors that are related to the given microstruture.
55
        Vecteurs de périodicité, en colonnes
56 57 58 59 60 61 62

    Returns
    -------
    repeated_ll : list
        Repeated lineloops that define the pattern over the whole domain associated with the given numbers of cells.
    """
    repeated_ll = cell_ll
63 64
    if gen_vect.shape != (3, 3):
        gen_vect_3D = np.zeros((3, 3))
65
        gen_vect_3D[: gen_vect.shape[0], : gen_vect.shape[1]] = gen_vect
66
    else:
67 68
        gen_vect_3D = gen_vect
    for k in range(len(nb_cells)):
69
        if nb_cells[k] > 1:
70 71
            new_contours = list()
            for i in range(1, int(nb_cells[k])):
72
                new_contours += [
73
                    geo.translation(ll, i * gen_vect_3D[:, k]) for ll in repeated_ll
74
                ]
75 76 77 78
            repeated_ll += new_contours
    repeated_ll = geo.remove_duplicates(repeated_ll)
    return repeated_ll

79

80 81 82 83 84 85 86 87 88 89 90 91
def offset_pattern(cell_ll, offset, cell_vect):
    """
    Translation of the lineloops that define the microstructure geometry of a unit cell.

    Parameters
    ----------
    cell_ll : list of instances of LineLoop
    offset : 1D array
        relative coordinates with respect to the unit-cell generating vectors of the point that will be moved to the origin
    gen_vect : 2D array
        The generating vectors that are related to the given unit-cell.
    """
92 93 94
    if cell_vect.shape != (3, 3):
        cell_vect_3D = np.zeros((3, 3))
        cell_vect_3D[: cell_vect.shape[0], : cell_vect.shape[1]] = cell_vect
95
    else:
96 97 98
        cell_vect_3D = cell_vect
    offset_vect_relat = np.zeros(3)
    for i, val in enumerate(offset):
99
        offset_vect_relat[i] = val % 1.0
100 101 102 103 104
    offset_vect_abs = np.dot(cell_vect_3D, offset_vect_relat)
    t_vect = -1 * offset_vect_abs
    shifted_ll = [geo.translation(ll, t_vect) for ll in cell_ll]
    return shifted_ll

105

106
class Gmsh2DRVE(object):
107 108 109
    def __init__(
        self, pattern_ll, cell_vect, nb_cells, offset, attractors, soft_mat, name
    ):
110 111
        """
        Contrat : Créer un maillage pour des RVE 2D, plans, comportant au plus 2 matériaux constitutifs et pouvant contenir plusieurs cellules.
112 113 114

        cell_vect : Vecteurs de périodicité en colonne, définit également le parallélogramme qui contient la cellule.
        #! La cellule doit être un parallélogramme.
115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131

        Parameters
        ----------
        pattern_ll : list
            Instances of LineLoop that define the contours of the microstructure.
        cell_vect : 2D array
            dimensions of the unit cell and directions of periodicity.
            (given in a 2D cartesian coordinate system)
        nb_cells : 1D array
            Numbers of cells in each direction of repetition/periodicity.
        offset : 1D array
            Relative position inside a cell of the point that will coincide with the origin of the global domain.
        attractors : list
            Instances of Point.
            Can also be = None or empty.
            It represent the points that will be used as attractors in the definition of the element characteristic length fields.
            Attractors are geometrical elements of the cell around which mesh refinement constraints will be set.
132
        name : string or Path
133
        """
134 135

        self.name = name.stem if isinstance(name, PurePath) else name
136 137 138 139
        model.add(self.name)
        model.setCurrent(self.name)

        if offset.any():
140 141 142 143
            nb_pattern = [
                math.ceil(val + 1) if offset[i] != 0 else math.ceil(val)
                for i, val in enumerate(nb_cells)
            ]
144 145 146 147
            nb_pattern = np.array(nb_pattern, dtype=np.int8)
            pattern_ll = offset_pattern(pattern_ll, offset, cell_vect)
        else:
            nb_pattern = np.int8(np.ceil(nb_cells))
148

149
        if not np.equal(nb_pattern, 1).all():
150 151
            duplicate_pattern(pattern_ll, nb_pattern, cell_vect)

152
        rve_vect = cell_vect * np.asarray(nb_cells)
153 154 155 156 157 158 159 160
        pt_o = np.zeros((3,))
        macro_vtcs = [
            pt_o,
            rve_vect[:, 0],
            rve_vect[:, 0] + rve_vect[:, 1],
            rve_vect[:, 1],
        ]
        macro_ll = geo.LineLoop([geo.Point(c) for c in macro_vtcs])
161 162 163 164 165 166
        macro_s = geo.PlaneSurface(macro_ll)

        if attractors:
            for entity in attractors:
                if not isinstance(entity, geo.Point):
                    raise TypeError(
167 168 169
                        """Use of curves as attractors for the refinement of the mesh
                    is not yet fully supported in our python library for gmsh."""
                    )
170 171 172 173 174
            if offset.any():
                attractors = offset_pattern(attractors, offset, cell_vect)
            if not np.equal(nb_pattern, 1).all():
                duplicate_pattern(attractors, nb_pattern, cell_vect)

175
        logger.info("Start boolean operations on surfaces")
176 177
        phy_surf = list()
        pattern_s = [geo.PlaneSurface(ll) for ll in pattern_ll]
178
        rve_s = geo.surface_bool_cut(macro_s, pattern_s)
179
        if len(rve_s) == 1:
180 181 182
            logger.info(
                "The main material domain of the RVE is connected (topological property)."
            )
183
        elif len(rve_s) == 0:
184 185 186
            logger.warning(
                "The boolean operation for creating the main material domain of the RVE return 0 surfaces."
            )
187
        else:
188 189 190
            logger.warning(
                "The main material domain of the RVE obtained by a boolean operation is disconnected (topological property)."
            )
191 192 193
        rve_s_phy = geo.PhysicalGroup(rve_s, 2, "microstruct_domain")
        phy_surf.append(rve_s_phy)
        if soft_mat:
194
            soft_s = geo.surface_bool_cut(macro_s, rve_s)
195 196
            soft_s_phy = geo.PhysicalGroup(soft_s, 2, "soft_domain")
            phy_surf.append(soft_s_phy)
197
        logger.info("Done boolean operations on surfaces")
198 199 200 201

        if attractors:
            need_sync = False
            for entity in attractors:
202 203 204
                if not entity.tag:
                    entity.add_gmsh()
                    need_sync = True
205
            if need_sync:
206
                factory.synchronize()  # ? Pourrait être enlevé ?
207 208

        for gp in phy_surf:
209 210 211 212
            gp.add_gmsh()
        factory.synchronize()

        data = model.getPhysicalGroups()
213 214 215 216 217 218 219
        details = [
            f"Physical group id : {dimtag[1]}, "
            + f"dimension : {dimtag[0]}, "
            + f"name : {model.getPhysicalName(*dimtag)}, "
            + f"nb of entitities {len(model.getEntitiesForPhysicalGroup(*dimtag))} \n"
            for dimtag in data
        ]
220 221
        logger.debug(f"All physical groups in the model : {data}")
        logger.debug(f"Physical groups details : \n {details}")
222
        logger.info("Done generating the gmsh geometrical model")
223
        if isinstance(name, PurePath):
224
            gmsh.write(str(name.with_suffix(".brep")))
225
        else:
226
            gmsh.write(f"{name}.brep")
227 228
        macro_bndry = macro_ll.sides
        if soft_mat:
229
            boundary = geo.AbstractSurface.get_surfs_boundary(rve_s + soft_s)
230 231 232
        else:
            try:
                s = one(rve_s)
233
                boundary = geo.AbstractSurface.get_surfs_boundary(s, recursive=False)
234
            except ValueError:
235 236 237
                boundary = geo.AbstractSurface.get_surfs_boundary(
                    rve_s, recursive=False
                )
Baptiste Durand's avatar
Baptiste Durand committed
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
        for b in boundary:
            b.get_boundary(get_coords=True)
        # factory.synchronize()
        micro_bndry = list()
        for macro_line in macro_bndry:
            fragments = geo.macro_line_fragments(boundary, macro_line)
            for c in fragments:
                c.gmsh_type = model.getType(1, c.tag)
            lines_only = [c for c in fragments if c.gmsh_type == "Line"]
            micro_bndry.append(lines_only)
        macro_dir = list()
        for i in range(len(macro_bndry) // 2):
            macro_line = macro_bndry[i]
            direction = macro_line.def_pts[-1].coord - macro_line.def_pts[0].coord
            macro_dir.append(direction)
253
        for i, crvs in enumerate(micro_bndry):
254
            msh.order_curves(crvs, macro_dir[i % 2], orientation=True)
255 256
        msh.set_periodicity_pairs(micro_bndry[0], micro_bndry[2])
        msh.set_periodicity_pairs(micro_bndry[1], micro_bndry[3])
257 258 259 260 261 262 263
        logger.info("Done defining a mesh periodicity constraint")
        tags = [
            "per_pair_1_slave",
            "per_pair_2_slave",
            "per_pair_1_mast",
            "per_pair_2_mast",
        ]
264 265 266 267 268
        per_pair_phy = list()
        for crvs, tag in zip(micro_bndry, tags):
            per_pair_phy.append(geo.PhysicalGroup(crvs, 1, tag))
        for gp in per_pair_phy:
            gp.add_gmsh()
269

270 271 272 273
        self.gen_vect = rve_vect
        self.nb_cells = nb_cells
        self.attractors = attractors if attractors else []
        self.phy_surf = phy_surf
274
        self.mesh_fields = []
275
        self.mesh_abs_path = ""
276 277 278

    def main_mesh_refinement(self, d_min_max, lc_min_max, sigmoid_interpol=False):
        model.setCurrent(self.name)
279 280 281 282
        attractors = {"points": self.attractors}
        logger.debug(
            f"When main_mesh_refinement(...) is called, physical groups in model : {model.getPhysicalGroups()}"
        )
283 284 285 286 287
        rve_s = self.phy_surf[0].entities
        for s in rve_s:
            if not s.boundary:
                s.get_boundary()
        rve_boundary = list(flatten([s.boundary for s in rve_s]))
288 289 290 291
        restrict_domain = {"surfaces": rve_s, "curves": rve_boundary}
        field = msh.set_mesh_refinement(
            d_min_max, lc_min_max, attractors, 10, sigmoid_interpol, restrict_domain
        )
292
        self.mesh_fields.append(field)
293 294 295

    def soft_mesh_refinement(self, d_min_max, lc_min_max, sigmoid_interpol=False):
        model.setCurrent(self.name)
296
        attractors = {"points": self.attractors}
297 298 299 300 301
        soft_s = self.phy_surf[1].entities
        for s in soft_s:
            if not s.boundary:
                s.get_boundary()
        soft_boundary = list(flatten([s.boundary for s in soft_s]))
302 303 304 305
        restrict_domain = {"surfaces": soft_s, "curves": soft_boundary}
        field = msh.set_mesh_refinement(
            d_min_max, lc_min_max, attractors, 1, sigmoid_interpol, restrict_domain
        )
306
        self.mesh_fields.append(field)
307

308
    def mesh_generate(self, mesh_field=None, directory: Path = None):
309
        """Generate a 2D mesh of the model which represent a RVE.
310

311 312 313 314 315
        Parameters
        ----------
        mesh_field : mesh_tools.Field, optional
            The characteristic length of the elements can be explicitly prescribe by means of this field.
            The default is None. In this case, the fields that have been created with the soft_mesh_refinement and main_mesh_refinement methods are used.
316 317 318 319
        directory : pathlib.Path, optional
            Indicate in which directory the .msh file must be created.
            If None (default), the .msh file is written in the current working directory.
            Ex: Path('/media/sf_VM_share/homog')
320
        """
321

322 323 324 325 326 327
        model.setCurrent(self.name)
        if not mesh_field:
            self.background_field = msh.set_background_mesh(self.mesh_fields)
        else:
            self.background_field = msh.set_background_mesh(mesh_field)
        data = model.getPhysicalGroups()
328
        logger.debug(f"Physical groups in model just before generating mesh : {data}")
329 330 331 332 333 334
        geo.PhysicalGroup.set_group_mesh(True)
        model.mesh.generate(1)
        gmsh.model.mesh.removeDuplicateNodes()
        model.mesh.generate(2)
        gmsh.model.mesh.removeDuplicateNodes()
        geo.PhysicalGroup.set_group_visibility(False)
335
        if directory:
336
            mesh_path = directory if not directory.suffix else directory.with_suffix("")
Baptiste Durand's avatar
Baptiste Durand committed
337
            if not mesh_path.exists():
338 339 340 341 342 343
                mesh_path.mkdir(mode=0o777, parents=True)
        else:
            mesh_path = Path.cwd()
        mesh_path = mesh_path.joinpath(f"{self.name}.msh")
        gmsh.write(str(mesh_path))
        self.mesh_abs_path = mesh_path
344

345 346 347 348 349 350 351 352 353

class Gmsh2DPart(object):
    def __init__(self, gen_vect, nb_cells, phy_surf, mesh_path: PurePath):
        self.gen_vect = gen_vect
        self.nb_cells = nb_cells
        self.phy_surf = phy_surf
        self.mesh_abs_path = mesh_path.resolve()


354
from .pantograph import pantograph_RVE, pantograph_offset_RVE, beam_pantograph_RVE
355
from .kagome import kagome_RVE
356

357

358 359
# from .other_2d_microstructures import auxetic_square_RVE

360

361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396
def Gmsh2DPartFromRVE(cell: Gmsh2DRVE, nb_cells, part_name=None):
    """[summary]

    Parameters
    ----------
    rve : Gmsh2DRVE
        [description]
    nb_cells : tuple, dimension 2 or 3
        Number of cells in each direction.
    part_name: str, optional
        Desired name for the mesh file
    Returns
    -------
    tuple
        Paths to the RVE mesh and the part mesh

    Remarques
    ----------
    Pour le moment, le RVE est composé d'une unique cellule.
    Pour le moment, le domaine macro est un parallélogramme aligné avec les axes
    de la cellule et contient un nombre entier de cellules.
    """

    name = cell.name
    cell_vect = cell.gen_vect
    # * 2D -> 3D
    if cell_vect.shape != (3, 3):
        cell_vect_3D = np.zeros((3, 3))
        cell_vect_3D[: cell_vect.shape[0], : cell_vect.shape[1]] = cell_vect
        cell_vect = cell_vect_3D
    if len(nb_cells) != 3:
        nb_cells = tuple(padded(nb_cells, 1, 3))

    # TODO : Activer le model gmsh correspondant au RVE
    model.setCurrent(name)
    # TODO : Créer un domaine macro
397
    part_vect = cell_vect * np.asarray(nb_cells)
398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453
    macro_vertices = [
        np.zeros((3,)),
        part_vect[0],
        part_vect[0] + part_vect[1],
        part_vect[1],
    ]
    macro_lloop = geo.LineLoop([geo.Point(c) for c in macro_vertices])
    macro_surf = geo.PlaneSurface(macro_lloop)

    translat_vectors = list()
    for translat_combination in product(*[range(i) for i in nb_cells]):
        if translat_combination == (0, 0, 0):
            continue  # Correspond à la cellule de base
        # * cell_vect : vectors in column
        t_vect = np.dot(cell_vect, np.array(translat_combination))
        translat_vectors.append(t_vect)
    # ? Exemple :
    # ? >>> nb_cells = (2,3)
    # ? >>> nb_cells = tuple(padded(nb_cells,1,3))
    # ? >>> nb_cells
    # * (2, 3, 1)
    # ? >>> translat_vectors = list()
    # ? >>> cell_vect = np.array(((4.,1.,0.),(3.,8.,0.),(0.,0.,0.)))
    # ? >>> for translat_combination in product(*[range(i) for i in nb_cells]):
    # ?         t_vect = np.dot(cell_vect, np.array(translat_combination))
    # ?         translat_vectors.append(t_vect)
    # ? >>> translat_vectors
    # * [array([0., 0., 0.]), array([1., 8., 0.]), array([ 2., 16.,  0.]), array([4., 3., 0.]), array([ 5., 11.,  0.]), array([ 6., 19.,  0.])] #noqa
    cell_surfaces_by_gp = [phy_surf.entities for phy_surf in cell.phy_surf]
    repeated_surfaces_by_gp = [list() for i in range(len(cell.phy_surf))]
    # * Structure de repeated_surfaces_by_gp :
    # * List avec :
    # *     pour chaque physical group, et pour chaque translation
    # *          la liste des surfaces (entitées) translatées

    for i, gp_surfaces in enumerate(cell_surfaces_by_gp):
        for t_vect in translat_vectors:
            dimTags = factory.copy([(2, s.tag) for s in gp_surfaces])
            factory.translate(dimTags, *(t_vect.tolist()))
            # ? Opération booléenne d'intersection ?
            # ? Pour détecter si surface entière : comparaison de boundingbox surface d'origine et bounding box resultat - vecteur translation
            this_translation_surfs = [geo.AbstractSurface(dt[1]) for dt in dimTags]
            repeated_surfaces_by_gp[i].append(this_translation_surfs)
    factory.synchronize()
    # TODO : Contraintes de périodicité
    for j, t_vect in enumerate(translat_vectors):
        master = list(chain.from_iterable(cell_surfaces_by_gp))
        all_surfs_this_transl = [surfs[j] for surfs in repeated_surfaces_by_gp]
        slaves = list(chain.from_iterable(all_surfs_this_transl))
        msh.set_periodicity_pairs(slaves, master, t_vect)
    # TODO : Extension des physical groups
    phy_surfaces = list()
    for i in range(len(cell.phy_surf)):
        all_surfaces = cell_surfaces_by_gp[i] + list(
            flatten(repeated_surfaces_by_gp[i])
        )
454 455 456 457 458 459
        tag = cell.phy_surf[i].tag + 1000
        name = cell.phy_surf[i].name
        phy_surfaces.append(geo.PhysicalGroup(all_surfaces, 2, name, tag))
    for gp in cell.phy_surf:
        gp.remove_gmsh()
    factory.synchronize()
460 461
    for gp in phy_surfaces:
        gp.add_gmsh()
manon222's avatar
manon222 committed
462 463 464 465
    all_gp = model.getPhysicalGroups()
    dimtags_part = [(gp.dim, gp.tag) for gp in phy_surfaces]
    remove_gp = [dt for dt in all_gp if not dt in dimtags_part]
    model.removePhysicalGroups(remove_gp)
466 467 468 469 470 471 472 473 474
    # ! Pour le moment, il semble impossible de réutiliser le tag d'un physical group
    # ! qui a été supprimé.
    # ! Voir : \Experimental\Test_traction_oct19\pb_complet\run_3\MWE_reuse_tag.py
    # ! Autre solution :
    # ! - Compléter les physical group existants ?
    # !      Impossible car groups déjà ajoutés au model
    # ! Utiliser un autre tag, avec une règle pour relier les 2.
    # !     Solution retenue. Règle choisie : le tag pour la part = 1000 + tag pour la cell

475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504
    # TODO : All mesh generation
    geo.PhysicalGroup.set_group_mesh(True)
    model.mesh.generate(1)
    model.mesh.generate(2)
    model.mesh.removeDuplicateNodes()
    gmsh.model.mesh.renumberNodes()
    gmsh.model.mesh.renumberElements()
    geo.PhysicalGroup.set_group_visibility(False)

    rve_path = cell.mesh_abs_path
    conversion = {
        "RVE": "PART",
        "rve": "part",
        "Rve": "Part",
        "cell": "part",
        "CELL": "PART",
        "Cell": "Part",
    }
    if part_name:
        part_path = rve_path.with_name(part_name).with_suffix(".msh")
    elif any(x in rve_path.name for x in conversion.keys()):
        name = rve_path.name
        for old, new in conversion.items():
            name = name.replace(old, new)
        part_path = rve_path.with_name(name).with_suffix(".msh")
    else:
        part_path = rve_path.with_name(rve_path.stem + "_part.msh")
    gmsh.write(str(part_path.with_suffix(".brep")))
    gmsh.write(str(part_path))
    return Gmsh2DPart(part_vect, nb_cells, phy_surfaces, part_path)
505

506

manon222's avatar
manon222 committed
507
from . import kagome
508
from . import pantograph
509 510 511

__all__ = [
    "pantograph_RVE",
manon222's avatar
manon222 committed
512
    "kagome_RVE",
513 514
    "pantograph_offset_RVE",
    "beam_pantograph_RVE",
515
    "pantograph_E11only_RVE",
516 517 518 519 520
    "auxetic_square_RVE",
    "Gmsh2DRVE",
    "Gmsh2DPart",
    "Gmsh2DPartFromRVE",
]