source: trunk/anuga_core/source/anuga/structures/structure.py @ 7984

Last change on this file since 7984 was 7983, checked in by habili, 14 years ago

structure classes to deal with irregular structures

File size: 3.9 KB
Line 
1from anuga.geometry.polygon import inside_polygon, polygon_area
2from anuga.config import g
3import anuga.utilities.log as log
4import inlet
5import numpy as num
6import math
7
8class Structure:
9    """Culvert flow - transfer water from one rectangular box to another.
10    Sets up the geometry of problem
11   
12    This is the base class for culverts. Inherit from this class (and overwrite
13    compute_discharge method for specific subclasses)
14   
15    Input: Two points, pipe_size (either diameter or width, height),
16    mannings_rougness,
17    """ 
18
19    def __init__(self,
20                 domain,
21                 end_points, 
22                 width=None,
23                 height=None,
24                 verbose=False):
25       
26        # Input check
27       
28        self.domain = domain
29
30        self.end_points = end_points
31 
32        self.width = width
33        self.height = height
34       
35        self.verbose=verbose
36       
37        # Create the fundamental culvert polygons and create inlet objects
38        self.__create_culvert_polygons()
39
40        #FIXME (SR) Put this into a foe loop to deal with more inlets
41
42        self.inlets = []
43        polygon0 = self.inlet_polygons[0]
44        outward_vector0 = self.culvert_vector
45        self.inlets.append(inlet.Inlet(self.domain, polygon0, outward_vector0))
46
47        polygon1 = self.inlet_polygons[1]
48        outward_vector1  = - self.culvert_vector
49        self.inlets.append(inlet.Inlet(self.domain, polygon1, outward_vector1))
50
51
52    def __call__(self):
53        msg = 'Method __call__ must be overridden by Culvert subclass'
54        raise Exception, msg
55
56    def __create_structure_polygons(self):
57
58        """Create polygons at the end of a culvert inlet and outlet.
59        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
60        for enquiring the total energy at both ends of the culvert and transferring flow.
61        """
62
63        # Calculate geometry
64        x0, y0 = self.end_points[0]
65        x1, y1 = self.end_points[1]
66
67        dx = x1 - x0
68        dy = y1 - y0
69
70        self.culvert_vector = num.array([dx, dy])
71        self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))
72        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
73
74        # Unit direction vector and normal
75        self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
76        self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
77
78        # Short hands
79        w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
80        h = 0.5*self.height*self.culvert_vector    # Vector of length=height in the
81                             # direction of the culvert
82
83        self.inlet_polygons = []
84
85        # Build exchange polygon and enquiry points 0 and 1
86        for i in [0, 1]:
87            i0 = (2*i-1)
88            p0 = self.end_points[i] + w
89            p1 = self.end_points[i] - w
90            p2 = p1 + i0*h
91            p3 = p0 + i0*h
92            self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
93
94        # Check that enquiry points are outside inlet polygons
95        for i in [0,1]:
96            polygon = self.inlet_polygons[i]
97            # FIXME (SR) Probably should calculate the area of all the triangles
98            # associated with this polygon, as there is likely to be some
99            # inconsistency between triangles and ploygon
100            area = polygon_area(polygon)
101           
102            msg = 'Polygon %s ' %(polygon)
103            msg += ' has area = %f' % area
104            assert area > 0.0, msg
105           
106   
107    def get_inlets(self):
108       
109        return self.inlets
110       
111       
112    def get_structure_length(self):
113       
114        return self.culvert_length
115       
116       
117    def get_structure_width(self):
118       
119        return self.width
120       
121       
122    def get_structure_height(self):
123   
124        return self.height
125       
Note: See TracBrowser for help on using the repository browser.