surfaces.py 17.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
# coding: utf8
"""
Created on 09/10/2018
@author: baptiste

Definition of the classes : LineLoop, PlaneSurface and AbstractSurface.

PlaneSurface:
Object designed to represent geometrical entitites of dimension two
and instantiate them in a gmsh model.
"""

manon222's avatar
manon222 committed
13 14
from .curves import Line, Arc
from .point import Point
15
from . import factory, np, logger, model
Baptiste Durand's avatar
black  
Baptiste Durand committed
16
from .tools import round_corner, offset, calcul_R, unit_vect, bisector, angle_between
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
from .curves import Line, AbstractCurve


def round_corner_2_sides(result_list):
    """ Permet de traiter les résultats d'une opération round_corner
    appliquée en série sur un ensemble de sommets.
    Une polyligne composée de segments et d'arc est composée.
    """
    sides = list()
    # ? à la bonne place ?
    for i, rslt in enumerate(result_list):
        new_line = rslt[0]
        new_line.def_pts[0] = result_list[i - 1][1].def_pts[-1]
        # Correction pour que le segment commence à la fin de l'arc précédent.
        new_arc = rslt[1]
        sides.extend([new_line, new_arc])
    return sides


class LineLoop(object):
    """
    Définit une courbe fermée, composée d'entitées géométriques 1D (Line, Arc...).
    """

    def __init__(self, elements, explicit=False):
        """
        La LineLoop peut être créée à partir :
            - d'une liste de sommets,
            - ou d'une liste d'objets Line/Arcs (explicit = True)
        """

        self.info_offset = False
        # ! A remplacer par quelque chose de mieux, comme l'utilisation de l'attribut "vertices" #noqa
        if explicit:
            self.sides = elements
            self.vertices = list()
        else:
            self.vertices = elements
            self.sides = list()
        self.tag = None

    def __eq__(self, other):
        """
         Opérateur de comparaison == surchargé pour les objets de la classe LineLoop
         Si la LineLoop n'est définie que par ses sommets :
            True ssi les listes de sommets sont égales, à un décalage d'indice près.
         Si la LineLoop est aussi définie par des Line/Arc :
            True ssi l'ensemble des éléments 1D qui composent la LineLoop est
            identique à celui de la LineLoop comparée.
         L'orientation est prise en compte.

         """
        if not isinstance(other, LineLoop):
            return False
        if self.sides or other.sides:
            # Si l'une des deux LineLoops est définie par des Line/Arc, comparaison au niveau de ces éléments. #noqa
            if len(self.sides) != len(other.sides):
                return False
            else:
                for elmt_1D in self.sides:
                    for other_elmt in other.sides:
                        test = elmt_1D == other_elmt
                        if test:
                            break
                    else:
Baptiste Durand's avatar
black  
Baptiste Durand committed
82 83
                        # Aucun break déclenché,
                        # i.e. si l'un des cote de la lineloop courante n'appartient pas à la LineLoop comparée
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106
                        return False
                else:
                    return True
        if len(self.vertices) != len(other.vertices):
            return False
        else:
            for shift in range(len(self.vertices)):
                if all(
                    p == other.vertices[i - shift] for i, p in enumerate(self.vertices)
                ):
                    return True
            else:
                return False

    def __ne__(self, other):
        return not self.__eq__(other)

    def plot2D(self, color="black"):
        """Représenter la polyligne dans un plot matplotlib.
        Disponible seulement en 2D pour l'instant."""
        if not self.sides:
            self.vertices_2_sides()
        for elmt in self.sides:
107
            elmt.plot2D(color)
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125

    def add_gmsh(self):
        if self.tag:
            return None
        if not self.sides:
            self.vertices_2_sides()
        for elmt in self.sides:
            if not elmt.tag:
                elmt.add_gmsh()
        self.tag = factory.addCurveLoop([elmt.tag for elmt in self.sides])
        return self.tag

    def reverse(self):
        self.sides.reverse()
        for elmt in self.sides:
            elmt.def_pts.reverse()
        self.vertices.reverse()

126
    def offset(self, t, method="vertex"):
127 128 129 130 131 132 133 134 135 136 137
        """ Opération d'offset appliquée sur tout les sommets de la LineLoop.

        Cette opération doit donc être faite assez tôt,
        avant que les Line/Arc composant la LineLoop soient créés.
        """
        assert not self.sides
        # Si on est déjà en présence de Lines, il est trop tard pour faire l'offset de cette façon #noqa
        new_vrtces = [None] * (len(self.vertices))
        self.info_offset = True
        for i in range(len(self.vertices)):
            new_vrtces[i - 1] = offset(
138 139 140 141 142
                self.vertices[i - 1],
                self.vertices[i - 2],
                self.vertices[i],
                t,
                method=method,
143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
            )
        self.offset_dpcmt = [
            np.linalg.norm(new.coord - prev.coord)
            for new, prev in zip(new_vrtces, self.vertices)
        ]
        # TODO  : regarder où c'est utiliser et si on peut avoir quelque chose de plus clair #noqa
        self.vertices = new_vrtces

    def round_corner_explicit(self, radii):
        """Opération d'arrondi des angles appliquée à tous les sommets du polygone.
        Les rayons sont indiqués de manière explicite, sous forme d'une liste.
        Liste de longueur 1 pour un rayon uniforme.
        """
        # TODOC
        if isinstance(radii, list):
            radii = [radii[i % len(radii)] for i in range(len(self.vertices))]
        else:
            radii = [radii] * (len(self.vertices))
        result_1D = list()
        for i in range(len(self.vertices)):
            result_1D.append(
                round_corner(
                    self.vertices[i - 1],
                    self.vertices[i - 2],
                    self.vertices[i],
                    radii[i - 1],
                    False,
                    False,
                )
            )
        self.sides = round_corner_2_sides(result_1D)

    def round_corner_incircle(self, radii):
        """ Opération d'arrondi des angles appliquée à tous les sommets du polygone.
        La méthode du cercle inscrit est utilisée.
        radii = liste de rayons à utiliser ou valeur (float) si rayon uniforme.
        Une liste de rayons de cercles inscrits peut être indiquée,
        liste de longueur 1 pour un rayon uniforme.
        Si la longueur de la liste de rayon est ni 1 ni égale au nombre de sommets,
        un modulo est utilisé.
        """
        # TODOC
        if isinstance(radii, list):
            effect_R = [radii[i % len(radii)] for i in range(len(self.vertices))]
        else:
            effect_R = [radii] * (len(self.vertices))
        if self.info_offset:
            effect_R = [
                R - offset_d for R, offset_d in zip(effect_R, self.offset_dpcmt)
            ]
            # ! ESSAI
        result_1D = list()
        for i in range(len(self.vertices)):
            result_1D.append(
                round_corner(
                    self.vertices[i - 1],
                    self.vertices[i - 2],
                    self.vertices[i],
                    effect_R[i - 1],
                    True,
                    False,
                )
            )
        self.sides = round_corner_2_sides(result_1D)

Baptiste Durand's avatar
black  
Baptiste Durand committed
208 209 210 211
    def round_corner_kagome(self, r, a, alpha):
        """ Opération d'arrondi des angles spécifique à la microstructure 'kagome',
        appliquée à tous les sommets du polygone.

manon222's avatar
manon222 committed
212
        alpha = paramètre d'ouverture du kagomé
Baptiste Durand's avatar
black  
Baptiste Durand committed
213
        a :
manon222's avatar
manon222 committed
214 215
        b : taille des éléments triangulaires
        theta : rotation du triangle à l'intérieur de la cellule de base
Baptiste Durand's avatar
black  
Baptiste Durand committed
216 217 218 219 220 221 222 223 224

        Parameters
        ----------
        r: float
            Rayon du cercle inscrit dans la jonction
        a: float
            taille de la cellule de base
        alpha: float
            paramètre d'ouverture de la microstructure
manon222's avatar
manon222 committed
225 226
        """

Baptiste Durand's avatar
black  
Baptiste Durand committed
227 228 229
        effect_R, phi_1, phi_2 = calcul_R(alpha, r, a)

        # ! ESSAI
manon222's avatar
manon222 committed
230 231 232 233 234
        result_1D = list()
        for i in range(len(self.vertices)):
            d2 = effect_R / np.sin(phi_2)
            d1 = effect_R / np.sin(phi_1)

Baptiste Durand's avatar
black  
Baptiste Durand committed
235 236 237
            dir_1 = self.vertices[i - 2].coord - self.vertices[i - 1].coord
            dir_2 = self.vertices[i].coord - self.vertices[i - 1].coord
            dir_1 = unit_vect(dir_1)
manon222's avatar
manon222 committed
238 239
            dir_2 = unit_vect(dir_2)

Baptiste Durand's avatar
black  
Baptiste Durand committed
240 241
            A = Point(self.vertices[i - 1].coord + effect_R * dir_1)
            C = Point(self.vertices[i - 1].coord + effect_R * dir_2)
manon222's avatar
manon222 committed
242 243

            alpha = angle_between(dir_1, dir_2, orient=True)
Baptiste Durand's avatar
black  
Baptiste Durand committed
244 245 246 247 248 249 250
            v_biss = bisector(dir_1, dir_2) if alpha >= 0 else - bisector(dir_1, dir_2)

            if abs(abs(angle_between(v_biss, dir_1)) - (np.pi / 2 - phi_2)) < 10e-14:
                # si on est du côté où l'angle vaut theta
                d = d2
            elif abs(abs(angle_between(v_biss, dir_1)) - (np.pi / 2 - phi_1)) < 10e-14:
                d = d1
manon222's avatar
manon222 committed
251 252
            else:
                raise ValueError("mauvaise gestion de d1 et d2")
Baptiste Durand's avatar
black  
Baptiste Durand committed
253
            B = Point(self.vertices[i - 1].coord + d * v_biss)
manon222's avatar
manon222 committed
254 255 256
            round_arc = Arc(A, B, C)
            racc_amt = Line(self.vertices[i - 2], A)
            racc_avl = Line(C, self.vertices[i])
Baptiste Durand's avatar
black  
Baptiste Durand committed
257 258
            curves_list = [racc_amt, round_arc, racc_avl]
            result_1D.append(curves_list)
manon222's avatar
manon222 committed
259 260
        self.sides = round_corner_2_sides(result_1D)

261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 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
    def vertices_2_sides(self):
        """ Méthode permettant de générer les segments reliant les sommets.
        Si une opération round_corner est utilisé, cette opération est inutile."""
        if self.sides:
            logger.warning(
                "Warning : attribut sides d'une LineLoop écrasé "
                "lors de l'utilisation de la méthode vertices_2_sides."
            )
        self.sides = [
            Line(self.vertices[i - 1], self.vertices[i])
            for i in range(len(self.vertices))
        ]


class PlaneSurface(object):
    """
    Calque de la fonction Plane Surface native de gmsh
    Créée à partir d'une LineLoop définissant le contour extérieur
    et, si nécessaire, de line loops définissant des trous internes
    """

    def __init__(self, ext_contour, holes=[]):
        self.ext_contour = ext_contour
        self.holes = holes
        self.tag = None
        self.boundary = ext_contour.sides + [crv for h in holes for crv in h.sides]
        # Pour favoriser le duck typing ?

    def __eq__(self, other):
        """
         Opérateur de comparaison == surchargé pour les objets de la classe Plane Surface
         Orientation prise en compte. Considérer le cas de surfaces non orientées ????
        """

        if not isinstance(other, PlaneSurface):
            return False
        if self.ext_contour != other.ext_contour:
            return False
        if len(self.holes) != len(other.holes):
            return False
        if len(self.holes) != 0:
            for contour in self.holes:
                for unContourOther in other.holes:
                    if contour == unContourOther:
                        break
                else:
                    return False
        return True

    def add_gmsh(self):
        if self.tag:
            return self.tag
        all_loops = (
            [self.ext_contour] if not self.holes else [self.ext_contour] + self.holes
        )
        for ll in all_loops:
            if not ll.tag:
                ll.add_gmsh()
        self.tag = factory.addPlaneSurface([ll.tag for ll in all_loops])


class AbstractSurface(object):
    """
    Surface dont on ne connait rien à part le tag.
    Une surface existante dans le modèle gmsh peut être identifiée à l'aide de l'API
    puis représentée par une instance de AbstractSurface.
    Il s'agit par exemple du résulat d'une opération booléenne.
    Par contre, ses bords sont a priori inconnus.
    """

    def __init__(self, tag):
        self.tag = tag
        self.boundary = []

    def get_boundary(self, recursive=True):
        """
        Récupérer les tags des entitées géométriques 1D qui composent le bord.

        Parameters
        ----------
        recursive : bool, optional
            If True, the boundaries of the 1-D entities that form the boundary
            of the AbstractSurface instance are also extracted from the gmsh model.
            Instances of Point are created to represent them.

        """
        self.boundary = AbstractSurface.get_surfs_boundary(self, recursive=recursive)

    @staticmethod
    def get_surfs_boundary(surfs, recursive=True):
        """
        Get the tags of all the 1D geometry entities that form the boundary of a surface
        or a group of surfaces.

        Parameters
        ----------
        recursive : bool, optional
            If True, the boundaries of the 1D entities are also extracted
            from the gmsh model.
            Instances of Point are created to represent them.

        """
        def_crv = []
        try:
            for s in surfs:
                if not s.tag:
                    s.add_gmsh()
            dim_tags = [(2, s.tag) for s in surfs]
        except TypeError:
            if isinstance(surfs, (PlaneSurface, AbstractSurface)):
                if not surfs.tag:
                    surfs.add_gmsh()
                dim_tags = (2, surfs.tag)
            else:
Baptiste Durand's avatar
black  
Baptiste Durand committed
375
                raise TypeError
376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392
        boundary_ = model.getBoundary(
            dim_tags, combined=True, oriented=False, recursive=False
        )
        for dimtag in boundary_:
            if dimtag[0] != 1:
                logger.warning(
                    "Unexpected type of geometrical entity "
                    f"in the boundary of surfaces {dim_tags}"
                )
                continue
            new_crv = AbstractCurve(dimtag[1])
            if recursive:
                new_crv.get_boundary()
            def_crv.append(new_crv)
        return def_crv


Baptiste Durand's avatar
Baptiste Durand committed
393
def surface_bool_cut(body, tool):
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 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454
    """
    Boolean operation of cutting performed on surfaces.

    Remove the aeras taken by the tool entities from the body surface.
    Removing a set of geometrical entities 'tools' at once is possible.
    The removeObject and removeTool parameters of the gmsh API function
    are set to False in order to keep the consistency between
    python geometrical objects and the gmsh geometrical model as far as possible.

    Parameters
    ----------
    body : PlaneSurface or AbstractSurface
        Main operand of the cut operation.
    tool : instance of PlaneSurface/AbstractSurface or list of instances of them
        Several tool areas can be removed to the body surface at once.
        To do this, the tool parameter must be a list.

    Return
    ----------
    cut_surf : AbstractSurface
        Python object that represents the surface obtained with the boolean operation.
        This will be a degenerate instance with only a tag attribut
        and a boundary attribut that can be evaluate later.
    """

    if not body.tag:
        body.add_gmsh()
    if not tool:  # * =True if empty list
        logger.warning(
            "No entity in the tool list for boolean cut operation."
            "The 'body' surface is returned."
        )
        return [body]
    try:
        _ = (element for element in tool)
    except TypeError:
        logger.debug("tool convert to list for boolean operation.")
        tool = [tool]
    for t in tool:
        if not t.tag:
            t.add_gmsh()
    output = factory.cut(
        [(2, body.tag)],
        [(2, t.tag) for t in tool],
        removeObject=False,
        removeTool=False,
    )
    logger.debug(f"Output of boolean operation 'cut' on surfaces : {output}")
    new_surf = list()
    for entity in output[0]:
        if entity[0] == 2:
            new_surf.append(AbstractSurface(entity[1]))
        else:
            logger.warning(
                "Some outputs of a cut boolean operation are not surfaces and"
                "therefore are not returned."
                f"\n Complete output from the API function : {output}"
            )
    return new_surf


Baptiste Durand's avatar
Baptiste Durand committed
455
def surface_bool_intersect(body, tool):
456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477
    """
    Boolean operation of intersection performed on surfaces.

    See the bool_cut_S() doc for more informations.
    """
    if not body.tag:
        body.add_gmsh()
    if isinstance(tool, PlaneSurface):
        tool = [tool]
    assert isinstance(tool, list)
    for t in tool:
        if not t.tag:
            t.add_gmsh()
    ops_output = []
    for t in tool:
        outpt = factory.intersect(
            [(2, body.tag)], [(2, t.tag)], removeObject=False, removeTool=False
        )
        if outpt[0]:
            ops_output.append(outpt)
        else:  # Tool entirely outside of body or entirely inside.
            t_copy_dimtag = factory.copy([(2, t.tag)])
478
            factory.synchronize()  # ! Peut être supprimé ?
479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495
            outpt = factory.intersect(
                [(2, body.tag)], t_copy_dimtag, removeObject=False, removeTool=True
            )
            if outpt[0]:
                ops_output.append(outpt)
    new_surf = []
    for outpt in ops_output:
        if outpt[0][0][0] == 2:
            new_surf.append(AbstractSurface(outpt[0][0][1]))
        else:
            warn_msg = (
                "Some entities that result from a intersection boolean operation "
                "are not surfaces and therefore are not returned. \n"
                f"Complete output from the API function : \n {ops_output}"
            )
            logger.warning(warn_msg)
    return new_surf