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