Function main

scripts/geodesic.py:121–166

()

    return div_X;

def main():
    args = parse_args();
    mesh = pymesh.load_mesh(args.input_mesh);
    assert(mesh.dim == 3);
    assert(mesh.vertex_per_face == 3);

    compute_edge_field(mesh);
    assert(mesh.has_attribute("edges"));
    compute_cotan_field(mesh);
    assert(mesh.has_attribute("cotan"));
    compute_area_and_normal_field(mesh);
    assert(mesh.has_attribute("vertex_area"));
    assert(mesh.has_attribute("face_normal"));

    assembler = pymesh.Assembler(mesh);
    L = assembler.assemble("laplacian") * -1;
    A = assembler.assemble("mass");
    G = assembler.assemble("gradient");

    t = np.mean(mesh.get_attribute("face_area").ravel());
    rhs = np.zeros(mesh.num_vertices);
    rhs[args.source] = 1.0;

    u = scipy.sparse.linalg.spsolve(A - t*L, rhs);
    mesh.add_attribute("u");
    mesh.set_attribute("u", u);

    grad_u = (G * u).reshape((-1, 3), order="C");
    grad_u_2 = compute_gradient(mesh, u);
    mesh.add_attribute("grad_u");
    mesh.set_attribute("grad_u", grad_u.ravel());

    X = -grad_u / norm(grad_u, axis=1)[...,np.newaxis];
    mesh.add_attribute("X");
    mesh.set_attribute("X", X.ravel());

    div_X = compute_divergence(mesh, X);
    mesh.add_attribute("div_X");
    mesh.set_attribute("div_X", div_X);

    phi = scipy.sparse.linalg.spsolve(L, div_X);
    phi = phi - np.amin(phi);
    mesh.add_attribute("phi");
    mesh.set_attribute("phi", phi);

    pymesh.save_mesh(args.output_mesh, mesh, "u", "grad_u", "X", "div_X", "phi");


if __name__ == "__main__":