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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1098 - (show 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 # $Id$
2
3 __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
9
10 from esys.escript import *
11 from esys.escript.modelframe import Model,ParameterSet
12 from esys import finley
13
14 class FinleyReader(ParameterSet):
15 """
16 reads finley mesh file.
17
18 @ivar source: mesh file in finley or gmsh format
19 @type source: C{DataSource}
20 @ivar intergrationOrder: integration order, default -1 (in).
21 @type intergrationOrder: C{int}
22 @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 """
27 def __init__(self,**kwargs):
28 """
29 initializes the object
30 """
31 super(FinleyReader,self).__init__(**kwargs)
32 self.declareParameter(source="none",
33 dim=None,
34 optimizeLabeling=True,
35 reducedIntegrationOrder=-1,
36 integrationOrder=-1)
37 self.__domain=None
38
39
40 def domain(self):
41 """
42 returns the domain
43
44 @return: the domain
45 @rtype: L{Domain}
46 """
47 if self.__domain == None:
48 if self.source.fileformat == "fly":
49 self.__domain=finley.ReadMesh(self.source.getLocalFileName(),self.integrationOrder)
50 elif self.source.fileformat == "gmsh":
51 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 else:
57 raise TypeError("unknown mesh file format %s."%self.source.fileformat)
58 self.trace("mesh read from %s in %s format."%(self.source.getLocalFileName(), self.source.fileformat))
59 return self.__domain
60 class RectangularDomain(ParameterSet):
61 """
62 Generates a mesh over a rectangular domain finley.
63
64 @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 @ivar intergrationOrder: integration order, default -1 (in).
75 @type intergrationOrder: C{int}
76 """
77 def __init__(self,**kwargs):
78 """
79 initializes the object
80 """
81 super(RectangularDomain,self).__init__(**kwargs)
82 self.declareParameter(dim=2,\
83 l=[1.,1.,1.],\
84 n=[10,10,10], \
85 order=1,\
86 periodic=[False,False,False],
87 integrationOrder=-1)
88 self.__domain=None
89
90 def domain(self):
91 """
92 returns the domain
93
94 @return: the domain
95 @rtype: L{Domain}
96 """
97 if self.__domain==None:
98 if self.dim==2:
99 self.__domain=finley.Rectangle(n0=self.n[0],\
100 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
121 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 def __init__(self,**kwargs):
131 """
132 set-up the object
133 """
134 super(UpdateGeometry, self).__init__(**kwargs)
135 self.declareParameter(domain=None,\
136 displacement=None)
137
138
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 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 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 class ScalarConstrainerOverBox(Model):
252 """
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
257 In the case that the spatial dimension is two, the arguments front and back are ignored.
258
259 @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 def __init__(self,**kwargs):
269 super(ScalarConstrainerOverBox, self).__init__(**kwargs)
270 self.declareParameter(domain=None, \
271 value=None, \
272 left=False, \
273 right=False, \
274 top=False, \
275 bottom=False, \
276 front=False, \
277 back=False, \
278 tol=1.e-8)
279 self.__value_of_constraint = None
280 self.__location_of_constraint=None
281 def location_of_constraint(self):
282 """
283 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 """
288 if self.__location_of_constraint == None: self.__setOutput()
289 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 if self.__location_of_constraint == None: self.__setOutput()
300 return self.__value_of_constraint
301
302 def __setOutput(self):
303 x=self.domain.getX()
304 self.__location_of_constraint=Scalar(0,x.getFunctionSpace())
305 if self.domain.getDim()==3:
306 x0,x1,x2=x[0],x[1],x[2]
307 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 else:
314 x0,x1=x[0],x[1]
315 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 if not self.value == None:
320 self.__value_of_constraint=self.__location_of_constraint*self.value
321
322 class VectorConstrainerOverBox(Model):
323 """
324 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
328 @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 """
343 def __init__(self, **kwargs):
344 super(VectorConstrainerOverBox, self).__init__(**kwargs)
345 self.declareParameter(domain=None, \
346 value=None, \
347 left=[0,0,0], \
348 right=[0,0,0], \
349 top=[0,0,0], \
350 bottom=[0,0,0], \
351 front=[0,0,0], \
352 back=[0,0,0], \
353 tol=1.e-8)
354 self.__value_of_constraint = None
355 self.__location_of_constraint=None
356
357 def location_of_constraint(self):
358 """
359 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 """
364 if self.__location_of_constraint == None: self.__setOutput()
365 return self.__location_of_constraint
366
367 def value_of_constraint(self):
368 """
369 return the values used to constrain a solution
370
371 @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 if self.__location_of_constraint == None: self.__setOutput()
376 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 if not self.value == None:
408 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 if not self.value == None:
424 self.__value_of_constraint=self.__location_of_constraint*self.value[:2]
425
426 # 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