source: inundation/pyvolution/domain.py @ 2569

Last change on this file since 2569 was 2569, checked in by ole, 18 years ago

Allowed None boundary objects (to accommodate ghost elements in the parallel implementation).
For sequential implementations, None would pave the way for a default boundary.

File size: 26.5 KB
Line 
1"""Class Domain - 2D triangular domains for finite-volume computations of
2   the shallow water wave equation
3
4
5   Copyright 2004
6   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
7   Geoscience Australia
8"""
9
10from mesh import Mesh
11from generic_boundary_conditions import Boundary
12from generic_boundary_conditions import File_boundary
13from generic_boundary_conditions import Dirichlet_boundary
14from generic_boundary_conditions import Time_boundary
15from generic_boundary_conditions import Transmissive_boundary
16
17import types
18
19class Domain(Mesh):
20
21    def __init__(self, coordinates, vertices, boundary = None,
22                 conserved_quantities = None, other_quantities = None,
23                 tagged_elements = None, geo_reference = None,
24                 use_inscribed_circle=False):
25
26        Mesh.__init__(self, coordinates, vertices, boundary,
27                      tagged_elements, geo_reference, use_inscribed_circle)
28
29        from Numeric import zeros, Float, Int
30        from quantity import Quantity, Conserved_quantity
31
32        #List of quantity names entering
33        #the conservation equations
34        #(Must be a subset of quantities)
35        if conserved_quantities is None:
36            self.conserved_quantities = []
37        else:
38            self.conserved_quantities = conserved_quantities
39
40        if other_quantities is None:
41            self.other_quantities = []
42        else:
43            self.other_quantities = other_quantities
44
45
46        #Build dictionary of Quantity instances keyed by quantity names
47        self.quantities = {}
48
49        #FIXME: remove later - maybe OK, though....
50        for name in self.conserved_quantities:
51            self.quantities[name] = Conserved_quantity(self)
52        for name in self.other_quantities:
53            self.quantities[name] = Quantity(self)
54
55        #Create an empty list for explicit forcing terms
56        self.forcing_terms = []
57
58
59        #Defaults
60        from config import max_smallsteps, beta_w, beta_h, epsilon, CFL
61        self.beta_w = beta_w
62        self.beta_h = beta_h
63        self.epsilon = epsilon
64
65        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
66        #Or maybe get rid of order altogether and use beta_w and beta_h
67        self.set_default_order(1)       
68        #self.default_order = 1
69        #self.order = self.default_order
70
71        self.smallsteps = 0
72        self.max_smallsteps = max_smallsteps
73        self.number_of_steps = 0
74        self.number_of_first_order_steps = 0
75        self.CFL = CFL
76
77        #Model time
78        self.time = 0.0
79        self.finaltime = None
80        self.min_timestep = self.max_timestep = 0.0
81        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
82
83        ######OBSOLETE
84        #Origin in UTM coordinates
85        #FIXME: This should be set if read by a msh file
86        #self.zone = zone
87        #self.xllcorner = xllcorner
88        #self.yllcorner = yllcorner
89
90
91        #Checkpointing and storage
92        from config import default_datadir
93        self.datadir = default_datadir
94        self.filename = 'domain'
95        self.checkpoint = False
96
97        #MH310505 To avoid calculating the flux across each edge twice, keep an integer (boolean) array,
98        #to be used during the flux calculation
99        N=self.number_of_elements
100        self.already_computed_flux = zeros((N, 3), Int)
101
102
103    def set_default_order(self, n):
104        """Set default (spatial) order to either 1 or 2
105        """
106
107        msg = 'Default order must be either 1 or 2. I got %s' %n
108        assert n in [1,2], msg
109       
110        self.default_order = n
111        self.order = self.default_order
112
113
114    #Public interface to Domain
115    def get_conserved_quantities(self, vol_id, vertex=None, edge=None):
116        """Get conserved quantities at volume vol_id
117
118        If vertex is specified use it as index for vertex values
119        If edge is specified use it as index for edge values
120        If neither are specified use centroid values
121        If both are specified an exeception is raised
122
123        Return value: Vector of length == number_of_conserved quantities
124
125        """
126
127        from Numeric import zeros, Float
128
129        if not (vertex is None or edge is None):
130            msg = 'Values for both vertex and edge was specified.'
131            msg += 'Only one (or none) is allowed.'
132            raise msg
133
134        q = zeros( len(self.conserved_quantities), Float)
135
136        for i, name in enumerate(self.conserved_quantities):
137            Q = self.quantities[name]
138            if vertex is not None:
139                q[i] = Q.vertex_values[vol_id, vertex]
140            elif edge is not None:
141                q[i] = Q.edge_values[vol_id, edge]
142            else:
143                q[i] = Q.centroid_values[vol_id]
144
145        return q
146
147    def set_time(self, time=0.0):
148        """Set the time"""
149
150        self.time = time
151
152    def set_quantity_vertices_dict(self, quantity_dict):
153        """Set values for named quantities.
154        The index is the quantity
155
156        name: Name of quantity
157        X: Compatible list, Numeric array, const or function (see below)
158
159        The values will be stored in elements following their
160        internal ordering.
161
162        """
163        for key in quantity_dict.keys():
164            self.set_quantity(key, quantity_dict[key], location='vertices')
165
166
167    def set_quantity(self, name, *args, **kwargs):
168        """Set values for named quantity
169
170
171        One keyword argument is documented here:
172        expression = None, # Arbitrary expression
173
174        expression:
175          Arbitrary expression involving quantity names
176
177        See Quantity.set_values for further documentation.
178        """
179
180        #FIXME (Ole): Allow new quantities here
181        #from quantity import Quantity, Conserved_quantity
182        #Create appropriate quantity object
183        ##if name in self.conserved_quantities:
184        ##    self.quantities[name] = Conserved_quantity(self)
185        ##else:
186        ##    self.quantities[name] = Quantity(self)
187
188
189        #Do the expression stuff
190        if kwargs.has_key('expression'):
191            expression = kwargs['expression']
192            del kwargs['expression']
193
194            Q = self.create_quantity_from_expression(expression)
195            kwargs['quantity'] = Q
196
197
198        #Assign values
199        self.quantities[name].set_values(*args, **kwargs)
200
201
202    def get_quantity(self, name, location='vertices', indices = None):
203        """Get quantity object.
204
205        name: Name of quantity
206
207        See methods inside the quantity object for more options
208        """
209
210        return self.quantities[name] #.get_values( location, indices = indices)
211
212
213    def get_quantity_object(self, name):
214        """Get object for named quantity
215
216        name: Name of quantity
217
218        FIXME: Obsolete
219        """
220
221        print 'get_quantity_object has been deprecated. Please use get_quantity'
222        return self.quantities[name]
223
224
225    def create_quantity_from_expression(self, expression):
226        """Create new quantity from other quantities using arbitrary expression
227
228        Combine existing quantities in domain using expression and return
229        result as a new quantity.
230
231        Note, the new quantity could e.g. be used in set_quantity
232
233        Valid expressions are limited to operators defined in class Quantity
234
235        Example:
236
237
238        """
239
240        from util import apply_expression_to_dictionary
241        return apply_expression_to_dictionary(expression, self.quantities)
242
243
244
245
246    def set_boundary(self, boundary_map):
247        """Associate boundary objects with tagged boundary segments.
248
249        Input boundary_map is a dictionary of boundary objects keyed
250        by symbolic tags to matched against tags in the internal dictionary
251        self.boundary.
252
253        As result one pointer to a boundary object is stored for each vertex
254        in the list self.boundary_objects.
255        More entries may point to the same boundary object
256
257        Schematically the mapping is from two dictionaries to one list
258        where the index is used as pointer to the boundary_values arrays
259        within each quantity.
260
261        self.boundary:          (vol_id, edge_id): tag
262        boundary_map (input):   tag: boundary_object
263        ----------------------------------------------
264        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
265
266
267        Pre-condition:
268          self.boundary has been built.
269
270        Post-condition:
271          self.boundary_objects is built
272
273        If a tag from the domain doesn't appear in the input dictionary an
274        exception is raised.
275        However, if a tag is not used to the domain, no error is thrown.
276        FIXME: This would lead to implementation of a
277        default boundary condition
278
279        Note: If a segment is listed in the boundary dictionary and if it is
280        not None, it *will* become a boundary -
281        even if there is a neighbouring triangle.
282        This would be the case for internal boundaries
283
284        Boundary objects that are None will be skipped.
285
286        FIXME: If set_boundary is called multiple times and if Boundary
287        object is changed into None, the neighbour structure will not be
288        restored!!!
289
290
291        """
292
293        self.boundary_objects = []
294        self.boundary_map = boundary_map  #Store for use with eg. boundary_stats.
295
296        #FIXME: Try to remove the sorting and fix test_mesh.py
297        x = self.boundary.keys()
298        x.sort()
299
300        #Loop through edges that lie on the boundary and associate them with
301        #callable boundary objects depending on their tags
302        for k, (vol_id, edge_id) in enumerate(x):
303            tag = self.boundary[ (vol_id, edge_id) ]
304
305            if boundary_map.has_key(tag):
306                B = boundary_map[tag]  #Get callable boundary object
307
308                if B is not None:
309                    self.boundary_objects.append( ((vol_id, edge_id), B) )
310                    self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
311                else:
312                    pass
313                    #FIXME: Check and perhaps fix neighbour structure
314
315
316            else:
317                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
318                msg += 'bound to a boundary object.\n'
319                msg += 'All boundary tags defined in domain must appear '
320                msg += 'in the supplied dictionary.\n'
321                msg += 'The tags are: %s' %self.get_boundary_tags()
322                raise msg
323
324
325    def set_region(self, functions):
326        # The order of functions in the list is used.
327        if type(functions) not in [types.ListType,types.TupleType]:
328            functions = [functions]
329        for function in functions:
330            for tag in self.tagged_elements.keys():
331                function(tag, self.tagged_elements[tag], self)
332
333                #Do we need to do this sort of thing?
334                #self = function(tag, self.tagged_elements[tag], self)
335
336    #MISC
337    def check_integrity(self):
338        Mesh.check_integrity(self)
339
340        for quantity in self.conserved_quantities:
341            msg = 'Conserved quantities must be a subset of all quantities'
342            assert quantity in self.quantities, msg
343
344        ##assert hasattr(self, 'boundary_objects')
345
346    def write_time(self):
347        print self.timestepping_statistics()
348
349        #Old version
350        #if self.min_timestep == self.max_timestep:
351        #    print 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
352        #          %(self.time, self.min_timestep, self.number_of_steps,
353        #            self.number_of_first_order_steps)
354        #elif self.min_timestep > self.max_timestep:
355        #    print 'Time = %.4f, steps=%d (%d)'\
356        #          %(self.time, self.number_of_steps,
357        #            self.number_of_first_order_steps)
358        #else:
359        #    print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
360        #          %(self.time, self.min_timestep,
361        #            self.max_timestep, self.number_of_steps,
362        #            self.number_of_first_order_steps)
363
364    def timestepping_statistics(self):
365        """Return string with time stepping statistics for printing or logging
366        """
367
368        msg = ''
369        if self.min_timestep == self.max_timestep:
370            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
371                   %(self.time, self.min_timestep, self.number_of_steps,
372                     self.number_of_first_order_steps)
373        elif self.min_timestep > self.max_timestep:
374            msg += 'Time = %.4f, steps=%d (%d)'\
375                   %(self.time, self.number_of_steps,
376                     self.number_of_first_order_steps)
377        else:
378            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
379                   %(self.time, self.min_timestep,
380                     self.max_timestep, self.number_of_steps,
381                     self.number_of_first_order_steps)
382
383        return msg
384
385
386    def write_boundary_statistics(self, quantities = None, tags = None):
387        print self.boundary_statistics(quantities, tags)
388
389    def boundary_statistics(self, quantities = None, tags = None):
390        """Output statistics about boundary forcing at each timestep
391
392
393        Input:
394          quantities: either None, a string or a list of strings naming the quantities to be reported
395          tags:       either None, a string or a list of strings naming the tags to be reported         
396
397
398        Example output:
399        Tag 'wall':
400            stage in [2, 5.5]
401            xmomentum in []
402            ymomentum in []
403        Tag 'ocean'
404
405
406        If quantities are specified only report on those. Otherwise take all conserved quantities.
407        If tags are specified only report on those, otherwise take all tags.
408
409        """
410
411        #Input checks
412        import types, string
413
414        if quantities is None:
415            quantities = self.conserved_quantities
416        elif type(quantities) == types.StringType:
417            quantities = [quantities] #Turn it into a list
418
419        msg = 'Keyword argument quantities must be either None, '
420        msg += 'string or list. I got %s' %str(quantities)
421        assert type(quantities) == types.ListType, msg
422
423
424        if tags is None:
425            tags = self.get_boundary_tags()
426        elif type(tags) == types.StringType:
427            tags = [tags] #Turn it into a list
428
429        msg = 'Keyword argument tags must be either None, '
430        msg += 'string or list. I got %s' %str(tags)
431        assert type(tags) == types.ListType, msg
432
433        #Determine width of longest quantity name (for cosmetic purposes)
434        maxwidth = 0
435        for name in quantities:
436            w = len(name)
437            if w > maxwidth:
438                maxwidth = w
439
440        #Output stats
441        msg = 'Boundary values at time %.4f:\n' %self.time
442        for tag in tags:
443            msg += '    %s:\n' %tag
444
445            for name in quantities:
446                q = self.quantities[name]
447
448                #Find range of boundary values for tag and q
449                maxval = minval = None
450                for i, ((vol_id, edge_id), B) in\
451                        enumerate(self.boundary_objects):
452                    if self.boundary[(vol_id, edge_id)] == tag:
453                        v = q.boundary_values[i]
454                        if minval is None or v < minval: minval = v
455                        if maxval is None or v > maxval: maxval = v
456
457                if minval is None or maxval is None:
458                    msg += '        Sorry no information available about' +\
459                           ' tag %s and quantity %s\n' %(tag, name)
460                else:
461                    msg += '        %s in [%12.8f, %12.8f]\n'\
462                           %(string.ljust(name, maxwidth), minval, maxval)
463
464
465        return msg
466
467
468    def get_name(self):
469        return self.filename
470
471    def set_name(self, name):
472        self.filename = name
473
474    def get_datadir(self):
475        return self.datadir
476
477    def set_datadir(self, name):
478        self.datadir = name
479
480
481
482    #def set_defaults(self):
483    #    """Set default values for uninitialised quantities.
484    #    Should be overridden or specialised by specific modules
485    #    """#
486    #
487    #    for name in self.conserved_quantities + self.other_quantities:
488    #        self.set_quantity(name, 0.0)
489
490
491    ###########################
492    #Main components of evolve
493
494    def evolve(self,
495               yieldstep = None,
496               finaltime = None,
497               duration = None,             
498               skip_initial_step = False):
499        """Evolve model through time starting from self.starttime.
500       
501       
502        yieldstep: Interval between yields where results are stored,
503                   statistics written and domain inspected or
504                   possibly modified. If omitted the internal predefined
505                   max timestep is used.
506                   Internally, smaller timesteps may be taken.               
507
508        duration: Duration of simulation
509
510        finaltime: Time where simulation should end
511
512        If both duration and finaltime are given an exception is thrown.
513
514
515        skip_initial_step: Boolean flag that decides whether the first
516        yield step is skipped or not. This is useful for example to avoid
517        duplicate steps when multiple evolve processes are dove tailed.
518
519
520        Evolve is implemented as a generator and is to be called as such, e.g.
521
522        for t in domain.evolve(yieldstep, finaltime):
523            <Do something with domain and t>
524
525
526        All times are given in seconds       
527
528        """
529
530        from config import min_timestep, max_timestep, epsilon
531
532        #FIXME: Maybe lump into a larger check prior to evolving
533        msg = 'Boundary tags must be bound to boundary objects before evolving system, '
534        msg += 'e.g. using the method set_boundary.\n'
535        msg += 'This system has the boundary tags %s '\
536               %self.get_boundary_tags()
537        assert hasattr(self, 'boundary_objects'), msg
538
539        ##self.set_defaults()
540
541        if yieldstep is None:
542            yieldstep = max_timestep
543        else:
544            yieldstep = float(yieldstep)
545
546        self.order = self.default_order
547
548
549        if finaltime is not None and duration is not None:
550            print 'F', finaltime, duration
551            msg = 'Only one of finaltime and duration may be specified'
552            raise msg
553        else:
554            if finaltime is not None:
555                self.finaltime = float(finaltime)
556            if duration is not None:
557                self.finaltime = self.starttime + float(duration)
558               
559
560       
561
562        self.yieldtime = 0.0 #Time between 'yields'
563
564        #Initialise interval of timestep sizes (for reporting only)
565        self.min_timestep = max_timestep
566        self.max_timestep = min_timestep
567        self.number_of_steps = 0
568        self.number_of_first_order_steps = 0
569
570        #update ghosts
571        self.update_ghosts()
572
573        #Initial update of vertex and edge values
574        self.distribute_to_vertices_and_edges()
575
576        #Initial update boundary values
577        self.update_boundary()
578
579        #Or maybe restore from latest checkpoint
580        if self.checkpoint is True:
581            self.goto_latest_checkpoint()
582
583        if skip_initial_step is False:
584            yield(self.time)  #Yield initial values
585
586        while True:
587
588            #Compute fluxes across each element edge
589            self.compute_fluxes()
590
591            #Update timestep to fit yieldstep and finaltime
592            self.update_timestep(yieldstep, finaltime)
593
594            #Update conserved quantities
595            self.update_conserved_quantities()
596
597            #update ghosts
598            self.update_ghosts()
599
600            #Update vertex and edge values
601            self.distribute_to_vertices_and_edges()
602
603            #Update boundary values
604            self.update_boundary()
605
606            #Update time
607            self.time += self.timestep
608            self.yieldtime += self.timestep
609            self.number_of_steps += 1
610            if self.order == 1:
611                self.number_of_first_order_steps += 1
612
613            #Yield results
614            if finaltime is not None and abs(self.time - finaltime) < epsilon:
615
616                #FIXME: There is a rare situation where the
617                #final time step is stored twice. Can we make a test?
618
619                # Yield final time and stop
620                yield(self.time)
621                break
622
623
624            if abs(self.yieldtime - yieldstep) < epsilon:
625                # Yield (intermediate) time and allow inspection of domain
626
627                if self.checkpoint is True:
628                    self.store_checkpoint()
629                    self.delete_old_checkpoints()
630
631                #Pass control on to outer loop for more specific actions
632                yield(self.time)
633
634                # Reinitialise
635                self.yieldtime = 0.0
636                self.min_timestep = max_timestep
637                self.max_timestep = min_timestep
638                self.number_of_steps = 0
639                self.number_of_first_order_steps = 0
640
641
642    def evolve_to_end(self, finaltime = 1.0):
643        """Iterate evolve all the way to the end
644        """
645
646        for _ in self.evolve(yieldstep=None, finaltime=finaltime):
647            pass
648
649
650
651    def update_boundary(self):
652        """Go through list of boundary objects and update boundary values
653        for all conserved quantities on boundary.
654        """
655
656        #FIXME: Update only those that change (if that can be worked out)
657        #FIXME: Boundary objects should not include ghost nodes.
658        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
659            if B is None:
660                print 'WARNING: Ignored boundary segment %d (None)'
661            else:               
662                q = B.evaluate(vol_id, edge_id)
663
664                for j, name in enumerate(self.conserved_quantities):
665                    Q = self.quantities[name]
666                    Q.boundary_values[i] = q[j]
667
668
669    def compute_fluxes(self):
670        msg = 'Method compute_fluxes must be overridden by Domain subclass'
671        raise msg
672
673
674    def update_timestep(self, yieldstep, finaltime):
675
676        from config import min_timestep
677
678        # self.timestep is calculated from speed of characteristics
679        # Apply CFL condition here
680        timestep = self.CFL*self.timestep
681
682        #Record maximal and minimal values of timestep for reporting
683        self.max_timestep = max(timestep, self.max_timestep)
684        self.min_timestep = min(timestep, self.min_timestep)
685
686        #Protect against degenerate time steps
687        if timestep < min_timestep:
688
689            #Number of consecutive small steps taken b4 taking action
690            self.smallsteps += 1
691
692            if self.smallsteps > self.max_smallsteps:
693                self.smallsteps = 0 #Reset
694
695                if self.order == 1:
696                    msg = 'WARNING: Too small timestep %.16f reached '\
697                          %timestep
698                    msg += 'even after %d steps of 1 order scheme'\
699                           %self.max_smallsteps
700                    print msg
701                    timestep = min_timestep  #Try enforcing min_step
702
703                    #raise msg
704                else:
705                    #Try to overcome situation by switching to 1 order
706                    self.order = 1
707
708        else:
709            self.smallsteps = 0
710            if self.order == 1 and self.default_order == 2:
711                self.order = 2
712
713
714        #Ensure that final time is not exceeded
715        if finaltime is not None and self.time + timestep > finaltime:
716            timestep = finaltime-self.time
717
718        #Ensure that model time is aligned with yieldsteps
719        if self.yieldtime + timestep > yieldstep:
720            timestep = yieldstep-self.yieldtime
721
722        self.timestep = timestep
723
724
725
726    def compute_forcing_terms(self):
727        """If there are any forcing functions driving the system
728        they should be defined in Domain subclass and appended to
729        the list self.forcing_terms
730        """
731
732        for f in self.forcing_terms:
733            f(self)
734
735
736
737    def update_conserved_quantities(self):
738        """Update vectors of conserved quantities using previously
739        computed fluxes specified forcing functions.
740        """
741
742        from Numeric import ones, sum, equal, Float
743
744        N = self.number_of_elements
745        d = len(self.conserved_quantities)
746
747        timestep = self.timestep
748
749        #Compute forcing terms
750        self.compute_forcing_terms()
751
752        #Update conserved_quantities
753        for name in self.conserved_quantities:
754            Q = self.quantities[name]
755            Q.update(timestep)
756
757            #Clean up
758            #Note that Q.explicit_update is reset by compute_fluxes
759
760            #MH090605 commented out the following since semi_implicit_update is now re-initialized
761            #at the end of the _update function in quantity_ext.c (This is called by the
762            #preceeding Q.update(timestep) statement above).
763            #For run_profile.py with N=128, the time of update_conserved_quantities is cut from 14.00 secs
764            #to 8.35 secs
765
766            #Q.semi_implicit_update[:] = 0.0
767
768    def update_ghosts(self):
769        pass
770
771    def distribute_to_vertices_and_edges(self):
772        """Extrapolate conserved quantities from centroid to
773        vertices and edge-midpoints for each volume
774
775        Default implementation is straight first order,
776        i.e. constant values throughout each element and
777        no reference to non-conserved quantities.
778        """
779
780        for name in self.conserved_quantities:
781            Q = self.quantities[name]
782            if self.order == 1:
783                Q.extrapolate_first_order()
784            elif self.order == 2:
785                Q.extrapolate_second_order()
786                Q.limit()
787            else:
788                raise 'Unknown order'
789            Q.interpolate_from_vertices_to_edges()
790
791
792    def centroid_norm(self, quantity, normfunc):
793        """Calculate the norm of the centroid values
794        of a specific quantity, using normfunc.
795
796        normfunc should take a list to a float.
797
798        common normfuncs are provided in the module utilities.norms
799        """
800        return normfunc(self.quantities[quantity].centroid_values)
801
802
803
804##############################################
805#Initialise module
806
807#Optimisation with psyco
808from config import use_psyco
809if use_psyco:
810    try:
811        import psyco
812    except:
813        import os
814        if os.name == 'posix' and os.uname()[4] == 'x86_64':
815            pass
816            #Psyco isn't supported on 64 bit systems, but it doesn't matter
817        else:
818            msg = 'WARNING: psyco (speedup) could not import'+\
819                  ', you may want to consider installing it'
820            print msg
821    else:
822        psyco.bind(Domain.update_boundary)
823        #psyco.bind(Domain.update_timestep)     #Not worth it
824        psyco.bind(Domain.update_conserved_quantities)
825        psyco.bind(Domain.distribute_to_vertices_and_edges)
826
827
828if __name__ == "__main__":
829    pass
Note: See TracBrowser for help on using the repository browser.