Changeset 4836


Ignore:
Timestamp:
Nov 20, 2007, 5:38:57 PM (16 years ago)
Author:
ole
Message:

Arranged for timestepping statistic for chosen triangle, e.g. one of the
gauges in the Okushiri example.
The underlying function is currently brute force, but OK for now.

Files:
6 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

    r4835 r4836  
    637637
    638638
    639     def timestepping_statistics(self, track_speeds=False):
     639    def timestepping_statistics(self,
     640                                track_speeds=False,
     641                                triangle_id=None):
    640642        """Return string with time stepping statistics for printing or logging
    641643
    642644        Optional boolean keyword track_speeds decides whether to report
    643645        location of smallest timestep as well as a histogram and percentile
    644         report.
     646        report.
     647
     648        Optional keyword triangle_id can be used to specify a particular
     649        triangle rather than the one with the largest speed.
    645650        """
    646651
     
    716721           
    717722            # Find index of largest computed flux speed
    718             k = self.k = argmax(self.max_speed)
     723            if triangle_id is None:
     724                k = self.k = argmax(self.max_speed)
     725            else:
     726                errmsg = 'Triangle_id %d does not exist in mesh: %s' %(triangle_id,
     727                                                                    str(self))
     728                assert 0 <= triangle_id < len(self), errmsg
     729                k = self.k = triangle_id
     730           
    719731
    720732            x, y = self.get_centroid_coordinates()[k]
     
    725737            msg += '  Triangle #%d with centroid (%.4f, %.4f), ' %(k, x, y)
    726738            msg += 'area = %.4f and radius = %.4f ' %(area, radius)
    727             msg += 'had the largest computed speed: %.6f m/s ' %(max_speed)
    728             msg += 'during last time interval. Quantities below '
    729             msg += 'are reported at their present value, and not what '
    730             msg += 'they were at the time the maximal speed was attained.'
    731             msg += 'To see this, rerun the model with yieldsteps smaller '
    732             msg += 'than the smallest internal timestep reported.'
     739            if triangle_id is None:
     740                msg += 'had the largest computed speed: %.6f m/s ' %(max_speed)
     741            else:
     742                msg += 'had computed speed: %.6f m/s ' %(max_speed)
     743               
    733744            if max_speed > 0.0:
    734745                msg += '(timestep=%.6f)\n' %(radius/max_speed)
  • anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

    r4808 r4836  
    262262    def get_node(self, i,
    263263                 absolute=False):
    264         """Return node coordinates.
    265 
    266         The nodes are ordered in an Nx2 array where N is the number of nodes.
    267         This is the same format they were provided in the constructor
    268         i.e. without any duplication.
     264        """Return node coordinates for triangle i.
    269265
    270266        Boolean keyword argument absolute determines whether coordinates
     
    597593
    598594       
     595       
  • anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

    r4780 r4836  
    860860        return str
    861861
     862    def get_triangle_containing_point(self, point):
     863        """Return triangle id for triangle containing specifiend point (x,y)
     864
     865        If point isn't within mesh, raise exception
     866
     867        """
     868       
     869        # FIXME(Ole): This function is currently brute force
     870        # because I needed it for diagnostics.
     871        # We should make it fast - probably based on the
     872        # quad tree structure.
     873        from anuga.utilities.polygon import is_outside_polygon,\
     874             is_inside_polygon
     875
     876        polygon = self.get_boundary_polygon()
     877        #print polygon, point
     878       
     879        if is_outside_polygon(point, polygon):
     880            msg = 'Point %s is outside mesh' %str(point)
     881            raise Exception, msg
     882
     883
     884        V = self.get_vertex_coordinates(absolute=True)
     885
     886        # FIXME: Horrible brute force
     887        for i, triangle in enumerate(self.triangles):
     888            poly = V[3*i:3*i+3]
     889            #print i, poly
     890
     891            if is_inside_polygon(point, poly, closed=True):
     892                return i
     893               
     894        return
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

    r4515 r4836  
    11951195
    11961196        assert allclose(absolute_points, boundary_polygon)
    1197        
     1197
     1198    def test_get_triangle_containing_point(self):
     1199
     1200        a = [0.0, 0.0]
     1201        b = [0.0, 2.0]
     1202        c = [2.0, 0.0]
     1203        d = [0.0, 4.0]
     1204        e = [2.0, 2.0]
     1205        f = [4.0, 0.0]
     1206
     1207        points = [a, b, c, d, e, f]
     1208        #bac, bce, ecf, dbe
     1209        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     1210        mesh = Mesh(points, vertices)
     1211       
     1212        mesh.check_integrity()
     1213
     1214
     1215        try:
     1216            id = mesh.get_triangle_containing_point([3.0, 5.0])
     1217        except:
     1218            pass
     1219        else:
     1220            msg = 'Should have caught point outside polygon (Non)'           
     1221            raise Exception, msg
     1222           
     1223        id = mesh.get_triangle_containing_point([0.5, 1.0])
     1224        assert id == 0
     1225
     1226        id = mesh.get_triangle_containing_point([1.0, 3.0])
     1227        assert id == 3       
     1228
     1229        for i, point in enumerate(mesh.get_centroid_coordinates()):
     1230            id = mesh.get_triangle_containing_point(point)
     1231            assert id == i       
     1232           
    11981233#-------------------------------------------------------------
    11991234if __name__ == "__main__":
    1200     #suite = unittest.makeSuite(Test_Mesh,'test_boundary_polygon_IIIa')
     1235    #suite = unittest.makeSuite(Test_Mesh,'test_get_triangle_containing_point')
    12011236    suite = unittest.makeSuite(Test_Mesh,'test')
    12021237    runner = unittest.TextTestRunner()
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r4827 r4836  
    462462
    463463       
    464     def timestepping_statistics(self, track_speeds=False):
     464    def timestepping_statistics(self,
     465                                track_speeds=False,
     466                                triangle_id=None):       
    465467        """Return string with time stepping statistics for printing or logging
    466468
     
    475477
    476478        # Call basic machinery from parent class
    477         msg = Generic_Domain.timestepping_statistics(self, track_speeds)
     479        msg = Generic_Domain.timestepping_statistics(self,
     480                                                     track_speeds,
     481                                                     triangle_id)
    478482
    479483        if track_speeds is True:
     
    482486            qwidth = self.qwidth
    483487       
    484             # Triangle with maximum speed
     488            # Selected triangle
    485489            k = self.k
    486490
  • anuga_validation/okushiri_2005/run_okushiri.py

    r4631 r4836  
    8282
    8383
     84# Select triangle containing ch5 for diagnostic output
     85# around known gauge
     86triangle_id = domain.get_triangle_containing_point([4.521, 1.196])
     87# This should get triangle id 32833 with centroid (4.5244, 1.1972)
     88
     89
    8490#-------------------------
    8591# Evolve through time
     
    8995
    9096for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
    91     domain.write_time(track_speeds=False)
     97    print domain.timestepping_statistics(track_speeds=True,
     98                                         triangle_id=triangle_id)
    9299
    93100print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracChangeset for help on using the changeset viewer.