Changeset 2705


Ignore:
Timestamp:
Apr 12, 2006, 10:58:17 AM (19 years ago)
Author:
jakeman
Message:

domain with working neighbour structure

Location:
development/pyvolution-1d
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed
  • development/pyvolution-1d/domain.py

    r2229 r2705  
    77   Geoscience Australia
    88"""
     9from generic_boundary_conditions import *
    910
    1011class Domain:
    1112
    12     def __init__(self, coordinates, conserved_quantities = None,
    13                  other_quantities = None):
     13    def __init__(self, coordinates, boundary = None,
     14                 conserved_quantities = None,
     15                 other_quantities = None,
     16                 tagged_elements = None):
    1417        """
    1518        Build 1D elements from x coordinates
     
    2528        self.number_of_elements = N = len(self.coordinates)-1
    2629
    27 
     30        #Allocate space for neighbour and boundary structures
     31        self.neighbours = zeros((N, 2), Int)
     32        self.neighbour_edges = zeros((N, 2), Int)
     33        self.number_of_boundaries = zeros(N, Int)
     34        self.surrogate_neighbours = zeros((N, 2), Int)
     35       
    2836        #Allocate space for geometric quantities
    2937        self.vertices  = zeros((N, 2), Float)
     
    4351
    4452            self.areas[i] = (xr-xl)
    45 
     53           
    4654##         print 'N', N
    4755##         print 'Centroid', self.centroids
    4856##         print 'Areas', self.areas
    4957##         print 'Vertex_Coordinates', self.vertices
     58           
     59            #Initialise Neighbours (-1 means that it is a boundary neighbour)
     60            self.neighbours[i, :] = [-1, -1]
     61            #Initialise edge ids of neighbours
     62            #In case of boundaries this slot is not used
     63            self.neighbour_edges[i, :] = [-1, -1]
     64
     65        #Build neighbour structure
     66        self.build_neighbour_structure()
     67
     68        #Build surrogate neighbour structure
     69        self.build_surrogate_neighbour_structure()         
     70
     71        #Build boundary dictionary mapping (id, edge) to symbolic tags
     72        self.build_boundary_dictionary(boundary)
     73
     74        #Build tagged element  dictionary mapping (tag) to array of elements
     75        self.build_tagged_elements_dictionary(tagged_elements)
     76
     77
     78       
     79        from quantity import Quantity, Conserved_quantity
     80       
     81        #List of quantity names entering
     82        #the conservation equations
     83        #(Must be a subset of quantities)
     84        if conserved_quantities is None:
     85            self.conserved_quantities = []
     86        else:
     87            self.conserved_quantities = conserved_quantities
     88
     89        if other_quantities is None:
     90            self.other_quantities = []
     91        else:
     92            self.other_quantities = other_quantities
     93
     94
     95        #Build dictionary of Quantity instances keyed by quantity names
     96        self.quantities = {}
     97
     98        #FIXME: remove later - maybe OK, though....
     99        for name in self.conserved_quantities:
     100            self.quantities[name] = Conserved_quantity(self)
     101        for name in self.other_quantities:
     102            self.quantities[name] = Quantity(self)
     103
     104        #Create an empty list for explicit forcing terms
     105        self.forcing_terms = []
     106
     107    def build_neighbour_structure(self):
     108        """Update all registered triangles to point to their neighbours.
     109
     110        Also, keep a tally of the number of boundaries for each triangle
     111
     112        Postconditions:
     113          neighbours and neighbour_edges is populated
     114          number_of_boundaries integer array is defined.
     115        """
     116
     117        #Step 1:
     118        #Build dictionary mapping from segments (2-tuple of points)
     119        #to left hand side edge (facing neighbouring triangle)
     120
     121        N = self.number_of_elements
     122        neighbourdict = {}
     123        l_edge = 0
     124        r_edge = 1
     125        for i in range(N):
     126
     127            #Register all segments as keys mapping to current triangle
     128            #and segment id
     129            #a = self.triangles[i, 0]
     130            #b = self.triangles[i, 1]
     131            #c = self.triangles[i, 2]
     132            a = self.vertices[i,0]
     133            b = self.vertices[i,1]
     134           
     135            """
     136            if neighbourdict.has_key((a,b)):
     137                    msg = "Edge 2 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[a,b][1],neighbourdict[a,b][0])
     138                    raise msg
     139            if neighbourdict.has_key((b,c)):
     140                    msg = "Edge 0 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[b,c][1],neighbourdict[b,c][0])
     141                    raise msg
     142            if neighbourdict.has_key((c,a)):
     143                    msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0])
     144                    raise msg
     145                    """
     146            #neighbourdict[a,b] = (i, 2) #(id, edge)
     147            #neighbourdict[b,c] = (i, 0) #(id, edge)
     148            #neighbourdict[c,a] = (i, 1) #(id, edge)
     149            #neighbourdict[a,b] = (i, 1) #(id, edge)
     150            #neighbourdict[b,a] = (i, 0) #(id, edge)
     151            neighbourdict[a,l_edge] = (i, 0) #(id, edge)
     152            neighbourdict[b,r_edge] = (i, 1) #(id, edge)
     153
     154
     155        #Step 2:
     156        #Go through triangles again, but this time
     157        #reverse direction of segments and lookup neighbours.
     158        for i in range(N):
     159            #a = self.triangles[i, 0]
     160            #b = self.triangles[i, 1]
     161            #c = self.triangles[i, 2]
     162           
     163            a = self.vertices[i,0]
     164            b = self.vertices[i,1]
     165           
     166            #self.number_of_boundaries[i] = 3
     167            self.number_of_boundaries[i] = 2
     168            if neighbourdict.has_key((b,l_edge)):
     169                self.neighbours[i, 1] = neighbourdict[b,l_edge][0]
     170                self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1]
     171                self.number_of_boundaries[i] -= 1
     172
     173            if neighbourdict.has_key((a,r_edge)):
     174                self.neighbours[i, 0] = neighbourdict[a,r_edge][0]
     175                self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1]
     176                self.number_of_boundaries[i] -= 1
     177            #if neighbourdict.has_key((b,a)):
     178            #    self.neighbours[i, 1] = neighbourdict[b,a][0]
     179            #    self.neighbour_edges[i, 1] = neighbourdict[b,a][1]
     180            #    self.number_of_boundaries[i] -= 1
     181
     182            #if neighbourdict.has_key((c,b)):
     183            #    self.neighbours[i, 0] = neighbourdict[c,b][0]
     184            #    self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
     185            #    self.number_of_boundaries[i] -= 1
     186
     187            #if neighbourdict.has_key((a,b)):
     188            #    self.neighbours[i, 0] = neighbourdict[a,b][0]
     189            #    self.neighbour_edges[i, 0] = neighbourdict[a,b][1]
     190            #    self.number_of_boundaries[i] -= 1
     191
     192    def build_surrogate_neighbour_structure(self):
     193        """Build structure where each triangle edge points to its neighbours
     194        if they exist. Otherwise point to the triangle itself.
     195
     196        The surrogate neighbour structure is useful for computing gradients
     197        based on centroid values of neighbours.
     198
     199        Precondition: Neighbour structure is defined
     200        Postcondition:
     201          Surrogate neighbour structure is defined:
     202          surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to
     203          triangles.
     204
     205        """
     206
     207        N = self.number_of_elements
     208        for i in range(N):
     209            #Find all neighbouring volumes that are not boundaries
     210            #for k in range(3):
     211            for k in range(2):   
     212                if self.neighbours[i, k] < 0:
     213                    self.surrogate_neighbours[i, k] = i #Point this triangle
     214                else:
     215                    self.surrogate_neighbours[i, k] = self.neighbours[i, k]
     216
     217    def build_boundary_dictionary(self, boundary = None):
     218        """Build or check the dictionary of boundary tags.
     219         self.boundary is a dictionary of tags,
     220         keyed by volume id and edge:
     221         { (id, edge): tag, ... }
     222
     223         Postconditions:
     224            self.boundary is defined.
     225        """
     226
     227        from config import default_boundary_tag
     228
     229        if boundary is None:
     230            boundary = {}
     231            for vol_id in range(self.number_of_elements):
     232                #for edge_id in range(0, 3):
     233                for edge_id in range(0, 2):
     234                    if self.neighbours[vol_id, edge_id] < 0:
     235                        boundary[(vol_id, edge_id)] = default_boundary_tag
     236        else:
     237            #Check that all keys in given boundary exist
     238            for vol_id, edge_id in boundary.keys():
     239                msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
     240                a, b = self.neighbours.shape
     241                assert vol_id < a and edge_id < b, msg
     242
     243                #FIXME: This assert violates internal boundaries (delete it)
     244                #msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id)
     245                #assert self.neighbours[vol_id, edge_id] < 0, msg
     246
     247            #Check that all boundary segments are assigned a tag
     248            for vol_id in range(self.number_of_elements):
     249                #for edge_id in range(0, 3):
     250                for edge_id in range(0, 2):   
     251                    if self.neighbours[vol_id, edge_id] < 0:
     252                        if not boundary.has_key( (vol_id, edge_id) ):
     253                            msg = 'WARNING: Given boundary does not contain '
     254                            msg += 'tags for edge (%d, %d). '\
     255                                   %(vol_id, edge_id)
     256                            msg += 'Assigning default tag (%s).'\
     257                                   %default_boundary_tag
     258
     259                            #FIXME: Print only as per verbosity
     260                            #print msg
     261
     262                            #FIXME: Make this situation an error in the future
     263                            #and make another function which will
     264                            #enable default boundary-tags where
     265                            #tags a not specified
     266                            boundary[ (vol_id, edge_id) ] =\
     267                                      default_boundary_tag
     268
     269
     270
     271        self.boundary = boundary
     272
     273    def build_tagged_elements_dictionary(self, tagged_elements = None):
     274        """Build the dictionary of element tags.
     275         self.tagged_elements is a dictionary of element arrays,
     276         keyed by tag:
     277         { (tag): [e1, e2, e3..] }
     278
     279         Postconditions:
     280            self.element_tag is defined
     281        """
     282        from Numeric import array, Int
     283
     284        if tagged_elements is None:
     285            tagged_elements = {}
     286        else:
     287            #Check that all keys in given boundary exist
     288            for tag in tagged_elements.keys():
     289                tagged_elements[tag] = array(tagged_elements[tag]).astype(Int)
     290
     291                msg = 'Not all elements exist. '
     292                assert max(tagged_elements[tag]) < self.number_of_elements, msg
     293        #print "tagged_elements", tagged_elements
     294        self.tagged_elements = tagged_elements
     295
     296    def get_boundary_tags(self):
     297        """Return list of available boundary tags
     298        """
     299
     300        tags = {}
     301        for v in self.boundary.values():
     302            tags[v] = 1
     303
     304        return tags.keys()
     305       
     306
     307    def get_conserved_quantities(self, vol_id, vertex=None, edge=None):
     308        """Get conserved quantities at volume vol_id
     309
     310        If vertex is specified use it as index for vertex values
     311        If edge is specified use it as index for edge values
     312        If neither are specified use centroid values
     313        If both are specified an exeception is raised
     314
     315        Return value: Vector of length == number_of_conserved quantities
     316
     317        """
     318
     319        from Numeric import zeros, Float
     320
     321        if not (vertex is None or edge is None):
     322            msg = 'Values for both vertex and edge was specified.'
     323            msg += 'Only one (or none) is allowed.'
     324            raise msg
     325
     326        q = zeros( len(self.conserved_quantities), Float)
     327
     328        for i, name in enumerate(self.conserved_quantities):
     329            Q = self.quantities[name]
     330            if vertex is not None:
     331                q[i] = Q.vertex_values[vol_id, vertex]
     332            elif edge is not None:
     333                q[i] = Q.edge_values[vol_id, edge]
     334            else:
     335                q[i] = Q.centroid_values[vol_id]
     336
     337        return q
     338       
    50339
    51340    def get_centroids(self):
     
    79368
    80369        return self.areas[elem_id]
     370
     371    def set_quantity(self, name, *args, **kwargs):
     372        """Set values for named quantity
     373
     374
     375        One keyword argument is documented here:
     376        expression = None, # Arbitrary expression
     377
     378        expression:
     379          Arbitrary expression involving quantity names
     380
     381        See Quantity.set_values for further documentation.
     382        """
     383
     384        #FIXME (Ole): Allow new quantities here
     385        #from quantity import Quantity, Conserved_quantity
     386        #Create appropriate quantity object
     387        ##if name in self.conserved_quantities:
     388        ##    self.quantities[name] = Conserved_quantity(self)
     389        ##else:
     390        ##    self.quantities[name] = Quantity(self)
     391
     392
     393        #Do the expression stuff
     394        if kwargs.has_key('expression'):
     395            expression = kwargs['expression']
     396            del kwargs['expression']
     397
     398            Q = self.create_quantity_from_expression(expression)
     399            kwargs['quantity'] = Q
     400
     401
     402        #Assign values
     403        self.quantities[name].set_values(*args, **kwargs)
     404
     405    def set_boundary(self, boundary_map):
     406        """Associate boundary objects with tagged boundary segments.
     407
     408        Input boundary_map is a dictionary of boundary objects keyed
     409        by symbolic tags to matched against tags in the internal dictionary
     410        self.boundary.
     411
     412        As result one pointer to a boundary object is stored for each vertex
     413        in the list self.boundary_objects.
     414        More entries may point to the same boundary object
     415
     416        Schematically the mapping is from two dictionaries to one list
     417        where the index is used as pointer to the boundary_values arrays
     418        within each quantity.
     419
     420        self.boundary:          (vol_id, edge_id): tag
     421        boundary_map (input):   tag: boundary_object
     422        ----------------------------------------------
     423        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
     424
     425
     426        Pre-condition:
     427          self.boundary has been built.
     428
     429        Post-condition:
     430          self.boundary_objects is built
     431
     432        If a tag from the domain doesn't appear in the input dictionary an
     433        exception is raised.
     434        However, if a tag is not used to the domain, no error is thrown.
     435        FIXME: This would lead to implementation of a
     436        default boundary condition
     437
     438        Note: If a segment is listed in the boundary dictionary and if it is
     439        not None, it *will* become a boundary -
     440        even if there is a neighbouring triangle.
     441        This would be the case for internal boundaries
     442
     443        Boundary objects that are None will be skipped.
     444
     445        FIXME: If set_boundary is called multiple times and if Boundary
     446        object is changed into None, the neighbour structure will not be
     447        restored!!!
     448        """
     449
     450        self.boundary_objects = []
     451       
     452
     453
     454
     455       
     456        self.boundary_map = boundary_map  #Store for use with eg. boundary_stats.
     457
     458        #FIXME: Try to remove the sorting and fix test_mesh.py
     459        x = self.boundary.keys()
     460        x.sort()
     461
     462        #Loop through edges that lie on the boundary and associate them with
     463        #callable boundary objects depending on their tags
     464        for k, (vol_id, edge_id) in enumerate(x):
     465            tag = self.boundary[ (vol_id, edge_id) ]
     466
     467            if boundary_map.has_key(tag):
     468                B = boundary_map[tag]  #Get callable boundary object
     469
     470                if B is not None:
     471                    self.boundary_objects.append( ((vol_id, edge_id), B) )
     472                    self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
     473                else:
     474                    pass
     475                    #FIXME: Check and perhaps fix neighbour structure
     476
     477            else:
     478                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
     479                msg += 'bound to a boundary object.\n'
     480                msg += 'All boundary tags defined in domain must appear '
     481                msg += 'in the supplied dictionary.\n'
     482                msg += 'The tags are: %s' %self.get_boundary_tags()
     483                raise msg
     484
     485
     486
     487    def check_integrity(self):
     488        #Mesh.check_integrity(self)
     489       
     490        for quantity in self.conserved_quantities:
     491            msg = 'Conserved quantities must be a subset of all quantities'
     492            assert quantity in self.quantities, msg
     493
     494        ##assert hasattr(self, 'boundary_objects')
     495
     496    def update_boundary(self):
     497        """Go through list of boundary objects and update boundary values
     498        for all conserved quantities on boundary.
     499        """
     500
     501        #FIXME: Update only those that change (if that can be worked out)
     502        #FIXME: Boundary objects should not include ghost nodes.
     503        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
     504            q = B.evaluate(vol_id, edge_id)
     505
     506            for j, name in enumerate(self.conserved_quantities):
     507                Q = self.quantities[name]
     508                Q.boundary_values[i] = q[j]
    81509
    82510
     
    101529    #D2 = Domain(points2)
    102530
    103 
    104 
    105 
    106 
    107 
Note: See TracChangeset for help on using the changeset viewer.