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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 953 by gross, Fri Feb 9 08:55:28 2007 UTC revision 1809 by ksteube, Thu Sep 25 06:43:44 2008 UTC
# Line 1  Line 1 
 # $Id$  
1    
2  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  ########################################################
3                      http://www.access.edu.au  #
4                  Primary Business: Queensland, Australia"""  # Copyright (c) 2003-2008 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7    #
8    # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11    #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
19               http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22  from esys.escript import *  from esys.escript import *
23  from esys.escript.modelframe import Model,ParameterSet  from esys.escript.modelframe import Model,ParameterSet
 from esys.pycad import TagMap  
24  from esys import finley  from esys import finley
25    
26  class FinleyReader(ParameterSet):  class FinleyReader(ParameterSet):
# Line 31  class FinleyReader(ParameterSet): Line 42  class FinleyReader(ParameterSet):
42            """            """
43            super(FinleyReader,self).__init__(**kwargs)            super(FinleyReader,self).__init__(**kwargs)
44            self.declareParameter(source="none",            self.declareParameter(source="none",
45                                  region_tag_map_source=None,                                  dim=None,
                                 surface_tag_map_source=None,  
46                                  optimizeLabeling=True,                                  optimizeLabeling=True,
47                                  reducedIntegrationOrder=-1,                                  reducedIntegrationOrder=-1,
48                                  integrationOrder=-1)                                  integrationOrder=-1)
49            self.__domain=None            self.__domain=None
           self.__region_tag_map=None  
           self.__surface_tag_map=None  
50    
51    
52         def domain(self):         def domain(self):
# Line 52  class FinleyReader(ParameterSet): Line 60  class FinleyReader(ParameterSet):
60               if  self.source.fileformat == "fly":               if  self.source.fileformat == "fly":
61                  self.__domain=finley.ReadMesh(self.source.getLocalFileName(),self.integrationOrder)                  self.__domain=finley.ReadMesh(self.source.getLocalFileName(),self.integrationOrder)
62               elif self.source.fileformat == "gmsh":               elif self.source.fileformat == "gmsh":
63                  self.__domain=finley.ReadGmsh(self.source.getLocalFileName(),self.integrationOrder,self.reducedIntegrationOrder, self.optimizeLabeling)                  if self.dim==None:
64                       dim=3
65                    else:
66                       dim=self.dim
67                    self.__domain=finley.ReadGmsh(self.source.getLocalFileName(),dim,self.integrationOrder,self.reducedIntegrationOrder, self.optimizeLabeling)
68               else:               else:
69                  raise TypeError("unknown mesh file format %s."%self.source.fileformat)                  raise TypeError("unknown mesh file format %s."%self.source.fileformat)
70               self.trace("mesh read from %s in %s format."%(self.source.getLocalFileName(), self.source.fileformat))                         self.trace("mesh read from %s in %s format."%(self.source.getLocalFileName(), self.source.fileformat))          
71            return self.__domain            return self.__domain
   
        def region_tag_map(self):  
           """  
           returns the map from regional tag names to tag integers used in the mesh  
   
           @return: the tag map  
           @rtype: L{TagMap}  
           """  
           if self.__region_tag_map == None:  
                self.__region_tag_map = TagMap()  
                if  self.region_tag_map_source != None:  
                    self.__region_tag_map.fillFromXML(open(self.region_tag_map_source.getLocalFileName()))  
                self.trace("region tag map read from %s in %s format."%(self.region_tag_map_source.getLocalFileName(), self.region_tag_map_source.fileformat))            
           return self.__region_tag_map  
   
        def surface_tag_map(self):  
           """  
           returns the map from surface tag names to tag integers used in the mesh  
   
           @return: the tag map  
           @rtype: L{TagMap}  
           """  
           if self.__surface_tag_map == None:  
                self.__surface_tag_map = TagMap()  
                if  self.surface_tag_map_source != None:  
                    self.__surface_tag_map.fillFromXML(open(self.surface_tag_map_source.getLocalFileName()))  
                self.trace("surface tag map read from %s in %s format."%(self.surface_tag_map_source.getLocalFileName(), self.surface_tag_map_source.fileformat))            
           return self.__surface_tag_map  
             
                         
72  class RectangularDomain(ParameterSet):  class RectangularDomain(ParameterSet):
73         """         """
74         Generates a mesh over a rectangular domain finley.         Generates a mesh over a rectangular domain finley.
# Line 127  class RectangularDomain(ParameterSet): Line 109  class RectangularDomain(ParameterSet):
109             if self.__domain==None:             if self.__domain==None:
110                if self.dim==2:                if self.dim==2:
111                     self.__domain=finley.Rectangle(n0=self.n[0],\                     self.__domain=finley.Rectangle(n0=self.n[0],\
112                                                  n1=self.n[1],\                                                  n1=self.n[2],\
113                                                  l0=self.l[0],\                                                  l0=self.l[0],\
114                                                  l1=self.l[1],\                                                  l1=self.l[2],\
115                                                  order=self.order, \                                                  order=self.order, \
116                                                  periodic0=self.periodic[0], \                                                  periodic0=self.periodic[0], \
117                                                  periodic1=self.periodic[1], \                                                  periodic1=self.periodic[2], \
118                                                  integrationOrder=self.integrationOrder)                                                  integrationOrder=self.integrationOrder)
119                else:                else:
120                     self.__domain=finley.Brick(n0=self.n[0],\                     self.__domain=finley.Brick(n0=self.n[0],\
# Line 248  class ConstrainerOverBox(Model): Line 230  class ConstrainerOverBox(Model):
230            return self.__value_of_constraint            return self.__value_of_constraint
231                    
232        def __setOutput(self):        def __setOutput(self):
233            x=self.domain.getX()            if self.__location_of_constraint == None:
234            val=self.value               x=self.domain.getX()
235            if isinstance(val, int) or isinstance(val, float):               val=self.value
236               shape=()               if isinstance(val, int) or isinstance(val, float):
237            elif isinstance(val, list) or isinstance(val, tuple) :                  shape=()
238               shape=(len(val),)               elif isinstance(val, list) or isinstance(val, tuple) :
239            elif isinstance(val, numarray.NumArray):                  shape=(len(val),)
240                shape=val.shape               elif isinstance(val, numarray.NumArray):
241            else:                   shape=val.shape
242                shape=val.getShape()               elif val == None:
243            self.__location_of_constraint=Data(0,shape,x.getFunctionSpace())                    shape=()
244            if self.domain.getDim()==3:               else:
245                  x0,x1,x2=x[0],x[1],x[2]                   shape=val.getShape()
246                  if self.left: self.__location_of_constraint+=whereZero(x0-inf(x0),self.tol)               self.__location_of_constraint=Data(0,shape,x.getFunctionSpace())
247                  if self.right: self.__location_of_constraint+=whereZero(x0-sup(x0),self.tol)               if self.domain.getDim()==3:
248                  if self.front: self.__location_of_constraint+=whereZero(x1-inf(x1),self.tol)                     x0,x1,x2=x[0],x[1],x[2]
249                  if self.back: self.__location_of_constraint+=whereZero(x1-sup(x1),self.tol)                     if self.left: self.__location_of_constraint+=whereZero(x0-inf(x0),self.tol)
250                  if self.bottom: self.__location_of_constraint+=whereZero(x2-inf(x2),self.tol)                     if self.right: self.__location_of_constraint+=whereZero(x0-sup(x0),self.tol)
251                  if self.top: self.__location_of_constraint+=whereZero(x2-sup(x2),self.tol)                     if self.front: self.__location_of_constraint+=whereZero(x1-inf(x1),self.tol)
252            else:                     if self.back: self.__location_of_constraint+=whereZero(x1-sup(x1),self.tol)
253                  x0,x1=x[0],x[1]                     if self.bottom: self.__location_of_constraint+=whereZero(x2-inf(x2),self.tol)
254                  if self.left: self.__location_of_constraint+=whereZero(x0-inf(x0),self.tol)                     if self.top: self.__location_of_constraint+=whereZero(x2-sup(x2),self.tol)
255                  if self.right: self.__location_of_constraint+=whereZero(x0-sup(x0),self.tol)               else:
256                  if self.bottom: self.__location_of_constraint+=whereZero(x1-inf(x1),self.tol)                     x0,x1=x[0],x[1]
257                  if self.top: self.__location_of_constraint+=whereZero(x1-sup(x1),self.tol)                     if self.left: self.__location_of_constraint+=whereZero(x0-inf(x0),self.tol)
258            self.__value_of_constraint=self.__location_of_constraint*self.value                     if self.right: self.__location_of_constraint+=whereZero(x0-sup(x0),self.tol)
259                       if self.bottom: self.__location_of_constraint+=whereZero(x1-inf(x1),self.tol)
260                       if self.top: self.__location_of_constraint+=whereZero(x1-sup(x1),self.tol)
261                 if not self.value == None:
262                       self.__value_of_constraint=self.__location_of_constraint*self.value
263  class ScalarConstrainerOverBox(Model):  class ScalarConstrainerOverBox(Model):
264        """        """
265        Creates a characteristic function for the location of constraints        Creates a characteristic function for the location of constraints
# Line 342  class ScalarConstrainerOverBox(Model): Line 328  class ScalarConstrainerOverBox(Model):
328                  if self.right: self.__location_of_constraint+=whereZero(x0-sup(x0),self.tol)                  if self.right: self.__location_of_constraint+=whereZero(x0-sup(x0),self.tol)
329                  if self.bottom: self.__location_of_constraint+=whereZero(x1-inf(x1),self.tol)                  if self.bottom: self.__location_of_constraint+=whereZero(x1-inf(x1),self.tol)
330                  if self.top: self.__location_of_constraint+=whereZero(x1-sup(x1),self.tol)                  if self.top: self.__location_of_constraint+=whereZero(x1-sup(x1),self.tol)
331            if self.value:            if not self.value == None:
332                self.__value_of_constraint=self.__location_of_constraint*self.value                self.__value_of_constraint=self.__location_of_constraint*self.value
333    
334  class VectorConstrainerOverBox(Model):  class VectorConstrainerOverBox(Model):
# Line 370  class VectorConstrainerOverBox(Model): Line 356  class VectorConstrainerOverBox(Model):
356             super(VectorConstrainerOverBox, self).__init__(**kwargs)             super(VectorConstrainerOverBox, self).__init__(**kwargs)
357             self.declareParameter(domain=None, \             self.declareParameter(domain=None, \
358                                   value=None,  \                                   value=None,  \
359                                   left=[0,0,0],  \                                   left=[False ,False ,False ],  \
360                                   right=[0,0,0],  \                                   right=[False ,False ,False ],  \
361                                   top=[0,0,0],  \                                   top=[False ,False ,False ],  \
362                                   bottom=[0,0,0],  \                                   bottom=[False ,False ,False ],  \
363                                   front=[0,0,0], \                                   front=[False ,False ,False ], \
364                                   back=[0,0,0], \                                   back=[False ,False ,False ], \
365                                   tol=1.e-8)                                   tol=1.e-8)
366             self.__value_of_constraint = None             self.__value_of_constraint = None
367             self.__location_of_constraint=None             self.__location_of_constraint=None
# Line 430  class VectorConstrainerOverBox(Model): Line 416  class VectorConstrainerOverBox(Model):
416               if self.top[0]: self.__location_of_constraint+=top_mask*[1.,0.,0.]               if self.top[0]: self.__location_of_constraint+=top_mask*[1.,0.,0.]
417               if self.top[1]: self.__location_of_constraint+=top_mask*[0.,1.,0.]               if self.top[1]: self.__location_of_constraint+=top_mask*[0.,1.,0.]
418               if self.top[2]: self.__location_of_constraint+=top_mask*[0.,0.,1.]               if self.top[2]: self.__location_of_constraint+=top_mask*[0.,0.,1.]
419               if self.value:               if not self.value == None:
420                  self.__value_of_constraint=self.__location_of_constraint*self.value                  self.__value_of_constraint=self.__location_of_constraint*self.value
421            else:            else:
422               x0,x1=x[0],x[1]               x0,x1=x[0],x[1]
# Line 446  class VectorConstrainerOverBox(Model): Line 432  class VectorConstrainerOverBox(Model):
432               top_mask=whereZero(x1-sup(x1),self.tol)               top_mask=whereZero(x1-sup(x1),self.tol)
433               if self.top[0]: self.__location_of_constraint+=top_mask*[1.,0.]               if self.top[0]: self.__location_of_constraint+=top_mask*[1.,0.]
434               if self.top[1]: self.__location_of_constraint+=top_mask*[0.,1.]               if self.top[1]: self.__location_of_constraint+=top_mask*[0.,1.]
435               if self.value:               if not self.value == None:
436                  self.__value_of_constraint=self.__location_of_constraint*self.value[:2]                  self.__value_of_constraint=self.__location_of_constraint*self.value[:2]
437    
438  # vim: expandtab shiftwidth=4:  class ConstrainerAtBoxVertex(Model):
439          """
440          Creates a characteristic function for the location of constraints
441          for all components of a value and selects the value from an initial value
442          ate these locations.
443    
444          In the case that the spatial dimension is two, the arguments front and back are ignored.
445    
446          @ivar domain: domain (in).
447          @ivar tol: absolute tolerance for "x=left, front, bottom vertex" condition, default 1.e-8 (in).
448          """
449          def __init__(self,**kwargs):
450               super(ConstrainerAtBoxVertex, self).__init__(**kwargs)
451               self.declareParameter(domain=None, \
452                                     value=None,  \
453                                     tol=1.e-8)
454               self.__value_of_constraint = None
455               self.__location_of_constraint=None
456          def location_of_constraint(self):
457              """
458              return the values used to constrain a solution
459    
460              @return: the mask marking the locations of the constraints
461              @rtype: L{escript.Scalar}
462              """
463              if self.__location_of_constraint == None: self.__setOutput()
464              return self.__location_of_constraint
465            
466          def value_of_constraint(self):
467              """
468              return the values used to constrain a solution
469    
470              @return: values to be used at the locations of the constraints. If
471                      L{value} is not given C{None} is rerturned.
472              @rtype: L{escript.Scalar}
473              """
474              if self.__location_of_constraint == None: self.__setOutput()
475              return self.__value_of_constraint
476            
477          def __setOutput(self):
478              if self.__location_of_constraint == None:
479                 x=self.domain.getX()
480                 val=self.value
481                 if isinstance(val, int) or isinstance(val, float):
482                    shape=()
483                 elif isinstance(val, list) or isinstance(val, tuple) :
484                    shape=(len(val),)
485                 elif isinstance(val, numarray.NumArray):
486                     shape=val.shape
487                 elif val == None:
488                      shape=()
489                 else:
490                     shape=val.getShape()
491                 if self.domain.getDim()==3:
492                       vertex=[inf(x[0]),inf(x[1]),inf(x[2])]
493                 else:
494                       vertex=[inf(x[0]),inf(x[1])]
495                 self.__location_of_constraint=whereZero(length(x-vertex),self.tol)*numarray.ones(shape)
496                 if not self.value == None:
497                       self.__value_of_constraint=self.__location_of_constraint*self.value
498    class ScalarConstrainerAtBoxVertex(Model):
499          """
500          Creates a characteristic function for the location of constraints
501          for a scalar value and selects the value from an initial value
502          ate these locations.
503    
504          In the case that the spatial dimension is two, the arguments front and back are ignored.
505    
506          @ivar domain: domain (in).
507          @ivar tol: absolute tolerance for "x=left, front, bottom vertex" condition, default 1.e-8 (in).
508          """
509          def __init__(self,**kwargs):
510               super(ScalarConstrainerAtBoxVertex, self).__init__(**kwargs)
511               self.declareParameter(domain=None, \
512                                     value=None,  \
513                                     tol=1.e-8)
514               self.__value_of_constraint = None
515               self.__location_of_constraint=None
516          def location_of_constraint(self):
517              """
518              return the values used to constrain a solution
519    
520              @return: the mask marking the locations of the constraints
521              @rtype: L{escript.Scalar}
522              """
523              if self.__location_of_constraint == None: self.__setOutput()
524              return self.__location_of_constraint
525            
526          def value_of_constraint(self):
527              """
528              return the values used to constrain a solution
529    
530              @return: values to be used at the locations of the constraints. If
531                      L{value} is not given C{None} is rerturned.
532              @rtype: L{escript.Scalar}
533              """
534              if self.__location_of_constraint == None: self.__setOutput()
535              return self.__value_of_constraint
536            
537          def __setOutput(self):
538              x=self.domain.getX()
539              self.__location_of_constraint=Scalar(0,x.getFunctionSpace())
540              if self.domain.getDim()==3:
541                       vertex=[inf(x[0]),inf(x[1]),inf(x[2])]
542              else:
543                     vertex=[inf(x[0]),inf(x[1])]
544              self.__location_of_constraint=whereZero(length(x-vertex),self.tol)
545              if not self.value == None:
546                  self.__value_of_constraint=self.__location_of_constraint*self.value
547    
548    class VectorConstrainerAtBoxVertex(Model):
549          """
550          Creates a characteristic function for the location of constraints vector value.
551          In the case that the spatial dimension is two, the arguments front and
552          back as well as the third component of each argument is ignored.
553    
554          @ivar domain: domain
555          @ivar comp_mask: list of three boolean. comp_mask[i]==True sets a constraint for the i-th component at the left, front, bottom vertex, default [False,False,False] (in).
556          @ivar tol: absolute tolerance for "x=left, front, bottom vertex" condition, default 1.e-8 (in).
557          """
558          def __init__(self, **kwargs):
559               super(VectorConstrainerAtBoxVertex, self).__init__(**kwargs)
560               self.declareParameter(domain=None, \
561                                     value=None,  \
562                                     comp_mask=[False, False, False],
563                                     tol=1.e-8)
564               self.__value_of_constraint = None
565               self.__location_of_constraint=None
566    
567          def location_of_constraint(self):
568              """
569              return the values used to constrain a solution
570    
571              @return: the mask marking the locations of the constraints
572              @rtype: L{escript.Vector}
573              """
574              if self.__location_of_constraint == None: self.__setOutput()
575              return self.__location_of_constraint
576            
577          def value_of_constraint(self):
578              """
579              return the values used to constrain a solution
580    
581              @return: values to be used at the locations of the constraints. If
582                      L{value} is not given C{None} is rerturned.
583              @rtype: L{escript.Vector}
584              """
585              if self.__location_of_constraint == None: self.__setOutput()
586              return self.__value_of_constraint
587            
588          def __setOutput(self):
589              x=self.domain.getX()
590              self.__location_of_constraint=Vector(0,x.getFunctionSpace())
591              if self.domain.getDim()==3:
592                 vertex=[inf(x[0]),inf(x[1]),inf(x[2])]
593                 msk=numarray.zeros((3,))
594                 if self.comp_mask[0]: msk[0]=1
595                 if self.comp_mask[1]: msk[1]=1
596                 if self.comp_mask[2]: msk[2]=1
597              else:
598                 vertex=[inf(x[0]),inf(x[1])]
599                 msk=numarray.zeros((2,))
600                 if self.comp_mask[0]: msk[0]=1
601                 if self.comp_mask[1]: msk[1]=1
602              self.__location_of_constraint=whereZero(length(x-vertex),self.tol)*numarray.ones(shape)
603              if not self.value == None:
604                    self.__value_of_constraint=self.__location_of_constraint*self.value
605    

Legend:
Removed from v.953  
changed lines
  Added in v.1809

  ViewVC Help
Powered by ViewVC 1.1.26