/[escript]/trunk/pycad/py_src/transformations.py
ViewVC logotype

Annotation of /trunk/pycad/py_src/transformations.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (hide annotations)
Mon Jan 20 03:37:18 2020 UTC (8 months ago) by uqaeller
File MIME type: text/x-python
File size: 5064 byte(s)
Updated the copyright header.


1 ksteube 1809
2 jfenwick 3981 ##############################################################################
3 ksteube 1312 #
4 uqaeller 6939 # Copyright (c) 2003-2020 by The University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 ksteube 1312 #
7 ksteube 1809 # Primary Business: Queensland, Australia
8 jfenwick 6112 # Licensed under the Apache License, version 2.0
9     # http://www.apache.org/licenses/LICENSE-2.0
10 ksteube 1312 #
11 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
13     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 uqaeller 6939 # Development from 2019 by School of Earth and Environmental Sciences
15 jfenwick 3981 #
16     ##############################################################################
17 ksteube 1312
18 sshaw 5706 from __future__ import print_function, division
19    
20 uqaeller 6939 __copyright__="""Copyright (c) 2003-2020 by The University of Queensland
21 jfenwick 3981 http://www.uq.edu.au
22 ksteube 1809 Primary Business: Queensland, Australia"""
23 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
24     http://www.apache.org/licenses/LICENSE-2.0"""
25 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
26 ksteube 1809
27 gross 905 """
28     transformations
29    
30 jfenwick 2625 :var __author__: name of author
31     :var __copyright__: copyrights
32     :var __license__: licence agreement
33     :var __url__: url entry point on documentation
34     :var __version__: version
35     :var __date__: date of the version
36     :var DEG: unit of degree
37     :var RAD: unit of radiant
38 gross 905 """
39    
40     __author__="Lutz Gross, l.gross@uq.edu.au"
41    
42 jfenwick 2455 import numpy
43 gross 905 import math
44    
45 jfenwick 2455 _TYPE=numpy.float64
46 gross 905 DEG=math.pi/180.
47     RAD=1.
48     class Transformation(object):
49     """
50 jfenwick 2625 General class to define an affine transformation *x->Ax+b*.
51 gross 905 """
52     def __init__(self):
53     """
54 caltinay 2180 Creates a linear transformation.
55 gross 905 """
56     pass
57 caltinay 2180
58 jfenwick 2455 def __call__(self,x=numpy.zeros((3,))):
59 gross 905 """
60 jfenwick 2625 Applies transformation to ``x``.
61 gross 905 """
62     raise NotImplementeError()
63 caltinay 2180
64 gross 905 class Translation(Transformation):
65     """
66 jfenwick 2625 Defines a translation *x->x+b*.
67 gross 905 """
68 jfenwick 2455 def __init__(self,b=numpy.zeros((3,),dtype=_TYPE)):
69 gross 905 """
70 jfenwick 2625 Creates the linear transformation *x->x+b*.
71 gross 905 """
72     super(Translation, self).__init__()
73 jfenwick 2455 self.__b=numpy.array(b,_TYPE)
74 caltinay 2180
75 jfenwick 2455 def __call__(self,x=numpy.zeros((3,))):
76 gross 905 """
77 jfenwick 2625 Applies translation to ``x``.
78 gross 905 """
79 jfenwick 2455 return numpy.array(x,_TYPE)+self.__b
80 caltinay 2180
81 jfenwick 3259 class Rotatation(Transformation):
82 gross 905 """
83 caltinay 2180 Defines a rotation.
84 gross 905 """
85 jfenwick 2455 def __init__(self,axis=numpy.ones((3,),dtype=_TYPE),point=numpy.zeros((3,),dtype=_TYPE),angle=0.*RAD):
86 gross 905 """
87 caltinay 2180 Creates a rotation using an axis and a point on the axis.
88 gross 905 """
89 jfenwick 2455 self.__axis=numpy.array(axis,dtype=_TYPE)
90     self.__point=numpy.array(point,dtype=_TYPE)
91     lax=numpy.dot(self.__axis,self.__axis)
92 gross 905 if not lax>0:
93     raise ValueError("points must be distinct.")
94     self.__axis/=math.sqrt(lax)
95     self.__angle=float(angle)
96 caltinay 2180
97 jfenwick 2455 def __call__(self,x=numpy.zeros((3,))):
98 gross 905 """
99 jfenwick 2625 Applies the rotation to ``x``.
100 gross 905 """
101 jfenwick 2455 x=numpy.array(x,_TYPE)
102 gross 916 z=x-self.__point
103 jfenwick 2455 z0=numpy.dot(z,self.__axis)
104 gross 905 z_per=z-z0*self.__axis
105 jfenwick 2455 lz_per=numpy.dot(z_per,z_per)
106 gross 907 if lz_per>0:
107     axis1=z_per/math.sqrt(lz_per)
108 gross 910 axis2=_cross(axis1,self.__axis)
109 jfenwick 2455 lax2=numpy.dot(axis2,axis2)
110 gross 910 if lax2>0:
111     axis2/=math.sqrt(lax2)
112 caltinay 2180 return z0*self.__axis+math.sqrt(lz_per)*(math.cos(self.__angle)*axis1-math.sin(self.__angle)*axis2)+self.__point
113 gross 910 else:
114     return x
115 gross 905 else:
116     return x
117 caltinay 2180
118 gross 907 def _cross(x, y):
119 gross 905 """
120 jfenwick 2625 Returns the cross product of ``x`` and ``y``.
121 gross 905 """
122 jfenwick 2455 return numpy.array([x[1] * y[2] - x[2] * y[1], x[2] * y[0] - x[0] * y[2], x[0] * y[1] - x[1] * y[0]], _TYPE)
123 gross 905
124     class Dilation(Transformation):
125     """
126 caltinay 2180 Defines a dilation.
127 gross 905 """
128 jfenwick 2455 def __init__(self,factor=1.,center=numpy.zeros((3,),dtype=_TYPE)):
129 gross 905 """
130 caltinay 2180 Creates a dilation with a center and a given expansion/contraction
131     factor.
132 gross 905 """
133 gross 915 if not abs(factor)>0:
134     raise ValueError("factor must be non-zero.")
135 gross 905 self.__factor=factor
136 jfenwick 2455 self.__center=numpy.array(center,dtype=_TYPE)
137 caltinay 2180
138 jfenwick 2455 def __call__(self,x=numpy.zeros((3,))):
139 gross 905 """
140 jfenwick 2625 Applies dilation to ``x``.
141 gross 905 """
142 jfenwick 2455 x=numpy.array(x,_TYPE)
143 gross 905 return self.__factor*(x-self.__center)+self.__center
144    
145     class Reflection(Transformation):
146     """
147 caltinay 2180 Defines a reflection on a plane.
148 gross 905 """
149 jfenwick 2455 def __init__(self,normal=numpy.ones((3,),dtype=_TYPE),offset=0.):
150 gross 905 """
151 caltinay 2180 Defines a reflection on a plane defined in normal form.
152 gross 905 """
153 jfenwick 2455 self.__normal=numpy.array(normal,dtype=_TYPE)
154     ln=math.sqrt(numpy.dot(self.__normal,self.__normal))
155 gross 905 if not ln>0.:
156     raise ValueError("normal must have positive length.")
157     self.__normal/=ln
158     if isinstance(offset,float) or isinstance(offset,int):
159     self.__offset=offset/ln
160     else:
161 jfenwick 2455 self.__offset=numpy.dot(numpy.array(offset,dtype=_TYPE),self.__normal)
162 caltinay 2180
163 jfenwick 2455 def __call__(self,x=numpy.zeros((3,))):
164 gross 905 """
165 jfenwick 2625 Applies reflection to ``x``.
166 gross 905 """
167 jfenwick 2455 x=numpy.array(x,_TYPE)
168     return x - 2*(numpy.dot(x,self.__normal)-self.__offset)*self.__normal
169 caltinay 2180

  ViewVC Help
Powered by ViewVC 1.1.26