kagome.py 6.77 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Nov  9 10:12:32 2020

@author: manon_thbaut
"""
import numpy as np
import logging
import ho_homog.geometry as geo
import gmsh
import ho_homog.mesh_tools as msh
13
from math import sin, cos, pi
14 15 16

from . import Gmsh2DRVE, logger

17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
model = gmsh.model
factory = model.occ

# * Logging
logger = logging.getLogger(__name__)  # http://sametmax.com/ecrire-des-logs-en-python/
logger.setLevel(logging.INFO)


formatter = logging.Formatter("%(asctime)s :: %(levelname)s :: %(message)s", "%H:%M")
stream_handler = logging.StreamHandler()
stream_handler.setLevel(logging.INFO)
stream_handler.setFormatter(formatter)
logger.addHandler(stream_handler)

xdmf_parameters = dict(
    flush_output=False, functions_share_mesh=True, rewrite_function_mesh=False
)


Baptiste Durand's avatar
Baptiste Durand committed
36
def kagome_RVE(alpha, r, a=None, b=None, nb_cells=(1, 1), offset=(0.0, 0.0), name=""):
37 38 39 40 41 42 43 44 45
    """[summary]

    Parameters
    ----------
    alpha : float
        Paramètre d'ouverture, intervalle [0, 0.5]
        0: configuration refermé
        0.5: configuration complétement ouverte
    name : str
Baptiste Durand's avatar
Baptiste Durand committed
46

47 48 49
    r : float
        junction thinness, rayon de jonction / côté d'un triangle

50 51 52 53
    Returns
    -------
    Instance of the Gmsh2DRVE class.

54 55 56
    La cellule est de taille constante,
    par contre les triangles ont une taille qui dépend de l'angle d'ouverture.
    """
57

58
    logger.info("Start defining the geometry")
59 60
    name = name if name else "kagome"

Baptiste Durand's avatar
Baptiste Durand committed
61 62 63 64 65 66 67 68 69
    if (not a) and (not b):
        raise ValueError("a or b (exclusive OR) must be imposed")
    elif a and b:
        raise ValueError("a and b cannot be imposed simultaneously")
    elif b:
        # si on choisit la taille des triangles b
        a = kagome_triangle_size_2_cell_size(alpha, b)
    elif a:
        b = kagome_cell_size_2_triangle_size(alpha, a)
70

Baptiste Durand's avatar
Baptiste Durand committed
71 72
    gen_vect = np.array(((a, a / 2), (0.0, np.sqrt(3) / 2 * a)))
    nb_cells, offset = np.asarray(nb_cells), np.asarray(offset)
73 74

    a0 = np.array((a, 0.0, 0.0))
75
    a1 = np.array((a / 2, np.sqrt(3) / 2 * a, 0.0))
76 77
    phi1 = geo.Point(alpha * a0)
    psi1 = geo.Point((1 - alpha) * a1)
Baptiste Durand's avatar
Baptiste Durand committed
78 79 80 81 82 83
    z_axis = [0, 0, 1]
    angle = 2 * np.pi / 3
    phi2 = geo.rotation(phi1, angle, z_axis)
    psi2 = geo.rotation(psi1, angle, z_axis)
    phi3 = geo.rotation(phi2, angle, z_axis)
    psi3 = geo.rotation(psi2, angle, z_axis)
84

85 86 87 88 89 90
    star_pts = [phi1, psi1, phi2, psi2, phi3, psi3]
    pattern_ll = [geo.LineLoop(star_pts, explicit=False)]
    trans = [a0, a1, a0 + a1]
    pattern_ll += [geo.translation(pattern_ll[0], t) for t in trans]

    pattern_ll = geo.remove_duplicates(pattern_ll)
Baptiste Durand's avatar
Baptiste Durand committed
91 92
    logger.info("Removing duplicate pattern line-loops: Done")
    logger.info(f"Number of pattern line-loops: {len(pattern_ll)}")
93
    pattern_ll = [round_corner_kagome(ll, r * b, a, alpha) for ll in pattern_ll]
Baptiste Durand's avatar
Baptiste Durand committed
94
    logger.info("Rounding all corners of pattern line-loops: Done")
95 96
    fine_pts = [pt for ll in pattern_ll for pt in ll.vertices]
    fine_pts = geo.remove_duplicates(fine_pts)
97

Baptiste Durand's avatar
Baptiste Durand committed
98
    return Gmsh2DRVE(pattern_ll, gen_vect, nb_cells, offset, fine_pts, False, name,)
99 100


Baptiste Durand's avatar
Baptiste Durand committed
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
def kagome_triangle_size_2_cell_size(alpha, b):
    """Passage de b (côté des triangles) à a (taille caractéristique cellule unitaire du kagome)."""
    _a = 1.0  # pour cellule de hauteur 1
    x = (1 - alpha) * _a / 2 - alpha * _a
    y = (1 - alpha) * _a * np.sqrt(3) / 2
    _b = np.sqrt(x ** 2 + y ** 2)
    theta = np.arcsin(np.sqrt(3) * alpha / (2 * _b))
    # compute theta for a unit cell

    a = 2 * b * np.cos(np.pi / 3 - theta)  # real cell
    return a


def kagome_cell_size_2_triangle_size(alpha, a):
    """Pour la microstructure "kagome", passage de a (taille caractéristique de la cellule unitaire) à  b (côté des triangles)

    Parameters
    ----------
    alpha : float
        paramètre d'ouverture de la microstructure
    a : float
        taille cellule unitaire parallélogramme.
    """
    t1 = (1 - alpha) * a / 2 - alpha * a
    t2 = (1 - alpha) * np.sqrt(3) * a / 2
    b = np.sqrt(t1 ** 2 + t2 ** 2)
    return b
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 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 208 209 210 211 212 213 214 215 216 217 218


def round_corner_kagome(lineloop, r, a, alpha):
    """ Opération d'arrondi des angles spécifique à la microstructure 'kagome',
    appliquée à tous les sommets du polygone.

    Parameters
    ----------
    lineloop : LineLoop
        Contour à modifier
    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
    """

    effect_r, phi_1, phi_2 = calcul_effective_r(alpha, r, a)
    vertices = lineloop.vertices
    # ! ESSAI
    results_1d = list()
    for i in range(len(vertices)):
        cur_pt = vertices[i - 1]
        d2 = effect_r / np.sin(phi_2)
        d1 = effect_r / np.sin(phi_1)

        dir_1 = vertices[i - 2].coord - cur_pt.coord
        dir_2 = vertices[i].coord - cur_pt.coord
        dir_1 = geo.unit_vect(dir_1)
        dir_2 = geo.unit_vect(dir_2)

        pt_amt = geo.translation(cur_pt, effect_r * dir_1)
        pt_avl = geo.translation(cur_pt, effect_r * dir_2)

        alpha = geo.angle_between(dir_1, dir_2, orient=True)
        v_biss = geo.bisector(dir_1, dir_2)
        if alpha < 0:
            v_biss = -v_biss

        if abs(abs(geo.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(geo.angle_between(v_biss, dir_1)) - (np.pi / 2 - phi_1)) < 10e-14:
            d = d1
        else:
            raise ValueError("mauvaise gestion de d1 et d2")
        center = geo.translation(cur_pt, d * v_biss)
        round_arc = geo.Arc(pt_amt, center, pt_avl)
        racc_amt = geo.Line(vertices[i - 2], pt_amt)
        racc_avl = geo.Line(pt_avl, vertices[i])
        curves_list = [racc_amt, round_arc, racc_avl]
        results_1d.append(curves_list)
    lineloop.sides = geo.surfaces.round_corner_2_sides(results_1d)
    return lineloop


def calcul_effective_r(alpha, r, a):
    """
    Méthode de construction des jonctions propre au kagomé.

    Parameters
    ----------
    alpha: float
        paramètre d'ouverture
    r: float
        rayon des jonctions
    a: float
        taille caractéristique de la cellule unitaire

    Returns
    -------
    tuple
        3 floats: Effective radius, phi_1, phi_2
    """

    b = kagome_cell_size_2_triangle_size(alpha, a)
    theta = np.arcsin(np.sqrt(3) * alpha * a / (2 * b))

    phi_2 = np.pi / 2 - theta
    if alpha < 1 / 3:
        phi_1 = np.pi / 6 - theta
        effect_r = (2 * r * sin(phi_1) * sin(phi_2)) / (
            sin(phi_2) * cos(phi_1) - sin(phi_1) * cos(phi_2) - sin(phi_2) + sin(phi_1)
        )
    else:
        phi_1 = theta - pi / 6
        effect_r = (2 * r * sin(phi_1) * sin(phi_2)) / (
            -sin(phi_2) * cos(phi_1) - sin(phi_1) * cos(phi_2) + sin(phi_2) + sin(phi_1)
        )
    return effect_r, phi_1, phi_2