Thanks, Carlo, with your help I managed to get what I wanted!
What helped me was the mention of the
method!
I believe, your
Code: Select all
bot_face = Part.Face(close_spline(sections[0]))
top_face = Part.Face(close_spline(sections[-1]))
works only for the special case of flat bottom and top sections and not in the general case of arbitrary section in the 3d space (as mentioned in my very first email).
Anyway, here is my code for anyone, who has the same (or a similar) problem. This code can even export the solid in four different file formats (see commented section near the bottom of the file):
Code: Select all
#!/usr/bin/env python3
import math
from OCC.Core import gp
from OCC.Core import TColgp
from OCC.Core import GeomAPI
from OCC.Core import GeomFill
from OCC.Core import GeomAbs
from OCC.Core import BRepFill
from OCC.Core import BRepBuilderAPI
from OCC.Core import GProp
from OCC.Core import BRepGProp
from OCC.Display.SimpleGui import init_display
NO_X = 10
NO_Z = 10
X = 1
Y = 0.1
Z = 2
def ellipse(x, a, b):
"""Calculate the y-coordinate of an ellipse with the one apex at (0,0) and one at (2a,0)."""
return b/a * math.sqrt(a**2 - (x-a)**2)
def parabola(x, a, b):
"""Calculate the y-coordinate of an parabola through the points (0,0), (a,b) and (2a,0)."""
return b * (1 - ((x-a)/a)**2)
def points():
"""Generate the points of our surface."""
array = TColgp.TColgp_Array2OfPnt(1, 2*NO_X-1, 1, NO_Z)
for j in range(NO_Z):
z = (j+10)/NO_Z * Z
b = Z/10 - 0.02*z
# First side (back)
for i in range(1, NO_X):
x = i/NO_X * X
y = parabola(x, X/2, 2*b) + ellipse(x, X/2, b)
point = gp.gp_Pnt(x, y, z)
array.SetValue(i, j+1, point)
# Second side (face)
k = i
for i in range(NO_X, 0, -1):
x = i/NO_X * X
y = parabola(x, X/2, 2*b) - ellipse(x, X/2, b)
point = gp.gp_Pnt(x, y, z)
k += 1
array.SetValue(k, j+1, point)
return array
# Fit a surface to our points
pts = points()
surface = GeomAPI.GeomAPI_PointsToBSplineSurface()
surface.Interpolate(pts, False)
surface = surface.Surface()
# Boundary lines to our surface
u_back, u_face, v_bottom, v_top = surface.Bounds()
bl_back = surface.UIso(u_back)
bl_face = surface.UIso(u_face)
bl_bottom = surface.VIso(v_bottom)
bl_top = surface.VIso(v_top)
# Fit a surface to the open aft end of our surface
surface_aft = GeomFill.geomfill_Surface(bl_back, bl_face)
# Boundary lines to the aft surface
u_bottom2, u_top2, v_face2, v_back2 = surface_aft.Bounds()
bl_bottom2 = surface_aft.UIso(u_bottom2)
bl_top2 = surface_aft.UIso(u_top2)
# Fit a surface to the open bottom of our surface
edge_bottom = BRepBuilderAPI.BRepBuilderAPI_MakeEdge(bl_bottom)
edge_bottom2 = BRepBuilderAPI.BRepBuilderAPI_MakeEdge(bl_bottom2)
filling_bottom = BRepFill.BRepFill_Filling()
filling_bottom.Add(edge_bottom.Edge(), GeomAbs.GeomAbs_C0)
filling_bottom.Add(edge_bottom2.Edge(), GeomAbs.GeomAbs_C0)
filling_bottom.Build()
# Fit a surface to the open top of our surface
edge_top = BRepBuilderAPI.BRepBuilderAPI_MakeEdge(bl_top)
edge_top2 = BRepBuilderAPI.BRepBuilderAPI_MakeEdge(bl_top2)
filling_top = BRepFill.BRepFill_Filling()
filling_top.Add(edge_top.Edge(), GeomAbs.GeomAbs_C0)
filling_top.Add(edge_top2.Edge(), GeomAbs.GeomAbs_C0)
filling_top.Build()
# Build a shell from all surfaces by sewing them together
shell = BRepBuilderAPI.BRepBuilderAPI_Sewing()
face_surface = BRepBuilderAPI.BRepBuilderAPI_MakeFace(surface, 0.001)
shell.Add(face_surface.Face())
face_aft = BRepBuilderAPI.BRepBuilderAPI_MakeFace(surface_aft, 0.001)
shell.Add(face_aft.Face())
shell.Add(filling_bottom.Face())
shell.Add(filling_top.Face())
shell.Perform()
# Reverse the shell, so that we get a positive volume
shell_shape = shell.SewedShape()
shell_shape.Complement()
# Build a solid
solid = BRepBuilderAPI.BRepBuilderAPI_MakeSolid(shell_shape)
solid.Build()
# Get the volume (which is the mass, if no density is defined)
gprop = GProp.GProp_GProps()
BRepGProp.brepgprop_VolumeProperties(solid.Shape(), gprop)
assert gprop.Mass() > 0, 'Shell is oriented the wrong way!'
print(f'Volume = {gprop.Mass()}')
print()
# # Export to .brep file
# from OCC.Core import BRepTools
# BRepTools.breptools_Write(solid, 'solid.brep')
# # Export other formats
# from OCC.Extend import DataExchange
# # Write STEP file
# DataExchange.write_step_file(solid.Shape(), 'solid.step')
# # Write STL file
# DataExchange.write_stl_file(
# solid.Shape(), 'solid.stl', linear_deflection=X/NO_X/10)
# # Write IGES file
# DataExchange.write_iges_file(solid.Shape(), 'solid.iges')
# print()
# Display the surfaces
display, start_display, _, _ = init_display()
display.DisplayShape(solid.Shape())
display.DisplayShape(gprop.CentreOfMass())
display.FitAll()
start_display()
Maybe all this is easier using FreeCAD's API, but I wanted to keep the number (and size) of the dependencies small.
