1 |
# $Id:$ |
2 |
|
3 |
""" |
4 |
Geometrical Primitives |
5 |
|
6 |
the concept is inspired by gmsh and very much focused on the fact that |
7 |
the classes are used to wrk with gmsh. |
8 |
|
9 |
@var __author__: name of author |
10 |
@var __copyright__: copyrights |
11 |
@var __license__: licence agreement |
12 |
@var __url__: url entry point on documentation |
13 |
@var __version__: version |
14 |
@var __date__: date of the version |
15 |
""" |
16 |
|
17 |
|
18 |
__author__="Lutz Gross, l.gross@uq.edu.au" |
19 |
__copyright__=""" Copyright (c) 2006 by ACcESS MNRF |
20 |
http://www.access.edu.au |
21 |
Primary Business: Queensland, Australia""" |
22 |
__license__="""Licensed under the Open Software License version 3.0 |
23 |
http://www.opensource.org/licenses/osl-3.0.php""" |
24 |
__url__="http://www.iservo.edu.au/esys/escript" |
25 |
__version__="$Revision:$" |
26 |
__date__="$Date:$" |
27 |
|
28 |
import numarray |
29 |
from transformations import _TYPE |
30 |
|
31 |
global global_primitive_id_counter |
32 |
global_primitive_id_counter=1 |
33 |
|
34 |
class Primitive(object): |
35 |
""" |
36 |
template for elementary geometrical object |
37 |
""" |
38 |
def __init__(self): |
39 |
""" |
40 |
""" |
41 |
global global_primitive_id_counter |
42 |
self.__ID=global_primitive_id_counter |
43 |
global_primitive_id_counter+=1 |
44 |
def getID(self): |
45 |
return self.__ID |
46 |
def __repr__(self): |
47 |
return "%s(%s)"%(self.__class__.__name__,self.getID()) |
48 |
def __cmp__(self,other): |
49 |
return cmp(self.getID(),other.getID()) |
50 |
def getPoints(self): |
51 |
""" |
52 |
returns the C{set} of points used to construct the primitive |
53 |
""" |
54 |
out=set() |
55 |
for i in self.getHistory(): |
56 |
if isinstance(i,Point): out.add(i) |
57 |
return out |
58 |
|
59 |
def setLocalScale(self,factor=1.): |
60 |
""" |
61 |
sets the local refinement factor |
62 |
""" |
63 |
for p in self.getPoints(): p.setLocalScale(factor) |
64 |
|
65 |
def isPoint(self): |
66 |
""" |
67 |
returns C{True} is the primitive is a L{Point} |
68 |
""" |
69 |
return False |
70 |
def isCurve(self): |
71 |
""" |
72 |
returns C{True} is the primitive is a L{Curve} |
73 |
""" |
74 |
return False |
75 |
def isSurface(self): |
76 |
""" |
77 |
returns C{True} is the primitive is a L{Surface} |
78 |
""" |
79 |
return False |
80 |
def isCurveLoop(self): |
81 |
""" |
82 |
returns C{True} is the primitive is a L{CurveLoop} |
83 |
""" |
84 |
return False |
85 |
def isSurfaceLoop(self): |
86 |
""" |
87 |
returns C{True} is the primitive is a L{SurfaceLoop} |
88 |
""" |
89 |
return False |
90 |
def getHistory(self): |
91 |
""" |
92 |
returns C{set} of primitive used to construct the primitive |
93 |
""" |
94 |
return set() |
95 |
|
96 |
|
97 |
#================================================== |
98 |
def __neg__(self): |
99 |
return ReversedPrimitive(self) |
100 |
def __pos__(self): |
101 |
return self.copy() |
102 |
def __add__(self,other): |
103 |
out=self.copy() |
104 |
out+=other |
105 |
return out |
106 |
def __iadd__(self,other): |
107 |
self.shift() |
108 |
def shift(self,shift): |
109 |
for p in self.getPoints(): p+=shift |
110 |
def copy(self): |
111 |
return Primitive() |
112 |
def getGmshCommand(self): |
113 |
raise NotImplementedError("getGmshCommand is not implemented for this class %s."%self.__class__.__name__) |
114 |
def translate(self,shift): |
115 |
raise NotImplementedError("translate is not implemented for this class %s."%self.__class__.__name__) |
116 |
|
117 |
class Point(Primitive): |
118 |
""" |
119 |
a three dimensional point |
120 |
""" |
121 |
def __init__(self,x=0.,y=0.,z=0.,local_scale=1.): |
122 |
""" |
123 |
creates a point with coorinates x,y,z with the local refinement factor local_scale |
124 |
""" |
125 |
super(Point, self).__init__() |
126 |
self.setCoordinates(x,y,z) |
127 |
self.setLocalScale(local_scale) |
128 |
def setLocalScale(self,factor=1.): |
129 |
""" |
130 |
sets the local refinement factor |
131 |
""" |
132 |
if factor<=0.: |
133 |
raise ValueError("scaling factor must be positive.") |
134 |
self.__local_scale=factor |
135 |
def getLocalScale(self): |
136 |
""" |
137 |
returns the local refinement factor |
138 |
""" |
139 |
return self.__local_scale |
140 |
def getCoordinates(self): |
141 |
""" |
142 |
returns the coodinates of the point as L{numarray.NumArray} object |
143 |
""" |
144 |
return self._x |
145 |
def setCoordinates(self,x,y,z): |
146 |
""" |
147 |
returns the coodinates of the point as L{numarray.NumArray} object |
148 |
""" |
149 |
self._x=numarray.array([x,y,z],_TYPE) |
150 |
def getHistory(self): |
151 |
""" |
152 |
returns C{set} of primitive used to construct the primitive |
153 |
""" |
154 |
return set([self]) |
155 |
|
156 |
def isColocated(self,point,tol=1.e-11): |
157 |
""" |
158 |
returns True if L{Point} point is colocation (same coordinates) |
159 |
that means if |self-point| <= tol * max(|self|,|point|) |
160 |
""" |
161 |
if isinstance(point,Point): |
162 |
point=point.getCoordinates() |
163 |
c=self.getCoordinates() |
164 |
d=c-point |
165 |
return numarray.dot(d,d)<=tol**2*max(numarray.dot(c,c),numarray.dot(point,point)) |
166 |
|
167 |
|
168 |
#============================================================= |
169 |
def copy(self): |
170 |
c=self.getCoordinates() |
171 |
return Point(c[0],c[1],c[2],local_scale=self.getLocalScale()) |
172 |
def isPoint(self): |
173 |
return True |
174 |
def getGmshCommand(self): |
175 |
c=self.getCoordinates() |
176 |
return "Point(%s) = {%e , %e, %e , %e * scale};"%(self.getID(),c[0],c[1],c[2], self.getLocalScale()) |
177 |
def shift(self,shift): |
178 |
""" |
179 |
shifts the point by a given shift |
180 |
""" |
181 |
self._x+=numarray.array(shift,numarray.Float64) |
182 |
def translate(self,shift): |
183 |
""" |
184 |
returns the point shifted by shift |
185 |
""" |
186 |
out=self.copy() |
187 |
out+=other |
188 |
return out |
189 |
|
190 |
|
191 |
class Curve(Primitive): |
192 |
""" |
193 |
a curve |
194 |
""" |
195 |
def __init__(self,*args): |
196 |
""" |
197 |
defines a curve form a set of control points |
198 |
""" |
199 |
super(Curve, self).__init__() |
200 |
l=len(args) |
201 |
for i in range(l): |
202 |
if not args[i].isPoint(): |
203 |
raise TypeError("%s-th argument is not a Point object."%i) |
204 |
self.__nodes=args |
205 |
def __len__(self): |
206 |
return len(self.__nodes) |
207 |
def isCurve(self): |
208 |
return True |
209 |
def getStart(self): |
210 |
""" |
211 |
returns start point |
212 |
""" |
213 |
return self.__nodes[0] |
214 |
|
215 |
def getEnd(self): |
216 |
""" |
217 |
returns end point |
218 |
""" |
219 |
return self.__nodes[-1] |
220 |
|
221 |
def getNodes(self): |
222 |
""" |
223 |
returns a list of the nodes |
224 |
""" |
225 |
return self.__nodes |
226 |
def getGmshCommand(self): |
227 |
out="" |
228 |
for i in self.getNodes(): |
229 |
if len(out)>0: |
230 |
out+=", %s"%i.getID() |
231 |
else: |
232 |
out="%s"%i.getID() |
233 |
return "Spline(%s) = {%s};"%(self.getID(),out) |
234 |
def getHistory(self): |
235 |
out=set([self]) |
236 |
for i in self.getNodes(): out|=i.getHistory() |
237 |
return out |
238 |
def getPoints(self): |
239 |
out=set() |
240 |
for i in self.getNodes(): out|=i.getPoints() |
241 |
return out |
242 |
|
243 |
|
244 |
class BezierCurve(Curve): |
245 |
""" |
246 |
a Bezier curve |
247 |
""" |
248 |
def __neg__(self): |
249 |
""" |
250 |
returns the line segment with swapped start and end points |
251 |
""" |
252 |
return BezierCurve(self.getNodes()[::-1]) |
253 |
def __add__(self,other): |
254 |
return BezierCurve([p+other for p in self.getNodes()]) |
255 |
def getGmshCommand(self): |
256 |
out="" |
257 |
for i in self.getNodes(): |
258 |
if len(out)>0: |
259 |
out+=", %s"%i.getID() |
260 |
else: |
261 |
out="%s"%i.getID() |
262 |
return "Bezier(%s) = {%s};"%(self.getID(),out) |
263 |
|
264 |
class BSplineCurve(Curve): |
265 |
""" |
266 |
a BSpline curve. Control points may be repeated. |
267 |
""" |
268 |
def __neg__(self): |
269 |
""" |
270 |
returns the line segment with swapped start and end points |
271 |
""" |
272 |
return BSplineCurve(self.getNodes()[::-1]) |
273 |
def __add__(self,other): |
274 |
return BSplineCurve([p+other for p in self.getNodes()]) |
275 |
def getGmshCommand(self): |
276 |
out="" |
277 |
for i in self.getNodes(): |
278 |
if len(out)>0: |
279 |
out+=", %s"%i.getID() |
280 |
else: |
281 |
out="%s"%i.getID() |
282 |
return "BSpline(%s) = {%s};"%(self.getID(),out) |
283 |
|
284 |
class Line(Curve): |
285 |
""" |
286 |
a line is defined by two L{Point}s |
287 |
""" |
288 |
def __init__(self,start,end): |
289 |
""" |
290 |
defines a curve form a set of control points |
291 |
""" |
292 |
super(Line, self).__init__(start,end) |
293 |
def __neg__(self): |
294 |
return ReversedPrimitive(self) |
295 |
def __add__(self,other): |
296 |
return Line(self.getEnd()+other,self.getStart()+other) |
297 |
def getGmshCommand(self): |
298 |
return "Line(%s) = {%s, %s};"%(self.getID(),self.getStart().getID(),self.getEnd().getID()) |
299 |
|
300 |
class Arc(Curve): |
301 |
""" |
302 |
defines an arc |
303 |
""" |
304 |
def __init__(self,center,start,end): |
305 |
""" |
306 |
creates an arc by the start point, end point and center |
307 |
""" |
308 |
if center.isPoint(): |
309 |
raise TypeError("center needs to be a Point object.") |
310 |
super(Arc, self).__init__(start,end) |
311 |
self.__center=center |
312 |
|
313 |
def getCenter(self): |
314 |
""" |
315 |
returns center |
316 |
""" |
317 |
return self.__center |
318 |
def __add__(self,other): |
319 |
return Arc(self.getCenter()+other,self.getStart()+other,self.getEnd()+other) |
320 |
|
321 |
def getHistory(self): |
322 |
out=set([self]) |
323 |
out|=self.getCenter().getHistory() |
324 |
for i in self.getNodes(): out|=i.getHistory() |
325 |
return out |
326 |
def getPoints(self): |
327 |
out=self.getCenter().getPoints() |
328 |
for i in self.getNodes(): out|=i.getPoints() |
329 |
return out |
330 |
def getGmshCommand(self): |
331 |
return "Circle(%s) = {%s, %s, %s};"%(self.getID(),self.getStart().getID(),self.getCenter().getID(),self.getEnd().getID()) |
332 |
|
333 |
class CurveLoop(Primitive): |
334 |
""" |
335 |
An oriented loop of curves. |
336 |
|
337 |
The loop must be closed and the L{Curves}s should be oriented consistently. |
338 |
""" |
339 |
def __init__(self,*curves): |
340 |
""" |
341 |
creates a polygon from a list of line curves. The curves must form a closed loop. |
342 |
""" |
343 |
super(CurveLoop, self).__init__() |
344 |
self.__curves=[] |
345 |
self.addCurve(*curves) |
346 |
def addCurve(self,*curves): |
347 |
for i in range(len(curves)): |
348 |
if not curves[i].isCurve(): |
349 |
raise TypeError("%s-th argument is not a Curve object."%i) |
350 |
self.__curves+=curves |
351 |
|
352 |
def isCurveLoop(self): |
353 |
return True |
354 |
def getCurves(self): |
355 |
return self.__curves |
356 |
def __add__(self,other): |
357 |
return CurveLoop(*tuple([c+other for c in self.getCurves()[::-1]])) |
358 |
def __len__(self): |
359 |
return len(self.__curves) |
360 |
def getHistory(self): |
361 |
out=set([self]) |
362 |
for i in self.getCurves(): out|=i.getHistory() |
363 |
return out |
364 |
def getPoints(self): |
365 |
out=set() |
366 |
for i in self.getCurves(): out|=i.getPoints() |
367 |
return out |
368 |
def getGmshCommand(self): |
369 |
out="" |
370 |
for i in self.getCurves(): |
371 |
if len(out)>0: |
372 |
out+=", %s"%i.getID() |
373 |
else: |
374 |
out="%s"%i.getID() |
375 |
return "Line Loop(%s) = {%s};"%(self.getID(),out) |
376 |
|
377 |
class Surface(Primitive): |
378 |
""" |
379 |
a surface |
380 |
""" |
381 |
def __init__(self,loop): |
382 |
""" |
383 |
creates a surface with boundary loop |
384 |
|
385 |
@param loop: L{CurveLoop} defining the boundary of the surface |
386 |
""" |
387 |
super(Surface, self).__init__() |
388 |
if not loop.isCurveLoop(): |
389 |
raise TypeError("argument loop needs to be a CurveLoop object.") |
390 |
self.__loop=loop |
391 |
def isSurface(self): |
392 |
return True |
393 |
def getBoundaryLoop(self): |
394 |
return self.__loop |
395 |
def __add__(self,other): |
396 |
return Surface(self.getBoundaryLoop()+other) |
397 |
def getHistory(self): |
398 |
out=set([self]) | self.getBoundaryLoop().getHistory() |
399 |
return out |
400 |
def getPoints(self): |
401 |
return self.getBoundaryLoop().getPoints() |
402 |
def getGmshCommand(self): |
403 |
return "Ruled Surface(%s) = {%s};"%(self.getID(),self.getBoundaryLoop().getID()) |
404 |
|
405 |
class PlaneSurface(Surface): |
406 |
""" |
407 |
a plane surface with holes |
408 |
""" |
409 |
def __init__(self,loop,holes=[]): |
410 |
""" |
411 |
creates a plane surface. |
412 |
|
413 |
@param loop: L{CurveLoop} defining the boundary of the surface |
414 |
@param holes: list of L{CurveLoop} defining holes in the surface. |
415 |
@note: A CurveLoop defining a hole should not have any lines in common with the exterior CurveLoop. |
416 |
A CurveLoop defining a hole should not have any lines in common with another CurveLoop defining a hole in the same surface. |
417 |
""" |
418 |
super(PlaneSurface, self).__init__(loop) |
419 |
for i in range(len(holes)): |
420 |
if not holes[i].inCurveLoop(): |
421 |
raise TypeError("%i th hole needs to be a CurveLoop object.") |
422 |
self.__holes=holes |
423 |
def getHoles(self): |
424 |
return self.__holes |
425 |
def __add__(self,other): |
426 |
return PlaneSurface(self.getBoundaryLoop()+other, holes=[h+other for h in self.getHoles()]) |
427 |
def getHistory(self): |
428 |
out=set([self]) | self.getBoundaryLoop().getHistory() |
429 |
for i in self.getHoles(): out|=i.getHistory() |
430 |
return out |
431 |
def getPoints(self): |
432 |
out=self.getBoundaryLoop().getPoints() |
433 |
for i in self.getHoles(): out|=i.getPoints() |
434 |
return out |
435 |
def getGmshCommand(self): |
436 |
out="" |
437 |
for i in self.getHoles(): |
438 |
if len(out)>0: |
439 |
out+=", %s"%i.getID() |
440 |
else: |
441 |
out="%s"%i.getID() |
442 |
if len(out)>0: |
443 |
return "Plane Surface(%s) = {%s, %s};"%(self.getID(),self.getBoundaryLoop().getID(), out) |
444 |
else: |
445 |
return "Plane Surface(%s) = {%s};"%(self.getID(),self.getBoundaryLoop().getID()) |
446 |
|
447 |
class RuledSurface(Surface): |
448 |
""" |
449 |
A ruled surface, i.e., a surface that can be interpolated using transfinite interpolation |
450 |
""" |
451 |
def __init__(self,loop): |
452 |
""" |
453 |
creates a ruled surface from a |
454 |
|
455 |
@param loop: L{CurveLoop} defining the boundary of the surface. There is a restriction of composed of either three or four L{Curve} objects. |
456 |
""" |
457 |
if not loop.isCurveLoop(): |
458 |
raise TypeError("argument loop needs to be a CurveLoop object.") |
459 |
if len(loop)<3: |
460 |
raise TypeError("the loop must contain at least three Curves.") |
461 |
super(RuledSurface, self).__init__(loop) |
462 |
def __add__(self,other): |
463 |
return RuledSurface(self.getBoundaryLoop()+other) |
464 |
def getGmshCommand(self): |
465 |
return "Ruled Surface(%s) = {%s};"%(self.getID(),self.getBoundaryLoop().getID()) |
466 |
|
467 |
class SurfaceLoop(Primitive): |
468 |
""" |
469 |
a surface loop. It defines the shell of a volume. |
470 |
|
471 |
The loop must represent a closed shell, and the L{Surface}s should be oriented consistently. |
472 |
""" |
473 |
def __init__(self,*surfaces): |
474 |
""" |
475 |
creates a surface loop |
476 |
""" |
477 |
super(SurfaceLoop, self).__init__() |
478 |
self.__surfaces=[] |
479 |
self.addSurface(*surfaces) |
480 |
def addSurface(self,*surfaces): |
481 |
for i in range(len(surfaces)): |
482 |
if not surfaces[i].isSurface(): |
483 |
raise TypeError("%s-th argument is not a Surface object."%i) |
484 |
self.__surfaces+=surfaces |
485 |
|
486 |
def isSurfaceLoop(self): |
487 |
return True |
488 |
def getSurfaces(self): |
489 |
return self.__surfaces |
490 |
def __add__(self,other): |
491 |
return SurfaceLoop([c+other for c in self.getSurfaces]) |
492 |
def __len__(self): |
493 |
return len(self.__surfaces) |
494 |
def getHistory(self): |
495 |
out=set([self]) |
496 |
for i in self.getSurfaces(): out|=i.getHistory() |
497 |
return out |
498 |
def getPoints(self): |
499 |
out=set() |
500 |
for i in self.getSurfaces(): out|=i.getPoints() |
501 |
return out |
502 |
def getGmshCommand(self): |
503 |
out="" |
504 |
for i in self.getSurfaces(): |
505 |
if len(out)>0: |
506 |
out+=", %s"%i.getID() |
507 |
else: |
508 |
out="%s"%i.getID() |
509 |
return "Surface Loop(%s) = {%s};"%(self.getID(),out) |
510 |
|
511 |
class Volume(Primitive): |
512 |
""" |
513 |
a volume with holes. |
514 |
""" |
515 |
def __init__(self,loop,holes=[]): |
516 |
""" |
517 |
creates a volume |
518 |
|
519 |
@param loop: L{SurfaceLoop} defining the boundary of the surface |
520 |
@param holes: list of L{SurfaceLoop} defining holes in the surface. |
521 |
@note: A SurfaceLoop defining a hole should not have any surfaces in common with the exterior SurfaceLoop. |
522 |
A SurfaceLoop defining a hole should not have any surfaces in common with another SurfaceLoop defining a hole in the same volume. |
523 |
""" |
524 |
super(Volume, self).__init__() |
525 |
if not loop.isSurfaceLoop(): |
526 |
raise TypeError("argument loop needs to be a SurfaceLoop object.") |
527 |
for i in range(len(holes)): |
528 |
if not holes[i].isSurfaceLoop(): |
529 |
raise TypeError("%i th hole needs to be a SurfaceLoop object.") |
530 |
self.__loop=loop |
531 |
self.__holes=holes |
532 |
def getHoles(self): |
533 |
return self.__holes |
534 |
def getSurfaceLoop(self): |
535 |
return self.__loop |
536 |
def __add__(self,other): |
537 |
return Volume(self.getSurfaceLoop()+other, holes=[h+other for h in self.getHoles()]) |
538 |
def getHistory(self): |
539 |
out=set([self]) | self.getSurfaceLoop().getHistory() |
540 |
for i in self.getHoles(): out|=i.getHistory() |
541 |
return out |
542 |
def getPoints(self): |
543 |
out=self.getSurfaceLoop().getPoints() |
544 |
for i in self.getHoles(): out|=i.Points() |
545 |
return out |
546 |
def getGmshCommand(self): |
547 |
out="" |
548 |
for i in self.getHoles(): |
549 |
if len(out)>0: |
550 |
out+=", %s"%i.getID() |
551 |
else: |
552 |
out="%s"%i.getID() |
553 |
if len(out)>0: |
554 |
return "Volume(%s) = {%s, %s};"%(self.getID(),self.getSurfaceLoop().getID(), out) |
555 |
else: |
556 |
return "Volume(%s) = {%s};"%(self.getID(),self.getSurfaceLoop().getID()) |
557 |
|
558 |
class ReversedPrimitive(object): |
559 |
def __init__(self,prim): |
560 |
self.__prim=prim |
561 |
def __getattr__(self,name): |
562 |
if name == "getID": |
563 |
return self.getReverseID |
564 |
else: |
565 |
return getattr(self.__prim,name) |
566 |
def getReverseID(self): |
567 |
return -self.__prim.getID() |
568 |
|
569 |
class PropertySet(Primitive): |
570 |
""" |
571 |
defines a group L{Primitive} objects. |
572 |
""" |
573 |
def __init__(self,tag=None,*items): |
574 |
super(PropertySet, self).__init__() |
575 |
self.__items=items |
576 |
self.__tag=tag |
577 |
def getHistory(self): |
578 |
out=set([self, self.getBoundaryLoop().getHistory()]) |
579 |
for i in self.getHoles(): out|=i.getHistory() |
580 |
return out |
581 |
|
582 |
class PrimitiveStack(object): |
583 |
def __init__(self,*items): |
584 |
self.__prims=set() |
585 |
for i in items: |
586 |
self.__prims|=i.getHistory() |
587 |
self.__prims=list(self.__prims) |
588 |
self.__prims.sort() |
589 |
|
590 |
def getGmshCommands(self): |
591 |
out="" |
592 |
for i in self.__prims: |
593 |
out+=i.getGmshCommand()+"\n" |
594 |
return out |