Ignore:
Timestamp:
Apr 24, 2009, 5:22:14 PM (16 years ago)
Author:
rwilson
Message:

Back-merge from Numeric trunk to numpy branch.

Location:
branches/numpy/anuga/abstract_2d_finite_volumes
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/numpy/anuga/abstract_2d_finite_volumes/domain.py

    r6553 r6902  
    9292        """
    9393
     94        number_of_full_nodes=None
     95        number_of_full_triangles=None
     96       
    9497        # Determine whether source is a mesh filename or coordinates
    9598        if type(source) == types.StringType:
     
    117120
    118121        # Expose Mesh attributes (FIXME: Maybe turn into methods)
     122        self.triangles = self.mesh.triangles       
    119123        self.centroid_coordinates = self.mesh.centroid_coordinates
    120124        self.vertex_coordinates = self.mesh.vertex_coordinates
     
    205209        # the ghost triangles.
    206210        if not num.allclose(self.tri_full_flag[:self.number_of_full_nodes], 1):
    207             print ('WARNING:  '
    208                    'Not all full triangles are store before ghost triangles')
     211            if self.numproc>1:
     212                print ('WARNING:  '
     213                       'Not all full triangles are store before ghost triangles')
    209214
    210215        # Defaults
     
    306311    def get_normal(self, *args, **kwargs):
    307312        return self.mesh.get_normal(*args, **kwargs)
     313    def get_triangle_containing_point(self, *args, **kwargs):
     314        return self.mesh.get_triangle_containing_point(*args, **kwargs)
    308315
    309316    def get_intersecting_segments(self, *args, **kwargs):
     
    17451752
    17461753    ##
    1747     # @brief ??
     1754    # @brief Sequential update of ghost cells
    17481755    def update_ghosts(self):
    1749         pass
    1750 
     1756        # We must send the information from the full cells and
     1757        # receive the information for the ghost cells
     1758        # We have a list with ghosts expecting updates
     1759
     1760        from Numeric import take,put
     1761
     1762
     1763        #Update of ghost cells
     1764        iproc = self.processor
     1765        if self.full_send_dict.has_key(iproc):
     1766
     1767            # now store full as local id, global id, value
     1768            Idf  = self.full_send_dict[iproc][0]
     1769
     1770            # now store ghost as local id, global id, value
     1771            Idg = self.ghost_recv_dict[iproc][0]
     1772
     1773            for i, q in enumerate(self.conserved_quantities):
     1774                Q_cv =  self.quantities[q].centroid_values
     1775                put(Q_cv,     Idg, take(Q_cv,     Idf))
     1776
     1777 
    17511778    ##
    17521779    # @brief Extrapolate conserved quantities from centroid to vertices
  • branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r6553 r6902  
    8282    """Time dependent boundary returns values for the
    8383    conserved quantities as a function of time.
    84     Must specify domain to get access to model time and a function f(t)
    85     which must return conserved quantities as a function time
     84    Must specify domain to get access to model time and a function of t
     85    which must return conserved quantities as a function time.
     86   
     87    Example:
     88      B = Time_boundary(domain,
     89                        function=lambda t: [(60<t<3660)*2, 0, 0])
     90     
     91      This will produce a boundary condition with is a 2m high square wave
     92      starting 60 seconds into the simulation and lasting one hour.
     93      Momentum applied will be 0 at all times.
     94                       
    8695    """
    8796
     
    98107        self.verbose = verbose
    99108
     109       
     110        # FIXME: Temporary code to deal with both f and function
     111        if function is not None and f is not None:
     112            raise Exception, 'Specify either function or f to Time_boundary'
     113           
     114        if function is None and f is None:
     115            raise Exception, 'You must specify a function to Time_boundary'
     116           
     117        if f is None:
     118            f = function
     119        #####
     120       
    100121        try:
    101122            q = f(0.0)
  • branches/numpy/anuga/abstract_2d_finite_volumes/mesh_factory.py

    r6304 r6902  
    201201
    202202    return points, elements, boundary
     203
     204
     205def rectangular_periodic(m_g, n_g, len1_g=1.0, len2_g=1.0, origin_g = (0.0, 0.0)):
     206
     207
     208    """Setup a rectangular grid of triangles
     209    with m+1 by n+1 grid points
     210    and side lengths len1, len2. If side lengths are omitted
     211    the mesh defaults to the unit square, divided between all the
     212    processors
     213
     214    len1: x direction (left to right)
     215    len2: y direction (bottom to top)
     216
     217    """
     218
     219    processor = 0
     220    numproc   = 1
     221
     222   
     223    n = n_g
     224    m_low  = -1
     225    m_high = m_g +1
     226
     227    m = m_high - m_low
     228
     229    delta1 = float(len1_g)/m_g
     230    delta2 = float(len2_g)/n_g
     231
     232    len1 = len1_g*float(m)/float(m_g)
     233    len2 = len2_g
     234    origin = ( origin_g[0]+float(m_low)/float(m_g)*len1_g, origin_g[1] )
     235
     236    #Calculate number of points
     237    Np = (m+1)*(n+1)
     238
     239    class VIndex:
     240
     241        def __init__(self, n,m):
     242            self.n = n
     243            self.m = m
     244
     245        def __call__(self, i,j):
     246            return j+i*(self.n+1)
     247
     248    class EIndex:
     249
     250        def __init__(self, n,m):
     251            self.n = n
     252            self.m = m
     253
     254        def __call__(self, i,j):
     255            return 2*(j+i*self.n)
     256
     257
     258    I = VIndex(n,m)
     259    E = EIndex(n,m)
     260
     261    points = num.zeros( (Np,2), num.Float)
     262
     263    for i in range(m+1):
     264        for j in range(n+1):
     265
     266            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
     267
     268    #Construct 2 triangles per rectangular element and assign tags to boundary
     269    #Calculate number of triangles
     270    Nt = 2*m*n
     271
     272
     273    elements = num.zeros( (Nt,3), num.Int)
     274    boundary = {}
     275    Idgl = []
     276    Idfl = []
     277    Idgr = []
     278    Idfr = []
     279
     280    full_send_dict = {}
     281    ghost_recv_dict = {}
     282    nt = -1
     283    for i in range(m):
     284        for j in range(n):
     285
     286            i1 = I(i,j+1)
     287            i2 = I(i,j)
     288            i3 = I(i+1,j+1)
     289            i4 = I(i+1,j)
     290
     291            #Lower Element
     292            nt = E(i,j)
     293            if i == 0:
     294                Idgl.append(nt)
     295
     296            if i == 1:
     297                Idfl.append(nt)
     298
     299            if i == m-2:
     300                Idfr.append(nt)
     301
     302            if i == m-1:
     303                Idgr.append(nt)
     304
     305            if i == m-1:
     306                if processor == numproc-1:
     307                    boundary[nt, 2] = 'right'
     308                else:
     309                    boundary[nt, 2] = 'ghost'
     310       
     311            if j == 0:
     312                boundary[nt, 1] = 'bottom'
     313            elements[nt,:] = [i4,i3,i2]
     314
     315            #Upper Element
     316            nt = E(i,j)+1
     317            if i == 0:
     318                Idgl.append(nt)
     319
     320            if i == 1:
     321                Idfl.append(nt)
     322
     323            if i == m-2:
     324                Idfr.append(nt)
     325
     326            if i == m-1:
     327                Idgr.append(nt)
     328
     329            if i == 0:
     330                if processor == 0:
     331                    boundary[nt, 2] = 'left'
     332                else:
     333                    boundary[nt, 2] = 'ghost'
     334            if j == n-1:
     335                boundary[nt, 1] = 'top'
     336            elements[nt,:] = [i1,i2,i3]
     337
     338    Idfl.extend(Idfr)
     339    Idgr.extend(Idgl)
     340
     341    Idfl = num.array(Idfl, num.Int)
     342    Idgr = num.array(Idgr, num.Int)
     343   
     344    full_send_dict[processor]  = [Idfl, Idfl]
     345    ghost_recv_dict[processor] = [Idgr, Idgr]
     346
     347
     348    return  points, elements, boundary, full_send_dict, ghost_recv_dict
     349
    203350
    204351def oblique(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):
  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py

    r6689 r6902  
    11781178        # Edges have already been deprecated in set_values, see changeset:5521,
    11791179        # but *might* be useful in get_values. Any thoughts anyone?
     1180        # YES (Ole): Edge values are necessary for volumetric balance check
    11801181
    11811182        if location not in ['vertices', 'centroids',
  • branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py

    r6553 r6902  
    397397                                     
    398398    def test_setting_timestepping_method(self):
    399         """test_set_quanitities_to_be_monitored
     399        """test_setting_timestepping_method
    400400        """
    401401
     
    803803                             [ 11.0,  11.0,  11.0]])
    804804                             
     805    def test_rectangular_periodic_and_ghosts(self):
     806
     807        from mesh_factory import rectangular_periodic
     808       
     809
     810        M=5
     811        N=2
     812        points, vertices, boundary, full_send_dict, ghost_recv_dict = rectangular_periodic(M, N)
     813
     814        assert num.allclose(ghost_recv_dict[0][0], [24, 25, 26, 27,  0,  1,  2,  3])
     815        assert num.allclose(full_send_dict[0][0] , [ 4,  5,  6,  7, 20, 21, 22, 23])
     816
     817        #print 'HERE'
     818       
     819        conserved_quantities = ['quant1', 'quant2']
     820        domain = Domain(points, vertices, boundary, conserved_quantities,
     821                        full_send_dict=full_send_dict,
     822                        ghost_recv_dict=ghost_recv_dict)
     823
     824
     825
     826
     827        assert num.allclose(ghost_recv_dict[0][0], [24, 25, 26, 27,  0,  1,  2,  3])
     828        assert num.allclose(full_send_dict[0][0] , [ 4,  5,  6,  7, 20, 21, 22, 23])
     829
     830        def xylocation(x,y):
     831            return 15*x + 9*y
     832
     833       
     834        domain.set_quantity('quant1',xylocation,location='centroids')
     835        domain.set_quantity('quant2',xylocation,location='centroids')
     836
     837
     838        assert num.allclose(domain.quantities['quant1'].centroid_values,
     839                            [  0.5,   1.,   5.,    5.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
     840                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
     841                               18.5,  19.,   23.,   23.5])
     842
     843
     844
     845        assert num.allclose(domain.quantities['quant2'].centroid_values,
     846                            [  0.5,   1.,   5.,    5.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
     847                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
     848                               18.5,  19.,   23.,   23.5])
     849
     850        domain.update_ghosts()
     851
     852
     853        assert num.allclose(domain.quantities['quant1'].centroid_values,
     854                            [  15.5,  16.,   20.,   20.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
     855                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
     856                                3.5,   4.,    8.,    8.5])
     857
     858
     859
     860        assert num.allclose(domain.quantities['quant2'].centroid_values,
     861                            [  15.5,  16.,   20.,   20.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
     862                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
     863                                3.5,   4.,    8.,    8.5])
     864
     865       
     866        assert num.allclose(domain.tri_full_flag, [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     867                                                   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0])
     868
     869        #Test that points are arranged in a counter clock wise order
     870        domain.check_integrity()
     871
     872
    805873    def test_that_mesh_methods_exist(self):
    806874        """test_that_mesh_methods_exist
     
    828896        domain.get_number_of_nodes()
    829897        domain.get_normal(0,0)
     898        domain.get_triangle_containing_point([0.4,0.5])
    830899        domain.get_intersecting_segments([[0.0, 0.0], [0.0, 1.0]])
    831900        domain.get_disconnected_triangles()
Note: See TracChangeset for help on using the changeset viewer.