/[escript]/trunk/modellib/py_src/geometry.py
ViewVC logotype

Annotation of /trunk/modellib/py_src/geometry.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1098 - (hide annotations)
Mon Apr 16 23:15:23 2007 UTC (12 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 20784 byte(s)
add a few #pragma ivdep which should speed up MV but cannot confiirm this on Pentium




1 jgs 127 # $Id$
2    
3 elspeth 628 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
4     http://www.access.edu.au
5     Primary Business: Queensland, Australia"""
6     __license__="""Licensed under the Open Software License version 3.0
7     http://www.opensource.org/licenses/osl-3.0.php"""
8 jgs 127
9 elspeth 628
10 jgs 149 from esys.escript import *
11     from esys.escript.modelframe import Model,ParameterSet
12     from esys import finley
13 jgs 127
14 gross 823 class FinleyReader(ParameterSet):
15 jgs 149 """
16 gross 819 reads finley mesh file.
17 jgs 147
18 gross 938 @ivar source: mesh file in finley or gmsh format
19     @type source: C{DataSource}
20 gross 823 @ivar intergrationOrder: integration order, default -1 (in).
21     @type intergrationOrder: C{int}
22 gross 938 @ivar reducedIntegrationOrder: reduced integration order, default -1 (in).
23     @type reducedIntegrationOrder: C{int}
24     @ivar optimizeLabeling: switches on optimization of the labeling of the nodes
25     @type optimizeLabeling: C{bool}
26 gross 398 """
27 gross 918 def __init__(self,**kwargs):
28 gross 823 """
29     initializes the object
30     """
31 gross 918 super(FinleyReader,self).__init__(**kwargs)
32 gross 823 self.declareParameter(source="none",
33 gross 997 dim=None,
34 gross 938 optimizeLabeling=True,
35     reducedIntegrationOrder=-1,
36     integrationOrder=-1)
37 gross 823 self.__domain=None
38 gross 398
39 gross 953
40 gross 823 def domain(self):
41     """
42     returns the domain
43    
44     @return: the domain
45     @rtype: L{Domain}
46     """
47 gross 911 if self.__domain == None:
48 gross 938 if self.source.fileformat == "fly":
49     self.__domain=finley.ReadMesh(self.source.getLocalFileName(),self.integrationOrder)
50     elif self.source.fileformat == "gmsh":
51 gross 997 if self.dim==None:
52     dim=3
53     else:
54     dim=self.dim
55     self.__domain=finley.ReadGmsh(self.source.getLocalFileName(),dim,self.integrationOrder,self.reducedIntegrationOrder, self.optimizeLabeling)
56 gross 938 else:
57     raise TypeError("unknown mesh file format %s."%self.source.fileformat)
58 gross 950 self.trace("mesh read from %s in %s format."%(self.source.getLocalFileName(), self.source.fileformat))
59 gross 823 return self.__domain
60     class RectangularDomain(ParameterSet):
61 gross 398 """
62     Generates a mesh over a rectangular domain finley.
63    
64 gross 819 @ivar dim: spatial dimension, default =2 (in).
65     @type dim: spatial dimension
66     @ivar l: spatial lengths, default [1.,1.,1.] (in).
67     @type l: C{list} of C{floats}s
68     @ivar n: number of elements, default [10,10,10] (in).
69     @type n: C{list} of C{int}s
70     @ivar order: element order, default 1 (in).
71     @type order: C{int}
72     @ivar periodic: flags for periodicity, default [False,False,False] (in).
73     @type periodic: C{list} of C{bool}s
74 gross 823 @ivar intergrationOrder: integration order, default -1 (in).
75     @type intergrationOrder: C{int}
76 jgs 147 """
77 gross 918 def __init__(self,**kwargs):
78 gross 819 """
79     initializes the object
80     """
81 gross 918 super(RectangularDomain,self).__init__(**kwargs)
82 jgs 147 self.declareParameter(dim=2,\
83 jgs 127 l=[1.,1.,1.],\
84     n=[10,10,10], \
85     order=1,\
86 gross 823 periodic=[False,False,False],
87     integrationOrder=-1)
88     self.__domain=None
89 jgs 127
90 gross 823 def domain(self):
91 gross 819 """
92 gross 823 returns the domain
93    
94     @return: the domain
95     @rtype: L{Domain}
96 gross 819 """
97 gross 911 if self.__domain==None:
98 gross 823 if self.dim==2:
99 gross 911 self.__domain=finley.Rectangle(n0=self.n[0],\
100 gross 823 n1=self.n[1],\
101     l0=self.l[0],\
102     l1=self.l[1],\
103     order=self.order, \
104     periodic0=self.periodic[0], \
105     periodic1=self.periodic[1], \
106     integrationOrder=self.integrationOrder)
107     else:
108     self.__domain=finley.Brick(n0=self.n[0],\
109     n1=self.n[1],\
110     n2=self.n[2],\
111     l0=self.l[0],\
112     l1=self.l[1],\
113     l2=self.l[2],\
114     order=self.order, \
115     periodic0=self.periodic[0], \
116     periodic1=self.periodic[1], \
117     periodic2=self.periodic[2], \
118     integrationOrder=self.integrationOrder)
119     return self.__domain
120 jgs 147
121 gross 823 class UpdateGeometry(Model):
122     """
123     applies a displacement field to a domain
124    
125     @ivar displacement: displacements applied to the original mesh coordinates (in).
126     @type displacement: L{escript.Vector}
127     @ivar domain: domain
128     @type domain: L{escript.Domain}
129     """
130 gross 918 def __init__(self,**kwargs):
131 gross 823 """
132     set-up the object
133     """
134 gross 918 super(UpdateGeometry, self).__init__(**kwargs)
135 gross 823 self.declareParameter(domain=None,\
136     displacement=None)
137 jgs 147
138 gross 823
139     def doInitialization(self):
140     """
141     initialize model
142     """
143     self.__x=self.domain.getX()
144     self.__reset=True
145    
146     def doStepPreprocessing(self,dt):
147     """
148     applies the current L{displacement} to mesh nodes if required.
149     """
150     if self.__reset:
151     self.trace("mesh nodes updated.")
152     self.domain.setX(self.__x+self.displacement)
153     self.__reset=False
154    
155     def doStep(self,dt):
156     """
157     applies the current L{displacement} to mesh nodes.
158     """
159     self.trace("mesh nodes updated.")
160     self.domain.setX(self.__x+self.displacement)
161     self.__reset=True
162    
163     def doStepPostprocessing(self,dt):
164     """
165     marks nodes as beeing updated.
166     """
167     self.__reset=False
168    
169 gross 953 class ConstrainerOverBox(Model):
170     """
171     Creates a characteristic function for the location of constraints
172     for all components of a value and selects the value from an initial value
173     ate these locations.
174    
175     In the case that the spatial dimension is two, the arguments front and back are ignored.
176    
177     @ivar domain: domain (in).
178     @ivar left: True to set a constraint at the left face of the domain (x[0]=min x[0]), default False (in).
179     @ivar right: True to set a constraint at the left face of the domain (x[0]=max x[0]), default False (in).
180     @ivar top: True to set a constraint at the left face of the domain (x[1]=min x[1]), default False (in).
181     @ivar bottom: True to set a constraint at the left face of the domain (x[1]=max x[1]), default False (in).
182     @ivar front: True to set a constraint at the left face of the domain (x[2]=min x[2]), default False (in).
183     @ivar back: True to set a constraint at the left face of the domain (x[2]=max x[2]), default False (in).
184     @ivar tol: absolute tolerance for "x=max x" condition, default 1.e-8 (in).
185     """
186     def __init__(self,**kwargs):
187     super(ConstrainerOverBox, self).__init__(**kwargs)
188     self.declareParameter(domain=None, \
189     value=None, \
190     left=False, \
191     right=False, \
192     top=False, \
193     bottom=False, \
194     front=False, \
195     back=False, \
196     tol=1.e-8)
197     self.__value_of_constraint = None
198     self.__location_of_constraint=None
199     def location_of_constraint(self):
200     """
201     return the values used to constrain a solution
202    
203     @return: the mask marking the locations of the constraints
204     @rtype: L{escript.Scalar}
205     """
206     if self.__location_of_constraint == None: self.__setOutput()
207     return self.__location_of_constraint
208    
209     def value_of_constraint(self):
210     """
211     return the values used to constrain a solution
212    
213     @return: values to be used at the locations of the constraints. If
214     L{value} is not given C{None} is rerturned.
215     @rtype: L{escript.Scalar}
216     """
217     if self.__location_of_constraint == None: self.__setOutput()
218     return self.__value_of_constraint
219    
220     def __setOutput(self):
221 gross 993 if self.__location_of_constraint == None:
222     x=self.domain.getX()
223     val=self.value
224     if isinstance(val, int) or isinstance(val, float):
225     shape=()
226     elif isinstance(val, list) or isinstance(val, tuple) :
227     shape=(len(val),)
228     elif isinstance(val, numarray.NumArray):
229     shape=val.shape
230     elif val == None:
231     shape=()
232     else:
233     shape=val.getShape()
234     self.__location_of_constraint=Data(0,shape,x.getFunctionSpace())
235     if self.domain.getDim()==3:
236     x0,x1,x2=x[0],x[1],x[2]
237     if self.left: self.__location_of_constraint+=whereZero(x0-inf(x0),self.tol)
238     if self.right: self.__location_of_constraint+=whereZero(x0-sup(x0),self.tol)
239     if self.front: self.__location_of_constraint+=whereZero(x1-inf(x1),self.tol)
240     if self.back: self.__location_of_constraint+=whereZero(x1-sup(x1),self.tol)
241     if self.bottom: self.__location_of_constraint+=whereZero(x2-inf(x2),self.tol)
242     if self.top: self.__location_of_constraint+=whereZero(x2-sup(x2),self.tol)
243     else:
244     x0,x1=x[0],x[1]
245     if self.left: self.__location_of_constraint+=whereZero(x0-inf(x0),self.tol)
246     if self.right: self.__location_of_constraint+=whereZero(x0-sup(x0),self.tol)
247     if self.bottom: self.__location_of_constraint+=whereZero(x1-inf(x1),self.tol)
248     if self.top: self.__location_of_constraint+=whereZero(x1-sup(x1),self.tol)
249     if not self.value == None:
250     self.__value_of_constraint=self.__location_of_constraint*self.value
251 gross 911 class ScalarConstrainerOverBox(Model):
252 gross 823 """
253     Creates a characteristic function for the location of constraints
254     for a scalar value and selects the value from an initial value
255     ate these locations.
256 jgs 127
257 gross 823 In the case that the spatial dimension is two, the arguments front and back are ignored.
258 jgs 127
259 gross 823 @ivar domain: domain (in).
260     @ivar left: True to set a constraint at the left face of the domain (x[0]=min x[0]), default False (in).
261     @ivar right: True to set a constraint at the left face of the domain (x[0]=max x[0]), default False (in).
262     @ivar top: True to set a constraint at the left face of the domain (x[1]=min x[1]), default False (in).
263     @ivar bottom: True to set a constraint at the left face of the domain (x[1]=max x[1]), default False (in).
264     @ivar front: True to set a constraint at the left face of the domain (x[2]=min x[2]), default False (in).
265     @ivar back: True to set a constraint at the left face of the domain (x[2]=max x[2]), default False (in).
266     @ivar tol: absolute tolerance for "x=max x" condition, default 1.e-8 (in).
267     """
268 gross 918 def __init__(self,**kwargs):
269     super(ScalarConstrainerOverBox, self).__init__(**kwargs)
270 jgs 127 self.declareParameter(domain=None, \
271 gross 823 value=None, \
272 jgs 147 left=False, \
273     right=False, \
274     top=False, \
275     bottom=False, \
276     front=False, \
277 gross 819 back=False, \
278 gross 823 tol=1.e-8)
279     self.__value_of_constraint = None
280     self.__location_of_constraint=None
281     def location_of_constraint(self):
282 jgs 149 """
283 gross 823 return the values used to constrain a solution
284    
285     @return: the mask marking the locations of the constraints
286     @rtype: L{escript.Scalar}
287 jgs 149 """
288 gross 908 if self.__location_of_constraint == None: self.__setOutput()
289 gross 823 return self.__location_of_constraint
290    
291     def value_of_constraint(self):
292     """
293     return the values used to constrain a solution
294    
295     @return: values to be used at the locations of the constraints. If
296     L{value} is not given C{None} is rerturned.
297     @rtype: L{escript.Scalar}
298     """
299 gross 908 if self.__location_of_constraint == None: self.__setOutput()
300 gross 823 return self.__value_of_constraint
301    
302     def __setOutput(self):
303 gross 819 x=self.domain.getX()
304 gross 823 self.__location_of_constraint=Scalar(0,x.getFunctionSpace())
305 gross 819 if self.domain.getDim()==3:
306     x0,x1,x2=x[0],x[1],x[2]
307 gross 823 if self.left: self.__location_of_constraint+=whereZero(x0-inf(x0),self.tol)
308     if self.right: self.__location_of_constraint+=whereZero(x0-sup(x0),self.tol)
309     if self.front: self.__location_of_constraint+=whereZero(x1-inf(x1),self.tol)
310     if self.back: self.__location_of_constraint+=whereZero(x1-sup(x1),self.tol)
311     if self.bottom: self.__location_of_constraint+=whereZero(x2-inf(x2),self.tol)
312     if self.top: self.__location_of_constraint+=whereZero(x2-sup(x2),self.tol)
313 gross 819 else:
314     x0,x1=x[0],x[1]
315 gross 823 if self.left: self.__location_of_constraint+=whereZero(x0-inf(x0),self.tol)
316     if self.right: self.__location_of_constraint+=whereZero(x0-sup(x0),self.tol)
317     if self.bottom: self.__location_of_constraint+=whereZero(x1-inf(x1),self.tol)
318     if self.top: self.__location_of_constraint+=whereZero(x1-sup(x1),self.tol)
319 gross 993 if not self.value == None:
320 gross 823 self.__value_of_constraint=self.__location_of_constraint*self.value
321 jgs 147
322 gross 911 class VectorConstrainerOverBox(Model):
323 jgs 149 """
324 gross 819 Creates a characteristic function for the location of constraints vector value.
325     In the case that the spatial dimension is two, the arguments front and
326     back as well as the third component of each argument is ignored.
327 jgs 127
328 gross 819 @ivar domain: domain
329     @ivar left: list of three boolean. left[i]==True sets a constraint for the i-th component at the left face of the domain (x[0]=min x[0]),
330     default [False,False,False] (in).
331     @ivar right: list of three boolean. left[i]==True sets a constraint for the i-th component at the right face of the domain (x[0]=max x[0]),
332     default [False,False,False] (in).
333     @ivar top: list of three boolean. left[i]==True sets a constraint for the i-th component at the top face of the domain (x[1]=min x[1]),
334     default [False,False,False] (in).
335     @ivar bottom: list of three boolean. left[i]==True sets a constraint for the i-th component at the bottom face of the domain (x[1]=min x[1]),
336     default [False,False,False] (in).
337     @ivar front: list of three boolean. left[i]==True sets a constraint for the i-th component at the front face of the domain (x[2]=min x[2]),
338     default [False,False,False] (in).
339     @ivar back: list of three boolean. left[i]==True sets a constraint for the i-th component at the back face of the domain (x[2]=max x[2]),
340     default [False,False,False] (in).
341     @ivar tol: absolute tolerance for "x=max x" condition, default 1.e-8 (in).
342 jgs 127 """
343 gross 918 def __init__(self, **kwargs):
344     super(VectorConstrainerOverBox, self).__init__(**kwargs)
345 jgs 127 self.declareParameter(domain=None, \
346 gross 823 value=None, \
347 jgs 147 left=[0,0,0], \
348     right=[0,0,0], \
349     top=[0,0,0], \
350     bottom=[0,0,0], \
351 gross 819 front=[0,0,0], \
352 gross 821 back=[0,0,0], \
353 gross 823 tol=1.e-8)
354     self.__value_of_constraint = None
355     self.__location_of_constraint=None
356    
357     def location_of_constraint(self):
358 jgs 149 """
359 gross 823 return the values used to constrain a solution
360    
361     @return: the mask marking the locations of the constraints
362     @rtype: L{escript.Vector}
363 jgs 149 """
364 gross 908 if self.__location_of_constraint == None: self.__setOutput()
365 gross 823 return self.__location_of_constraint
366    
367     def value_of_constraint(self):
368     """
369     return the values used to constrain a solution
370 jgs 149
371 gross 823 @return: values to be used at the locations of the constraints. If
372     L{value} is not given C{None} is rerturned.
373     @rtype: L{escript.Vector}
374     """
375 gross 908 if self.__location_of_constraint == None: self.__setOutput()
376 gross 823 return self.__value_of_constraint
377    
378     def __setOutput(self):
379     x=self.domain.getX()
380     self.__location_of_constraint=Vector(0,x.getFunctionSpace())
381     if self.domain.getDim()==3:
382     x0,x1,x2=x[0],x[1],x[2]
383     left_mask=whereZero(x0-inf(x0),self.tol)
384     if self.left[0]: self.__location_of_constraint+=left_mask*[1.,0.,0.]
385     if self.left[1]: self.__location_of_constraint+=left_mask*[0.,1.,0.]
386     if self.left[2]: self.__location_of_constraint+=left_mask*[0.,0.,1.]
387     right_mask=whereZero(x0-sup(x0),self.tol)
388     if self.right[0]: self.__location_of_constraint+=right_mask*[1.,0.,0.]
389     if self.right[1]: self.__location_of_constraint+=right_mask*[0.,1.,0.]
390     if self.right[2]: self.__location_of_constraint+=right_mask*[0.,0.,1.]
391     front_mask=whereZero(x1-inf(x1),self.tol)
392     if self.front[0]: self.__location_of_constraint+=front_mask*[1.,0.,0.]
393     if self.front[1]: self.__location_of_constraint+=front_mask*[0.,1.,0.]
394     if self.front[2]: self.__location_of_constraint+=front_mask*[0.,0.,1.]
395     back_mask=whereZero(x1-sup(x1),self.tol)
396     if self.back[0]: self.__location_of_constraint+=back_mask*[1.,0.,0.]
397     if self.back[1]: self.__location_of_constraint+=back_mask*[0.,1.,0.]
398     if self.back[2]: self.__location_of_constraint+=back_mask*[0.,0.,1.]
399     bottom_mask=whereZero(x2-inf(x2),self.tol)
400     if self.bottom[0]: self.__location_of_constraint+=bottom_mask*[1.,0.,0.]
401     if self.bottom[1]: self.__location_of_constraint+=bottom_mask*[0.,1.,0.]
402     if self.bottom[2]: self.__location_of_constraint+=bottom_mask*[0.,0.,1.]
403     top_mask=whereZero(x2-sup(x2),self.tol)
404     if self.top[0]: self.__location_of_constraint+=top_mask*[1.,0.,0.]
405     if self.top[1]: self.__location_of_constraint+=top_mask*[0.,1.,0.]
406     if self.top[2]: self.__location_of_constraint+=top_mask*[0.,0.,1.]
407 gross 993 if not self.value == None:
408 gross 823 self.__value_of_constraint=self.__location_of_constraint*self.value
409     else:
410     x0,x1=x[0],x[1]
411     left_mask=whereZero(x0-inf(x0),self.tol)
412     if self.left[0]: self.__location_of_constraint+=left_mask*[1.,0.]
413     if self.left[1]: self.__location_of_constraint+=left_mask*[0.,1.]
414     right_mask=whereZero(x0-sup(x0),self.tol)
415     if self.right[0]: self.__location_of_constraint+=right_mask*[1.,0.]
416     if self.right[1]: self.__location_of_constraint+=right_mask*[0.,1.]
417     bottom_mask=whereZero(x1-inf(x1),self.tol)
418     if self.bottom[0]: self.__location_of_constraint+=bottom_mask*[1.,0.]
419     if self.bottom[1]: self.__location_of_constraint+=bottom_mask*[0.,1.]
420     top_mask=whereZero(x1-sup(x1),self.tol)
421     if self.top[0]: self.__location_of_constraint+=top_mask*[1.,0.]
422     if self.top[1]: self.__location_of_constraint+=top_mask*[0.,1.]
423 gross 993 if not self.value == None:
424 gross 823 self.__value_of_constraint=self.__location_of_constraint*self.value[:2]
425    
426 jgs 149 # vim: expandtab shiftwidth=4:

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26