Summing squared area changes from Voronoi cells given area of triangles in 3D?

Question

I have a list of triangles in 3D that form a surface (ie a triangulation). The structure is a deformed triangular lattice. I want to know the change in area of the deformed hexagons of the voronoi tessalation of the lattice with respect to the rest area of the undeformed lattice cells (ie with respect to a regular hexagon). In fact, I really want the sum of the squared change in area of the hexagonal unit cells associated with those triangles.

Background/Math details: I'm approximating a curved elastic sheet by a triangular lattice. One way to tune the poisson ratio (elastic constant) of the sheet is by adding a 'volumetric' strain energy term to the energy. I'm trying to compute a 'volumetric' strain energy of a deformed, elastic, triangular lattice, defined as: U_volumetric = 1/2 T (e_v)^2, where e_v=deltaV/V is determined by the change in area of a voronoi cell with respect to its reference area, which is a known constant.

Reference: https://www.researchgate.net/publication/265853755_Finite_element_implementation_of_a_non-local_particle_method_for_elasticity_and_fracture_analysis

Want:

Sum[ (DeltaA/ A).^2 ] over all hexagonal cells.

My data is stored in the variables:

xyz = [ x1,y1,z1; x2,y2,z2; etc] %the vertices/particles in 3D

TRI = [ vertex0, vertex1, vertex2; etc] % where vertex0 is the row of xyz for the particle sitting at vertex 0 of the first triangle.

NeighborList = [ p1n1, p1n2, p1n3, p1n4, p1n5,p1n6 ; p2n1...] % where p1n1 is particle 1's first nearest neighbor as a row index for xyz. For example, xyz(NL(1,1),:) returns the xyz location of particle 1's first neighbor.

AreaTRI = [ areaTRI1; areaTRI2; etc]

I am writing this in MATLAB.

As of now, I am approximating the amount of area attributed to each vertex as 1/3 of the triangle's area, then summing over the 6 nearest neighbor triangles. But a voronoi cell area will NOT be exactly equal to Sum_(i=0,1,...5) 1/3* areaTRI_i, so this is a bad approximation. See the image in the link above, which I think makes this clearer.

Answer 1

You can do this using the DUALMESH-submission on the file exchange:

DUALMESH is a toolbox of mesh processing routines that allow the construction of "dual" meshes based on underlying simplicial triangulations. Support is provided for various planar and surface triangulation types, including non-Delaunay and non-manifold types.

Simply use the following commands to generate a vector areas of all the dual elements' areas. The ordering will correspond to the nodes xyz.

[cp,ce,pv,ev] = makedual2(xyz, TRI);
[~,areas(cp(:,1))] = geomdual2(cp,ce,pv,ev);

You might want to have a look at the boundary areas using:

trisurf(TRI, xyz(:,1), xyz(:,2), areas);

The dual cells of boundary nodes theoretically are unbounded and thus should have infinite area. This submission handles it differently however: Instead of an unbounded cell it will return the intersection of the unbounded cell with the original mesh.

Also mind that your question is not well defined if the mesh you are working with is not planar, as the dual mesh cells will be planar and won't scale the same way as the triangles. So this solution will probably only work correctly if your mesh is really 2D. (From what I can tell, the paper you mention is also only for the 2D-case.)