diff --git a/src/bvpy/domains/abstract.py b/src/bvpy/domains/abstract.py
index 56c1c57d652e966f323331cb207695285c3e0a05..dc3b1a5a3f9a5e99a827b161ca5d809fb251da88 100644
--- a/src/bvpy/domains/abstract.py
+++ b/src/bvpy/domains/abstract.py
@@ -22,7 +22,10 @@
import gmsh
import sys
+import os
+import shutil
from bvpy import logger as log
+from bvpy import MPI_COMM_WORLD
import meshio as mio
from abc import ABC, abstractmethod
import tempfile
@@ -157,7 +160,7 @@ class AbstractDomain(ABC):
AbstractDomain._domain_instances += 1
if AbstractDomain._domain_instances == 1:
- gmsh.initialize()
+ gmsh.initialize() # for all processes
if clear:
gmsh.clear()
@@ -167,8 +170,15 @@ class AbstractDomain(ABC):
self._cell_size = cell_size
self._algorithm = algorithm
- self._mesh_dir = tempfile.TemporaryDirectory()
- self._dir_name = self._mesh_dir.name
+ if MPI_COMM_WORLD.Get_rank() == 0:
+ # Creation of the temporary folder only for the process of rank 0
+ self._mesh_dir = tempfile.TemporaryDirectory() # keep a reference to avoid suppression of the temp file
+ dir_name = self._mesh_dir .name
+ else:
+ dir_name = None
+
+ # Diffusion of the folder to the other processess
+ self._dir_name = MPI_COMM_WORLD.bcast(dir_name, root=0)
self.set_cell_type(cell_type)
self.set_verbosity(verbosity)
@@ -194,10 +204,16 @@ class AbstractDomain(ABC):
"""Clears the temporary file containing the domain.
"""
- self._mesh_dir.cleanup()
+ #self._mesh_dir.cleanup()
+
+ if MPI_COMM_WORLD.Get_rank() == 0:
+ # Remove the temporary folder and its content using process of rank 0
+ if self._dir_name and os.path.isdir(self._dir_name):
+ shutil.rmtree(self._dir_name)
+
AbstractDomain._domain_instances -= 1
if AbstractDomain._domain_instances == 0:
- gmsh.finalize()
+ gmsh.finalize() # for all processes
def __repr__(self):
"""Returns a short description of the considered domain.
@@ -300,6 +316,13 @@ class AbstractDomain(ABC):
"""
pass
+ def create_geometry(self, *args, **kwargs):
+ """ Routine to control the update of the geometry if multi-threading is used"""
+
+ if MPI_COMM_WORLD.Get_rank() == 0:
+ log.debug('Create geometry of the domain!')
+ # - only the process of rank 0 updates the geometry
+ self.geometry(*args, **kwargs)
def run_gui(self):
"""Generates a quick gmsh gui with the current model in it.
@@ -433,7 +456,7 @@ class AbstractDomain(ABC):
if self.mesh is not None:
self.model.mesh.clear()
- self._mesh_conversion()
+ self._mesh_conversion_mpi()
self.ds = fe.Measure('ds', domain=self.mesh, subdomain_data=self.bdata)
self.dx = fe.Measure('dx', domain=self.mesh, subdomain_data=self.cdata)
@@ -463,6 +486,7 @@ class AbstractDomain(ABC):
ele_dim = CELL_TYPE[self._cell_type]
log.debug(f"Convert a mesh with topological element of dimension {ele_dim}")
+ # Seul le processus de rang 0 génère et écrit le maillage
self.model.mesh.generate(ele_dim)
# - Assert that we will get all the entites of the gmsh mesh during the conversion by catching the number of
@@ -618,6 +642,204 @@ class AbstractDomain(ABC):
self.bdata = _generate_default_boundary_data(mesh)
self.sub_boundary_names = {0: 'boundary'}
+ def _mesh_conversion_mpi(self):
+ """Converts a temporary gmsh mesh into a fenics (dolphin) one.
+
+ Returns
+ -------
+ mesh : :class:`Mesh<dolfin.cpp.mesh.Mesh>`
+ Fenics mesh corresponding to the domain.
+ bdata : :class:`MeshFunctionSizet<dolfin.cpp.mesh.MeshFunction>`
+ MeshFunction that labels the borders (i.e. outermost facets)
+ of the mesh of the domain.
+ cdata : :class:`MeshFunctionSizet<dolfin.cpp.mesh.MeshFunction>`
+ MeshFunction that labels the inner elements (i.e. cells)
+ of the mesh of the domain.
+
+ """
+ # - Get the finite element topological dimension
+ ele_dim = CELL_TYPE[self._cell_type]
+ log.debug(f"Convert a mesh with topological element of dimension {ele_dim}")
+
+ # - Set default values for the domain & boundaries labels/name
+ has_bdata = False # flag to fill automatically labels to boundary
+ sub_boundary_names = {0: 'boundary'} # default naming for the boundary
+ has_cdata = False # flag to fill automatically labels if needed
+ sub_domain_names = {0: 'domain'} # default naming for the domain
+
+ ### From GMSH to MESHIO
+ # Only the rank 0 process generate the gmsh mesh, convert it to MeshIO and save it as XDMF file.
+ # It will avoid duplicates
+ if MPI_COMM_WORLD.Get_rank() == 0:
+ log.debug('Generate GMSH mesh and convert to MeshIO')
+ self.model.mesh.generate(ele_dim)
+
+ # - Assert that we will get all the entites of the gmsh mesh during the conversion by catching the number of
+ # entities
+ ELEMENT_GMSH_TYPES = {'line': 1, 'triangle': 2, 'tetra': 4} # Gmsh code for a given element type
+ element_type, element_tags, _ = self.model.mesh.getElements(ele_dim) # get the entities of dimension `ele_dim`
+ N = len(element_tags[np.where(element_type == ELEMENT_GMSH_TYPES[self._cell_type])[0][0]]) # number of entities
+
+ gmsh.write(self._dir_name+'/mesh.msh')
+
+ path = self._dir_name + '/mesh.msh'
+ mesh_io = mio.read(path, file_format='gmsh')
+
+ # - Get the coordinates of the vertices
+ ## Note : results if embedded space is != 3 ?
+ mesh_io.points = mesh_io.points[:, :self.dim]
+
+ if self._cell_type not in mesh_io.cells_dict:
+ raise RuntimeError(f'None elements of the required topological dimension (element_dim = {ele_dim} for '
+ f'{self._cell_type}) have been found during the conversion from Gmsh to MeshIo.'
+ f'\nKnown issue 1: no Entities of the required topological dimension ({ele_dim}) '
+ f'are defined in your Gmsh mesh. Use YourAbstractObject.factory.getEntities(dim={ele_dim}) '
+ f'to check if it is the case.\nKnown issue 2: some PhysicalGroups have been set for '
+ f'Entities of topological dimension < {ele_dim}. If PhysicalGroups are defined, be sure '
+ f'to define them for all the Entities at topological dimension = {ele_dim}. If not, '
+ f'all the unmarked Entities will be ignored.\nKnown issue 3: ???')
+ elif len(mesh_io.cells_dict[self._cell_type]) != N:
+ raise Warning(f'Some elements have been lost during the conversion of the gmsh mesh, one possible '
+ f'reason is the partial labelling of entities at topological dimension {ele_dim}. Make '
+ f'sure that you label all the entities ! Unlabelled entities at dim {ele_dim} will be'
+ f' ignored...')
+
+ # - Extract physical properties (and possible name) for the domain
+ if ('gmsh:physical' in mesh_io.cell_data_dict.keys()
+ and self._cell_type in mesh_io.cell_data_dict["gmsh:physical"].keys()):
+ log.debug(f"Found labels at the dimension {ele_dim}")
+ has_cdata = True
+
+ # Notes:
+ # 'gmsh: physical' lists the tags given to group of entities of same dimension. Here, the dimension will
+ # be defined by ele_dim (from the _cell_type).
+ # There tags can be defined by addPhysicalGroup(ele_dim, ele_tags, tagOfGroup) in gmsh.
+ # It is possible to set a name for a given PhysicalGroup using setPhysicalName in gmsh.
+ # Gmsh does not allowed intersection of PhysicalGroup, the first tag for a given entity
+ # will be the one you will find here.
+ # Ex: if elements [0, 1] are tagged with 4,
+ # then 1 is tagged with 2,
+ # you will see ele. 1 with tag 4 here.
+ # Important warning : In the case where some entities have a PhysicialGroup, if it exists entities without
+ # any PhysicialGroup, they will be ignored !
+ # TO BE CLEAR : you will loose part of your mesh !
+ # Conclusion : using PhysicalGroup, the resulting mesh will always contains PhysicialGroups that form a
+ # partition of the mesh (at the given topological dimension)
+ # However, it is possible that some of the PhysicialGroup do not have a PhysicalName, we will
+ # deal with this situation by setting a default name.
+
+ labels = mesh_io.cell_data_dict["gmsh:physical"][self._cell_type]
+
+ # - Catch PhysicalName associated with PhysicalGroup
+ # - reformat field_data as {tag : name} and filter PhysicalName (`ele_dim` topological dimension)
+ if mesh_io.field_data:
+ sub_domain_names = {lab_dim[0]: lab_name
+ for lab_name, lab_dim in mesh_io.field_data.items()
+ if lab_dim[1] == (ele_dim)} # tag: name
+ log.debug(f"domains : field data => {sub_domain_names} found")
+ else:
+ sub_domain_names = {i: f'sub_domain_{i}' for i in np.unique(labels)}
+ log.debug(f"domains : default field data => {sub_domain_names}")
+
+ if set(labels) != set(sub_domain_names):
+ log.warning('Domain : the number of PhysicalGroup does not correspond to the number of PhysicalName !')
+ log.warning('Domain : some of your PhysicalGroup will be automatically named...')
+
+ # - set automatic names for the other labels
+ sub_domain_names.update({lab: f'sub_domain_{i}'
+ for i, lab in enumerate(set(labels) - set(sub_domain_names))})
+
+ mio.write(self._dir_name + '/' + 'mesh.xdmf',
+ mio.Mesh(points=mesh_io.points,
+ cells={self._cell_type:
+ mesh_io.cells_dict[self._cell_type]},
+ cell_data={f"{self._cell_type}_labels": [labels]}))
+ else:
+ log.debug(f"No labels found at the dimension {ele_dim}, use default.")
+
+ mio.write(self._dir_name + '/' + 'mesh.xdmf',
+ mio.Mesh(points=mesh_io.points,
+ cells={self._cell_type:
+ mesh_io.cells_dict[self._cell_type]}))
+
+ # - Extract physical properties (and possible names) for the boundaries (ele_dim - 1)
+ bnd_cell_type = [cell_type for cell_type, d in CELL_TYPE.items() if d == (ele_dim - 1)][0]
+
+ if ('gmsh:physical' in mesh_io.cell_data_dict.keys()
+ and bnd_cell_type in mesh_io.cell_data_dict["gmsh:physical"].keys()):
+ has_bdata = True
+ labels = mesh_io.cell_data_dict["gmsh:physical"][bnd_cell_type]
+
+ # - Catch PhysicalName associated with PhysicalGroup
+ # - reformat field_data as {tag : name} and filter PhysicalName (`ele_dim` - 1 topological dimension)
+ if mesh_io.field_data:
+ sub_boundary_names = {lab_dim[0]: lab_name
+ for lab_name, lab_dim in mesh_io.field_data.items()
+ if lab_dim[1] == (ele_dim - 1)} # tag: name
+ log.debug(f"boundary : field data => {sub_boundary_names} found")
+ else:
+ sub_boundary_names = {i: f'sub_boundary_{i}' for i in np.unique(labels)}
+ log.debug(f"boundary : default field data => {sub_boundary_names}")
+
+ if set(labels) != set(sub_boundary_names):
+ log.warning('Boundary: the number of PhysicalGroup does not correspond to the number of PhysicalName !')
+ log.warning('Boundary: some of your PhysicalGroup will be automatically named...')
+
+ # - set automatic names for the other labels
+ sub_boundary_names.update({lab: f'sub_boundary_{i}'
+ for i, lab in enumerate(set(labels)-set(sub_boundary_names))})
+
+ bnd_mesh = mio.Mesh(points=mesh_io.points,
+ cells={bnd_cell_type:
+ mesh_io.cells_dict[bnd_cell_type]},
+ cell_data={f"boundary": [labels]})
+ mio.xdmf.write(self._dir_name + '/mesh_boundary.xdmf', bnd_mesh)
+
+ # Update state of other processes from rank 0 (labels and sub-domain names)
+ has_cdata = MPI_COMM_WORLD.bcast(has_cdata, root=0)
+ has_bdata = MPI_COMM_WORLD.bcast(has_bdata, root=0)
+ sub_domain_names = MPI_COMM_WORLD.bcast(sub_domain_names, root=0)
+ sub_boundary_names = MPI_COMM_WORLD.bcast(sub_boundary_names, root=0)
+
+ # Synchronize the MPI processes (other processes will wait for the rank 0 for the mesh creation
+ MPI_COMM_WORLD.Barrier()
+
+ ### Reading of the XDMF mesh in parallel (all the processes will get the complete mesh)
+ log.debug('Generate FeniCS mesh !')
+ log.debug(f'{self._cell_type} for {MPI_COMM_WORLD.Get_rank()}')
+
+ # Init FeniCS mesh
+ mesh = fe.Mesh()
+
+ # Read the saved mesh
+ with fe.XDMFFile(MPI_COMM_WORLD, self._dir_name + '/' + 'mesh.xdmf') as f:
+ f.read(mesh)
+
+ # Fill the cdata
+ if has_cdata:
+ mvc = fe.MeshValueCollection("size_t", mesh, ele_dim)
+ with fe.XDMFFile(MPI_COMM_WORLD, self._dir_name + '/' + 'mesh.xdmf') as infile:
+ infile.read(mvc, f'{self._cell_type}_labels')
+ cdata = fe.cpp.mesh.MeshFunctionSizet(mesh, mvc)
+ else:
+ cdata = _generate_default_cell_data(mesh)
+
+ # Fill the bdata
+ if has_bdata:
+ mvc = fe.MeshValueCollection("size_t", mesh, ele_dim - 1)
+ with fe.XDMFFile(MPI_COMM_WORLD, self._dir_name + '/mesh_boundary.xdmf') as infile:
+ infile.read(mvc, 'boundary')
+ bdata = fe.cpp.mesh.MeshFunctionSizet(mesh, mvc)
+ else:
+ bdata = _generate_default_boundary_data(mesh)
+
+ # Stored the mesh information
+ self.mesh = mesh
+ self.cdata = cdata
+ self.sub_domain_names = sub_domain_names
+ self.bdata = bdata
+ self.sub_boundary_names = sub_boundary_names
+
def info(self):
"""Prints the descrition of the domain.
@@ -628,6 +850,20 @@ class AbstractDomain(ABC):
"""
print(self._details)
+class GmshProxy:
+ def __init__(self, factory):
+ self._factory = factory
+
+ def __getattr__(self, name):
+ """Redirige les appels de méthode vers factory si le rang est 0."""
+ def method(*args, **kwargs):
+ log.debug(f'Call from {MPI_COMM_WORLD.Get_rank()} of the method {name}')
+ if MPI_COMM_WORLD.Get_rank() == 0:
+ log.debug(f'Ok for {MPI_COMM_WORLD.Get_rank()} of the method {name}')
+ return getattr(self._factory, name)(*args, **kwargs)
+ else:
+ log.debug(f'NO for {MPI_COMM_WORLD.Get_rank()} of the method {method}')
+ return method
class BuiltInModel(object):
"""A model based on the Gmsh `geo` factory.
@@ -672,7 +908,8 @@ class OccModel(object):
"""
super().__init__()
- self.factory = gmsh.model.occ
+ #self.factory = gmsh.model.occ
+ self.factory = GmshProxy(gmsh.model.occ)
def __sub__(self, other):
"""Cutting the left hand side Domain with the right hand side domain.
diff --git a/src/bvpy/domains/primitives.py b/src/bvpy/domains/primitives.py
index 7b8e18d959d19fd976a7a9e26716ac5403731a06..8c759f0ced54089499fdba466f86e24671ded081 100644
--- a/src/bvpy/domains/primitives.py
+++ b/src/bvpy/domains/primitives.py
@@ -124,7 +124,7 @@ class Rectangle(AbstractDomain, OccModel):
super(Rectangle, self).__init__(**kwargs)
self._characteristic_size = {'Length': length,
'Width': width}
- self.geometry(x, y, length, width)
+ self.create_geometry(x, y, length, width)
def geometry(self, x, y, length, width):
"""Adds a rectangle to the Gmsh factory computing the domain.
@@ -492,7 +492,7 @@ class Cylinder(AbstractDomain, OccModel):
super(Cylinder, self).__init__(**kwargs)
self._characteristic_size = {'Length': max(dx, dy, dz),
'Width': min(dx, dy, dz)}
- self.geometry(x, y, z, dx, dy, dz, r, open)
+ self.create_geometry(x, y, z, dx, dy, dz, r, open)
def geometry(self, x, y, z, dx, dy, dz, r, open):
self.volumes[0] = self.factory.addCylinder(x, y, z, dx, dy, dz, r)