"""
NRV-:class:`.MshCreator` handling.
"""
import math
import os
import gmsh
import numpy as np
from ....backend._file_handler import rmv_ext
from ....backend._log_interface import pass_info, rise_error, rise_warning
from ....backend._NRV_Class import NRV_class
from ....backend._parameters import parameters
dir_path = parameters.nrv_path + "/_misc"
###############
## Constants ##
###############
pi = np.pi
[docs]
def is_MshCreator(object):
"""
Check if an object is a MshCreator, return True if yes, else False
Parameters
----------
result : object
object to test
Returns
-------
bool
True it the type is a MshCreator object
"""
return isinstance(object, MshCreator)
[docs]
def clear_gmsh():
if gmsh.isInitialized():
gmsh.finalize()
[docs]
class MshCreator(NRV_class):
"""
Class handling the creation of a gmsh mesh (https://gmsh.info/doc/texinfo/gmsh.html)
Contains methodes dealing with the mesh geometries, physical domains and feilds
"""
[docs]
def __init__(self, D=3, ver_level=None):
"""
initialisation of the MshCreator
Parameters
----------
D : int(1,2,3)
mesh dimension
ver_level : int(0,1,2,3,4,5,99)
verbosity level of gmsh (see MshCreator.set_verbosity), by default 2
"""
super().__init__()
self.type = "GMSH"
self.D = D
self.entities = {}
self.volumes = []
self.volumes_com = []
self.volumes_bd = []
self.faces = []
self.faces_com = []
self.faces_bd = []
self.id_domains = np.array([], dtype=int)
self.dim_domains = np.array([], dtype=int)
self.name_domain = []
self.file = ""
self.Nfeild = 0
self.res = 1
clear_gmsh()
gmsh.initialize()
self.verbosity_level = ver_level
self.set_verbosity(ver_level)
self.model = gmsh.model
self.model.add("DFG 3D")
# Mesh properties
self.is_generated = False
self.N_entities = 0
self.N_nodes = 0
self.N_elements = 0
self.Ncore = parameters.GMSH_Ncores
gmsh.option.setNumber("General.NumThreads", self.Ncore)
if self.Ncore > 1:
gmsh.option.set_number("Mesh.Algorithm3D", 10)
else:
gmsh.option.set_number("Mesh.Algorithm3D", 1)
#####################
## special methods ##
#####################
[docs]
def get_obj(self):
"""
update and return list of mesh entities
Returns
-------
self.entities :dict
"""
return self.entities
[docs]
def get_volumes(self, com=False, bd=False):
"""
update and return list of mesh volumes (optional: with their center of mass)
Returns
-------
self.faces :list[tuple]
"""
self.model.occ.synchronize()
self.volumes = self.model.getEntities(dim=3)
self.volumes_com = []
self.volumes_bd = []
for i in range(len(self.volumes)):
volume = self.volumes[i]
self.volumes_com += [
np.round(self.model.occ.getCenterOfMass(volume[0], volume[1]), 4)
]
self.volumes_bd += [
np.round(self.model.occ.getBoundingBox(volume[0], volume[1]), 4)
]
if com:
if bd:
return self.volumes, self.volumes_com, self.volumes_bd
else:
return self.volumes, self.volumes_com
else:
if bd:
return self.volumes, self.volumes_bd
return self.volumes
[docs]
def get_info(self, verbose=False):
entities = self.model.getEntities()
nodeTags = self.model.mesh.getNodes()[0]
elemTags = self.model.mesh.getElements()[1]
self.N_entities = len(entities)
self.N_nodes = len(nodeTags)
self.N_elements = sum(len(i) for i in elemTags)
if verbose:
pass_info("Mesh properties:")
pass_info("Number of processes : " + str(self.Ncore))
pass_info("Number of entities : " + str(self.N_entities))
pass_info("Number of nodes : " + str(self.N_nodes))
pass_info("Number of elements : " + str(self.N_elements))
return self.Ncore, self.N_entities, self.N_nodes, self.N_elements
[docs]
def get_mesh_info(self, verbose=False):
rise_warning("DEPRECATED method use get_info instead")
self.get_info(verbose=verbose)
[docs]
def get_faces(self, com=False, bd=False):
"""
update and return list of mesh face (optional: with their center of mass)
Returns
-------
self.faces :list[tuple]
"""
self.model.occ.synchronize()
self.faces = self.model.getEntities(dim=2)
self.faces_com = []
self.faces_bd = []
for i in range(len(self.faces)):
face = self.faces[i]
self.faces_com += [
np.round(self.model.occ.getCenterOfMass(face[0], face[1]), 4)
]
self.faces_bd += [
np.round(self.model.occ.getBoundingBox(face[0], face[1]), 4)
]
if com:
if bd:
return self.faces, self.faces_com, self.faces_bd
else:
return self.faces, self.faces_com
else:
if bd:
return self.faces, self.faces_com
return self.faces
[docs]
def get_res(self):
"""
return the global resolution saved (usefull when no field are set)
Returns
-------
res :float
global resolution saved in the object
"""
return self.res
[docs]
def set_res(self, new_res):
"""
set the global resolution saved (usefull when no field are set)
Parameters
----------
new_res :float
global resolution to set the object
"""
self.res = new_res
[docs]
def set_verbosity(self, i=None):
"""
from gmsh: Level of information printed on the terminal and the message console.
Parameters
----------
i : int (1, 2, 3, 4, 5 or 99)
- 0, silent except for fatal errors
- 1, +errors
- 2, +warnings
- 3, +direct
- 4, +information
- 5, +status
- 99, +debug
"""
if i is None:
i = parameters.VERBOSITY_LEVEL
self.verbosity_level = i
gmsh.option.setNumber("General.Verbosity", self.verbosity_level)
[docs]
def set_chara_blen(self, i=0):
"""
from gmsh: Extend characteristic lengths from the boundaries inside the surface/volume
Parameters
----------
i : int, float, bool
Parameter value, by default 0
"""
i = float(i)
gmsh.option.set_number("Mesh.CharacteristicLengthExtendFromBoundary", i)
##############################################################################################
####################################### geometry methods ##################################
##############################################################################################
[docs]
def add_point(self, x=0, y=0, z=0):
"""
add a point to the mesh
Parameters
----------
x : float
x position of the first face center
y : float
y position of the first face center
z : float
z position of the first face center
Returns
-------
point : int
id of the added object
"""
point = self.model.occ.addPoint(x, y, z)
self.model.occ.synchronize()
return point
[docs]
def add_line(self, X0, X1):
"""
add a line to the mesh
Parameters
----------
X0 : int or tupple(3)
x position of the first face center
X1 : int or tupple(3)
y position of the first face center
Returns
-------
line : int
id of the added object
"""
if isinstance(X0, int):
ix0 = X0
elif np.iterable(X0):
if len(X0) != self.D:
rise_error("not enough dimension given in add_line")
else:
ix0 = self.add_point(X0[0], X0[1], X0[2])
if isinstance(X1, int):
ix0 = X0
elif np.iterable(X1):
if len(X1) != self.D:
rise_error("not enough dimension given in add_line")
else:
ix1 = self.add_point(X1[0], X1[1], X1[2])
line = self.model.occ.addLine(ix0, ix1)
self.model.occ.synchronize()
return line
[docs]
def add_box(self, x=0, y=0, z=0, ax=5, ay=1, az=1):
"""
add a box to the mesh
Parameters
----------
x : float
x position of the first face center
y : float
y position of the first face center
z : float
z position of the first face center
ax : float
Box length along x
ay : float
Box length along y
az : float
Box length along z
Returns
-------
box : int
id of the added object
"""
if self.D == 3:
parameters = {"x": x, "y": y, "z": z, "ax": ax, "ay": ay, "az": az}
box = self.model.occ.addBox(x, y, z, ax, ay, az)
self.model.occ.synchronize()
bounds = self.model.getEntities(dim=2)[-6:]
self.entities[box] = {
"type": "box",
"parameters": parameters,
"bounds": bounds,
"dim": 3,
}
return box
else:
rise_warning("Not added : add_cylinder requiere 3D mesh")
return None
[docs]
def add_cylinder(self, x=0, y=0, z=0, L=5, R=1):
"""
add a x-oriented cylinder to mesh entities
Parameters
----------
x : float
x position of the x-min face center
y : float
y position of the x-min face center
z : float
z position of the x-min face center
L : float
Cylinder length along x
R : float
Cylinder radius
Returns
-------
cyl : int
id of the added object
"""
if self.D == 3:
parameters = {"x": x, "y": y, "z": z, "L": L, "R": R}
cyl = self.model.occ.addCylinder(x, y, z, L, 0, 0, R)
self.model.occ.synchronize()
bounds = self.model.getEntities(dim=2)[-3:]
self.entities[cyl] = {
"type": "cylinder",
"parameters": parameters,
"bounds": bounds,
"dim": 3,
}
return cyl
else:
rise_warning("Not added : add_cylinder requiere 3D mesh")
return None
[docs]
def add_cone(self, x=0, y=0, z=0, L=5, R1=1, R2=0):
"""
add a x-oriented cone to mesh entities
Parameters
----------
x : float
x position of the x-min face center.
y : float
y position of the x-min face center.
z : float
z position of the x-min face center.
L : float
Cone length along x.
R1 : float
Cone x-min face radius.
R2 : float
Cone x-max face radius.
Returns
-------
cone : int
id of the added object
"""
if self.D == 3:
parameters = {"x": x, "y": y, "z": z, "L": L, "R1": R1, "R2": R2}
cone = self.model.occ.addCone(x, y, z, L, 0, 0, R1, R2)
self.model.occ.synchronize()
bounds = self.model.getEntities(dim=2)[-3:]
self.entities[cone] = {
"type": "cone",
"parameters": parameters,
"bounds": bounds,
"dim": 3,
}
return cone
else:
rise_warning("Not added : add_cylinder requiere 3D mesh")
return None
[docs]
def rotate(self, volume, angle, x=0, y=0, z=0, ax=0, ay=0, az=0, rad=True):
"""
rotate volume
Parameters
----------
volume : int
gmsh id of the volume
angle: : float
angle of the rotation
x : float
x-position of the center of the rotation, by default 0
y : float
y-position of the center of the rotation, by default 0
z : float
z-position of the center of the rotation, by default 0
ax : float
x-coefficient of the rotation axis direction vector, by default 0
ay : float
y-coefficient of the rotation axis direction vector, by default 0
az : float
z-coefficient of the rotation axis direction vector, by default 0
rad : bool
if true angle considered in rad, else in degree, by default 0
"""
if not rad:
angle = math.radians(angle)
self.model.occ.rotate([(3, volume)], x, y, z, ax, ay, az, angle)
[docs]
def fragment(self, IDs=None, dim=3, verbose=True):
"""
Fragmentation of the mesh important to link entities to each other
Parameters
----------
IDs : list or None
list of IDs of the element to fragments, if None all dim dimension elements are
fragmented, by default None
dim : int
dimension of the elements considerated, by default 3
verbose : bool
print or not the verbose on the temrminal, by default False
"""
if IDs is None:
if dim == 2:
list_obj = self.get_faces(com=False)
elif dim == 3:
list_obj = self.get_volumes(com=False)
else:
list_obj = [(dim, k) for k in self.entities]
elif not np.iterable(IDs) or len(IDs) < 2:
rise_warning("Need at least 2 entities to fragment")
return -1
else:
list_obj = [(dim, k) for k in IDs]
frag = self.model.occ.fragment([list_obj[0]], [k for k in list_obj[1:]])
new_entities = {}
if verbose:
pass_info(
"Warning: New volume generated by fragmentation, bounds are no longer up to date"
)
for i in frag[0][:]:
mask = [i in k for k in frag[1]]
p_id = (np.array(list_obj)[mask][:, 1]).tolist()
p_types = [self.entities[k]["type"] for k in p_id]
dim = min([self.entities[k]["dim"] for k in p_id])
com = np.round(self.model.occ.getCenterOfMass(i[0], i[1]), 3)
parameters = {"p_id": p_id, "p_types": p_types, "com": com}
new_entities[i[1]] = {
"type": "fragment",
"parameters": parameters,
"dim": dim,
}
self.entities.update(new_entities)
##############################################################################################
####################################### domains methods ###################################
##############################################################################################
[docs]
def add_domains(self, obj_IDs, phys_ID, dim=None, name=None):
"""
add domains (ID + name) to a goupe of entities
Caution: as to be used after all entities are placed
Parameters
----------
fname : str
path and name of saving file. If ends with ".msh" only save in ".msh" file
"""
if not np.iterable(obj_IDs):
obj_IDs = [obj_IDs]
if dim is None:
dim = max([self.entities[k]["dim"] for k in obj_IDs])
self.model.addPhysicalGroup(dim, obj_IDs, phys_ID)
if name is None:
name = "domain " + str(self.n_domains)
self.id_domains = np.append(self.id_domains, phys_ID)
self.dim_domains = np.append(self.dim_domains, dim)
self.name_domain += [name]
@property
def n_domains(self):
return len(self.id_domains)
@property
def domains_1D(self):
I = np.where(self.dim_domains == 1)
return self.id_domains[I]
@property
def domains_2D(self):
I = np.where(self.dim_domains == 2)
return self.id_domains[I]
@property
def domains_3D(self):
I = np.where(self.dim_domains == 3)
return self.id_domains[I]
##############################################################################################
####################################### feilds methods ####################################
##############################################################################################
[docs]
def refine_entities(self, ent_ID, res_in, dim, res_out=None, IncludeBoundary=True):
"""
refine mesh resolution in a list of faces or volumes IDs
Parameters
----------
ent_ID : list[int] or int
ID or list of ID of the entities where the resolution should be changed
res_in : float
resolution inside the entities
dim : int (2 or 3)
dimention of the considerated entities
res_out : float
resolution outside the entities if None take self.res, default None
IncludeBoundary : bool
include the boundaries for the refinment, default True
"""
self.Nfeild += 1
if not np.iterable(ent_ID):
ent_ID = [ent_ID]
if dim == 2:
typelist = "SurfacesList"
else:
typelist = "VolumesList"
if res_out is None:
res_out = self.res
self.model.mesh.field.add("Constant", self.Nfeild)
self.model.mesh.field.setNumbers(self.Nfeild, typelist, ent_ID)
self.model.mesh.field.setNumber(self.Nfeild, "IncludeBoundary", IncludeBoundary)
self.model.mesh.field.setNumber(self.Nfeild, "VIn", res_in)
self.model.mesh.field.setNumber(self.Nfeild, "VOut", res_out)
return self.Nfeild
[docs]
def refine_threshold(
self, ent_ID, dim, dist_min, dist_max, res_min, res_max=None, N_samples=100
):
"""
refine mesh resolution in a list of faces or volumes IDs
Parameters
----------
ent_ID : list[int] or int
ID or list of ID of the entities where the resolution should be changed
res_in : float
resolution inside the entities
dim : int (2 or 3)
dimention of the considerated entities
res_out : float
resolution outside the entities if None take self.res, default None
IncludeBoundary : bool
include the boundaries for the refinment, default True
"""
self.Nfeild += 1
if not np.iterable(ent_ID):
ent_ID = [ent_ID]
if dim == 2:
typelist = "SurfacesList"
elif dim == 1:
typelist = "CurvesList"
else:
typelist = "PointsList"
if res_max is None:
res_max = self.res
self.Nfeild += 1
i_Dist_field = self.Nfeild
self.model.mesh.field.add("Distance", i_Dist_field)
self.model.mesh.field.setNumbers(i_Dist_field, typelist, ent_ID)
self.model.mesh.field.setNumber(i_Dist_field, "Sampling", N_samples)
self.Nfeild += 1
self.model.mesh.field.add("Threshold", self.Nfeild)
self.model.mesh.field.setNumber(self.Nfeild, "InField", i_Dist_field)
self.model.mesh.field.setNumber(self.Nfeild, "SizeMin", res_min)
self.model.mesh.field.setNumber(self.Nfeild, "SizeMax", res_max)
self.model.mesh.field.setNumber(self.Nfeild, "DistMin", dist_min)
self.model.mesh.field.setNumber(self.Nfeild, "DistMax", dist_max)
self.model.mesh.field.setNumber(self.Nfeild, "Sigmoid", 1)
return self.Nfeild
[docs]
def refine_min(self, feild_IDs):
"""
refine mesh resolution taking the minimum value for a list of refinment fields
Parameters
----------
feild_IDs : list[int]
list of field from wich the minimum should be taken
"""
self.Nfeild += 1
self.model.mesh.field.add("Min", self.Nfeild)
self.model.mesh.field.setNumbers(self.Nfeild, "FieldsList", feild_IDs)
return self.Nfeild
[docs]
def refinement_callback(self, meshSizeCallback):
"""
Add a call back function which is apply to the mesh size defined by fields and return
the final mesh size
Parameters
----------
meshSizeCallback : nrv.utils.MeshCallBack
"""
self.model.mesh.setSizeCallback(meshSizeCallback)
##############################################################################################
################################ generate and saving methods ##############################
##############################################################################################
[docs]
def generate(self):
"""
genetate the mesh
"""
if not self.is_generated:
if self.Nfeild > 0:
self.model.mesh.field.setAsBackgroundMesh(self.Nfeild)
else:
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.1 * self.res)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", self.res)
self.model.occ.synchronize()
self.model.mesh.generate(self.D)
self.is_generated = True
[docs]
def save(self, fname, generate=True):
"""
Save mesh in fname in ".msh"
Parameters
----------
fname : str
path and name of saving file. If ends with ".msh" only save in ".msh" file
"""
if generate == True:
self.generate()
fname = rmv_ext(fname)
gmsh.write(fname + ".msh")
self.file = fname
[docs]
def load(self, fname):
"""
load mesh from ".msh" file
"""
if ".msh" not in fname:
fname += ".msh"
if not os.path.isfile(fname):
rise_warning(fname + " not found: Mesh could not be loaded")
else:
gmsh.merge(fname)
self.file = fname
[docs]
def visualize(self, fname=None):
if fname is None:
self.generate()
gmsh.fltk.run()
elif fname is not None:
self.save(fname=fname)
os.system("gmsh " + self.file + ".msh")