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

Contents of /trunk/pycad/py_src/primitives.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 917 - (show annotations)
Tue Jan 2 02:46:53 2007 UTC (12 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 24138 byte(s)
some EsysXML input files. rebuild still fails for these files
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, Translation, Dilation, Transformation
30
31
32 def resetGlobalPrimitiveIdCounter():
33 """
34 initializes the global primitive ID counter
35 """
36 global global_primitive_id_counter
37 global_primitive_id_counter=1
38
39 def setToleranceForColocation(tol=1.e-11):
40 """
41 set the global tolerance for colocation checks to tol
42 """
43 global global_tolerance_for_colocation
44 global_tolerance_for_colocation=tol
45
46 def getToleranceForColocation():
47 """
48 returns the global tolerance for colocation checks
49 """
50 return global_tolerance_for_colocation
51
52 resetGlobalPrimitiveIdCounter()
53 setToleranceForColocation()
54
55 class Primitive(object):
56 """
57 template for elementary geometrical object
58 """
59 def __init__(self):
60 """
61
62 """
63 global global_primitive_id_counter
64 self.__ID=global_primitive_id_counter
65 global_primitive_id_counter+=1
66
67 def getID(self):
68 return self.__ID
69
70 def __repr__(self):
71 return "%s(%s)"%(self.__class__.__name__,self.getID())
72 def __cmp__(self,other):
73 return cmp(self.getID(),other.getID())
74
75 def getConstructionPoints(self):
76 """
77 returns the points used to construct the primitive
78 """
79 out=set()
80 for i in self.getPrimitives():
81 if isinstance(i,Point): out.add(i)
82 return list(out)
83
84 def getPrimitives(self):
85 """
86 returns primitives used to construct the primitive
87 """
88 return []
89
90 def copy(self):
91 """
92 returns a deep copy of the object
93 """
94 return Primitive()
95
96
97 def modifyBy(self,transformation):
98 """
99 modifies the coordinates by applying a transformation
100 """
101 for p in self.getConstructionPoints(): p.modifyBy(transformation)
102
103
104 def __add__(self,other):
105 """
106 returns a new object shifted by other
107 """
108 return self.apply(Translation(numarray.array(other,_TYPE)))
109
110 def __sub__(self,other):
111 """
112 returns a new object shifted by other
113 """
114 return self.apply(Translation(-numarray.array(other,_TYPE)))
115
116 def __iadd__(self,other):
117 """
118 shifts the point by other
119 """
120 self.modifyBy(Translation(numarray.array(other,_TYPE)))
121 return self
122
123 def __isub__(self,other):
124 """
125 shifts the point by -other
126 """
127 self.modifyBy(Translation(-numarray.array(other,_TYPE)))
128 return self
129
130 def __imul__(self,other):
131 """
132 modifies object by applying L{Transformation} other. If other is not a L{Transformation} it will try convert it.
133 """
134 if isinstance(other,int) or isinstance(other,float):
135 trafo=Dilation(other)
136 elif isinstance(other,numarray.NumArray):
137 trafo=Translation(other)
138 elif isinstance(other,Transformation):
139 trafo=other
140 else:
141 raise TypeError, "cannot convert argument to Trnsformation class object."
142 self.modifyBy(trafo)
143 return self
144
145 def __rmul__(self,other):
146 """
147 applies L{Transformation} other to object. If other is not a L{Transformation} it will try convert it.
148 """
149 if isinstance(other,int) or isinstance(other,float):
150 trafo=Dilation(other)
151 elif isinstance(other,numarray.NumArray):
152 trafo=Translation(other)
153 elif isinstance(other,Transformation):
154 trafo=other
155 else:
156 raise TypeError, "cannot convert argument to Transformation class object."
157 return self.apply(trafo)
158
159 def __neg__(self):
160 return ReversedPrimitive(self)
161
162 def setLocalScale(self,factor=1.):
163 """
164 sets the local refinement factor
165 """
166 for p in self.getConstructionPoints(): p.setLocalScale(factor)
167
168 def getGmshCommand(self, local_scaling_factor=1.):
169 """
170 returns the Gmsh command(s) to create the primitive
171 """
172 raise NotImplementedError("getGmshCommand is not implemented for this class %s."%self.__class__.__name__)
173
174 def apply(self,transformation):
175 """
176 returns a new object by applying the transformation
177 """
178 raise NotImplementedError("apply is not implemented for this class %s."%self.__class__.__name__)
179
180 def isColocated(self,primitive):
181 """
182 returns True is the two primitives are located at the smae position
183 """
184 raise NotImplementedError("isColocated is not implemented for this class %s."%self.__class__.__name__)
185
186 class Point(Primitive):
187 """
188 a three dimensional point
189 """
190 def __init__(self,x=0.,y=0.,z=0.,local_scale=1.):
191 """
192 creates a point with coorinates x,y,z with the local refinement factor local_scale
193 """
194 super(Point, self).__init__()
195 self.setCoordinates(numarray.array([x,y,z],_TYPE))
196 self.setLocalScale(local_scale)
197
198 def setLocalScale(self,factor=1.):
199 """
200 sets the local refinement factor
201 """
202 if factor<=0.:
203 raise ValueError("scaling factor must be positive.")
204 self.__local_scale=factor
205 def getLocalScale(self):
206 """
207 returns the local refinement factor
208 """
209 return self.__local_scale
210 def getCoordinates(self):
211 """
212 returns the coodinates of the point as L{numarray.NumArray} object
213 """
214 return self._x
215 def setCoordinates(self,x):
216 """
217 returns the coodinates of the point as L{numarray.NumArray} object
218 """
219 if not isinstance(x, numarray.NumArray):
220 self._x=numarray.array(x,_TYPE)
221 else:
222 self._x=x
223
224 def getPrimitives(self):
225 """
226 returns primitives used to construct the primitive
227 """
228 return [self]
229
230 def isColocated(self,primitive):
231 """
232 returns True if L{Point} primitive is colocation (same coordinates)
233 that means if |self-primitive| <= tol * max(|self|,|primitive|)
234 """
235 if isinstance(primitive,Point):
236 primitive=primitive.getCoordinates()
237 c=self.getCoordinates()
238 d=c-primitive
239 return numarray.dot(d,d)<=getToleranceForColocation()**2*max(numarray.dot(c,c),numarray.dot(primitive,primitive))
240 else:
241 return False
242
243 def copy(self):
244 """
245 returns a deep copy of the point
246 """
247 c=self.getCoordinates()
248 return Point(c[0],c[1],c[2],local_scale=self.getLocalScale())
249
250 def modifyBy(self,transformation):
251 """
252 modifies the coordinates by applying a transformation
253 """
254 self.setCoordinates(transformation(self.getCoordinates()))
255
256 def apply(self,transformation):
257 """
258 returns a new L{Point} by applying the transformation
259 """
260 new_p=self.copy()
261 new_p.modifyBy(transformation)
262 return new_p
263
264 def getGmshCommand(self, local_scaling_factor=1.):
265 """
266 returns the Gmsh command(s) to create the primitive
267 """
268 c=self.getCoordinates()
269 return "Point(%s) = {%s , %s, %s , %s };"%(self.getID(),c[0],c[1],c[2], self.getLocalScale()*local_scaling_factor)
270
271 class Primitive1D(Primitive):
272 """
273 general one-dimensional primitive
274 """
275 def __init__(self,*args):
276 """
277 create a one-dimensional primitive
278 """
279 super(Primitive1D, self).__init__()
280
281 class Curve(Primitive1D):
282 """
283 a curve defined through a list of control points.
284 """
285 def __init__(self,*points):
286 """
287 defines a curve form control points
288 """
289 if len(points)<2:
290 raise TypeError("Curve needs at least two points")
291 super(Curve, self).__init__()
292 i=0
293 for p in points:
294 i+=1
295 if not isinstance(p,Point): raise TypeError("%s-th argument is not a Point object."%i)
296 self.__points=points
297
298 def __len__(self):
299 """
300 returns the number of control points
301 """
302 return len(self.__points)
303
304 def getStartPoint(self):
305 """
306 returns start point
307 """
308 return self.__points[0]
309
310 def getEndPoint(self):
311 """
312 returns end point
313 """
314 return self.__points[-1]
315
316 def getControlPoints(self):
317 """
318 returns a list of the points
319 """
320 return self.__points
321
322 def getPrimitives(self):
323 """
324 returns primitives used to construct the Curve
325 """
326 out=set()
327 for p in self.getControlPoints(): out|=set(p.getPrimitives())
328 out.add(self)
329 return list(out)
330
331 def copy(self):
332 """
333 returns a deep copy
334 """
335 new_p=[]
336 for p in self.getControlPoints(): new_p.append(p.copy())
337 return self.__class__(*tuple(new_p))
338
339
340 def apply(self,transformation):
341 """
342 applies transformation
343 """
344 new_p=[]
345 for p in self.getControlPoints(): new_p.append(p.apply(transformation))
346 return self.__class__(*tuple(new_p))
347
348 def isColocated(self,primitive):
349 """
350 returns True curves are on the same position
351 """
352 if isinstance(primitive,self.__class__):
353 if len(primitive) == len(self):
354 cp0=self.getControlPoints()
355 cp1=primitive.getControlPoints()
356 for i in range(len(cp0)):
357 if not cp0[i].isColocated(cp1[i]):
358 return False
359 return True
360 else:
361 return False
362 else:
363 return False
364
365 class Spline(Curve):
366 """
367 a spline curve defined through a list of control points.
368 """
369 def getGmshCommand(self):
370 """
371 returns the Gmsh command(s) to create the Curve
372 """
373 out=""
374 for i in self.getControlPoints():
375 if len(out)>0:
376 out+=", %s"%i.getID()
377 else:
378 out="%s"%i.getID()
379 return "Spline(%s) = {%s};"%(self.getID(),out)
380
381
382 class BezierCurve(Curve):
383 """
384 a Bezier curve
385 """
386 def getGmshCommand(self):
387 """
388 returns the Gmsh command(s) to create the Curve
389 """
390 out=""
391 for i in self.getControlPoints():
392 if len(out)>0:
393 out+=", %s"%i.getID()
394 else:
395 out="%s"%i.getID()
396 return "Bezier(%s) = {%s};"%(self.getID(),out)
397
398 class BSpline(Curve):
399 """
400 a BSpline curve. Control points may be repeated.
401 """
402 def getGmshCommand(self):
403 """
404 returns the Gmsh command(s) to create the Curve
405 """
406 out=""
407 for i in self.getControlPoints():
408 if len(out)>0:
409 out+=", %s"%i.getID()
410 else:
411 out="%s"%i.getID()
412 return "BSpline(%s) = {%s};"%(self.getID(),out)
413
414 class Line(Curve):
415 """
416 a line is defined by two L{Point}s
417 """
418 def __init__(self,*points):
419 """
420 defines a line with start and end point
421 """
422 if len(points)!=2:
423 raise TypeError("Line needs two points")
424 super(Line, self).__init__(*points)
425 def getGmshCommand(self):
426 """
427 returns the Gmsh command(s) to create the Curve
428 """
429 return "Line(%s) = {%s, %s};"%(self.getID(),self.getStartPoint().getID(),self.getEndPoint().getID())
430
431
432 class Arc(Primitive1D):
433 """
434 defines an arc which is strictly, smaller than Pi
435 """
436 def __init__(self,center,start,end):
437 """
438 creates an arc by the start point, end point and center
439 """
440 if not isinstance(center,Point): raise TypeError("center needs to be a Point object.")
441 if not isinstance(end,Point): raise TypeError("end needs to be a Point object.")
442 if not isinstance(start,Point): raise TypeError("start needs to be a Point object.")
443 # TODO: check length of circle.
444 super(Arc, self).__init__()
445 self.__center=center
446 self.__start=start
447 self.__end=end
448
449 def getStartPoint(self):
450 """
451 returns start point
452 """
453 return self.__start
454
455 def getEndPoint(self):
456 """
457 returns end point
458 """
459 return self.__end
460
461 def getCenterPoint(self):
462 """
463 returns center
464 """
465 return self.__center
466
467 def getPrimitives(self):
468 """
469 returns the primitives used to construct the Curve
470 """
471 out=set()
472 out|=set(self.getStartPoint().getPrimitives())
473 out|=set(self.getEndPoint().getPrimitives())
474 out|=set(self.getCenterPoint().getPrimitives())
475 out.add(self)
476 return list(out)
477
478 def getGmshCommand(self):
479 """
480 returns the Gmsh command(s) to create the primitive
481 """
482 return "Circle(%s) = {%s, %s, %s};"%(self.getID(),self.getStartPoint().getID(),self.getCenterPoint().getID(),self.getEndPoint().getID())
483
484 def copy(self):
485 """
486 returns a deep copy
487 """
488 return Arc(self.getCenterPoint().copy(),self.getStartPoint().copy(),self.getEndPoint().copy())
489
490 def apply(self,transformation):
491 """
492 applies transformation
493 """
494 return Arc(self.getCenterPoint().apply(transformation),self.getStartPoint().apply(transformation),self.getEndPoint().apply(transformation))
495
496 def isColocated(self,primitive):
497 """
498 returns True curves are on the same position
499 """
500 if isinstance(primitive,Arc):
501 if not self.getCenterPoint().isColocated(primitive.getCenterPoint()): return False
502 if not self.getEndPoint().isColocated(primitive.getEndPoint()): return False
503 if not self.getStartPoint().isColocated(primitive.getStartPoint()): return False
504 return True
505 else:
506 return False
507
508 #=================================================================================================================================
509 class CurveLoop(Primitive):
510 """
511 An oriented loop of curves.
512
513 The loop must be closed and the L{Curves}s should be oriented consistently.
514 """
515 def __init__(self,*curves):
516 """
517 creates a polygon from a list of line curves. The curves must form a closed loop.
518 """
519 super(CurveLoop, self).__init__()
520 self.__curves=[]
521 self.addCurve(*curves)
522 def addCurve(self,*curves):
523 for i in range(len(curves)):
524 if not curves[i].isCurve():
525 raise TypeError("%s-th argument is not a Curve object."%i)
526 self.__curves+=curves
527
528 def isCurveLoop(self):
529 return True
530 def getCurves(self):
531 return self.__curves
532 def __add__(self,other):
533 return CurveLoop(*tuple([c+other for c in self.getCurves()[::-1]]))
534 def __len__(self):
535 return len(self.__curves)
536 def getPrimitives(self):
537 out=set([self])
538 for i in self.getCurves(): out|=i.getPrimitives()
539 return out
540 def getConstructionPoints(self):
541 out=set()
542 for i in self.getCurves(): out|=i.getConstructionPoints()
543 return out
544 def getGmshCommand(self):
545 out=""
546 for i in self.getCurves():
547 if len(out)>0:
548 out+=", %s"%i.getID()
549 else:
550 out="%s"%i.getID()
551 return "Line Loop(%s) = {%s};"%(self.getID(),out)
552
553 class Surface(Primitive):
554 """
555 a surface
556 """
557 def __init__(self,loop):
558 """
559 creates a surface with boundary loop
560
561 @param loop: L{CurveLoop} defining the boundary of the surface
562 """
563 super(Surface, self).__init__()
564 if not loop.isCurveLoop():
565 raise TypeError("argument loop needs to be a CurveLoop object.")
566 self.__loop=loop
567 def isSurface(self):
568 return True
569 def getBoundaryLoop(self):
570 return self.__loop
571 def __add__(self,other):
572 return Surface(self.getBoundaryLoop()+other)
573 def getPrimitives(self):
574 out=set([self]) | self.getBoundaryLoop().getPrimitives()
575 return out
576 def getConstructionPoints(self):
577 return self.getBoundaryLoop().getConstructionPoints()
578 def getGmshCommand(self):
579 return "Ruled Surface(%s) = {%s};"%(self.getID(),self.getBoundaryLoop().getID())
580
581 class PlaneSurface(Surface):
582 """
583 a plane surface with holes
584 """
585 def __init__(self,loop,holes=[]):
586 """
587 creates a plane surface.
588
589 @param loop: L{CurveLoop} defining the boundary of the surface
590 @param holes: list of L{CurveLoop} defining holes in the surface.
591 @note: A CurveLoop defining a hole should not have any lines in common with the exterior CurveLoop.
592 A CurveLoop defining a hole should not have any lines in common with another CurveLoop defining a hole in the same surface.
593 """
594 super(PlaneSurface, self).__init__(loop)
595 for i in range(len(holes)):
596 if not holes[i].inCurveLoop():
597 raise TypeError("%i th hole needs to be a CurveLoop object.")
598 self.__holes=holes
599 def getHoles(self):
600 return self.__holes
601 def __add__(self,other):
602 return PlaneSurface(self.getBoundaryLoop()+other, holes=[h+other for h in self.getHoles()])
603 def getPrimitives(self):
604 out=set([self]) | self.getBoundaryLoop().getPrimitives()
605 for i in self.getHoles(): out|=i.getPrimitives()
606 return out
607 def getConstructionPoints(self):
608 out=self.getBoundaryLoop().getConstructionPoints()
609 for i in self.getHoles(): out|=i.getConstructionPoints()
610 return out
611 def getGmshCommand(self):
612 out=""
613 for i in self.getHoles():
614 if len(out)>0:
615 out+=", %s"%i.getID()
616 else:
617 out="%s"%i.getID()
618 if len(out)>0:
619 return "Plane Surface(%s) = {%s, %s};"%(self.getID(),self.getBoundaryLoop().getID(), out)
620 else:
621 return "Plane Surface(%s) = {%s};"%(self.getID(),self.getBoundaryLoop().getID())
622
623 class RuledSurface(Surface):
624 """
625 A ruled surface, i.e., a surface that can be interpolated using transfinite interpolation
626 """
627 def __init__(self,loop):
628 """
629 creates a ruled surface from a
630
631 @param loop: L{CurveLoop} defining the boundary of the surface. There is a restriction of composed of either three or four L{Curve} objects.
632 """
633 if not loop.isCurveLoop():
634 raise TypeError("argument loop needs to be a CurveLoop object.")
635 if len(loop)<3:
636 raise TypeError("the loop must contain at least three Curves.")
637 super(RuledSurface, self).__init__(loop)
638 def __add__(self,other):
639 return RuledSurface(self.getBoundaryLoop()+other)
640 def getGmshCommand(self):
641 return "Ruled Surface(%s) = {%s};"%(self.getID(),self.getBoundaryLoop().getID())
642
643 class SurfaceLoop(Primitive):
644 """
645 a surface loop. It defines the shell of a volume.
646
647 The loop must represent a closed shell, and the L{Surface}s should be oriented consistently.
648 """
649 def __init__(self,*surfaces):
650 """
651 creates a surface loop
652 """
653 super(SurfaceLoop, self).__init__()
654 self.__surfaces=[]
655 self.addSurface(*surfaces)
656 def addSurface(self,*surfaces):
657 for i in range(len(surfaces)):
658 if not surfaces[i].isSurface():
659 raise TypeError("%s-th argument is not a Surface object."%i)
660 self.__surfaces+=surfaces
661
662 def isSurfaceLoop(self):
663 return True
664 def getSurfaces(self):
665 return self.__surfaces
666 def __add__(self,other):
667 return SurfaceLoop([c+other for c in self.getSurfaces])
668 def __len__(self):
669 return len(self.__surfaces)
670 def getPrimitives(self):
671 out=set([self])
672 for i in self.getSurfaces(): out|=i.getPrimitives()
673 return out
674 def getConstructionPoints(self):
675 out=set()
676 for i in self.getSurfaces(): out|=i.getConstructionPoints()
677 return out
678 def getGmshCommand(self):
679 out=""
680 for i in self.getSurfaces():
681 if len(out)>0:
682 out+=", %s"%i.getID()
683 else:
684 out="%s"%i.getID()
685 return "Surface Loop(%s) = {%s};"%(self.getID(),out)
686
687 class Volume(Primitive):
688 """
689 a volume with holes.
690 """
691 def __init__(self,loop,holes=[]):
692 """
693 creates a volume
694
695 @param loop: L{SurfaceLoop} defining the boundary of the surface
696 @param holes: list of L{SurfaceLoop} defining holes in the surface.
697 @note: A SurfaceLoop defining a hole should not have any surfaces in common with the exterior SurfaceLoop.
698 A SurfaceLoop defining a hole should not have any surfaces in common with another SurfaceLoop defining a hole in the same volume.
699 """
700 super(Volume, self).__init__()
701 if not loop.isSurfaceLoop():
702 raise TypeError("argument loop needs to be a SurfaceLoop object.")
703 for i in range(len(holes)):
704 if not holes[i].isSurfaceLoop():
705 raise TypeError("%i th hole needs to be a SurfaceLoop object.")
706 self.__loop=loop
707 self.__holes=holes
708 def getHoles(self):
709 return self.__holes
710 def getSurfaceLoop(self):
711 return self.__loop
712 def __add__(self,other):
713 return Volume(self.getSurfaceLoop()+other, holes=[h+other for h in self.getHoles()])
714 def getPrimitives(self):
715 out=set([self]) | self.getSurfaceLoop().getPrimitives()
716 for i in self.getHoles(): out|=i.getPrimitives()
717 return out
718 def getConstructionPoints(self):
719 out=self.getSurfaceLoop().getConstructionPoints()
720 for i in self.getHoles(): out|=i.Points()
721 return out
722 def getGmshCommand(self):
723 out=""
724 for i in self.getHoles():
725 if len(out)>0:
726 out+=", %s"%i.getID()
727 else:
728 out="%s"%i.getID()
729 if len(out)>0:
730 return "Volume(%s) = {%s, %s};"%(self.getID(),self.getSurfaceLoop().getID(), out)
731 else:
732 return "Volume(%s) = {%s};"%(self.getID(),self.getSurfaceLoop().getID())
733
734 class ReversedPrimitive(object):
735 def __init__(self,prim):
736 self.__prim=prim
737 def __getattr__(self,name):
738 if name == "getID":
739 return self.getReverseID
740 else:
741 return getattr(self.__prim,name)
742 def getReverseID(self):
743 return -self.__prim.getID()
744
745 class PropertySet(Primitive):
746 """
747 defines a group L{Primitive} objects.
748 """
749 def __init__(self,tag=None,*items):
750 super(PropertySet, self).__init__()
751 self.__items=items
752 self.__tag=tag
753 def getPrimitives(self):
754 out=set([self, self.getBoundaryLoop().getPrimitives()])
755 for i in self.getHoles(): out|=i.getPrimitives()
756 return out
757
758 class PrimitiveStack(object):
759 def __init__(self,*items):
760 self.__prims=set()
761 for i in items:
762 self.__prims|=i.getPrimitives()
763 self.__prims=list(self.__prims)
764 self.__prims.sort()
765
766 def getGmshCommands(self):
767 out=""
768 for i in self.__prims:
769 out+=i.getGmshCommand()+"\n"
770 return out

  ViewVC Help
Powered by ViewVC 1.1.26