Source code for plugins.pyformex_gts

#
##
##  This file is part of pyFormex 2.0  (Mon Sep 14 12:29:05 CEST 2020)
##  pyFormex is a tool for generating, manipulating and transforming 3D
##  geometrical models by sequences of mathematical operations.
##  Home page: http://pyformex.org
##  Project page:  http://savannah.nongnu.org/projects/pyformex/
##  Copyright 2004-2020 (C) Benedict Verhegghe (benedict.verhegghe@ugent.be)
##  Distributed under the GNU General Public License version 3 or later.
##
##  This program is free software: you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation, either version 3 of the License, or
##  (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program.  If not, see http://www.gnu.org/licenses/.
##
"""Operations on triangulated surfaces using GTS functions.

This module provides access to GTS from insisde pyFormex.
"""

import pyformex as pf
from pyformex.lazy import *
from pyformex import utils
from pyformex.multi import multitask

# No warning when in sphinx/docmodule mode
if not pf.sphinx:
    if not utils.External.has('gts-bin'):
        utils.warn("error_no_gts_bin")
    if not utils.External.has('gts-extra'):
        utils.warn("error_no_gts_extra")


#
# gts commands used:
#   in Debian package: stl2gts gts2stl gtscheck
#   not in Debian package: gtssplit gtscoarsen gtsrefine gtssmooth gtsinside
#


[docs]def read_gts_intersectioncurve(fn): """Read the intersection curve of a boolean operation Returns ------- Mesh A Mesh of eltype Line2 containing the line segments on the intersection curve. """ import re from pyformex.formex import Formex RE = re.compile("^VECT 1 2 0 2 0 (?P<data>.*)$") r = [] for line in open(fn, 'r'): m = RE.match(line) if m: r.append(m.group('data')) nelems = len(r) x = fromstring('\n'.join(r), sep=' ').reshape(-1, 2, 3) return Formex(x).toMesh()
[docs]def boolean(self, surf, op, check=False, verbose=False): """Perform a boolean operation with another surface. Boolean operations between surfaces are a basic operation in free surface modeling. Both surfaces should be closed orientable non-intersecting manifolds. Use the :meth:`check` method to find out. The boolean operations are set operations on the enclosed volumes: union('+'), difference('-') or intersection('*'). Parameters ---------- surf: TriSurface Another TriSurface that is a closed manifold surface. op: '+', '-' or '*' The boolean operation to perform: union('+'), difference('-') or intersection('*'). check: bool If True, a check is done that the surfaces are not self-intersecting; if one of them is, the set of self-intersecting faces is written (as a GtsSurface) on standard output verbose: bool If True, print statistics about the surface. Returns ------- TriSurface A closed manifold TriSurface that is the volume union, difference or intersection of self with surf. Note ---- This method uses the external command 'gtsset' and will not run if it is not installed (available from pyformex/extras). """ return self.gtsset(surf, op, filt = '', ext='.gts', check=check, verbose=verbose).fuse().compact()
[docs]def intersection(self, surf, check=False, verbose=False): """Return the intersection curve(s) of two surfaces. Boolean operations between surfaces are a basic operation in free surface modeling. Both surfaces should be closed orientable non-intersecting manifolds. Use the :meth:`check` method to find out. Parameters ---------- surf: TriSurface A closed manifold surface. check: bool, optional If True, a check is made that the surfaces are not self-intersecting; if one of them is, the set of self-intersecting faces is written (as a GtsSurface) on standard output verbose: bool, optional If True, statistics about the surface are printed on stdout. Returns ------- Mesh A Mesh with eltype Line2 holding all the line segments of the intersection curve(s). """ return self.gtsset(surf, op='*', ext='.list', curve=True, check=check, verbose=verbose)
def gtsset(self, surf, op, filt='', ext='', curve=False, check=False, verbose=False): """_Perform the boolean/intersection methods. See the boolean/intersection methods for more info. Parameters not explained there: - `filt`: a filter command to be executed on the gtsset output - `ext`: extension of the result file - curve: if True, an intersection curve is computed, else the surface. Returns the resulting TriSurface (curve=False), or Mesh with eltype Line2 (curve=True), or None if the input surfaces do not intersect. """ # import here to avoid circular import from pyformex.trisurface import TriSurface op = { '+': 'union', '-': 'diff', '*': 'inter', 'a': 'all', }[op] options = '' if curve: options += '-i' if check: options += ' -s' if verbose: options += ' -v' with utils.TempDir() as tmpdir: tmp = tmpdir / 'surface1.gts' tmp1 = tmpdir / 'surface2.gts' tmp2 = tmpdir / ('result'+ext) #print("Writing temp file %s" % tmp) self.write(tmp, 'gts') #print("Writing temp file %s" % tmp1) surf.write(tmp1, 'gts') print("Performing boolean operation") cmd = f"cd {tmpdir} && gtsset {options} {op} {tmp} {tmp1} {filt}" P = utils.system(cmd, stdout=open(tmp2, 'w'), shell=True, verbose=verbose) if P.returncode or verbose: print(P.stdout) if P.returncode: print(P.stderr) return None if curve: res = read_gts_intersectioncurve(tmp2) elif op=='all': res = {} for k in ['s1in2', 's1out2', 's2in1', 's2out1']: res[k] = TriSurface.read(tmpdir / (k+'.gts')) else: res = TriSurface.read(tmp2) return res # TODO: change pts to a Coords def gtsinside(self, pts, dir=0): """_Test whether points are inside the surface. pts is a plex-1 Formex. dir is the shooting direction. Returns a list of point numbers that are inside. This is not intended to be used directly. Use inside instead """ S = self.rollAxes(dir) P = pts.rollAxes(dir) with utils.TempDir() as tmpdir: tmp = tmpdir / 'surface.gts' tmp1 = tmpdir / 'points.dta' tmp2 = tmpdir / 'result.out' #print("Writing temp file %s" % tmp) S.write(tmp, 'gts') #print("Writing temp file %s" % tmp1) with open(tmp1, 'w') as f: P.coords.tofile(f, sep=' ') f.write('\n') #print("Performing inside testing") cmd = "gtsinside %s %s" % (tmp, tmp1) P = utils.command(cmd, stdout=open(tmp2, 'w')) if P.returncode: #print("An error occurred during the testing.\nSee file %s for more details." % tmp2) print(P.stdout) return None #print("Reading results from %s" % tmp2) ind = fromfile(tmp2, sep=' ', dtype=Int) return ind
[docs]def inside(self, pts, atol='auto', multi=True): """Test which of the points pts are inside the surface. Parameters: - `pts`: a (usually 1-plex) Formex or a data structure that can be used to initialize a Formex. Returns an integer array with the indices of the points that are inside the surface. The indices refer to the onedimensional list of points as obtained from pts.points(). """ from pyformex.formex import Formex if not isinstance(pts, Formex): pts = Formex(pts) if atol == 'auto': atol = pts.dsize()*0.001 # determine bbox of common space of surface and points bb = bboxIntersection(self, pts) if (bb[0] > bb[1]).any(): # No bbox intersection: no points inside return array([], dtype=Int) # Limit the points to the common part # Add point numbers as property, to allow return of original numbers pts.setProp(arange(pts.nelems())) pts = pts.clip(pts.testBbox(bb, atol=atol)) ins = zeros((3, pts.nelems()), dtype=bool) if pf.scriptMode == 'script' or not multi: for i in range(3): dirs = roll(arange(3), -i)[1:] # clip the surface perpendicular to the shooting direction # !! gtsinside seems to fail sometimes when using clipping S = self # .clip(testBbox(self,bb,dirs=dirs,atol=atol),compact=False) # find inside points shooting in direction i ok = gtsinside(S, pts, dir=i) ins[i, ok] = True else: tasks = [(gtsinside, (self, pts, i)) for i in range(3)] ind = multitask(tasks, 3) for i in range(3): ins[i, ind[i]] = True ok = where(ins.sum(axis=0) > 1)[0] return pts.prop[ok]
# End