Merge branch 'aux_triangle_microstructure'

.gitignore

@@ -9,4 +9,4 @@ demo/loc_E12_u.pdf
 *.pyc.*
 *.egg-info/
 ho_homog/.idea
-.idea/HO_homog.iml
+*.iml

ho_homog/geometry/curves.py

@@ -208,7 +208,7 @@ class AbstractCurve(Curve):
         """
         bndry_logger.debug(f"Abstract Curve -> get_boundary. self.tag : {self.tag}")
         def_pts = []
-        boundary = model.getBoundary((1, self.tag), combined=False)
+        boundary = model.getBoundary((1, self.tag), False, False, False)
         bndry_logger.debug(
             f"Abstract Curve -> get_boundary. raw API return : {boundary}"
         )

ho_homog/geometry/surfaces.py

@@ -10,6 +10,8 @@ Object designed to represent geometrical entitites of dimension two
 and instantiate them in a gmsh model.
 """
 
+from .curves import Line, Arc
+from .point import Point
 from . import factory, np, logger, model
 from .tools import round_corner, offset
 from .curves import Line, AbstractCurve
@@ -77,7 +79,8 @@ class LineLoop(object):
                         if test:
                             break
                     else:
-                        # Aucun break déclenché, i.e. si l'un des cote de la lineloop courante n'appartient pas à la LineLoop comparée #noqa
+                        # Aucun break déclenché,
+                        # i.e. si l'un des cote de la lineloop courante n'appartient pas à la LineLoop comparée
                         return False
                 else:
                     return True
@@ -316,7 +319,7 @@ class AbstractSurface(object):
                     surfs.add_gmsh()
                 dim_tags = (2, surfs.tag)
             else:
-                raise (TypeError)
+                raise TypeError
         boundary_ = model.getBoundary(
             dim_tags, combined=True, oriented=False, recursive=False
         )

ho_homog/geometry/tools.py

@@ -118,7 +118,15 @@ def round_corner(inp_pt, pt_amt, pt_avl, r, junction_raduis=False, plot=False):
     if alpha < 0:
         # corriger le cas des angles \in ]-pi;0[ = ]pi;2pi[]. L'arrondi est toujours dans le secteur <180° #noqa
         v_biss = -v_biss
-
+        
+   # if abs(alpha-math.pi)<10E-4:
+       # pt_amt_milieu=Point((pt_amt.coord+inp_pt.coord)/2)
+       # pt_avl_milieu=Point((pt_avl.coord+inp_pt.coord)/2)
+       # racc_amt = Line(pt_amt_milieu, inp_pt)
+        #racc_avl = Line(inp_pt, pt_avl_milieu)
+       # geoList = [racc_amt, racc_avl]
+        #raise ValueError("angle quasi plat")
+        #return geoList
     if junction_raduis:
         if plot:
             R = copy.deepcopy(r)
@@ -167,6 +175,33 @@ def round_corner(inp_pt, pt_amt, pt_avl, r, junction_raduis=False, plot=False):
     return geoList
 
 
+def fine_point_jonction(inp_pt, pt_amt, pt_avl, r,sharp_r):
+    """
+    Retourne le point au centre d'une jonction défini par deux angles aigus "dans le même sens"
+
+    """
+    from .transformations import translation
+
+    # Direction en amont et en aval
+    v_amt = unit_vect(pt_amt.coord - inp_pt.coord)
+    v_avl = unit_vect(pt_avl.coord - inp_pt.coord)
+
+    alpha = angle_between(v_amt, v_avl, orient=True)
+    v_biss = bisector(v_amt, v_avl)
+    if alpha < 0:
+        # corriger le cas des angles \in ]-pi;0[ = ]pi;2pi[]. L'arrondi est toujours dans le secteur <180° #noqa
+        v_biss = -v_biss
+        
+   
+
+    # Calcul des distances centre - sommet de l'angle et sommet - point de raccordement
+    dist_center = float(r/2+sharp_r) / abs(math.sin(alpha / 2.0))
+
+
+    pt_ctr = translation(inp_pt, dist_center * v_biss)
+
+    return pt_ctr
+
 def offset(pt, pt_dir1, pt_dir2, distance, method="vertex"):
     """[summary]
 

ho_homog/mesh_generate/__init__.py

@@ -31,6 +31,7 @@ logger = logging.getLogger(__name__)  # http://sametmax.com/ecrire-des-logs-en-p
 
 __all__ = [
     "pantograph",
+    "kagome",
     "duplicate_pattern",
     "offset_pattern",
     "Gmsh2DRVE",
@@ -351,6 +352,8 @@ class Gmsh2DPart(object):
 
 
 from .pantograph import pantograph_RVE, pantograph_offset_RVE, beam_pantograph_RVE
+from .kagome import kagome_RVE
+
 
 # from .other_2d_microstructures import auxetic_square_RVE
 
@@ -456,6 +459,10 @@ def Gmsh2DPartFromRVE(cell: Gmsh2DRVE, nb_cells, part_name=None):
     factory.synchronize()
     for gp in phy_surfaces:
         gp.add_gmsh()
+    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)
     # ! 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
@@ -497,10 +504,12 @@ def Gmsh2DPartFromRVE(cell: Gmsh2DRVE, nb_cells, part_name=None):
     return Gmsh2DPart(part_vect, nb_cells, phy_surfaces, part_path)
 
 
+from . import kagome
 from . import pantograph
 
 __all__ = [
     "pantograph_RVE",
+    "kagome_RVE",
     "pantograph_offset_RVE",
     "beam_pantograph_RVE",
     "pantograph_E11only_RVE",

ho_homog/mesh_generate/kagome.py

+#!/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
+from math import sin, cos, pi
+
+from . import Gmsh2DRVE, logger
+
+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
+)
+
+
+def kagome_RVE(alpha, r, a=None, b=None, nb_cells=(1, 1), offset=(0.0, 0.0), name=""):
+    """[summary]
+
+    Parameters
+    ----------
+    alpha : float
+        Paramètre d'ouverture, intervalle [0, 0.5]
+        0: configuration refermé
+        0.5: configuration complétement ouverte
+    name : str
+
+    r : float
+        junction thinness, rayon de jonction / côté d'un triangle
+
+    Returns
+    -------
+    Instance of the Gmsh2DRVE class.
+
+    La cellule est de taille constante,
+    par contre les triangles ont une taille qui dépend de l'angle d'ouverture.
+    """
+
+    logger.info("Start defining the geometry")
+    name = name if name else "kagome"
+
+    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)
+
+    gen_vect = np.array(((a, a / 2), (0.0, np.sqrt(3) / 2 * a)))
+    nb_cells, offset = np.asarray(nb_cells), np.asarray(offset)
+
+    a0 = np.array((a, 0.0, 0.0))
+    a1 = np.array((a / 2, np.sqrt(3) / 2 * a, 0.0))
+    phi1 = geo.Point(alpha * a0)
+    psi1 = geo.Point((1 - alpha) * a1)
+    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)
+
+    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)
+    logger.info("Removing duplicate pattern line-loops: Done")
+    logger.info(f"Number of pattern line-loops: {len(pattern_ll)}")
+    pattern_ll = [round_corner_kagome(ll, r * b, a, alpha) for ll in pattern_ll]
+    logger.info("Rounding all corners of pattern line-loops: Done")
+    fine_pts = [pt for ll in pattern_ll for pt in ll.vertices]
+    fine_pts = geo.remove_duplicates(fine_pts)
+
+    return Gmsh2DRVE(pattern_ll, gen_vect, nb_cells, offset, fine_pts, False, name,)
+
+
+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
+
+
+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

test/test_mesh_generate_kagome.py

+# coding: utf-8
+
+import gmsh
+import ho_homog.geometry as geo
+import ho_homog.mesh_generate as mg
+
+geo.init_geo_tools()
+geo.set_gmsh_option("Mesh.MshFileVersion", 4.1)
+
+
+def test_kagome():
+    """ Generate a mesh that corresponds to one unit cell of 'kagome' microstructure.
+    Two cases : alpha = 0.5, fully open, alpha = 0.1 triangles almost close"""
+    a = 1
+    junction_thinness = 1e-2
+    for alpha in [0.1, 0.5]:
+        rve = mg.kagome.kagome_RVE(
+            0.5, junction_thinness, a=a, name=f"kagome_{alpha:.2f}"
+        )
+        lc_ratio = 1 / 2
+        lc_min_max = (lc_ratio * junction_thinness * a, lc_ratio * a)
+        d_min_max = (2 * junction_thinness * a, a)
+        rve.main_mesh_refinement(d_min_max, lc_min_max, False)
+        rve.mesh_generate()
+        gmsh.model.mesh.renumberNodes()
+        gmsh.model.mesh.renumberElements()
+        gmsh.write(str(rve.mesh_abs_path))