source: development/momentum_sink/rotate_buildings.py @ 2379

import Numeric
import math
# add inundation dir to your pythonpath
from pmesh.mesh import Mesh
from coordinate_transforms.geo_reference import Geo_reference


def create_mesh(maximum_triangle_area,
                mesh_file=None,
                triangles_in_name = False):
    """
    triangles_in_name, if True is used to append the number of
    triangles in the mesh to the mesh file name.
    """
    # create a mesh instance of class Mesh
    m = Mesh()

    # Boundary of problem
    WidtH = 100 # width of boudary in metres
    #W = WidtH/8
    #L = W
    outer_polygon = [[0,0],[1.5*WidtH,0],[1.5*WidtH,WidtH],[0,WidtH]]
    print outer_polygon
    m.add_region_from_polygon(outer_polygon, tags={'wall':[0,1,2], 'wave':[3]})

    # inner polygons => building boundaries
    depth = 15 # depth of building side to oncoming wave
    whs = depth/2 
    breadth = 20 # breadth of building
    lhs = breadth/2 

    # Sample rotation Matrix
    Thd = 30 # Degrees
    Th = Thd * 3.14259/180 # Deg to Rad.

    wh1 = (-whs) * math.cos(Th) + (-lhs) * math.sin(Th)
    lh1 = (-lhs) * math.cos(Th) - (-whs) * math.sin(Th)

    wh2 = (+whs) * math.cos(Th) + (-lhs) * math.sin(Th)
    lh2 = (-lhs) * math.cos(Th) - (+whs) * math.sin(Th)

    wh3 = (+whs) * math.cos(Th) + (+lhs) * math.sin(Th)
    lh3 = (+lhs) * math.cos(Th) - (+whs) * math.sin(Th)

    wh4 = (-whs) * math.cos(Th) + (+lhs) * math.sin(Th)
    lh4 = (+lhs) * math.cos(Th) - (-whs) * math.sin(Th)

    print "building footprint"
    print depth * breadth , "m^2"
    block = 625
    BL = block**0.5
    ForDep = (0.2*WidtH) + (BL/2)
    RearDep = 1.2*WidtH 
    porosity = breadth/BL
    print porosity, " Building porosity"
    
    Breadths = Numeric.arrayrange( (BL/2), WidtH, (BL))
    print Breadths, "Breadths"
    Depths = Numeric.arrayrange( ForDep, RearDep, BL )
    print Depths, "Depths"

    for D in Depths:
        #Breadths = Breadths + BL/2
        for B in Breadths:
            polygon = [[D+wh1,B+lh1],[D+wh2,B+lh2],[D+wh3,B+lh3],[D+wh4,B+lh4]]
            m.add_hole_from_polygon(polygon, tags={'wall':[0,1,2,3]})# Adds holes with reflective boundaries.        
        

            #print polygon
    m.generate_mesh(maximum_triangle_area=maximum_triangle_area)
    triangle_count = m.get_triangle_count()
    
    if mesh_file is None:    
        return m, triangle_count
    else:
        if triangles_in_name is True:
            mesh_file = mesh_file[:-4] + '_' + str(triangle_count) \
                        + mesh_file[-4:]
        m.export_mesh_file(mesh_file)
        return mesh_file, triangle_count

#-------------------------------------------------------------
if __name__ == "__main__":
    _, triangle_count = create_mesh(10,mesh_file="test.tsh")
    print "triangle_count",triangle_count