/src/examples/Level1/Geometry/geomplate.py

  1#!/usr/bin/env python
  2
  3# TODO:
  4# * need examples where the tangency to constraining faces is respected
  5
  6##Copyright 2009-2011 Jelle Ferina (jelleferinga@gmail.com)
  7##
  8##This file is part of pythonOCC.
  9##
 10##pythonOCC is free software: you can redistribute it and/or modify
 11##it under the terms of the GNU Lesser General Public License as published by
 12##the Free Software Foundation, either version 3 of the License, or
 13##(at your option) any later version.
 14##
 15##pythonOCC is distributed in the hope that it will be useful,
 16##but WITHOUT ANY WARRANTY; without even the implied warranty of
 17##MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 18##GNU Lesser General Public License for more details.
 19##
 20##You should have received a copy of the GNU Lesser General Public License
 21##along with pythonOCC.  If not, see <http://www.gnu.org/licenses/>.
 22import os
 23from OCC.Utils.Construct import make_closed_polygon, make_n_sided, make_vertex, make_face, make_wire
 24
 25from OCC.gp import *
 26from OCC.BRepBuilderAPI import *
 27
 28from OCC.Utils.Topology import WireExplorer, Topo
 29from OCC.BRepAdaptor import *
 30from OCC.BRep import *
 31from OCC.ShapeAnalysis import *
 32from OCC.GeomLProp import *
 33
 34import types, sys, time
 35
 36from OCC.Display.SimpleGui import init_display
 37from OCC.Utils.DataExchange.IGES import IGESImporter
 38from OCC.BRepFill import *
 39from OCC.GeomPlate import *
 40from OCC.GEOMAlgo import *
 41
 42display, start_display, add_menu, add_function_to_menu = init_display()
 43
 44try:
 45    from scipy import arange
 46    from scipy.optimize import fsolve
 47except ImportError:
 48    print 'scipy not installed, will not be able to run the solve radius example'
 49
 50
 51def geom_plate(event=None):
 52    display.EraseAll()
 53    p1,p2,p3,p4,p5 = gp_Pnt(0,0,0),gp_Pnt(0,10,0),gp_Pnt(0,10,10),gp_Pnt(0,0,10),gp_Pnt(5,5,5)
 54    poly = make_closed_polygon([p1,p2,p3,p4])
 55    edges = [i for i in Topo(poly).edges()]
 56    face = make_n_sided(edges, [p5])
 57    display.DisplayShape(edges)
 58    display.DisplayShape(make_vertex(p5))
 59    display.DisplayShape(face, update=True)
 60
 61#===============================================================================
 62# Find a surface such that the radius at the vertex is n
 63#===============================================================================
 64
 65def build_plate(polygon, points):
 66    '''
 67    build a surface from a constraining polygon(s) and point(s)
 68    @param polygon:     list of polygons ( TopoDS_Shape)
 69    @param points:      list of points ( gp_Pnt ) 
 70    '''
 71    # plate surface
 72    bpSrf = GeomPlate_BuildPlateSurface(3,15,2)
 73    
 74    # add curve constraints
 75    for poly in polygon:
 76        for edg in WireExplorer(poly).ordered_edges():
 77            c = BRepAdaptor_HCurve()
 78            c.ChangeCurve().Initialize(edg)
 79            constraint = BRepFill_CurveConstraint(c.GetHandle(), 0)
 80            bpSrf.Add(constraint.GetHandle())
 81        
 82    # add point constraint
 83    for pt in points:
 84        bpSrf.Add(GeomPlate_PointConstraint(pt, 0).GetHandle())
 85        bpSrf.Perform()
 86    
 87    maxSeg, maxDeg, critOrder = 9,8,0
 88    tol  = 1e-4
 89    dmax = max([tol,10*bpSrf.G0Error()])
 90    
 91    srf = bpSrf.Surface()
 92    plate = GeomPlate_MakeApprox(srf, tol, maxSeg, maxDeg, dmax, critOrder)
 93    
 94    uMin, uMax, vMin, vMax = srf.GetObject().Bounds()
 95    
 96    return make_face(plate.Surface(), uMin, uMax, vMin, vMax)
 97
 98def radius_at_uv(face, u, v):
 99    '''
100    returns the mean radius at a u,v coordinate
101    @param face:    surface input
102    @param u,v:     u,v coordinate
103    '''
104    h_srf = BRep_Tool().Surface(face)
105    uv_domain = GeomLProp_SurfaceTool().Bounds(h_srf)
106    curvature = GeomLProp_SLProps(h_srf, u, v, 1, 1e-6)
107    try:
108        _crv_min = 1./curvature.MinCurvature()
109    except ZeroDivisionError:
110        _crv_min = 0.
111    
112    try:
113        _crv_max = 1./curvature.MaxCurvature()
114    except ZeroDivisionError:
115        _crv_max = 0.
116    return abs((_crv_min+_crv_max)/2.)
117
118def uv_from_projected_point_on_face(face, pt):
119    '''
120    returns the uv coordinate from a projected point on a face
121    '''
122    srf = BRep_Tool().Surface(face)
123    sas = ShapeAnalysis_Surface(srf)
124    uv  = sas.ValueOfUV(pt, 1e-2)
125    print 'distance',sas.Value(uv).Distance(pt)    
126    return uv.Coord()
127
128
129class RadiusConstrainedSurface():
130    '''
131    returns a surface that has `radius` at `pt`
132    '''
133    def __init__(self, display, poly, pnt, targetRadius):
134        self.display = display
135        self.targetRadius = targetRadius
136        self.poly = poly
137        self.pnt = pnt
138        self.plate = self.build_surface()
139        
140    def build_surface(self):
141        '''
142        builds and renders the plate
143        '''
144        self.plate = build_plate([self.poly], [self.pnt])
145        self.display.EraseAll()
146        self.display.DisplayShape(self.plate)
147        vert = make_vertex(self.pnt)
148        self.display.DisplayShape(vert, update=True)
149    
150    def radius(self, z):
151        '''
152        sets the height of the point constraining the plate, returns the radius at this point
153        '''
154        if isinstance(z, types.FloatType):
155            self.pnt.SetX(z)
156        else:
157            self.pnt.SetX(float(z[0]))
158        self.build_surface()
159        uv = uv_from_projected_point_on_face(self.plate, self.pnt)
160        radius = radius_at_uv(self.plate, uv[0], uv[1])
161        print 'z:',z, 'radius:', radius
162        self.curr_radius = radius
163        return self.targetRadius-abs(radius)
164
165    def solve(self):
166        fsolve(self.radius, 1, maxfev=1000)
167        return self.plate
168
169
170def solve_radius(event=None):
171    display.EraseAll()
172    p1,p2,p3,p4,p5 = gp_Pnt(0,0,0),gp_Pnt(0,10,0),gp_Pnt(0,10,10),gp_Pnt(0,0,10),gp_Pnt(5,5,5)
173    poly = make_closed_polygon([p1,p2,p3,p4])
174    for i in arange(0.1,3.,0.2).tolist():
175        rcs = RadiusConstrainedSurface(display, poly, p5, i )
176        face = rcs.solve()
177        print 'Goal: %s radius: %s' % ( i, rcs.curr_radius )
178        time.sleep(0.5)
179
180#===============================================================================
181# 
182#===============================================================================
183
184def build_geom_plate(edges):
185    bpSrf = GeomPlate_BuildPlateSurface(3,9,12)
186    
187    # add curve constraints
188    for edg in edges:
189        c = BRepAdaptor_HCurve()
190        print 'edge:',edg
191        c.ChangeCurve().Initialize(edg)
192        constraint = BRepFill_CurveConstraint(c.GetHandle(), 0)
193        bpSrf.Add(constraint.GetHandle())
194    
195    # add point constraint
196    try:
197        bpSrf.Perform()
198    except RuntimeError:
199        print 'failed to build the geom plate surface '
200    
201    maxSeg, maxDeg, critOrder = 9,8,0
202    tol  = 1e-4
203    dmax = max([tol,10*bpSrf.G0Error()])
204    
205    srf = bpSrf.Surface()
206#    plate = GeomPlate_MakeApprox(srf, tol, maxSeg, maxDeg, dmax, critOrder)
207    plate = GeomPlate_MakeApprox(srf, 1e-04, 100, 9, 1e-03, 0)
208    
209    uMin, uMax, vMin, vMax = srf.GetObject().Bounds()
210    face = make_face(plate.Surface(), uMin, uMax, vMin, vMax)
211    return face
212
213
214def build_curve_network(event=None):
215    '''
216    mimic the curve network surfacing command from rhino
217    '''
218    print 'Importing IGES file...',
219    pth = os.path.dirname(os.path.abspath(__file__))
220    pth = os.path.abspath(os.path.join(pth, '../../data/IGES/curve_geom_plate.igs'))
221    iges = IGESImporter(pth)
222    iges.read_file()
223    print 'done.'
224    # GetShapes returns 36 TopoDS_Shape, while the TopoDS_Compound contains
225    # just the 6 curves I made in rhino... hmmm...
226    print 'Building geomplate...',
227    topo = Topo(iges.get_compound())
228    edges_list = list(topo.edges())
229    face = build_geom_plate(edges_list)
230    print 'done.'
231    display.EraseAll()
232    display.DisplayShape(edges_list)
233   
234    print 'Cutting out of edges...',
235    # Make a wire from outer edges
236    _edges = [edges_list[2], edges_list[3], edges_list[4], edges_list[5]]
237    outer_wire = make_wire(_edges)
238
239    ################ --- ################
240
241#    for i in edges_list:
242#        if i not in _edges:
243#            _edges.append(i)
244#
245#    import ipdb; ipdb.set_trace()
246#    f = make_n_sided(_edges)
247
248
249    # TODO: cut out the shape from the edges in _edges
250
251#    if outer_wire:
252#        splitter = GEOMAlgo_Splitter()
253#        splitter.AddTool(outer_wire)
254#        splitter.AddShape(face)
255#        splitter.Perform()
256#        shp = splitter.Shape()
257#        display.DisplayShape(shp)
258
259def exit(event=None):
260    sys.exit()
261
262if __name__ == "__main__":
263    build_curve_network()
264    add_menu('geom plate')
265    add_function_to_menu('geom plate', geom_plate)
266    add_function_to_menu('geom plate', solve_radius)
267    add_function_to_menu('geom plate', build_curve_network)
268    add_function_to_menu('geom plate', exit)
269    start_display()