Skip to content

Commit

Permalink
Add: function for closed trajectories
Browse files Browse the repository at this point in the history
  • Loading branch information
@OberGue
OberGue committed Jul 2, 2024
1 parent 4387404 commit 6737fb8
Show file tree
Hide file tree
Showing 2 changed files with 179 additions and 14 deletions.
131 changes: 121 additions & 10 deletions examples/show_swept.py
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
import sys

import gustaf as gus
import numpy as np

Expand All @@ -6,25 +8,87 @@
if __name__ == "__main__":

### TRAJECTORY ###

# arbitrary trajectory
# dict_trajectory = {
# "degrees": [2],
# "knot_vectors": [[0.0, 0.0, 0.0, 0.333, 0.666, 1.0, 1.0, 1.0]],
# "control_points": np.array(
# [
# [0.0, 0.0, 0.0],
# [0.0, 0.0, 5.0],
# [10.0, 5.0, 0.0],
# [15.0, 0.0, -5.0],
# [20.0, 0.0, 0.0],
# ]
# ),
# }

# 2D questionmark
# dict_trajectory = {
# "degrees": [3],
# "knot_vectors": [[0.0, 0.0, 0.0, 0.0, 0.2, 0.4,
# 0.6, 0.8, 0.9, 1.0, 1.0, 1.0, 1.0]],
# "control_points": np.array([
# [0.5, 0], # Startpunkt
# [0.5, 2],
# [1.0, 3],
# [2.0, 4],
# [2.15, 5],
# [1.8, 5.9],
# [1.0, 6.2],
# [-0.25, 6],
# [-0.5, 5],])
# }
# init trajectory as bspline

# closed 3D questionmark
dict_trajectory = {
"degrees": [2],
"knot_vectors": [[0.0, 0.0, 0.0, 0.333, 0.666, 1.0, 1.0, 1.0]],
"degrees": [3],
"knot_vectors": [
[
0.0,
0.0,
0.0,
0.0,
0.1,
0.2,
0.3,
0.4,
0.5,
0.6,
0.8,
0.9,
1.0,
1.0,
1.0,
1.0,
]
],
"control_points": np.array(
[
[0.0, 0.0, 0.0],
[0.0, 0.0, 5.0],
[10.0, 5.0, 0.0],
[15.0, 0.0, -5.0],
[20.0, 0.0, 0.0],
[0.5, 0, 0],
[0.5, 2, 0.3],
[1.0, 3, 0.1],
[2.0, 4, -0.1],
[2.15, 5, -0.2],
[1.8, 5.9, -0.4],
[1.0, 6.2, -0.3],
[-0.25, 6, -0.1],
[-0.5, 5.0, 0.1],
[-2.0, 4.0, 0.2],
[-1, 3.0, 0.1],
[0.5, 0.0, 0.0],
]
),
}

# init trajectory as bspline
trajectory = splinepy.BSpline(**dict_trajectory)

# alternatively, use helpme to create a trajectory
# trajectory = splinepy.helpme.create.circle(10)

# insert knots and control points
trajectory.uniform_refine(0, 2)
trajectory.uniform_refine(0, 1)

### CROSS SECTION ###
dict_cross_section = {
Expand Down Expand Up @@ -69,3 +133,50 @@
["Swept Surface", swept_surface],
resolution=101,
)

sys.exit()

### EXPORT A SWEPT SPLINE ###
dict_export_cs = {
"degrees": [1],
"knot_vectors": [[0.0, 0.0, 1.0, 1.0]],
"control_points": np.array(
[
[0.0, 0.0],
[1.0, 0.0],
]
),
}
export_cs = splinepy.BSpline(**dict_export_cs)

dict_export_traj = {
"degrees": [3],
"knot_vectors": [[0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0, 1.0]],
"control_points": np.array(
[
[0.0, 0.0],
[1.0, 2.0],
[2.0, 0.0],
[3.0, -2.0],
[4.0, 0.0],
]
),
}
export_traj = splinepy.BSpline(**dict_export_traj)

swept_surface = splinepy.helpme.create.swept(
trajectory=export_traj,
cross_section=export_cs,
cross_section_normal=[-1, 0, 0],
)

gus.show(
["Trajectory", export_traj],
["Cross Section", export_cs],
["Swept Surface", swept_surface],
resolution=101,
)

projection = swept_surface.create.embedded(2)
gus.show(["Projection", projection], resolution=101)
splinepy.io.mfem.export("testmeshmesh.mesh", projection)
62 changes: 58 additions & 4 deletions splinepy/helpme/create.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -308,14 +308,16 @@ def swept(
receives rotation into the direction of the trajectory tangent
vector and is then placed at the evaluation points of the
trajectory's knots.
Fundamental ideas can be found in the NURBS Book, Piegl & Tiller,
2nd edition, chapter 10.4 Swept Surfaces.
Parameters
----------
cross_section : Spline
Cross-section to be swept
trajectory : Spline
Trajectory along which the cross-section is swept
cross_section_normal : np.ndarray
cross_section_normal : np.array
Normal vector of the cross-section
Default is [0, 0, 1]
Expand Down Expand Up @@ -344,6 +346,9 @@ def swept(
cross_section = cross_section.create.embedded(3)

def transformation_matrices(traj, par_value):

### TRANSFORMATION MATRICES ###

# tangent vector x on trajectory at parameter value 0
x = traj.derivative([par_value[0]], [1])
x = (x / _np.linalg.norm(x)).ravel()
Expand All @@ -358,17 +363,19 @@ def transformation_matrices(traj, par_value):
vec = [x[2], -x[1], x[0]]
B.append(_np.cross(x, vec) / _np.linalg.norm(_np.cross(x, vec)))

# initialize transformation matrices
# initialize transformation matrices and x-collection
T = []
A = []
x_collection = []

# evaluating transformation matrices for each trajectory point
for i in range(len(par_value)):
# tangent vector x on trajectory at parameter value i
x = traj.derivative([par_value[i]], [1])
x = (x / _np.linalg.norm(x)).ravel()
x_collection.append(x)

# projecting B_(i-1) onto the plane normal to x
# projecting B_(i) onto the plane normal to x
B.append(B[i] - _np.dot(B[i], x) * x)
B[i + 1] = B[i + 1] / _np.linalg.norm(B[i + 1])

Expand All @@ -382,7 +389,49 @@ def transformation_matrices(traj, par_value):
# array of transformation matrices from local to global coordinates
A.append(_np.linalg.inv(T[i]))

# rotation matrix around y
# own procedure, if trajectory is closed and B[0] != B[-1]
# according to NURBS Book, Piegl & Tiller, 2nd edition, p. 483
if _np.array_equal(
traj.evaluate([[0]]), traj.evaluate([par_value[-1]])
) and not _np.array_equal(B[0], B[-1]):
# reset transformation matrices
T = []
A = []
B_reverse = [None] * len(B)
B_reverse[0] = B[-1]
for i in range(len(par_value)):
# redo the calculation of B using x_collection from before
B_reverse[i + 1] = (
B_reverse[i]
- _np.dot(B_reverse[i], x_collection[i]) * x_collection[i]
)
B_reverse[i + 1] = B_reverse[i + 1] / _np.linalg.norm(
B_reverse[i + 1]
)
# middle point between B and B_reverse
B_reverse[i + 1] = (B[i + 1] + B_reverse[i + 1]) * 0.5
B_reverse[i + 1] = B_reverse[i + 1] / _np.linalg.norm(
B_reverse[i + 1]
)
# defining y and z axis-vectors
z = B_reverse[i + 1]
y = _np.cross(z, x_collection[i])

# array of transformation matrices from global to local coordinates
T.append(_np.vstack((x_collection[i], y, z)))

# array of transformation matrices from local to global coordinates
A.append(_np.linalg.inv(T[i]))

# check if the beginning and the end of the B-vector are the same
if not _np.allclose(B_reverse[0], B_reverse[-1]):
_log.warning(
"Vector calculation is not exact due to the "
"trajectory being closed in an uncommon way."
)

### ROTATION MATRIX AROUND Y ###

angle_of_cs_normal = _np.arctan2(
cross_section_normal[2], cross_section_normal[0]
)
Expand All @@ -394,8 +443,11 @@ def transformation_matrices(traj, par_value):
]
)
R = _np.where(_np.abs(R) < 1e-10, 0, R)

return A, R

### REFINEMENT OF TRAJECTORY ###

par_value = trajectory.greville_abscissae()
par_value = par_value.reshape(-1, 1)

Expand Down Expand Up @@ -440,6 +492,8 @@ def transformation_matrices(traj, par_value):
par_value = _np.sort(par_value)
par_value = par_value.reshape(-1, 1)

### SWEEPING PROCESS ###

# evaluate trajectory at the parameter values
evals = trajectory.evaluate(par_value)

Expand Down

0 comments on commit 6737fb8

Please sign in to comment.