# Allow user-provided vertex/edge areas in CollisionMesh (e.g. for axisymmetric models) #266

Motivation

We use ipctk (1.6.0, Python bindings) in fedoo, a finite element library, including for 2D axisymmetric models. There, the collision mesh is the planar (r, z) outline of a solid of revolution, and each planar collision stands for a ring of contact. The area a vertex represents is therefore 2*pi*R * ds rather than the planar barycentric length ds computed by CollisionMesh.

The same need arises whenever the quadrature weight of a surface vertex is not its barycentric area in the collision mesh: plane strain/stress with a thickness, periodic or symmetric sectors, shells with a custom lumping, etc.

Current workaround

After each NormalCollisions.build, we rescale the weight of every collision from Python:

collisions.build(mesh, vertices, dhat, broad_phase=broad_phase)
edges, faces = mesh.edges, mesh.faces
for i in range(len(collisions)):
    c = collisions[i]
    ids = c.vertex_ids(edges, faces)
    n = c.num_vertices()
    R = radius[ids[0]] if n == 3 else radius[ids[:n]].mean()
    c.weight *= 2 * np.pi * R

This gives a consistent potential, gradient and hessian, and the tangential collisions inherit the weight. It has three drawbacks:

1. Cost: a Python loop over all collisions at every Newton iteration and line-search trial (about 16 µs per collision; 3 s for 170k collisions, versus 0.04 s for the build itself).
2. OGC: ogc.TrustRegion.update_if_needed rebuilds the collision set internally, so the rescaled weights are lost and we cannot use OGC in this setting.
3. Merged terms: with use_area_weighting and IMPROVED_MAX_APPROX, a vertex-vertex term merges contributions from both vertices (-(A_p/2 + A_v/2)). From the merged weight alone, the exact value -(w_p*A_p/2 + w_v*A_v/2) cannot be recovered, so we use a mean radius.

Proposal

Let the user provide (or scale) the vertex and edge areas of a CollisionMesh, for instance:

mesh.vertex_areas = 2 * np.pi * R * np.asarray(mesh.vertex_areas)   # setter
# or
mesh.set_vertex_areas(areas); mesh.set_edge_areas(areas)

vertex_areas / edge_areas are currently read-only (AttributeError: property of 'CollisionMesh' object has no setter). Since NormalCollisions.build already takes the weights from these areas when use_area_weighting is on, a setter would solve the three points above without any axisymmetric-specific code in ipctk.

Related question

initial_barrier_stiffness computes its lower bound from average_mass, bbox_diagonal and the barrier only. When collision weights are not of order one (area weighting, or the scaling above), the effective stiffness kappa * weight of that lower bound depends on the weights and on the unit system. Is scaling average_mass by 1 / mean(weight) the intended way to handle this, or would a weight-aware variant make sense?

I am happy to open a PR for the setters if this direction suits you.