#!/usr/bin/env python
# 05.10.2005, c
"""
Given a mesh file, this script extracts its surface and prints it to stdout in
form of a list where each row is [element, face, component]. A component
corresponds to a contiguous surface region - for example, a cubical mesh with a
spherical hole has two surface components. Two surface faces sharing a single
node belong to one component.
With '-m' option, a mesh of the surface is created and saved in
'<original path>/surf_<original mesh file name>.mesh'.
"""
import sys
sys.path.append('.')
from optparse import OptionParser
import numpy as nm
import scipy.sparse as sp
import sfepy
from sfepy.base.base import output
from sfepy.base.ioutils import edit_filename
from sfepy.discrete.fem import Mesh, FEDomain
from sfepy.discrete.common.extmods.cmesh import (create_mesh_graph,
graph_components)
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 sp.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
usage = """%prog [options] filename_in|- filename_out|-
'-' is for stdin, stdout
""" + __doc__.rstrip()
[docs]def main():
parser = OptionParser(usage=usage, version="%prog " + sfepy.__version__)
parser.add_option("-m", "--mesh",
action="store_true", dest="save_mesh",
default=False,
help="save surface mesh")
parser.add_option("-n", "--no-surface",
action="store_true", dest="no_surface",
default=False,
help="do not output surface [default: %default]")
(options, args) = parser.parse_args()
if (len(args) == 2):
filename_in = args[0];
filename_out = args[1];
else:
parser.print_help(),
return
if (filename_in == '-'):
file_in = sys.stdin
else:
file_in = open(filename_in, "r");
mesh = Mesh.from_file(filename_in)
if (filename_in != '-'):
file_in.close()
domain = FEDomain('domain', mesh)
if options.save_mesh:
region = domain.create_region('surf', 'vertices of surface', 'facet')
surf_mesh = Mesh.from_region(region, mesh,
localize=True, is_surface=True)
aux = edit_filename(filename_in, prefix='surf_', new_ext='.mesh')
surf_mesh.write(aux, io='auto')
if domain.has_faces():
domain.fix_element_orientation()
lst, surf_faces = get_surface_faces(domain)
if options.no_surface:
return
gr_s = surface_graph(surf_faces, mesh.n_nod)
n_comp, comps = surface_components(gr_s, surf_faces)
output('number of surface components:', n_comp)
ccs, comps = comps, nm.zeros((0,1), nm.int32)
for cc in ccs:
comps = nm.concatenate((comps, cc[:,nm.newaxis]), 0)
out = nm.concatenate((lst, comps), 1)
if (filename_out == '-'):
file_out = sys.stdout
else:
file_out = open(filename_out, "w");
for row in out:
file_out.write('%d %d %d\n' % (row[0], row[1], row[2]))
if (filename_out != '-'):
file_out.close()
if __name__=='__main__':
main()