How are types defined in input files and scripts?

[IsaacPincus]

Hi!

I've been trying to mess around with adding local area constraints to a membrane simulation through potentialConstantAreaTriangle, and I've become a bit confused about how vertex types are defined. In the first example, the faces_N*.inp and vertices_N*.inp files denote every vertex and face with type 1, which seems to be read in via __add_vertex(...), meaning that each vertex will have vertex.type = 1.

However, the scripts seem to have the forces added to type 0 vertices, for example in disclination.py:
evolver.add_force("Mesh>Harmonic", {"k":{"0":k}, "l0":{"0":l0}})

When I step through potentialHarmonic.cpp with a debugger, the type of the edges does indeed seem to be 0 during runtime, but then I'm not sure exactly how this relates to the properties read in from the input files. When I tried the same with potentialConstantAreaTriangle.cpp, there the face type seems to be 1 despite apparently defining things in a similar way:
evolver.add_force("Mesh>Constant Area", { "sigma": {"0": str(1e2)}, "target_area": {"0": str(target_area_per_face)} })
and at this point I'm fairly confused and not sure if I'm making simple mistakes.

Is there some easier way to understand how types are defined in pyMembrane? Happy to provide some minimal working examples if that makes things clearer! Thanks in advance for any help.

[fdmatoz]

Dear Isaac,

Thank you for your message. Types are defined in each object, i.e, vertices, edges and faces. When you load the mesh, I set-up the types to a given value. You can check them like so,

Example-Issac.zip

Now, Mesh>Constant Area uses the faces as such:

for (int face_index = 0; face_index < Numfaces; face_index++)
{
    int type = faces[face_index].type;
    int v1 = faces[face_index].v1;
    int v2 = faces[face_index].v2;
    int v3 = faces[face_index].v3;

    real g_now[3];
    pymemb::compute_form_factor_triangle(g_now, vertices[v1].r, vertices[v2].r, vertices[v3].r, _box);
    real energy = ComputeVertexConstantAreaTriangleEnergy_lambda(g_now, _sigma[type], _target_area[type]);

    /// Add energy to that face
    faces[face_index].energy += energy;
}

The reason for which I didn't want to read all of the information from the file is that people might want to do it by themselves but perhaps I should been more clear on the tutorials.

Please let me know..

Regards

D.

[IsaacPincus]

Thanks Daniel, this is very helpful! So if I'm understanding correctly, the types of the faces and vertices are '1', but the types of the edges are still '0' since we haven't explicitly set them in the input files? Then, the 'Mesh>Bending>Dihedral' potential uses the edges for computation of energies, so we set it for type '0', while the 'Mesh>Constant Area' potential uses the faces, so we feed in type '1'? As in:

# bending potential
kappa = str(1.0)
evolver.add_force("Mesh>Bending>Dihedral", {"kappa":{"0":kappa}})

# constant area potential
initial_area = compute.area()
target_area_per_face = initial_area/len(system.faces)
kD = 1e5
print(str(target_area_per_face))
evolver.add_force("Mesh>Constant Area", {
    "sigma": {"1": str(kD)},
    "target_area": {"1": str(target_area_per_face)}
})