Ignore:
Timestamp:
May 6, 2008, 6:03:26 PM (16 years ago)
Author:
ole
Message:

Work on ticket:175 with special consideration to mesh intersections
Also reinstated disabled Okada test which currently fails but is faster.

File:
1 edited

Legend:

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

    r5276 r5286  
    1010from general_mesh import General_mesh
    1111from math import pi, sqrt
     12from Numeric import array, allclose
    1213       
    1314
     
    230231            if ratio <= min_ratio: min_ratio = ratio
    231232        return max_ratio,min_ratio
     233
     234
    232235
    233236    def build_neighbour_structure(self):
     
    628631
    629632        N = len(self)
    630         #Get x,y coordinates for all vertices for all triangles
     633
     634        # Get x,y coordinates for all vertices for all triangles
    631635        V = self.get_vertex_coordinates()
    632         #Check each triangle
     636
     637        # Check each triangle
    633638        for i in range(N):
    634639
     
    637642            x2, y2 = V[3*i+2, :]
    638643           
    639             #Check that area hasn't been compromised
     644            # Check that area hasn't been compromised
    640645            area = self.areas[i]
    641646            ref = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
     
    644649            assert abs((area - ref)/area) < epsilon, msg
    645650
    646             #Check that points are arranged in counter clock-wise order
     651            # Check that points are arranged in counter clock-wise order
    647652            v0 = [x1-x0, y1-y0]
    648653            v1 = [x2-x1, y2-y1]
     
    656661            assert a0 < pi and a1 < pi and a2 < pi, msg
    657662
    658             #Check that normals are orthogonal to edge vectors
    659             #Note that normal[k] lies opposite vertex k
     663            # Check that normals are orthogonal to edge vectors
     664            # Note that normal[k] lies opposite vertex k
    660665
    661666            normal0 = self.normals[i, 0:2]
     
    665670            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
    666671
    667                 #Normalise
     672                # Normalise
    668673                l_u = sqrt(u[0]*u[0] + u[1]*u[1])
    669674                l_v = sqrt(v[0]*v[0] + v[1]*v[1])               
    670675
    671                 x = (u[0]*v[0] + u[1]*v[1])/l_u/l_v #Inner product
     676                x = (u[0]*v[0] + u[1]*v[1])/l_u/l_v # Inner product
    672677               
    673678                msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
     
    896901
    897902
     903    def _get_intersecting_segments(self, line, triangle_intersections={}):
     904      """Find edges intersected by line
     905
     906      Input:
     907          line - list of two points forming a segmented line
     908          triangle_intersections to be updated
     909      Output:
     910          list of instances of class Triangle_intersection
     911
     912      This method is used by the public method
     913      get_intersecting_segments(self, polyline) which also contains
     914      more documentation.
     915      """
     916
     917      from anuga.utilities.polygon import intersection
     918     
     919      msg = 'Line segment must contain exactly two points'
     920      assert len(line) == 2, msg
     921     
     922     
     923      # Check intersection with edge segments for all triangles
     924      # FIXME (Ole): This should be implemented in C
     925      V = self.get_vertex_coordinates()
     926      N = len(self)
     927
     928      for i in range(N):
     929          # Get nodes and edge segments for each triangle
     930          x0, y0 = V[3*i, :]
     931          x1, y1 = V[3*i+1, :]
     932          x2, y2 = V[3*i+2, :]
     933           
     934
     935          edge_segments = [[[x0,y0], [x1, y1]],
     936                            [[x1,y1], [x2, y2]],
     937                            [[x2,y2], [x0, y0]]]
     938
     939          # Find segments that are intersected by line
     940          intersections = []
     941          for edge in edge_segments:
     942
     943              status, value = intersection(line, edge)
     944              if status == 1:
     945                  # Normal intersection of one edge
     946
     947                  # Exclude singular intersections with vertices
     948                  if not(allclose(value, edge[0]) or\
     949                         allclose(value, edge[1])):
     950                      intersections.append(value)
     951
     952              if status == 2:
     953                  # Edge is sharing a segment with line
     954
     955                  # Add identified line segment
     956                  for j in range(value.shape[0]):
     957                      intersections.append(value[j,:])
     958
     959
     960          msg = 'There can be only two or no intersections'
     961          assert len(intersections) in [0,2], msg
     962
     963
     964          if len(intersections) == 2:
     965
     966              # Calculate attributes for this segment
     967
     968              # Origin of intersecting line to be used for direction
     969              xi0 = line[0][0]
     970              eta0 = line[0][1]
     971
     972              # End points of intersecting segment
     973              x0, y0 = intersections[0]
     974              x1, y1 = intersections[1]
     975
     976
     977              # Determine which end point is closer to the origin of the line
     978              # This is necessary for determining the direction of
     979              # the line and the normals
     980
     981              # Distances from line origin to the two intersections
     982              z0 = array([x0 - xi0, y0 - eta0])
     983              z1 = array([x1 - xi0, y1 - eta0])             
     984              d0 = sqrt(sum(z0**2))
     985              d1 = sqrt(sum(z1**2))
     986                 
     987              if d1 < d0:
     988                  # Swap
     989                  xi, eta = x0, y0
     990                  x0, y0 = x1, y1
     991                  x1, y1 = xi, eta
     992
     993              # (x0,y0) is now the origin of the intersecting segment
     994                 
     995
     996              # Normal direction (right hand side relative to direction of line)
     997              vector = array([x1 - x0, y1 - y0]) # Segment vector
     998              length = sqrt(sum(vector**2))      # Segment length
     999              normal = array([vector[1], -vector[0]])/length
     1000
     1001
     1002              segment = ((x0,y0), (x1, y1))   
     1003              T = Triangle_intersection(segment=segment,
     1004                                        normal=normal,
     1005                                        length=length,
     1006                                        triangle_id=i)
     1007
     1008
     1009              # Add segment unless it was done earlier
     1010              if not triangle_intersections.has_key(segment):   
     1011                  triangle_intersections[segment] = T
     1012
     1013
     1014      #return triangle_intersections
     1015
     1016
     1017
    8981018    def get_intersecting_segments(self, polyline):
    8991019      """Find edges intersected by polyline
     
    9031023
    9041024      Output:
    905           dictionary where
    906               keys are triangle ids for those elements that are intersected
    907               values are a list of segments each represented as a triplet:
    908                   length of the intersecting segment
    909                   right hand normal
    910                   midpoint of intersecting segment
    911 
    912       The polyline may break inside any triangle causing multiple segments per triabgle.
     1025          list of instances of class Triangle_intersection
     1026
     1027      The polyline may break inside any triangle causing multiple
     1028      segments per triangle - consequently the same triangle may
     1029      appear in several entries.
     1030
     1031      If a polyline segment coincides with a triangle edge,
     1032      the the entire shared segment will be used.
     1033
     1034      Onle one of the triangles thus intersected will be used and that
     1035      is the first one encoutered.
     1036
     1037      intersections with single vertices are ignored.
     1038
     1039      Resulting segments are unsorted
    9131040      """
    9141041
    915       pass
     1042      msg = 'Polyline must contain at least two points'
     1043      assert len(polyline) >= 2, msg
     1044
     1045      # For all segments in polyline
     1046      triangle_intersections = {}
     1047      for i, point0 in enumerate(polyline[:-1]):
     1048          point1 = polyline[i+1]
     1049         
     1050          line = [point0, point1]
     1051
     1052          self._get_intersecting_segments(line, triangle_intersections)
     1053
     1054
     1055      return triangle_intersections.values()
     1056 
    9161057
    9171058 
     
    9311072            return []
    9321073       
     1074
     1075
     1076class Triangle_intersection:
     1077    """Store information about line segments intersecting a triangle
     1078   
     1079    Attributes are
     1080
     1081        segment: Line segment intersecting triangle [[x0,y0], [x1, y1]]
     1082        normal: [a,b] right hand normal to segment
     1083        length: Length of intersecting segment
     1084        triangle_id: id (in mesh) of triangle being intersected
     1085
     1086    """
     1087
     1088
     1089    def __init__(self,
     1090                 segment=None,
     1091                 normal=None,
     1092                 length=None,
     1093                 triangle_id=None):
     1094        self.segment = segment             
     1095        self.normal = normal
     1096        self.length = length
     1097        self.triangle_id = triangle_id
     1098       
     1099
     1100    def __repr__(self):
     1101        s = 'Triangle_intersection('
     1102        s += 'segment=%s, normal=%s, length=%s, triangle_id=%s)'\
     1103             %(self.segment,
     1104               self.normal,
     1105               self.length,
     1106               self.triangle_id)
     1107   
     1108        return s
Note: See TracChangeset for help on using the changeset viewer.