property_topomesh_optimization.py 文件源码

def property_topomesh_area_smoothing_force(topomesh,target_areas=None):

    compute_topomesh_property(topomesh,'vertices',degree=2)
    compute_topomesh_property(topomesh,'length',degree=1)

    triangle_vertices = topomesh.wisp_property('vertices',degree=2).values()
    rotated_triangle_vertices = np.transpose([triangle_vertices[:,2],triangle_vertices[:,0],triangle_vertices[:,1]])
    antirotated_triangle_vertices = np.transpose([triangle_vertices[:,1],triangle_vertices[:,2],triangle_vertices[:,0]])
    triangle_vertices = np.append(np.append(triangle_vertices,rotated_triangle_vertices,axis=0),antirotated_triangle_vertices,axis=0)

    edge_index_list = np.array([[1, 2],[0, 1],[0, 2]])
    triangle_edge_vertices = triangle_vertices[:,edge_index_list]

    triangle_edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,1]) - topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,0])
    #triangle_edge_lengths = np.power(np.sum(np.power(triangle_edge_vectors,2.0),axis=2),0.5)
    triangle_edge_lengths = np.linalg.norm(triangle_edge_vectors,axis=2)
    triangle_edge_directions = triangle_edge_vectors/triangle_edge_lengths[...,np.newaxis]

    triangle_perimeters = np.sum(triangle_edge_lengths,axis=1)
    triangle_areas = np.sqrt((triangle_perimeters/2.0)*(triangle_perimeters/2.0-triangle_edge_lengths[:,0])*(triangle_perimeters/2.0-triangle_edge_lengths[:,1])*(triangle_perimeters/2.0-triangle_edge_lengths[:,2]))
    triangle_areas[np.where(triangle_areas==0.0)] = 0.001

    if target_areas is None:
        target_areas = triangle_areas.mean()*np.ones(len(triangle_areas))
    else:
        target_areas = np.tile(target_areas,(3))

    area_edge_1_force = (target_areas-triangle_areas)*(triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,1]**2)*triangle_edge_lengths[:,1]
    area_edge_2_force = (target_areas-triangle_areas)*(triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2)*triangle_edge_lengths[:,2]

    area_unitary_force = -(area_edge_1_force[:,np.newaxis]*triangle_edge_directions[:,1] + area_edge_2_force[:,np.newaxis]*triangle_edge_directions[:,2])

    triangle_force = np.transpose([nd.sum(area_unitary_force[:,d],triangle_vertices[:,0],index=list(topomesh.wisps(0))) for d in xrange(3)])

    return triangle_force