# source:development/pyvolution-1d/domain.py@2716

Last change on this file since 2716 was 2716, checked in by jakeman, 18 years ago

File size: 29.6 KB
Line
1"""Class Domain - 1D domains for finite-volume computations of
2   the shallow water wave equation
3
4
6   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
7   Geoscience Australia
8"""
9from generic_boundary_conditions import *
10from coordinate_transforms.geo_reference import Geo_reference
11
12class Domain:
13
14    def __init__(self, coordinates, boundary = None,
15                 conserved_quantities = None, other_quantities = None,
16                 tagged_elements = None, geo_reference = None):
17        """
18        Build 1D elements from x coordinates
19        """
20
21        from Numeric import array, zeros, Float, Int
22
23        self.beta = 1.0
24        #Store Points
25        self.coordinates = array(coordinates)
26
27        if geo_reference is None:
28            self.geo_reference = Geo_reference() #Use defaults
29        else:
30            self.geo_reference = geo_reference
31
32        #Register number of Elements
33        self.number_of_elements = N = len(self.coordinates)-1
34
35        #Allocate space for neighbour and boundary structures
36        self.neighbours = zeros((N, 2), Int)
37        #self.neighbour_edges = zeros((N, 2), Int)
38        self.neighbour_vertices = zeros((N, 2), Int)
39        self.number_of_boundaries = zeros(N, Int)
40        self.surrogate_neighbours = zeros((N, 2), Int)
41
42        #Allocate space for geometric quantities
43        self.vertices  = zeros((N, 2), Float)
44        self.centroids = zeros(N, Float)
45        self.areas     = zeros(N, Float)
46        for i in range(N):
47            xl = self.coordinates[i]
48            xr = self.coordinates[i+1]
49            self.vertices[i,0] = xl
50            self.vertices[i,1] = xr
51
52            centroid = (xl+xr)/2
53            self.centroids[i] = centroid
54
55            msg = 'Coordinates should be ordered, smallest to largest'
56            assert xr>xl, msg
57
58            self.areas[i] = (xr-xl)
59
60##         print 'N', N
61##         print 'Centroid', self.centroids
62##         print 'Areas', self.areas
63##         print 'Vertex_Coordinates', self.vertices
64
65            #Initialise Neighbours (-1 means that it is a boundary neighbour)
66            self.neighbours[i, :] = [-1, -1]
67            #Initialise edge ids of neighbours
68            #Initialise vertex ids of neighbours
69            #In case of boundaries this slot is not used
70            #self.neighbour_edges[i, :] = [-1, -1]
71            self.neighbour_vertices[i, :] = [-1, -1]
72
73        self.build_vertexlist()
74
75        #Build neighbour structure
76        self.build_neighbour_structure()
77
78        #Build surrogate neighbour structure
79        self.build_surrogate_neighbour_structure()
80
81        #Build boundary dictionary mapping (id, edge) to symbolic tags
82        #Build boundary dictionary mapping (id, vertex) to symbolic tags
83        self.build_boundary_dictionary(boundary)
84
85        #Build tagged element  dictionary mapping (tag) to array of elements
86        self.build_tagged_elements_dictionary(tagged_elements)
87
88        from quantity import Quantity, Conserved_quantity
89
90        #List of quantity names entering
91        #the conservation equations
92        #(Must be a subset of quantities)
93        if conserved_quantities is None:
94            self.conserved_quantities = []
95        else:
96            self.conserved_quantities = conserved_quantities
97
98        if other_quantities is None:
99            self.other_quantities = []
100        else:
101            self.other_quantities = other_quantities
102
103
104        #Build dictionary of Quantity instances keyed by quantity names
105        self.quantities = {}
106
107        #FIXME: remove later - maybe OK, though....
108        for name in self.conserved_quantities:
109            self.quantities[name] = Conserved_quantity(self)
110        for name in self.other_quantities:
111            self.quantities[name] = Quantity(self)
112
113        #Create an empty list for explicit forcing terms
114        self.forcing_terms = []
115
116        #Defaults
117        from config import max_smallsteps, beta_w, beta_h, epsilon, CFL
118        self.beta_w = beta_w
119        self.beta_h = beta_h
120        self.epsilon = epsilon
121
122        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
123        #Or maybe get rid of order altogether and use beta_w and beta_h
124        self.default_order = 1
125        self.order = self.default_order
126        self.smallsteps = 0
127        self.max_smallsteps = max_smallsteps
128        self.number_of_steps = 0
129        self.number_of_first_order_steps = 0
130        self.CFL = CFL
131
132        #Model time
133        self.time = 0.0
134        self.finaltime = None
135        self.min_timestep = self.max_timestep = 0.0
136        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
137        #Checkpointing and storage
140        self.filename = 'domain'
141        self.checkpoint = False
142
143    def __len__(self):
144        return self.number_of_elements
145
146    def build_vertexlist(self):
147        """Build vertexlist index by vertex ids and for each entry (point id)
148        build a list of (triangles, vertex_id) pairs that use the point
149        as vertex.
150
151        Preconditions:
152          self.coordinates and self.triangles are defined
153
154        Postcondition:
155          self.vertexlist is built
156        """
157        from Numeric import array
158
159        vertexlist = [None]*len(self.coordinates)
160        for i in range(self.number_of_elements):
161
162            #a = self.triangles[i, 0]
163            #b = self.triangles[i, 1]
164            #c = self.triangles[i, 2]
165            a = i
166            b = i + 1
167
168            #Register the vertices v as lists of
169            #(triangle_id, vertex_id) tuples associated with them
170            #This is used for smoothing
171            #for vertex_id, v in enumerate([a,b,c]):
172            for vertex_id, v in enumerate([a,b]):
173                if vertexlist[v] is None:
174                    vertexlist[v] = []
175
176                vertexlist[v].append( (i, vertex_id) )
177
178        self.vertexlist = vertexlist
179
180
181    def build_neighbour_structure(self):
182        """Update all registered triangles to point to their neighbours.
183
184        Also, keep a tally of the number of boundaries for each triangle
185
186        Postconditions:
187          neighbours and neighbour_edges is populated
188          neighbours and neighbour_vertices is populated
189          number_of_boundaries integer array is defined.
190        """
191
192        #Step 1:
193        #Build dictionary mapping from segments (2-tuple of points)
194        #to left hand side edge (facing neighbouring triangle)
195
196        N = self.number_of_elements
197        neighbourdict = {}
198        #l_edge = 0
199        #r_edge = 1
200        l_vertex = 0
201        r_vertex = 1
202        for i in range(N):
203
204            #Register all segments as keys mapping to current triangle
205            #and segment id
206            #a = self.triangles[i, 0]
207            #b = self.triangles[i, 1]
208            #c = self.triangles[i, 2]
209            a = self.vertices[i,0]
210            b = self.vertices[i,1]
211
212            """
213            if neighbourdict.has_key((a,b)):
214                    msg = "Edge 2 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[a,b][1],neighbourdict[a,b][0])
215                    raise msg
216            if neighbourdict.has_key((b,c)):
217                    msg = "Edge 0 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[b,c][1],neighbourdict[b,c][0])
218                    raise msg
219            if neighbourdict.has_key((c,a)):
220                    msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0])
221                    raise msg
222                    """
223            #neighbourdict[a,b] = (i, 2) #(id, edge)
224            #neighbourdict[b,c] = (i, 0) #(id, edge)
225            #neighbourdict[c,a] = (i, 1) #(id, edge)
226            #neighbourdict[a,b] = (i, 1) #(id, edge)
227            #neighbourdict[b,a] = (i, 0) #(id, edge)
228            #neighbourdict[a,l_edge] = (i, 0) #(id, edge)
229            #neighbourdict[b,r_edge] = (i, 1) #(id, edge)
230            neighbourdict[a,l_vertex] = (i, 0) #(id, vertex)
231            neighbourdict[b,r_vertex] = (i, 1) #(id, vertex)
232
233
234        #Step 2:
235        #Go through triangles again, but this time
236        #reverse direction of segments and lookup neighbours.
237        for i in range(N):
238            #a = self.triangles[i, 0]
239            #b = self.triangles[i, 1]
240            #c = self.triangles[i, 2]
241
242            a = self.vertices[i,0]
243            b = self.vertices[i,1]
244
245            #self.number_of_boundaries[i] = 3
246            self.number_of_boundaries[i] = 2
247            #if neighbourdict.has_key((b,l_edge)):
248            if neighbourdict.has_key((b,l_vertex)):
249                #self.neighbours[i, 1] = neighbourdict[b,l_edge][0]
250                #self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1]
251                self.neighbours[i, 1] = neighbourdict[b,l_vertex][0]
252                self.neighbour_vertices[i, 1] = neighbourdict[b,l_vertex][1]
253                self.number_of_boundaries[i] -= 1
254
255            #if neighbourdict.has_key((a,r_edge)):
256            if neighbourdict.has_key((a,r_vertex)):
257                #self.neighbours[i, 0] = neighbourdict[a,r_edge][0]
258                #self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1]
259                self.neighbours[i, 0] = neighbourdict[a,r_vertex][0]
260                self.neighbour_vertices[i, 0] = neighbourdict[a,r_vertex][1]
261                self.number_of_boundaries[i] -= 1
262
263            #if neighbourdict.has_key((b,a)):
264            #    self.neighbours[i, 1] = neighbourdict[b,a][0]
265            #    self.neighbour_edges[i, 1] = neighbourdict[b,a][1]
266            #    self.number_of_boundaries[i] -= 1
267
268            #if neighbourdict.has_key((c,b)):
269            #    self.neighbours[i, 0] = neighbourdict[c,b][0]
270            #    self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
271            #    self.number_of_boundaries[i] -= 1
272
273            #if neighbourdict.has_key((a,b)):
274            #    self.neighbours[i, 0] = neighbourdict[a,b][0]
275            #    self.neighbour_edges[i, 0] = neighbourdict[a,b][1]
276            #    self.number_of_boundaries[i] -= 1
277
278    def build_surrogate_neighbour_structure(self):
279        """Build structure where each triangle edge points to its neighbours
280        if they exist. Otherwise point to the triangle itself.
281
282        The surrogate neighbour structure is useful for computing gradients
283        based on centroid values of neighbours.
284
285        Precondition: Neighbour structure is defined
286        Postcondition:
287          Surrogate neighbour structure is defined:
288          surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to
289          triangles.
290
291        """
292
293        N = self.number_of_elements
294        for i in range(N):
295            #Find all neighbouring volumes that are not boundaries
296            #for k in range(3):
297            for k in range(2):
298                if self.neighbours[i, k] < 0:
299                    self.surrogate_neighbours[i, k] = i #Point this triangle
300                else:
301                    self.surrogate_neighbours[i, k] = self.neighbours[i, k]
302
303    def build_boundary_dictionary(self, boundary = None):
304        """Build or check the dictionary of boundary tags.
305         self.boundary is a dictionary of tags,
306         keyed by volume id and edge:
307         { (id, edge): tag, ... }
308
309         Postconditions:
310            self.boundary is defined.
311        """
312
313        from config import default_boundary_tag
314
315        if boundary is None:
316            boundary = {}
317            for vol_id in range(self.number_of_elements):
318                #for edge_id in range(0, 3):
319                #for edge_id in range(0, 2):
320                for vertex_id in range(0, 2):
321                    #if self.neighbours[vol_id, edge_id] < 0:
322                    if self.neighbours[vol_id, vertex_id] < 0:
323                        #boundary[(vol_id, edge_id)] = default_boundary_tag
324                        boundary[(vol_id, vertex_id)] = default_boundary_tag
325        else:
326            #Check that all keys in given boundary exist
327            #for vol_id, edge_id in boundary.keys():
328            for vol_id, vertex_id in boundary.keys():
329                #msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
330                msg = 'Segment (%d, %d) does not exist' %(vol_id, vertex_id)
331                a, b = self.neighbours.shape
332                #assert vol_id < a and edge_id < b, msg
333                assert vol_id < a and vertex_id < b, msg
334
335                #FIXME: This assert violates internal boundaries (delete it)
336                #msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id)
337                #assert self.neighbours[vol_id, edge_id] < 0, msg
338
339            #Check that all boundary segments are assigned a tag
340            for vol_id in range(self.number_of_elements):
341                #for edge_id in range(0, 3):
342                #for edge_id in range(0, 2):
343                for vertex_id in range(0, 2):
344                    #if self.neighbours[vol_id, edge_id] < 0:
345                    if self.neighbours[vol_id, vertex_id] < 0:
346                        #if not boundary.has_key( (vol_id, edge_id) ):
347                        if not boundary.has_key( (vol_id, vertex_id) ):
348                            msg = 'WARNING: Given boundary does not contain '
349                            #msg += 'tags for edge (%d, %d). '\
350                            #       %(vol_id, edge_id)
351                            msg += 'tags for vertex (%d, %d). '\
352                                   %(vol_id, vertex_id)
353                            msg += 'Assigning default tag (%s).'\
354                                   %default_boundary_tag
355
356                            #FIXME: Print only as per verbosity
357                            #print msg
358
359                            #FIXME: Make this situation an error in the future
360                            #and make another function which will
361                            #enable default boundary-tags where
362                            #tags a not specified
363                            #boundary[ (vol_id, edge_id) ] =\
364                            boundary[ (vol_id, vertex_id) ] =\
365                                      default_boundary_tag
366
367
368
369        self.boundary = boundary
370
371    def build_tagged_elements_dictionary(self, tagged_elements = None):
372        """Build the dictionary of element tags.
373         self.tagged_elements is a dictionary of element arrays,
374         keyed by tag:
375         { (tag): [e1, e2, e3..] }
376
377         Postconditions:
378            self.element_tag is defined
379        """
380        from Numeric import array, Int
381
382        if tagged_elements is None:
383            tagged_elements = {}
384        else:
385            #Check that all keys in given boundary exist
386            for tag in tagged_elements.keys():
387                tagged_elements[tag] = array(tagged_elements[tag]).astype(Int)
388
389                msg = 'Not all elements exist. '
390                assert max(tagged_elements[tag]) < self.number_of_elements, msg
391        #print "tagged_elements", tagged_elements
392        self.tagged_elements = tagged_elements
393
394    def get_boundary_tags(self):
395        """Return list of available boundary tags
396        """
397
398        tags = {}
399        for v in self.boundary.values():
400            tags[v] = 1
401
402        return tags.keys()
403
404    def get_vertex_coordinates(self, obj = False):
405        """Return all vertex coordinates.
406        Return all vertex coordinates for all triangles as an Nx6 array
407        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
408
409        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
410        FIXME, we might make that the default.
411        FIXME Maybe use keyword: continuous = False for this condition?
412
413
414        """
415
416        if obj is True:
417            from Numeric import concatenate, reshape
418            #V = self.vertex_coordinates
419            V = self.vertices
420            #return concatenate( (V[:,0:2], V[:,2:4], V[:,4:6]), axis=0)
421
422            N = V.shape[0]
423            #return reshape(V, (3*N, 2))
424            return reshape(V, (N, 2))
425        else:
426            #return self.vertex_coordinates
427            return self.vertices
428
429    def get_conserved_quantities(self, vol_id, vertex=None):#, edge=None):
430        """Get conserved quantities at volume vol_id
431
432        If vertex is specified use it as index for vertex values
433        If edge is specified use it as index for edge values
434        If neither are specified use centroid values
435        If both are specified an exeception is raised
436
437        Return value: Vector of length == number_of_conserved quantities
438
439        """
440
441        from Numeric import zeros, Float
442
443        #if not (vertex is None):# or edge is None):
444        #    msg = 'Values for both vertex and edge was specified.'
445        #    msg += 'Only one (or none) is allowed.'
446       #     raise msg
447
448        q = zeros( len(self.conserved_quantities), Float)
449
450        for i, name in enumerate(self.conserved_quantities):
451            Q = self.quantities[name]
452            if vertex is not None:
453                q[i] = Q.vertex_values[vol_id, vertex]
454            #elif edge is not None:
455            #    q[i] = Q.edge_values[vol_id, edge]
456            else:
457                q[i] = Q.centroid_values[vol_id]
458
459        return q
460
461
462    def get_centroids(self):
463        """Return all coordinates of centroids
464        Return x coordinate of centroid for each element as a N array
465        """
466
467        return self.centroids
468
469    def get_vertices(self):
470        """Return all coordinates of centroids
471        Return x coordinate of centroid for each element as a N array
472        """
473
474        return self.vertices
475
476    def get_coordinate(self, elem_id, vertex=None):
477        """Return coordinate of centroid,
478        or left or right vertex.
479        Left vertex (vertex=0). Right vertex (vertex=1)
480        """
481
482        if vertex is None:
483            return self.centroids[elem_id]
484        else:
485            return self.vertices[elem_id,vertex]
486
487    def get_area(self, elem_id):
488        """Return area of element id
489        """
490
491        return self.areas[elem_id]
492
493    def get_quantity(self, name, location='vertices', indices = None):
494        """Get values for named quantity
495
496        name: Name of quantity
497
498        In case of location == 'centroids' the dimension values must
499        be a list of a Numerical array of length N, N being the number
500        of elements. Otherwise it must be of dimension Nx3.
501
502        Indices is the set of element ids that the operation applies to.
503
504        The values will be stored in elements following their
505        internal ordering.
506        """
507
508        return self.quantities[name].get_values( location, indices = indices)
509
510    def set_quantity(self, name, *args, **kwargs):
511        """Set values for named quantity
512
513
514        One keyword argument is documented here:
515        expression = None, # Arbitrary expression
516
517        expression:
518          Arbitrary expression involving quantity names
519
520        See Quantity.set_values for further documentation.
521        """
522
523        #FIXME (Ole): Allow new quantities here
524        #from quantity import Quantity, Conserved_quantity
525        #Create appropriate quantity object
526        ##if name in self.conserved_quantities:
527        ##    self.quantities[name] = Conserved_quantity(self)
528        ##else:
529        ##    self.quantities[name] = Quantity(self)
530
531
532        #Do the expression stuff
533        if kwargs.has_key('expression'):
534            expression = kwargs['expression']
535            del kwargs['expression']
536
537            Q = self.create_quantity_from_expression(expression)
538            kwargs['quantity'] = Q
539
540
541        #Assign values
542        self.quantities[name].set_values(*args, **kwargs)
543
544    def set_boundary(self, boundary_map):
545        """Associate boundary objects with tagged boundary segments.
546
547        Input boundary_map is a dictionary of boundary objects keyed
548        by symbolic tags to matched against tags in the internal dictionary
549        self.boundary.
550
551        As result one pointer to a boundary object is stored for each vertex
552        in the list self.boundary_objects.
553        More entries may point to the same boundary object
554
555        Schematically the mapping is from two dictionaries to one list
556        where the index is used as pointer to the boundary_values arrays
557        within each quantity.
558
559        self.boundary:          (vol_id, edge_id): tag
560        boundary_map (input):   tag: boundary_object
561        ----------------------------------------------
562        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
563
564
565        Pre-condition:
566          self.boundary has been built.
567
568        Post-condition:
569          self.boundary_objects is built
570
571        If a tag from the domain doesn't appear in the input dictionary an
572        exception is raised.
573        However, if a tag is not used to the domain, no error is thrown.
574        FIXME: This would lead to implementation of a
575        default boundary condition
576
577        Note: If a segment is listed in the boundary dictionary and if it is
578        not None, it *will* become a boundary -
579        even if there is a neighbouring triangle.
580        This would be the case for internal boundaries
581
582        Boundary objects that are None will be skipped.
583
584        FIXME: If set_boundary is called multiple times and if Boundary
585        object is changed into None, the neighbour structure will not be
586        restored!!!
587        """
588
589        self.boundary_objects = []
590
591
592
593
594
595        self.boundary_map = boundary_map  #Store for use with eg. boundary_stats.
596
597        #FIXME: Try to remove the sorting and fix test_mesh.py
598        x = self.boundary.keys()
599        x.sort()
600
601        #Loop through edges that lie on the boundary and associate them with
602        #callable boundary objects depending on their tags
603        #for k, (vol_id, edge_id) in enumerate(x):
604        for k, (vol_id, vertex_id) in enumerate(x):
605            #tag = self.boundary[ (vol_id, edge_id) ]
606            tag = self.boundary[ (vol_id, vertex_id) ]
607
608            if boundary_map.has_key(tag):
609                B = boundary_map[tag]  #Get callable boundary object
610
611                if B is not None:
612                    #self.boundary_objects.append( ((vol_id, edge_id), B) )
613                    #self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
614                    self.boundary_objects.append( ((vol_id, vertex_id), B) )
615                    self.neighbours[vol_id, vertex_id] = -len(self.boundary_objects)
616                else:
617                    pass
618                    #FIXME: Check and perhaps fix neighbour structure
619
620            else:
621                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
622                msg += 'bound to a boundary object.\n'
623                msg += 'All boundary tags defined in domain must appear '
624                msg += 'in the supplied dictionary.\n'
625                msg += 'The tags are: %s' %self.get_boundary_tags()
626                raise msg
627
628
629
630    def check_integrity(self):
631        #Mesh.check_integrity(self)
632
633        for quantity in self.conserved_quantities:
634            msg = 'Conserved quantities must be a subset of all quantities'
635            assert quantity in self.quantities, msg
636
637        ##assert hasattr(self, 'boundary_objects')
638
639    def get_name(self):
640        return self.filename
641
642    def set_name(self, name):
643        self.filename = name
644
647
650
651#Main components of evolve
652
653    def evolve(self, yieldstep = None, finaltime = None,
654               skip_initial_step = False):
655        """Evolve model from time=0.0 to finaltime yielding results
656        every yieldstep.
657
658        Internally, smaller timesteps may be taken.
659
660        Evolve is implemented as a generator and is to be called as such, e.g.
661
662        for t in domain.evolve(timestep, yieldstep, finaltime):
663            <Do something with domain and t>
664
665        """
666
667        from config import min_timestep, max_timestep, epsilon
668
669        #FIXME: Maybe lump into a larger check prior to evolving
670        msg = 'Boundary tags must be bound to boundary objects before evolving system, '
671        msg += 'e.g. using the method set_boundary.\n'
672        msg += 'This system has the boundary tags %s ' %self.get_boundary_tags()
673        assert hasattr(self, 'boundary_objects'), msg
674
675        ##self.set_defaults()
676
677        if yieldstep is None:
678            yieldstep = max_timestep
679        else:
680            yieldstep = float(yieldstep)
681
682        self.order = self.default_order
683
684
685        self.yieldtime = 0.0 #Time between 'yields'
686
687        #Initialise interval of timestep sizes (for reporting only)
688        self.min_timestep = max_timestep
689        self.max_timestep = min_timestep
690        self.finaltime = finaltime
691        self.number_of_steps = 0
692        self.number_of_first_order_steps = 0
693
694        #update ghosts
695        #self.update_ghosts()
696
697        #Initial update of vertex and edge values
698        self.distribute_to_vertices_and_edges()
699
700        #Initial update boundary values
701        self.update_boundary()
702
703        #Or maybe restore from latest checkpoint
704        if self.checkpoint is True:
705            self.goto_latest_checkpoint()
706
707        if skip_initial_step is False:
708            yield(self.time)  #Yield initial values
709
710        while True:
711
712            #Compute fluxes across each element edge
713            self.compute_fluxes()
714
715            #Update timestep to fit yieldstep and finaltime
716            self.update_timestep(yieldstep, finaltime)
717
718            #Update conserved quantities
719            self.update_conserved_quantities()
720
721            #update ghosts
722            #self.update_ghosts()
723
724            #Update vertex and edge values
725            self.distribute_to_vertices_and_edges()
726
727            #Update boundary values
728            self.update_boundary()
729
730            #Update time
731            self.time += self.timestep
732            self.yieldtime += self.timestep
733            self.number_of_steps += 1
734            if self.order == 1:
735                self.number_of_first_order_steps += 1
736
737            #Yield results
738            if finaltime is not None and abs(self.time - finaltime) < epsilon:
739
740                #FIXME: There is a rare situation where the
741                #final time step is stored twice. Can we make a test?
742
743                # Yield final time and stop
744                yield(self.time)
745                break
746
747
748            if abs(self.yieldtime - yieldstep) < epsilon:
749                # Yield (intermediate) time and allow inspection of domain
750
751                if self.checkpoint is True:
752                    self.store_checkpoint()
753                    self.delete_old_checkpoints()
754
755                #Pass control on to outer loop for more specific actions
756                yield(self.time)
757
758                # Reinitialise
759                self.yieldtime = 0.0
760                self.min_timestep = max_timestep
761                self.max_timestep = min_timestep
762                self.number_of_steps = 0
763                self.number_of_first_order_steps = 0
764
765    def distribute_to_vertices_and_edges(self):
766        """Extrapolate conserved quantities from centroid to
767        vertices and edge-midpoints for each volume
768
769        Default implementation is straight first order,
770        i.e. constant values throughout each element and
771        no reference to non-conserved quantities.
772        """
773
774        for name in self.conserved_quantities:
775            Q = self.quantities[name]
776            if self.order == 1:
777                Q.extrapolate_first_order()
778            elif self.order == 2:
779                Q.extrapolate_second_order()
780                Q.limit()
781            else:
782                raise 'Unknown order'
783            Q.interpolate_from_vertices_to_edges()
784
785
786    def update_boundary(self):
787        """Go through list of boundary objects and update boundary values
788        for all conserved quantities on boundary.
789        """
790
791        #FIXME: Update only those that change (if that can be worked out)
792        #FIXME: Boundary objects should not include ghost nodes.
793        #for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
794        #    q = B.evaluate(vol_id, edge_id)
795        for i, ((vol_id, vertex_id), B) in enumerate(self.boundary_objects):
796            q = B.evaluate(vol_id, vertex_id)
797
798            for j, name in enumerate(self.conserved_quantities):
799                Q = self.quantities[name]
800                Q.boundary_values[i] = q[j]
801
802    def compute_forcing_terms(self):
803        """If there are any forcing functions driving the system
804        they should be defined in Domain subclass and appended to
805        the list self.forcing_terms
806        """
807
808        for f in self.forcing_terms:
809            f(self)
810
811
812if __name__ == "__main__":
813
814    points1 = [0.0, 1.0, 2.0, 3.0]
815    D1 = Domain(points1)
816
817    print D1.get_coordinate(0)
818    print D1.get_coordinate(0,1)
819    print 'Number of Elements = ',D1.number_of_elements
820
821    try:
822        print D1.get_coordinate(3)
823    except:
824        pass
825    else:
826        msg =  'Should have raised an out of bounds exception'
827        raise msg
828
829    #points2 = [0.0, 1.0, 2.0, 3.0, 2.5]
830    #D2 = Domain(points2)
831
Note: See TracBrowser for help on using the repository browser.