Shapely preprocessor: For discussion, not to merge

[connorferster]

Hi @robbievanleeuwen,

I have begun the "overhaul" of the pre-processor to run on shapely. I think it is working quite beautifully! I have not run any of the tests but have just been playing with the sections and plotting them to review them. In the process, I have made some design changes that I wanted to discuss with you, and get your thoughts on, before I proceed further:

1. Geometry class has a new attribute, geom, which contains the shapely object (Polygon)
2. Instantiating a new Geometry object now only requires a shapely object and no other parameters
3. Most of the Geometry methods have been re-implemented, relying on shapely's own methods. Some methods are no longer required (e.g. Geometry.calculate_facet_length(...) because shapely can produce the length of the perimeter from Polygon.length).
4. .control_points and .holes are now automatically generated, with no input required from the user, using Polygon.representative_point() which returns a "cheaply computed point that is guaranteed to be within the shape". Holes in Polygon instances are stored as LinearRing instances within Polygon.interiors. Converting these LinearRing instances to Polygons allows me to use the .representative_point() method to also generate the "hole" control points.
5. The addition of a new method (with some helper functions outside of the Geometry class) called create_facets_and_control_points(...). This method converts the Geometry.geom shapely object from "shapely format" into "sectionproperties format". This method is run whenever either Geometry.create_mesh(...) or Geometry.plot_geometry(...) is called. This new method populates self.points, self.facets, self.control_points, self.holes, and self.perimeter. Once these attributes are populated, then .create_mesh(...) and .plot_geometry(...) work exactly as normal.
6. All section geometries are now instances of Geometry class instead of instances of one of the several custom classes that inherit from Geometry and are essential a Geometry object but with a custom constructor. The custom constructors are now standalone functions e.g. no longer a Rhs class but a rectangular_hollow_section(...) function which returns a new Geometry instance.

Based on my work so far, here are my thoughts on how I might proceed and how things might roll out. Please let me know your thoughts on these ideas, also:

• Reduce the amount of custom classes in the preprocessor down to just Geometry, Material, and Section (currently CrossSection). Section is composed of one or more Geometry instances and one or more Material instances.
• Geometry "constructor" functions simply return a new Geometry instance (with appropriate shape) instead of returning the instance and "shift"-ing them. Translation, rotation, and mirroring are methods that only exist within Geometry class.
• Move Geometry.create_mesh(...) into Section.create_mesh(...). It seems that meshing occurs immediately before analysis and analysis only occurs on complete sections (i.e. geometry + material).
• Remove the .add_point(...), .add_facet(...), .add_hole(...), etc. methods from Geometry. The idea being that the user either, a) creates sections entirely with the constructor functions; or b) the user manipulates the shapely object(s) (within Geometry.geom) directly until they are happy with it. It would be awkward to add points (etc.) on-the-fly to shapely objects when it is likely easier to create a list of vertices all at once and feed that list into shapely.

Things have been going really smoothly on this project for me. I am eager to continue working on it and getting this implemented and tested in the next few weeks.

[connorferster]

See the below gist for a quick demo:

https://gist.github.com/connorferster/f009671157b3108a9f6494f89274c5ab

[robbievanleeuwen]

Overall this is looking great! Love your work 👍 👍 👍 Very happy with the direction this is going - I think this will unlock a lot of potential and make a lot of the code simpler and potentially faster!

I'd be interested to see how the implementation works with the CustomSection - it would be good to maintain a method to implement a cross-section with a list of points, facets and holes and a bonus to have the control points and perimeter generated by shapely! I think this will also tie nicely into the work that is currently being done on the dxf import which is currently essentially just populating points, facets & holes, but will also need to somehow accomodate multiple regions.

Responses to your talking points:

• Reduce the amount of custom classes in the preprocessor down to just Geometry, Material, and Section (currently CrossSection). Section is composed of one or more Geometry instances and one or more Material instances.

Sounds like a good idea. I imagine this will be a bit of a rework with the Section class and managing multiple instances of Geometry and Material but will be a solid approach.

• Geometry "constructor" functions simply return a new Geometry instance (with appropriate shape) instead of returning the instance and "shift"-ing them. Translation, rotation, and mirroring are methods that only exist within Geometry class.

👍

• Move Geometry.create_mesh(...) into Section.create_mesh(...). It seems that meshing occurs immediately before analysis and analysis only occurs on complete sections (i.e. geometry + material).

Good idea. The only thing to keep in mind is that meshing is used in an algorithm for the plastic analysis but I don't see this as being a problem as long as it's thought about.

• Remove the .add_point(...), .add_facet(...), .add_hole(...), etc. methods from Geometry. The idea being that the user either, a) creates sections entirely with the constructor functions; or b) the user manipulates the shapely object(s) (within Geometry.geom) directly until they are happy with it. It would be awkward to add points (etc.) on-the-fly to shapely objects when it is likely easier to create a list of vertices all at once and feed that list into shapely.

Happy for these to be removed. These no longer serve a function.

[connorferster]

Thanks for taking the time for the great feedback. I look forward to working on this more after Christmas and trying out a strategy for the MergedSection and CrossSection.

[connorferster]

New Gists:

https://gist.github.com/connorferster/43aa50e8c0c9eda971275e7adcf211fd

https://gist.github.com/connorferster/d393791e3af7ae885bc00c8553db73fc

Currently working on doc strings.

Are you open to switching your test suite from unittest to pytest? If so, I can take care of that.

[connorferster]

Also: happy new year!

[connorferster]

A couple of new gists:

https://gist.github.com/connorferster/9c822b4716f028a1c595dd8809a91f86

https://gist.github.com/connorferster/0f592c376d21112238a9d7b2f4f41895

I feel good about where everything has gotten to. I am going to focus on completing documentation and creating tests.

Since sectionproperties does not have many tests currently, are you open to switching your testing suite to pytest? I am going to be writing lots of new tests for the new Geometry class and it could be a good opportunity to switch to a modern testing framework. Essentially, your tests become easier to write (pure functions with assert statements) and you can get more information following each test run. I believe most modern python projects use pytest these days. I highly recommend it.

Let me know your thoughts on this

[Spectre5]

I haven't been involved at all with this PR, but I've been considering adding more tests for section-properties as well (it needs MANY more, I think) and I would also highly recommend switching to pytest, for what it's worth.

[robbievanleeuwen]

Happy for the testing switch! Great work 🎆 Tests have been on the to-do list for a while but never really got around to it!

[robbievanleeuwen]

Looking great @connorferster, loving the gists!

Re: reinforced concrete - I had thought about this a little bit in the past and wasn't sure how best to proceed. Problem - meshing circles (reinforcement) with straight lines underestimates the area of the circle. You can get close to the correct area by using a lot of subdivisions but this seems unecessarily expensive. One option would be to implement finite elements with curved edges but this is beyond my knowledge. Another option would be to use a fixed number of subdivisions e.g. 16 or 32 and set the diameter slightly bigger than the reinforcement bar but this would give you the exact area (I'm sure there's an analytical relationship between the 'increased' diameter required and the number of subdivisons). Not sure if you have any thoughts on this?

Another thought - one of my plans was to switch from meshpy to triangle (#24). Triangle appears to have the same functionality as meshpy but is well maintained and isn't such a pain to install (particularly for windows users). I don't think it would be much work to do the switch as the architecture is largely the same. Not sure if you have time or are interested in incorporating this into your massive overhaul 😄 ? If not, it's something I can tackle in the not too distant future.

Thanks again!

Robbie

[connorferster]

Sounds good, @robbievanleeuwen!

1. I started playing around with triangle but decided to revert my changes to "stay on track" and just focus on the preprocessing. I think moving to triangle would be good. This crashing situation with meshpy is not reliable and I do not know why it is happening because no exceptions or error codes are raised from the OS, just a silent exit. In Jupyter, I re-run once, twice, thrice and by then it usually works. However, for testing it is an untenable situation because it causes the test suite to exit early. I am currently unable to complete a test suite run (on either my local or on the Travis CI) because of this crashing issue. Have you experienced this crash (i.e. segmentation fault) issue on Travis CI before?

2a. Re: concrete - If that is the case, then it seems that the analysis will tend to be conservative, which is good. Given that FEA is almost always going to be an approximation of a physical system, I think that this discrepancy between meshed area vs actual area is within the realm of expectation of a "sophisticated" user (i.e. a user at least familiar with FE theory and its limitations).

2b. You mention on your main docs page that there are plans for a "concrete module". What did you have in mind for that? I am wondering if, after the preprocessing overhaul is complete, the next steps would be for the post-processing to consider both a compressive yield stress and separate tensile yield stress. I am new to the world of concrete analysis but I am under the impression that cracked section analysis requires an iterative procedure to account for the stress redistribution during the cracking process.

3. One thing that is tripping me up right now are the NASTRAN sections. Currently, I have left them as distinct classes instead of replacing them with functions that generate the required geometries. The reason for this is because of the method that is unique to the NASTRAN sections: .getStressPoints()

In my recent reading about NASTRAN, it seems that reporting stress values at these corner points on the section is characteristic of NASTRAN and should not be omitted. I had hoped to be able to dynamically calculate the locations of these stress points for an arbitrary section (allowing for stress point reporting on any shaped section) but it seems that, even for the NASTRAN sections, they are not consistently in the same place (e.g. in HAT section, two points lie in the mid-depth of the bottom flange instead of the bottom corners).

I like the cleanliness of having just two geometry classes that can be used to represent any geometry but I am unsure about the best way to handle these NASTRAN sections with their additional method. Currently, I am thinking of defining the stress points as an @property for all Geometry instances but the values would only be set by the .setter when the geometry was created with a create_nastran_section...() function. What are your thoughts on this? Once created, the recovery points could be dynamic with the section, e.g. if the shapely geometry is rotated, the recovery points would also be rotated and reported from that new rotated location.

[Czarified]

I just wanted to throw my two cents in here. I've been using your gists to learn shapely and follow along, and I'm really loving the re-write! Building up sections is pretty straight-forward, but I think it will deserve it's own section in the docs working through the operations with examples. For now, the gists are easy enough to follow. I'd be happy to help flesh out those docs from the examples you've already written up, once we're ready.

2. When I've worked with composite structures in the past, I'd have compressive and tensile properties, and change the one used based on the applied loads. And then the same for the allowables, using the correct one for the final local stress state. My experience, however, is with carbon-fiber composites, where the fiber-volume ratios are much higher than reinforced concrete. I agree with Connor on 2a, some slight conservatism is good.

3. I think a common stressPoints property is a good move. It may be beneficial to allow users to define their own stress recovery points, if they don't need the defaults or have a unique case. As an example, maybe I know a certain location is prone to damage, and I want to make sure I have higher margins at this spot, even if it isn't the most critical or at one of the default corners/intersections. I could run the analysis twice with a reduced area there and guess what my section could look like after damage, or I could run it once and pull results from this spot specifically to verify I have sufficient capability.

   • Either way, it would really benefit us Nastran users to be able to have loads on a bar element, plug them into a more refined cross-section here, and then directly compare stresses as a sanity check.

[Spectre5]

I think all cross sections should have definable stress recovery points. By default, for non-Nastran sections I'd recommend that we use the convex hull of the vertices. This may be more that 4 points. But the user should also be able to add extra ones and/or disable the default ones. Sometimes you want the stresses at the center of a flange, such as with shear dominated loading.

[aegis1980]

@connorferster

Re. your comment before...

1. I started playing around with triangle but decided to revert my changes to "stay on track" and just focus on the preprocessing. I think moving to triangle would be good. This crashing situation with meshpy is not reliable and I do not know why it is happening because no exceptions or error codes are raised from the OS, just a silent exit. In Jupyter, I re-run once, twice, thrice and by then it usually works. However, for testing it is an untenable situation because it causes the test suite to exit early. I am currently unable to complete a test suite run (on either my local or on the Travis CI) because of this crashing issue. Have you experienced this crash (i.e. segmentation fault) issue on Travis CI before?

This sounds like the problem I was having in my (attempt) [https://github.com//pull/40]
Never got to bottom of it. A couple of valid runs then exiting with no trace or anything. Tried everything to chase it down.

[robbievanleeuwen]

Hi all,

Keeping the number system we've got going:

1. I have run into the segmentation fault issue before and agree it is particularly frustrating. I seem to remember narrowing it down to, or at least attributing it to problems with the input geometry resulting in the mesh algorithm not being able to generate a mesh. From memory this was caused by duplicate or 'relatively' close nodes or facets that didn't make sense e.g. intersected other facets. Maybe have a careful look at the input geometry you are trying to mesh? Otherwise I can have a look at the suspect geometry and see if I can figure it out? The seeming randomness of the errors is puzzling though...

2. Agree with the approach to modelling circles and the conservatism with smaller areas. I guess the inner perfectionist in me wanted to get something more exact. Regarding a concrete module, was thinking of the following:

• Some pre-processor functions to easily define common reinforced cross-sections.
• Compute some code related properties, e.g. ultimate flexural strength/shear strength etc/stiffness etc.
• Some more advanced analysis methods e.g. M-N interaction diagram, parabolic stress profiles, moment-curvature stiffness diagrams etc.

3. Can't really comment a whole lot on the NASTRAN business as I have never used it and don't really understand the relevance of the stress points (this feature was added by another user). Happy to defer to others on this one, but your approach seems sound.

[connorferster]

@robbievanleeuwen

Thanks for that. Have taken a break from working on it this week. I feel like I am pretty darn close and want to get this complete so it can be reviewed and merged!

Road map for me:

1. Finish writing the nastran_...() functions for geometry generation.
2. Update the docstrings
3. Figure out what's happening with meshpy on your original tests and get them running reliably
4. Write some basic unit tests for Geometry and CompoundGeometry

Once this is complete, I am interested in writing a biaxial interaction diagram module. But that's for a little bit later. First things are first.

[connorferster]

I am going to cancel this pull request. Getting close to finishing. Will keep working on the branch and re-submit the pull request when it is ready.

[Spectre5]

For the rectangle, do you have a minimum working example for it that consistently fails? I'll try the same one to see if it fails for me or not too.

[Spectre5]

You can also set a PR to be "in work" in which case it can't be merged

[connorferster]

It's okay actually! I adjusted the way I generate Geometry.points and Geometry.facets so it does not include the first point again at the end (to close the loop). This seems to have fixed some of the problems I was having. I am able to get pytest to run consistently on test_rectangle.py but I am getting funny results on some properties, most notably qx and qy.

[connorferster]

And on my circular geometry, I am getting funny values for ixx_c and iyy_c, for example. Not sure why yet. Going to compare point orders but I don't know if that should create any difference since we have the facets lists to order the control the point connectivity.