From 7320d61fefaa2a657c72efb02f0fb0be3b0f9e9a Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Thu, 18 Apr 2024 21:07:07 +0200 Subject: [PATCH 01/53] Add: show_swept-file and swept-function in create-file --- examples/show_swept.py | 69 +++++++++++++++++++++++++++++++++++++++ splinepy/helpme/create.py | 36 ++++++++++++++++++++ 2 files changed, 105 insertions(+) create mode 100644 examples/show_swept.py diff --git a/examples/show_swept.py b/examples/show_swept.py new file mode 100644 index 000000000..5236c45fc --- /dev/null +++ b/examples/show_swept.py @@ -0,0 +1,69 @@ +import gustaf as gus +import numpy as np + +import splinepy + +if __name__ == "__main__": + # trajectory + # define degrees + ds1 = [2] + # define knot vectors + kvs1 = [[0.0, 0.0, 0.0, 1.0, 1.0, 1.0]] + # define control points + cps1 = np.array( + [ + [0.0, 0.0, 0.0], + [3.0, 1.0, 0.0], + [4.0, 0.0, 0.0], + ] + ) + # init trajectory as bspline + trajectory = splinepy.BSpline( + degrees=ds1, + knot_vectors=kvs1, + control_points=cps1, + ) + # show + # trajectory.show() + + # cross section + # define degrees + ds1 = [1] + + # define knot vectors + kvs1 = [[0.0, 0.0, 1.0, 1.0]] + + # define control points + cps1 = np.array( + [ + [0.0, 0.0, -1.0], + [0.0, 0.0, 1.0], + ] + ) + + # init cross section as bspline + cross_section = splinepy.BSpline( + degrees=ds1, + knot_vectors=kvs1, + control_points=cps1, + ) + + def der(data, on): + return data.derivative(on, [1]) + + trajectory.spline_data["der"] = splinepy.SplineDataAdaptor( + trajectory, function=der + ) + trajectory.show_options["arrow_data"] = "der" + trajectory.show_options["arrow_data_on"] = np.linspace(0, 1, 20).reshape( + -1, 1 + ) + + derivative = trajectory.derivative([[0.5], [1]], [1]) + print(derivative) + + gus.show( + ["Trajectory", trajectory], + ["Cross section", cross_section], + resolution=50, + ) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 1e3d4ee3d..bc319ce53 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -298,6 +298,42 @@ def revolved( return type(spline)(**spline_dict) +def swept(cross_section, trajectory, nsections): + """Sweeps a cross-section along a trajectory + + Parameters + ---------- + crossection : Spline + Cross-section to be swept + trajectory : Spline + Trajectory along which the cross-section is swept + nsections : int + Number of sections trajectory is divided into + + Returns + ------- + swept_spline : Spline + Spline resulting from the sweep + """ + + from splinepy.spline import Spline as _Spline + + # Check input type + if not isinstance(cross_section, _Spline): + raise NotImplementedError("Sweeps only works for splines") + if not isinstance(trajectory, _Spline): + raise NotImplementedError("Sweeps only works for splines") + + # Check if trajectory is a curve + if trajectory.dim != 1: + raise ValueError("Trajectory must be a curve") + + # dummy variable for nsections + _ = nsections + + # IMPLEMENTATION OF SWEEPING SURFACE + + def from_bounds(parametric_bounds, physical_bounds): """Creates a minimal spline with given parametric bounds, physical bounds. Physical bounds can have less or equal number of From 8e544a11dee97ae11ee515a12673f7a930980f7e Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Mon, 29 Apr 2024 17:21:59 +0200 Subject: [PATCH 02/53] Add: first kinda swept surface implementation in show_swept-file --- examples/show_swept.py | 60 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 49 insertions(+), 11 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 5236c45fc..c411fad27 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -23,18 +23,27 @@ knot_vectors=kvs1, control_points=cps1, ) - # show - # trajectory.show() + # define sections along trajectory + nsect = len(kvs1[0])-ds1[0]-1 + + # compute parameter values for inserting cross sections + par_value = np.zeros((nsect,1)) + for i in range(nsect): + par_value[i][0] = (kvs1[0][i+1]+kvs1[0][i+2])/2 + + # evaluate trajectory at these parameter values + evals = trajectory.evaluate(par_value) + # cross section # define degrees - ds1 = [1] + ds_cs = [1] # define knot vectors - kvs1 = [[0.0, 0.0, 1.0, 1.0]] + kvs_cs = [[0.0, 0.0, 1.0, 1.0]] # define control points - cps1 = np.array( + cps_cs = np.array( [ [0.0, 0.0, -1.0], [0.0, 0.0, 1.0], @@ -43,11 +52,38 @@ # init cross section as bspline cross_section = splinepy.BSpline( - degrees=ds1, - knot_vectors=kvs1, - control_points=cps1, + degrees=ds_cs, + knot_vectors=kvs_cs, + control_points=cps_cs, + ) + + # set cross section CPs along trajectory + cps_eval = [] + for eval_point in evals: + # for-loop representing the nr. of cross section CPs + # note: only working for this special case + for i in range(len(cps_cs)): + current_cps = [] + current_cps.append(eval_point[0]) + current_cps.append(eval_point[1]) + current_cps.append(cps_cs[i][2]) + cps_eval.append(current_cps) + + fitting_points = np.array(cps_eval) + + # fit surface + interpolated_surface, _ = splinepy.helpme.fit.surface( + fitting_points=fitting_points, + size=[2, 3], + n_control_points=[2, 3], + degrees=[1, 2], ) + + # testing derivative and evaluation + # derivative = trajectory.derivative([[0.5], [1]], [1]) + # evals = trajectory.evaluate([[0.5]]) + # show def der(data, on): return data.derivative(on, [1]) @@ -59,11 +95,13 @@ def der(data, on): -1, 1 ) - derivative = trajectory.derivative([[0.5], [1]], [1]) - print(derivative) - gus.show( ["Trajectory", trajectory], ["Cross section", cross_section], resolution=50, ) + + gus.show( + ["Surface", interpolated_surface], + resolution=50, + ) From f63b2d5c4ab445a7770af62a1b15ab55347750d0 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Thu, 2 May 2024 14:12:41 +0200 Subject: [PATCH 03/53] Fix: Moved computation from show_swept-file to create-file --- examples/show_swept.py | 32 ++++-------------------------- splinepy/helpme/create.py | 41 ++++++++++++++++++++++++++++++--------- 2 files changed, 36 insertions(+), 37 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index c411fad27..d1246eafe 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -27,14 +27,6 @@ # define sections along trajectory nsect = len(kvs1[0])-ds1[0]-1 - # compute parameter values for inserting cross sections - par_value = np.zeros((nsect,1)) - for i in range(nsect): - par_value[i][0] = (kvs1[0][i+1]+kvs1[0][i+2])/2 - - # evaluate trajectory at these parameter values - evals = trajectory.evaluate(par_value) - # cross section # define degrees ds_cs = [1] @@ -57,26 +49,10 @@ control_points=cps_cs, ) - # set cross section CPs along trajectory - cps_eval = [] - for eval_point in evals: - # for-loop representing the nr. of cross section CPs - # note: only working for this special case - for i in range(len(cps_cs)): - current_cps = [] - current_cps.append(eval_point[0]) - current_cps.append(eval_point[1]) - current_cps.append(cps_cs[i][2]) - cps_eval.append(current_cps) - - fitting_points = np.array(cps_eval) - - # fit surface - interpolated_surface, _ = splinepy.helpme.fit.surface( - fitting_points=fitting_points, - size=[2, 3], - n_control_points=[2, 3], - degrees=[1, 2], + interpolated_surface = splinepy.helpme.create.swept( + trajectory=trajectory, + cross_section=cross_section, + nsections=nsect, ) # testing derivative and evaluation diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index bc319ce53..cec1be6bb 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -6,7 +6,7 @@ from splinepy import settings as _settings from splinepy.utils import log as _log from splinepy.utils.data import make_matrix as _make_matrix - +from splinepy.helpme.fit import surface def embedded(spline, new_dimension): """Embeds Spline in given dimension. @@ -324,15 +324,38 @@ def swept(cross_section, trajectory, nsections): if not isinstance(trajectory, _Spline): raise NotImplementedError("Sweeps only works for splines") - # Check if trajectory is a curve - if trajectory.dim != 1: - raise ValueError("Trajectory must be a curve") - - # dummy variable for nsections - _ = nsections - - # IMPLEMENTATION OF SWEEPING SURFACE + # compute parameter values for inserting cross sections + par_value = _np.zeros((nsections,1)) + for i in range(nsections): + par_value[i][0] = (trajectory.knot_vectors[0][i+1]+trajectory.knot_vectors[0][i+2])/2 + + # evaluate trajectory at these parameter values + evals = trajectory.evaluate(par_value) + + # set cross section CPs along trajectory + cps_eval = [] + for eval_point in evals: + # for-loop representing the nr. of cross section CPs + # note: only working for this special case + for i in range(len(cross_section.control_points)): + current_cps = [] + current_cps.append(eval_point[0]) + current_cps.append(eval_point[1]) + current_cps.append(cross_section.control_points[i][2]) + cps_eval.append(current_cps) + + fitting_points = _np.array(cps_eval) + + # fit surface + interpolated_surface, _ = surface( + fitting_points=fitting_points, + size=[2, 3], + n_control_points=[2, 3], + degrees=[1, 2], + ) + return interpolated_surface + def from_bounds(parametric_bounds, physical_bounds): """Creates a minimal spline with given parametric bounds, physical bounds. From d871575f03869909bcab1d119b071650cf477578 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Thu, 2 May 2024 14:51:15 +0200 Subject: [PATCH 04/53] Add: function for transformation matrix --- splinepy/helpme/create.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index cec1be6bb..7272e4c9d 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -324,6 +324,23 @@ def swept(cross_section, trajectory, nsections): if not isinstance(trajectory, _Spline): raise NotImplementedError("Sweeps only works for splines") + # compute transformation matrix + # parameters: trajectory traj and parametric value v + def transformation_matrix(traj, v): + # origin of local coordinate system + o = traj.evaluate([[v]]) + # local directions in global coordinates + x = traj.derivative([[v]], [1]) / _np.linalg.norm(traj.derivative([[v]], [1]) ) + z = _np.cross(traj.derivative([[v]], [1]), traj.derivative([[v]], [2])) / _np.linalg.norm(_np.cross(traj.derivative([[v]], [1]), traj.derivative([[v]], [2]))) + y = _np.cross(x,z) + # transformation matrix + A = [] + A = _np.vstack((x, y, z)) + return A + + A = transformation_matrix(trajectory, 0) + print(A) + # compute parameter values for inserting cross sections par_value = _np.zeros((nsections,1)) for i in range(nsections): From f3bcd76cb28a6950505aaee7be40c4697f27db1f Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Thu, 2 May 2024 17:44:35 +0200 Subject: [PATCH 05/53] Add: show_options changed in show_swept-file --- examples/show_swept.py | 5 ++++- splinepy/helpme/create.py | 14 ++++++++------ 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index d1246eafe..2f1cbfa07 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -54,7 +54,7 @@ cross_section=cross_section, nsections=nsect, ) - + # testing derivative and evaluation # derivative = trajectory.derivative([[0.5], [1]], [1]) # evals = trajectory.evaluate([[0.5]]) @@ -70,6 +70,9 @@ def der(data, on): trajectory.show_options["arrow_data_on"] = np.linspace(0, 1, 20).reshape( -1, 1 ) + trajectory.show_options["control_mesh"] = False + cross_section.show_options["control_mesh"] = False + interpolated_surface.show_options["control_mesh"] = False gus.show( ["Trajectory", trajectory], diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 7272e4c9d..89b34a97c 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -332,12 +332,14 @@ def transformation_matrix(traj, v): # local directions in global coordinates x = traj.derivative([[v]], [1]) / _np.linalg.norm(traj.derivative([[v]], [1]) ) z = _np.cross(traj.derivative([[v]], [1]), traj.derivative([[v]], [2])) / _np.linalg.norm(_np.cross(traj.derivative([[v]], [1]), traj.derivative([[v]], [2]))) - y = _np.cross(x,z) - # transformation matrix - A = [] - A = _np.vstack((x, y, z)) + y = _np.cross(z,x) + # transformation matrix from global to local coordinates + T = [] + T = _np.vstack((x, y, z)) + # transformation matrix from local to global coordinates + A = _np.linalg.inv(T) return A - + A = transformation_matrix(trajectory, 0) print(A) @@ -349,7 +351,7 @@ def transformation_matrix(traj, v): # evaluate trajectory at these parameter values evals = trajectory.evaluate(par_value) - # set cross section CPs along trajectory + # set cross section evaluation points along trajectory cps_eval = [] for eval_point in evals: # for-loop representing the nr. of cross section CPs From a6c3eba9324c809a824ae8eb04764944b868c3cc Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Mon, 6 May 2024 12:45:54 +0200 Subject: [PATCH 06/53] Add: positioning of cross-section control points --- examples/show_swept.py | 12 ++++--- splinepy/helpme/create.py | 69 ++++++++++++++++++++++++--------------- 2 files changed, 50 insertions(+), 31 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 2f1cbfa07..d1f675a7c 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -25,20 +25,22 @@ ) # define sections along trajectory - nsect = len(kvs1[0])-ds1[0]-1 + nsect = len(kvs1[0]) - ds1[0] - 1 # cross section # define degrees - ds_cs = [1] + ds_cs = [2] # define knot vectors - kvs_cs = [[0.0, 0.0, 1.0, 1.0]] + kvs_cs = [[0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0]] # define control points cps_cs = np.array( [ - [0.0, 0.0, -1.0], - [0.0, 0.0, 1.0], + [0.0, 0.0, 0.0], + [1.0, 1.0, 0.0], + [2.0, 0.0, 0.0], + [3.0, 1.0, 0.0], ] ) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 89b34a97c..202df67e0 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -4,9 +4,10 @@ from gustaf.utils import arr as _arr from splinepy import settings as _settings +from splinepy.helpme.fit import surface from splinepy.utils import log as _log from splinepy.utils.data import make_matrix as _make_matrix -from splinepy.helpme.fit import surface + def embedded(spline, new_dimension): """Embeds Spline in given dimension. @@ -324,57 +325,73 @@ def swept(cross_section, trajectory, nsections): if not isinstance(trajectory, _Spline): raise NotImplementedError("Sweeps only works for splines") - # compute transformation matrix + # compute transformation matrix # parameters: trajectory traj and parametric value v def transformation_matrix(traj, v): - # origin of local coordinate system - o = traj.evaluate([[v]]) # local directions in global coordinates - x = traj.derivative([[v]], [1]) / _np.linalg.norm(traj.derivative([[v]], [1]) ) - z = _np.cross(traj.derivative([[v]], [1]), traj.derivative([[v]], [2])) / _np.linalg.norm(_np.cross(traj.derivative([[v]], [1]), traj.derivative([[v]], [2]))) - y = _np.cross(z,x) + x = traj.derivative([[v]], [1]) / _np.linalg.norm( + traj.derivative([[v]], [1]) + ) + z = _np.cross( + traj.derivative([[v]], [1]), traj.derivative([[v]], [2]) + ) / _np.linalg.norm( + _np.cross(traj.derivative([[v]], [1]), traj.derivative([[v]], [2])) + ) + y = _np.cross(z, x) # transformation matrix from global to local coordinates T = [] T = _np.vstack((x, y, z)) # transformation matrix from local to global coordinates A = _np.linalg.inv(T) + A = _np.negative(A) return A - - A = transformation_matrix(trajectory, 0) - print(A) # compute parameter values for inserting cross sections - par_value = _np.zeros((nsections,1)) + par_value = _np.zeros((nsections, 1)) for i in range(nsections): - par_value[i][0] = (trajectory.knot_vectors[0][i+1]+trajectory.knot_vectors[0][i+2])/2 + par_value[i][0] = ( + trajectory.knot_vectors[0][i + 1] + + trajectory.knot_vectors[0][i + 2] + ) / 2 # evaluate trajectory at these parameter values evals = trajectory.evaluate(par_value) # set cross section evaluation points along trajectory cps_eval = [] - for eval_point in evals: - # for-loop representing the nr. of cross section CPs - # note: only working for this special case - for i in range(len(cross_section.control_points)): - current_cps = [] - current_cps.append(eval_point[0]) - current_cps.append(eval_point[1]) - current_cps.append(cross_section.control_points[i][2]) - cps_eval.append(current_cps) + for index, eval_point in enumerate(evals): + # evaluate transformation matrix for each trajectory point + A = transformation_matrix(trajectory, par_value[index][0]) + for cscp in cross_section.control_points: + # rotation matrix for 90 degrees around y-axis + R = _np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) + # rotate cross section in trajectory direction + normal_cscp = _np.dot(R, cscp) + # transform cross section to global coordinates + normal_cscp = _np.dot(A, normal_cscp) + # translate cross section to trajectory point + normal_cscp = normal_cscp + eval_point + # append control point to list + cps_eval.append(normal_cscp) fitting_points = _np.array(cps_eval) - + # fit surface interpolated_surface, _ = surface( fitting_points=fitting_points, - size=[2, 3], - n_control_points=[2, 3], - degrees=[1, 2], + size=[ + len(cross_section.control_points), + len(trajectory.control_points), + ], + n_control_points=[ + len(cross_section.control_points), + len(trajectory.control_points), + ], + degrees=[int(cross_section.degrees), int(trajectory.degrees)], ) return interpolated_surface - + def from_bounds(parametric_bounds, physical_bounds): """Creates a minimal spline with given parametric bounds, physical bounds. From 78ea5a2ebb676dd19651fd18aea216f253f16f5b Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Mon, 6 May 2024 14:05:42 +0200 Subject: [PATCH 07/53] Fix: now, evaluation-points are fitted instead of CPs --- splinepy/helpme/create.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 202df67e0..341db4f81 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -3,6 +3,7 @@ import numpy as _np from gustaf.utils import arr as _arr +import splinepy from splinepy import settings as _settings from splinepy.helpme.fit import surface from splinepy.utils import log as _log @@ -359,6 +360,7 @@ def transformation_matrix(traj, v): # set cross section evaluation points along trajectory cps_eval = [] + cs_evals = [] for index, eval_point in enumerate(evals): # evaluate transformation matrix for each trajectory point A = transformation_matrix(trajectory, par_value[index][0]) @@ -374,7 +376,23 @@ def transformation_matrix(traj, v): # append control point to list cps_eval.append(normal_cscp) - fitting_points = _np.array(cps_eval) + cross_sec_placed_cps = _np.array(cps_eval) + number_of_cps_we_take = len(cross_section.control_points) + new_cross_section = splinepy.BSpline( + degrees=cross_section.degrees, + knot_vectors=cross_section.knot_vectors, + control_points=cross_sec_placed_cps[ + index + * number_of_cps_we_take : (index + 1) + * number_of_cps_we_take + ], + ) + + # take care, this has to be the same as the nr of CPs each CS for now + evaluations = new_cross_section.evaluate([[0], [0.33], [0.66], [1]]) + cs_evals.append(evaluations) + + fitting_points = _np.array(cs_evals).reshape(-1, 3) # fit surface interpolated_surface, _ = surface( From 5f1666900fe2f40ce42972b18c1ae2c297cf0430 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Tue, 7 May 2024 14:30:55 +0200 Subject: [PATCH 08/53] Add: rotation matrix and centered cross section --- examples/show_swept.py | 5 +++++ splinepy/helpme/create.py | 42 +++++++++++++++++++++++++++++++++------ 2 files changed, 41 insertions(+), 6 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index d1f675a7c..d238b4385 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -51,10 +51,15 @@ control_points=cps_cs, ) + # user should define the normal vector of the cross section + cs_nv = np.array([0, 0, 1]) + + # create interpolated surface interpolated_surface = splinepy.helpme.create.swept( trajectory=trajectory, cross_section=cross_section, nsections=nsect, + cross_section_normal=cs_nv, ) # testing derivative and evaluation diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 341db4f81..ad91b1857 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -300,7 +300,12 @@ def revolved( return type(spline)(**spline_dict) -def swept(cross_section, trajectory, nsections): +def swept( + cross_section, + trajectory, + nsections, + cross_section_normal=None, +): """Sweeps a cross-section along a trajectory Parameters @@ -311,6 +316,9 @@ def swept(cross_section, trajectory, nsections): Trajectory along which the cross-section is swept nsections : int Number of sections trajectory is divided into + cross_section_normal : np.ndarray + Normal vector of the cross-section + Default is [0, 0, 1] Returns ------- @@ -326,6 +334,9 @@ def swept(cross_section, trajectory, nsections): if not isinstance(trajectory, _Spline): raise NotImplementedError("Sweeps only works for splines") + if cross_section_normal is None: + cross_section_normal = _np.array([0, 0, 1]) + # compute transformation matrix # parameters: trajectory traj and parametric value v def transformation_matrix(traj, v): @@ -338,14 +349,31 @@ def transformation_matrix(traj, v): ) / _np.linalg.norm( _np.cross(traj.derivative([[v]], [1]), traj.derivative([[v]], [2])) ) + if z[0][2] < 0: + z = -z y = _np.cross(z, x) # transformation matrix from global to local coordinates T = [] T = _np.vstack((x, y, z)) + # transformation matrix from local to global coordinates A = _np.linalg.inv(T) - A = _np.negative(A) - return A + + # rotation matrix around y + angle_of_x = _np.arctan2(x[0][2], x[0][0]) + angle_of_cs_normal = _np.arctan2( + cross_section_normal[2], cross_section_normal[0] + ) + rotation_angle = angle_of_cs_normal - angle_of_x + R = _np.array( + [ + [_np.cos(rotation_angle), 0, _np.sin(rotation_angle)], + [0, 1, 0], + [-_np.sin(rotation_angle), 0, _np.cos(rotation_angle)], + ] + ) + + return A, R # compute parameter values for inserting cross sections par_value = _np.zeros((nsections, 1)) @@ -358,15 +386,17 @@ def transformation_matrix(traj, v): # evaluate trajectory at these parameter values evals = trajectory.evaluate(par_value) + # find a center of the cross section + cs_center = cross_section.evaluate([[0.5]]).flatten() + cross_section.control_points = cross_section.control_points - cs_center + # set cross section evaluation points along trajectory cps_eval = [] cs_evals = [] for index, eval_point in enumerate(evals): # evaluate transformation matrix for each trajectory point - A = transformation_matrix(trajectory, par_value[index][0]) + A, R = transformation_matrix(trajectory, par_value[index][0]) for cscp in cross_section.control_points: - # rotation matrix for 90 degrees around y-axis - R = _np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) # rotate cross section in trajectory direction normal_cscp = _np.dot(R, cscp) # transform cross section to global coordinates From afa5569acc30bab0e0346d28ca60fc34cbcab86e Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Tue, 7 May 2024 16:36:47 +0200 Subject: [PATCH 09/53] Fix: getting fitting flexible for CS eval points; Note: WIP --- splinepy/helpme/create.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index ad91b1857..73ffa93a7 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -334,6 +334,7 @@ def swept( if not isinstance(trajectory, _Spline): raise NotImplementedError("Sweeps only works for splines") + # setting default value for cross_section_normal if cross_section_normal is None: cross_section_normal = _np.array([0, 0, 1]) @@ -349,6 +350,7 @@ def transformation_matrix(traj, v): ) / _np.linalg.norm( _np.cross(traj.derivative([[v]], [1]), traj.derivative([[v]], [2])) ) + # check if z is pointing in the right direction if z[0][2] < 0: z = -z y = _np.cross(z, x) @@ -372,7 +374,6 @@ def transformation_matrix(traj, v): [-_np.sin(rotation_angle), 0, _np.cos(rotation_angle)], ] ) - return A, R # compute parameter values for inserting cross sections @@ -386,7 +387,7 @@ def transformation_matrix(traj, v): # evaluate trajectory at these parameter values evals = trajectory.evaluate(par_value) - # find a center of the cross section + # translate cross section to origin cs_center = cross_section.evaluate([[0.5]]).flatten() cross_section.control_points = cross_section.control_points - cs_center @@ -418,22 +419,24 @@ def transformation_matrix(traj, v): ], ) - # take care, this has to be the same as the nr of CPs each CS for now - evaluations = new_cross_section.evaluate([[0], [0.33], [0.66], [1]]) + # evaluate cross section at parameter values + evaluations = new_cross_section.evaluate( + [[0], [0.25], [0.5], [0.75], [1]] + ) cs_evals.append(evaluations) fitting_points = _np.array(cs_evals).reshape(-1, 3) - # fit surface + # fit surface - take care of size and n_control_points --> not sure yet interpolated_surface, _ = surface( fitting_points=fitting_points, size=[ - len(cross_section.control_points), - len(trajectory.control_points), + len(evaluations), + nsections, ], n_control_points=[ - len(cross_section.control_points), - len(trajectory.control_points), + len(evaluations), + nsections, ], degrees=[int(cross_section.degrees), int(trajectory.degrees)], ) From 3a4cb530b3af9874e301a7f2747d5fe57adb7902 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Tue, 14 May 2024 17:28:20 +0200 Subject: [PATCH 10/53] Add: Transformation matrix now with projection; Note: WIP --- examples/show_swept.py | 8 ++--- splinepy/helpme/create.py | 74 ++++++++++++++++++++++++++++++--------- 2 files changed, 61 insertions(+), 21 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index d238b4385..339b96b4b 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -13,8 +13,8 @@ cps1 = np.array( [ [0.0, 0.0, 0.0], - [3.0, 1.0, 0.0], - [4.0, 0.0, 0.0], + [3.0, 1, 0.0], + [4.0, 3, 0.0], ] ) # init trajectory as bspline @@ -29,10 +29,10 @@ # cross section # define degrees - ds_cs = [2] + ds_cs = [3] # define knot vectors - kvs_cs = [[0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0]] + kvs_cs = [[0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]] # define control points cps_cs = np.array( diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 73ffa93a7..7f10abfd4 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -330,9 +330,9 @@ def swept( # Check input type if not isinstance(cross_section, _Spline): - raise NotImplementedError("Sweeps only works for splines") + raise NotImplementedError("Sweep only works for splines") if not isinstance(trajectory, _Spline): - raise NotImplementedError("Sweeps only works for splines") + raise NotImplementedError("Sweep only works for splines") # setting default value for cross_section_normal if cross_section_normal is None: @@ -340,29 +340,67 @@ def swept( # compute transformation matrix # parameters: trajectory traj and parametric value v - def transformation_matrix(traj, v): + # def transformation_matrix(traj, v): + # # local directions in global coordinates + # x = traj.derivative([[v]], [1]) / _np.linalg.norm( + # traj.derivative([[v]], [1]) + # ) + # z = _np.cross( + # traj.derivative([[v]], [1]), traj.derivative([[v]], [2]) + # ) / _np.linalg.norm( + # _np.cross(traj.derivative([[v]], [1]), traj.derivative([[v]], [2]))) + + # # check if z is pointing in the right direction + # if z[0][2] < 0: + # z = -z + # y = _np.cross(z, x) + # # transformation matrix from global to local coordinates + # T = [] + # 3 + # T = _np.vstack((x, y, z)) + + # # transformation matrix from local to global coordinates + # A = _np.linalg.inv(T) + + # # rotation matrix around y + # angle_of_x = _np.arctan2(x[0][2], x[0][0]) + # angle_of_cs_normal = _np.arctan2( + # cross_section_normal[2], cross_section_normal[0] + # ) + # rotation_angle = angle_of_cs_normal - angle_of_x + # R = _np.array( + # [ + # [_np.cos(rotation_angle), 0, _np.sin(rotation_angle)], + # [0, 1, 0], + # [-_np.sin(rotation_angle), 0, _np.cos(rotation_angle)], + # ] + # ) + # return A, R + + def transformation_matrix_with_projection(traj, v): # local directions in global coordinates - x = traj.derivative([[v]], [1]) / _np.linalg.norm( + x = ( traj.derivative([[v]], [1]) - ) - z = _np.cross( - traj.derivative([[v]], [1]), traj.derivative([[v]], [2]) - ) / _np.linalg.norm( - _np.cross(traj.derivative([[v]], [1]), traj.derivative([[v]], [2])) - ) - # check if z is pointing in the right direction - if z[0][2] < 0: - z = -z + / _np.linalg.norm(traj.derivative([[v]], [1])) + ).flatten() + # arbitrary vector B_0 normal to x + vec = [1, 1, 1] + B_0 = _np.cross(x, vec) + # projecting B_0 onto the plane normal to x + b_i = B_0 - _np.dot(B_0, x) * x + z = b_i / _np.linalg.norm(b_i) + y = _np.cross(z, x) # transformation matrix from global to local coordinates T = [] + T = _np.vstack((x, y, z)) # transformation matrix from local to global coordinates A = _np.linalg.inv(T) # rotation matrix around y - angle_of_x = _np.arctan2(x[0][2], x[0][0]) + angle_of_x = _np.arctan2(x[2], x[0]) angle_of_cs_normal = _np.arctan2( cross_section_normal[2], cross_section_normal[0] ) @@ -396,7 +434,9 @@ def transformation_matrix(traj, v): cs_evals = [] for index, eval_point in enumerate(evals): # evaluate transformation matrix for each trajectory point - A, R = transformation_matrix(trajectory, par_value[index][0]) + A, R = transformation_matrix_with_projection( + trajectory, par_value[index][0] + ) for cscp in cross_section.control_points: # rotate cross section in trajectory direction normal_cscp = _np.dot(R, cscp) @@ -435,8 +475,8 @@ def transformation_matrix(traj, v): nsections, ], n_control_points=[ - len(evaluations), - nsections, + len(cross_section.control_points), + len(trajectory.control_points), ], degrees=[int(cross_section.degrees), int(trajectory.degrees)], ) From 9f07a672216b01da80dc3b43b7dd5da85d918bcf Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 15 May 2024 16:17:49 +0200 Subject: [PATCH 11/53] Fix: calculation of transformation matrix - still working on projection; Note: WIP --- examples/show_swept.py | 8 ++--- splinepy/helpme/create.py | 72 ++++++++++++++++++++++++--------------- 2 files changed, 49 insertions(+), 31 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 339b96b4b..d238b4385 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -13,8 +13,8 @@ cps1 = np.array( [ [0.0, 0.0, 0.0], - [3.0, 1, 0.0], - [4.0, 3, 0.0], + [3.0, 1.0, 0.0], + [4.0, 0.0, 0.0], ] ) # init trajectory as bspline @@ -29,10 +29,10 @@ # cross section # define degrees - ds_cs = [3] + ds_cs = [2] # define knot vectors - kvs_cs = [[0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]] + kvs_cs = [[0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0]] # define control points cps_cs = np.array( diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 7f10abfd4..a756b3084 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -356,7 +356,6 @@ def swept( # y = _np.cross(z, x) # # transformation matrix from global to local coordinates # T = [] - # 3 # T = _np.vstack((x, y, z)) # # transformation matrix from local to global coordinates @@ -377,41 +376,54 @@ def swept( # ) # return A, R - def transformation_matrix_with_projection(traj, v): - # local directions in global coordinates - x = ( - traj.derivative([[v]], [1]) - / _np.linalg.norm(traj.derivative([[v]], [1])) - ).flatten() + def transformation_matrix_with_projection(traj, index, par_value): + # for i in range(len(traj.knot_vectors[0])): + x = traj.derivative([[0]], [1]) / _np.linalg.norm( + traj.derivative([[0]], [1]) + ) # arbitrary vector B_0 normal to x - vec = [1, 1, 1] - B_0 = _np.cross(x, vec) - # projecting B_0 onto the plane normal to x - b_i = B_0 - _np.dot(B_0, x) * x - z = b_i / _np.linalg.norm(b_i) + vec = [1.5, 0.3, 1] + B = [] + B.append(_np.cross(x, vec) / _np.linalg.norm(_np.cross(x, vec))) + z = B[0] + + for i in range(1, index + 1): + # local directions in global coordinates + x = ( + traj.derivative([[par_value[i][0]]], [1]) + / _np.linalg.norm(traj.derivative([[par_value[i][0]]], [1])) + ).flatten() + + # projecting B_i onto the plane normal to x + B.append(B[i - 1] - _np.dot(B[i - 1], x) * x) + B[i] = B[i] / _np.linalg.norm(B[i]) + z = B[i] y = _np.cross(z, x) - # transformation matrix from global to local coordinates - T = [] + # transformation matrix from global to local coordinates T = _np.vstack((x, y, z)) # transformation matrix from local to global coordinates A = _np.linalg.inv(T) # rotation matrix around y - angle_of_x = _np.arctan2(x[2], x[0]) - angle_of_cs_normal = _np.arctan2( - cross_section_normal[2], cross_section_normal[0] - ) - rotation_angle = angle_of_cs_normal - angle_of_x - R = _np.array( - [ - [_np.cos(rotation_angle), 0, _np.sin(rotation_angle)], - [0, 1, 0], - [-_np.sin(rotation_angle), 0, _np.cos(rotation_angle)], - ] - ) + # angle_of_x = _np.arctan2(x[2], x[0]) + # angle_of_cs_normal = _np.arctan2( + # cross_section_normal[2], cross_section_normal[0] + # ) + # rotation_angle = angle_of_cs_normal - angle_of_x + # R = [] + # R.append(_np.array( + # [ + # [_np.cos(rotation_angle), 0, _np.sin(rotation_angle)], + # [0, 1, 0], + # [-_np.sin(rotation_angle), 0, _np.cos(rotation_angle)], + # ] + # )) + + # rotation matrix for rotation 90° around z-axis + R = _np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]]) return A, R # compute parameter values for inserting cross sections @@ -429,16 +441,22 @@ def transformation_matrix_with_projection(traj, v): cs_center = cross_section.evaluate([[0.5]]).flatten() cross_section.control_points = cross_section.control_points - cs_center + # # evaluate transformation matrix for each trajectory point + # A, R, B = transformation_matrix_with_projection( + # trajectory + # ) + # set cross section evaluation points along trajectory cps_eval = [] cs_evals = [] for index, eval_point in enumerate(evals): # evaluate transformation matrix for each trajectory point A, R = transformation_matrix_with_projection( - trajectory, par_value[index][0] + trajectory, index, par_value ) for cscp in cross_section.control_points: # rotate cross section in trajectory direction + # cscp = cscp.reshape(-1, 1) normal_cscp = _np.dot(R, cscp) # transform cross section to global coordinates normal_cscp = _np.dot(A, normal_cscp) From d4f1b8f34da5b40497590d3c26c0dfa17c9221d3 Mon Sep 17 00:00:00 2001 From: Jaewook Lee Date: Tue, 21 May 2024 15:37:30 +0200 Subject: [PATCH 12/53] minor clean up during discussion --- splinepy/helpme/create.py | 41 ++++++++++++++++++++++----------------- 1 file changed, 23 insertions(+), 18 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index a756b3084..8f4930967 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -303,7 +303,7 @@ def revolved( def swept( cross_section, trajectory, - nsections, + nsections=None, cross_section_normal=None, ): """Sweeps a cross-section along a trajectory @@ -338,6 +338,10 @@ def swept( if cross_section_normal is None: cross_section_normal = _np.array([0, 0, 1]) + # make copies so we can work on it inplace + trajectory = trajectory.copy() + cross_section = cross_section.copy() + # compute transformation matrix # parameters: trajectory traj and parametric value v # def transformation_matrix(traj, v): @@ -378,8 +382,8 @@ def swept( def transformation_matrix_with_projection(traj, index, par_value): # for i in range(len(traj.knot_vectors[0])): - x = traj.derivative([[0]], [1]) / _np.linalg.norm( - traj.derivative([[0]], [1]) + x = traj.derivative([par_value[0]], [1]) / _np.linalg.norm( + traj.derivative([par_value[0]], [1]) ) # arbitrary vector B_0 normal to x vec = [1.5, 0.3, 1] @@ -390,8 +394,8 @@ def transformation_matrix_with_projection(traj, index, par_value): for i in range(1, index + 1): # local directions in global coordinates x = ( - traj.derivative([[par_value[i][0]]], [1]) - / _np.linalg.norm(traj.derivative([[par_value[i][0]]], [1])) + traj.derivative([par_value[i]], [1]) + / _np.linalg.norm(traj.derivative([par_value[i]], [1])) ).flatten() # projecting B_i onto the plane normal to x @@ -427,19 +431,22 @@ def transformation_matrix_with_projection(traj, index, par_value): return A, R # compute parameter values for inserting cross sections - par_value = _np.zeros((nsections, 1)) - for i in range(nsections): - par_value[i][0] = ( - trajectory.knot_vectors[0][i + 1] - + trajectory.knot_vectors[0][i + 2] - ) / 2 + if nsections is None: + par_value = trajectory.greville_abscssae() + else: + bounds = trajectory.parametric_bounds.ravel() + par_value = _np.linspace(*bounds, nsections) + par_value = par_value.reshape(-1, 1) # evaluate trajectory at these parameter values evals = trajectory.evaluate(par_value) # translate cross section to origin - cs_center = cross_section.evaluate([[0.5]]).flatten() - cross_section.control_points = cross_section.control_points - cs_center + cross_para_center = _np.mean(cross_section.parametric_bounds, axis=0) + cs_center = cross_section.evaluate( + cross_para_center.reshape(-1, 1) + ).ravel() + cross_section.control_points -= cs_center # # evaluate transformation matrix for each trajectory point # A, R, B = transformation_matrix_with_projection( @@ -457,9 +464,9 @@ def transformation_matrix_with_projection(traj, index, par_value): for cscp in cross_section.control_points: # rotate cross section in trajectory direction # cscp = cscp.reshape(-1, 1) - normal_cscp = _np.dot(R, cscp) + normal_cscp = _np.matmul(R, cscp) # transform cross section to global coordinates - normal_cscp = _np.dot(A, normal_cscp) + normal_cscp = _np.matmul(A, normal_cscp) # translate cross section to trajectory point normal_cscp = normal_cscp + eval_point # append control point to list @@ -478,9 +485,7 @@ def transformation_matrix_with_projection(traj, index, par_value): ) # evaluate cross section at parameter values - evaluations = new_cross_section.evaluate( - [[0], [0.25], [0.5], [0.75], [1]] - ) + evaluations = new_cross_section.sample(7) cs_evals.append(evaluations) fitting_points = _np.array(cs_evals).reshape(-1, 3) From 7d07e79780fc50997f54fdd9b873ed62fa165f68 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 22 May 2024 20:28:49 +0200 Subject: [PATCH 13/53] Fix: create-file transformation matrices now calculated once; rotation works in this case --- splinepy/helpme/create.py | 68 ++++++++++++++++++++------------------- 1 file changed, 35 insertions(+), 33 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 8f4930967..26c25d3f6 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -344,7 +344,7 @@ def swept( # compute transformation matrix # parameters: trajectory traj and parametric value v - # def transformation_matrix(traj, v): + # def transformation_matrix_without_projection(traj, v): # # local directions in global coordinates # x = traj.derivative([[v]], [1]) / _np.linalg.norm( # traj.derivative([[v]], [1]) @@ -380,36 +380,42 @@ def swept( # ) # return A, R - def transformation_matrix_with_projection(traj, index, par_value): - # for i in range(len(traj.knot_vectors[0])): + def transformation_matrices(traj, par_value): + # tangent vector x on trajectory at parameter value 0 x = traj.derivative([par_value[0]], [1]) / _np.linalg.norm( traj.derivative([par_value[0]], [1]) ) - # arbitrary vector B_0 normal to x + + # evaluating arbitrary vector B_0 normal to x vec = [1.5, 0.3, 1] B = [] B.append(_np.cross(x, vec) / _np.linalg.norm(_np.cross(x, vec))) - z = B[0] - for i in range(1, index + 1): - # local directions in global coordinates + # initializing transformation matrices + T = [] + A = [] + + # 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]) / _np.linalg.norm(traj.derivative([par_value[i]], [1])) ).flatten() - # projecting B_i onto the plane normal to x - B.append(B[i - 1] - _np.dot(B[i - 1], x) * x) - B[i] = B[i] / _np.linalg.norm(B[i]) - z = B[i] + # projecting B_(i-1) 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]) - y = _np.cross(z, x) + # defining y and z axis-vectors + z = B[i + 1] + y = _np.cross(z, x) - # transformation matrix from global to local coordinates - T = _np.vstack((x, y, z)) + # array of transformation matrices from global to local coordinates + T.append(_np.vstack((x, y, z))) - # transformation matrix from local to global coordinates - A = _np.linalg.inv(T) + # array of transformation matrices from local to global coordinates + A.append(_np.linalg.inv(T[i])) # rotation matrix around y # angle_of_x = _np.arctan2(x[2], x[0]) @@ -426,8 +432,8 @@ def transformation_matrix_with_projection(traj, index, par_value): # ] # )) - # rotation matrix for rotation 90° around z-axis - R = _np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]]) + # rotation matrix for rotation 90° around z-axis and 90° around y-axis + R = _np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]]) return A, R # compute parameter values for inserting cross sections @@ -441,40 +447,35 @@ def transformation_matrix_with_projection(traj, index, par_value): # evaluate trajectory at these parameter values evals = trajectory.evaluate(par_value) - # translate cross section to origin + # evaluate center of cross section and translate to origin cross_para_center = _np.mean(cross_section.parametric_bounds, axis=0) cs_center = cross_section.evaluate( cross_para_center.reshape(-1, 1) ).ravel() cross_section.control_points -= cs_center - # # evaluate transformation matrix for each trajectory point - # A, R, B = transformation_matrix_with_projection( - # trajectory - # ) + # evaluate transformation matrices for every trajectory point + A, R = transformation_matrices(trajectory, par_value) # set cross section evaluation points along trajectory cps_eval = [] cs_evals = [] for index, eval_point in enumerate(evals): - # evaluate transformation matrix for each trajectory point - A, R = transformation_matrix_with_projection( - trajectory, index, par_value - ) + # place every control point of cross section separately for cscp in cross_section.control_points: # rotate cross section in trajectory direction - # cscp = cscp.reshape(-1, 1) normal_cscp = _np.matmul(R, cscp) # transform cross section to global coordinates - normal_cscp = _np.matmul(A, normal_cscp) + normal_cscp = _np.matmul(A[index], normal_cscp) # translate cross section to trajectory point normal_cscp = normal_cscp + eval_point # append control point to list cps_eval.append(normal_cscp) + # create spline of new cross section cross_sec_placed_cps = _np.array(cps_eval) number_of_cps_we_take = len(cross_section.control_points) - new_cross_section = splinepy.BSpline( + temp_cross_section = splinepy.BSpline( degrees=cross_section.degrees, knot_vectors=cross_section.knot_vectors, control_points=cross_sec_placed_cps[ @@ -485,16 +486,17 @@ def transformation_matrix_with_projection(traj, index, par_value): ) # evaluate cross section at parameter values - evaluations = new_cross_section.sample(7) - cs_evals.append(evaluations) + temp_evaluations = temp_cross_section.sample(7) + cs_evals.append(temp_evaluations) + # defining points for fitting routine fitting_points = _np.array(cs_evals).reshape(-1, 3) # fit surface - take care of size and n_control_points --> not sure yet interpolated_surface, _ = surface( fitting_points=fitting_points, size=[ - len(evaluations), + len(temp_evaluations), nsections, ], n_control_points=[ From fc9b911d0976c396bb2409b91c265f4aa8110e71 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 22 May 2024 21:56:24 +0200 Subject: [PATCH 14/53] Fix: changed arbitrary vector vec and rotation matrix in create-file; changed trajectory in show_swept-file --- examples/show_swept.py | 9 ++++----- splinepy/helpme/create.py | 6 +++--- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index d238b4385..f128ed44e 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -13,8 +13,8 @@ cps1 = np.array( [ [0.0, 0.0, 0.0], - [3.0, 1.0, 0.0], - [4.0, 0.0, 0.0], + [3.0, 1.0, 0.5], + [4.0, 0.0, 2.0], ] ) # init trajectory as bspline @@ -32,15 +32,14 @@ ds_cs = [2] # define knot vectors - kvs_cs = [[0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0]] + kvs_cs = [[0.0, 0.0, 0.0, 1.0, 1.0, 1.0]] # define control points cps_cs = np.array( [ [0.0, 0.0, 0.0], - [1.0, 1.0, 0.0], + [1.0, 2.0, 0.0], [2.0, 0.0, 0.0], - [3.0, 1.0, 0.0], ] ) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 26c25d3f6..ec120b6eb 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -387,7 +387,7 @@ def transformation_matrices(traj, par_value): ) # evaluating arbitrary vector B_0 normal to x - vec = [1.5, 0.3, 1] + vec = [0, 1, 0] B = [] B.append(_np.cross(x, vec) / _np.linalg.norm(_np.cross(x, vec))) @@ -432,8 +432,8 @@ def transformation_matrices(traj, par_value): # ] # )) - # rotation matrix for rotation 90° around z-axis and 90° around y-axis - R = _np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]]) + # rotation matrix for rotation 90° around y-axis + R = _np.array([[0, 0, 1], [0, 1, 0], [1, 0, 0]]) return A, R # compute parameter values for inserting cross sections From 8e6c00669db0cc0b43128bce66d9b8941696d5bf Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Thu, 23 May 2024 12:09:55 +0200 Subject: [PATCH 15/53] Fix: random vector vec now dependent of tangent --- examples/show_swept.py | 2 +- splinepy/helpme/create.py | 11 ++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index f128ed44e..cd90ffe52 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -13,7 +13,7 @@ cps1 = np.array( [ [0.0, 0.0, 0.0], - [3.0, 1.0, 0.5], + [2.0, 1.0, 1.0], [4.0, 0.0, 2.0], ] ) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index ec120b6eb..777f33f44 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -339,7 +339,7 @@ def swept( cross_section_normal = _np.array([0, 0, 1]) # make copies so we can work on it inplace - trajectory = trajectory.copy() + trajectory = trajectory.create.embedded(3) cross_section = cross_section.copy() # compute transformation matrix @@ -382,12 +382,13 @@ def swept( def transformation_matrices(traj, par_value): # tangent vector x on trajectory at parameter value 0 - x = traj.derivative([par_value[0]], [1]) / _np.linalg.norm( + x = ( traj.derivative([par_value[0]], [1]) - ) + / _np.linalg.norm(traj.derivative([par_value[0]], [1])) + ).ravel() - # evaluating arbitrary vector B_0 normal to x - vec = [0, 1, 0] + # evaluating a vector B_0 normal to x + vec = [-x[1], x[0], -x[2]] B = [] B.append(_np.cross(x, vec) / _np.linalg.norm(_np.cross(x, vec))) From fb3ef9563133e2dfa2a8ceef4b6957c2a9215ef6 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Sun, 9 Jun 2024 19:59:41 +0200 Subject: [PATCH 16/53] Add: new function in fit.py-file; Fix: adaptions in create.py-file acc to new func --- examples/show_swept.py | 19 +++++--- splinepy/helpme/create.py | 44 ++++++++++-------- splinepy/helpme/fit.py | 96 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 132 insertions(+), 27 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index cd90ffe52..d63186afd 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -8,31 +8,34 @@ # define degrees ds1 = [2] # define knot vectors - kvs1 = [[0.0, 0.0, 0.0, 1.0, 1.0, 1.0]] + kvs1 = [[0.0, 0.0, 0.0, 0.25, 0.75, 1.0, 1.0, 1.0]] # define control points cps1 = np.array( [ [0.0, 0.0, 0.0], - [2.0, 1.0, 1.0], - [4.0, 0.0, 2.0], + [5.0, 2.0, 4.0], + [10.0, 5.0, 8.0], + [15.0, 2.0, 4.0], + [20.0, 5.0, 0.0], ] ) + # init trajectory as bspline trajectory = splinepy.BSpline( degrees=ds1, knot_vectors=kvs1, - control_points=cps1, - ) + control_points=cps1,) + # define sections along trajectory nsect = len(kvs1[0]) - ds1[0] - 1 # cross section # define degrees - ds_cs = [2] + ds_cs = [3] # define knot vectors - kvs_cs = [[0.0, 0.0, 0.0, 1.0, 1.0, 1.0]] + kvs_cs = [[0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0, 1.0]] # define control points cps_cs = np.array( @@ -40,6 +43,8 @@ [0.0, 0.0, 0.0], [1.0, 2.0, 0.0], [2.0, 0.0, 0.0], + [3.0, -2.0, 0.0], + [4.0, 0.0, 0.0], ] ) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 777f33f44..287f66fd2 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -5,7 +5,7 @@ import splinepy from splinepy import settings as _settings -from splinepy.helpme.fit import surface +from splinepy.helpme.fit import surface_from_cross_sections from splinepy.utils import log as _log from splinepy.utils.data import make_matrix as _make_matrix @@ -459,9 +459,10 @@ def transformation_matrices(traj, par_value): A, R = transformation_matrices(trajectory, par_value) # set cross section evaluation points along trajectory - cps_eval = [] - cs_evals = [] + cps = [] + cross_section_splines = [] for index, eval_point in enumerate(evals): + temp_csp = [] # place every control point of cross section separately for cscp in cross_section.control_points: # rotate cross section in trajectory direction @@ -471,33 +472,36 @@ def transformation_matrices(traj, par_value): # translate cross section to trajectory point normal_cscp = normal_cscp + eval_point # append control point to list - cps_eval.append(normal_cscp) + temp_csp.append(normal_cscp) - # create spline of new cross section - cross_sec_placed_cps = _np.array(cps_eval) - number_of_cps_we_take = len(cross_section.control_points) + # create temp spline of new cross section + cross_sec_placed_cps = _np.array(temp_csp) temp_cross_section = splinepy.BSpline( degrees=cross_section.degrees, knot_vectors=cross_section.knot_vectors, - control_points=cross_sec_placed_cps[ - index - * number_of_cps_we_take : (index + 1) - * number_of_cps_we_take - ], + control_points=cross_sec_placed_cps ) + + # append cross section spline to list + cross_section_splines.append(temp_cross_section) + cps.append(cross_sec_placed_cps) - # evaluate cross section at parameter values - temp_evaluations = temp_cross_section.sample(7) - cs_evals.append(temp_evaluations) + cps = _np.array(cps) + degrees=[int(cross_section.degrees), int(trajectory.degrees)] + knot_vectors=[cross_section.knot_vectors[0][:], trajectory.knot_vectors[0][:]] - # defining points for fitting routine - fitting_points = _np.array(cs_evals).reshape(-1, 3) + fitting_surface = splinepy.BSpline( + degrees, + knot_vectors, + control_points=_np.array(cps).reshape(-1,3),) # fit surface - take care of size and n_control_points --> not sure yet - interpolated_surface, _ = surface( - fitting_points=fitting_points, + interpolated_surface, _ = surface_from_cross_sections( + fitting_spline=fitting_surface, + cross_section_splines=cross_section_splines, + cps=cps, size=[ - len(temp_evaluations), + len(cross_section.control_points), nsections, ], n_control_points=[ diff --git a/splinepy/helpme/fit.py b/splinepy/helpme/fit.py index 3450dd96a..01206153a 100644 --- a/splinepy/helpme/fit.py +++ b/splinepy/helpme/fit.py @@ -5,6 +5,7 @@ from splinepy.utils import log as _log from splinepy.utils.data import has_scipy as _has_scipy from splinepy.utils.data import make_matrix as _make_matrix +import splinepy if _has_scipy: from scipy.sparse.linalg import spsolve as _spsolve @@ -589,3 +590,98 @@ def surface( "residual": residual, } return fitted_spline, residual + +def surface_from_cross_sections( + cross_section_splines, + cps, + size, + degrees=None, + n_control_points=None, + knot_vectors=None, + fitting_spline=None, + associated_queries=None, + centripetal=True, + interpolate_endpoints=True, + verbose_output=False, +): + """ + ... + + Parameters + ---------- + ... + Returns + ------- + ... + """ + + # dim = fitting_points.shape[1] + dim = 3 + + # get evaluation locations + u_k = [None, None] + if associated_queries is not None and ( + knot_vectors is not None or fitting_spline is not None + ): + # check dimensions of queries + if ( + associated_queries[0].shape[1] != 1 + or associated_queries[1].shape[1] != 1 + ): + raise ValueError( + "Associated queries in each direction must have dimension 1!" + ) + u_k[0] = _np.asanyarray(associated_queries[0]) + u_k[1] = _np.asanyarray(associated_queries[1]) + + if fitting_spline.para_dim != 2: + raise ValueError( + "If fitting spline is given as one spline, " + "parametric dimension must be 2!" + ) + + # extract spline in each direction (knot_vectors, weights etc. + # of original spline), works for all type of splines! + fitting_splines = fitting_spline.extract.boundaries([2, 0]) + fitted_spline = fitting_spline.copy() + + n_control_points = fitted_spline.control_mesh_resolutions + + # index helper + mi_interim_cps = MultiIndex((n_control_points[0], size[1])) + + # loop first dim + # curve fit for every j in n_points_v + interim_control_points = _np.empty((n_control_points[0] * size[1], dim)) + # residual_u = _np.empty(size[1]) + # fitted_spline_u = fitting_splines[0] + for v in range(len(cross_section_splines)): + interim_control_points[mi_interim_cps[:, v]] = ( + cross_section_splines[v].control_points + ) + + # loop trajectory dim + # curve fit for every k in n_control_points_u + residual_v = _np.empty(n_control_points[0]) + fitted_spline_v = fitting_splines[1] + for u in range(n_control_points[0]): + fitted_spline_v, residual_v[u] = curve( + fitting_points=cps[:,u], + fitting_spline=fitted_spline_v, + #associated_queries=u_k[1], + centripetal=centripetal, + interpolate_endpoints=interpolate_endpoints, + ) + + interim_control_points[mi_interim_cps[u, : n_control_points[1]]] = ( + fitted_spline_v.control_points + ) + + # copy + fitted_spline.control_points = interim_control_points[ + mi_interim_cps[: n_control_points[0], : n_control_points[1]] + ] + + residual = _np.linalg.norm(residual_v) + + return fitted_spline, residual From 17b9cc9ae1479c442710966a10b43bf20bd333d8 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Sun, 9 Jun 2024 20:04:04 +0200 Subject: [PATCH 17/53] Fix: pre-commit changes --- examples/show_swept.py | 4 ++-- splinepy/helpme/create.py | 20 ++++++++++++-------- splinepy/helpme/fit.py | 20 ++++++++++---------- 3 files changed, 24 insertions(+), 20 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index d63186afd..4c8f131e4 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -24,8 +24,8 @@ trajectory = splinepy.BSpline( degrees=ds1, knot_vectors=kvs1, - control_points=cps1,) - + control_points=cps1, + ) # define sections along trajectory nsect = len(kvs1[0]) - ds1[0] - 1 diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 287f66fd2..81c309155 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -479,21 +479,25 @@ def transformation_matrices(traj, par_value): temp_cross_section = splinepy.BSpline( degrees=cross_section.degrees, knot_vectors=cross_section.knot_vectors, - control_points=cross_sec_placed_cps + control_points=cross_sec_placed_cps, ) - + # append cross section spline to list cross_section_splines.append(temp_cross_section) cps.append(cross_sec_placed_cps) cps = _np.array(cps) - degrees=[int(cross_section.degrees), int(trajectory.degrees)] - knot_vectors=[cross_section.knot_vectors[0][:], trajectory.knot_vectors[0][:]] + degrees = [int(cross_section.degrees), int(trajectory.degrees)] + knot_vectors = [ + cross_section.knot_vectors[0][:], + trajectory.knot_vectors[0][:], + ] fitting_surface = splinepy.BSpline( - degrees, - knot_vectors, - control_points=_np.array(cps).reshape(-1,3),) + degrees, + knot_vectors, + control_points=_np.array(cps).reshape(-1, 3), + ) # fit surface - take care of size and n_control_points --> not sure yet interpolated_surface, _ = surface_from_cross_sections( @@ -508,7 +512,7 @@ def transformation_matrices(traj, par_value): len(cross_section.control_points), len(trajectory.control_points), ], - degrees=[int(cross_section.degrees), int(trajectory.degrees)], + # degrees=[int(cross_section.degrees), int(trajectory.degrees)], ) return interpolated_surface diff --git a/splinepy/helpme/fit.py b/splinepy/helpme/fit.py index 01206153a..2c6f59cf0 100644 --- a/splinepy/helpme/fit.py +++ b/splinepy/helpme/fit.py @@ -5,7 +5,6 @@ from splinepy.utils import log as _log from splinepy.utils.data import has_scipy as _has_scipy from splinepy.utils.data import make_matrix as _make_matrix -import splinepy if _has_scipy: from scipy.sparse.linalg import spsolve as _spsolve @@ -591,18 +590,19 @@ def surface( } return fitted_spline, residual + def surface_from_cross_sections( cross_section_splines, cps, size, - degrees=None, + # degrees=None, n_control_points=None, knot_vectors=None, fitting_spline=None, associated_queries=None, centripetal=True, interpolate_endpoints=True, - verbose_output=False, + # verbose_output=False, ): """ ... @@ -614,7 +614,7 @@ def surface_from_cross_sections( ------- ... """ - + # dim = fitting_points.shape[1] dim = 3 @@ -639,7 +639,7 @@ def surface_from_cross_sections( "If fitting spline is given as one spline, " "parametric dimension must be 2!" ) - + # extract spline in each direction (knot_vectors, weights etc. # of original spline), works for all type of splines! fitting_splines = fitting_spline.extract.boundaries([2, 0]) @@ -656,9 +656,9 @@ def surface_from_cross_sections( # residual_u = _np.empty(size[1]) # fitted_spline_u = fitting_splines[0] for v in range(len(cross_section_splines)): - interim_control_points[mi_interim_cps[:, v]] = ( - cross_section_splines[v].control_points - ) + interim_control_points[mi_interim_cps[:, v]] = cross_section_splines[ + v + ].control_points # loop trajectory dim # curve fit for every k in n_control_points_u @@ -666,9 +666,9 @@ def surface_from_cross_sections( fitted_spline_v = fitting_splines[1] for u in range(n_control_points[0]): fitted_spline_v, residual_v[u] = curve( - fitting_points=cps[:,u], + fitting_points=cps[:, u], fitting_spline=fitted_spline_v, - #associated_queries=u_k[1], + # associated_queries=u_k[1], centripetal=centripetal, interpolate_endpoints=interpolate_endpoints, ) From 3c82d3e040abf2c9a3e029d238578f1cbf449998 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 12 Jun 2024 11:21:52 +0200 Subject: [PATCH 18/53] Add: changed show-settings to better possibility of comparison --- examples/show_swept.py | 7 +++++-- splinepy/helpme/create.py | 2 +- splinepy/helpme/fit.py | 1 - 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 4c8f131e4..fde1869e2 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -59,7 +59,7 @@ cs_nv = np.array([0, 0, 1]) # create interpolated surface - interpolated_surface = splinepy.helpme.create.swept( + interpolated_surface, fitting_surface = splinepy.helpme.create.swept( trajectory=trajectory, cross_section=cross_section, nsections=nsect, @@ -84,6 +84,7 @@ def der(data, on): trajectory.show_options["control_mesh"] = False cross_section.show_options["control_mesh"] = False interpolated_surface.show_options["control_mesh"] = False + fitting_surface.show_options["control_mesh"] = False gus.show( ["Trajectory", trajectory], @@ -92,6 +93,8 @@ def der(data, on): ) gus.show( - ["Surface", interpolated_surface], + ["Trajectory", trajectory], + ["Fitting Surface", fitting_surface], + ["Fitted Surface", interpolated_surface], resolution=50, ) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 81c309155..2f996c0ec 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -515,7 +515,7 @@ def transformation_matrices(traj, par_value): # degrees=[int(cross_section.degrees), int(trajectory.degrees)], ) - return interpolated_surface + return interpolated_surface, fitting_surface def from_bounds(parametric_bounds, physical_bounds): diff --git a/splinepy/helpme/fit.py b/splinepy/helpme/fit.py index 2c6f59cf0..8c6002e97 100644 --- a/splinepy/helpme/fit.py +++ b/splinepy/helpme/fit.py @@ -649,7 +649,6 @@ def surface_from_cross_sections( # index helper mi_interim_cps = MultiIndex((n_control_points[0], size[1])) - # loop first dim # curve fit for every j in n_points_v interim_control_points = _np.empty((n_control_points[0] * size[1], dim)) From 142d2265f58856a43fd413f395ae4c2a5c147f53 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 12 Jun 2024 17:07:23 +0200 Subject: [PATCH 19/53] Big update; Add: nurbs, traj refinement; Del: fitting routine --- examples/show_swept.py | 120 ++++++++++---------------- splinepy/helpme/create.py | 171 +++++++++++++++----------------------- splinepy/helpme/fit.py | 95 --------------------- 3 files changed, 113 insertions(+), 273 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index fde1869e2..9a0d9cc55 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -4,97 +4,69 @@ import splinepy if __name__ == "__main__": - # trajectory - # define degrees - ds1 = [2] - # define knot vectors - kvs1 = [[0.0, 0.0, 0.0, 0.25, 0.75, 1.0, 1.0, 1.0]] - # define control points - cps1 = np.array( - [ - [0.0, 0.0, 0.0], - [5.0, 2.0, 4.0], - [10.0, 5.0, 8.0], - [15.0, 2.0, 4.0], - [20.0, 5.0, 0.0], - ] - ) - - # init trajectory as bspline - trajectory = splinepy.BSpline( - degrees=ds1, - knot_vectors=kvs1, - control_points=cps1, - ) - # define sections along trajectory - nsect = len(kvs1[0]) - ds1[0] - 1 + ### 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], + [5.0, 5.0, 0.0], + [10.0, 7.5, 0.0], + [10.0, 5.0, 0.0], + [10.0, 0.0, 0.0], + ] + ), + } - # cross section - # define degrees - ds_cs = [3] + # init trajectory as bspline + trajectory = splinepy.BSpline(**dict_trajectory) - # define knot vectors - kvs_cs = [[0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0, 1.0]] + # insert knots and control points + # trajectory.uniform_refine([0], 1) - # define control points - cps_cs = np.array( - [ - [0.0, 0.0, 0.0], - [1.0, 2.0, 0.0], - [2.0, 0.0, 0.0], - [3.0, -2.0, 0.0], - [4.0, 0.0, 0.0], - ] - ) + ### CROSS SECTION ### + dict_cross_section = { + "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, 0.0], + [1.0, 2.0, 0.0], + [2.0, 0.0, 0.0], + [3.0, -2.0, 0.0], + [4.0, 0.0, 0.0], + ] + ), + } # init cross section as bspline - cross_section = splinepy.BSpline( - degrees=ds_cs, - knot_vectors=kvs_cs, - control_points=cps_cs, - ) + cross_section = splinepy.BSpline(**dict_cross_section) + + # alternatively, use helpme to create a cross section + # cross_section = splinepy.helpme.create.surface_circle(.5).nurbs - # user should define the normal vector of the cross section + # user can define the normal vector of the cross section, in case + # the cross section is not planar in the x-y plane (default) cs_nv = np.array([0, 0, 1]) - # create interpolated surface - interpolated_surface, fitting_surface = splinepy.helpme.create.swept( + ### SWEEP ### + swept_surface = splinepy.helpme.create.swept( trajectory=trajectory, cross_section=cross_section, - nsections=nsect, + nsections=20, cross_section_normal=cs_nv, ) - # testing derivative and evaluation - # derivative = trajectory.derivative([[0.5], [1]], [1]) - # evals = trajectory.evaluate([[0.5]]) - - # show - def der(data, on): - return data.derivative(on, [1]) - - trajectory.spline_data["der"] = splinepy.SplineDataAdaptor( - trajectory, function=der - ) - trajectory.show_options["arrow_data"] = "der" - trajectory.show_options["arrow_data_on"] = np.linspace(0, 1, 20).reshape( - -1, 1 - ) + ### VISUALIZATION ### trajectory.show_options["control_mesh"] = False cross_section.show_options["control_mesh"] = False - interpolated_surface.show_options["control_mesh"] = False - fitting_surface.show_options["control_mesh"] = False - - gus.show( - ["Trajectory", trajectory], - ["Cross section", cross_section], - resolution=50, - ) + swept_surface.show_options["control_mesh"] = False gus.show( ["Trajectory", trajectory], - ["Fitting Surface", fitting_surface], - ["Fitted Surface", interpolated_surface], - resolution=50, + ["Cross Section", cross_section], + ["Swept Surface", swept_surface], + resolution=51, ) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 2f996c0ec..0c28a531e 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -5,7 +5,6 @@ import splinepy from splinepy import settings as _settings -from splinepy.helpme.fit import surface_from_cross_sections from splinepy.utils import log as _log from splinepy.utils.data import make_matrix as _make_matrix @@ -328,7 +327,7 @@ def swept( from splinepy.spline import Spline as _Spline - # Check input type + # check input type if not isinstance(cross_section, _Spline): raise NotImplementedError("Sweep only works for splines") if not isinstance(trajectory, _Spline): @@ -340,45 +339,7 @@ def swept( # make copies so we can work on it inplace trajectory = trajectory.create.embedded(3) - cross_section = cross_section.copy() - - # compute transformation matrix - # parameters: trajectory traj and parametric value v - # def transformation_matrix_without_projection(traj, v): - # # local directions in global coordinates - # x = traj.derivative([[v]], [1]) / _np.linalg.norm( - # traj.derivative([[v]], [1]) - # ) - # z = _np.cross( - # traj.derivative([[v]], [1]), traj.derivative([[v]], [2]) - # ) / _np.linalg.norm( - # _np.cross(traj.derivative([[v]], [1]), traj.derivative([[v]], [2]))) - - # # check if z is pointing in the right direction - # if z[0][2] < 0: - # z = -z - # y = _np.cross(z, x) - # # transformation matrix from global to local coordinates - # T = [] - # T = _np.vstack((x, y, z)) - - # # transformation matrix from local to global coordinates - # A = _np.linalg.inv(T) - - # # rotation matrix around y - # angle_of_x = _np.arctan2(x[0][2], x[0][0]) - # angle_of_cs_normal = _np.arctan2( - # cross_section_normal[2], cross_section_normal[0] - # ) - # rotation_angle = angle_of_cs_normal - angle_of_x - # R = _np.array( - # [ - # [_np.cos(rotation_angle), 0, _np.sin(rotation_angle)], - # [0, 1, 0], - # [-_np.sin(rotation_angle), 0, _np.cos(rotation_angle)], - # ] - # ) - # return A, R + cross_section = cross_section.create.embedded(3) def transformation_matrices(traj, par_value): # tangent vector x on trajectory at parameter value 0 @@ -392,7 +353,7 @@ def transformation_matrices(traj, par_value): B = [] B.append(_np.cross(x, vec) / _np.linalg.norm(_np.cross(x, vec))) - # initializing transformation matrices + # initialize transformation matrices T = [] A = [] @@ -419,30 +380,39 @@ def transformation_matrices(traj, par_value): A.append(_np.linalg.inv(T[i])) # rotation matrix around y - # angle_of_x = _np.arctan2(x[2], x[0]) - # angle_of_cs_normal = _np.arctan2( - # cross_section_normal[2], cross_section_normal[0] - # ) - # rotation_angle = angle_of_cs_normal - angle_of_x - # R = [] - # R.append(_np.array( - # [ - # [_np.cos(rotation_angle), 0, _np.sin(rotation_angle)], - # [0, 1, 0], - # [-_np.sin(rotation_angle), 0, _np.cos(rotation_angle)], - # ] - # )) - - # rotation matrix for rotation 90° around y-axis - R = _np.array([[0, 0, 1], [0, 1, 0], [1, 0, 0]]) + angle_of_x = _np.arctan2(x[2], x[0]) + angle_of_cs_normal = _np.arctan2( + cross_section_normal[2], cross_section_normal[0] + ) + rotation_angle = angle_of_cs_normal - angle_of_x + R = _np.array( + [ + [_np.cos(rotation_angle), 0, _np.sin(rotation_angle)], + [0, 1, 0], + [_np.sin(rotation_angle), 0, _np.cos(rotation_angle)], + ] + ) + R = _np.where(_np.abs(R) < 1e-10, 0, R) + return A, R # compute parameter values for inserting cross sections if nsections is None: - par_value = trajectory.greville_abscssae() + par_value = trajectory.greville_abscissae() + v_dir_kv = trajectory.knot_vectors[0] else: bounds = trajectory.parametric_bounds.ravel() - par_value = _np.linspace(*bounds, nsections) + par_value = _np.linspace(*bounds, nsections + 1) + # compute new knot vector for traj-direction of swept spline + traj_deg = trajectory.degrees[0] + v_dir_kv = _np.empty(traj_deg + len(par_value) + 1) + v_dir_kv[: (traj_deg + 1)] = 0.0 + v_dir_kv[-(traj_deg + 1) :] = 1.0 + v_dir_kv[(traj_deg + 1) : -(traj_deg + 1)] = ( + _np.convolve(par_value.ravel(), _np.ones(traj_deg), "valid")[1:-1] + / traj_deg + ) + par_value = par_value.reshape(-1, 1) # evaluate trajectory at these parameter values @@ -451,16 +421,15 @@ def transformation_matrices(traj, par_value): # evaluate center of cross section and translate to origin cross_para_center = _np.mean(cross_section.parametric_bounds, axis=0) cs_center = cross_section.evaluate( - cross_para_center.reshape(-1, 1) + cross_para_center.reshape(-1, cross_section.para_dim) ).ravel() cross_section.control_points -= cs_center # evaluate transformation matrices for every trajectory point A, R = transformation_matrices(trajectory, par_value) - # set cross section evaluation points along trajectory - cps = [] - cross_section_splines = [] + # set cross section control points along trajectory + swept_spline_cps = [] for index, eval_point in enumerate(evals): temp_csp = [] # place every control point of cross section separately @@ -470,52 +439,46 @@ def transformation_matrices(traj, par_value): # transform cross section to global coordinates normal_cscp = _np.matmul(A[index], normal_cscp) # translate cross section to trajectory point - normal_cscp = normal_cscp + eval_point + normal_cscp += eval_point # append control point to list temp_csp.append(normal_cscp) - # create temp spline of new cross section - cross_sec_placed_cps = _np.array(temp_csp) - temp_cross_section = splinepy.BSpline( - degrees=cross_section.degrees, - knot_vectors=cross_section.knot_vectors, - control_points=cross_sec_placed_cps, - ) - - # append cross section spline to list - cross_section_splines.append(temp_cross_section) - cps.append(cross_sec_placed_cps) - - cps = _np.array(cps) - degrees = [int(cross_section.degrees), int(trajectory.degrees)] - knot_vectors = [ - cross_section.knot_vectors[0][:], - trajectory.knot_vectors[0][:], - ] - - fitting_surface = splinepy.BSpline( - degrees, - knot_vectors, - control_points=_np.array(cps).reshape(-1, 3), - ) + # collect all control points + swept_spline_cps.append(_np.array(temp_csp)) - # fit surface - take care of size and n_control_points --> not sure yet - interpolated_surface, _ = surface_from_cross_sections( - fitting_spline=fitting_surface, - cross_section_splines=cross_section_splines, - cps=cps, - size=[ - len(cross_section.control_points), - nsections, + # create spline dictionary + dict_swept_spline = { + "degrees": [*cross_section.degrees, *trajectory.degrees], + "knot_vectors": [ + *cross_section.knot_vectors, + v_dir_kv, ], - n_control_points=[ - len(cross_section.control_points), - len(trajectory.control_points), - ], - # degrees=[int(cross_section.degrees), int(trajectory.degrees)], - ) + "control_points": _np.asarray(swept_spline_cps).reshape( + -1, cross_section.dim + ), + } + + # check if spline is rational + if cross_section.is_rational or trajectory.is_rational: + + def weights(spline): + if spline.is_rational: + return spline.weights + return _np.ones(spline.control_points.shape[0]) + + trajectory_weights = weights(trajectory) + cross_section_weights = weights(cross_section) + dict_swept_spline["weights"] = _np.outer( + trajectory_weights, cross_section_weights + ).reshape(-1, 1) + spline_type = splinepy.NURBS + else: + spline_type = splinepy.BSpline + + # create swept spline + swept_spline = spline_type(**dict_swept_spline) - return interpolated_surface, fitting_surface + return swept_spline def from_bounds(parametric_bounds, physical_bounds): diff --git a/splinepy/helpme/fit.py b/splinepy/helpme/fit.py index 8c6002e97..3450dd96a 100644 --- a/splinepy/helpme/fit.py +++ b/splinepy/helpme/fit.py @@ -589,98 +589,3 @@ def surface( "residual": residual, } return fitted_spline, residual - - -def surface_from_cross_sections( - cross_section_splines, - cps, - size, - # degrees=None, - n_control_points=None, - knot_vectors=None, - fitting_spline=None, - associated_queries=None, - centripetal=True, - interpolate_endpoints=True, - # verbose_output=False, -): - """ - ... - - Parameters - ---------- - ... - Returns - ------- - ... - """ - - # dim = fitting_points.shape[1] - dim = 3 - - # get evaluation locations - u_k = [None, None] - if associated_queries is not None and ( - knot_vectors is not None or fitting_spline is not None - ): - # check dimensions of queries - if ( - associated_queries[0].shape[1] != 1 - or associated_queries[1].shape[1] != 1 - ): - raise ValueError( - "Associated queries in each direction must have dimension 1!" - ) - u_k[0] = _np.asanyarray(associated_queries[0]) - u_k[1] = _np.asanyarray(associated_queries[1]) - - if fitting_spline.para_dim != 2: - raise ValueError( - "If fitting spline is given as one spline, " - "parametric dimension must be 2!" - ) - - # extract spline in each direction (knot_vectors, weights etc. - # of original spline), works for all type of splines! - fitting_splines = fitting_spline.extract.boundaries([2, 0]) - fitted_spline = fitting_spline.copy() - - n_control_points = fitted_spline.control_mesh_resolutions - - # index helper - mi_interim_cps = MultiIndex((n_control_points[0], size[1])) - # loop first dim - # curve fit for every j in n_points_v - interim_control_points = _np.empty((n_control_points[0] * size[1], dim)) - # residual_u = _np.empty(size[1]) - # fitted_spline_u = fitting_splines[0] - for v in range(len(cross_section_splines)): - interim_control_points[mi_interim_cps[:, v]] = cross_section_splines[ - v - ].control_points - - # loop trajectory dim - # curve fit for every k in n_control_points_u - residual_v = _np.empty(n_control_points[0]) - fitted_spline_v = fitting_splines[1] - for u in range(n_control_points[0]): - fitted_spline_v, residual_v[u] = curve( - fitting_points=cps[:, u], - fitting_spline=fitted_spline_v, - # associated_queries=u_k[1], - centripetal=centripetal, - interpolate_endpoints=interpolate_endpoints, - ) - - interim_control_points[mi_interim_cps[u, : n_control_points[1]]] = ( - fitted_spline_v.control_points - ) - - # copy - fitted_spline.control_points = interim_control_points[ - mi_interim_cps[: n_control_points[0], : n_control_points[1]] - ] - - residual = _np.linalg.norm(residual_v) - - return fitted_spline, residual From 17b100c339f0d0d72bd9bfcb6f0a7634e947b2a2 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 19 Jun 2024 15:42:31 +0200 Subject: [PATCH 20/53] Fix: calculation of rotation matrix --- examples/show_swept.py | 12 ++++++------ splinepy/helpme/create.py | 24 +++++++++++++++--------- 2 files changed, 21 insertions(+), 15 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 9a0d9cc55..f55996f95 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -12,10 +12,10 @@ "control_points": np.array( [ [0.0, 0.0, 0.0], - [5.0, 5.0, 0.0], - [10.0, 7.5, 0.0], + [0.0, 0.0, 5.0], [10.0, 5.0, 0.0], - [10.0, 0.0, 0.0], + [15.0, 0.0, -5.0], + [20.0, 0.0, 0.0], ] ), } @@ -24,7 +24,7 @@ trajectory = splinepy.BSpline(**dict_trajectory) # insert knots and control points - # trajectory.uniform_refine([0], 1) + trajectory.uniform_refine([0], 3) ### CROSS SECTION ### dict_cross_section = { @@ -55,7 +55,7 @@ swept_surface = splinepy.helpme.create.swept( trajectory=trajectory, cross_section=cross_section, - nsections=20, + nsections=None, cross_section_normal=cs_nv, ) @@ -68,5 +68,5 @@ ["Trajectory", trajectory], ["Cross Section", cross_section], ["Swept Surface", swept_surface], - resolution=51, + resolution=101, ) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 0c28a531e..e890175bf 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -305,7 +305,11 @@ def swept( nsections=None, cross_section_normal=None, ): - """Sweeps a cross-section along a trajectory + """ + Sweeps a cross-section along a trajectory. The cross-section + receives rotation into the direction of the trajectory tangent + vector and is then placed at the evaluation points of the + trajectory's knots. Parameters ---------- @@ -348,10 +352,15 @@ def transformation_matrices(traj, par_value): / _np.linalg.norm(traj.derivative([par_value[0]], [1])) ).ravel() - # evaluating a vector B_0 normal to x + # evaluating a vector normal to x vec = [-x[1], x[0], -x[2]] B = [] - B.append(_np.cross(x, vec) / _np.linalg.norm(_np.cross(x, vec))) + # avoid dividing by zero + if _np.linalg.norm(_np.cross(x, vec)) > 1e-10: + B.append(_np.cross(x, vec) / _np.linalg.norm(_np.cross(x, vec))) + else: + vec = [x[2], -x[1], x[0]] + B.append(_np.cross(x, vec) / _np.linalg.norm(_np.cross(x, vec))) # initialize transformation matrices T = [] @@ -363,7 +372,7 @@ def transformation_matrices(traj, par_value): x = ( traj.derivative([par_value[i]], [1]) / _np.linalg.norm(traj.derivative([par_value[i]], [1])) - ).flatten() + ).ravel() # projecting B_(i-1) onto the plane normal to x B.append(B[i] - _np.dot(B[i], x) * x) @@ -380,20 +389,17 @@ def transformation_matrices(traj, par_value): A.append(_np.linalg.inv(T[i])) # rotation matrix around y - angle_of_x = _np.arctan2(x[2], x[0]) angle_of_cs_normal = _np.arctan2( cross_section_normal[2], cross_section_normal[0] ) - rotation_angle = angle_of_cs_normal - angle_of_x R = _np.array( [ - [_np.cos(rotation_angle), 0, _np.sin(rotation_angle)], + [_np.cos(angle_of_cs_normal), 0, _np.sin(angle_of_cs_normal)], [0, 1, 0], - [_np.sin(rotation_angle), 0, _np.cos(rotation_angle)], + [_np.sin(angle_of_cs_normal), 0, _np.cos(angle_of_cs_normal)], ] ) R = _np.where(_np.abs(R) < 1e-10, 0, R) - return A, R # compute parameter values for inserting cross sections From 42ffcbfbb59702939e049cbae3d9e836d7805109 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Thu, 20 Jun 2024 11:44:18 +0200 Subject: [PATCH 21/53] Add: refinement routine; Fix: no more derivative redundance --- examples/show_swept.py | 3 +- splinepy/helpme/create.py | 77 +++++++++++++++++++++++++-------------- 2 files changed, 50 insertions(+), 30 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index f55996f95..88a1a775e 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -24,7 +24,7 @@ trajectory = splinepy.BSpline(**dict_trajectory) # insert knots and control points - trajectory.uniform_refine([0], 3) + trajectory.uniform_refine([0], 2) ### CROSS SECTION ### dict_cross_section = { @@ -55,7 +55,6 @@ swept_surface = splinepy.helpme.create.swept( trajectory=trajectory, cross_section=cross_section, - nsections=None, cross_section_normal=cs_nv, ) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index e890175bf..ad8a4f9bc 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -302,7 +302,6 @@ def revolved( def swept( cross_section, trajectory, - nsections=None, cross_section_normal=None, ): """ @@ -317,8 +316,6 @@ def swept( Cross-section to be swept trajectory : Spline Trajectory along which the cross-section is swept - nsections : int - Number of sections trajectory is divided into cross_section_normal : np.ndarray Normal vector of the cross-section Default is [0, 0, 1] @@ -345,12 +342,18 @@ def swept( trajectory = trajectory.create.embedded(3) cross_section = cross_section.create.embedded(3) + # print(trajectory.knot_vectors[0]) + # print(len(trajectory.control_points)) + + # # insert knots + # trajectory.insert_knots(0, [0.69, 0.7, 0.1]) + # print(trajectory.knot_vectors[0]) + # print(len(trajectory.control_points)) + def transformation_matrices(traj, par_value): # tangent vector x on trajectory at parameter value 0 - x = ( - traj.derivative([par_value[0]], [1]) - / _np.linalg.norm(traj.derivative([par_value[0]], [1])) - ).ravel() + x = traj.derivative([par_value[0]], [1]) + x = (x / _np.linalg.norm(x)).ravel() # evaluating a vector normal to x vec = [-x[1], x[0], -x[2]] @@ -369,10 +372,8 @@ def transformation_matrices(traj, par_value): # 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]) - / _np.linalg.norm(traj.derivative([par_value[i]], [1])) - ).ravel() + x = traj.derivative([par_value[i]], [1]) + x = (x / _np.linalg.norm(x)).ravel() # projecting B_(i-1) onto the plane normal to x B.append(B[i] - _np.dot(B[i], x) * x) @@ -402,22 +403,42 @@ def transformation_matrices(traj, par_value): R = _np.where(_np.abs(R) < 1e-10, 0, R) return A, R - # compute parameter values for inserting cross sections - if nsections is None: - par_value = trajectory.greville_abscissae() - v_dir_kv = trajectory.knot_vectors[0] - else: - bounds = trajectory.parametric_bounds.ravel() - par_value = _np.linspace(*bounds, nsections + 1) - # compute new knot vector for traj-direction of swept spline - traj_deg = trajectory.degrees[0] - v_dir_kv = _np.empty(traj_deg + len(par_value) + 1) - v_dir_kv[: (traj_deg + 1)] = 0.0 - v_dir_kv[-(traj_deg + 1) :] = 1.0 - v_dir_kv[(traj_deg + 1) : -(traj_deg + 1)] = ( - _np.convolve(par_value.ravel(), _np.ones(traj_deg), "valid")[1:-1] - / traj_deg - ) + par_value = trajectory.greville_abscissae() + par_value = par_value.reshape(-1, 1) + + curv = [] + for i in par_value: + # calculate derivate of trajectory at parametric value i + curv.append(round(_np.linalg.norm(trajectory.derivative([i], [2])), 2)) + + # insert knots and control points in trajectory area with highest curvature + max_curv = max(curv) + max_indices = [i for i, x in enumerate(curv) if x == max_curv] + par_insertion_values = [] + knot_insertion_values = [] + par_valuesss = par_value.ravel() + p = trajectory.degrees[0] + for i in max_indices: + if i > len(trajectory.knot_vectors[0]) - p - 1: + break + else: + par_insertion_values.append( + (par_valuesss[i] + par_valuesss[i + 1]) / 2 + ) + knot_insertion_values.append( + ( + trajectory.knot_vectors[0][i + p] + + trajectory.knot_vectors[0][i + p + 1] + ) + / 2 + ) + trajectory.insert_knots(0, knot_insertion_values) + par_val_reshape = par_value.reshape(-1) + new_insertion_values_reshape = _np.asanyarray(par_insertion_values) + par_value = _np.concatenate( + (par_val_reshape, new_insertion_values_reshape) + ) + par_value = _np.sort(par_value) par_value = par_value.reshape(-1, 1) @@ -457,7 +478,7 @@ def transformation_matrices(traj, par_value): "degrees": [*cross_section.degrees, *trajectory.degrees], "knot_vectors": [ *cross_section.knot_vectors, - v_dir_kv, + *trajectory.knot_vectors, ], "control_points": _np.asarray(swept_spline_cps).reshape( -1, cross_section.dim From 98bde93b5bce72957061f14d025921dd273a23ac Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Thu, 27 Jun 2024 14:50:50 +0200 Subject: [PATCH 22/53] Fix: no more global splinepy import --- examples/show_swept.py | 6 +++--- splinepy/helpme/create.py | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 88a1a775e..b599c1a67 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -22,9 +22,9 @@ # init trajectory as bspline trajectory = splinepy.BSpline(**dict_trajectory) - + trajectory = splinepy.helpme.create.circle(10) # insert knots and control points - trajectory.uniform_refine([0], 2) + trajectory.uniform_refine([0], 4) ### CROSS SECTION ### dict_cross_section = { @@ -45,7 +45,7 @@ cross_section = splinepy.BSpline(**dict_cross_section) # alternatively, use helpme to create a cross section - # cross_section = splinepy.helpme.create.surface_circle(.5).nurbs + cross_section = splinepy.helpme.create.surface_circle(0.5).nurbs # user can define the normal vector of the cross section, in case # the cross section is not planar in the x-y plane (default) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index ad8a4f9bc..d6b52c6c3 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -3,7 +3,6 @@ import numpy as _np from gustaf.utils import arr as _arr -import splinepy from splinepy import settings as _settings from splinepy.utils import log as _log from splinepy.utils.data import make_matrix as _make_matrix @@ -326,6 +325,8 @@ def swept( Spline resulting from the sweep """ + from splinepy import NURBS as _NURBS + from splinepy import BSpline as _BSpline from splinepy.spline import Spline as _Spline # check input type @@ -439,7 +440,6 @@ def transformation_matrices(traj, par_value): (par_val_reshape, new_insertion_values_reshape) ) par_value = _np.sort(par_value) - par_value = par_value.reshape(-1, 1) # evaluate trajectory at these parameter values @@ -498,9 +498,9 @@ def weights(spline): dict_swept_spline["weights"] = _np.outer( trajectory_weights, cross_section_weights ).reshape(-1, 1) - spline_type = splinepy.NURBS + spline_type = _NURBS else: - spline_type = splinepy.BSpline + spline_type = _BSpline # create swept spline swept_spline = spline_type(**dict_swept_spline) From d79096ff56e573abedf2581fe70ce12e1567f31f Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Mon, 1 Jul 2024 12:07:22 +0200 Subject: [PATCH 23/53] Add: self-refinement now working; added function in bspline.py (note: not ideal yet) --- examples/show_swept.py | 6 ++-- splinepy/bspline.py | 58 +++++++++++++++++++++++++++++++++++- splinepy/helpme/create.py | 62 +++++++++++++++++++-------------------- 3 files changed, 90 insertions(+), 36 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index b599c1a67..573ad39a6 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -22,9 +22,9 @@ # init trajectory as bspline trajectory = splinepy.BSpline(**dict_trajectory) - trajectory = splinepy.helpme.create.circle(10) + # trajectory = splinepy.helpme.create.circle(10) # insert knots and control points - trajectory.uniform_refine([0], 4) + trajectory.uniform_refine(0, 2) ### CROSS SECTION ### dict_cross_section = { @@ -45,7 +45,7 @@ cross_section = splinepy.BSpline(**dict_cross_section) # alternatively, use helpme to create a cross section - cross_section = splinepy.helpme.create.surface_circle(0.5).nurbs + # cross_section = splinepy.helpme.create.surface_circle(0.5).nurbs # user can define the normal vector of the cross section, in case # the cross section is not planar in the x-y plane (default) diff --git a/splinepy/bspline.py b/splinepy/bspline.py index e0d6f7320..a71efe4d2 100644 --- a/splinepy/bspline.py +++ b/splinepy/bspline.py @@ -179,6 +179,63 @@ def determine_new_knots(kv_unique, n_knots): ) self.insert_knots(para_dim, new_knots) + def uniform_refine_fixed_knots(self, para_dims=None, total_knots=4): + """ + Uniformly refines the knot vector in given direction(s) by adding a fixed + total number of knots, which are evenly distributed. + + Parameters + ---------- + para_dims : int or list + list of parametric dimensions to be refined (default None -> all) + total_knots : int + total number of new knots to be added + + Returns + -------- + None + """ + # if no para_dim is given - assume that each dimension is refined + if para_dims is None: + para_dims = range(self.para_dim) + + # if an integer is given, make it a list + elif isinstance(para_dims, int): + para_dims = [para_dims] + + def determine_new_knots_fixed(kv_unique, total_knots): + if total_knots == 0: + return [] + # Calculate the number of new knots to be added between each + # pair of existing knots + num_spans = len(kv_unique) - 1 + new_knots_per_span = total_knots // num_spans + additional_knots = total_knots % num_spans + + new_knots = [] + for i in range(num_spans): + span = kv_unique[i + 1] - kv_unique[i] + num_knots_in_span = new_knots_per_span + ( + 1 if i < additional_knots else 0 + ) + if num_knots_in_span > 0: + new_knots_in_span = kv_unique[i] + span * _np.linspace( + 1 / (num_knots_in_span + 1), + 1 - 1 / (num_knots_in_span + 1), + num_knots_in_span, + ) + new_knots.extend(new_knots_in_span) + return _np.array(new_knots) + + # determine new knots for each para_dim and insert the knots + for para_dim in para_dims: + new_knots = determine_new_knots_fixed( + # recompute unique to allow duplicating para_dims. + kv_unique=self.unique_knots[para_dim], + total_knots=total_knots, + ) + self.insert_knots(para_dim, new_knots) + def knot_insertion_matrix( self, parametric_dimension=None, knots=None, beziers=False ): @@ -319,7 +376,6 @@ def remove_knots(self, parametric_dimension, knots, tolerance=None): self._logd(f"Tried to remove {len(knots)} knot(s).") self._logd(f"Actually removed {sum(removed)} knot(s).") - return removed def normalize_knot_vectors(self): diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index d6b52c6c3..9a66502a9 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -311,7 +311,7 @@ def swept( Parameters ---------- - crossection : Spline + cross_section : Spline Cross-section to be swept trajectory : Spline Trajectory along which the cross-section is swept @@ -343,14 +343,6 @@ def swept( trajectory = trajectory.create.embedded(3) cross_section = cross_section.create.embedded(3) - # print(trajectory.knot_vectors[0]) - # print(len(trajectory.control_points)) - - # # insert knots - # trajectory.insert_knots(0, [0.69, 0.7, 0.1]) - # print(trajectory.knot_vectors[0]) - # print(len(trajectory.control_points)) - def transformation_matrices(traj, par_value): # tangent vector x on trajectory at parameter value 0 x = traj.derivative([par_value[0]], [1]) @@ -407,42 +399,48 @@ def transformation_matrices(traj, par_value): par_value = trajectory.greville_abscissae() par_value = par_value.reshape(-1, 1) + ### insert knots in trajectory area with highest curvature ### curv = [] for i in par_value: - # calculate derivate of trajectory at parametric value i + # calculate curvature of trajectory at parametric value i curv.append(round(_np.linalg.norm(trajectory.derivative([i], [2])), 2)) - - # insert knots and control points in trajectory area with highest curvature + # evaluate the par_values-vector indices of the maximum curvature points max_curv = max(curv) max_indices = [i for i, x in enumerate(curv) if x == max_curv] - par_insertion_values = [] - knot_insertion_values = [] - par_valuesss = par_value.ravel() - p = trajectory.degrees[0] - for i in max_indices: - if i > len(trajectory.knot_vectors[0]) - p - 1: - break + # prepare two matrices for the insertion + insertion_values = [] + # compute the new insertion values + par_values = par_value.ravel() + + for maxi in max_indices: + if maxi == 0: + insertion_values.append( + (par_values[maxi] + par_values[maxi + 1]) / 2 + ) + elif maxi == len(par_values) - 1: + insertion_values.append( + (par_values[maxi] + par_values[maxi - 1]) / 2 + ) else: - par_insertion_values.append( - (par_valuesss[i] + par_valuesss[i + 1]) / 2 + insertion_values.append( + (par_values[maxi] + par_values[maxi - 1]) / 2 ) - knot_insertion_values.append( - ( - trajectory.knot_vectors[0][i + p] - + trajectory.knot_vectors[0][i + p + 1] - ) - / 2 + insertion_values.append( + (par_values[maxi] + par_values[maxi + 1]) / 2 ) - trajectory.insert_knots(0, knot_insertion_values) - par_val_reshape = par_value.reshape(-1) - new_insertion_values_reshape = _np.asanyarray(par_insertion_values) + + # insert knots into the trajectory's knot vector + insertion_values = _np.unique(insertion_values) + trajectory.uniform_refine_fixed_knots(0, len(insertion_values)) + # insert knots into the parameter values par_value = _np.concatenate( - (par_val_reshape, new_insertion_values_reshape) + (par_value.reshape(-1), _np.asanyarray(insertion_values)) ) + # sort parameter values par_value = _np.sort(par_value) par_value = par_value.reshape(-1, 1) - # evaluate trajectory at these parameter values + # evaluate trajectory at the parameter values evals = trajectory.evaluate(par_value) # evaluate center of cross section and translate to origin From 35961fcba5f6a6fbae55f7a4b94720fe67e9aaf3 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Tue, 2 Jul 2024 10:58:59 +0200 Subject: [PATCH 24/53] Add: function for closed trajectories --- examples/show_swept.py | 131 +++++++++++++++++++++++++++++++++++--- splinepy/helpme/create.py | 62 ++++++++++++++++-- 2 files changed, 179 insertions(+), 14 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 573ad39a6..9e3a82e50 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -1,3 +1,5 @@ +import sys + import gustaf as gus import numpy as np @@ -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 = { @@ -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) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 9a66502a9..f17001788 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -308,6 +308,8 @@ 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 ---------- @@ -315,7 +317,7 @@ def swept( 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] @@ -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() @@ -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]) @@ -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] ) @@ -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) @@ -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) From f08140b1ffebb5e98ef4f71c72db7ba5c37159ab Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Tue, 2 Jul 2024 13:18:55 +0200 Subject: [PATCH 25/53] Add: some checks were added --- splinepy/helpme/create.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index f17001788..912aa6779 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -336,6 +336,10 @@ def swept( raise NotImplementedError("Sweep only works for splines") if not isinstance(trajectory, _Spline): raise NotImplementedError("Sweep only works for splines") + if not trajectory.para_dim == 1: + raise NotImplementedError("Trajectory must be 1D") + if not len(cross_section_normal) == 3: + raise ValueError("Cross section normal must be 3D") # setting default value for cross_section_normal if cross_section_normal is None: @@ -459,7 +463,7 @@ def transformation_matrices(traj, par_value): # evaluate the par_values-vector indices of the maximum curvature points max_curv = max(curv) max_indices = [i for i, x in enumerate(curv) if x == max_curv] - # prepare two matrices for the insertion + # prepare matrix for the insertion insertion_values = [] # compute the new insertion values par_values = par_value.ravel() @@ -484,7 +488,7 @@ def transformation_matrices(traj, par_value): # insert knots into the trajectory's knot vector insertion_values = _np.unique(insertion_values) trajectory.uniform_refine_fixed_knots(0, len(insertion_values)) - # insert knots into the parameter values + # add insertion values to the existing parameter values par_value = _np.concatenate( (par_value.reshape(-1), _np.asanyarray(insertion_values)) ) @@ -551,6 +555,7 @@ def weights(spline): trajectory_weights, cross_section_weights ).reshape(-1, 1) spline_type = _NURBS + else: spline_type = _BSpline From 9a5c99d4da6e8232d06ea9cbdd2aa35f47d47a94 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Tue, 2 Jul 2024 13:23:48 +0200 Subject: [PATCH 26/53] Add: reference to nurbs book --- splinepy/helpme/create.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 912aa6779..2f90750c5 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -374,6 +374,7 @@ def transformation_matrices(traj, par_value): # evaluating transformation matrices for each trajectory point for i in range(len(par_value)): + # calculation according to NURBS Book, eq. 10.27 # tangent vector x on trajectory at parameter value i x = traj.derivative([par_value[i]], [1]) x = (x / _np.linalg.norm(x)).ravel() @@ -405,6 +406,7 @@ def transformation_matrices(traj, par_value): B_reverse[0] = B[-1] for i in range(len(par_value)): # redo the calculation of B using x_collection from before + # according to NURBS Book, eq. 10.27 B_reverse[i + 1] = ( B_reverse[i] - _np.dot(B_reverse[i], x_collection[i]) * x_collection[i] From e1d46a9aebbfb9fc72b33d036b46f3bd6f8046a9 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Tue, 2 Jul 2024 13:36:38 +0200 Subject: [PATCH 27/53] Add: auto_refinement variable --- examples/show_swept.py | 1 + splinepy/helpme/create.py | 91 ++++++++++++++++++++++----------------- 2 files changed, 52 insertions(+), 40 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 9e3a82e50..e4007169e 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -120,6 +120,7 @@ trajectory=trajectory, cross_section=cross_section, cross_section_normal=cs_nv, + auto_refinement=True, ) ### VISUALIZATION ### diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 2f90750c5..ed0fa2639 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -302,6 +302,7 @@ def swept( cross_section, trajectory, cross_section_normal=None, + auto_refinement=False, ): """ Sweeps a cross-section along a trajectory. The cross-section @@ -320,6 +321,9 @@ def swept( cross_section_normal : np.array Normal vector of the cross-section Default is [0, 0, 1] + auto_refinement : bool + If True, the trajectory will be refined at points of + highest curvature. Default is False. Returns ------- @@ -340,6 +344,8 @@ def swept( raise NotImplementedError("Trajectory must be 1D") if not len(cross_section_normal) == 3: raise ValueError("Cross section normal must be 3D") + if not isinstance(auto_refinement, bool): + raise ValueError("auto_refinement must be a boolean") # setting default value for cross_section_normal if cross_section_normal is None: @@ -349,6 +355,10 @@ def swept( trajectory = trajectory.create.embedded(3) cross_section = cross_section.create.embedded(3) + # initialize parameter values + par_value = trajectory.greville_abscissae() + par_value = par_value.reshape(-1, 1) + def transformation_matrices(traj, par_value): ### TRANSFORMATION MATRICES ### @@ -454,49 +464,50 @@ def transformation_matrices(traj, par_value): ### REFINEMENT OF TRAJECTORY ### - par_value = trajectory.greville_abscissae() - par_value = par_value.reshape(-1, 1) + if auto_refinement: + ## inserts knots in trajectory area with highest curvature # - ### insert knots in trajectory area with highest curvature ### - curv = [] - for i in par_value: - # calculate curvature of trajectory at parametric value i - curv.append(round(_np.linalg.norm(trajectory.derivative([i], [2])), 2)) - # evaluate the par_values-vector indices of the maximum curvature points - max_curv = max(curv) - max_indices = [i for i, x in enumerate(curv) if x == max_curv] - # prepare matrix for the insertion - insertion_values = [] - # compute the new insertion values - par_values = par_value.ravel() - - for maxi in max_indices: - if maxi == 0: - insertion_values.append( - (par_values[maxi] + par_values[maxi + 1]) / 2 - ) - elif maxi == len(par_values) - 1: - insertion_values.append( - (par_values[maxi] + par_values[maxi - 1]) / 2 - ) - else: - insertion_values.append( - (par_values[maxi] + par_values[maxi - 1]) / 2 - ) - insertion_values.append( - (par_values[maxi] + par_values[maxi + 1]) / 2 + curv = [] + for i in par_value: + # calculate curvature of trajectory at parametric value i + curv.append( + round(_np.linalg.norm(trajectory.derivative([i], [2])), 2) ) + # evaluate the par_values-vector indices of the maximum curvature points + max_curv = max(curv) + max_indices = [i for i, x in enumerate(curv) if x == max_curv] + # prepare matrix for the insertion + insertion_values = [] + # compute the new insertion values + par_values = par_value.ravel() + + for maxi in max_indices: + if maxi == 0: + insertion_values.append( + (par_values[maxi] + par_values[maxi + 1]) / 2 + ) + elif maxi == len(par_values) - 1: + insertion_values.append( + (par_values[maxi] + par_values[maxi - 1]) / 2 + ) + else: + insertion_values.append( + (par_values[maxi] + par_values[maxi - 1]) / 2 + ) + insertion_values.append( + (par_values[maxi] + par_values[maxi + 1]) / 2 + ) - # insert knots into the trajectory's knot vector - insertion_values = _np.unique(insertion_values) - trajectory.uniform_refine_fixed_knots(0, len(insertion_values)) - # add insertion values to the existing parameter values - par_value = _np.concatenate( - (par_value.reshape(-1), _np.asanyarray(insertion_values)) - ) - # sort parameter values - par_value = _np.sort(par_value) - par_value = par_value.reshape(-1, 1) + # insert knots into the trajectory's knot vector + insertion_values = _np.unique(insertion_values) + trajectory.uniform_refine_fixed_knots(0, len(insertion_values)) + # add insertion values to the existing parameter values + par_value = _np.concatenate( + (par_value.reshape(-1), _np.asanyarray(insertion_values)) + ) + # sort parameter values + par_value = _np.sort(par_value) + par_value = par_value.reshape(-1, 1) ### SWEEPING PROCESS ### From 384a0d22e1d69c5f464d36a6eb83546fc4597036 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Tue, 2 Jul 2024 17:23:14 +0200 Subject: [PATCH 28/53] Fix: positioning of cs-cps now acc to traj-cps; Fix: auto_refinement routine --- examples/show_swept.py | 132 ++++++++++++++++++++------------------ splinepy/helpme/create.py | 45 +++++++++---- 2 files changed, 100 insertions(+), 77 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index e4007169e..ba22802a1 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -25,63 +25,67 @@ # } # 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": [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, - ] + [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, 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], + [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": [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.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], + # ] + # ), + # } trajectory = splinepy.BSpline(**dict_trajectory) # alternatively, use helpme to create a trajectory @@ -91,25 +95,25 @@ trajectory.uniform_refine(0, 1) ### CROSS SECTION ### - dict_cross_section = { - "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, 0.0], - [1.0, 2.0, 0.0], - [2.0, 0.0, 0.0], - [3.0, -2.0, 0.0], - [4.0, 0.0, 0.0], - ] - ), - } + # dict_cross_section = { + # "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, 0.0], + # [1.0, 2.0, 0.0], + # [2.0, 0.0, 0.0], + # [3.0, -2.0, 0.0], + # [4.0, 0.0, 0.0], + # ] + # ), + # } - # init cross section as bspline - cross_section = splinepy.BSpline(**dict_cross_section) + # # init cross section as bspline + # cross_section = splinepy.BSpline(**dict_cross_section) # alternatively, use helpme to create a cross section - # cross_section = splinepy.helpme.create.surface_circle(0.5).nurbs + cross_section = splinepy.helpme.create.surface_circle(0.5).nurbs # user can define the normal vector of the cross section, in case # the cross section is not planar in the x-y plane (default) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index ed0fa2639..45b50ec68 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -440,7 +440,7 @@ def transformation_matrices(traj, par_value): 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]): + if not _np.allclose(B_reverse[0], B_reverse[-1], rtol=1e-3): _log.warning( "Vector calculation is not exact due to the " "trajectory being closed in an uncommon way." @@ -500,20 +500,39 @@ def transformation_matrices(traj, par_value): # insert knots into the trajectory's knot vector insertion_values = _np.unique(insertion_values) - trajectory.uniform_refine_fixed_knots(0, len(insertion_values)) - # add insertion values to the existing parameter values - par_value = _np.concatenate( - (par_value.reshape(-1), _np.asanyarray(insertion_values)) - ) - # sort parameter values - par_value = _np.sort(par_value) + # convert knot vector to list + kv_list = trajectory.knot_vectors[0].numpy() + + # insert knots into the trajectory's knot vector + if any(value in insertion_values for value in kv_list): + add = _np.concatenate((insertion_values, kv_list)) + add = _np.unique(add) + # remove the existing knots + add = add[~_np.isin(add, kv_list)] + # check if there are any knots to insert + if len(add) == 0: + _log.warning("Auto Refinement couldn't insert knots.") + else: + trajectory.insert_knots(0, add) + # give information about the inserted knots + _log.info( + f"Auto Refinement inserted {len(add)} " + "knots into the trajectory." + ) + else: + trajectory.insert_knots(0, insertion_values) + # give information about the inserted knots + _log.info( + f"Auto Refinement inserted {len(insertion_values)} " + "knots into the trajectory." + ) + + # recalculate parameter values + par_value = trajectory.greville_abscissae() par_value = par_value.reshape(-1, 1) ### SWEEPING PROCESS ### - # evaluate trajectory at the parameter values - evals = trajectory.evaluate(par_value) - # evaluate center of cross section and translate to origin cross_para_center = _np.mean(cross_section.parametric_bounds, axis=0) cs_center = cross_section.evaluate( @@ -526,7 +545,7 @@ def transformation_matrices(traj, par_value): # set cross section control points along trajectory swept_spline_cps = [] - for index, eval_point in enumerate(evals): + for index in range(len(par_value)): temp_csp = [] # place every control point of cross section separately for cscp in cross_section.control_points: @@ -535,7 +554,7 @@ def transformation_matrices(traj, par_value): # transform cross section to global coordinates normal_cscp = _np.matmul(A[index], normal_cscp) # translate cross section to trajectory point - normal_cscp += eval_point + normal_cscp += trajectory.control_points[index] # append control point to list temp_csp.append(normal_cscp) From 46c8f5473bb235baedb5034a87840f6b7a2cccf4 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Tue, 2 Jul 2024 17:32:50 +0200 Subject: [PATCH 29/53] Add: test functions in create-test-file --- splinepy/helpme/create.py | 3 +- tests/helpme/test_create.py | 76 +++++++++++++++++++++++++++++++++++++ 2 files changed, 78 insertions(+), 1 deletion(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 45b50ec68..f2a11d547 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -342,7 +342,8 @@ def swept( raise NotImplementedError("Sweep only works for splines") if not trajectory.para_dim == 1: raise NotImplementedError("Trajectory must be 1D") - if not len(cross_section_normal) == 3: + + if cross_section_normal is not None and not len(cross_section_normal) == 3: raise ValueError("Cross section normal must be 3D") if not isinstance(auto_refinement, bool): raise ValueError("auto_refinement must be a boolean") diff --git a/tests/helpme/test_create.py b/tests/helpme/test_create.py index 265f3c103..1c79c7b23 100644 --- a/tests/helpme/test_create.py +++ b/tests/helpme/test_create.py @@ -343,3 +343,79 @@ def test_determinant_spline( det_spl.evaluate(queries=rnd_queries).ravel(), np.linalg.det(sp_i.jacobian(queries=rnd_queries)), ), f"{sp_i.whatami} at index {idx} failed determinant spline" + + +def test_swept_basic_functionality(): + cross_section = splinepy.BSpline( + degrees=[2], + control_points=[[0, 0], [0.5, 1], [1, 0]], + knot_vectors=[[0, 0, 0, 1, 1, 1]], + ) + trajectory = splinepy.BSpline( + degrees=[3], + control_points=[[0, 0], [1, 2], [2, 3], [3, 3]], + knot_vectors=[[0, 0, 0, 0, 1, 1, 1, 1]], + ) + result = splinepy.helpme.create.swept(cross_section, trajectory) + assert result is not None + assert ( + result.control_points.shape[0] + == cross_section.control_points.shape[0] + * trajectory.control_points.shape[0] + ) + + +def test_swept_with_custom_normal(): + cross_section = splinepy.BSpline( + degrees=[2], + control_points=[[0, 0], [0.5, 1], [1, 0]], + knot_vectors=[[0, 0, 0, 1, 1, 1]], + ) + trajectory = splinepy.BSpline( + degrees=[3], + control_points=[[0, 0], [1, 2], [2, 3], [3, 3]], + knot_vectors=[[0, 0, 0, 0, 1, 1, 1, 1]], + ) + custom_normal = np.array([0, 1, 0]) + result = splinepy.helpme.create.swept( + cross_section, trajectory, cross_section_normal=custom_normal + ) + assert result is not None + + +def test_swept_invalid_inputs(): + cross_section = splinepy.BSpline( + degrees=[2], + control_points=[[0, 0], [0.5, 1], [1, 0]], + knot_vectors=[[0, 0, 0, 1, 1, 1]], + ) + invalid_trajectory = "invalid_trajectory" + with pytest.raises(NotImplementedError): + splinepy.helpme.create.swept(cross_section, invalid_trajectory) + + invalid_cross_section = "invalid_cross_section" + trajectory = splinepy.BSpline( + degrees=[3], + control_points=[[0, 0], [1, 2], [2, 3], [3, 3]], + knot_vectors=[[0, 0, 0, 0, 1, 1, 1, 1]], + ) + with pytest.raises(NotImplementedError): + splinepy.helpme.create.swept(invalid_cross_section, trajectory) + + +def test_swept_rational_splines(): + cross_section = splinepy.NURBS( + degrees=[2], + control_points=[[0, 0], [0.5, 1], [1, 0]], + weights=[1, 0.5, 1], + knot_vectors=[[0, 0, 0, 1, 1, 1]], + ) + trajectory = splinepy.NURBS( + degrees=[3], + control_points=[[0, 0], [1, 2], [2, 3], [3, 3]], + weights=[1, 0.5, 0.5, 1], + knot_vectors=[[0, 0, 0, 0, 1, 1, 1, 1]], + ) + result = splinepy.helpme.create.swept(cross_section, trajectory) + assert result is not None + assert result.is_rational From 75a5bb61296820dfedce18cfc98cfc752d81b919 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 3 Jul 2024 10:18:15 +0200 Subject: [PATCH 30/53] Rm: removed unnecessary insertion function in bspline.py --- splinepy/bspline.py | 57 --------------------------------------------- 1 file changed, 57 deletions(-) diff --git a/splinepy/bspline.py b/splinepy/bspline.py index a71efe4d2..b932b27df 100644 --- a/splinepy/bspline.py +++ b/splinepy/bspline.py @@ -179,63 +179,6 @@ def determine_new_knots(kv_unique, n_knots): ) self.insert_knots(para_dim, new_knots) - def uniform_refine_fixed_knots(self, para_dims=None, total_knots=4): - """ - Uniformly refines the knot vector in given direction(s) by adding a fixed - total number of knots, which are evenly distributed. - - Parameters - ---------- - para_dims : int or list - list of parametric dimensions to be refined (default None -> all) - total_knots : int - total number of new knots to be added - - Returns - -------- - None - """ - # if no para_dim is given - assume that each dimension is refined - if para_dims is None: - para_dims = range(self.para_dim) - - # if an integer is given, make it a list - elif isinstance(para_dims, int): - para_dims = [para_dims] - - def determine_new_knots_fixed(kv_unique, total_knots): - if total_knots == 0: - return [] - # Calculate the number of new knots to be added between each - # pair of existing knots - num_spans = len(kv_unique) - 1 - new_knots_per_span = total_knots // num_spans - additional_knots = total_knots % num_spans - - new_knots = [] - for i in range(num_spans): - span = kv_unique[i + 1] - kv_unique[i] - num_knots_in_span = new_knots_per_span + ( - 1 if i < additional_knots else 0 - ) - if num_knots_in_span > 0: - new_knots_in_span = kv_unique[i] + span * _np.linspace( - 1 / (num_knots_in_span + 1), - 1 - 1 / (num_knots_in_span + 1), - num_knots_in_span, - ) - new_knots.extend(new_knots_in_span) - return _np.array(new_knots) - - # determine new knots for each para_dim and insert the knots - for para_dim in para_dims: - new_knots = determine_new_knots_fixed( - # recompute unique to allow duplicating para_dims. - kv_unique=self.unique_knots[para_dim], - total_knots=total_knots, - ) - self.insert_knots(para_dim, new_knots) - def knot_insertion_matrix( self, parametric_dimension=None, knots=None, beziers=False ): From 65a320f359246cf755f8447a11b201da6b779df9 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 3 Jul 2024 11:03:23 +0200 Subject: [PATCH 31/53] Add: opportunity to set cross-sec CPs on traj eval points or traj CPs --- examples/show_swept.py | 105 +++++++++++++++++++------------------- splinepy/helpme/create.py | 24 +++++++-- 2 files changed, 72 insertions(+), 57 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index ba22802a1..003899b89 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -25,67 +25,67 @@ # } # 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": [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, - # ] + # [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, 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], + # [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": [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.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], + ] + ), + } trajectory = splinepy.BSpline(**dict_trajectory) # alternatively, use helpme to create a trajectory @@ -125,6 +125,7 @@ cross_section=cross_section, cross_section_normal=cs_nv, auto_refinement=True, + set_on_trajectory=True, ) ### VISUALIZATION ### diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index f2a11d547..0eca124bf 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -303,6 +303,7 @@ def swept( trajectory, cross_section_normal=None, auto_refinement=False, + set_on_trajectory=False, ): """ Sweeps a cross-section along a trajectory. The cross-section @@ -324,6 +325,11 @@ def swept( auto_refinement : bool If True, the trajectory will be refined at points of highest curvature. Default is False. + set_on_trajectory : bool + If True, the cross-section will be placed at the evaluation + points of the trajectory's knots. If False, the cross-section + will be placed at the control points of the trajectory. + Default is False. Returns ------- @@ -466,7 +472,7 @@ def transformation_matrices(traj, par_value): ### REFINEMENT OF TRAJECTORY ### if auto_refinement: - ## inserts knots in trajectory area with highest curvature # + ## inserts knots in trajectory area with highest curvature ## curv = [] for i in par_value: @@ -546,16 +552,24 @@ def transformation_matrices(traj, par_value): # set cross section control points along trajectory swept_spline_cps = [] - for index in range(len(par_value)): + for i, par_val in enumerate(par_value): temp_csp = [] + if set_on_trajectory: + evals = trajectory.evaluate([par_val]).ravel() # place every control point of cross section separately for cscp in cross_section.control_points: # rotate cross section in trajectory direction normal_cscp = _np.matmul(R, cscp) # transform cross section to global coordinates - normal_cscp = _np.matmul(A[index], normal_cscp) - # translate cross section to trajectory point - normal_cscp += trajectory.control_points[index] + normal_cscp = _np.matmul(A[i], normal_cscp) + # check whether user wants to place cross section + # at evaluation points or control points of trajectory + if set_on_trajectory: + # translate cross section to trajectory evaluation point + normal_cscp += evals + else: + # translate cross section to trajectory control point + normal_cscp += trajectory.control_points[i] # append control point to list temp_csp.append(normal_cscp) From 9d34f7c35ea670ba98d87728011c9385e3cbf31e Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 3 Jul 2024 11:38:22 +0200 Subject: [PATCH 32/53] Add: swept description; Fix: using settings.TOLERANCE --- examples/show_swept.py | 144 +++++++++++++++++++------------------- splinepy/helpme/create.py | 15 ++-- 2 files changed, 82 insertions(+), 77 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 003899b89..25212c14c 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -10,19 +10,19 @@ ### 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], - # ] - # ), - # } + 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 = { @@ -46,78 +46,78 @@ # } # init trajectory as bspline - # closed 3D questionmark - dict_trajectory = { - "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.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], - ] - ), - } + # # closed 3D questionmark + # dict_trajectory = { + # "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.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], + # ] + # ), + # } 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, 1) + trajectory.uniform_refine(0, 3) ### CROSS SECTION ### - # dict_cross_section = { - # "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, 0.0], - # [1.0, 2.0, 0.0], - # [2.0, 0.0, 0.0], - # [3.0, -2.0, 0.0], - # [4.0, 0.0, 0.0], - # ] - # ), - # } + dict_cross_section = { + "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, 0.0], + [1.0, 2.0, 0.0], + [2.0, 0.0, 0.0], + [3.0, -2.0, 0.0], + [4.0, 0.0, 0.0], + ] + ), + } - # # init cross section as bspline - # cross_section = splinepy.BSpline(**dict_cross_section) + # init cross section as bspline + cross_section = splinepy.BSpline(**dict_cross_section) # alternatively, use helpme to create a cross section - cross_section = splinepy.helpme.create.surface_circle(0.5).nurbs + # cross_section = splinepy.helpme.create.surface_circle(0.5).nurbs # user can define the normal vector of the cross section, in case # the cross section is not planar in the x-y plane (default) - cs_nv = np.array([0, 0, 1]) + cs_nv = np.array([0, 1, 1]) ### SWEEP ### swept_surface = splinepy.helpme.create.swept( @@ -125,7 +125,7 @@ cross_section=cross_section, cross_section_normal=cs_nv, auto_refinement=True, - set_on_trajectory=True, + set_on_trajectory=False, ) ### VISUALIZATION ### diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 0eca124bf..4a46b0242 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -308,9 +308,14 @@ def swept( """ Sweeps a cross-section along a trajectory. The cross-section 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, + vector and is then placed either at the evaluation points of the + trajectory's knots or at the trajectory's control points. This + depends on the value of the set_on_trajectory parameter. + + The sweeping process has some limitations, since the cross-section + cannot be preserved exactly along the whole trajectory. + + The maths behind can be found in the NURBS Book, Piegl & Tiller, 2nd edition, chapter 10.4 Swept Surfaces. Parameters @@ -378,7 +383,7 @@ def transformation_matrices(traj, par_value): vec = [-x[1], x[0], -x[2]] B = [] # avoid dividing by zero - if _np.linalg.norm(_np.cross(x, vec)) > 1e-10: + if _np.linalg.norm(_np.cross(x, vec)) > _settings.TOLERANCE: B.append(_np.cross(x, vec) / _np.linalg.norm(_np.cross(x, vec))) else: vec = [x[2], -x[1], x[0]] @@ -465,7 +470,7 @@ def transformation_matrices(traj, par_value): [_np.sin(angle_of_cs_normal), 0, _np.cos(angle_of_cs_normal)], ] ) - R = _np.where(_np.abs(R) < 1e-10, 0, R) + R = _np.where(_np.abs(R) < _settings.TOLERANCE, 0, R) return A, R From 04b75e1cdcd05d8af5c66c15fec5fb5ad66062e8 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 3 Jul 2024 12:57:40 +0200 Subject: [PATCH 33/53] Fix: get rid of x,y,z --- splinepy/helpme/create.py | 114 +++++++++++++++++++++----------------- 1 file changed, 62 insertions(+), 52 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 4a46b0242..a85a24d79 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -352,12 +352,16 @@ def swept( if not isinstance(trajectory, _Spline): raise NotImplementedError("Sweep only works for splines") if not trajectory.para_dim == 1: - raise NotImplementedError("Trajectory must be 1D") + raise NotImplementedError( + "Trajectory must be of parametric dimension 1" + ) if cross_section_normal is not None and not len(cross_section_normal) == 3: - raise ValueError("Cross section normal must be 3D") + raise ValueError("Cross section normal must be a 3D vector") if not isinstance(auto_refinement, bool): raise ValueError("auto_refinement must be a boolean") + if not isinstance(set_on_trajectory, bool): + raise ValueError("set_on_trajectory must be a boolean") # setting default value for cross_section_normal if cross_section_normal is None: @@ -375,90 +379,96 @@ 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() + # tangent vector 'e1' of trajectory at parameter value 0 + e1 = traj.derivative([par_value[0]], [1]) + e1 = (e1 / _np.linalg.norm(e1)).ravel() - # evaluating a vector normal to x - vec = [-x[1], x[0], -x[2]] + # evaluating a vector normal to e1 + vec = [-e1[1], e1[0], -e1[2]] B = [] # avoid dividing by zero - if _np.linalg.norm(_np.cross(x, vec)) > _settings.TOLERANCE: - B.append(_np.cross(x, vec) / _np.linalg.norm(_np.cross(x, vec))) + if _np.linalg.norm(_np.cross(e1, vec)) > _settings.TOLERANCE: + B.append(_np.cross(e1, vec) / _np.linalg.norm(_np.cross(e1, vec))) else: - vec = [x[2], -x[1], x[0]] - B.append(_np.cross(x, vec) / _np.linalg.norm(_np.cross(x, vec))) + vec = [e1[2], -e1[1], e1[0]] + B.append(_np.cross(e1, vec) / _np.linalg.norm(_np.cross(e1, vec))) - # initialize transformation matrices and x-collection + # initialize transformation matrices and e1-collection T = [] A = [] - x_collection = [] + tang_collection = [] # evaluating transformation matrices for each trajectory point for i in range(len(par_value)): # calculation according to NURBS Book, eq. 10.27 - # 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) onto the plane normal to x - B.append(B[i] - _np.dot(B[i], x) * x) + # tangent vector e1 on trajectory at parameter value i + e1 = traj.derivative([par_value[i]], [1]) + e1 = (e1 / _np.linalg.norm(e1)).ravel() + # collecting tangent vectors for later use + tang_collection.append(e1) + + # projecting B_(i) onto the plane normal to e1 + B.append(B[i] - _np.dot(B[i], e1) * e1) B[i + 1] = B[i + 1] / _np.linalg.norm(B[i + 1]) - # defining y and z axis-vectors - z = B[i + 1] - y = _np.cross(z, x) + # defining e2 and e3 vectors + e3 = B[i + 1] + e2 = _np.cross(e3, e1) # array of transformation matrices from global to local coordinates - T.append(_np.vstack((x, y, z))) + T.append(_np.vstack((e1, e2, e3))) # array of transformation matrices from local to global coordinates A.append(_np.linalg.inv(T[i])) - # own procedure, if trajectory is closed and B[0] != B[-1] + # separate procedure, if trajectory is closed and B[0] != B[-1] + # recalculate B-vector and middle the values between B and B_rec # according to NURBS Book, Piegl & Tiller, 2nd edition, p. 483 - if _np.array_equal( + is_trajectory_closed = _np.array_equal( traj.evaluate([[0]]), traj.evaluate([par_value[-1]]) - ) and not _np.array_equal(B[0], B[-1]): + ) + is_B_start_equal_B_end = _np.array_equal(B[0], B[-1]) + + if is_trajectory_closed and not is_B_start_equal_B_end: # reset transformation matrices T = [] A = [] - B_reverse = [None] * len(B) - B_reverse[0] = B[-1] + # preallocate B_rec + B_rec = [None] * len(B) + # make sure start of B_rec is equal to the end of B + B_rec[0] = B[-1] + # redo the calculation of B using tang_collection from before + # in order to avoid recalculating the tangent vectors; + # calculation according to NURBS Book, eq. 10.27 for i in range(len(par_value)): - # redo the calculation of B using x_collection from before - # according to NURBS Book, eq. 10.27 - 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] + B_rec[i + 1] = ( + B_rec[i] + - _np.dot(B_rec[i], tang_collection[i]) + * tang_collection[i] ) - # 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]) + B_rec[i + 1] = B_rec[i + 1] / _np.linalg.norm(B_rec[i + 1]) + # middle point between B and B_rec + B_rec[i + 1] = (B[i + 1] + B_rec[i + 1]) * 0.5 + # normalizing B_rec + B_rec[i + 1] = B_rec[i + 1] / _np.linalg.norm(B_rec[i + 1]) + # defining e2 and e3 axis-vectors + e3 = B_rec[i + 1] + e2 = _np.cross(e3, tang_collection[i]) # array of transformation matrices from global to local coordinates - T.append(_np.vstack((x_collection[i], y, z))) + T.append(_np.vstack((tang_collection[i], e2, e3))) # 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], rtol=1e-3): + if not _np.allclose(B_rec[0], B_rec[-1], rtol=1e-3): _log.warning( "Vector calculation is not exact due to the " "trajectory being closed in an uncommon way." ) - ### ROTATION MATRIX AROUND Y ### + ### ROTATION MATRIX AROUND e2 VECTOR ### angle_of_cs_normal = _np.arctan2( cross_section_normal[2], cross_section_normal[0] @@ -559,6 +569,7 @@ def transformation_matrices(traj, par_value): swept_spline_cps = [] for i, par_val in enumerate(par_value): temp_csp = [] + # evaluate trajectory if user wants to set cross section on traj. if set_on_trajectory: evals = trajectory.evaluate([par_val]).ravel() # place every control point of cross section separately @@ -567,8 +578,8 @@ def transformation_matrices(traj, par_value): normal_cscp = _np.matmul(R, cscp) # transform cross section to global coordinates normal_cscp = _np.matmul(A[i], normal_cscp) - # check whether user wants to place cross section - # at evaluation points or control points of trajectory + # check if user wants to place cross section at + # evaluation points or control points of trajectory if set_on_trajectory: # translate cross section to trajectory evaluation point normal_cscp += evals @@ -593,7 +604,7 @@ def transformation_matrices(traj, par_value): ), } - # check if spline is rational + # add weights properly if spline is rational if cross_section.is_rational or trajectory.is_rational: def weights(spline): @@ -607,7 +618,6 @@ def weights(spline): trajectory_weights, cross_section_weights ).reshape(-1, 1) spline_type = _NURBS - else: spline_type = _BSpline From fa153f6e8429de14b5ae2ab41e0316baee1b30dc Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 3 Jul 2024 13:29:00 +0200 Subject: [PATCH 34/53] Fix: embedded calc of transf matrices in normal code - without function now --- splinepy/helpme/create.py | 218 ++++++++++++++++++-------------------- 1 file changed, 105 insertions(+), 113 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index a85a24d79..c3131bee1 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -375,115 +375,6 @@ def swept( par_value = trajectory.greville_abscissae() par_value = par_value.reshape(-1, 1) - def transformation_matrices(traj, par_value): - - ### TRANSFORMATION MATRICES ### - - # tangent vector 'e1' of trajectory at parameter value 0 - e1 = traj.derivative([par_value[0]], [1]) - e1 = (e1 / _np.linalg.norm(e1)).ravel() - - # evaluating a vector normal to e1 - vec = [-e1[1], e1[0], -e1[2]] - B = [] - # avoid dividing by zero - if _np.linalg.norm(_np.cross(e1, vec)) > _settings.TOLERANCE: - B.append(_np.cross(e1, vec) / _np.linalg.norm(_np.cross(e1, vec))) - else: - vec = [e1[2], -e1[1], e1[0]] - B.append(_np.cross(e1, vec) / _np.linalg.norm(_np.cross(e1, vec))) - - # initialize transformation matrices and e1-collection - T = [] - A = [] - tang_collection = [] - - # evaluating transformation matrices for each trajectory point - for i in range(len(par_value)): - # calculation according to NURBS Book, eq. 10.27 - # tangent vector e1 on trajectory at parameter value i - e1 = traj.derivative([par_value[i]], [1]) - e1 = (e1 / _np.linalg.norm(e1)).ravel() - # collecting tangent vectors for later use - tang_collection.append(e1) - - # projecting B_(i) onto the plane normal to e1 - B.append(B[i] - _np.dot(B[i], e1) * e1) - B[i + 1] = B[i + 1] / _np.linalg.norm(B[i + 1]) - - # defining e2 and e3 vectors - e3 = B[i + 1] - e2 = _np.cross(e3, e1) - - # array of transformation matrices from global to local coordinates - T.append(_np.vstack((e1, e2, e3))) - - # array of transformation matrices from local to global coordinates - A.append(_np.linalg.inv(T[i])) - - # separate procedure, if trajectory is closed and B[0] != B[-1] - # recalculate B-vector and middle the values between B and B_rec - # according to NURBS Book, Piegl & Tiller, 2nd edition, p. 483 - is_trajectory_closed = _np.array_equal( - traj.evaluate([[0]]), traj.evaluate([par_value[-1]]) - ) - is_B_start_equal_B_end = _np.array_equal(B[0], B[-1]) - - if is_trajectory_closed and not is_B_start_equal_B_end: - # reset transformation matrices - T = [] - A = [] - # preallocate B_rec - B_rec = [None] * len(B) - # make sure start of B_rec is equal to the end of B - B_rec[0] = B[-1] - # redo the calculation of B using tang_collection from before - # in order to avoid recalculating the tangent vectors; - # calculation according to NURBS Book, eq. 10.27 - for i in range(len(par_value)): - B_rec[i + 1] = ( - B_rec[i] - - _np.dot(B_rec[i], tang_collection[i]) - * tang_collection[i] - ) - B_rec[i + 1] = B_rec[i + 1] / _np.linalg.norm(B_rec[i + 1]) - # middle point between B and B_rec - B_rec[i + 1] = (B[i + 1] + B_rec[i + 1]) * 0.5 - # normalizing B_rec - B_rec[i + 1] = B_rec[i + 1] / _np.linalg.norm(B_rec[i + 1]) - # defining e2 and e3 axis-vectors - e3 = B_rec[i + 1] - e2 = _np.cross(e3, tang_collection[i]) - - # array of transformation matrices from global to local coordinates - T.append(_np.vstack((tang_collection[i], e2, e3))) - - # 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_rec[0], B_rec[-1], rtol=1e-3): - _log.warning( - "Vector calculation is not exact due to the " - "trajectory being closed in an uncommon way." - ) - - ### ROTATION MATRIX AROUND e2 VECTOR ### - - angle_of_cs_normal = _np.arctan2( - cross_section_normal[2], cross_section_normal[0] - ) - R = _np.array( - [ - [_np.cos(angle_of_cs_normal), 0, _np.sin(angle_of_cs_normal)], - [0, 1, 0], - [_np.sin(angle_of_cs_normal), 0, _np.cos(angle_of_cs_normal)], - ] - ) - R = _np.where(_np.abs(R) < _settings.TOLERANCE, 0, R) - - return A, R - ### REFINEMENT OF TRAJECTORY ### if auto_refinement: @@ -553,6 +444,110 @@ def transformation_matrices(traj, par_value): par_value = trajectory.greville_abscissae() par_value = par_value.reshape(-1, 1) + ### TRANSFORMATION MATRICES ### + + # tangent vector 'e1' of trajectory at parameter value 0 + e1 = trajectory.derivative([par_value[0]], [1]) + e1 = (e1 / _np.linalg.norm(e1)).ravel() + + # evaluating a vector normal to e1 + vec = [-e1[1], e1[0], -e1[2]] + B = [] + # avoid dividing by zero + if _np.linalg.norm(_np.cross(e1, vec)) > _settings.TOLERANCE: + B.append(_np.cross(e1, vec) / _np.linalg.norm(_np.cross(e1, vec))) + else: + vec = [e1[2], -e1[1], e1[0]] + B.append(_np.cross(e1, vec) / _np.linalg.norm(_np.cross(e1, vec))) + + # initialize transformation matrices and tangent-vector-collection + T = [] + A = [] + tang_collection = [] + + # evaluating transformation matrices for each trajectory point + for i in range(len(par_value)): + # calculation according to NURBS Book, eq. 10.27 + # tangent vector e1 on trajectory at parameter value i + e1 = trajectory.derivative([par_value[i]], [1]) + e1 = (e1 / _np.linalg.norm(e1)).ravel() + # collecting tangent vectors for later use + tang_collection.append(e1) + + # projecting B_(i) onto the plane normal to e1 + B.append(B[i] - _np.dot(B[i], e1) * e1) + B[i + 1] = B[i + 1] / _np.linalg.norm(B[i + 1]) + + # defining e2 and e3 vectors + e3 = B[i + 1] + e2 = _np.cross(e3, e1) + + # array of transformation matrices from global to local coordinates + T.append(_np.vstack((e1, e2, e3))) + + # array of transformation matrices from local to global coordinates + A.append(_np.linalg.inv(T[i])) + + # separate procedure, if trajectory is closed and B[0] != B[-1] + # recalculate B-vector and middle the values between B and B_rec + # according to NURBS Book, Piegl & Tiller, 2nd edition, p. 483 + is_trajectory_closed = _np.array_equal( + trajectory.evaluate([[0]]), trajectory.evaluate([par_value[-1]]) + ) + is_B_start_equal_B_end = _np.array_equal(B[0], B[-1]) + + if is_trajectory_closed and not is_B_start_equal_B_end: + # reset transformation matrices + T = [] + A = [] + # preallocate B_rec + B_rec = [None] * len(B) + # make sure start of B_rec is equal to the end of B + B_rec[0] = B[-1] + # redo the calculation of B using tang_collection from before + # in order to avoid recalculating the tangent vectors; + # calculation according to NURBS Book, eq. 10.27 + for i in range(len(par_value)): + B_rec[i + 1] = ( + B_rec[i] + - _np.dot(B_rec[i], tang_collection[i]) * tang_collection[i] + ) + B_rec[i + 1] = B_rec[i + 1] / _np.linalg.norm(B_rec[i + 1]) + # middle point between B and B_rec + B_rec[i + 1] = (B[i + 1] + B_rec[i + 1]) * 0.5 + # normalizing B_rec + B_rec[i + 1] = B_rec[i + 1] / _np.linalg.norm(B_rec[i + 1]) + # defining e2 and e3 axis-vectors + e3 = B_rec[i + 1] + e2 = _np.cross(e3, tang_collection[i]) + + # array of transformation matrices from global to local coordinates + T.append(_np.vstack((tang_collection[i], e2, e3))) + + # 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_rec[0], B_rec[-1], rtol=1e-3): + _log.warning( + "Vector calculation is not exact due to the " + "trajectory being closed in an uncommon way." + ) + + ### ROTATION MATRIX AROUND e2 VECTOR ### + + angle_of_cs_normal = _np.arctan2( + cross_section_normal[2], cross_section_normal[0] + ) + R = _np.array( + [ + [_np.cos(angle_of_cs_normal), 0, _np.sin(angle_of_cs_normal)], + [0, 1, 0], + [_np.sin(angle_of_cs_normal), 0, _np.cos(angle_of_cs_normal)], + ] + ) + R = _np.where(_np.abs(R) < _settings.TOLERANCE, 0, R) + ### SWEEPING PROCESS ### # evaluate center of cross section and translate to origin @@ -562,14 +557,11 @@ def transformation_matrices(traj, par_value): ).ravel() cross_section.control_points -= cs_center - # evaluate transformation matrices for every trajectory point - A, R = transformation_matrices(trajectory, par_value) - # set cross section control points along trajectory swept_spline_cps = [] for i, par_val in enumerate(par_value): temp_csp = [] - # evaluate trajectory if user wants to set cross section on traj. + # evaluate trajectory if user wants to set cross section on trajectory. if set_on_trajectory: evals = trajectory.evaluate([par_val]).ravel() # place every control point of cross section separately From 058957acfd57ec5d9e302746a6f1b6af2116d494 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 3 Jul 2024 13:39:05 +0200 Subject: [PATCH 35/53] Add: debug log at cs_normal_vector --- examples/show_swept.py | 2 +- splinepy/helpme/create.py | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 25212c14c..5139c7036 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -117,7 +117,7 @@ # user can define the normal vector of the cross section, in case # the cross section is not planar in the x-y plane (default) - cs_nv = np.array([0, 1, 1]) + cs_nv = np.array([0, 0, 1]) ### SWEEP ### swept_surface = splinepy.helpme.create.swept( diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index c3131bee1..e08b6d50b 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -366,6 +366,8 @@ def swept( # setting default value for cross_section_normal if cross_section_normal is None: cross_section_normal = _np.array([0, 0, 1]) + # add debug message + _log.debug("No cross section normal given. Defaulting to [0, 0, 1].") # make copies so we can work on it inplace trajectory = trajectory.create.embedded(3) From c59a26b8f0a7ab9101e8bd274de44595214254be Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 3 Jul 2024 13:43:49 +0200 Subject: [PATCH 36/53] Add: debug log for B division by zero --- splinepy/helpme/create.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index e08b6d50b..526fa069e 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -461,6 +461,10 @@ def swept( else: vec = [e1[2], -e1[1], e1[0]] B.append(_np.cross(e1, vec) / _np.linalg.norm(_np.cross(e1, vec))) + # add debug message + _log.debug( + "Division by zero occurred. Using alternative vector for B." + ) # initialize transformation matrices and tangent-vector-collection T = [] From 2c73c8728e769ffe00d59bd3ca0bdbab1b3f6dce Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Mon, 8 Jul 2024 10:41:55 +0200 Subject: [PATCH 37/53] Add: functionality to rotate cross_section around traj-tangent-vector --- splinepy/helpme/create.py | 37 ++++++++++++++++++++++++++++++------- 1 file changed, 30 insertions(+), 7 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 526fa069e..b3fd7a67d 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -304,6 +304,7 @@ def swept( cross_section_normal=None, auto_refinement=False, set_on_trajectory=False, + rotation_adaption=None, ): """ Sweeps a cross-section along a trajectory. The cross-section @@ -335,6 +336,18 @@ def swept( points of the trajectory's knots. If False, the cross-section will be placed at the control points of the trajectory. Default is False. + rotation_adaption : float + Angle in radians by which the cross-section is rotated around + the trajectory tangent vector. This is an additional rotation + if the user wants to adapt the cross-section rotation. + Example with rectangular crossection: + x + x x x x x x x + x x x x + x x --> rotation around pi/4 --> x x + x x x x + x x x x x x x + x Returns ------- @@ -540,18 +553,28 @@ def swept( "trajectory being closed in an uncommon way." ) - ### ROTATION MATRIX AROUND e2 VECTOR ### + ### ROTATION MATRIX ### + # calculate angle of cross section normal vector around e2-axis angle_of_cs_normal = _np.arctan2( cross_section_normal[2], cross_section_normal[0] ) - R = _np.array( - [ - [_np.cos(angle_of_cs_normal), 0, _np.sin(angle_of_cs_normal)], - [0, 1, 0], - [_np.sin(angle_of_cs_normal), 0, _np.cos(angle_of_cs_normal)], - ] + + # calculate rotation matrix for cross section normal vector + R = _arr.rotation_matrix_around_axis( + axis=[0, 1, 0], rotation=angle_of_cs_normal, degree=False ) + + # rotate cross section around trajectory tangent vector if wanted + if rotation_adaption is not None: + R = _np.matmul( + _arr.rotation_matrix_around_axis( + axis=[1, 0, 0], rotation=rotation_adaption, degree=False + ), + R, + ) + + # remove numerical noise R = _np.where(_np.abs(R) < _settings.TOLERANCE, 0, R) ### SWEEPING PROCESS ### From c2ecb5b15639850f104bafcefbb1d17ce7e396ad Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Mon, 8 Jul 2024 10:42:46 +0200 Subject: [PATCH 38/53] Add: test function to compare swept and extruded --- tests/helpme/test_create.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/helpme/test_create.py b/tests/helpme/test_create.py index 1c79c7b23..9200a2734 100644 --- a/tests/helpme/test_create.py +++ b/tests/helpme/test_create.py @@ -419,3 +419,24 @@ def test_swept_rational_splines(): result = splinepy.helpme.create.swept(cross_section, trajectory) assert result is not None assert result.is_rational + + +def test_swept_to_extruded(): + cross_section = splinepy.NURBS( + degrees=[2], + control_points=[[-0.5, -1 / 3], [0, 2 / 3], [0.5, -1 / 3]], + weights=[1, 0.5, 1], + knot_vectors=[[0, 0, 0, 1, 1, 1]], + ) + trajectory = splinepy.NURBS( + degrees=[1], + control_points=[[0, 0, 0], [0, 0, 3]], + weights=[1, 1], + knot_vectors=[[0, 0, 1, 1]], + ) + + result = splinepy.helpme.create.swept( + cross_section, trajectory, rotation_adaption=np.pi / 2 + ) + extruded_result = splinepy.helpme.create.extruded(cross_section, [0, 0, 3]) + assert np.allclose(result.control_points, extruded_result.control_points) From ad3a125621a3a7a56eddac1064a9c260fd0165ac Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Mon, 8 Jul 2024 11:50:54 +0200 Subject: [PATCH 39/53] Add: test_swept_with_costum_normal now with derivative check --- tests/helpme/test_create.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tests/helpme/test_create.py b/tests/helpme/test_create.py index 9200a2734..bfee65a84 100644 --- a/tests/helpme/test_create.py +++ b/tests/helpme/test_create.py @@ -380,7 +380,11 @@ def test_swept_with_custom_normal(): result = splinepy.helpme.create.swept( cross_section, trajectory, cross_section_normal=custom_normal ) - assert result is not None + trajectory = splinepy.helpme.create.embedded(trajectory, 3) + swept_derivative = result.derivative([[0.5, 0]], [0, 1]) + traj_derivative = trajectory.derivative([[0]], [1]) + + assert np.allclose(swept_derivative, traj_derivative) def test_swept_invalid_inputs(): From 249eb09878efb0a3da2f04065ef315b34230b7b0 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 10 Jul 2024 10:54:49 +0200 Subject: [PATCH 40/53] Rm: auto_refinement; Rm: clean up show_swept file --- examples/show_swept.py | 165 +++++--------------------------------- splinepy/helpme/create.py | 75 ----------------- 2 files changed, 22 insertions(+), 218 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 5139c7036..c3f047e85 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -1,5 +1,3 @@ -import sys - import gustaf as gus import numpy as np @@ -9,126 +7,54 @@ ### TRAJECTORY ### - # arbitrary trajectory + # define a questionmark-trajectory 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.2, 0.4, 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], # 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], ] ), } - # 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": [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.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], - # ] - # ), - # } + # create spline of trajectory dict 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, 3) + # refine trajectory by inserting knots and control points + trajectory.uniform_refine(0, 1) ### CROSS SECTION ### - dict_cross_section = { - "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, 0.0], - [1.0, 2.0, 0.0], - [2.0, 0.0, 0.0], - [3.0, -2.0, 0.0], - [4.0, 0.0, 0.0], - ] - ), - } - - # init cross section as bspline - cross_section = splinepy.BSpline(**dict_cross_section) - # alternatively, use helpme to create a cross section - # cross_section = splinepy.helpme.create.surface_circle(0.5).nurbs + # helpme to create a circular cross section + cross_section = splinepy.helpme.create.surface_circle(0.5).nurbs # user can define the normal vector of the cross section, in case # the cross section is not planar in the x-y plane (default) cs_nv = np.array([0, 0, 1]) ### SWEEP ### + + # create swept surface swept_surface = splinepy.helpme.create.swept( trajectory=trajectory, cross_section=cross_section, cross_section_normal=cs_nv, - auto_refinement=True, set_on_trajectory=False, ) ### VISUALIZATION ### + trajectory.show_options["control_mesh"] = False cross_section.show_options["control_mesh"] = False swept_surface.show_options["control_mesh"] = False @@ -139,50 +65,3 @@ ["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) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index b3fd7a67d..299ffaf98 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -302,7 +302,6 @@ def swept( cross_section, trajectory, cross_section_normal=None, - auto_refinement=False, set_on_trajectory=False, rotation_adaption=None, ): @@ -328,9 +327,6 @@ def swept( cross_section_normal : np.array Normal vector of the cross-section Default is [0, 0, 1] - auto_refinement : bool - If True, the trajectory will be refined at points of - highest curvature. Default is False. set_on_trajectory : bool If True, the cross-section will be placed at the evaluation points of the trajectory's knots. If False, the cross-section @@ -371,8 +367,6 @@ def swept( if cross_section_normal is not None and not len(cross_section_normal) == 3: raise ValueError("Cross section normal must be a 3D vector") - if not isinstance(auto_refinement, bool): - raise ValueError("auto_refinement must be a boolean") if not isinstance(set_on_trajectory, bool): raise ValueError("set_on_trajectory must be a boolean") @@ -390,75 +384,6 @@ def swept( par_value = trajectory.greville_abscissae() par_value = par_value.reshape(-1, 1) - ### REFINEMENT OF TRAJECTORY ### - - if auto_refinement: - ## inserts knots in trajectory area with highest curvature ## - - curv = [] - for i in par_value: - # calculate curvature of trajectory at parametric value i - curv.append( - round(_np.linalg.norm(trajectory.derivative([i], [2])), 2) - ) - # evaluate the par_values-vector indices of the maximum curvature points - max_curv = max(curv) - max_indices = [i for i, x in enumerate(curv) if x == max_curv] - # prepare matrix for the insertion - insertion_values = [] - # compute the new insertion values - par_values = par_value.ravel() - - for maxi in max_indices: - if maxi == 0: - insertion_values.append( - (par_values[maxi] + par_values[maxi + 1]) / 2 - ) - elif maxi == len(par_values) - 1: - insertion_values.append( - (par_values[maxi] + par_values[maxi - 1]) / 2 - ) - else: - insertion_values.append( - (par_values[maxi] + par_values[maxi - 1]) / 2 - ) - insertion_values.append( - (par_values[maxi] + par_values[maxi + 1]) / 2 - ) - - # insert knots into the trajectory's knot vector - insertion_values = _np.unique(insertion_values) - # convert knot vector to list - kv_list = trajectory.knot_vectors[0].numpy() - - # insert knots into the trajectory's knot vector - if any(value in insertion_values for value in kv_list): - add = _np.concatenate((insertion_values, kv_list)) - add = _np.unique(add) - # remove the existing knots - add = add[~_np.isin(add, kv_list)] - # check if there are any knots to insert - if len(add) == 0: - _log.warning("Auto Refinement couldn't insert knots.") - else: - trajectory.insert_knots(0, add) - # give information about the inserted knots - _log.info( - f"Auto Refinement inserted {len(add)} " - "knots into the trajectory." - ) - else: - trajectory.insert_knots(0, insertion_values) - # give information about the inserted knots - _log.info( - f"Auto Refinement inserted {len(insertion_values)} " - "knots into the trajectory." - ) - - # recalculate parameter values - par_value = trajectory.greville_abscissae() - par_value = par_value.reshape(-1, 1) - ### TRANSFORMATION MATRICES ### # tangent vector 'e1' of trajectory at parameter value 0 From dcc4773bd267c025c27428db8d2d25d6a0d92a32 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Fri, 12 Jul 2024 23:27:30 +0200 Subject: [PATCH 41/53] Fix: new test functions --- splinepy/helpme/create.py | 19 ++-- tests/helpme/test_create.py | 196 ++++++++++++++++++++++++++---------- 2 files changed, 155 insertions(+), 60 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 299ffaf98..7907342e7 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -355,20 +355,21 @@ def swept( from splinepy import BSpline as _BSpline from splinepy.spline import Spline as _Spline - # check input type + ### INPUT CHECKS ### + if not isinstance(cross_section, _Spline): - raise NotImplementedError("Sweep only works for splines") + raise TypeError("Sweep only works for splines") if not isinstance(trajectory, _Spline): - raise NotImplementedError("Sweep only works for splines") + raise TypeError("Sweep only works for splines") if not trajectory.para_dim == 1: - raise NotImplementedError( - "Trajectory must be of parametric dimension 1" - ) + raise TypeError("Trajectory must be of parametric dimension 1") + if not isinstance(set_on_trajectory, bool): + raise TypeError("set_on_trajectory must be a boolean") + if not isinstance(rotation_adaption, (float, int, type(None))): + raise TypeError("rotation_adaption must be a float or int") if cross_section_normal is not None and not len(cross_section_normal) == 3: raise ValueError("Cross section normal must be a 3D vector") - if not isinstance(set_on_trajectory, bool): - raise ValueError("set_on_trajectory must be a boolean") # setting default value for cross_section_normal if cross_section_normal is None: @@ -376,6 +377,8 @@ def swept( # add debug message _log.debug("No cross section normal given. Defaulting to [0, 0, 1].") + ### STARTING CALCULATIONS ### + # make copies so we can work on it inplace trajectory = trajectory.create.embedded(3) cross_section = cross_section.create.embedded(3) diff --git a/tests/helpme/test_create.py b/tests/helpme/test_create.py index bfee65a84..c92d23483 100644 --- a/tests/helpme/test_create.py +++ b/tests/helpme/test_create.py @@ -1,5 +1,6 @@ import numpy as np import pytest +from gustaf.utils import arr import splinepy @@ -346,101 +347,192 @@ def test_determinant_spline( def test_swept_basic_functionality(): - cross_section = splinepy.BSpline( - degrees=[2], - control_points=[[0, 0], [0.5, 1], [1, 0]], - knot_vectors=[[0, 0, 0, 1, 1, 1]], + cross_section = splinepy.NURBS( + degrees=[2, 1], + knot_vectors=[ + [0, 0, 0, 1, 1, 1], + [0, 0, 1, 1], + ], + control_points=[ + [-1.0, 0.0], + [-1.0, 1.0], + [0.0, 1.0], + [-2.0, 0.0], + [-2.0, 2.0], + [0.0, 2.0], + ], + weights=[ + [1.0], + [2**-0.5], + [1.0], + [1.0], + [2**-0.5], + [1.0], + ], ) trajectory = splinepy.BSpline( degrees=[3], control_points=[[0, 0], [1, 2], [2, 3], [3, 3]], knot_vectors=[[0, 0, 0, 0, 1, 1, 1, 1]], ) + + # create result --> 3D body should be created result = splinepy.helpme.create.swept(cross_section, trajectory) + + # test result's type assert result is not None + assert result.is_rational + + # test invalid input + invalid_trajectory = "invalid_trajectory" + with pytest.raises(TypeError): + splinepy.helpme.create.swept(cross_section, invalid_trajectory) + + # test result's shape assert ( result.control_points.shape[0] == cross_section.control_points.shape[0] * trajectory.control_points.shape[0] ) + # test if result is 3D + assert result.dim == 3 -def test_swept_with_custom_normal(): - cross_section = splinepy.BSpline( + +def test_swept_to_extruded(): + # create linear swept surface and compare it to extruded surface + + cross_section = splinepy.NURBS( degrees=[2], - control_points=[[0, 0], [0.5, 1], [1, 0]], + control_points=[[-0.5, -1 / 3], [0, 2 / 3], [0.5, -1 / 3]], + weights=[1, 0.5, 1], knot_vectors=[[0, 0, 0, 1, 1, 1]], ) - trajectory = splinepy.BSpline( - degrees=[3], - control_points=[[0, 0], [1, 2], [2, 3], [3, 3]], - knot_vectors=[[0, 0, 0, 0, 1, 1, 1, 1]], + trajectory = splinepy.NURBS( + degrees=[1], + control_points=[[0, 0, 0], [0, 0, 3]], + weights=[1, 1], + knot_vectors=[[0, 0, 1, 1]], ) - custom_normal = np.array([0, 1, 0]) + result = splinepy.helpme.create.swept( - cross_section, trajectory, cross_section_normal=custom_normal + cross_section, trajectory, rotation_adaption=np.pi / 2 ) - trajectory = splinepy.helpme.create.embedded(trajectory, 3) - swept_derivative = result.derivative([[0.5, 0]], [0, 1]) - traj_derivative = trajectory.derivative([[0]], [1]) + extruded_result = splinepy.helpme.create.extruded(cross_section, [0, 0, 3]) + assert np.allclose(result.control_points, extruded_result.control_points) - assert np.allclose(swept_derivative, traj_derivative) +def test_swept_control_point_placing(np_rng): + # check if cross-section's control points are always placed in the correct angle -def test_swept_invalid_inputs(): cross_section = splinepy.BSpline( degrees=[2], control_points=[[0, 0], [0.5, 1], [1, 0]], knot_vectors=[[0, 0, 0, 1, 1, 1]], ) - invalid_trajectory = "invalid_trajectory" - with pytest.raises(NotImplementedError): - splinepy.helpme.create.swept(cross_section, invalid_trajectory) - invalid_cross_section = "invalid_cross_section" trajectory = splinepy.BSpline( degrees=[3], - control_points=[[0, 0], [1, 2], [2, 3], [3, 3]], + control_points=[[0, 0, 0], [1, 0, 1], [0, 0, 2], [0, 0, 3]], knot_vectors=[[0, 0, 0, 0, 1, 1, 1, 1]], ) - with pytest.raises(NotImplementedError): - splinepy.helpme.create.swept(invalid_cross_section, trajectory) + result = splinepy.helpme.create.swept(cross_section, trajectory) + def calculate_cs_normal_vector(cs_cp, res_cp, rand): -def test_swept_rational_splines(): - cross_section = splinepy.NURBS( + # take the control points of the chosen cross-section + taken_cs_cp = res_cp[ + rand * len(cs_cp) : rand * len(cs_cp) + len(cs_cp) + ] + P1 = taken_cs_cp[0] + P2 = taken_cs_cp[1] + P3 = taken_cs_cp[2] + + # calculate normal vector of cross-section + v1 = P2 - P1 + v2 = P3 - P1 + return np.cross(v1, v2) + + cs_cp = cross_section.control_points + res_cp = result.control_points + + # choose id of cross-section to take from result + rand = np_rng.integers(0, len(res_cp) / len(cs_cp), 1)[0] + + normal_vector = calculate_cs_normal_vector(cs_cp, res_cp, rand) + + # evaluate trajectory parameter value at the position of the cross-section + rand_traj_pos = trajectory.greville_abscissae()[rand] + # calculate tangent vector of trajectory at the position of the cross-section + rand_traj_tangent = trajectory.derivative([rand_traj_pos], 1).ravel() + + # check if normal vector is parallel to tangent vector + # --> then transformation of cross-section is correct + assert np.allclose(np.cross(rand_traj_tangent, normal_vector).all(), 0) + + ## TEST ALSO WITH CUSTOM NORMAL VECTOR ## + + # make sure that the angle between trajectory's tangent and cross-section's + # normal vector is always the same + custom_normal = np.array([1, 0, 1]) + + angle_of_custom_cs_normal = np.arctan2(custom_normal[2], custom_normal[0]) + # calculate rotation matrix for cross-section normal vector + R = arr.rotation_matrix_around_axis( + axis=[0, 1, 0], rotation=angle_of_custom_cs_normal, degree=False + ) + # in sweeping, the cross-section is rotated, we mimic this here + test_normal = np.matmul(R, custom_normal) + + start_traj_tang = trajectory.derivative([[0]], 1).ravel() + angle_at_start = np.arccos( + np.dot(start_traj_tang, test_normal) + / (np.linalg.norm(start_traj_tang) * np.linalg.norm(test_normal)) + ) + + result_with_cus_normal = splinepy.helpme.create.swept( + cross_section, trajectory, cross_section_normal=custom_normal + ) + + result_with_cus_normal_cps = result_with_cus_normal.control_points + normal_vector_custom_version = ( + calculate_cs_normal_vector(cs_cp, result_with_cus_normal_cps, rand) + * -1 + ) + + angle_at_rand = np.arccos( + np.dot(rand_traj_tangent, normal_vector_custom_version) + / ( + np.linalg.norm(rand_traj_tangent) + * np.linalg.norm(normal_vector_custom_version) + ) + ) + + assert np.allclose(angle_at_start, angle_at_rand) + + +def test_swept_cs_centering(np_rng): + # check if cross-section's center lays on the trajectory + + cross_section = splinepy.BSpline( degrees=[2], control_points=[[0, 0], [0.5, 1], [1, 0]], - weights=[1, 0.5, 1], knot_vectors=[[0, 0, 0, 1, 1, 1]], ) - trajectory = splinepy.NURBS( + trajectory = splinepy.BSpline( degrees=[3], - control_points=[[0, 0], [1, 2], [2, 3], [3, 3]], - weights=[1, 0.5, 0.5, 1], + control_points=[[0, 0, 0], [1, 2, 0], [2, 3, 0], [3, 3, 0]], knot_vectors=[[0, 0, 0, 0, 1, 1, 1, 1]], ) result = splinepy.helpme.create.swept(cross_section, trajectory) - assert result is not None - assert result.is_rational - -def test_swept_to_extruded(): - cross_section = splinepy.NURBS( - degrees=[2], - control_points=[[-0.5, -1 / 3], [0, 2 / 3], [0.5, -1 / 3]], - weights=[1, 0.5, 1], - knot_vectors=[[0, 0, 0, 1, 1, 1]], - ) - trajectory = splinepy.NURBS( - degrees=[1], - control_points=[[0, 0, 0], [0, 0, 3]], - weights=[1, 1], - knot_vectors=[[0, 0, 1, 1]], - ) + # choose random parameter value of trajectory + r_par_val = np_rng.random(1) + # evaluate swept spline at the chosen position of the trajectory + # --> cross-section parameter value 0.5 ensures that the center of the + # cross-section is evaluated + coords = result.evaluate([[0.5, r_par_val[0]]]).ravel() + # evaluate trajectory at the chosen position of the trajectory + ref_coords = trajectory.evaluate([r_par_val]).ravel() - result = splinepy.helpme.create.swept( - cross_section, trajectory, rotation_adaption=np.pi / 2 - ) - extruded_result = splinepy.helpme.create.extruded(cross_section, [0, 0, 3]) - assert np.allclose(result.control_points, extruded_result.control_points) + assert np.allclose(coords, ref_coords) From 27c60943468e46d7d443b96127d57eb6996624f0 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Fri, 12 Jul 2024 23:28:32 +0200 Subject: [PATCH 42/53] Rm: german comment --- examples/show_swept.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index c3f047e85..2cafeb431 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -15,7 +15,7 @@ ], "control_points": np.array( [ - [0.5, 0], # Startpunkt + [0.5, 0], [0.5, 2], [1.0, 3], [2.0, 4], From bc627a8526e64a7f39293c911d89bd956ca55b0e Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Fri, 12 Jul 2024 23:32:08 +0200 Subject: [PATCH 43/53] Add: capabilities of sweeping in the docstring --- examples/show_swept.py | 2 +- splinepy/helpme/create.py | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 2cafeb431..00979a049 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -15,7 +15,7 @@ ], "control_points": np.array( [ - [0.5, 0], + [0.5, 0], [0.5, 2], [1.0, 3], [2.0, 4], diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 7907342e7..dd1885b68 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -310,7 +310,9 @@ def swept( receives rotation into the direction of the trajectory tangent vector and is then placed either at the evaluation points of the trajectory's knots or at the trajectory's control points. This - depends on the value of the set_on_trajectory parameter. + depends on the value of the set_on_trajectory parameter. It can + create both a surface or a solid, depending on the dimension of + the cross-section. The sweeping process has some limitations, since the cross-section cannot be preserved exactly along the whole trajectory. From 08d903d937859e3240860924c9585ba632f437af Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Fri, 12 Jul 2024 23:57:17 +0200 Subject: [PATCH 44/53] Fix: changed indentation in docstring --- splinepy/helpme/create.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index dd1885b68..e0025eeb4 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -335,17 +335,17 @@ def swept( will be placed at the control points of the trajectory. Default is False. rotation_adaption : float - Angle in radians by which the cross-section is rotated around - the trajectory tangent vector. This is an additional rotation - if the user wants to adapt the cross-section rotation. - Example with rectangular crossection: - x - x x x x x x x - x x x x - x x --> rotation around pi/4 --> x x - x x x x - x x x x x x x - x + Angle in radians by which the cross-section is rotated around + the trajectory tangent vector. This is an additional rotation + if the user wants to adapt the cross-section rotation. + Example with rectangular crossection: + x + x x x x x x x + x x x x + x x --> rotation around pi/4 --> x x + x x x x + x x x x x x x + x Returns ------- From 86e00f40b8c0348a19f46c74eab45582cdc57927 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Sat, 13 Jul 2024 00:11:23 +0200 Subject: [PATCH 45/53] Fix: added points in docstring --- splinepy/helpme/create.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index e0025eeb4..18676ba26 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -339,13 +339,13 @@ def swept( the trajectory tangent vector. This is an additional rotation if the user wants to adapt the cross-section rotation. Example with rectangular crossection: - x - x x x x x x x - x x x x - x x --> rotation around pi/4 --> x x - x x x x - x x x x x x x - x + . x + . x x x x x x x + . x x x x + . x x --> rotation around pi/4 --> x x + . x x x x + . x x x x x x x + . x Returns ------- From 016b7c6636f7ea51826a8053ce2181c984788f12 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Mon, 9 Sep 2024 16:25:51 +0200 Subject: [PATCH 46/53] Add: show_swept-file now shows more examples --- examples/show_swept.py | 99 +++++++++++++++++++++++++++++++++------ splinepy/helpme/create.py | 2 +- 2 files changed, 86 insertions(+), 15 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 00979a049..eb509054c 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -7,7 +7,7 @@ ### TRAJECTORY ### - # define a questionmark-trajectory + # define a hook-trajectory dict_trajectory = { "degrees": [3], "knot_vectors": [ @@ -34,34 +34,105 @@ # refine trajectory by inserting knots and control points trajectory.uniform_refine(0, 1) - ### CROSS SECTION ### + ### CROSS SECTIONS ### - # helpme to create a circular cross section - cross_section = splinepy.helpme.create.surface_circle(0.5).nurbs + # 1. create a circular line-cross-section + cross_section_circle1D = splinepy.helpme.create.circle(0.5).nurbs + + # 2. create a circular surface-cross-section + cross_section_circle2D = splinepy.helpme.create.surface_circle(0.5).nurbs + + # 3. create a rectangular surface-cross-section + cross_section_rect2D = splinepy.helpme.create.box(1, 1).nurbs # user can define the normal vector of the cross section, in case - # the cross section is not planar in the x-y plane (default) - cs_nv = np.array([0, 0, 1]) + # the cross section is not planar in the x-y plane (or the user + # wants crooked sweeping) + cs_nv = np.array([1, 0, 1]) ### SWEEP ### # create swept surface - swept_surface = splinepy.helpme.create.swept( + swept_surface_1D_circle = splinepy.helpme.create.swept( + trajectory=trajectory, + cross_section=cross_section_circle1D, + cross_section_normal=None, + set_on_trajectory=False, + rotation_adaption=None, + ) + + # create swept solid (circular nr. 1) + # the cross-sections are set on the trajectory's control points (default) + swept_surface_2D_circle_1 = splinepy.helpme.create.swept( trajectory=trajectory, - cross_section=cross_section, + cross_section=cross_section_circle2D, + cross_section_normal=None, + set_on_trajectory=False, + rotation_adaption=None, + ) + + # create crooked swept solid (circular nr. 2) + # the cross-sections are set on the trajectory's evaluation points + swept_surface_2D_circle_2 = splinepy.helpme.create.swept( + trajectory=trajectory, + cross_section=cross_section_circle2D, + cross_section_normal=None, + set_on_trajectory=True, + rotation_adaption=None, + ) + + # create swept solid set on trajectory (circular nr. 3) + # the cross-section's normal vector is not default; crooked sweeping + swept_surface_2D_circle_3 = splinepy.helpme.create.swept( + trajectory=trajectory, + cross_section=cross_section_circle2D, cross_section_normal=cs_nv, set_on_trajectory=False, + rotation_adaption=None, ) - ### VISUALIZATION ### + # create swept solid (rectangular nr. 1) + swept_surface_2D_rect_1 = splinepy.helpme.create.swept( + trajectory=trajectory, + cross_section=cross_section_rect2D, + cross_section_normal=None, + set_on_trajectory=False, + rotation_adaption=None, + ) + + # create swept solid (rectangular nr. 2) + # rotation adaption with 45 degrees + swept_surface_2D_rect_2 = splinepy.helpme.create.swept( + trajectory=trajectory, + cross_section=cross_section_rect2D, + cross_section_normal=None, + set_on_trajectory=False, + rotation_adaption=45 * np.pi / 180, + ) - trajectory.show_options["control_mesh"] = False - cross_section.show_options["control_mesh"] = False - swept_surface.show_options["control_mesh"] = False + ### VISUALIZATION ### gus.show( ["Trajectory", trajectory], - ["Cross Section", cross_section], - ["Swept Surface", swept_surface], + ["1D Cross Section", cross_section_circle1D], + ["Swept Surface", swept_surface_1D_circle], + resolution=101, + control_mesh=False, + ) + + gus.show( + ["2D Cross Section", cross_section_circle2D], + ["Swept Solid - Set on Control Points", swept_surface_2D_circle_1], + ["Swept Solid - Set on Evaluation Points", swept_surface_2D_circle_2], + ["Swept Solid - Crooked Sweeping", swept_surface_2D_circle_3], + resolution=101, + control_mesh=False, + ) + + gus.show( + ["New Cross Section", cross_section_rect2D], + ["Swept Solid without Rotation", swept_surface_2D_rect_1, trajectory], + ["Swept Solid with 45° Rotation", swept_surface_2D_rect_2, trajectory], resolution=101, + control_mesh=False, ) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 18676ba26..9fb9a3314 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -338,7 +338,7 @@ def swept( Angle in radians by which the cross-section is rotated around the trajectory tangent vector. This is an additional rotation if the user wants to adapt the cross-section rotation. - Example with rectangular crossection: + Example with rectangular cross-section: . x . x x x x x x x . x x x x From 43e48d4cc73755e2f2d7e37a8391b95426e210dc Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Tue, 10 Sep 2024 18:54:40 +0200 Subject: [PATCH 47/53] Rm: .all()-function in test_create-file --- tests/helpme/test_create.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/helpme/test_create.py b/tests/helpme/test_create.py index c92d23483..38f09b74a 100644 --- a/tests/helpme/test_create.py +++ b/tests/helpme/test_create.py @@ -468,7 +468,7 @@ def calculate_cs_normal_vector(cs_cp, res_cp, rand): # check if normal vector is parallel to tangent vector # --> then transformation of cross-section is correct - assert np.allclose(np.cross(rand_traj_tangent, normal_vector).all(), 0) + assert np.allclose(np.cross(rand_traj_tangent, normal_vector), 0) ## TEST ALSO WITH CUSTOM NORMAL VECTOR ## From 255b8743a76433c4521d6d3e07e344fd92831478 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Wed, 11 Sep 2024 14:34:54 +0200 Subject: [PATCH 48/53] Fix: show_options and variable names in show_swept --- examples/show_swept.py | 99 +++++++++++++++++++++++++----------------- 1 file changed, 59 insertions(+), 40 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index eb509054c..983ac6dbb 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -36,14 +36,14 @@ ### CROSS SECTIONS ### - # 1. create a circular line-cross-section - cross_section_circle1D = splinepy.helpme.create.circle(0.5).nurbs + # 1. create a circular 1D-line-cross-section + cross_section_circle = splinepy.helpme.create.circle(0.5).nurbs - # 2. create a circular surface-cross-section - cross_section_circle2D = splinepy.helpme.create.surface_circle(0.5).nurbs + # 2. create a circular 2D-surface-cross-section + cross_section_disk = splinepy.helpme.create.surface_circle(0.5).nurbs - # 3. create a rectangular surface-cross-section - cross_section_rect2D = splinepy.helpme.create.box(1, 1).nurbs + # 3. create a rectangular 2D-surface-cross-section + cross_section_plate = splinepy.helpme.create.box(1, 1).nurbs # user can define the normal vector of the cross section, in case # the cross section is not planar in the x-y plane (or the user @@ -53,48 +53,48 @@ ### SWEEP ### # create swept surface - swept_surface_1D_circle = splinepy.helpme.create.swept( + swept_surface_circle = splinepy.helpme.create.swept( trajectory=trajectory, - cross_section=cross_section_circle1D, + cross_section=cross_section_circle, cross_section_normal=None, set_on_trajectory=False, rotation_adaption=None, ) - # create swept solid (circular nr. 1) - # the cross-sections are set on the trajectory's control points (default) - swept_surface_2D_circle_1 = splinepy.helpme.create.swept( + # create crooked swept solid (circular nr. 1) + # the cross-section's normal vector is not default; crooked sweeping + swept_solid_disk_1 = splinepy.helpme.create.swept( trajectory=trajectory, - cross_section=cross_section_circle2D, - cross_section_normal=None, + cross_section=cross_section_disk, + cross_section_normal=cs_nv, set_on_trajectory=False, rotation_adaption=None, ) - # create crooked swept solid (circular nr. 2) - # the cross-sections are set on the trajectory's evaluation points - swept_surface_2D_circle_2 = splinepy.helpme.create.swept( + # create swept solid (circular nr. 2) + # the cross-sections are set on the trajectory's control points (default) + swept_solid_disk_2 = splinepy.helpme.create.swept( trajectory=trajectory, - cross_section=cross_section_circle2D, + cross_section=cross_section_disk, cross_section_normal=None, - set_on_trajectory=True, + set_on_trajectory=False, rotation_adaption=None, ) - # create swept solid set on trajectory (circular nr. 3) - # the cross-section's normal vector is not default; crooked sweeping - swept_surface_2D_circle_3 = splinepy.helpme.create.swept( + # create swept solid (circular nr. 3) + # the cross-sections are set on the trajectory's evaluation points + swept_solid_disk_3 = splinepy.helpme.create.swept( trajectory=trajectory, - cross_section=cross_section_circle2D, - cross_section_normal=cs_nv, - set_on_trajectory=False, + cross_section=cross_section_disk, + cross_section_normal=None, + set_on_trajectory=True, rotation_adaption=None, ) # create swept solid (rectangular nr. 1) - swept_surface_2D_rect_1 = splinepy.helpme.create.swept( + swept_solid_plate_1 = splinepy.helpme.create.swept( trajectory=trajectory, - cross_section=cross_section_rect2D, + cross_section=cross_section_plate, cross_section_normal=None, set_on_trajectory=False, rotation_adaption=None, @@ -102,9 +102,9 @@ # create swept solid (rectangular nr. 2) # rotation adaption with 45 degrees - swept_surface_2D_rect_2 = splinepy.helpme.create.swept( + swept_solid_plate_2 = splinepy.helpme.create.swept( trajectory=trajectory, - cross_section=cross_section_rect2D, + cross_section=cross_section_plate, cross_section_normal=None, set_on_trajectory=False, rotation_adaption=45 * np.pi / 180, @@ -112,27 +112,46 @@ ### VISUALIZATION ### + # first window: swept surface gus.show( ["Trajectory", trajectory], - ["1D Cross Section", cross_section_circle1D], - ["Swept Surface", swept_surface_1D_circle], - resolution=101, + ["1D Cross Section", cross_section_circle], + ["Swept Surface", swept_surface_circle], + resolution=51, control_mesh=False, + control_point_ids=False, ) + # adjust show options + swept_solid_disk_2.show_options["alpha"] = 0.3 + swept_solid_disk_3.show_options["alpha"] = 0.3 + trajectory.show_options["control_points"] = False + + # second window: swept solids (circular) gus.show( - ["2D Cross Section", cross_section_circle2D], - ["Swept Solid - Set on Control Points", swept_surface_2D_circle_1], - ["Swept Solid - Set on Evaluation Points", swept_surface_2D_circle_2], - ["Swept Solid - Crooked Sweeping", swept_surface_2D_circle_3], - resolution=101, + ["2D Cross Section", cross_section_disk], + ["Swept Solid - Crooked Sweeping", swept_solid_disk_1], + [ + "Swept Solid - Set on Control Points", + swept_solid_disk_2, + trajectory, + ], + [ + "Swept Solid - Set on Evaluation Points", + swept_solid_disk_3, + trajectory, + ], + resolution=51, control_mesh=False, + control_point_ids=False, ) + # third window: swept solids (rectangular) gus.show( - ["New Cross Section", cross_section_rect2D], - ["Swept Solid without Rotation", swept_surface_2D_rect_1, trajectory], - ["Swept Solid with 45° Rotation", swept_surface_2D_rect_2, trajectory], - resolution=101, + ["New Cross Section", cross_section_plate], + ["Swept Solid without Rotation", swept_solid_plate_1], + ["Swept Solid with 45° Rotation", swept_solid_plate_2], + resolution=51, control_mesh=False, + control_point_ids=False, ) From 4207da0730995e0f700bc4610cf815ea7b12f846 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Mon, 16 Sep 2024 10:55:02 +0200 Subject: [PATCH 49/53] Fix: changed comments of cs_nv --- examples/show_swept.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/show_swept.py b/examples/show_swept.py index 983ac6dbb..ecc16ecc6 100644 --- a/examples/show_swept.py +++ b/examples/show_swept.py @@ -45,9 +45,9 @@ # 3. create a rectangular 2D-surface-cross-section cross_section_plate = splinepy.helpme.create.box(1, 1).nurbs - # user can define the normal vector of the cross section, in case - # the cross section is not planar in the x-y plane (or the user - # wants crooked sweeping) + # Define a custom normal vector for the cross-section when: + # a) The cross-section is not planar in the x-y plane, or + # b) Crooked sweeping is desired cs_nv = np.array([1, 0, 1]) ### SWEEP ### From 57517044d9f983a0b1f5d923035b2cfcf79f64e9 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Mon, 16 Sep 2024 11:27:58 +0200 Subject: [PATCH 50/53] Add: assertions and comments in test_basic_functionality-function in test_create-file --- tests/helpme/test_create.py | 31 ++++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/tests/helpme/test_create.py b/tests/helpme/test_create.py index 38f09b74a..93167f1ec 100644 --- a/tests/helpme/test_create.py +++ b/tests/helpme/test_create.py @@ -346,6 +346,10 @@ def test_determinant_spline( ), f"{sp_i.whatami} at index {idx} failed determinant spline" +### TESTS FOR SWEEPING ### + + +# test the basic functionality of the swept-function def test_swept_basic_functionality(): cross_section = splinepy.NURBS( degrees=[2, 1], @@ -398,9 +402,30 @@ def test_swept_basic_functionality(): # test if result is 3D assert result.dim == 3 + # test if result's degrees are correct + assert np.allclose( + result.degrees, + [ + cross_section.degrees[0], + cross_section.degrees[1], + trajectory.degrees[0], + ], + ) + + # test if result's knot vectors are correct + assert np.allclose( + result.knot_vectors[0], cross_section.knot_vectors[0] + ), "Knot vector of first cross-section dimension does not match." + assert np.allclose( + result.knot_vectors[1], cross_section.knot_vectors[1] + ), "Knot vector of second cross-section dimension does not match." + assert np.allclose( + result.knot_vectors[2], trajectory.knot_vectors[0] + ), "Knot vector of trajectory does not match." + +# create linear swept surface and compare it to extruded surface def test_swept_to_extruded(): - # create linear swept surface and compare it to extruded surface cross_section = splinepy.NURBS( degrees=[2], @@ -422,8 +447,8 @@ def test_swept_to_extruded(): assert np.allclose(result.control_points, extruded_result.control_points) +# check if cross-section's control points are always placed in the correct angle def test_swept_control_point_placing(np_rng): - # check if cross-section's control points are always placed in the correct angle cross_section = splinepy.BSpline( degrees=[2], @@ -511,8 +536,8 @@ def calculate_cs_normal_vector(cs_cp, res_cp, rand): assert np.allclose(angle_at_start, angle_at_rand) +# check if cross-section's center lays on the trajectory def test_swept_cs_centering(np_rng): - # check if cross-section's center lays on the trajectory cross_section = splinepy.BSpline( degrees=[2], From 53851c7a5923b0db170d2e57135c7d3e9ea28579 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Mon, 16 Sep 2024 12:07:24 +0200 Subject: [PATCH 51/53] Fix: error messages in create-file --- splinepy/helpme/create.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 9fb9a3314..7c197b176 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -371,13 +371,13 @@ def swept( raise TypeError("rotation_adaption must be a float or int") if cross_section_normal is not None and not len(cross_section_normal) == 3: - raise ValueError("Cross section normal must be a 3D vector") + raise ValueError("cross_section_normal must be a 3D vector") # setting default value for cross_section_normal if cross_section_normal is None: cross_section_normal = _np.array([0, 0, 1]) # add debug message - _log.debug("No cross section normal given. Defaulting to [0, 0, 1].") + _log.debug("No cross_section_normal given. Defaulting to [0, 0, 1].") ### STARTING CALCULATIONS ### From ba5264df9c0bf38c4b4643f467a6528f09031d7f Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Mon, 16 Sep 2024 12:49:53 +0200 Subject: [PATCH 52/53] Fix: clarified error messages at input checks of create-file --- splinepy/helpme/create.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index 7c197b176..ae5eff5a6 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -355,16 +355,17 @@ def swept( from splinepy import NURBS as _NURBS from splinepy import BSpline as _BSpline - from splinepy.spline import Spline as _Spline ### INPUT CHECKS ### - if not isinstance(cross_section, _Spline): - raise TypeError("Sweep only works for splines") - if not isinstance(trajectory, _Spline): - raise TypeError("Sweep only works for splines") + if not isinstance(cross_section, (_BSpline, _NURBS)): + raise TypeError( + "cross_section must be an instance of BSpline or NURBS" + ) + if not isinstance(trajectory, (_BSpline, _NURBS)): + raise TypeError("trajectory must be an instance of BSpline or NURBS") if not trajectory.para_dim == 1: - raise TypeError("Trajectory must be of parametric dimension 1") + raise TypeError("trajectory must be of parametric dimension 1") if not isinstance(set_on_trajectory, bool): raise TypeError("set_on_trajectory must be a boolean") if not isinstance(rotation_adaption, (float, int, type(None))): From b2a3d5e2071fa9eb93f4b4f579e36b99e0bfeb35 Mon Sep 17 00:00:00 2001 From: Guenther Obermair Date: Tue, 17 Sep 2024 18:15:11 +0200 Subject: [PATCH 53/53] Fix: rotation matrix calculation now also for e2 entries --- splinepy/helpme/create.py | 52 +++++++++++++++++++++++++-------------- 1 file changed, 33 insertions(+), 19 deletions(-) diff --git a/splinepy/helpme/create.py b/splinepy/helpme/create.py index ae5eff5a6..3a74ad5bf 100644 --- a/splinepy/helpme/create.py +++ b/splinepy/helpme/create.py @@ -373,7 +373,6 @@ def swept( if cross_section_normal is not None and not len(cross_section_normal) == 3: raise ValueError("cross_section_normal must be a 3D vector") - # setting default value for cross_section_normal if cross_section_normal is None: cross_section_normal = _np.array([0, 0, 1]) @@ -485,18 +484,33 @@ def swept( ) ### ROTATION MATRIX ### + # rotates cross-section normal vector into global e1 direction - # calculate angle of cross section normal vector around e2-axis - angle_of_cs_normal = _np.arctan2( - cross_section_normal[2], cross_section_normal[0] - ) + # skip calculation if cross-section normal vector is default + if _np.array_equal(cross_section_normal, _np.array([0, 0, 1])): + R = _np.array([[0.0, 0.0, 1.0], [0.0, 1.0, 0.0], [-1.0, 0.0, 0.0]]) + # else, calculate rotation matrix for given cross-section normal vector + else: + # calculate angle of cross-section normal vector around e2-axis + angle_of_cs_normal_1 = _np.arctan2( + cross_section_normal[2], cross_section_normal[0] + ) - # calculate rotation matrix for cross section normal vector - R = _arr.rotation_matrix_around_axis( - axis=[0, 1, 0], rotation=angle_of_cs_normal, degree=False - ) + # calculate angle of cross-section normal vector around e3-axis + angle_of_cs_normal_2 = _np.arctan2( + cross_section_normal[1], cross_section_normal[2] + ) + + # calculate rotation matrix for cross-section normal vector + R1 = _arr.rotation_matrix_around_axis( + axis=[0, 1, 0], rotation=angle_of_cs_normal_1, degree=False + ) + R2 = _arr.rotation_matrix_around_axis( + axis=[0, 0, 1], rotation=angle_of_cs_normal_2, degree=False + ) + R = _np.matmul(R2, R1) - # rotate cross section around trajectory tangent vector if wanted + # rotate cross-section around trajectory tangent vector (e1) if wanted if rotation_adaption is not None: R = _np.matmul( _arr.rotation_matrix_around_axis( @@ -510,33 +524,33 @@ def swept( ### SWEEPING PROCESS ### - # evaluate center of cross section and translate to origin + # evaluate center of cross-section and translate to origin cross_para_center = _np.mean(cross_section.parametric_bounds, axis=0) cs_center = cross_section.evaluate( cross_para_center.reshape(-1, cross_section.para_dim) ).ravel() cross_section.control_points -= cs_center - # set cross section control points along trajectory + # set cross-section control points along trajectory swept_spline_cps = [] for i, par_val in enumerate(par_value): temp_csp = [] - # evaluate trajectory if user wants to set cross section on trajectory. + # evaluate trajectory if user wants to set cross-section on trajectory. if set_on_trajectory: evals = trajectory.evaluate([par_val]).ravel() - # place every control point of cross section separately + # place every control point of cross-section separately for cscp in cross_section.control_points: - # rotate cross section in trajectory direction + # rotate cross-section in trajectory direction normal_cscp = _np.matmul(R, cscp) - # transform cross section to global coordinates + # transform cross-section to global coordinates normal_cscp = _np.matmul(A[i], normal_cscp) - # check if user wants to place cross section at + # check if user wants to place cross-section at # evaluation points or control points of trajectory if set_on_trajectory: - # translate cross section to trajectory evaluation point + # translate cross-section to trajectory evaluation point normal_cscp += evals else: - # translate cross section to trajectory control point + # translate cross-section to trajectory control point normal_cscp += trajectory.control_points[i] # append control point to list temp_csp.append(normal_cscp)