Non-perpendicular hinge axis possible for multibody solver #299

kccwing · 2024-09-10T03:21:40Z

Changes to hinge_for and hinge_node_for, and unit tests for multibody beam formulation using the augmented lagrangian method - cases with dominating scaling vs penalty factor checked, as well as a global rotation of the double pendulum with non-perpendicular axis. 3 new tests take in total of 3 minutes to run (3*30 time steps)

Second derivatives of constraint equations evaluated manually - to be replaced with AD version in due course.

(+ temporarily commented code for dynamic trim)

fixes to both constraint types, for augmented lagrange multiplier method

1. der_Tan_by_xv: limit for small crv 2. der_skewp_v and der_skewpT_v

ben-l-p · 2024-09-10T09:04:25Z

sharpy/postproc/aerogridplot.py

+                        if struct_tstep.mb_dict is None:
+                            coords[counter, :] += struct_tstep.for_pos[0:3]
+                        else:
+                            #TODO: uncomment for dynamic trim


Does this need to be uncommented for dynamic trim, or can this be removed? (for this code and other similar comments in commit)

nope, these fix the issue that multibody simulations used to have outputs in the local A frame. The commented out portion will be useful for dynamic trim, but I haven't quite decided how best to proceed with it so I left the pointers there

sharpy/structure/utils/lagrangeconstraints.py

wong-hl · 2024-09-10T11:53:50Z

sharpy/structure/utils/lagrangeconstraints.py

+            mat = ag.multiply_matrices(-np.eye(3))
+
+            vec = ag.multiply_matrices(-node_cga, node_dot_Ra)                                    


Is it possible to replace the multiply_matrices() method with a numpy one or @ as Ben mentioned? If the numpy method is done on an instance of a numpy array it should be fairly easy to replace with the find and replace. Otherwise, it might be more involved using sed?? (I'm not 100% sure how feasible it is)

I am definitely a fan of the @ symbol for matrix multiplication - makes matrix algebra look much cleaner! Also, there's an identity matrix being multiplied by nothing here...

mat and vec structure is left intact for easier debugging and future work on these derivatives (until we replace the current approach) - open to replacing all with @ - will look into this and do it at least for the code involved in this PR.

wong-hl · 2024-09-10T11:54:48Z

sharpy/utils/generate_cases.py

@@ -1979,7 +1979,7 @@ def check(self):
                raise RuntimeError(("'%s' parameter required in '%s' lagrange constraint" % (param, self.behaviour)))
        has_behaviour = False
        for param, value in self.__dict__.items():
-            if not param in ['behaviour', 'scalingFactor', 'penaltyFactor']:
+            if not param in ['behaviour', 'scalingFactor', 'penaltyFactor', 'rot_axisA2']:


Rename rot_axisA2 to be consistent with the other params, e.g. rotAxisA2

rot_axisA2 for backwards compatibility for existing job scripts, as well as consistency with rot_axisB in lagrangeconstraints.py

Ah the joys of consistently inconsistent naming conventions 😂

wong-hl · 2024-09-10T11:58:34Z

tests/coupled/multibody/test_double_pendulum_slanted.py

+        LC2.scalingFactor = 1e6
+        LC2.penaltyFactor = 0.
+
+        LC = []
+        LC.append(LC1)
+        LC.append(LC2)
+
+        MB1 = gc.BodyInformation()
+        MB1.body_number = 0
+        MB1.FoR_position = np.zeros((6,),)
+        MB1.FoR_velocity = np.zeros((6,),)
+        MB1.FoR_acceleration = np.zeros((6,),)
+        MB1.FoR_movement = 'free'
+        MB1.quat = ag.rotation2quat(ag.multiply_matrices(ag.rotation3d_z(lateral_ini),ag.quat2rotation(np.array([1.0,0.0,0.0,0.0]))))
+
+        MB2 = gc.BodyInformation()
+        MB2.body_number = 1
+        MB2.FoR_position = np.array([node_pos2[0, 0], node_pos2[0, 1], node_pos2[0, 2], 0.0, 0.0, 0.0])
+        MB2.FoR_velocity = np.zeros((6,),)
+        MB2.FoR_acceleration = np.zeros((6,),)
+        MB2.FoR_movement = 'free'
+        MB2.quat = ag.rotation2quat(ag.multiply_matrices(ag.rotation3d_z(lateral_ini),ag.quat2rotation(np.array([1.0,0.0,0.0,0.0]))))
+
+        MB = []
+        MB.append(MB1)
+        MB.append(MB2)
+
+        # Write files
+        gc.clean_test_files(SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])
+        SimInfo.generate_solver_file()
+        SimInfo.generate_dyn_file(numtimesteps)
+        beam1.generate_h5_files(SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])
+        gc.generate_multibody_file(LC, MB,SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])
+
+        # sharpy_output = sharpy.sharpy_main.main(['',
+        #                                     SimInfo.solvers['SHARPy']['route'] +
+        #                                     SimInfo.solvers['SHARPy']['case'] +
+        #                                     '.sharpy'])
+
+        # Same case with penalty weights
+        global name_hinge_slanted_pen
+        name_hinge_slanted_pen = 'name_hinge_slanted_pen'
+        SimInfo.solvers['SHARPy']['case'] = name_hinge_slanted_pen
+
+        LC1.scalingFactor = 1e-24
+        LC1.penaltyFactor = 1e0
+        LC2.scalingFactor = 1e-24
+        LC2.penaltyFactor = 1e0
+
+        gc.clean_test_files(SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])
+        SimInfo.generate_solver_file()
+        SimInfo.generate_dyn_file(numtimesteps)
+        beam1.generate_h5_files(SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])
+        gc.generate_multibody_file(LC, MB,SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])
+
+        # Same case with rotated global reference
+        global name_hinge_slanted_lateralrot
+        name_hinge_slanted_lateralrot = 'name_hinge_slanted_lateralrot'
+        SimInfo.solvers['SHARPy']['case'] = name_hinge_slanted_lateralrot
+
+        LC2.rot_axisB = np.array([np.sin(45.0*deg2rad),np.cos(45.0*deg2rad),0.0])
+        LC2.rot_axisA2 = np.array([np.sin(45.0*deg2rad),np.cos(45.0*deg2rad),0.0])
+
+        gc.clean_test_files(SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])
+        SimInfo.generate_solver_file()
+        SimInfo.generate_dyn_file(numtimesteps)
+        beam1.generate_h5_files(SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])
+        gc.generate_multibody_file(LC, MB,SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])
+
+        #     # Same case with spherical joints
+        #     global name_spherical
+        #     name_spherical = 'dpg_spherical'
+        #     SimInfo.solvers['SHARPy']['case'] = name_spherical
+
+        #     SimInfo.solvers['NonLinearDynamicMultibody']['time_integrator'] = 'NewmarkBeta'
+        #     SimInfo.solvers['NonLinearDynamicMultibody']['time_integrator_settings'] = {'newmark_damp': 0.15,
+        #                                                                                     'dt': dt}
+
+        #     LC1 = gc.LagrangeConstraint()
+        #     LC1.behaviour = 'spherical_FoR'
+        #     LC1.body_FoR = 0
+        #     LC1.scalingFactor = 1e6
+
+        #     LC2 = gc.LagrangeConstraint()
+        #     LC2.behaviour = 'spherical_node_FoR'
+        #     LC2.node_in_body = nnodes1-1
+        #     LC2.body = 0
+        #     LC2.body_FoR = 1
+        #     LC2.scalingFactor = 1e6
+
+        gc.clean_test_files(SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])
+        SimInfo.generate_solver_file()
+        SimInfo.generate_dyn_file(numtimesteps)
+        beam1.generate_h5_files(SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])
+        gc.generate_multibody_file(LC, MB,SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])
+
+        # sharpy_output = sharpy.sharpy_main.main(['',
+        #                                      SimInfo.solvers['SHARPy']['route'] +
+        #                                      SimInfo.solvers['SHARPy']['case'] +
+        #                                      '.sharpy'])


Is it possible to break up the _set_up() method into smaller functions so that the method reads as a sequence of methods that explains the overall set-up is instead of needing to know all the individual steps?

Again, this is how the case generation class works as of current, which isn't part of this commit. It would maybe benefit from more documentation (I don't personally like using it as it is only convenient for simple cases), however as we already have a fair few cases defined using this, modifying it might not be a great idea.

If this is how they are all defined, is it worth having a base test class? Or are they sufficiently different?

There are two methods of creating a model. There is the approach of writing every parameter separately (chord, elastic axis for each element node, structure connectivity etc etc) which can be a bit tedious at times, but that's what I always use. This class creates all this for you for simple configurations, for example a single cantilever wing without taper, or in this case a pendulum. Both will create the same H5 files which are passed to the solver. I personally feel that this case is fine - it preserves the structure of a "normal" case so hypothetically if you wanted to extend this case for a full simulation, for example a triple pendulum, then you very easily could. I don't see much benefit here in a special definition just for tests.

ideally i would create a class for the pendulum setup just like the flying wings for the model build, then another higher level of abstraction for case setup like what @ACea15 had done, but yes I do see value in this, especially the latter - I should have a good version of it by the time I am back in February. Thanks!

wong-hl · 2024-09-10T12:03:20Z

tests/coupled/multibody/test_double_pendulum_slanted.py

+import shutil                 
+from sharpy.utils.constants import deg2rad
+
+class TestDoublePendulumSlanted(unittest.TestCase):


I know this is how we do unit tests in aerospace. But, this is not a unit test by the normal definition, as it should use mocking for the other components. This is technically an integration test. I'm ngl, I still haven't figured out what's an effective method for testing scientific code is 🤷‍♀️

This is how other tests are defined, so even if there is a better way to do this, I would rather they be consistent for now. We can look into swapping them all to a better definition in a future release.

I agree. I don't think we should change it, and I'm not 100% certain it would be worth it if we did????

Have you come across the hypothesis testing package? It's something I've briefly looked into and seems like it might be useful

indeed @wong-hl - I came across this when I was looking for a way to have a parametrised setup() within uniitest - haven't heard about mocking so I'll read up on it!

Mainly trimmed down lagrangeconstraints.py (numerics should be as before). Also cut some unneccesary comments I found along the way, and lastly cut the pesky multiply_matrices() function.

ben-l-p

I cleaned up some of the code last week to try and make it a little bit less intimidating. I have been using this solver for a while already and have found it seemes to work as intended - again, great contribution Kelvin!

kccwing · 2024-10-09T16:59:38Z

Appreciate the help with the PR @ben-l-p and @wong-hl! Thank you again :)

Shall I merge this one as well and then it goes in the latest (minor) version?

ben-l-p · 2024-10-09T17:04:01Z

Appreciate the help with the PR @ben-l-p and @wong-hl! Thank you again :)

Shall I merge this one as well and then it goes in the latest (minor) version?

I'm happy with it - if Hui Ling is too then we should be good to go. I would say we should leave this in develop now and merge to main in due course, as it shouldn't take me too long (don't quote me on this) to add my multibody edits too, so it can all be in one version.

kccwing · 2024-10-09T17:06:26Z

Brilliant! No worries, and thanks again for the help and verification 😊

wong-hl

LGTM

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

kccwing added 6 commits September 10, 2024 05:07

Test for non-perpen hinge axis added

ab832f8

Post processing fixes for for_pos

d6e471c

(+ temporarily commented code for dynamic trim)

Hinge_for and hinge_node_for fixes for derivatives

1cea00b

fixes to both constraint types, for augmented lagrange multiplier method

Algebra functions for hinge axis

648ebd3

1. der_Tan_by_xv: limit for small crv 2. der_skewp_v and der_skewpT_v

Maintaining backwards compatibility for rot_axisA2

6372f17

Maintaining backwards compatibility for rot_axisA2

78a283c

kccwing requested review from wong-hl and ben-l-p September 10, 2024 04:13

kccwing added the enhancement label Sep 10, 2024

ben-l-p reviewed Sep 10, 2024

View reviewed changes

sharpy/structure/utils/lagrangeconstraints.py Outdated Show resolved Hide resolved

wong-hl reviewed Sep 10, 2024

View reviewed changes

kccwing and others added 2 commits September 11, 2024 17:48

Removing commented code for LM_C - B^TB terms

3b5ecad

Cleaning up some code

369c2dd

Mainly trimmed down lagrangeconstraints.py (numerics should be as before). Also cut some unneccesary comments I found along the way, and lastly cut the pesky multiply_matrices() function.

ben-l-p self-requested a review October 9, 2024 15:43

ben-l-p approved these changes Oct 9, 2024

View reviewed changes

wong-hl approved these changes Oct 9, 2024

View reviewed changes

kccwing and others added 3 commits October 17, 2024 12:39

Big fix to read rot_axisA2 as optional argument indeed

b662a7f

avoid checking for equality of a list and string for custom yb_FoR input

c89cd2e

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Merge branch 'develop' into dev_hingeaxis

7f36d54

ben-l-p merged commit 69a4cdf into develop Oct 30, 2024
6 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Non-perpendicular hinge axis possible for multibody solver #299

Non-perpendicular hinge axis possible for multibody solver #299

kccwing commented Sep 10, 2024 •

edited

Loading

ben-l-p Sep 10, 2024 •

edited

Loading

kccwing Sep 10, 2024

wong-hl Sep 10, 2024

ben-l-p Sep 10, 2024

kccwing Sep 11, 2024

wong-hl Sep 10, 2024

kccwing Sep 11, 2024

wong-hl Sep 11, 2024

wong-hl Sep 10, 2024

ben-l-p Sep 10, 2024

wong-hl Sep 10, 2024

ben-l-p Sep 10, 2024

kccwing Sep 11, 2024 •

edited

Loading

wong-hl Sep 10, 2024

ben-l-p Sep 10, 2024

wong-hl Sep 10, 2024

kccwing Sep 11, 2024

ben-l-p left a comment

kccwing commented Oct 9, 2024

ben-l-p commented Oct 9, 2024

kccwing commented Oct 9, 2024

wong-hl left a comment

		mat = ag.multiply_matrices(-np.eye(3))

		vec = ag.multiply_matrices(-node_cga, node_dot_Ra)

Non-perpendicular hinge axis possible for multibody solver #299

Non-perpendicular hinge axis possible for multibody solver #299

Conversation

kccwing commented Sep 10, 2024 • edited Loading

ben-l-p Sep 10, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kccwing Sep 11, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ben-l-p left a comment

Choose a reason for hiding this comment

kccwing commented Oct 9, 2024

ben-l-p commented Oct 9, 2024

kccwing commented Oct 9, 2024

wong-hl left a comment

Choose a reason for hiding this comment

kccwing commented Sep 10, 2024 •

edited

Loading

ben-l-p Sep 10, 2024 •

edited

Loading

kccwing Sep 11, 2024 •

edited

Loading