Meshpy to Triangle switch

[normanrichardson]

I've seen a good amount of discussion about switching away from using Meshpy and using Triangle for meshing. This got me curious about how sectionproperties uses Meshpy and what this switch would entail. I saw this is planned to happen post v2 so I thought I would make a discussion for this switch. Here is what I have learnt so far:

1. Both Meshpy and Triangle reference the same C library, Quake, for tri meshing. I thought this was great as we should be able to confidently expect the exact same mesh results provided we get the setup correct. I would suggest going down this route with testing.
2. Quake takes its meshing arguments as switches.
3. The main geometric data structure of Meshpy and Triangle are similar. Meshpy has a MeshInfo object and Triangle uses a dictionary. Unfortunately, the naming convention of the attributes (in the Meshpy case) and keys (in the Triangle case) are different. In a sectionproperties context, here is a mapping of the important ones:

Meshpy MeshInfo	Triangle dictionary	sectionproperties
points	vertices	points
facets	segments	facets
holes	holes	holes
regions (column 0 and 1)	regions (column 0 and 1)	control_points
regions (column 2)	regions (column 2)	region_id
regions (column 3)	regions (column 3)	mesh_sizes
elements	triangles	elements
element_attributes	triangle_attributes	element_attributes

4. Ultimately, both Meshpy and Triangle use Quakes switches (from 2.). In Meshpy case it uses a series of arguments in meshpy.triangle.build, and it builds a switch string that it passes to Quake. In Triangle's case it takes the switch string directly as an argument in triangle.triangulate. Given this, I thought I would try find the equivalent calls between Meshpy and Triangle. Sectionproperties calls Mephpy at two locations with different mesh requirements:

sectionproperties	Meshpy Implimentation	Triangle Implementation	comments
pre.create_mesh	meshpy.triangle.build(mesh, min_angle=30, mesh_order=2, quality_meshing=True, attributes=True, volume_constraints=True)	triangle.triangulate(data, 'Qzpq30Aao2')	Tested here
PlasticSection.create_plastic_mesh	meshpy.triangle.build(mesh, mesh_order=2, quality_meshing=False, attributes=True)	triangle.triangulate(data, 'QzpAao2')	Not tested

Some things that come to mind are:

• The sectionproperties notation is taken from Meshpy, is this going to stay the same? To me I think this is a fair option.
• What is the best way to test the two libraries? Are gists sufficient?
• What other testing can be done to make sure we set traingle's opts correctly? I was thinking comparing meshes of some typical sections in sectionproperties.pre.sections, but this time just running these directly threw Meshpy and Triangle ie. pull the data points, facets etc. out, build meshInfo and a matching dictionary, run the function in 4, and comapre the meshes.

[robbievanleeuwen]

Thanks for the comprehensive rundown @normanrichardson! This will be super helpful when we make the switch!

[normanrichardson]

No worries, it's interesting stuff. I will also note that on the v2 branch PlasticSection.create_plastic_mesh has been updated/removed and no longer needs meshing. This makes sense to me as it looks like shapely can provide all the geometric operations and data needed for that root-finding process.

[connorferster]

Now that this has been merged into v2 on #153 (beautifully done, Robbie!) it will likely make installations easier for some folks.

For interest, here is a casual benchmark I ran between triangle and meshpy. There seems to be a slight improvement (~15%) after switching to triangle:
https://gist.github.com/connorferster/74c39da3da2710315921d95bc67eb03d