# Source code for sfepy.mesh.mesh_tools

```from sfepy.discrete.fem import Mesh, FEDomain
import scipy.sparse as sps
import numpy as nm
from sfepy.base.compat import factorial
from sfepy.base.base import output
from sfepy.discrete.common.extmods.cmesh import (create_mesh_graph,
graph_components)

[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.identity(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_matrix((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_matrix((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 expand2d(mesh2d, dist, rep):
"""
Expand 2D planar mesh into 3D volume,

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.
"""
if len(mesh2d.descs) > 1:
raise ValueError('More than one cell type (%s). Not supported!'
% ', '.join(mesh2d.descs))

nel = mesh2d.n_el
nnd = mesh2d.n_nod
et = mesh2d.descs[0]
coors = mesh2d.coors
conn = mesh2d.get_conn(et)

zcoor = nm.arange(rep + 1) * dist
coors3d = nm.hstack([nm.tile(coors, (rep + 1, 1)),
nm.tile(zcoor, (nnd,1)).T.flatten()[:,nm.newaxis]])
ngroups = nm.tile(mesh2d.cmesh.vertex_groups, (rep + 1,))

if et == '2_4':
descs3d = '3_8'
conn3d = nm.zeros((nel * rep, 8), dtype=nm.int32)
mats3d = nm.tile(mesh2d.cmesh.cell_groups, (1, rep)).squeeze()

elif et == '2_3':
descs3d = '3_4'
conn3d = nm.zeros((3 * nel * rep, 4), dtype=nm.int32)
mats3d = nm.tile(mesh2d.cmesh.cell_groups, (1, 3 * rep)).squeeze()

for ii in range(rep):
bgn0 = nnd * ii
bgn1 = bgn0 + nnd
if et == '2_4':
bge0 = nel * ii
bge1 = bge0 + nel
conn3d[bge0:bge1,:4] = conn + bgn0
conn3d[bge0:bge1,4:] = conn + bgn1

elif et == '2_3':
# 0 1 2 5
bge0 = 3 * nel * ii
bge1 = bge0 + nel
conn3d[bge0:bge1,:] = nm.array([conn[:,0] + bgn0,
conn[:,1] + bgn0,
conn[:,2] + bgn0,
conn[:,2] + bgn1]).T
# 0 1 5 4
bge0 += nel
bge1 += nel
conn3d[bge0:bge1,:] = nm.array([conn[:,0] + bgn0,
conn[:,1] + bgn0,
conn[:,2] + bgn1,
conn[:,1] + bgn1]).T
# 0 4 5 3
bge0 += nel
bge1 += nel
conn3d[bge0:bge1,:] = nm.array([conn[:,0] + bgn0,
conn[:,1] + bgn1,
conn[:,2] + bgn1,
conn[:,0] + bgn1]).T

mesh3d = Mesh.from_data('mesh', coors3d, ngroups, [conn3d],
[mats3d], [descs3d])

return mesh3d

[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_matrix((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

[docs]
def extract_edges(mesh, eps=1e-16):
"""
Extract outline edges of a given mesh.
The outline edge is an edge for which norm(nvec_1 - nvec_2) < eps,
where nvec_1 and nvec_2 are the normal vectors of the incident facets.

Parameters
----------
mesh : Mesh
The 3D or 2D mesh.
eps : float
The tolerance parameter of the outline edge searching algorithm.

Returns
-------
mesh_out : tuple
The data of the outline mesh, Mesh.from_data() format, i.e.
(coors, ngroups, ed_conns, mat_ids, descs).
"""
domain = FEDomain('domain', mesh)
cmesh = domain.cmesh

output('Mesh - dimension: %d, vertices: %d, elements: %d'
% (mesh.dim, mesh.n_nod, mesh.n_el))

if mesh.dim == 2:
oedges = cmesh.get_surface_facets()
mesh_coors = nm.hstack([cmesh.coors,
nm.zeros((cmesh.coors.shape[0], 1))])

elif mesh.dim == 3:
cmesh.setup_connectivity(1, 2)
cmesh.setup_connectivity(3, 2)

sfaces = cmesh.get_surface_facets()
_, idxs = nm.unique(cmesh.get_conn(3, 2).indices, return_index=True)

normals = cmesh.get_facet_normals()[idxs, :]
se_map, se_off = cmesh.get_incident(1, sfaces, 2, ret_offsets=True)
sedges = nm.unique(se_map)
n_se = sedges.shape[0]

# remap surface edges to continuous range
se_remap = -nm.ones(sedges.max() + 1)
se_remap[sedges] = nm.arange(n_se)
se_map0 = se_remap[se_map]

# surface face/edge connectivity matrix (n_surf x n_edge)
n_ef = nm.diff(se_off)[0]  # = 2
n_sf = se_map.shape[0] // n_ef
row = nm.repeat(nm.arange(n_sf), n_ef)
sf2e = sps.coo_matrix((nm.ones((n_sf * n_ef,), dtype=bool),
(row , se_map0)), shape=(n_sf, n_se))
# edge to face map (n_edge x 2)
se2f = sf2e.tocsc().indices.reshape((sedges.shape[0], 2))

snormals = normals[sfaces]
err = nm.linalg.norm(snormals[se2f[:, 0]] - snormals[se2f[:, 1]],
axis=1)
oedges = sedges[nm.where(err > eps)[0]]
mesh_coors = cmesh.coors

else:
raise NotImplementedError

# save outline mesh
if oedges.shape[0] > 0:
ec_idxs = nm.unique(cmesh.get_incident(0, oedges, 1))
ed_coors = mesh_coors[ec_idxs, :]
ngroups = nm.zeros((ed_coors.shape[0],), dtype=nm.int16)

aux = cmesh.get_conn(1, 0).indices

ed_conns = aux.reshape((aux.shape[0] // 2, 2))[oedges, :]
ec_remap = -nm.ones((ec_idxs.max() + 1, ), dtype=nm.int64)
ec_remap[ec_idxs] = nm.arange(ec_idxs.shape[0])
ed_conns = ec_remap[ed_conns]
mat_ids = nm.ones((ed_conns.shape[0],), dtype=nm.int16)

mesh_out = ed_coors, ngroups, [ed_conns], [mat_ids], ['3_2']

return mesh_out

else:
raise ValueError('no outline edges found (eps=%e)!' % eps)

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):
nnz, prow, icol = create_mesh_graph(n_nod, n_nod, len(surf_faces),
surf_faces, surf_faces)
data = nm.empty((nnz,), dtype=nm.int32)
data.fill(2)
return sps.csr_matrix((data, icol, prow), (n_nod, n_nod))

[docs]
def surface_components(gr_s, surf_faces):
"""
Determine surface components given surface mesh connectivity graph.
"""
n_nod = gr_s.shape[0]
n_comp, flag = graph_components(n_nod, gr_s.indptr, gr_s.indices)

comps = []
for ii, face in enumerate(surf_faces):
comp = flag[face[:,0]]
comps.append(comp)

return n_comp, comps

```