Commit 1c15376a authored by pacome's avatar pacome
Browse files

Merge branch 'clean_shapes'

parents 2dd480d9 486b9d30
......@@ -344,7 +344,31 @@ class Shape(ABC):
return General(func=new_func)
def translate(self, vect):
vect = np.asarray(vect)
func = self.geometry
def newfunc(x):
return func(translate(x, -vect))
self.geometry = newfunc
def rotate(self, angle, axis=None, center='self'):
'''3D not tested !!!'''
func = self.geometry
try:
angle = -angle
except TypeError:
angle = -np.asarray(angle)
def newfunc(x):
return func(rotate(x, angle, axis=axis, center=center))
self.geometry = newfunc
def reflect(self, axis):
func = self.geometry
def newfunc(x):
return func(reflect(x, axis))
self.geometry = newfunc
(),
class PiecewiseFunction(object):
'''
When called apllies self.evaluate_coordinates to a vector of
......@@ -461,13 +485,12 @@ class General(Shape):
def geometry(self, x):
return self.func(x)
class Rectangle(General):
'''
Define a rectangle
'''
def __init__(self, length, corner=None):
def __init__(self, length, corner=None, center=None):
'''
Class that defines a rectangular shape.
When called and a numpy array is given as parameter it returns a boolean
......@@ -482,6 +505,11 @@ class Rectangle(General):
the lengths of the sides
if Ly (Lz) is None, Ly (Lz) will be equal to Lx
(usefull to quickly draw cubes)
center = (x0, y0, z0) : numbers
overwrites the corner argument
if center=None the corner will be used
the lower left corner of the (if z0 si None will be 2d)
Returns:
--------
func: function
......@@ -491,9 +519,9 @@ class Rectangle(General):
super().__init__()
self.length=length
self.corner=corner
self.center=center
self.prepare_param()
def prepare_param(self):
if isinstance(self.length, (int, float)):
......@@ -508,6 +536,9 @@ class Rectangle(General):
self.corner = np.asanyarray(self.corner)
if self.center is not None:
self.center = np.asarray(self.center)
assert len(self.length) == len(self.corner), (
'The corner does not have the same'
+ ' size as length / or a corner has'
......@@ -520,37 +551,25 @@ class Rectangle(General):
It returns True if the point is within the shape defined
by the function or False otherwise
'''
x = np.asarray(x)
if len(x.shape) == 1:
ndim, npt = x.shape[0], -1
x = x[None, :]
else:
npt, ndim = x.shape
assert (ndim == 2 or ndim == 3)
x = np.asanyarray(x)
if len(x.shape) > 1:
if len(self.length) == 3:
return ((self.corner[0] <= x[:, 0])
* (x[:, 0] <= self.corner[0] + self.length[0])
* (self.corner[1] <= x[:, 1])
* (x[:, 1] <= self.corner[1] + self.length[1])
* (self.corner[2] <= x[:, 2])
* (x[:, 2] <= self.corner[2] + self.length[2]))
else:
return ((self.corner[0] <= x[:, 0])
* (x[:, 0] <= self.corner[0] + self.length[0])
* (self.corner[1] <= x[:, 1])
* (x[:, 1] <= self.corner[1] + self.length[1]))
if self.center is None:
isin = np.all(x <= self.corner + self.length, axis=1) \
* np.all(x >= self.corner, axis=1)
else:
if len(self.length) == 3:
return ((self.corner[0] <= x[0]
<= self.corner[0] + self.length[0])
and (self.corner[1] <= x[1]
<= self.corner[1] + self.length[1])
and (self.corner[2] <= x[2]
<= self.corner[2] + self.length[2]))
else:
return ((self.corner[0] <= x[0]
<= self.corner[0] + self.length[0])
and (self.corner[1] <= x[1]
<= self.corner[1] + self.length[1]))
isin = np.all(x <= self.center + self.length/2, axis=1) \
* np.all(x >= self.center - self.length/2, axis=1)
if npt is -1:
return isin[0]
else:
return isin
class Ellipsoid(General):
'''
......@@ -560,7 +579,7 @@ class Ellipsoid(General):
def __init__(self, radius, center=None):
'''
Class that defines a ellipsoid shape.
When called and a numpy array is given as parameter it returns a
When called and a numpy array is given as parameter it retu(),rns a
boolean numpy array of the same length of the latter.
Returns true if inside an 2d of 3d ellipsoid
......@@ -694,20 +713,10 @@ class InHull(Shape):
'''
if self.points_coordinates is None:
self.points_coordinates = self.generate_points()
vect_parent = vect
def fun():
vect = np.asarray(vect_parent)
if len(self.points_coordinates.shape) < 2:
# translate a point
return self.points_coordinates + np.asarray(vect)
elif len(self.points_coordinates.shape) == 2:
return (np.repeat([vect], len(self.points_coordinates), axis=0)
+ self.points_coordinates)
self.points_coordinates = fun()
self.points_operations.append(fun)
self.points_coordinates = translate(self.points_coordinates, vect)
self.points_operations.append(
lambda : translate(self.points_coordinates, vect))
def rotate(self, angle, axis=None, center='self'):
'''Rotate a point or a list of points in 2 or 3d around a center
......@@ -720,7 +729,7 @@ class InHull(Shape):
angle of rotation.
if 3d, tx is either the angle around the x axis or
the angle to rotate around the specified axis
ty: number
ty: numbernp.repeat([rot_mat], len(x), axis=0)
if axis is None, the angle of the rotation around y
will be ignored if axis is not None
tz: number (opt, defaults to None)
......@@ -739,49 +748,229 @@ class InHull(Shape):
if self.points_coordinates is None:
self.points_coordinates = self.generate_points()
if not isinstance(angle, (list, tuple)):
angle = (angle,)
angle = np.asarray(angle)
center_parent = center
def fun():
if len(self.points_coordinates.shape) == 2:
npt, ndim = self.points_coordinates.shape
else: raise ValueError
if ndim is 2:
rot_mat = np.array([[np.cos(angle[0]), -np.sin(angle[0])],
[np.sin(angle[0]), np.cos(angle[0])]])
elif ndim is 3:
rot_mat = _rotation_matrix(angle, axis)
else: raise ValueError
if center_parent is None:
center = np.zeros(ndim)
elif center_parent is 'self':
center = np.mean(self.points_coordinates, axis=0)
elif isinstance(center_parent, (tuple, list, np.ndarray)):
center = np.asarray(center_parent)
else: raise ValueError
center_mat = np.repeat([center],
len(self.points_coordinates), axis=0)
ret = (np.einsum(
'kij, kj -> ki',
np.repeat([rot_mat], len(self.points_coordinates), axis=0),
self.points_coordinates - center_mat)
+ center_mat)
if npt is 1:
return ret[0]
else:
return ret
self.points_coordinates = rotate(self.points_coordinates, angle,
axis=axis, center=center)
self.points_operations.append(
lambda : rotate(self.points_coordinates, angle,
axis=axis, center=center))
def reflect(self, axis):
'''returns reflections of point(s) x wrt to axis
Parameters:
-----------
axis: array like
2 (3) points that generate the line (plan)
wrt which to reflect in 2d (3d)
'''
if self.points_coordinates is None:
self.points_coordinates = self.generate_points()
self.points_coordinates = reflect(self.points_coordinates, axis)
self.points_operations.append(lambda : reflect(self.points_coordinates,
axis))
#####################
def box(length, corner=None, center=None):
'''
Returns points of a 2d of 3d box
Parameters:
-----------
length = (Lx, Ly, Lz): numbers
the lengths of the sides
if Ly (Lz) is None, Ly (Lz) will be equal to Lx
(usefull to quickly draw cubes)
corner = (x0, y0, z0) : numbers
the lower left corner of the (if z0 si None will be 2d)
center = (x0, y0, z0) : numbers
overwrites the corner argument
if center=None the corner will be used
the lower left corner of the (if z0 si None will be 2d)
Returns:
--------
points: np.ndarray
the points defining the box
'''
if isinstance(length, (int, float)):
length = (length, length)
self.points_coordinates = fun()
self.points_operations.append(fun)
length = np.asarray(length)
ndim = len(length)
if corner is None:
corner = np.zeros(ndim)
elif isinstance(corner, int):
corner = corner * np.ones(ndim)
corner = np.asarray(corner)
if center is not None:
center = np.asarray(center)
assert len(length) == len(corner), (
'The corner does not have the same'
+ ' size as length / or a corner has'
+ ' been given and length is an int')
points = np.array([[0, 0], [0, length[1]], length[:2], [length[0], 0]])
if ndim == 3:
points = np.vstack((np.c_[points, np.zeros(4)],
np.c_[points, np.ones(4) * length[2]]))
if center is None:
points = points + corner
else:
points = points - length / 2 + center
return points
def translate(x, vect):
'''
Translate a point or a list of points by a vector vect
Parameters:
-----------
x: tuple or liste of tuple
the points to translate
vect: list, tuple or numpy array
the translation vector
'''
x = np.asarray(x)
vect = np.asarray(vect)
if len(x.shape) < 2:
# translate a point
return x + np.asarray(vect)
elif len(x.shape) == 2:
return np.repeat([vect], len(x), axis=0) + x
def rotate(x, angle, axis=None, center='self'):
'''Rotate a point or a list of points in 2 or 3d around a center
of an angle or a combination of angles or around an axis
Parameters:
-----------
angle: (tx, ty, tz) or number
tx: number
angle of rotation.
if 3d, tx is either the angle around the x axis or
the angle to rotate around the specified axis
ty: number
if axis is None, the angle of the rotation around y
will be ignored if axis is not None
tz: number (opt, defaults to None)
if axis is None, the angle of the rotation around z
will be ignored if axis is not None
axis: arraylike (opt, defaults to None)
if axis is None, rotation around main axes
if not None, the axis of rotation
center: None, arraylike or 'self'(opt, defaults to None)
the center of the rotation
if None, will be the origin according to the dimension
if self, will be the average of the points
(center of regular polygon for 'inplace' rotation)
'''
x = np.asarray(x)
center_parent = center
if not isinstance(angle, (list, tuple)):
angle = (angle,)
angle = np.asarray(angle)
if len(x.shape) == 2:
npt, ndim = x.shape
elif len(x.shape) == 1:
npt, ndim = -1, x.shape[0]
x = x[None, :]
else: raise ValueError
if ndim is 2:
rot_mat = np.array([[np.cos(angle[0]), -np.sin(angle[0])],
[np.sin(angle[0]), np.cos(angle[0])]])
elif ndim is 3:
rot_mat = _rotation_matrix(angle, axis)
else: raise ValueError
if center_parent is None:
center = np.zeros(ndim)
elif center_parent is 'self':
center = np.mean(x, axis=0)
elif isinstance(center_parent, (tuple, list, np.ndarray)):
center = np.asarray(center_parent)
else: raise ValueError
center_mat = np.repeat([center],
len(x), axis=0)
ret = (np.einsum(
'kij, kj -> ki',
np.repeat([rot_mat], len(x), axis=0), x - center_mat)
+ center_mat)
if npt is -1:
return ret[0]
else:
return ret
def reflect(x, axis):
'''returns reflections of point(s) x wrt to axis
Parameters:
-----------
x: array like
point or list of points in 2 or 3d
axis: array like
2 (3) points that generate the line (plan)
wrt which to reflect in 2d (3d)
Returns:
--------
ret: array like (shape of x)
the reflected point(s)
'''
x = np.asarray(x)
if len(x.shape) == 1:
ndim, npt = x.shape[0], -1
x = x[None, :]
else:
npt, ndim = x.shape
assert (ndim == 2 or ndim == 3)
assert isinstance(axis, (list, np.ndarray))
axis = np.asarray(axis)
assert axis.shape[1] == ndim
A = axis[0]
Amat = np.repeat([A], len(x), axis=0)
if ndim == 2:
u = (axis[1:] - A)[0]
v = np.array([u[1], -u[0]])
# u, v is othogonal basis
mirror = np.array([[1, 0],[0, -1]])
Mat = np.array([u, v])
elif ndim == 3:
u, C = axis[1:] - A
w = np.cross(u, C)
v = np.cross(w, u)
# u, v, w is othogonal basis
mirror = np.array([[1, 0, 0],[0, 1, 0], [0 ,0 ,-1]])
Mat = np.array([u, v, w])
Matinv = np.linalg.inv(Mat)
transf = Mat.T @ mirror @ Matinv.T
ret = np.einsum('kij, kj -> ki', np.repeat([transf], len(x), axis=0),
(x - Amat)) + Amat
if npt is -1:
return ret[0]
else:
return ret
#################
class Delaunay(InHull):
'''
......@@ -833,7 +1022,7 @@ class RegularPolygon(InHull):
Parameters:
-----------
n: integer (opt, defaults to 3)
npoints: integer (opt, defaults to 3)
the number of sides (points) of the regular polygon
center: 2d tuple (opt, defaults to [0, 0])
the center of the regular polygon
......@@ -845,20 +1034,17 @@ class RegularPolygon(InHull):
flat: bool (opt, defaults to True)
it True, the polygon is rotated to have his lowest
edge parallel to x axis
'''
super().__init__()
self.radius = radius
self.npoints = npoints
self.center = center
self.clockwise = clockwise
self.flat = flat
super().__init__()
def generate_points(self):
'''
Returns a numpy array containing the coordinates of the points
defining the shape.
'''
......@@ -918,14 +1104,13 @@ class ExtrudedRegularPolygon(InHull):
'''
super().__init__()
self.radius = radius
self.npoints = npoints
self.center = center
self.clockwise = clockwise
self.flat = flat
self.Lz = Lz
super().__init__()
def generate_points(self):
'''
......@@ -961,11 +1146,11 @@ class ExtrudedRegularPolygon(InHull):
##########################################################
# test plots
def __test_plot_ellipse():
ellipse = Ellipsoid(radius=2)
mm = continuous.grid._2d_mesh(*continuous.grid._points_from_2dbounding_box([-3, 3, -3, 3],
0.1))
mm = continuous.grid._2d_mesh(
*continuous.grid._points_from_2dbounding_box([-3, 3, -3, 3], 0.1))
x, y = list(zip(*mm))
xx, yy = list(zip(*continuous.grid._apply_stencil(mm, ellipse)))
plt.scatter(x, y, label='bbox')
......@@ -994,6 +1179,26 @@ def __test_plot_rectangle():
plt.scatter(x, y, c=mesh.point_label)
plt.show()
def __test_plot_rectangle_center():
rect1 = Rectangle(length=4, center=[2, 2])
rect2 = Rectangle(length=2, center=[2, 2])
rect3 = Rectangle(length=1, center=[2, 2])
bbox1 = [0, 4, 0, 4]
bbox2 = [1, 3, 1, 3]
mesh1 = (bbox1, 0.15, rect1, 0)
mesh2 = (bbox2, 0.1, rect2, 1)
hole1 = (rect3)
points = ([[2, 2]], 2)
mesh = GridBuilder(meshs=[mesh1, mesh2], holes=[hole1],
points=[points])
x, y = list(zip(*mesh.points))
plt.scatter(x, y, c=mesh.point_label)
plt.show()
def __test_plot_ellipse_circ():
......@@ -1062,6 +1267,34 @@ def __test_plot_rotation_translation():
plt.axes().set_aspect('equal', 'datalim')
plt.show()
def __test_reflection():
'''test if two reflexions is equal to identity'''
np.random.seed(112)
ndim = 2
axis = np.random.rand(ndim, ndim)
pt = np.random.rand(150, ndim)
assert(np.allclose(pt, reflect(reflect(pt, axis), axis)))
ndim = 3
axis = np.random.rand(ndim, ndim)
pt = np.random.rand(150, ndim)
assert(np.allclose(pt, reflect(reflect(pt, axis), axis)))
print('test reflexion passed')
def __test_plot_reflection():
ndim = 2
axis = np.random.rand(ndim, ndim) - 0.5
pt = np.random.rand(10, ndim)
newpt = reflect(pt, axis)
plt.figure(figsize=(12,12))
plt.plot(*pt.T, marker='o')
plt.plot(*axis.T, lw=4)
plt.plot(*newpt.T, marker='^')
plt.axes().set_aspect('equal')
plt.show()
#test_plot_ellipse()
#test_plot_rectangle()
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment