Commit 30837335 authored by Baptiste Durand's avatar Baptiste Durand

Split geometry into several modules

parent feb75de4
This diff is collapsed.
# coding: utf8
"""
Created on 09/10/2018
@author: baptiste
Definition of the classes : Curve, Line, Arc and AbstractCurve.
Objects designed to represent geometrical entitites of dimension one
and instantiate them in a gmsh model.
"""
import logging
import warnings
from logging.handlers import RotatingFileHandler
from pathlib import Path
from . import factory, logger, model, np, plt
from .point import Point
bndry_logger = logging.getLogger("bndry")
bndry_logger.setLevel(logging.DEBUG)
bndry_logger.propagate = False
log_path = Path("~/ho_homog_log/gmsh_getBoundary_output.log").expanduser()
if not log_path.parent.exists():
log_path.parent.mkdir(mode=0o777, parents=True)
file_handler = RotatingFileHandler(str(log_path), "a", 1000000, 1)
file_handler.setLevel(logging.DEBUG)
file_handler.setFormatter(
logging.Formatter("%(asctime)s :: %(levelname)s :: %(message)s")
)
bndry_logger.addHandler(file_handler)
bndry_logger.warning("***Output of the getBoundary function of the gmsh API***")
bndry_logger.warning(
"""
*****
This logger and the associated handler are a hacky workaround to avoid errors
related to gmsh.model.occ.getBoundary(). Therefore they should not be removed
or deactivated.
*****"""
)
class Curve(object):
"""Superclass that is used to define the other classes : Line, Arc and AbstractCurve.
It is designed to represent geometrical entities of dimension one.
"""
def __init__(self, def_pts_list, gmsh_api_add_function):
self.def_pts = def_pts_list
self.tag = None
self.gmsh_constructor = gmsh_api_add_function
Curve.all_instances.append(self)
def __eq__(self, other):
"""
Return True if and only if :
- both self and other are instances of the same subclass,
AND
- The coordinates of the points that are used to define
these two Lines (or Arcs) are equal.
"""
if not type(other) is type(self):
return False
return all(p == q for p, q in zip(self.def_pts, other.def_pts))
def __ne__(self, other):
return not self.__eq__(other)
def add_gmsh(self):
"""Instantiate a geometrical entity in the current gmsh model"""
if self.tag:
return None
for pt in self.def_pts:
if not pt.tag:
pt.add_gmsh()
self.tag = self.gmsh_constructor(*[p.tag for p in self.def_pts])
return self.tag
# IDEA
# ? @property
# ? def tag(self):
# ? if self.__tag is None:
# ? self.add_gmsh()
# ? return self.__tag
def reverse(self):
""" Peut-être utile au moment de définir les contraites de périodicité.
# TODO À tester !!!
"""
self.def_pts.reverse()
if self.tag:
self.tag *= -1
class Line(Curve):
"""A line is defined by 2 instances of Point (start + end)"""
def __init__(self, start_pt, end_pt):
Curve.__init__(self, [start_pt, end_pt], factory.addLine)
def __str__(self):
"""Affichage plus clair des coordonnées des points de départ et d'arrivée."""
prt_str = (
f"Line {self.tag if self.tag else '--'}, "
f"Point tags : start {self.def_pts[0].tag}, end {self.def_pts[1].tag}"
)
return prt_str
def length(self):
return np.linalg.norm(self.def_pts[0].coord - self.def_pts[1].coord)
def direction(self):
""" Renvoie un vecteur unitaire correspondant à la direction de la ligne"""
return 1 / self.length() * (self.def_pts[1].coord - self.def_pts[0].coord)
def plot2D(self, color="black"):
"""En 2D seulement. Tracer la ligne dans un plot matplotlib. """
x = [pt.coord[0] for pt in self.def_pts]
y = [pt.coord[1] for pt in self.def_pts]
plt.plot(x, y, color=color)
class Arc(Curve):
"""Classe définissant un arc de cercle caractérisé par :
- son point de départ
- son point d'arrivée
- son centre
- son rayon
"""
def __init__(self, start_pt, center_pt, end_pt):
""" Crée un arc de cercle.
Vérification préalable de l'égalité des distances centre<->extrémités.
"""
d1 = np.linalg.norm(start_pt.coord - center_pt.coord)
d2 = np.linalg.norm(end_pt.coord - center_pt.coord)
np.testing.assert_almost_equal(d1, d2, decimal=10)
Curve.__init__(self, [start_pt, center_pt, end_pt], factory.addCircleArc)
self.radius = (d1 + d2) / 2
def __str__(self):
prt_str = (
f"Circle arc {self.tag if self.tag else '--'}, "
f"Point tags : start {self.def_pts[0].tag}, center {self.def_pts[1].tag}, "
f"end {self.def_pts[2].tag}"
)
return prt_str
def plot2D(
self,
circle_color="Green",
end_pts_color="Blue",
center_color="Orange",
pt_size=5,
):
"""Représenter l'arc de cercle dans un plot matplotlib.
Disponible seulement en 2D pour l'instant."""
self.def_pts[0].plot(end_pts_color, pt_size)
self.def_pts[2].plot(end_pts_color, pt_size)
self.def_pts[1].plot(center_color, pt_size)
circle = plt.Circle(
(self.def_pts[1].coord[0], self.def_pts[1].coord[1]),
self.radius,
color=circle_color,
fill=False,
)
ax = plt.gca()
ax.add_patch(circle)
class AbstractCurve(Curve):
# TODOC
@staticmethod
def empty_constructor(*args):
warnings.warn(
"Adding an AbstractCurve instance to the gmsh model is not supported. "
"\n These python objects actually represent unknown geometrical entities"
" that already exist in the gmsh model.",
UserWarning,
)
def __init__(self, tag):
"""
Créer une représentation d'une courbe existant dans le modèle Gmsh.
# ! A corriger
# ! Lors de l'instantiation, les points extrémités peuvent être donnés,
# ! soit explicitement, soit par leurs tags.
"""
Curve.__init__(self, [], AbstractCurve.empty_constructor)
self.tag = tag
def get_boundary(self, get_coords=True):
"""
Récupérer les points correspondants aux extrémités de la courbe dans le modèle Gmsh.
Parameters
----------
coords : bool, optional
If true, les coordonnées des points extrémités sont aussi récupérés.
"""
bndry_logger.debug(f"Abstract Curve -> get_boundary. self.tag : {self.tag}")
def_pts = []
boundary = model.getBoundary((1, self.tag), False, False, False)
bndry_logger.debug(
f"Abstract Curve -> get_boundary. raw API return : {boundary}"
)
for pt_dimtag in boundary:
if not pt_dimtag[0] == 0:
raise TypeError(
f"The boundary of the geometrical entity {self.tag} are not points."
)
coords = model.getValue(0, pt_dimtag[1], []) if get_coords else []
logger.debug(repr(coords))
new_pt = Point(np.array(coords))
new_pt.tag = pt_dimtag[1]
def_pts.append(new_pt)
self.def_pts = def_pts
def plot2D(self, color="black"):
"""En 2D seulement. Tracé des points de def_pts, reliés en pointillés"""
x = [pt.coord[0] for pt in self.def_pts]
y = [pt.coord[1] for pt in self.def_pts]
plt.plot(x, y, color=color, linestyle="dashed")
# coding: utf8
"""
Created on 09/10/2018
@author: baptiste
Definition of a class designed to represent gmsh physical entities.
"""
from . import logger, factory, model, set_gmsh_option
class PhysicalGroup(object):
"""
Create and manage groups of geometrical entities in the gmsh model.
Physical groups can be visible in the exported mesh.
"Groups of elementary geometrical entities can also be defined
and are called “physical” groups.
These physical groups cannot be modified by geometry commands: their only
purpose is to assemble elementary entities into larger groups so that they can be
referred to by the mesh module as single entities." (gmsh reference manual 26/09/2018)
"""
all_groups = dict()
def __init__(self, entities, geo_dim, name=None):
"""
Gather multiple instances of one of the geometrical Python classes
(Point, Curve, PlaneSurface) to form a single object in the gmsh model.
Parameters
----------
entities : list
The instances that will compose the physical group.
They must have the same geometrical dimension
0 for points,
1 for line, arcs, instances of AbstractCurve,
2 for surfaces,
3 for volumes
geo_dim : int
Geometrical dimension of the entities that are gathered.
name : string, optional
name of the group.
"""
try:
_ = (item for item in entities)
except TypeError:
logger.debug("Physical group single entity -> Input conversion to list")
entities = [entities]
self.entities = entities
self.dim = geo_dim
self.name = name
self.tag = None
try:
PhysicalGroup.all_groups[self.dim].append(self)
except KeyError:
PhysicalGroup.all_groups[self.dim] = [self]
def add_gmsh(self):
factory.synchronize()
if self.tag:
return self.tag
tags = list()
for item in self.entities:
if not item.tag:
item.add_gmsh()
tags.append(item.tag)
factory.synchronize()
self.tag = model.addPhysicalGroup(self.dim, tags)
# ! TEMPORAIRE, un appel à synchronize devrait pouvoir être enlevé.
logger.info(f"Physical group {self.tag} of dim {self.dim} added to gmsh")
phy_before = model.getPhysicalGroups()
factory.synchronize()
phy_after = model.getPhysicalGroups()
dbg_msg = (
"Call of factory.synchronize(). Physical groups in the model \n"
f"before : \n {phy_before} \n "
f"and after : \n {phy_after}"
)
logger.debug(dbg_msg)
factory.synchronize()
logger.debug(
f"And after a 2nd call to synchronize() : {model.getPhysicalGroups()}"
)
if self.name:
model.setPhysicalName(self.dim, self.tag, self.name)
def add_to_group(self, entities):
"""
Add geometrical entities to an existing physical group.
entitites :
A geometrical entity or a list of geometrical entities.
The appended items must be of the same geometrical dimension.
"""
if self.tag:
raise AttributeError(
"The physical group has been already defined in the gmsh model."
"It is too late to add entities to this group."
)
if isinstance(entities, list):
self.entities += entities
else:
self.entities.append(entities)
def set_color(self, rgba_color, recursive=False):
"""
Choisir la couleur du maillage des éléments géométriques de l'entité physique.
Parameters
----------
rgba_color : list of 4 integers between 0 and 255.
RGBA code of the desired color.
recursive : bool, optional
Apply the color setting to the parent geometrical entities as well.
"""
dimtags = [(self.dim, e.tag) for e in self.entities]
model.setVisibility(dimtags, 1)
model.setColor(dimtags, *rgba_color, recursive=recursive)
def set_visibility(self, val: bool, recursive: bool = False):
"""Show or hide entities that belong to the physical group."""
dim_tags = [(self.dim, entity.tag) for entity in self.entities]
visibility = 1 if bool(val) else 0
model.setVisibility(dim_tags, visibility, bool(recursive))
@classmethod
def set_group_visibility(cls, val):
"""
Make only entities that belong to at least one physical group visible,
or make all geometrical entities visibles.
Parameters
----------
val : bool
If True, only entities that to at least one physical group visible.
If False, all geometrical entities will be visible.
"""
if val:
model.setVisibility(model.getEntities(), 0)
dimtags = list()
for gps in cls.all_groups.values():
for gp in gps:
dimtags += [(gp.dim, ent.tag) for ent in gp.entities]
model.setVisibility(dimtags, 1, recursive=True)
else:
model.setVisibility(model.getEntities(), 1)
return None
@classmethod
def set_group_mesh(cls, val):
"""
Mesh only the entities that belong to at least one physical group.
Based on the EXPERIMENTAL MeshOnlyVisible option.
Parameters
----------
val : bool
If True, only entities that to at least one physical group will be mesh.
If False, a mesh will be generate for all geometrical entities.
"""
if val:
cls.set_group_visibility(1)
set_gmsh_option("Mesh.MeshOnlyVisible", 1)
else:
cls.set_group_visibility(0)
set_gmsh_option("Mesh.MeshOnlyVisible", 0)
# coding: utf8
"""
Created on 09/10/2018
@author: baptiste
Definition of the Point class.
Base geometrical object for model definition in gmsh.
"""
from . import factory, plt, np
class Point(object):
# TODOC
def __init__(self, coord=np.array((0.0, 0.0)), mesh_size: float = 0.0):
"""Create a new point.
Parameters
----------
coord : 1D ndarray (vector), optional
coordinates of the point, in 2D or 3D.
By default : np.array((0., 0.))
mesh_size : float, optional
Characteristic length imposed for the mesh generation at this point.
By default : 0
"""
coord = np.asarray(coord) # * existing arrays are not copied
dim = coord.shape[0]
self.coord = coord
if dim == 2:
self.coord = np.append(self.coord, [0.0])
# ? Choix : on utilise toujours des coordonnés en 3D.
self.tag = None
self.mesh_size = mesh_size
def __repr__(self):
"""Represent a Point object with the coordinates of this point."""
return f"Pt {self.tag} ({str(self.coord)}) "
def __eq__(self, other):
"""Opérateur de comparaison == redéfini pour les objets de la classe Point.
Renvoie True ssi les coordonnées sont égales,
sans prendre en compte la signature de l'objet.
"""
if not isinstance(other, Point):
return False
if np.array_equal(self.coord, other.coord):
return True
elif np.allclose(self.coord, other.coord):
return True
else:
return False
def __ne__(self, other):
return not self.__eq__(other)
# NOTE Adapter à la 3D ?
def plot2D(self, color="red", size=5):
plt.plot(self.coord[0], self.coord[1], marker="o", markersize=size, color=color)
def add_gmsh(self):
""" Add the point to the current gmsh model via the python API."""
if self.tag:
# CAVEAT : True if the geometrical entity has already been instantiated.
return None # for information purposes only.
self.tag = factory.addPoint(*self.coord, self.mesh_size)
return self.tag
# 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.
"""
from . import factory, np, logger, model
from .tools import round_corner, offset
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:
# Aucun break déclenché, i.e. si l'un des cote de la lineloop courante n'appartient pas à la LineLoop comparée #noqa
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:
elmt.plot(color)
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()
def offset(self, t):
""" 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(
self.vertices[i - 1], self.vertices[i - 2], self.vertices[i], t
)
self.offset_dpcmt = [