__init__.py 18.8 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 34 35 36 37 38 39 40 41
__all__ = [
    "pantograph",
    "duplicate_pattern",
    "offset_pattern",
    "Gmsh2DRVE",
    "Gmsh2DPart",
    "Gmsh2DPartFromRVE",
]


42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
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.

    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
61 62
    if gen_vect.shape != (3, 3):
        gen_vect_3D = np.zeros((3, 3))
63
        gen_vect_3D[: gen_vect.shape[0], : gen_vect.shape[1]] = gen_vect
64
    else:
65 66
        gen_vect_3D = gen_vect
    for k in range(len(nb_cells)):
67
        if nb_cells[k] > 1:
68 69
            new_contours = list()
            for i in range(1, int(nb_cells[k])):
70 71 72
                new_contours += [
                    geo.translation(ll, i * gen_vect_3D[k]) for ll in repeated_ll
                ]
73 74 75 76
            repeated_ll += new_contours
    repeated_ll = geo.remove_duplicates(repeated_ll)
    return repeated_ll

77

78 79 80 81 82 83 84 85 86 87 88 89
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.
    """
90 91 92
    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
93
    else:
94 95 96
        cell_vect_3D = cell_vect
    offset_vect_relat = np.zeros(3)
    for i, val in enumerate(offset):
97
        offset_vect_relat[i] = val % 1.0
98 99 100 101 102
    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

103

104
class Gmsh2DRVE(object):
105 106 107
    def __init__(
        self, pattern_ll, cell_vect, nb_cells, offset, attractors, soft_mat, name
    ):
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
        """
        Contrat : Créer un maillage pour des RVE 2D, plans, comportant au plus 2 matériaux constitutifs et pouvant contenir plusieurs cellules.
        #! La cellule est un parallélogramme.

        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.
128
        name : string or Path
129
        """
130 131

        self.name = name.stem if isinstance(name, PurePath) else name
132 133 134 135
        model.add(self.name)
        model.setCurrent(self.name)

        if offset.any():
136 137 138 139
            nb_pattern = [
                math.ceil(val + 1) if offset[i] != 0 else math.ceil(val)
                for i, val in enumerate(nb_cells)
            ]
140 141 142 143
            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))
144

145
        if not np.equal(nb_pattern, 1).all():
146 147
            duplicate_pattern(pattern_ll, nb_pattern, cell_vect)

148
        rve_vect = cell_vect * nb_cells[:, np.newaxis]
149
        O = np.zeros((3,))
150
        macro_vtx = [O, rve_vect[0], rve_vect[0] + rve_vect[1], rve_vect[1]]
151 152 153 154 155 156 157
        macro_ll = geo.LineLoop([geo.Point(c) for c in macro_vtx])
        macro_s = geo.PlaneSurface(macro_ll)

        if attractors:
            for entity in attractors:
                if not isinstance(entity, geo.Point):
                    raise TypeError(
158 159 160
                        """Use of curves as attractors for the refinement of the mesh
                    is not yet fully supported in our python library for gmsh."""
                    )
161 162 163 164 165
            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)

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

        if attractors:
            need_sync = False
            for entity in attractors:
193 194 195
                if not entity.tag:
                    entity.add_gmsh()
                    need_sync = True
196
            if need_sync:
197
                factory.synchronize()  # ? Pourrait être enlevé ?
198 199

        for gp in phy_surf:
200 201 202 203
            gp.add_gmsh()
        factory.synchronize()

        data = model.getPhysicalGroups()
204 205 206 207 208 209 210
        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
        ]
211 212
        logger.debug(f"All physical groups in the model : {data}")
        logger.debug(f"Physical groups details : \n {details}")
213
        logger.info("Done generating the gmsh geometrical model")
214
        if isinstance(name, PurePath):
215
            gmsh.write(str(name.with_suffix(".brep")))
216
        else:
217
            gmsh.write(f"{name}.brep")
218 219
        macro_bndry = macro_ll.sides
        if soft_mat:
220
            boundary = geo.AbstractSurface.get_surfs_boundary(rve_s + soft_s)
221 222 223 224 225 226 227 228
        else:
            try:
                s = one(rve_s)
                boundary = geo.AbstractSurface.get_surfs_boundary(s)
            except ValueError:
                boundary = geo.AbstractSurface.get_surfs_boundary(rve_s)
        factory.synchronize()
        micro_bndry = [geo.macro_line_fragments(boundary, M_ln) for M_ln in macro_bndry]
229 230 231 232
        macro_dir = [
            macro_bndry[i].def_pts[-1].coord - macro_bndry[i].def_pts[0].coord
            for i in range(len(macro_bndry) // 2)
        ]
233
        for i, crvs in enumerate(micro_bndry):
234
            msh.order_curves(crvs, macro_dir[i % 2], orientation=True)
235 236
        msh.set_periodicity_pairs(micro_bndry[0], micro_bndry[2])
        msh.set_periodicity_pairs(micro_bndry[1], micro_bndry[3])
237 238 239 240 241 242 243
        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",
        ]
244 245 246 247 248
        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()
249

250 251 252 253
        self.gen_vect = rve_vect
        self.nb_cells = nb_cells
        self.attractors = attractors if attractors else []
        self.phy_surf = phy_surf
254
        self.mesh_fields = []
255
        self.mesh_abs_path = ""
256 257 258

    def main_mesh_refinement(self, d_min_max, lc_min_max, sigmoid_interpol=False):
        model.setCurrent(self.name)
259 260 261 262
        attractors = {"points": self.attractors}
        logger.debug(
            f"When main_mesh_refinement(...) is called, physical groups in model : {model.getPhysicalGroups()}"
        )
263 264 265 266 267
        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]))
268 269 270 271
        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
        )
272
        self.mesh_fields.append(field)
273 274 275

    def soft_mesh_refinement(self, d_min_max, lc_min_max, sigmoid_interpol=False):
        model.setCurrent(self.name)
276
        attractors = {"points": self.attractors}
277 278 279 280 281
        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]))
282 283 284 285
        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
        )
286
        self.mesh_fields.append(field)
287

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

291 292 293 294 295
        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.
296 297 298 299
        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')
300
        """
301

302 303 304 305 306 307
        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()
308
        logger.debug(f"Physical groups in model just before generating mesh : {data}")
309 310 311 312 313 314
        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)
315
        if directory:
316
            mesh_path = directory if not directory.suffix else directory.with_suffix("")
Baptiste Durand's avatar
Baptiste Durand committed
317
            if not mesh_path.exists():
318 319 320 321 322 323
                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
324

325 326 327 328 329 330 331 332 333

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


334
from .pantograph import pantograph_RVE, pantograph_offset_RVE, beam_pantograph_RVE
335

336 337
# from .other_2d_microstructures import auxetic_square_RVE

338

339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 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 397 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
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
    part_vect = cell_vect * np.asarray(nb_cells)[:, np.newaxis]
    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])
        )
432 433 434 435 436 437 438
        tag = cell.phy_surf[i].tag + 1000
        name = cell.phy_surf[i].name
        phy_surfaces.append(geo.PhysicalGroup(all_surfaces, 2, name, tag))
    # gmsh.fltk.run()
    for gp in cell.phy_surf:
        gp.remove_gmsh()
    factory.synchronize()
439 440
    for gp in phy_surfaces:
        gp.add_gmsh()
Baptiste Durand's avatar
Baptiste Durand committed
441 442 443 444
    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)
445 446 447 448 449 450 451 452 453
    # ! 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

454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483
    # 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)
484

485

486
from . import pantograph
487 488 489 490 491

__all__ = [
    "pantograph_RVE",
    "pantograph_offset_RVE",
    "beam_pantograph_RVE",
492
    "pantograph_E11only_RVE",
493 494 495 496 497
    "auxetic_square_RVE",
    "Gmsh2DRVE",
    "Gmsh2DPart",
    "Gmsh2DPartFromRVE",
]