Changeset 12888 for issm/trunk-jpl/src/m/utils/Geometry/FlagElements.py

Timestamp:

08/03/12 10:21:22 (13 years ago)

Author:

jschierm

Message:

Running python versions of setflowequation (and more).

File:

• issm/trunk-jpl/src/m/utils/Geometry/FlagElements.py (modified) (1 diff)

issm/trunk-jpl/src/m/utils/Geometry/FlagElements.py

@@ -1 @@
-from numpy import * 
+import numpy
+import os
+#from basinzoom import *
+#from ContourToMesh import *
+from MatlabFuncs import *
+
 def FlagElements(md,region):
-#FLAGELEMENTS - flag the elements in an region
-#
-#   The region can be given with an exp file, a list of elements.
-#
-#   Usage: 
-#      flag=FlagElements(md,region);
-#
-#   Example:
-#      flag=FlagElements(md,'all');
-#      flag=FlagElements(md,'');
-#      flag=FlagElements(md,'Domain.exp');
-#      flag=FlagElements(md,'~Domain.exp');
-#      flag=FlagElements(md,md.mask.elementongroundedice);
+        """
+        FLAGELEMENTS - flag the elements in an region
 
-        if isinstance(region,basestring):
-                if not(region):
-                        flag=zeros(md.mesh.numberofelements,'bool')
-                        invert=0;
-                elif region=='all':
-                        flag=ones(md.mesh.numberofelements,'bool')
-                        invert=0;
+            The region can be given with an exp file, a list of elements.
+
+            Usage: 
+               flag=FlagElements(md,region);
+
+            Example:
+               flag=FlagElements(md,'all');
+               flag=FlagElements(md,'');
+               flag=FlagElements(md,'Domain.exp');
+               flag=FlagElements(md,'~Domain.exp');
+               flag=FlagElements(md,md.mask.elementongroundedice);
+        """
+
+        if   isinstance(region,str):
+                if   not region:
+                        flag=numpy.zeros(md.mesh.numberofelements,'bool')
+                        invert=0
+                elif strcmpi(region,'all'):
+                        flag=numpy.ones(md.mesh.numberofelements,'bool')
+                        invert=0
                 else:
                         #make sure that we actually don't want the elements outside the domain outline!
-                        if region[0]=='~':
+                        if strcmpi(region[0],'~'):
                                 region=region[1:]
-                                invert=1;
+                                invert=1
                         else:
-                                invert=0;
-                        
+                                invert=0
+
                         #does the region domain outline exist or do we have to look for xlim,ylim in basinzoom?
-                        if not os.path.isfile(region):
-                                [xlim,ylim]=basinzoom('basin',region);
-                                flag_nodes=double(md.mesh.x<xlim(2) & md.mesh.x>xlim(1) &  md.mesh.y<ylim(2) & md.mesh.y>ylim(1));
-                                flag=prod(flag_nodes(md.mesh.elements),2);
+                        if not os.path.exists(region):
+                                if len(region)>3 and not strcmp(region[-4:],'.exp'):
+                                        raise IOError("Error: File 'region' not found!" % region)
+                                xlim,ylim=basinzoom('basin',region)
+                                flag_nodes=numpy.logical_and(numpy.logical_and(md.mesh.x<xlim[1],md.mesh.x>xlim[0]),numpy.logical_and(md.mesh.y<ylim[1],md.mesh.y>ylim[0])).astype(float)
+                                flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1)
                         else:
                                 #ok, flag elements
-                                flag=ContourToMesh(md.mesh.elements[:,0:3],md.mesh.x,md.mesh.y,region,'element',1);
-                
+                                flag=ContourToMesh(md.mesh.elements[:,0:3],md.mesh.x,md.mesh.y,region,'element',1)
+
                 if invert:
-                        flag=~flag;
-        
-        elif isinstance(region,nparray):
-                if len(region)!=md.mesh.numberofelements:
-                        print FlagElements.__doc__
-                        print 'Flaglist for region must be of same size as number of elements in model'
-                        return []
-                flag=region;
+                        flag=numpy.logical_not(flag)
+
+        elif isinstance(region,numpy.nparray) or isinstance(region,bool):
+                if not numpy.size(region,0)==md.mesh.numberofelements:
+                        raise TypeError("Flaglist for region must be of same size as number of elements in model.")
+                flag=region
+
         else:
-                print 'Invalid region option'
-                return []
+                raise TypeError("Invalid region option")
 
-        return flag;
+        return flag
+