Changeset 659


Ignore:
Timestamp:
Dec 2, 2004, 2:11:36 PM (20 years ago)
Author:
ole
Message:

Added functionality for settimng quantities according to polygons

Location:
inundation/ga/storm_surge
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/examples/beach.py

    r653 r659  
    2121from pmesh2domain import pmesh_to_domain_instance
    2222from config import data_dir
     23from util import read_polygon, Polygon_function
    2324
    2425######################
    … …  
    4950
    5051
    51 def inside(p,polygon):
    52     return True
    5352
    54 def friction(x, y):
    55     """Read friction from polygon file and set accordingly
    56     """
    5753
    58     from Numeric import array
    59     x = array(x).astype('f')
    60     y = array(y).astype('f')   
     54p0 = [[20,20], [40,20], [40,40], [20,40]]         
    6155   
    62     #Get polygon
    63     fid = open('beach.poly')
    64     lines = fid.readlines()
    65     fid.close()
    66     polygon = []
    67     for line in lines:
    68         fields = line.split(',')
    69         polygon.append( [float(fields[0]), float(fields[1])] )
    70 
    71     print polygon
    72 
    73     z = 0*x
    74     for i in range(len(x)):
    75         if inside( (x[i], y[i]), polygon):
    76             z[i] = 1
    77         else:
    78             z[i] = 0
    79 
    80 friction([0],[1])   
    81 
    82 domain.set_quantity('elevation', x_slope)
    83 domain.set_quantity('level', x_slope)
     56domain.set_quantity('elevation', Polygon_function( [(p0, 50)]), 'centroids' )
     57domain.set_quantity('level', Polygon_function( [(p0, 50)] ), 'centroids' )
    8458domain.set_quantity('friction', 0.07)
    8559
    8660print domain.get_extent()
    87 import sys; sys.exit()
    8861
    8962######################

  • inundation/ga/storm_surge/pyvolution/domain.py

    r648 r659  
    121121            self.set_quantity(key, quantity_dict[key], location='vertices')
    122122       
     123
    123124    def set_quantity(self, name, X, location='vertices', indexes = None):
    124125        """Set values for named quantity

  • inundation/ga/storm_surge/pyvolution/flatbed.py

    r505 r659  
    1010from mesh_factory import rectangular
    1111from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
    12      Constant_height, wind_Stress
     12     Constant_height
    1313from Numeric import array
     14from util import Polygon_function, read_polygon
    1415
    1516#Create basic mesh
    16 N = 150
    17 points, vertices, boundary = rectangular(N, N, 1000, 1000)
     17N = 20
     18points, vertices, boundary = rectangular(N, N, 100, 100)
    1819
    1920#Create shallow water domain
    2021domain = Domain(points, vertices, boundary)
    21 domain.smooth = True
    22 domain.visualise = False
    2322domain.store = True
    24 domain.default_order=1
     23domain.set_name('polygons')
     24domain.default_order=2
    2525
    2626#Set driving forces
    … …  
    2929inflow_level = 10.0
    3030domain.set_quantity('friction', manning)
     31
     32
     33
     34
     35#Define polynomials
     36p0 = [[20,27], [30,25], [40,40], [20,40]]         
     37p1 = [[80,19], [90,20], [85,50], [80,55], [75,58], [70,60], [60,24]]         
     38p2 = read_polygon('testpoly.txt')
     39
     40
     41#Set elevation 
     42domain.set_quantity('elevation',
     43        Polygon_function([(p0, 23), (p1,10), (p2,15)]))
     44
     45
     46         
    3147
    3248
    … …  
    4359######################
    4460#Evolution
    45 for t in domain.evolve(yieldstep = 10, finaltime = 500):
     61for t in domain.evolve(yieldstep = 1, finaltime = 1):
    4662    domain.write_time()
    4763

  • inundation/ga/storm_surge/pyvolution/quantity.py

    r593 r659  
    147147            #Intialise centroid and edge_values
    148148            self.interpolate()
    149            
     149           
     150        if location == 'centroids':
     151            #Extrapolate 1st order - to capture notion of area being specified
     152            self.extrapolate_first_order()           
    150153         
    151154    def get_values(self, location='vertices', indexes = None):
    … …  
    479482
    480483           
     484    def extrapolate_first_order(self):
     485        """Extrapolate conserved quantities from centroid to
     486        vertices for each volume using
     487        first order scheme.
     488        """
     489       
     490        qc = self.centroid_values
     491        qv = self.vertex_values
     492
     493        for i in range(3):
     494            qv[:,i] = qc
    481495           
    482496
    … …  
    486500
    487501    boundary values, storage and method for updating, and
    488     methods for extrapolation from centropid to vertices inluding
     502    methods for (second order) extrapolation from centroid to vertices inluding
    489503    gradients and limiters
    490504    """
    … …  
    525539       
    526540       
    527     def extrapolate_first_order(self):
    528         """Extrapolate conserved quantities from centroid to
    529         vertices for each volume using
    530         first order scheme.
    531         """
    532        
    533         qc = self.centroid_values
    534         qv = self.vertex_values
    535 
    536         for i in range(3):
    537             qv[:,i] = qc
    538            
    539            
    540541    def extrapolate_second_order(self):
    541542        #Call correct module function

  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r623 r659  
    8181        bed = self.quantities['elevation']       
    8282
    83         msg = 'All water levels must be greater than the bed elevation'
    84         assert alltrue( greater_equal(
    85             level.vertex_values, bed.vertex_values )), msg
    86 
    87         assert alltrue( greater_equal(
    88             level.edge_values, bed.edge_values )), msg
    89 
    90         assert alltrue( greater_equal(
    91             level.centroid_values, bed.centroid_values )), msg       
     83        #msg = 'All water levels must be greater than the bed elevation'
     84        #assert alltrue( greater_equal(
     85        #    level.vertex_values, bed.vertex_values )), msg
     86        #
     87        #assert alltrue( greater_equal(
     88        #    level.edge_values, bed.edge_values )), msg
     89        #
     90        #assert alltrue( greater_equal(
     91        #    level.centroid_values, bed.centroid_values )), msg       
    9292
    9393

  • inundation/ga/storm_surge/pyvolution/test_util.py

    r656 r659  
    356356       
    357357    def test_point_on_line(self):
    358                
     358               
    359359        #Endpoints first       
    360360        assert point_on_line( (0, 0), [(0,0), (1,0)] ) 
    … …  
    371371        assert not point_on_line( (0, 0.5), [(0,0), (1,1)] )
    372372       
     373        #From real example that failed
     374        assert not point_on_line( (40,50), [(40,20), (40,40)] )
     375       
     376       
     377        #From real example that failed
     378        assert not point_on_line( (40,19), [(40,20), (40,40)] )
     379                       
     380       
    373381       
    374382       
    375383    def test_inside_polygon(self):
    376    
     384
     385       
    377386        #Simplest case: Polygon is the unit square
    378387        polygon = [[0,0], [1,0], [1,1], [0,1]]
    … …  
    394403        assert not inside_polygon( (0.5, 0.), polygon, closed=False)   
    395404        assert not inside_polygon( (1., 0.5), polygon, closed=False)   
    396        
     405
     406
     407   
     408        #From real example (that failed)
     409        polygon = [[20,20], [40,20], [40,40], [20,40]]         
     410        points = [ [40, 50] ]
     411        res = inside_polygon(points, polygon)
     412        assert len(res) == 0
     413       
     414        polygon = [[20,20], [40,20], [40,40], [20,40]]         
     415        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
     416        res = inside_polygon(points, polygon)
     417        assert len(res) == 2
     418        assert allclose(res, [0,1])
     419       
     420               
    397421             
    398422        #More convoluted and non convex polygon
    … …  
    455479        #assert not inside_polygon( (0.5, 0.5), polygon )       
    456480       
    457         points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]                   
    458        
    459481
    460482                     

  • inundation/ga/storm_surge/pyvolution/util.py

    r656 r659  
    7070    y1 = line[1][1]
    7171
    72    
    73     a = array([x - x0, y - y0])                 
     72    a = array([x - x0, y - y0])
     73    a_normal = array([a[1], -a[0]])
     74                       
    7475    b = array([x1 - x0, y1 - y0])
    7576   
    76     len_a = sqrt(sum(a**2))           
    77     len_b = sqrt(sum(b**2))                 
    78                
    79     if len_a == 0 or len_a == len_b:
    80        return True
     77    if dot(a_normal, b) == 0:
     78        #Point is somewhere on the infinite extension of the line
     79
     80        len_a = sqrt(sum(a**2))               
     81        len_b = sqrt(sum(b**2))                                 
     82        if dot(a, b) >= 0 and len_a <= len_b:
     83           return True
     84        else:   
     85           return False
    8186    else:
    82        x = dot(a, b)/len_b       
    83        return allclose(x, len_a)
    84 
    85        
    86    
     87      return False 
     88   
     89
    8790def inside_polygon(point, polygon, closed = True):
    8891    """Determine whether points are inside or outside a polygon
    … …  
    155158    py = polygon[:,1]   
    156159
    157                        
     160   
    158161    #Begin algorithm
    159162    indices = []
    … …  
    162165        x = points[k, 0]
    163166        y = points[k, 1]       
     167
    164168       
    165169        inside = False
    … …  
    257261        for polygon, value in self.regions:
    258262            indices = inside_polygon(points, polygon)
    259            
    260263            for i in indices:
    261264                z[i] = value
    262265
     266        #from sys import exit; exit()   
    263267        return z                 
     268
     269def read_polygon(filename):
     270    """Read points assumed to form a polygon
     271       There must be exactly two numbers in each line
     272    """             
     273   
     274    #Get polygon
     275    fid = open(filename)
     276    lines = fid.readlines()
     277    fid.close()
     278    polygon = []
     279    for line in lines:
     280        fields = line.split(',')
     281        polygon.append( [float(fields[0]), float(fields[1])] )
     282
     283    return polygon     
     284             
    264285             
    265286class File_function:
Note: See TracChangeset for help on using the changeset viewer.