import dolfinx
import numpy as np
import ufl
from ..io import *
from mpi4py import MPI
from petsc4py import PETSc
[docs]
class Mesh():
"""
A base class for creating and managing computational meshes using Dolfinx.
Parameters
-------------
comm: MPI communicator, default is MPI.COMM_WORLD
mesh: dolfinx mesh object, default is None
mesh_file: mesh file name, default is None
boundary: dolfinx mesh tags for boundary, default is None
subdomain: dolfinx mesh tags for subdomain, default is None
gdim: geometric dimension of the mesh, optional, used when loading from file
"""
def __init__(self, **kwargs):
self.comm = kwargs.get('comm', MPI.COMM_WORLD)
# --- Get the mesh -- - #
_msh = kwargs.get('mesh', None)
_msh_file = kwargs.get('mesh_file', None)
# --- Load fenics mesh objects --- #
if _msh is not None:
# Assigm dolfinx mesh object to flationx mesh object
self.msh = _msh
# Set empty boundary and subdomain if not provided, otherwise use the provided ones
_boundary = kwargs.get('boundary', None)
_subdomain = kwargs.get('subdomain', None)
if _boundary is None:
self.boundary = dolfinx.mesh.meshtags(self.msh, self.get_fdim(), np.array([]), np.array([]))
else:
self.boundary = _boundary
if _subdomain is None:
self.subdomain = dolfinx.mesh.meshtags(self.msh, self.get_tdim(), np.array([]), np.array([]))
else:
self.subdomain = _subdomain
elif _msh_file is not None:
# Option to provide gdim instead of having the mesh reader guess it
_gdim = kwargs.get('gdim', None)
self.msh, self.subdomain, self.boundary = io.read_mesh(_msh_file, gdim=_gdim, comm=self.comm)
[docs]
def get_tdim(self):
"""
Returns the topological dimension of the mesh.
"""
return self.msh.topology.dim
[docs]
def get_fdim(self):
"""
Returns the facet dimension of the mesh; the dimension of the mesh boundary.
"""
return self.msh.topology.dim - 1
[docs]
def get_gdim(self):
"""
Returns the geometric dimension of the mesh.
"""
return self.msh.geometry.dim
[docs]
def get_cell_diameter(self):
"""
Returns the ufl cell diameter of the mesh.
"""
return ufl.CellDiameter(self.msh)
[docs]
def get_facet_normal(self):
"""
Returns the ufl facet normal of the mesh.
"""
return ufl.FacetNormal(self.msh)
def _mark_entities(self, marking_dict, entity_dim):
"""
Mark entities (facets or cells) of the mesh based on user-defined markers. Method called by mark_boundary and mark_subdomain.
Parameters
-----------
marking_dict : dict, A dictionary where keys are marker ids and values are functions that return boolean arrays
indicating which entities to mark.
entity_dim : int, The dimension of the entities to mark (e.g., facets or cells).
Returns
-------
dolfinx.mesh.MeshTags: A MeshTags object containing the marked entities.
"""
# Build connectivity
if entity_dim != self.get_tdim():
self.msh.topology.create_connectivity(entity_dim, self.get_tdim())
# Create a list of entities and marking ids
entity_ids = []
all_marking_ids = []
# for marker, idx in zip(marking_function_lst, marking_ids):
for idx, marker in marking_dict.items():
this_bnd = dolfinx.mesh.locate_entities(self.msh, entity_dim, marker)
this_ids = [idx for i in range(len(this_bnd))]
entity_ids.extend(this_bnd)
all_marking_ids.extend(this_ids)
# Sort marking and entities by entity ids
entity_ids = np.array(entity_ids)
all_marking_ids = np.array(all_marking_ids)
sorted_ids = np.argsort(entity_ids)
entity_ids = entity_ids[sorted_ids]
all_marking_ids = all_marking_ids[sorted_ids]
return dolfinx.mesh.meshtags(self.msh, entity_dim, entity_ids, all_marking_ids)
[docs]
def mark_boundary(self, marking_dict):
"""
Mark the boundary of the mesh with user-defined markers.
Parameters
-----------
marking_dict : dict, A dictionary where keys are marker ids and values are functions that return boolean arrays
indicating which facets to mark.
"""
self.boundary = self._mark_entities(marking_dict, self.get_fdim())
[docs]
def mark_subdomain(self, marking_dict):
"""
Mark the subdomains of the mesh with user-defined markers.
Parameters
-----------
marking_dict : dict, A dictionary where keys are marker ids and values are functions that return boolean arrays
indicating which cells to mark.
"""
self.subdomain = self._mark_entities(marking_dict, self.get_tdim())
[docs]
def get_num_facets_local(self):
"""
Returns the number of facets in the local mesh partition.
"""
return self.msh.topology.index_map(self.get_tdim()-1).size_local
[docs]
def get_mean_boundary_normal(self, boundary_id):
"""
Calculate the mean normal vector of a given boundary.
Parameters
------------
boundary_id (int): The id of the boundary for which to compute the mean normal.
Returns
------------
np.ndarray: A unit vector representing the mean outward normal of the boundary.
"""
n = self.get_facet_normal()
ds = ufl.Measure('ds', domain=self.msh, subdomain_data=self.boundary)
local_normal_array = []
for i in range(self.get_gdim()):
form_i = dolfinx.fem.form(n[i] * ds(boundary_id))
assemble_i = dolfinx.fem.assemble_scalar(form_i)
local_normal_array.append(assemble_i)
local_normal_array = np.array(local_normal_array)
global_normal_array = self.msh.comm.allreduce(local_normal_array, op=MPI.SUM)
norm = np.linalg.norm(global_normal_array, 2)
normal_array = global_normal_array / norm
return normal_array
[docs]
def get_boundary_centroid(self, boundary_id):
"""
Calculate the centroid of a given boundary.
Parameters
------------
boundary_id (int): The id of the boundary for which to compute the centroid.
Returns
------------
np.ndarray: A point representing the centroid of the boundary.
"""
# Get the facets of the boundary with the given id
boundary_facets = self.boundary.find(boundary_id)
facet_midpoints = []
x = self.msh.geometry.x
# Loop through each facet and compute the midpoint
for facet in boundary_facets:
facet_vertices = self.msh.topology.connectivity(self.get_fdim(), 0).links(facet)
facet_coords = x[facet_vertices]
facet_midpoints.append(np.mean(facet_coords, axis=0))
# Gather midpoints from all processes
if len(facet_midpoints) > 0:
midpoints_local = np.array(facet_midpoints)
else:
midpoints_local = np.empty((0, 3))
all_midpoints = self.msh.comm.allgather(midpoints_local)
centroid = None
# If root process, compute the centroid
if self.msh.comm.rank == 0:
# Concatenate all midpoints from all processes
all_midpoints_flat = np.concatenate(all_midpoints)
# Compute the mean to get the centroid
if len(all_midpoints_flat) > 0:
centroid = np.mean(all_midpoints_flat, axis=0)
# Broadcast the result to all processes
centroid = np.array(self.msh.comm.bcast(centroid, root=0))
return centroid
[docs]
def get_boundary_area(self, boundary_id):
"""
Calculate the area of a given boundary.
Parameters
------------
boundary_id (int): The id of the boundary for which to compute the area.
Returns
------------
float: The area of the boundary.
"""
one = dolfinx.fem.Constant(self.msh, PETSc.ScalarType(1.0))
ds = ufl.Measure('ds', domain=self.msh, subdomain_data=self.boundary)
form = dolfinx.fem.form(one * ds(boundary_id))
area = dolfinx.fem.assemble_scalar(form)
# Sum the area across all processes
area = self.msh.comm.allreduce(area, op=MPI.SUM)
return area
[docs]
def get_mean_cell_diameter(self):
"""
Get mean cell diameter.
Returns
---------
float: the average length of the cell diameters within the mesh.
"""
cell_diameter = self.get_cell_diameter()
integrated_cell_diameter = dolfinx.fem.assemble_scalar(dolfinx.fem.form(cell_diameter * ufl.dx))
volume = dolfinx.fem.assemble_scalar(dolfinx.fem.form(dolfinx.fem.Constant(self.msh, dolfinx.default_scalar_type(1.0)) * ufl.dx))
mean_cell_diameter = integrated_cell_diameter / volume
return mean_cell_diameter
[docs]
class Boundary:
"""
A simple class to represent and store properties of a mesh boundary.
Parameters
----------
mesh : flatiron_tk.Mesh
The mesh object containing the boundary.
boundary_id : int
The identifier for the specific boundary.
**kwargs : dict
Additional attributes to be set for the boundary.
Default Attributes
------------------
id : int
The boundary identifier.
area : float
The area of the boundary.
radius : float
The equivalent radius of the boundary (assuming circular shape).
centroid : np.ndarray
The centroid coordinates of the boundary.
normal : np.ndarray
The mean outward normal vector of the boundary.
"""
def __init__(self, mesh, boundary_id, **kwargs):
self.id = boundary_id
self.area = mesh.get_boundary_area(boundary_id)
self.radius = None
if mesh.msh.geometry.dim == 3: self.radius = np.sqrt(self.area / np.pi)
elif mesh.msh.geometry.dim == 2: self.radius = self.area / 2
self.centroid = mesh.get_boundary_centroid(boundary_id)
self.normal = mesh.get_mean_boundary_normal(boundary_id)
# Set additional attributes from kwargs
for key, value in kwargs.items():
setattr(self, key, value)
[docs]
def set_new_attributes(self, **kwargs):
"""
Sets new attributes or updates existing ones for the Boundary instance.
Parameters
----------
**kwargs : dict
Key-value pairs representing the attributes to be set or updated.
"""
for key, value in kwargs.items():
if hasattr(self, key):
print(f'\033[91m[WARNING]\033[0m Attribute \"{key}\" already exists and will be overwritten.')
setattr(self, key, value)