Overhaul pre-processor to include shapely

[robbievanleeuwen]

Is your feature request related to a problem? Please describe. No specific problems, but would increase functionality and simplify a fair bit of the code.

Describe the solution you'd like Use shapely in all it's glory.

Describe alternatives you've considered Writing the code yourself - tried that, no fun.

Additional context Suggesting in issue #9 by Agent6-6-6.

[aegis1980]

Hi, Am sure I saw a 'DWG import' todo somewhere - this is more related to that, but I am using shapely. A WIP but (as part of another project) I made this: https://github.com/aegis1980/cad_to_shapely. I am having a few issues with making is work on PyPi (imports mostly: Python not my goto) but at some point it will be a pip-able thing. Currently only DXF (I ve written code for IGES and SVG - need to add it). Anyway, this gist imports a (Rhino3D generated DXF) aluminium extrusion (file attached) and spins it around for sectionproperties processing. Need to check I get the same answers as I get using Strand7 so unverified. @robbievanleeuwen Good work BTW ;) And a strong antipodean influence too. section_holes_complex.zip

[connorferster]

Hi @robbievanleeuwen,

I have some time before January and I would like to work on this. I am experienced with shapely and have previously started looking into how to integrate it into section-properties.

I will make a fork and start working on a branch so you can review the progress. As I go, let me know if you would like to advise me further or chat about the direction it's going.

Thanks for this excellent library! I think using shapely will be an excellent addition to it.

[robbievanleeuwen]

Hi @connorferster,

I initially had a few tries at implementing shapely but ran into a few hurdles and didn't delve into the docs deep enough to pull through. If you think you could make it work I think it would be a fantastic addition and potentially open up a few other possibilities with sectionproperties, so I'm more than happy for you to give it a try!

Give me a shout if you have any questions or would like to discuss anything 👍

Cheers, Robbie

[connorferster]

Now, I am able to get test_rectangle.py to consistently run in pytest and I no longer have my issues when creating the mesh in Jupyter when also using matplotlib to .plot_geometry(). Great news.

Additionally, my tests were crashing because of PlasticSection: I had not reviewed it or updated it yet to see if there were inconsistencies with the new pre-processor. However, my intention is to fit the new pre-processor into the existing CrossSection and analysis interface so it should not require updating.

I am however getting some funny values when calculating some geometric properties such as Ixx_c/Iyy_c and qx/qy. Looking into it now....

[robbievanleeuwen]

The plastic section calculation involves regenerating geometry and remeshing so I'm not surpirsed this isn't working yet! Give me a shout if you need some help understanding how it works 😄

[connorferster]

Thanks @robbievanleeuwen, will let you know as I get to it (next week?)

I found what was failing my tests with the qx, qy, Ixx_c, Iyy_c, etc. values. It seems that the reported values are actually Eqx, EIxx_c, etc. and not yet pure geometric value. It seems that these are values reported from the MeshPy elements and that sectionproperties would have to divide these values by E. I was running the tests with a "steel" material of E=200e3 and would not have expected the E value to affect the geometric properties. I realized this when I ran the display_values function and clued in that you do actually report them as E.qx, etc. Once I changed to E=1, I started getting exact match values.

Are you open to altering this behaviour so that the geometric properties are reported as pure geometric properties while the current reported properties being relabelled as something like Eqx?

[robbievanleeuwen]

The current behaviour switches to outputting 'material affected' properties as soon as a material property is specified. So if you wanted a purely geometric analysis with one material, there would be no need to specify the material.

I get that this is a little confusing and the reference to this behaviour is buried in the docs (probably not in the most obvious section). I think you've probably highlighted that this is not really a very good way of going about things and can cause a lot of confusion and wasted time (sorry!).

Definitely open to altering the behaviour, noting that as soon as you have more than one material a lot of the section propreties become undefined and therefore need to be modulus weighted.

[connorferster]

Aaah...I see it now. I take it back, I think your existing approach is best. I was thinking of the "default material" being a bit of an after-thought just because I am so used to having to specify a material to do any analysis in commercial software. However, it's a good solution: a gift to the user that I was foolishly rejecting!

[robbievanleeuwen]

I think it probably needs a bit more of a highlight in the docs - opened issue #53.

[connorferster]

Was modelling a funny-shaped hollow column section today in some commercial software. It was so painful. I wished I could have been using sectionproperties with the new pre-processor functions!!! Would have been so much easier! I really want to get this live...

[connorferster]

Hi @robbievanleeuwen,

I am noticing some funny behaviour with the geometric properties. They seem to be dependent upon the section's absolute location within the cartesian plane. For example, using the rectangular_section() function to create a box centered on the origin, I generate the following results: Note: qx, qy, specifically.

If I shift the section so that the origin is the bottom left corner, I generate the following results (which would be considered the "correct" results):

If I shift the section further still:

Note: Ixx_g, Iyy_g, and Ixy_g also get out of whack.

Have you encountered this behaviour before? Is this expected and I just missed it in the docs?

Generating code:

[robbievanleeuwen]

Hi @connorferster this is expected behaviour. Ixx_g and qx are taken about the global axis, which is (0,0). So these values should change as you move the location of the section around the cartesian plane. I would expect qx=qy=0 when the global axis coincides with the centroid in example 1. Ixx_c is taken about the centroidal axis so remains unchanged in the various examples. I think the various axes are described in the docs, e.g. here.

[connorferster]

It seems I changed the implementation of rectangular_section() (at some point...for some reason...?) to return a rectangle centered at the origin, instead of your original implementation centered on b/2, d/2. No idea why I did that, especially since I want to keep all of your shape generating procedures the same as you made them. This explains why the rectangular_section() results were inconsistent...solved!

[connorferster]

Hi @robbievanleeuwen,

I am having trouble with the plastic properties: some of the values with the shapely pre-processor ("shapely") are incorrect when compared with the main branch ("main"). I am investigating it by using the same CircularSection used in the examples. Here are some of the things I have found so far:

1. Point coordinates and facets are identical between "shapely" and "main"
2. Number of mesh elements are the same between "shapely" and "main"
3. Plastic properties for each element are identical between "shapely" and "main" for some elements but different for other elements
4. Number of iterations to converge with the brent root finding algorithm (via CrossSection.pc_algorithm) is just over double (four iterations in "main", 8-10 iterations in "shapely")

My plan is to next to check is to see if every Tri6 element is identical between "shapely" and "main". However, I am wondering if these symptoms are familiar to you during your development of the library.

Edit:

It seems that just under 50% of the plastic mesh elements are different between "shapely" and "main". I am confused by this given that the points and facets used as the input are identical and I am using the same mesh size for each.

Edit 2:

The original CrossSection mesh nodes and mesh elements are identical between "shapely" and "main". Only the plastic mesh are different (and only ~50% of them).

Edit 3:

Aha! I have narrowed the problem down to the part where the second mesh is generated based off of the shifted section (the original mesh shifted so the centroid is at the origin). The elements are identical in the first mesh but then get messed up in the second, shifted mesh....

Edit 4: So it seems that the values for the centroid are slightly different between "shapely" and "main" but that, to me, the difference should be immaterial to the meshing especially since the generating code for this geometry is creating a geometry centered on the origin. Also, given that the original CrossSection points, facets, meshes, and Tri6 elements are identical in both "shapely" and "main", why would there be this tiny difference in the calculation of the centroid?

• cx, cy ("main"): -1.6133929953056836e-15 -3.044650639835494e-16
• cx, cy ("shapely"): 5.436876142563384e-16, -1.388215708401184e-15

To confirm that the Tri6 elements in both "shapely" and "main" are identical, I set up logging in PlasticSection.calculate_centroid() to log the contents of the incoming elements argument. Using Notepad++ to run the comparison, they are identical and the code in calculate_centroid() is untouched (by me, I mean).

Edit 5: Image showing comparisons of Tri6 elements (showing apparently identical elements) and individual centroid calculation results for those elements (which shows discrepancies in the 12th to 14th decimal places:

Having stepped away from the problem this afternoon, obviously the input elements must be different in order to get different calculation results. The __repr__() method auto-created by dataclasses.dataclass only shows me eight decimal digits but the discrepancies are occurring in the 12th-14th decimal places. To accurately compare in the logs, I will need to write a new __repr__, etc.

I think I will start at the top of the chain and work down until I start to find discrepancies.

[connorferster]

Hi @robbievanleeuwen, @Agent6-6-6, @Czarified, @BenjaminFraser, @Spectre5

Further to the chain of information above:

Edit 6: The Section/CrossSection mesh nodes and elements and attributes are all identical to the 16th digit at time of instantiation. Additionally, the Tri6 elements between "main" and "shapely" are also identical to the 16th digit.

Edit 7: Once the Section/CrossSection is instantiated, the next step I take is to do .calculate_geometric_properties. Almost all of the individual calculated values in each element (e.g. area, qx, qy, etc) are divergent from "main" at either the 12th or 13th decimal place. The occasional element has a matching value. I thought that the problem may have been the incremental sum but the problem is that each individual element's properties are slightly different.

I have found the very first instance of divergence in the calculation: it is the calculation for the Jacobian. When it is calculated in "shapely", using the exact same inputs (exact match to 16th decimal place, anyway), the result is divergent in the 15th and 16th decimal place (sometimes the 14th place).

However, I am stumped on what to do about this and am in need of some new ideas from you or other in this small community before pursuing testing further.

[connorferster]

Shapely log file: https://drive.google.com/file/d/1c2tlHwUAV7uZeiXkZkRFBS9Gp3BhTRjo/view?usp=sharing

"Main" log file ("legacy"): https://drive.google.com/file/d/16MiJFBQHEJwQc_FboirErJYlkFGZrqGw/view?usp=sharing

[Agent6-6-6]

@connorferster, are you possibly just running into numerical precision issues here?

While the numbers differ, they are still pretty darn close to zero!

I wonder if you're taking more iterations (presumably using Brent's solver in scipy.optimize.brent), then there's something to look at there. Because it seems odd that this would be the case for what is the exact same problem as you describe it.

Does the same occur with other section types?

I must admit I've only partly kept up with this thread and not taken the shapely version for a spin yet, nor looked at the code specifically.

The only time I've seen Brent algorithm spit out slightly different answers (and also generally take more iterations) is if you set the tolerance to large, for example having for example 5 decimal places of precision might take more iterations than 8 decimal places. It sounds counter intuitive but going smaller sometimes results in less iterations, it's probably something to do with how many times you're evoking each type of interpolation in the algorithm or some other weirdness.

I'm sure you've checked you're using exactly the same input parameters for scipy.optimize.brent.

[connorferster]

Definitely running into floating point precision errors here :)

What is really stumping me is the question of where the error originates? I am using Robbie’s exact same code to generate point coordinates with the only real difference being they first go into a shapely object before going into Geometry.points. So is something happening within the shapely object somewhere beyond the 16 decimal that is preventing me from getting matching results?

That is my thought for my next test. But, gawrsh...

Sent from my iPhone

On Feb 8, 2021, at 23:30, Agent6-6-6 notifications@github.com wrote:

 @connorferster, are you possibly just running into numerical precision issues here?

While the numbers differ, they are still pretty darn close to zero!

I wonder if you're taking more iterations (presumably using Brent's solver in scipy.optimize.brent), then there's something to look at there. Because it seems odd that this would be the case for what is the exact same problem as you describe it.

Does the same occur with other section types?

I must admit I've only partly kept up with this thread and not taken the shapely version for a spin yet, nor looked at the code specifically.

The only time I've seen Brent algorithm spit out slightly different answers (and also generally take more iterations) is if you set the tolerance to large, for example having for example 5 decimal places of precision might take more iterations than 8 decimal places. It sounds counter intuitive but going smaller sometimes results in less iterations, it's probably something to do with how many times you're evoking each type of interpolation in the algorithm or some other weirdness.

I'm sure you've checked you're using exactly the same input parameters for scipy.optimize.brent.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[Spectre5]

I am having trouble with the plastic properties: some of the values with the shapely pre-processor ("shapely") are incorrect when compared with the main branch ("main").

How different, in the end, are your results compared to current master ("main")?

[Czarified]

This may be a silly question, but does the exact same geometry input (in master) always generate the exact same output to those 14th-16th decimal places? It may be some chaos introduced with the optimization algorithm. Since the differences are so small, I would tend to note the difference, and ignore them for now; relying more on an overall test for similarity.

Yes, it's very curious that they're different, since the shapely object should be exactly the same input. Is the modified code with your logging on your branch already?

The fea.shape_function function doesn't seem to be doing anything that would cause J to differ. It's all based on the gauss points, which seem to be hard-coded anyway. So I'm not sure where the deviation on the Jacobian is even coming from. The only difference is the coords np.array passed to shape_function. But you confirmed already that the points are exact out to 16 digits, right? So coords should match exactly between master and shapely-branch. 🤔 😕

EDIT: Another difference, I don't completely understand: Tri6.geometric_properties uses 6-point Gaussian integration, while Tri6.plastic_properties uses 3-point. Is that intentional? I don't know the numerical methods as well as Robbie does.

[connorferster]

@Spectre5

How different, in the end, are your results compared to current master ("main")?

For the geometric properties: a very small fraction of a percent with discrepancies for final results occurring in the 6th to 8th significant figure. For the plastic properties: the results are just wrong.

Observing the incorrect results (and incorrect intermediary values) started me on this chase. I noticed the slight discrepancy in the geometric properties earlier on but, because they were so small, they did not particularly bother me. However, once I noticed that the plastic properties were completely incorrect, I started wondering if there was a relationship between these two things.

@Czarified

This may be a silly question, but does the exact same geometry input (in master) always generate the exact same output to those 14th-16th decimal places?

I don't think that's silly at all! That was something I checked early on, too. All of my runs on the "main" branch with the circular section geometry exactly matches the results published in the first example in the Docs. So, I scratched that possibility out.

Is the modified code with your logging on your branch already?

Just pushed it. You can find it on my repo under the "shapely..." branch

[robbievanleeuwen]

Hi all,

Unpacking the above, there seems to be two issues:

1. Precision errors - I don't feel like I have a deep enough understanding of computer science to be able to help on this one. And perhaps this is my naivety/lack of compsci knowledge, but are we particularly concerned about variations in the 12th decimal place? My guess would be that there are some extra operations within shapely that mean precision comes into the equation but that's just an uneducated hunch.

2. Plastic properties - I wouldn't expect the algorithm to take twice as long to converge. There are some useful debugging tools built into the plastic section analysis which might help you get to the bottom of this. They visualise each step of the algorithm and show the trial plastic neutral axis geometry and the generated mesh as well as print some useful info to the terminal. I would try setting verbose=True and debug=True - docs here. This may help you compare what the plastic algorithm is doing in each step for the main branch and the shapely branch.

@Czarified to answer your question: plastic analysis only computes areas - when we compute areas 3 point integration gives exact results (in fact one point integration would be enough as the Jacobian is constant in the Tri6 element!) The geometric analysis requires higher order integration as we are integrating higher order functions which vary over the Tri6 element, e.g. the computation of second moments of area has quadratic and mixed terms - this needs six point integration to get an exact result. I talk about this briefly in section 4 of this.

[Agent6-6-6]

I'm not concerned about the 12th decimal place being inaccurate. If there are minor differences and it's decided that 12 significant figures is where it's at in terms of accuracy, then just state that in the docs and you could even give an option to hard round to that many significant figures if desired so answers are always the same. The fact that you have to define curved regions as a discrete number of straight facets kind of throws even achieving a few decimal places of accuracy out the window anyway for practical uses.

Perhaps some of the plastic properties can be worked out using shapely, taking the section perimeter and splitting it using a horizontal line. Work out plastic modulus from centroid of each divided half once you have the dividing line in the right location to make area each side equal (or if you have different yield stresses use the values weighted by yield stress). I say this as an alternative method to just to give you another result to see if it is the same or not as one of the results to understand if there is an issue. You don't need to mesh the full section necessarily to determine the plastic axis and plastic section modulus unless I'm mistaken.

[connorferster]

Hi @robbievanleeuwen, @Agent6-6-6, @Spectre5, @Czarified,

First, thank you for your support and suggestions. I really appreciate it!

Update: I believe I found the problem with the plastic section properties. I had had a previous exchange with Robbie where I noticed that there were methods in Geometry for adding individual lines and points to the section and I had taken them out after we both thought that would no longer be necessary with a shapely interface. As such, I began commenting out lines that referenced these methods (or lines that I thought were referencing those methods in CrossSection). As such, I broke a critical piece of the plastic section calculation where the section is cut along the principal axes.

However, I think that shapely can have an easier time of this with shapely.ops.split(geometry, splitter) and I can shorten the code for PlasticSection.add_line.

Going through this preprocessor update, I am really appreciating the amount of work Robbie did to build a geometry engine from scratch. It's excellent.

[robbievanleeuwen]

Hi @connorferster,

Sorry I remember this exchange and it seems like I gave you some dodgy advice - looks like I did use these methods after all!

I think that both your suggestion of using the split method in shapely and @Agent6-6-6's suggestion of not requiring meshing are great! The only reason I need to remesh in the plastic algorithm is to recompute the area above and below the trial plastic neutral axis - however if there is a way to get this info from shapely we can remove the constant remeshing from this algorithm completely and will simplify the whole thing!

Thanks again :)

[Spectre5]

Great investigative work! I like the idea of potentially looking into other solutions for getting area "above" and "below". It would probably be faster too instead of meshing and re-meshing many times.

[connorferster]

@robbievanleeuwen

I was scratching my head for the last couple days wondering why my warping properties tests were not passing the test suite where yours did! Then I noticed that you did not actually test the warping properties ;) When I actually connected the test for warping properties in test_rectangular.py then the results from the main branch failed the tests also.

When I checked my results against the main branch results, they were the same but neither (completely) match the analytical results from the test suite (within the prescribed tolerance, anyway).

I will change the test suite expected results to match the results calculated from the main branch, taking the main branch results as correct (for now, anyway).

[connorferster]

Hi @robbievanleeuwen,

The original test suite now runs and passes, which is great. However, I wanted to a create a section that is a more elaborate test but the darned thing keeps crashing in the last iterations of the plastic section calculation. It is crashing Triangle. Here is the error message I can generate when I run a my test script from the console:

Internal error in segmentintersection():
  Topological inconsistency after splitting a segment.
  Please report this bug to jrs@cs.berkeley.edu
  Include the message above, your input data set, and the exact
    command line you used to run Triangle.

Here is a gist of the same script to demonstrate what I am doing: https://gist.github.com/connorferster/46d46881b10e844a8bd90ed926ec7d58

It seems that some funny geometry is being created which is crashing Triangle. However, when I use the debug=True flag in .calculate_plastic_properties() I cannot discern where the broken geometry is.

I have gone line-by-line through the small amount of modifications I have made to cross_section.py, fea.py, and pre.py to try and find something that is breaking it. Of course, when I output the points, facets, holes, and ctrl points from my "compiled" geometry and feed it into the main branch as a CustomSection, it works. I have been doing line-by-line comparisons of the main branch to try and find what the problem is and I cannot find it.

Would you be willing to do a web conference chat sometime in the next week to review together? I am hoping that you may be able to see something that I cannot see.

I am so close to getting this closed out and a pull request submitted. Hoping that we can get this fixed together.

Many thanks,

Connor

[Spectre5]

I can understand wanting to get to the bottom of this issue, even just to be sure there isn't some deeper or bigger bug. I know I'd want to. But that said, once figured out, are you still thinking about looking into switching to shapely to do this plastic computation part instead of using triangle to continually remesh it? Or do you feel that is better saved for a separate PR?

[connorferster]

I think it is a separate PR. I started looking into it over the last week and I realized that there would be significant code changes to cross_section.py to make it work. It is well beyond the scope of what I want to contribute on this PR (which itself is a pretty big scope!). While I have previously been playing around with the mechanics in shapely that can perform this calculation, I have not gotten the process 100% clear in my mind yet.

The shapely preprocessor is connected so that the data coming in to the analysis modules should appear no different than if it were coming from the Geometry objects on the main branch. I fixed the previous plastic sections problem by removing a “modification” (aka naively breaking the code) I made to cross_section.py. So I am surprised now that this new issue has occurred.

On Feb 15, 2021, at 22:53, Spectre5 notifications@github.com wrote:

 I can understand wanting to get to the bottom of this issue, even just to be sure there isn't some deeper or bigger bug. I know I'd want to. But that said, once figured out, are you still thinking about looking into switching to shapely to do this plastic computation part instead of using triangle to continually remesh it? Or do you feel that is better saved for a separate PR?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[Spectre5]

Makes sense. Have you tried rounding the points location that Triangle gets? It might at least tell you if the issue is another numerical precision issue

[robbievanleeuwen]

Hi @connorferster, sorry I've been MIA the last 2 weeks. I'm back at uni (again!) and it's just kicked off so I've been a bit preoccupied!

Good find with the testing, it seems like I never finished that file - very strange!!

Very happy to have a chat via a video call, would be nice to meet as well! Maybe the easiest way to organise would be to flick me an email at robbie.vanleeuwen@gmail.com and we can find a time that suits both our time zones ;)

[connorferster]

@robbievanleeuwen @aegis1980 @Czarified @Spectre5 @Agent6-6-6

The new shapely pre-processor is finished and 100% operational. I have just made a push to my fork and have posted this (lengthy) gist to showcase most of the new features:

https://gist.github.com/connorferster/29a93ab1d07741d9ab4306a8e458a083

As I understand, Robbie is quite busy on a new venture right now. A community has been gathering around sectionproperties which I think is really special. I would really appreciate for all interested parties to review the above gist and to review the code in the new sectionproperties.pre.sections file, where the lion's share of the pre-processor resides. How do you feel about the code? Is it reasonably easy to understand? Do you think it looks maintainable?

Additionally, if people have time and interest, please try cloning my fork onto your machine and install it with setup.py and try it out for yourself!

At the end of the gist, are three action items I would love, love, love to have some help with. Particularly updating the documentation. Please drop me a line if you are interested in assisting with updating examples in the documentation with me. If there are more of us pitching in, we could probably get it done in a short weekend session!

Thanks for everyone's support and interest in this new pre-processor! I hope that we can continue building sectionproperties for years to come so that there is a free software option for this fundamental piece of mechanical and structural engineering.

Edit: After a conversation with Robbie, I did end up implementing the plastic section calculations with shapely. So much easier (and faster to calculate)! Works beautifully.

[Agent6-6-6]

@connorferster, finally taking this for a spin to test....

I thought it would be useful to start a bit of a list on things that are different to enable people to convert existing scripts once things get merged, working through converting some older scripts to the shapely version I've noted the following so far:-

• CustomSection no longer exists, similar functionality is provided now by .Geometry.from_points
• All the previous built-in geometry definitions have been renamed slightly (typically the 'shift' parameter has also been removed). Can still shift using shift_section after defining geometry.
• It might be useful to have a before/after table correlating the above name changes in the documentation to help people mix and match the changes
• mirror_section, shift_section, etc works differently than previous, you can no longer apply directly to a created geometry, before you could go geo1.mirror_section(….), now you have to go geo1=geo1.mirror_section(….) (as far as I can tell if I'm interpreting things correctly?)
• Mirror_section works differently (see point below as have not quite figured it out....)
• shift_section works differently, now accepting two separate variables for x and y offset instead of a list. Can no longer be applied like noted above for mirror_section example directly to a geometry instance
• MergedGeometry is now renamed CompoundGeometry

Couple of other points as I progress through converting a single script that I've noted:-

• I had to manually change the shapely.log path in sections.py as it was hardcoded to your desktop path (probably not intended for final code obviously)
• I required to install icecream package, this isn't a requirement in the setup? Include in requirements if intended for final code?
• I required to update my matplotlib (had very old version (like 3.1.2)) that apparently didn't have matplotlib.colors which is required somewhere, might need to enforce some minimum version here?

This is probably me being dumb, but the example given in the docs for mirror_section doesn't provide an example of mirroring about an arbitrary point, and since things have changed I'm not entirely sure how to get it to work.

Previously this following worked geometry5.mirror_section(axis='x', mirror_point=[0, 0])

Looking at the docs it has Union[list[float, float], str] as the type, what is the str if I specify a specific point? I'm missing something but not sure what?

I can do geometry5 = geometry5.mirror_section(axis='x', mirror_point=[0, 0, 0]) and it gets past the line but simply trying anything to try get things to work ... :), not sure what I'm supposed to put in the position of the third zero? Throws an error if you only have [x, y] coordinates as implied in the docs? TypeError: can only concatenate list (not "tuple") to list

Can you provide a small example of how mirroring might work about an arbitrary point as I'm clearly missing something (probably obvious)! Thanks

[Agent6-6-6]

@connorferster One other thing I just noted, in plot_mesh, materials=True is the default. This puts the label on the plot, but if you have no materials defined I'd probably prefer this to be set to false as the default option. Presumably you can check if there are materials defined, and if none are defined set this argument to False? This was the previous behaviour, a plot of the mesh produced with no legend.

Theres also a debugging print(legend_list) that is probably left over in plot_mesh? https://github.com/connorferster/section-properties/blob/77cd4a8b36a4d335b41a74f70a3a894d8771f68b/sectionproperties/analysis/cross_section.py#L1072

One other possible improvement:- Is it possible when you have multiple geometry elements and if one mesh_size is specified with create_mesh to just to apply to all elements irrespective of number as a default behaviour? Currently you need a list with the same number of mesh_size entries as geometry elements. But I find I'm usually just applying the same mesh size throughout.

[connorferster]

@Agent6-6-6

Agreed! A comparative list would be a good plan.

You have certainly caught some of old instrumentation there (logs, icecream). Thanks for catching those. Will remove.

More tomorrow...

[connorferster]

Hi @Agent6-6-6

Thanks for pointing out those gaffs. I have fixed .mirror_section() and removed the instrumentation and print statements.

Regarding a list of old vs new, do you think it would suffice to redo all of the current examples in the documentation with the new preprocessor, showing them side-by-side? Or with old code commented out with new code below each line?

What do you think the best format may be for presenting the comparison?

Thanks!

[Agent6-6-6]

@connorferster, I was imagining just a table or list describing the differences as a separate page in the docs would suffice. It's really just to give users a bit of a headstart in making the required adjustments to their previous scripts. But an example that uses/demonstrates the majority of the changes to old methods and demonstrates the different approaches/changes compared side by side has merit with some appropriate code comments explaining the differences. I don't think you need to explain any new features here, just what might need to change in older code to update for shapely pre-processor.

I'm not sure how @robbievanleeuwen intended to treat the new shapely version, as next version 1.09 or whatever its up to, or a completely new package? That might have some bearing on how the docs are approached. I'd imagine it's just a natural progression to next version within the same package. So all the old methods become redundant and if someone wants to find out about them they can install an older version or review the older versions of the docs.

I don't really see any benefit in providing both versions in the docs (not withstanding it would be a significant amount of work to do both I'd imagine). If someone wants to learn about the previous code nomenclature they can drill back to an earlier version of the docs anyway.

I'll let you know if I note anything else if you've fixed mirror_section 👍

[connorferster]

Thank you @Agent6-6-6 for reviewing the code. I really appreciate it!

That sounds good. I will play around with a couple of ideas to see what is perhaps more communicative and comprehensible. The documentation will need to be updated for the new version anyway. I can add a page of "API changes from v1.x.x" or something to that affect.

Because this will be introducing several breaking changes to the code, I think the v2.0.0 will be most appropriate (if following semver, which I think is a good idea).

[Agent6-6-6]

Hi @connorferster, I've only just got back onto testing your latest changes and I note in fea.py there is the same hardcoded log path that was in sections.py at the start. Probably needs some tidying up as well as currently errors on the following code:-

https://github.com/connorferster/section-properties/blob/3266c21313d63673b6b3c9c8048bb5fec7fdfb2c/sectionproperties/analysis/fea.py#L11

[Agent6-6-6]

@connorferster, the following lines put into create_mesh routine after mesh.set_holes(holes) # set holes line will achieve this request:-

One other possible improvement:- Is it possible when you have multiple geometry elements and if one mesh_size is specified with create_mesh to just to apply to all elements irrespective of number as a default behaviour? Currently you need a list with the same number of mesh_size entries as geometry elements. But I find I'm usually just applying the same mesh size throughout.

    if len(mesh_sizes) == 1: # if single mesh size value provided, create meshes for all regions at this size
        mesh_sizes = mesh_sizes * len(control_points)

[connorferster]

Agreed! Easy fix.

[Czarified]

@connorferster , at this point, do you want to submit it as a PR, so we can give a more "formal" review? GitHub has built-in functionality for this, and it will be documented in the same PR that implements the changes, for record-keeping. I'm going to take an initial look today, but I don't want to build/install the branch on my work computer right now. I'll dig in deeper at home.

[Agent6-6-6]

Yeah will do.

[connorferster]

Confirmed! Will do this eve.

[connorferster]

PR Submitted

[Agent6-6-6]

Hi @connorferster This just tripped me up as another subtle change, argument mesh_sizes has changed to mesh_size in the Geometry Class create_mesh function.

This is an issue because if you have code that could result in either a single or multiple geometry elements, then the expected/required kwarg seems to change... for example adding or not adding a strengthening plate to a I-section based on some code logic requires different arguments when you create the mesh. mesh = geometry.create_mesh(mesh_sizes=[mesh_area]) vs mesh = geometry.create_mesh(mesh_size=[mesh_area])

Then if you pass a single geometry element it fails using mesh_sizes as an argument, but if you pass a compound geometry it works using mesh_sizes as an argument. Conversely it will work with a single geometry element using mesh_size as an argument, but fail if there is compound geometry using mesh_size

If that made no sense I can write an example to demonstrate. But I believe it should universally be one or the other here if there is no reason to change to mesh_size, then keep as it was?

Either have it one way or the other

(appreciate the documents not fully updated, but in most locations it would need updating to reflect the final decision)

[Agent6-6-6]

@connorferster , also the external perimeter output seems to be incorrect for compound geometry. I think it's returning the bounding box perimeter rather than the true external perimeter of a section (see example below). Still seeing if I can figure out why in the code, but you might know why if more familiar with shapely and how it's implemented.

For example:- perimeter = 1788mm (section tables put 310UC97 as 1790mm)

compared with adding 2 330x16 side plates, external perimeter should be 1378mm but the code is returning 1334mm which is the perimeter of a bounding box

what it seems to be reporting back is this, vs the true perimeter:-

Compared to v1.08 I don't think the perimeter output worked at all for merged geometry, just reporting zero.

Seeing almost double the speed executing code which is really nice 👍

[Agent6-6-6]

Regarding the perimeter, finally found it in the code, it's because section-properties is using the convex_hull function from shapely. What we require in this situation is the concave hull of the compound geometry.

Unfortunately, shapely does not include a concave hull formulation. The problem being there is no single answer to the problem, it's really dependant on the concave hull algorithm you're using.

I found this which may be of some use. https://gist.github.com/dwyerk/10561690

Will do some more reading....

[Czarified]

What results require the perimeter? This is good to know, but not really necessary for any stress calculations, right? Maybe we can just note it in the docs. Since you've outlined it so well here, we can copy most of it over.