1

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.

Area of deformed hexagon

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.

NPMitchell
  • 193
  • 1
  • 10
  • Change with respect to what? Also: could you please provide us with an image? – knedlsepp Feb 20 '15 at 16:44
  • *summing squared area changes* vs. *I want to know the change in volume*, What is it you want to know? Maybe provide us with the underlying math also. – knedlsepp Feb 20 '15 at 16:51
  • @knedlsepp: I posted an image here: http://s1167.photobucket.com/user/npmitchell/media/Slide1_zpsyhqg9tvb.jpg.html – NPMitchell Feb 20 '15 at 19:03
  • @knedlsepp I want a 'volumetric strain' of a deformed lattice with respect to a reference lattice of equilateral triangles. I've updated the question to reflect this. – NPMitchell Feb 20 '15 at 19:05
  • @knedlsepp: Updated question to explain the context of the problem. Thank you! – NPMitchell Feb 20 '15 at 19:18
  • I suspect *[this](http://www.mathworks.com/matlabcentral/fileexchange/48609-dualmesh-polygonal-mesh-construction)* could be a starting point at least. From what I can tell at first sight, it includes a function `triadual2` that is said to be able to triangulate the dual mesh cells. Computing areas of triangles is easy, so you will simply need to use this function, find out how to reference the dual cells and compute the new triangles' areas. – knedlsepp Feb 20 '15 at 19:26
  • One thing: You say your triangles are in 3D, which makes this a bit more difficult. I imagine that the dual mesh should usually have planar hexagons. How would you define your deformed hexagon if the middle vertex was not in the same plane, as the others? Do you really have 3D data? The paper only mentions 2D... – knedlsepp Feb 21 '15 at 09:54
  • I have a weakly curved surface. I took your suggestion for the dualmesh toolbox, and it works very well and seems extremely efficient. For each surface, I check that the approximation by planar surfaces is sufficiently accurate, but this is straightforward to do, as the accuracy scales with the density (or the number of points in my lattice). I may experiment with accounting for non-planar deformations, but it looks like this isn't much of a problem if the unit cell width is sufficiently smaller than either radius of curvature of the surface. – NPMitchell Feb 22 '15 at 06:20
  • Ok. I just wanted to point out that the paper only mentions planar surfaces and the underlying mathematics may not be accurate for curved surfaces. Especially as at points of high curvature the surface of the planar dual voronoi cell will *de*crease using this computation, whereas it really should *in*crease if one would assume the cells to be on the triangular surface. – knedlsepp Feb 22 '15 at 11:12

1 Answers1

0

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.)

knedlsepp
  • 6,065
  • 3
  • 20
  • 41