import numpy as nm
import scipy.sparse as sps
from scipy.sparse.csgraph import connected_components
from scipy.spatial.transform import Rotation
from scipy.spatial import cKDTree
from sfepy.base.compat import factorial
from sfepy.base.base import output
from sfepy.discrete.equations import create_dof_graph
from sfepy.discrete.fem import Mesh, FEDomain
from sfepy.discrete.fem.utils import prepare_translate
[docs]
def elems_q2t(el):
nel, nnd = el.shape
if nnd > 4:
q2t = nm.array([[0, 2, 3, 6],
[0, 3, 7, 6],
[0, 7, 4, 6],
[0, 5, 6, 4],
[1, 5, 6, 0],
[1, 6, 2, 0]])
else:
q2t = nm.array([[0, 1, 2],
[0, 2, 3]])
ns, nn = q2t.shape
nel *= ns
out = nm.zeros((nel, nn), dtype=nm.int32)
for ii in range(ns):
idxs = nm.arange(ii, nel, ns)
out[idxs,:] = el[:, q2t[ii,:]]
return nm.ascontiguousarray(out)
[docs]
def triangulate(mesh, verbose=False):
"""
Triangulate a 2D or 3D tensor product mesh: quadrilaterals->triangles,
hexahedrons->tetrahedrons.
Parameters
----------
mesh : Mesh
The input mesh.
Returns
-------
mesh : Mesh
The triangulated mesh.
"""
conns = None
for k, new_desc in [('3_8', '3_4'), ('2_4', '2_3')]:
if k in mesh.descs:
conns = mesh.get_conn(k)
break
if conns is not None:
nelo = conns.shape[0]
output('initial mesh: %d elements' % nelo, verbose=verbose)
new_conns = elems_q2t(conns)
nn = new_conns.shape[0] // nelo
new_cgroups = nm.repeat(mesh.cmesh.cell_groups, nn)
output('new mesh: %d elements' % new_conns.shape[0], verbose=verbose)
mesh = Mesh.from_data(mesh.name, mesh.coors,
mesh.cmesh.vertex_groups,
[new_conns], [new_cgroups], [new_desc])
return mesh
[docs]
def smooth_mesh(mesh, n_iter=4, lam=0.6307, mu=-0.6347,
weights=None, bconstr=True,
volume_corr=False):
"""
FE mesh smoothing.
Based on:
[1] Steven K. Boyd, Ralph Muller, Smooth surface meshing for automated
finite element model generation from 3D image data, Journal of
Biomechanics, Volume 39, Issue 7, 2006, Pages 1287-1295,
ISSN 0021-9290, 10.1016/j.jbiomech.2005.03.006.
(http://www.sciencedirect.com/science/article/pii/S0021929005001442)
Parameters
----------
mesh : mesh
FE mesh.
n_iter : integer, optional
Number of iteration steps.
lam : float, optional
Smoothing factor, see [1].
mu : float, optional
Unshrinking factor, see [1].
weights : array, optional
Edge weights, see [1].
bconstr: logical, optional
Boundary constraints, if True only surface smoothing performed.
volume_corr: logical, optional
Correct volume after smoothing process.
Returns
-------
coors : array
Coordinates of mesh nodes.
"""
def laplacian(coors, weights):
n_nod = coors.shape[0]
displ = (weights - sps.eye_array(n_nod)) @ coors
return displ
def taubin(coors0, weights, lam, mu, n_iter):
coors = coors0.copy()
for ii in range(n_iter):
displ = laplacian(coors, weights)
if nm.mod(ii, 2) == 0:
coors += lam * displ
else:
coors += mu * displ
return coors
def get_volume(el, nd):
from sfepy.linalg.utils import dets_fast
dim = nd.shape[1]
nnd = el.shape[1]
etype = '%d_%d' % (dim, nnd)
if etype == '2_4' or etype == '3_8':
el = elems_q2t(el)
nel = el.shape[0]
#bc = nm.zeros((dim, ), dtype=nm.double)
mul = 1.0 / factorial(dim)
if dim == 3:
mul *= -1.0
mtx = nm.ones((nel, dim + 1, dim + 1), dtype=nm.double)
mtx[:,:,:-1] = nd[el,:]
vols = mul * dets_fast(mtx)
vol = vols.sum()
bc = nm.dot(vols, mtx.sum(1)[:,:-1] / nnd)
bc /= vol
return vol, bc
from sfepy.base.timing import Timer
output('smoothing...')
timer = Timer(start=True)
if weights is None:
n_nod = mesh.n_nod
domain = FEDomain('mesh', mesh)
cmesh = domain.cmesh
# initiate all vertices as inner - hierarchy = 2
node_group = nm.ones((n_nod,), dtype=nm.int16) * 2
# boundary vertices - set hierarchy = 4
if bconstr:
# get "vertices of surface"
facets = cmesh.get_surface_facets()
f_verts = cmesh.get_incident(0, facets, cmesh.dim - 1)
node_group[f_verts] = 4
# generate costs matrix
e_verts = cmesh.get_conn(1, 0).indices
fc1, fc2 = e_verts[0::2], e_verts[1::2]
idxs = nm.where(node_group[fc2] >= node_group[fc1])
rows1 = fc1[idxs]
cols1 = fc2[idxs]
idxs = nm.where(node_group[fc1] >= node_group[fc2])
rows2 = fc2[idxs]
cols2 = fc1[idxs]
crows = nm.concatenate((rows1, rows2))
ccols = nm.concatenate((cols1, cols2))
costs = sps.coo_array((nm.ones_like(crows), (crows, ccols)),
shape=(n_nod, n_nod),
dtype=nm.double)
# generate weights matrix
idxs = list(range(n_nod))
aux = sps.coo_array((1.0 / nm.asarray(costs.sum(1)).squeeze(),
(idxs, idxs)),
shape=(n_nod, n_nod),
dtype=nm.double)
#aux.setdiag(1.0 / costs.sum(1))
weights = (aux.tocsc() @ costs.tocsc()).tocsr()
coors = taubin(mesh.coors, weights, lam, mu, n_iter)
output('...done in %.2f s' % timer.stop())
if volume_corr:
output('rescaling...')
timer.start()
volume0, bc = get_volume(mesh.conns[0], mesh.coors)
volume, _ = get_volume(mesh.conns[0], coors)
scale = volume0 / volume
output('scale factor: %.2f' % scale)
coors = (coors - bc) * scale + bc
output('...done in %.2f s' % timer.stop())
return coors
[docs]
def merge_lines(mesh, eps=1e-18):
"""
Merge edges of an edge-only mesh that are in the same direction w.r.t. the
tolerance `eps`.
"""
coors, ngroups, conns, mat_ids, ctype = mesh
conns = conns[0]
# vertices to edges map
n_v = coors.shape[0]
n_e = conns.shape[0]
row = nm.repeat(nm.arange(n_e), 2)
aux = sps.coo_array((nm.ones((n_e * 2,), dtype=bool),
(row, conns.flatten())), shape=(n_e, n_v))
v2e = aux.tocsc()
n_epv = nm.diff(v2e.indptr)
# directional vectors of edges
de = coors[conns[:, 1], :] - coors[conns[:, 0], :]
de = de / nm.linalg.norm(de, axis=1)[:, nm.newaxis]
eflag = nm.ones((n_e, ), dtype=bool)
valid_e = nm.where(eflag)[0]
e_remove = []
while len(valid_e) > 0:
ie = valid_e[0]
d = de[ie]
buff = [(ie, conns[ie, 0]), (ie, conns[ie, 1])]
eflag[ie] = False # invalidate edge
while len(buff) > 0:
e, v = buff.pop(-1)
if n_epv[v] == 2:
idx = v2e.indptr[v]
aux = v2e.indices[idx]
next_e = v2e.indices[idx + 1] if aux == e else aux
if not eflag[next_e]: # valid edge?
continue
if nm.linalg.norm(de[next_e] - d) < eps\
or nm.linalg.norm(de[next_e] + d) < eps:
next_ec = conns[next_e, :]
new_v = next_ec[0] if next_ec[1] == v else next_ec[1]
idx = 0 if conns[e, 0] == v else 1
conns[e, idx] = new_v # reconnect edge
idx = v2e.indptr[new_v]
aux = v2e.indices[idx]
idx += 0 if aux == next_e else 1
v2e.indices[idx] = e # update v2e map
buff.append((e, new_v)) # continue in searching
eflag[next_e] = False # invalidate edge
e_remove.append(next_e)
valid_e = nm.where(eflag)[0]
if len(e_remove) > 0:
# remove unused edges and vertices
eflag.fill(True)
eflag[nm.asarray(e_remove)] = False
remap = -nm.ones((n_v, ), dtype=nm.int64)
remap[conns[eflag, :]] = 1
vidx = nm.where(remap > 0)[0]
remap[vidx] = nm.arange(len(vidx))
conns_new = remap[conns[eflag, :]]
return coors[vidx, :], ngroups[vidx],\
[conns_new], [mat_ids[0][eflag]], ctype
else:
return mesh
def _get_facets(vertices, offsets, ii, n_fp):
facets = []
for ic in range(n_fp):
facets.append(vertices[offsets[ii] + ic][:, None])
facets = nm.concatenate(facets, axis=1)
return nm.ascontiguousarray(facets.astype(nm.int32))
[docs]
def get_surface_faces(domain):
cmesh = domain.cmesh
faces = cmesh.get_surface_facets()
vertices_f, offs_f = cmesh.get_incident(0, faces,
cmesh.dim - 1, ret_offsets=True)
n_fp = nm.diff(offs_f)
surf_faces = []
itri = nm.where(n_fp == 3)[0]
if itri.size:
surf_faces.append(_get_facets(vertices_f, offs_f, itri, 3))
itet = nm.where(n_fp == 4)[0]
if itet.size:
surf_faces.append(_get_facets(vertices_f, offs_f, itet, 4))
cells_c, offs_c = cmesh.get_incident(cmesh.dim, faces, cmesh.dim - 1,
ret_offsets=True)
ids = cmesh.get_local_ids(faces, cmesh.dim - 1, cells_c, offs_c,
cmesh.dim)
lst = nm.c_[cells_c, ids]
return lst, surf_faces
[docs]
def surface_graph(surf_faces, n_nod):
graph = 0
for conn in surf_faces:
graph += create_dof_graph(conn, conn, (n_nod, n_nod)).tocsr()
return graph
[docs]
def surface_components(gr_s, surf_faces):
"""
Determine surface components given surface mesh connectivity graph.
"""
n_comp, flag = connected_components(gr_s)
comps = []
for ii, face in enumerate(surf_faces):
comp = flag[face[:,0]]
comps.append(comp)
return n_comp, comps
[docs]
def get_cell_vertices_only(mesh):
"""
Remove vertices not used in any cell.
Parameters
----------
mesh: Mesh
FE mesh
Returns
-------
out: Mesh
FE mesh
"""
data = list(mesh._get_io_data())
vertices = nm.unique([nm.unique(conn.ravel()) for conn in data[2]])
remap = prepare_translate(vertices, nm.arange(len(vertices)))
coors = data[0][vertices]
ngroups = data[1][vertices]
conns = [remap[conn] for conn in data[2]]
return Mesh.from_data(mesh.name + '_cleaned', coors, ngroups,
conns, data[3], data[4])
[docs]
def get_mesh_by_cgroup(mesh, value, cell_vertices_only=True):
"""
Extract mesh cells using cell_group value(s).
Parameters
----------
mesh: Mesh
FE mesh
value: int, list, or tuple
Select cells with a given value, list, or range: (value[0], value[1])
cell_vertices_only: bool
If True, remove free vertices
Returns
-------
out: Mesh
FE mesh
"""
new_conns, new_mat_ids, new_descs = [], [], []
for desc in mesh.descs:
conns, cidxs = mesh.get_conn(desc, ret_cells=True)
mat_ids = mesh.cmesh.cell_groups[cidxs]
if isinstance(value, tuple):
ncidxs = nm.logical_and(mat_ids >= value[0], mat_ids <= value[1])
elif isinstance(value, list):
ncidxs = nm.logical_or.reduce([mat_ids == k for k in value])
else:
ncidxs = mat_ids == value
if ncidxs.sum() > 0:
new_conns.append(conns[ncidxs])
new_mat_ids.append(mat_ids[ncidxs])
new_descs.append(desc)
out = Mesh.from_data(mesh.name, mesh.coors,
mesh.cmesh.vertex_groups,
new_conns, new_mat_ids, new_descs)
if cell_vertices_only:
out = get_cell_vertices_only(out)
return out
[docs]
def get_mesh_by_ngroup(mesh, value, cell_vertices_only=True):
"""
Extract mesh cells using cell_group value(s).
Parameters
----------
mesh: Mesh
FE mesh
value: int, list, or tuple
Select nodes with a given value, list, or range: (value[0], value[1])
cell_vertices_only: bool
If True, remove free vertices
Returns
-------
out: Mesh
FE mesh
"""
vgroups = mesh.cmesh.vertex_groups
if isinstance(value, tuple):
vidxs = nm.logical_and(vgroups >= value[0], vgroups <= value[1])
elif isinstance(value, list):
vidxs = nm.logical_or.reduce([vgroups == k for k in value])
else:
vidxs = vgroups == value
new_vgroups = vgroups[vidxs]
new_coors = mesh.coors[vidxs]
remap = -nm.ones((len(vidxs),), dtype=nm.int64)
remap[vidxs] = nm.arange(vidxs.sum())
new_conns, new_mat_ids, new_descs = [], [], []
for desc in mesh.descs:
conns, cidxs = mesh.get_conn(desc, ret_cells=True)
mat_ids = mesh.cmesh.cell_groups[cidxs]
ncidxs = (remap[conns] >= 0).sum(axis=1) == conns.shape[1]
if ncidxs.sum() > 0:
new_conns.append(remap[conns[ncidxs]])
new_mat_ids.append(mat_ids[ncidxs])
new_descs.append(desc)
out = Mesh.from_data(mesh.name, new_coors, new_vgroups,
new_conns, new_mat_ids, new_descs)
if cell_vertices_only:
out = get_cell_vertices_only(out)
return out
[docs]
def stack_descs(conns, mat_ids, descs):
sdescs = set(descs)
sconns = {k: [] for k in sdescs}
smat_ids = {k: [] for k in sdescs}
for k, desc in enumerate(descs):
sconns[desc].append(conns[k])
smat_ids[desc].append(mat_ids[k])
return ([nm.vstack(sconns[k]) for k in sdescs],
[nm.hstack(smat_ids[k]) for k in sdescs],
list(sdescs))
[docs]
def extrude(mesh, cline, twist=None, scale=None, nvec=None,
wedge_to_tetra=True):
"""
Create a solid 3D mesh from a given planar 2D mesh by extruding it.
The new points in each layer lie in a plane perpendicular to one of the
cline segments. Each layer can be twisted and scaled.
Parameters
----------
mesh: Mesh
2D planar FE mesh (tri or quad elements)
cline: list of coordinates
Points of central line.
twist: float, or array
Angle of twist in each layer.
scale: float, or array
Scale factor in each layer.
nvec: array
Normal vectors of layers.
wedge_to_tetra: bool
If True, convert wedge elements to tetrahedrons.
Returns
-------
out: Mesh
3D FE mesh
"""
def get_rotation(v0, v1):
cross = nm.cross(v0, v1)
dot = nm.dot(v0, v1)
if nm.linalg.norm(cross) < 1e-18:
# Already aligned or opposite
if dot > 0:
rotation = Rotation.identity()
else:
# pi/2 rotation around any perpendicular axis
axis = nm.array([0, 1, 0])
if nm.allclose(v0, axis):
axis = nm.array([1, 0, 0])
axis = nm.cross(v0, axis)
axis /= nm.linalg.norm(axis)
rotation = Rotation.from_rotvec(nm.pi * axis)
else:
axis = cross / nm.linalg.norm(cross)
angle = nm.arccos(nm.clip(dot, -1.0, 1.0))
rotation = Rotation.from_rotvec(angle * axis)
return rotation
coors0 = mesh.coors
if coors0.shape[1] == 2:
coors0 = nm.hstack([coors0, coors0[:, [0]] * 0])
ccoor = nm.mean(coors0, axis=0)
coors0 = coors0 - ccoor[None, :]
nvec0 = nm.array([0, 0, 1.])
cline = nm.array(cline)
if len(cline.shape) == 1:
dz, nz = cline[0], int(cline[1])
cline = nm.tile(ccoor, (nz + 1, 1))
cline[:, 2] = nm.arange(nz + 1) * dz
if nvec is not None:
nvec = nm.array(nvec)
if len(nvec.shape) == 0:
nvec = nm.repeat(nvec, len(cline) - 1)
else:
nvec = cline[1:] - cline[:-1]
nvec = nvec / nm.linalg.norm(nvec, axis=1)[:, None]
if twist is not None:
twist = nm.array(twist)
if len(twist.shape) == 0:
twist = nm.cumsum(nm.repeat(twist, len(cline)))
if scale is not None:
scale = nm.array(scale)
if len(scale.shape) == 0:
scale = nm.cumsum(nm.repeat(scale, len(cline)))
new_coors = [coors0 + cline[[0], :]]
new_conns, new_mat_ids, new_descs = [], [], []
mat_ids = mesh.cmesh.cell_groups
nnd = coors0.shape[0]
new_vgroups = nm.hstack([mesh.cmesh.vertex_groups] * len(cline))
for k in range(len(cline) - 1):
rotation = get_rotation(nvec0, nvec[k])
coors = rotation.apply(coors0)
if twist is not None:
rotation = Rotation.from_rotvec(nvec[k] * twist[k])
coors = rotation.apply(coors)
if scale is not None:
coors *= scale[k]
new_coors.append(coors + cline[[k + 1], :])
for desc in mesh.descs:
conns, cidxs = mesh.get_conn(desc, ret_cells=True)
econns = nm.hstack([conns + nnd * k, conns + nnd * (k + 1)])
emat_ids = mat_ids[cidxs]
if econns.shape[1] == 6:
if wedge_to_tetra: # wedges -> 3 x tetrahedron
econns1 = nm.empty((3* econns.shape[0], 4),
dtype=econns.dtype)
econns1[0::3] = econns[:, [0, 1, 2, 5]]
econns1[1::3] = econns[:, [0, 1, 5, 4]]
econns1[2::3] = econns[:, [0, 4, 5, 3]]
econns = econns1
emat_ids = nm.repeat(mat_ids, 3)
new_conns.append(econns)
new_mat_ids.append(emat_ids)
new_descs.append(f'3_{econns.shape[1]}')
out = Mesh.from_data(mesh.name,
nm.vstack(new_coors), new_vgroups,
*stack_descs(new_conns, new_mat_ids, new_descs))
return out
[docs]
def expand2d(mesh2d, dist, rep):
"""
Expand a 2D planar mesh into a 3D volume,
convert triangular/quad elements to tetrahedrons/hexahedrons.
Parameters
----------
mesh2d : Mesh
The 2D mesh.
dist : float
The elements size in the 3rd direction.
rep : int
The number of elements in the 3rd direction.
Returns
-------
mesh3d : Mesh
The 3D mesh.
"""
return extrude(mesh2d, (dist, rep))
[docs]
def get_unique_coor_map(coors, eps=1e-18):
tree = cKDTree(coors)
pairs = nm.array(list(tree.query_pairs(r=eps)))
if len(pairs) > 0:
n, row, col = len(coors), pairs[:, 0], pairs[:, 1]
graph = sps.csr_array((nm.ones(len(row)), (row, col)), shape=(n, n))
_, remap = connected_components(graph + graph.T)
return remap
[docs]
def merge_nodes(mesh, eps=1e-12):
"""
Merge duplicate mesh nodes.
Parameters
----------
mesh: Mesh
FE mesh.
eps: float
Tolerance for duplicity search.
Returns
-------
mesh: Mesh
FE mesh.
"""
remap = get_unique_coor_map(mesh.coors, eps=eps)
if remap is not None:
_, nidxs = nm.unique(remap, return_index=True)
new_conns, new_mat_ids = [], []
for desc in mesh.descs:
conns, cidxs = mesh.get_conn(desc, ret_cells=True)
new_mat_ids.append(mesh.cmesh.cell_groups[cidxs])
new_conns.append(remap[conns])
return Mesh.from_data(mesh.name, mesh.coors[nidxs],
mesh.cmesh.vertex_groups[nidxs],
new_conns, new_mat_ids, mesh.descs)
else:
return mesh
[docs]
def revolve(mesh, nphi, p, v=(1, 0, 0), phi_max=None, wedge_to_tetra=True):
"""
Create a solid 3D mesh from a given planar 2D mesh by revolving it
around the axis defined by a vector.
Parameters
----------
mesh: Mesh
2D planar FE mesh.
nphi: int
Number of elements in the circumferential direction
p: list, tuple, or numpy.ndarray
Coordinates of a point on the axis of rotation
v: list, tuple, or numpy.ndarray
Directional vector of the axis of rotation
phi_max: float or None
Angle of the revolution, if None, 360 degrees is used (closed loop)
wedge_to_tetra: bool
If True, convert wedge elements to tetrahedrons.
Returns
-------
out: Mesh
3D FE mesh.
"""
if phi_max is None:
phi_max = 360
phi = nm.linspace(0, 2*nm.pi, nphi + 1)
else:
phi = nm.linspace(0, nm.deg2rad(phi_max), nphi + 1)
p = nm.array(p)[None, :]
v = nm.array(v)
coors = mesh.coors
ccoor = nm.mean(coors, axis=0)
if coors.shape[1] == 2:
ccoor = nm.array([ccoor[0], ccoor[1], 0])
cline = nm.vstack([Rotation.from_rotvec(v * ph).apply(ccoor - p)
for ph in phi]) + p
nvec0 = nm.array([0, 0, 1])
nvec = nm.vstack([Rotation.from_rotvec(v * ph).apply(nvec0)
for ph in phi[1:]])
out = extrude(mesh, cline, nvec=nvec, wedge_to_tetra=wedge_to_tetra)
if phi_max == 360:
out = merge_nodes(out)
return out
[docs]
def mirror(mesh, p, v):
"""
Duplicate the mesh by mirroring it. The mirror plane is defined by a point
in the plane and by a normal vector to that plane.
Parameters
----------
mesh: Mesh
FE mesh.
p: list, tuple, or numpy.ndarray
Coordinates of a point at the mirror plane.
v: list, tuple, or numpy.ndarray
Normal vector of the mirror plane.
Returns
-------
out: Mesh
FE mesh.
"""
mirror_map = {
'2_3': nm.array([1, 0, 2]),
'2_4': nm.array([1, 0, 3, 2]),
'3_4': nm.array([1, 0, 2, 3]),
'3_8': nm.array([1, 0, 3, 2, 5, 4, 7, 6]),
}
coors = mesh.coors
if coors.shape[1] == 2:
coors = nm.hstack([coors, coors[:, [0]] * 0])
p = nm.asarray(p)[None, :]
v = nm.asarray(v)
v = v / nm.linalg.norm(v)
vc = coors - p
dist = nm.dot(vc, v)
new_coors = nm.vstack([
coors,
coors - 2 * dist[:, None] * v
])
new_vgroups = nm.hstack([mesh.cmesh.vertex_groups] * 2)
new_conns, new_mat_ids, new_descs = [], [], []
i1 = coors.shape[0]
for desc in mesh.descs:
conns, cidxs = mesh.get_conn(desc, ret_cells=True)
new_mat_ids.append(nm.hstack([mesh.cmesh.cell_groups[cidxs]] * 2))
new_descs.append(desc)
new_conns.append(nm.vstack([conns, conns[:, mirror_map[desc]] + i1]))
out = Mesh.from_data(mesh.name, new_coors, new_vgroups,
new_conns, new_mat_ids, mesh.descs)
return merge_nodes(out)