updated nonlinearstatic and dependencies to accept optional load ramp…

…ing (with test files)

sharpy/routines/static.py

@@ -21,6 +21,7 @@ def sol_101(self,
                 s_delta_curved=1e-3,
                 s_load_steps=1,
                 s_relaxation=0.01,
+                l_ramp=False,
                 modify_settings=None,              
                 **kwargs):
 
@@ -43,6 +44,7 @@ def sol_101(self,
                                              min_delta=s_tolerance,
                                              delta_curved=s_delta_curved,
                                              num_load_steps=s_load_steps,
+                                             load_ramping=l_ramp,
                                              relaxation_factor=s_relaxation)
         if primary:
             if modify_settings is not None:

sharpy/solvers/nonlinearstatic.py

@@ -28,6 +28,18 @@ class NonLinearStatic(_BaseStructural):
     settings_default = _BaseStructural.settings_default.copy()
     settings_description = _BaseStructural.settings_description.copy()
 
+    # additional setting for load ramping - disabled by default, only enable with structural only simulations
+    settings_types['load_ramping'] = 'bool'
+    settings_default['load_ramping'] = False
+    settings_description['load_ramping'] = 'Flag to enable load ramping'
+
+    # additional parameter for load ramping - only used if load ramping enabled to decide ramp steps
+    # initial value does not matter -> value will be set and modified at runtime
+    settings_types['load_ramping_conv'] = 'bool'
+    settings_default['load_ramping_conv'] = False
+    settings_description['load_ramping_conv'] = 'Flag for convergence and subsequent ramping'
+    # end of load ramping edit
+
     settings_types['initial_position'] = 'list(float)'
     settings_default['initial_position'] = np.array([0.0, 0.0, 0.0])
 
@@ -47,10 +59,100 @@ def initialise(self, data, custom_settings=None):
         settings.to_custom_types(self.settings, self.settings_types, self.settings_default)
 
     def run(self):
-        self.data.structure.timestep_info[self.data.ts].for_pos[0:3] = self.settings['initial_position']
-        xbeamlib.cbeam3_solv_nlnstatic(self.data.structure, self.settings, self.data.ts)
-        self.extract_resultants()
-        return self.data
+        # print(self.settings)
+        # load ramping edits - unchanged if flag = False for compatibility with FSI simulations
+        if not self.settings['load_ramping']:
+            self.data.structure.timestep_info[self.data.ts].for_pos[0:3] = self.settings['initial_position']
+            xbeamlib.cbeam3_solv_nlnstatic(self.data.structure, self.settings, self.data.ts)
+            self.extract_resultants()
+            return self.data
+        else:
+            self.data.ts = 0
+            self.settings['load_ramping_conv'] = 0
+            load_ramping_marker = 0
+            while not self.settings['load_ramping_conv']:
+
+                if load_ramping_marker >= 5:
+                    print('Maximum load step doubling reached (5 times). Check your case to determine cause. Exiting...')
+                    break
+
+                load_ramping_marker += 1
+
+                print('Solution running with %d load step(s). ' % self.settings['num_load_steps'])
+
+                for i_step in range(self.settings['num_load_steps'] + 1):
+#                   print(i_step)
+                    if (i_step == self.settings['num_load_steps'] and
+                    self.settings['num_load_steps'] > 0):
+                        break
+
+
+                    # load step coefficient
+                    if not self.settings['num_load_steps'] == 0:
+                        load_step_multiplier = (i_step + 1.0)/self.settings['num_load_steps']
+                    else:
+                        load_step_multiplier = 1.0
+
+
+
+
+                    # new storage every load step
+                    if i_step > 0:
+                        self.data.ts += 1
+                        self.next_step()
+
+
+
+                    n_node, _ = self.data.structure.timestep_info[self.data.ts].pos.shape
+                    struct_forces = np.full([n_node,6], 0)
+                    # print(struct_forces)
+
+                    # copy force in beam
+                    old_g = self.settings['gravity']
+                    self.settings['gravity'] = old_g*load_step_multiplier
+
+#                   print(self.settings['gravity'])
+#                   print(load_step_multiplier)
+
+                    temp1 = load_step_multiplier*(struct_forces + self.data.structure.ini_info.steady_applied_forces)
+                    self.data.structure.timestep_info[self.data.ts].steady_applied_forces[:] = temp1
+                    # run beam
+                    self.data.structure.timestep_info[self.data.ts].for_pos[0:3] = self.settings['initial_position']
+
+#                   print('Before solver: '+str(self.settings['load_ramping_conv']))
+#                   print('Before solver: '+str(self.data.structure.timestep_info[self.data.ts].pos))
+                    xbeamlib.cbeam3_solv_nlnstatic(self.data.structure, self.settings, self.data.ts)
+#                   print('After solver: '+str(self.settings['load_ramping_conv']))
+
+#                   print('After solver: '+str(self.data.structure.timestep_info[self.data.ts].pos))
+                    self.extract_resultants()
+
+
+                    self.settings['gravity'] = old_g
+#                     (self.data.structure.timestep_info[self.data.ts].total_forces[0:3],
+#                         self.data.structure.timestep_info[self.data.ts].total_forces[3:6]) = (
+#                             self.extract_resultants(self.data.structure.timestep_info[self.data.ts]))
+
+
+                    # convergence
+
+                if self.settings['load_ramping_conv']:
+                    self.cleanup_timestep_info()
+                    print('Converged. Exiting...')
+                elif not self.settings['load_ramping_conv']:        
+                    print('Solution did not converge with %d load step(s). Doubling number of load steps...' % self.settings['num_load_steps'])
+                    self.settings['num_load_steps'] *= 2  
+                    self.del_timestep_info()
+                    self.run()
+                else:
+                    self.settings['load_ramping_conv'] = True
+                    self.cleanup_timestep_info()
+                    print('Unknown error 14239, exiting...')
+            return self.data
+
+
+
+
 
     def next_step(self):
         self.data.structure.next_step()
@@ -81,5 +183,23 @@ def create_q_vector(self, tstep=None):
             tstep = self.data.structure.timestep_info[-1]
 
         xb.xbeam_solv_disp2state(self.data.structure, tstep)
+
+    def cleanup_timestep_info(self):
+        if len(self.data.structure.timestep_info) > 1:
+            # copy last info to first
+            self.data.structure.timestep_info[0] = self.data.structure.timestep_info[-1].copy()
+            # delete all the rest
+            while len(self.data.structure.timestep_info) - 1:
+                del self.data.structure.timestep_info[-1]
+        self.data.ts = 0
+
+    def del_timestep_info(self):
+        if len(self.data.structure.timestep_info) >= 1:
+            # delete all the rest
+            while len(self.data.structure.timestep_info)-1 :
+                del self.data.structure.timestep_info[-1]
+        self.data.ts = 0
+        self.data.structure.timestep_info[self.data.ts].for_pos[0:3] = self.settings['initial_position']
+
 
 

sharpy/structure/utils/xbeamlib.py

@@ -35,6 +35,7 @@ class Xbopts(ct.Structure):
                 ("gravity_dir_z", ct.c_double),
                 ("balancing", ct.c_bool),
                 ("relaxation_factor", ct.c_double),
+                ("load_ramping_conv", ct.c_bool)
                 ]
 
     def __init__(self):
@@ -59,6 +60,7 @@ def __init__(self):
         self.gravity_dir_z = ct.c_double(1.0)
         self.balancing = ct.c_bool(False)
         self.relaxation_factor = ct.c_double(0.3)
+        self.load_ramping_conv = ct.c_bool(False)
 
 
 xbeamlib = ct_utils.import_ctypes_lib(xbeam_lib_path, 'libxbeam')
@@ -67,7 +69,7 @@ def __init__(self):
 doubleP = ct.POINTER(ct.c_double)
 intP = ct.POINTER(ct.c_int)
 charP = ct.POINTER(ct.c_char_p)
-
+boolP = ct.POINTER(ct.c_bool)
 
 def cbeam3_solv_nlnstatic(beam, settings, ts):
     """@brief Python wrapper for f_cbeam3_solv_nlnstatic
@@ -81,6 +83,8 @@ def cbeam3_solv_nlnstatic(beam, settings, ts):
     n_mass = ct.c_int(beam.n_mass)
     n_stiff = ct.c_int(beam.n_stiff)
 
+
+
     xbopts = Xbopts()
     xbopts.PrintInfo = ct.c_bool(settings['print_info'])
     xbopts.Solution = ct.c_int(112)
@@ -95,6 +99,7 @@ def cbeam3_solv_nlnstatic(beam, settings, ts):
     xbopts.gravity_dir_y = ct.c_double(gravity_vector[1])
     xbopts.gravity_dir_z = ct.c_double(gravity_vector[2])
     xbopts.relaxation_factor = ct.c_double(settings['relaxation_factor'])
+    xbopts.load_ramping_conv = ct.c_bool(settings['load_ramping_conv'])
 
     # here we only need to set the flags at True, all the forces are follower
     xbopts.FollowerForce = ct.c_bool(True)
@@ -127,6 +132,9 @@ def cbeam3_solv_nlnstatic(beam, settings, ts):
                             beam.timestep_info[ts].steady_applied_forces.ctypes.data_as(doubleP),
                             beam.timestep_info[ts].gravity_forces.ctypes.data_as(doubleP)
                             )
+
+#   print(xbopts.load_ramping_conv)
+    settings['load_ramping_conv']=xbopts.load_ramping_conv
 
 
 def cbeam3_loads(beam, ts):