class Point:
    def __init__(self, x=None, y=None):
        self.x = x or 0.0
        self.y = y or 0.0

    def __str__(self):
        return '(' + str(self.x) + ', ' + str(self.y) + ')'

    def __hash__(self):
        return self.__str__().__hash__()

    def __eq__(self, other):
        return self.__str__() == other.__str__()


class Edge:
    def __init__(self, a=None, b=None):
        self.a = a or Point()
        self.b = b or Point()

    def __str__(self):
        return '[' + str(self.a) + '; ' + str(self.b) + ']'


def _point_on_segment(edge, p):
    if min(edge.a.x, edge.b.x) <= p.x <= max(edge.a.x, edge.b.x) and min(edge.a.y, edge.b.y) <= p.y <= max(edge.a.y, edge.b.y):
        return True
    return False

def _intersection_point(E1, E2):
    P1 = E1.a
    P2 = E1.b
    P3 = E2.a
    P4 = E2.b

    A1 = (P1.y - P2.y) / (P1.x - P2.x) if P1.x != P2.x else P1.x
    A2 = (P3.y - P4.y) / (P3.x - P4.x) if P3.x != P4.x else P3.x
    B1 = P1.y - A1 * P1.x
    B2 = P3.y - A2 * P3.x

    # parallel segments
    if A1 == A2:
        return True

    Xa = (B2 - B1) / (A1 - A2)
    Ya = A1 * Xa + B1

    # no intersections
    if Xa < max(min(P1.x, P2.x), min(P3.x, P4.x)) or Xa > min(max(P1.x, P2.x), max(P3.x, P4.x)):
        return False

    return Point(Xa, Ya)


class Shape:
    def __init__(self, edges=None, vertices=None):
        self.edges = edges or []
        self.vertices = vertices or []

    def area(self):
        S = 0.0
        for i, j in zip(self.vertices, self.vertices[1:]) + [(self.vertices[-1], self.vertices[0])]:
            S += (i.x + j.x) * (i.y - j.y)
        return 0.5 * abs(S)

    def point_in_polygon(self, A):
        N = len(self.vertices)
        inside = False

        P1 = self.vertices[0]
        for i in xrange(N+1):
            P2 = self.vertices[i % N]
            if min(P1.y, P2.y) < A.y <= max(P1.y, P2.y):
                if A.x <= max(P1.x, P2.x):
                    if P1.y != P2.y:
                        xints = P1.x + (A.y - P1.y)*(P2.x - P1.x)/(P2.y - P1.y)
                    if P1.x == P2.x or A.x <= xints:
                        inside = not inside
            P1 = P2

        return inside

    def __str__(self):
        return '{ ' + ' | '.join([str(edge) for edge in self.edges]) + ' }'


def _intersect_poly(poly1, poly2):
    points = {}
    res = Shape()
    extra = {}

    for edge1 in poly1.edges:
        points[edge1] = []

        for edge2 in poly2.edges:
            new_vertex = _intersection_point(edge1, edge2)
            if new_vertex not in (True, False):
                points[edge1].append(new_vertex)
                if edge2 in extra:
                    extra[edge2].append(new_vertex)
                else:
                    extra[edge2] = [new_vertex]

        # Edge doesn't cross
        if len(points[edge1]) == 0:
            continue

        elif len(points[edge1]) == 1:
            a = edge1.a
            b = edge1.b

            res.vertices.append(points[edge1][0])
            if poly2.point_in_polygon(a):
                res.edges.append(Edge(a, points[edge1][0]))
                res.vertices.append(a)
            else:
                res.edges.append(Edge(b, points[edge1][0]))
                res.vertices.append(b)

        elif len(points[edge1]) == 2:
            a = points[edge1][0]
            b = points[edge1][1]

            res.edges.append(Edge(a, b))
            res.vertices.append(a)
            res.vertices.append(b)

    for k, v in extra.iteritems():
        if len(v) == 2:
            res.edges.append(Edge(v[0], v[1]))

    res.vertices = list(set(res.vertices))
    return res if len(res.vertices) != 0 else None


# TRIANGLES
print('Getting info from 515_flap_verh.bdf')

triangle_grid_file = open("515_flap_verh.bdf", "r")
triangle_grid_strings = triangle_grid_file.readlines()
triangle_grid_file.close()

# Nodes coords
nodes_strings = filter(lambda s: s.startswith('GRID*'), triangle_grid_strings)

triangle_grid_nodes = {}
for s in nodes_strings:
    temp_str = s.split()
    N = int(temp_str[1])
    x = float(temp_str[3])
    y = float(temp_str[4])

    triangle_grid_nodes[N] = Point(x, y)

# Triangles
elems_strings = filter(lambda s: s.startswith('CTRIA3'), triangle_grid_strings)

triangle_grid_elements = {}
for s in elems_strings:
    temp_str = s.split()
    N = int(temp_str[1])
    v1 = int(temp_str[3])
    v2 = int(temp_str[4])
    v3 = int(temp_str[5])

    triangle_grid_elements[N] = (v1, v2, v3)


# Pressures
press_strings = filter(lambda s: s.startswith('PLOAD2*'), triangle_grid_strings)

triangle_grid_pressures = {}
for s in press_strings:
    temp_str = s.split()
    N = int(temp_str[3])
    P = float(temp_str[2])

    triangle_grid_pressures[N] = P


# QUADS
print('Getting info from 000_Flaperon.bdf')

quad_grid_file = open("000_Flaperon.bdf", "r")
quad_grid_strings = quad_grid_file.readlines()
quad_grid_file.close()

# Nodes
nodes_strings = filter(lambda s: s.startswith('GRID'), quad_grid_strings)

quad_grid_nodes = {}
for s in nodes_strings:
    temp_str = s.split()
    N = int(temp_str[1])
    x = float(temp_str[2]) / 1000.0
    y = float(temp_str[3]) / 1000.0

    quad_grid_nodes[N] = Point(x, y)


# Elems
elems_strings = filter(lambda s: s.startswith('CQUAD4'), quad_grid_strings)

quad_grid_elements = {}
for s in elems_strings:
    temp_str = s.split()
    N = int(temp_str[1])
    v1 = int(temp_str[3])
    v2 = int(temp_str[4])
    v3 = int(temp_str[5])
    v4 = int(temp_str[6])

    quad_grid_elements[N] = (v1, v2, v3, v4)


# BUILD SHAPES
print('Building triangles...')

TRIANGLES = {}
for k, v in triangle_grid_elements.iteritems():
    v1 = triangle_grid_nodes[v[0]]
    v2 = triangle_grid_nodes[v[1]]
    v3 = triangle_grid_nodes[v[2]]

    e1 = Edge(v1, v2)
    e2 = Edge(v2, v3)
    e3 = Edge(v1, v3)

    S = Shape(edges=[e1,e2,e3], vertices=[v1,v2,v3])
    area = S.area()
    pressure = triangle_grid_pressures[k]

    TRIANGLES[k] = (S, area, pressure)

# BUILD QUADS
print('Building quads...')

QUADS = {}
for k, v in quad_grid_elements.iteritems():
    v1 = quad_grid_nodes[v[0]]
    v2 = quad_grid_nodes[v[1]]
    v3 = quad_grid_nodes[v[2]]
    v4 = quad_grid_nodes[v[3]]

    e1 = Edge(v1, v2)
    e2 = Edge(v2, v3)
    e3 = Edge(v3, v4)
    e4 = Edge(v1, v4)

    S = Shape(edges=[e1,e2,e3,e4], vertices=[v1,v2,v3,v4])

    QUADS[k] = S


# MAIN LOOP
print('Calculation...')

res_pressures = {}

for triangle_n, triangle_info in TRIANGLES.iteritems():
    triangle = triangle_info[0]
    triangle_area = triangle_info[1]
    triangle_pressure = triangle_info[2]

    for quad_n, quad in QUADS.iteritems():
        intersection = _intersect_poly(triangle, quad)
        if intersection:
            areas_ratio = intersection.area() / triangle_area
            pressure = triangle_pressure * areas_ratio
            if quad_n in res_pressures:
                res_pressures[quad_n].append(pressure)
            else:
                res_pressures[quad_n] = [pressure]


print('Done!')

res_file = open('res.txt', 'w')
RES = {}
for k, v in res_pressures.iteritems():
    RES[k] = sum(v) / len(v)
    res_file.write(str(k) + ': ' + str(RES[k]) + '\n')

res_file.close()