diff --git a/poisson/continuous/shapes.py b/poisson/continuous/shapes.py
index 1a322fa9ee0bd5daa04d6250c6766cbc31667cd4..d1420ef2e8422af2407dde6cab82573d6717cba7 100644
--- a/poisson/continuous/shapes.py
+++ b/poisson/continuous/shapes.py
@@ -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()