Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# Institute: Norwegian University of Science and Technology (NTNU)
# Date: March 2016
#
from sys import path
path.append('../')
from splipy import *
import splipy.curve_factory as curves
import splipy.surface_factory as surfaces
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from math import pi, cos, sin
# create the three sides of the triangle, each consisting of a circle segment
c1 = curves.circle_segment(pi/3)
c2 = c1.clone().rotate(2*pi/3) + [1,0]
c3 = c1.clone().rotate(4*pi/3) + [cos(pi/3), sin(pi/3)]
# merge the three cirlces into one, and center it at the origin
c = c1.append(c2).append(c3)
c -= c.center()
# plot the reuleaux triangle
t = np.linspace(c.start(0), c.end(0), 151) # 151 parametric evaluation points
x = c(t) # evaluate (x,y)-coordinates
plt.plot(x[:,0], x[:,1])
plt.axis('equal')
plt.show()
# split the triangle in two, and align this with the y-axis
two_parts = c.split((c.start(0) + c.end(0)) / 2.0)
cprime = +(np.sqrt(abs(rx**2*ry**2 - rx**2*xp[1]**2 - ry**2*xp[0]**2) /
(rx**2*xp[1]**2 + ry**2*xp[0]**2)) *
np.array([rx*xp[1]/ry, -ry*xp[0]/rx]))
center = np.linalg.solve(R.T, cprime) + (startpoint+xend)/2
def arccos(vec1, vec2):
return (np.sign(vec1[0]*vec2[1] - vec1[1]*vec2[0]) *
np.arccos(vec1.dot(vec2)/np.linalg.norm(vec1)/np.linalg.norm(vec2)))
tmp1 = np.divide( xp - cprime, [rx,ry])
tmp2 = np.divide(-xp - cprime, [rx,ry])
theta1 = arccos(np.array([1,0]), tmp1)
delta_t= arccos(tmp1, tmp2) % (2*np.pi)
if not sweep_flag and delta_t > 0:
delta_t -= 2*np.pi
elif sweep_flag and delta_t < 0:
delta_t += 2*np.pi
curve_piece = (curve_factory.circle_segment(delta_t)*[rx,ry]).rotate(theta1) + center
# curve_piece = curve_factory.circle_segment(delta_t)
elif piece[0] == 'z' or piece[0] == 'Z':
# periodic curve
# curve_piece = Curve(BSplineBasis(2), [startpoint, last_curve[0]])
# curve_piece.reparam([0, curve_piece.length()])
# last_curve.append(curve_piece).make_periodic(0)
last_curve.make_periodic(0)
result.append(last_curve)
last_curve = None
continue
else:
raise RuntimeError('Unknown path parameter:' + piece)
if(curve_piece.length()>state.controlpoint_absolute_tolerance):
curve_piece.reparam([0, curve_piece.length()])
:param float theta: Angle to revolve, in radians
:param array-like axis: Axis of rotation
:return: The revolved surface
:rtype: Volume
"""
surf = surf.clone() # clone input surface, throw away old reference
surf.set_dimension(3) # add z-components (if not already present)
surf.force_rational() # add weight (if not already present)
# align axis with the z-axis
normal_theta = atan2(axis[1], axis[0])
normal_phi = atan2(sqrt(axis[0]**2 + axis[1]**2), axis[2])
surf.rotate(-normal_theta, [0,0,1])
surf.rotate(-normal_phi, [0,1,0])
path = CurveFactory.circle_segment(theta=theta)
n = len(surf) # number of control points of the surface
m = len(path) # number of control points of the sweep
cp = np.zeros((m * n, 4))
dt = np.sign(theta)*(path.knots(0)[1] - path.knots(0)[0]) / 2.0
for i in range(m):
weight = path[i,-1]
cp[i * n:(i + 1) * n, :] = np.reshape(surf.controlpoints.transpose(1, 0, 2), (n, 4))
cp[i * n:(i + 1) * n, 2] *= weight
cp[i * n:(i + 1) * n, 3] *= weight
surf.rotate(dt)
result = Volume(surf.bases[0], surf.bases[1], path.bases[0], cp, True)
# rotate it back again
result.rotate(normal_phi, [0,1,0])
result.rotate(normal_theta, [0,0,1])
:param float theta: Angle to revolve, in radians
:param array-like axis: Axis of rotation
:return: The revolved surface
:rtype: Surface
"""
curve = curve.clone() # clone input curve, throw away input reference
curve.set_dimension(3) # add z-components (if not already present)
curve.force_rational() # add weight (if not already present)
# align axis with the z-axis
normal_theta = atan2(axis[1], axis[0])
normal_phi = atan2(sqrt(axis[0]**2 + axis[1]**2), axis[2])
curve.rotate(-normal_theta, [0,0,1])
curve.rotate(-normal_phi, [0,1,0])
circle_seg = CurveFactory.circle_segment(theta)
n = len(curve) # number of control points of the curve
m = len(circle_seg) # number of control points of the sweep
cp = np.zeros((m * n, 4))
# loop around the circle and set control points by the traditional 9-point
# circle curve with weights 1/sqrt(2), only here C0-periodic, so 8 points
dt = 0
t = 0
for i in range(m):
x,y,w = circle_seg[i]
dt = atan2(y,x) - t
t += dt
curve.rotate(dt)
cp[i * n:(i + 1) * n, :] = curve[:]
cp[i * n:(i + 1) * n, 2] *= w
def sphere(r=1, center=(0,0,0), zaxis=(0,0,1), xaxis=(1,0,0)):
""" Create a spherical shell.
:param float r: Radius
:param array-like center: Local origin of the sphere
:param array-like zaxis: direction of the north/south pole of the parametrization
:param array-like xaxis: direction of the longitudal sem
:return: The spherical shell
:rtype: Surface
"""
circle = CurveFactory.circle_segment(pi, r)
circle.rotate(-pi / 2)
circle.rotate(pi / 2, (1, 0, 0)) # flip up into xz-plane
result = revolve(circle)
result.rotate(rotate_local_x_axis(xaxis, zaxis))
return flip_and_move_plane_geometry(result, center, zaxis)