Source code for extract_surface

#!/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'.
from __future__ import absolute_import
import sys
from six.moves import range
from argparse import ArgumentParser, RawDescriptionHelpFormatter

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,

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
[docs]def main(): parser = ArgumentParser(description=__doc__, formatter_class=RawDescriptionHelpFormatter) parser.add_argument("--version", action="version", version="%(prog)s " + sfepy.__version__) parser.add_argument("-m", "--mesh", action="store_true", dest="save_mesh", default=False, help="save surface mesh") parser.add_argument("-n", "--no-surface", action="store_true", dest="no_surface", default=False, help="do not output surface [default: %(default)s]") parser.add_argument('filename_in', help="'-' is for stdin") parser.add_argument('filename_out', help="'-' is for stdout") options = parser.parse_args() filename_in = options.filename_in filename_out = options.filename_out 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()