From 6c58c4ee1e2272875a33059471fac4b6f477af43 Mon Sep 17 00:00:00 2001 From: Justin Ngan <75423952+justinngan92@users.noreply.github.com> Date: Wed, 16 Apr 2025 22:58:07 +1000 Subject: [PATCH 1/3] resolve numpy 2.0 compatibility issue * version that works under numpy 2.0 * Not yet check with numpy 1.0 or with other dependency combinations (e.g., scipy) --- src/ospgrillage/__init__.py | 2 +- src/ospgrillage/load.py | 74 ++--- src/ospgrillage/members.py | 10 +- src/ospgrillage/mesh.py | 6 +- src/ospgrillage/osp_grillage.py | 350 ++++++++++++++++++------ src/ospgrillage/postprocessing.py | 21 +- src/ospgrillage/{static.py => utils.py} | 0 tests/test_load.py | 64 ++++- 8 files changed, 385 insertions(+), 142 deletions(-) rename src/ospgrillage/{static.py => utils.py} (100%) diff --git a/src/ospgrillage/__init__.py b/src/ospgrillage/__init__.py index 4bc2fd02..f8330e82 100644 --- a/src/ospgrillage/__init__.py +++ b/src/ospgrillage/__init__.py @@ -1,4 +1,4 @@ -from ospgrillage.static import * +from ospgrillage.utils import * from ospgrillage.mesh import * from ospgrillage.load import * from ospgrillage.material import * diff --git a/src/ospgrillage/load.py b/src/ospgrillage/load.py index ffe5b10f..d22b75a1 100644 --- a/src/ospgrillage/load.py +++ b/src/ospgrillage/load.py @@ -667,35 +667,32 @@ def _define_patch_edge_lines(self): # create equation of plane from four straight lines # create interpolate object f - p = np.array( - [ - [self.load_point_1.p, self.load_point_2.p], - [self.load_point_3.p, self.load_point_4.p], - ] - ) + x = np.array( [ - [self.load_point_1.x, self.load_point_2.x], - [self.load_point_3.x, self.load_point_4.x], + self.load_point_1.x, + self.load_point_2.x, + self.load_point_3.x, + self.load_point_4.x, ] ) z = np.array( [ - [self.load_point_1.z, self.load_point_2.z], - [self.load_point_3.z, self.load_point_4.z], + self.load_point_1.z, + self.load_point_2.z, + self.load_point_3.z, + self.load_point_4.z, ] ) - - # create function to get interpolation of p - self.patch_mag_interpolate = interpolate.interp2d(x, z, p) - mod_list = [ls for ls in self.point_list if ls is not None] - self.patch_min_dim = min( + p = np.array( [ - get_distance(p1, p2) - for (p1, p2) in zip(mod_list, mod_list[1:] + [mod_list[0]]) - if all([p1 is not None, p2 is not None]) + self.load_point_1.p, + self.load_point_2.p, + self.load_point_3.p, + self.load_point_4.p, ] ) + elif self.load_point_4 is None: # update line 3 self.line_3 = LineLoading( @@ -704,35 +701,28 @@ def _define_patch_edge_lines(self): # TODO create equation of plane from 3 points # create interpolate object f - p = np.array( - [ - [self.load_point_1.p, self.load_point_2.p], - [self.load_point_3.p, self.load_point_4.p], - ] - ) x = np.array( [ - [self.load_point_1.x, self.load_point_2.x], - [self.load_point_3.x, self.load_point_4.x], + self.load_point_1.x, + self.load_point_2.x, + self.load_point_3.x, ] ) z = np.array( [ - [self.load_point_1.z, self.load_point_2.z], - [self.load_point_3.z, self.load_point_4.z], + self.load_point_1.z, + self.load_point_2.z, + self.load_point_3.z, ] ) - - # create function to get interpolation of p - self.patch_mag_interpolate = interpolate.interp2d(x, z, p) - mod_list = [ls for ls in self.point_list if ls is not None] - self.patch_min_dim = min( + p = np.array( [ - get_distance(p1, p2) - for (p1, p2) in zip(mod_list, mod_list[1:] + [mod_list[0]]) - if all([p1 is not None, p2 is not None]) + self.load_point_1.p, + self.load_point_2.p, + self.load_point_3.p, ] ) + elif self.load_point_8 is not None: # TODO # point 1 2 3 @@ -745,6 +735,18 @@ def _define_patch_edge_lines(self): raise ValueError( "Patch load points for patch load {} not valid".format(self.name) ) + # create function to get interpolation of p + self.patch_mag_interpolate = interpolate.LinearNDInterpolator( + list(zip(x, z)), p + ) + mod_list = [ls for ls in self.point_list if ls is not None] + self.patch_min_dim = min( + [ + get_distance(p1, p2) + for (p1, p2) in zip(mod_list, mod_list[1:] + [mod_list[0]]) + if all([p1 is not None, p2 is not None]) + ] + ) # --------------------------------------------------------------------------------------------------------------- diff --git a/src/ospgrillage/members.py b/src/ospgrillage/members.py index 274bf81b..fbba51e5 100644 --- a/src/ospgrillage/members.py +++ b/src/ospgrillage/members.py @@ -9,6 +9,8 @@ # if TYPE_CHECKING: # from material import Material +g = 9.81 +mm = 1e3 def create_section(**kwargs): @@ -233,7 +235,7 @@ def __init__( self.section_command_flag = False # # determine mass from section area and material density - self.mass = self.section.A * self.material.density + self.mass = self.section.A * self.material.density * mm # kg/mm elif any( [ self.section.op_ele_type == "ShellMITC4", @@ -423,7 +425,7 @@ def get_element_command_str( ele_str = None if self.section.op_ele_type == "ElasticTimoshenkoBeam": section_input = self.get_member_prop_arguments(ele_width) - ele_str = 'ops.element("{type}", {tag}, *{node_tag_list}, *{memberprop}, {transftag}, {mass})\n'.format( + ele_str = 'ops.element("{type}", {tag}, *{node_tag_list}, *{memberprop}, {transftag},"-mass", {mass})\n'.format( type=self.section.op_ele_type, tag=ele_tag, node_tag_list=node_tag_list, @@ -433,7 +435,7 @@ def get_element_command_str( ) elif self.section.op_ele_type == "elasticBeamColumn": section_input = self.get_member_prop_arguments(ele_width) - ele_str = 'ops.element("{type}", {tag}, *{node_tag_list}, *{memberprop}, {transftag}, {mass})\n'.format( + ele_str = 'ops.element("{type}", {tag}, *{node_tag_list}, *{memberprop}, {transftag}, "-mass", {mass})\n'.format( type=self.section.op_ele_type, tag=ele_tag, node_tag_list=node_tag_list, @@ -442,7 +444,7 @@ def get_element_command_str( mass=self.mass, ) elif self.section.op_ele_type == "nonlinearBeamColumn": - ele_str = 'ops.element("{type}",{tag},*{node_tag_list},{num_int_pt},{sectag},{transftag},{mass})\n'.format( + ele_str = 'ops.element("{type}",{tag},*{node_tag_list},{num_int_pt},{sectag},{transftag},"-mass", {mass})\n'.format( type=self.section.op_ele_type, tag=ele_tag, node_tag_list=node_tag_list, diff --git a/src/ospgrillage/mesh.py b/src/ospgrillage/mesh.py index 562866ae..0f30952e 100644 --- a/src/ospgrillage/mesh.py +++ b/src/ospgrillage/mesh.py @@ -7,7 +7,7 @@ """ import math -from ospgrillage.static import * +from ospgrillage.utils import * from collections import namedtuple @@ -1486,7 +1486,7 @@ def _get_geo_transform_tag(self, ele_nodes: list, offset=None): node_i = self.node_spec[ele_nodes[0]]["coordinate"] node_j = self.node_spec[ele_nodes[1]]["coordinate"] vxz = self._get_vector_xz(node_i, node_j) - vxz = [np.round(num, decimals=self.decimal_lim) for num in vxz] + vxz = [np.round(num, decimals=self.decimal_lim).tolist() for num in vxz] tag_value = self.transform_dict.setdefault( repr(vxz) + "|" + repr(offset), self.transform_counter + 1 ) @@ -2178,7 +2178,7 @@ def _get_geo_transform_tag(self, ele_nodes, offset=None): global_offset_i = [a + b for a, b in zip(node_i, local_offset)] global_offset_j = [a - b for a, b in zip(node_j, local_offset)] global_offset = [global_offset_i, global_offset_j] - vxz = [np.round(num, decimals=self.decimal_lim) for num in vxz] + vxz = [np.round(num, decimals=self.decimal_lim).tolist() for num in vxz] tag_value = self.transform_dict.setdefault( repr(vxz) + "|" + repr(global_offset), self.transform_counter + 1 ) diff --git a/src/ospgrillage/osp_grillage.py b/src/ospgrillage/osp_grillage.py index 3cf8c128..ce5b63ef 100644 --- a/src/ospgrillage/osp_grillage.py +++ b/src/ospgrillage/osp_grillage.py @@ -339,6 +339,26 @@ def _create_mesh(self, **kwargs): print("Meshing complete") return mesh_obj + def _write_imports(self): + """write the import functions""" + with open(self.filename, "w") as file_handle: + # create py file or overwrite existing + # writing headers and description at top of file + file_handle.write( + "# Grillage generator wizard\n# Model name: {}\n".format( + self.model_name + ) + ) + # time + now = datetime.now() + dt_string = now.strftime("%d/%m/%Y %H:%M:%S") + file_handle.write("# Constructed on:{}\n".format(dt_string)) + # necessary imports + file_handle.write( + "import numpy as np\nimport math\nimport openseespy.opensees as ops" + "\nimport vfo.vfo as opsplt\n" + ) + # interface function def create_osp_model(self, pyfile: bool = False): """ @@ -351,23 +371,7 @@ def create_osp_model(self, pyfile: bool = False): self.pyfile = pyfile # if output mode, create the py file if self.pyfile: - with open(self.filename, "w") as file_handle: - # create py file or overwrite existing - # writing headers and description at top of file - file_handle.write( - "# Grillage generator wizard\n# Model name: {}\n".format( - self.model_name - ) - ) - # time - now = datetime.now() - dt_string = now.strftime("%d/%m/%Y %H:%M:%S") - file_handle.write("# Constructed on:{}\n".format(dt_string)) - # necessary imports - file_handle.write( - "import numpy as np\nimport math\nimport openseespy.opensees as ops" - "\nimport vfo.vfo as opsplt\n" - ) + self._write_imports() self._write_op_model() # run model generation in OpenSees or write generation command to py file @@ -385,10 +389,9 @@ def _run_mesh_generation(self): # write / execute model commands self._write_op_node(self.Mesh_obj) # write node() commands self._write_op_fix(self.Mesh_obj) # write fix() command for support nodes + self._write_mass() self._write_geom_transf(self.Mesh_obj) # x dir members - # write / execute variable definition command - # write / execute material and sections for mat_str in self.material_command_list: if self.pyfile: @@ -409,7 +412,6 @@ def _run_mesh_generation(self): self.model_command_list.append(sec_str) # write /execute element commands - for ele_tag, ele_str in self.element_command_list.items(): if self.pyfile: with open(self.filename, "a") as file_handle: @@ -579,6 +581,24 @@ def _write_op_fix(self, mesh_obj): eval(fix_str) self.model_command_list.append(fix_str) + # TODO + def _write_mass(self): + """write or evaluates the mass commands""" + + # d = {node number: [1,2,3,0,0,0] } # dx dy dz, rx, ry, rz + # + mass_dict = {} + for node, mass_list in mass_dict.items(): + mass_str = "ops.mass({nodetag},*{mass_list})".format( + nodetag=node, mass_list=mass_list + ) + if self.pyfile: + with open(self.filename, "a") as file_handle: + file_handle.write(mass_str) + else: + eval(mass_str) + self.model_command_list.append(mass_str) + def _write_equal_dof(self, node_tag_list: list, dof: list = None): """ Write OpenseesPy's equalDOF command. @@ -719,9 +739,9 @@ def _create_standard_element_list(self): self.common_grillage_element_z_group.update({key: val}) # populate start edge and end edge entries self.common_grillage_element_z_group[self.common_grillage_element_keys[4]] = [0] - self.common_grillage_element_z_group[ - self.common_grillage_element_keys[5] - ] = list(range(1, self.Mesh_obj.global_edge_count)) + self.common_grillage_element_z_group[self.common_grillage_element_keys[5]] = ( + list(range(1, self.Mesh_obj.global_edge_count)) + ) self.common_grillage_element_z_group[self.common_grillage_element_keys[6]] = [ 0 ] # proxy 0 for set_member() loop @@ -1690,9 +1710,9 @@ def _assign_patch_load(self, patch_load_obj: PatchLoading) -> list: p_list = [] for tag in nodes: coord = self.Mesh_obj.node_spec[tag]["coordinate"] - p = patch_load_obj.patch_mag_interpolate(coord[0], coord[2])[ - 0 - ] # object function returns array like + p = patch_load_obj.patch_mag_interpolate( + coord[0], coord[2] + ) # object function returns array like p_list.append(LoadPoint(coord[0], coord[1], coord[2], p)) # get centroid of patch on grid xc, yc, zc = get_patch_centroid(p_list) @@ -1764,16 +1784,16 @@ def _assign_patch_load(self, patch_load_obj: PatchLoading) -> list: ) # object function returns array like # p is array object, extract p_list.append( - LoadPoint(int_point[0], int_point[1], int_point[2], p[0]) + LoadPoint(int_point[0], int_point[1], int_point[2], p) if int_point != [] else [] ) # loop each node in grid points for items in node_in_grid: coord = self.Mesh_obj.node_spec[items]["coordinate"] - p = patch_load_obj.patch_mag_interpolate(coord[0], coord[2])[ - 0 - ] # object function returns array like + p = patch_load_obj.patch_mag_interpolate( + coord[0], coord[2] + ) # object function returns array like p_list.append(LoadPoint(coord[0], coord[1], coord[2], p)) # Loop each p_list object to find duplicates if any, remove duplicate for count, point in enumerate(p_list): @@ -1945,9 +1965,9 @@ def add_load_case( "load_factor": load_factor, } list_of_incr_load_case_dict.append(increment_load_case_dict) - self.moving_load_case_dict[ - moving_load_obj.name - ] = list_of_incr_load_case_dict + self.moving_load_case_dict[moving_load_obj.name] = ( + list_of_incr_load_case_dict + ) if self.diagnostics: print("Moving load case: {} created".format(moving_load_obj.name)) @@ -1966,7 +1986,7 @@ def analyze(self, **kwargs) -> None: * all (`bool`): If True, runs all load cases. If not provided, default to True. * load_case ('list' or 'str'): String or list of name strings for selected load case to be analyzed. * set_verbose(`bool`): If True, incremental load case report is not printed to terminal (default True) - + * analysis_type ('str'): The type of analysis. Default is "Static". :except: raise ValueError if missing arguments for either load_case=, or all= """ @@ -1975,7 +1995,8 @@ def analyze(self, **kwargs) -> None: # get run options from kwargs all_flag = True # Default true selected_load_case: list = kwargs.get("load_case", None) # - + analysis_type = kwargs.get("analysis_type", "Static") + step = kwargs.get("step", 1) # default 1 # check if any load cases are defined if self.load_case_list == [] and self.moving_load_case_dict == {}: raise Exception("No load cases were defined") @@ -2034,6 +2055,8 @@ def analyze(self, **kwargs) -> None: load_factor = load_case_dict["load_factor"] load_case_analysis = Analysis( analysis_name=load_case_obj.name, + analysis_type=analysis_type, + step=step, ops_grillage_name=self.model_name, pyfile=self.pyfile, time_series_counter=self.global_time_series_counter, @@ -2042,6 +2065,7 @@ def analyze(self, **kwargs) -> None: ele_counter=self.Mesh_obj.element_counter, constraint_type=self.constraint_type, load_case=load_case_obj, + time_increment=kwargs.get("time_increment", 0.01), ) load_case_analysis.add_load_command(load_command, load_factor=load_factor) # run the Analysis object, collect results, and store Analysis object in the list for Analysis load case @@ -2070,6 +2094,8 @@ def analyze(self, **kwargs) -> None: load_factor = load_case_dict["load_factor"] incremental_analysis = Analysis( analysis_name=load_case_obj.name, + analysis_type=analysis_type, + step=step, ops_grillage_name=self.model_name, pyfile=self.pyfile, time_series_counter=self.global_time_series_counter, @@ -2078,6 +2104,7 @@ def analyze(self, **kwargs) -> None: ele_counter=self.Mesh_obj.element_counter, constraint_type=self.constraint_type, load_case=load_case_obj, + time_increment=kwargs.get("time_increment", 0.01), ) incremental_analysis.add_load_command( load_command, load_factor=load_factor @@ -2418,6 +2445,79 @@ def clear_load_cases(self, **kwargs): # remove all results self.results = Results(self.Mesh_obj) # reset results + @staticmethod + def get_MCK(): + """Returns the mass stiffness and damping matrices. + Note model needs to be pyfile=False mode + + ref:https://portwooddigital.com/2020/05/17/gimme-all-your-damping-all-your-mass-and-stiffness-too/ + """ + ops.wipeAnalysis() + ops.system("FullGeneral") + ops.analysis("Transient") + + # Mass + ops.integrator("GimmeMCK", 1.0, 0.0, 0.0) + ops.analyze(1, 0.0) + + # Number of equations in the model + N = ops.systemSize() # Has to be done after analyze + + M = ops.printA("-ret") # Or use ops.printA('-file','M.out') + M = np.array(M) # Convert the list to an array + M.shape = (N, N) # Make the array an NxN matrix + + # Stiffness + ops.integrator("GimmeMCK", 0.0, 0.0, 1.0) + ops.analyze(1, 0.0) + K = ops.printA("-ret") + K = np.array(K) + K.shape = (N, N) + + # Damping + ops.integrator("GimmeMCK", 0.0, 1.0, 0.0) + ops.analyze(1, 0.0) + C = ops.printA("-ret") + C = np.array(C) + C.shape = (N, N) + + massDOFs = [] + for nd in ops.getNodeTags(): + + for j in range(6): # NDF is number of DOFs/node + if ops.nodeMass(nd, j + 1) > 0.0: + massDOFs.append(ops.nodeDOFs(nd)[j]) + + return M, C, K + + @staticmethod + def store_state(): + """Use for transient analysis to store previous state""" + save = {} + for tag in ops.getNodeTags(): + d = ops.nodeDisp(tag) + v = ops.nodeVel(tag) + a = ops.nodeAccel(tag) + save[tag] = {"d": d, "v": v, "a": a} + return save + + @staticmethod + def set_previous_state(save): + """Set the disp, vel, and acc of all nodes in the last analysis state to the current analysis""" + for node, val in save.items(): + d = val["d"] + v = val["v"] + a = val["a"] + + for i, magnitude in enumerate(d): + ops.setNodeDisp(node, i + 1, magnitude, "-commit") + for i, magnitude in enumerate(v): + ops.setNodeVel(node, i + 1, magnitude, "-commit") + for i, magnitude in enumerate(a): + ops.setNodeAccel(node, i + 1, magnitude, "-commit") + + ... + # --------------------------------------------------------------------------------------------------------------------- class Analysis: @@ -2450,6 +2550,7 @@ def __init__( step: int = 1, **kwargs, ): + self.analysis_name = analysis_name self.ops_grillage_name = ops_grillage_name self.time_series_tag = None @@ -2459,7 +2560,28 @@ def __init__( self.pyfile = pyfile self.analysis_file_name = ( self.analysis_name + "of" + self.ops_grillage_name + ".py" - ) # py file name + ) + + # default newmark parameters for Average Accelerations + gamma = 0.5 + beta = 0.25 + time_increment = kwargs.get("time_increment", 0.01) + if analysis_type == "Transient": + if kwargs.get("linear_acceleration"): + beta = 1 / 6 + + # Preset dict for different analysis + analysis_arguments = { + "Static": { + "analyse": "ops.analyze({})\n".format(self.step), + "integrator": "ops.integrator('LoadControl', 1)\n", + }, + "Transient": { + "analyse": "ops.analyze({},{})\n".format(self.step, time_increment), + "integrator": "ops.integrator('Newmark', {},{})\n".format(gamma, beta), + }, + } + # list recording load commands, time series and pattern for the input load case self.load_cases_dict_list = ( [] @@ -2471,6 +2593,8 @@ def __init__( self.constraint_type = kwargs.get("constraint_type", "Plain") # Default plain # Variables recording results of analysis self.node_disp = dict() # key node tag, val list of dof + self.node_vel = dict() # key node tag, val list of dof + self.node_accel = dict() # key node tag, val list of dof self.ele_force = ( dict() ) # key ele tag, val list of forces on nodes of ele[ order according to ele tag] @@ -2489,9 +2613,11 @@ def __init__( type=self.constraint_type ) # default plain self.algorithm_command = "ops.algorithm('Linear')\n" # default linear - self.analyze_command = "ops.analyze({})\n".format(self.step) # default 1 step + self.analyze_command = analysis_arguments[analysis_type][ + "analyse" + ] # default 1 step self.analysis_command = 'ops.analysis("{}")\n'.format(analysis_type) - self.intergrator_command = "ops.integrator('LoadControl', 1)\n" + self.intergrator_command = analysis_arguments[analysis_type]["integrator"] self.sensitivity_integrator_command = "ops." self.mesh_node_counter = node_counter # set node counter based on current Mesh self.mesh_ele_counter = ele_counter # set ele counter based on current Mesh @@ -2625,7 +2751,11 @@ def extract_grillage_responses(self): # first loop extract node displacements for node_tag in ops.getNodeTags(): disp_list = ops.nodeDisp(node_tag) + vel_list = ops.nodeVel(node_tag) + accel_list = ops.nodeAccel(node_tag) self.node_disp.setdefault(node_tag, disp_list) + self.node_vel.setdefault(node_tag, vel_list) + self.node_accel.setdefault(node_tag, accel_list) # loop through all elements in Mesh, extract local forces for ele_tag in ops.getEleTags(): @@ -2659,13 +2789,29 @@ def __init__(self, mesh_obj: Mesh): self.mesh_obj = mesh_obj # coordinates for dimensions self.displacement_component = [ - "dx", - "dy", - "dz", + "x", + "y", + "z", "theta_x", "theta_y", "theta_z", ] + self.vel_component = [ + "dx", + "dy", + "dz", + "dtheta_x", + "dtheta_y", + "dtheta_z", + ] + self.acc_component = [ + "ddx", + "ddy", + "ddz", + "ddtheta_x", + "ddtheta_y", + "ddtheta_z", + ] self.force_component = [ "Vx_i", "Vy_i", @@ -2720,6 +2866,8 @@ def extract_analysis( if analysis_obj: # compile ele forces for each node node_disp = analysis_obj.node_disp + node_vel = analysis_obj.node_vel + node_accel = analysis_obj.node_accel ele_force_dict = dict.fromkeys( list(ops.getEleTags()) ) # dict key is element tag, value is ele force from @@ -2746,7 +2894,13 @@ def extract_analysis( ele_nodes_dict.update({ele_num: ele_nodes}) self.basic_load_case_record_global_forces.setdefault( analysis_obj.analysis_name, - [node_disp, global_ele_force_dict, ele_nodes_dict], + [ + node_disp, + node_vel, + node_accel, + global_ele_force_dict, + ele_nodes_dict, + ], ) # if moving load, input is a list of analysis obj elif list_of_inc_analysis: @@ -2758,6 +2912,8 @@ def extract_analysis( for inc_analysis_obj in list_of_inc_analysis: # compile ele forces for each node node_disp = inc_analysis_obj.node_disp + node_vel = inc_analysis_obj.node_vel + node_accel = inc_analysis_obj.node_accel ele_force_dict = dict.fromkeys(list(ops.getEleTags())) ele_force_global_dict = dict.fromkeys(list(ops.getEleTags())) ele_nodes_dict = dict.fromkeys(list(ops.getEleTags())) @@ -2776,7 +2932,13 @@ def extract_analysis( ele_force_global_dict.update({ele_num: ele_forces}) inc_load_case_global_force_record.setdefault( inc_analysis_obj.analysis_name, - [node_disp, ele_force_global_dict, ele_nodes_dict], + [ + node_disp, + node_vel, + node_accel, + ele_force_global_dict, + ele_nodes_dict, + ], ) self.moving_load_case_record.append(inc_load_case_record) @@ -2793,6 +2955,8 @@ def compile_data_array(self, local_force_option=True, main_ele_tags=None): # Sort data for dataArrays # for basic load case {loadcasename:[{1:,2:...},{1:,2:...}], ... , loadcasename:[{1:,2:...},{1:,2:...} } basic_node_disp_list = [] + basic_node_vel_list = [] + basic_node_accel_list = [] basic_load_case_coord = [] basic_ele_force_list = [] extracted_ele_nodes_list = False # a 2D array of ele node i and ele node j @@ -2814,33 +2978,39 @@ def compile_data_array(self, local_force_option=True, main_ele_tags=None): basic_node_disp_list.append( [a for a in list(resp_list_of_2_dict[0].values())] ) # list index 0 is disp + basic_node_vel_list.append( + [a for a in list(resp_list_of_2_dict[1].values())] + ) + basic_node_accel_list.append( + [a for a in list(resp_list_of_2_dict[2].values())] + ) # extract force basic_ele_force_list.append( - [a for a in list(resp_list_of_2_dict[1].values())] + [a for a in list(resp_list_of_2_dict[3].values())] ) # list index 1 is force # extract based on element type base_ele_force_list_beam.append( [ a - for a in list(resp_list_of_2_dict[1].values()) + for a in list(resp_list_of_2_dict[3].values()) if len(a) == len(self.force_component) ] ) base_ele_force_list_shell.append( [ a - for a in list(resp_list_of_2_dict[1].values()) + for a in list(resp_list_of_2_dict[3].values()) if len(a) == len(self.force_component_shell) ] ) if not extracted_ele_nodes_list: ele_nodes_list = list( - resp_list_of_2_dict[2].values() + resp_list_of_2_dict[4].values() ) # list index 2 is ele nodes variable extracted_ele_nodes_list = ( True # set to true, only extract if its the first time extracting ) - ele_tag = list(resp_list_of_2_dict[2].keys()) + ele_tag = list(resp_list_of_2_dict[4].keys()) # for section forces of each element # Coordinate of Load Case dimension basic_load_case_coord.append(load_case_name) @@ -2858,19 +3028,24 @@ def compile_data_array(self, local_force_option=True, main_ele_tags=None): basic_node_disp_list.append( [a for a in list(inc_resp_list_of_2_dict[0].values())] ) - + basic_node_vel_list.append( + [a for a in list(inc_resp_list_of_2_dict[1].values())] + ) + basic_node_accel_list.append( + [a for a in list(inc_resp_list_of_2_dict[2].values())] + ) if local_force_option: basic_ele_force_list.append( [ a - for a in list(inc_resp_list_of_2_dict[1].values()) + for a in list(inc_resp_list_of_2_dict[3].values()) if len(a) == len(self.force_component) ] ) else: # global force basic_ele_force_list.append( - [a for a in list(inc_resp_list_of_2_dict[1].values())] + [a for a in list(inc_resp_list_of_2_dict[3].values())] ) # lists for shell model output @@ -2878,8 +3053,8 @@ def compile_data_array(self, local_force_option=True, main_ele_tags=None): [ a for key, a in zip( - list(inc_resp_list_of_2_dict[1].keys()), - list(inc_resp_list_of_2_dict[1].values()), + list(inc_resp_list_of_2_dict[3].keys()), + list(inc_resp_list_of_2_dict[3].values()), ) if len(a) == len(self.force_component) if key < main_ele_tags @@ -2888,7 +3063,7 @@ def compile_data_array(self, local_force_option=True, main_ele_tags=None): base_ele_force_list_shell.append( [ a - for a in list(inc_resp_list_of_2_dict[1].values()) + for a in list(inc_resp_list_of_2_dict[3].values()) if len(a) == len(self.force_component_shell) ] ) @@ -2896,11 +3071,13 @@ def compile_data_array(self, local_force_option=True, main_ele_tags=None): # inc_moving_load_case_coord.append(increment_load_case_name) basic_load_case_coord.append(increment_load_case_name) if not extracted_ele_nodes_list: - ele_nodes_list = list(inc_resp_list_of_2_dict[2].values()) - ele_tag = list(inc_resp_list_of_2_dict[2].keys()) + ele_nodes_list = list(inc_resp_list_of_2_dict[4].values()) + ele_tag = list(inc_resp_list_of_2_dict[4].keys()) extracted_ele_nodes_list = True # convert to np array format - basic_array = np.array(basic_node_disp_list, dtype=object) + basic_array_disp = np.array(basic_node_disp_list, dtype=object) + basic_array_vel = np.array(basic_node_vel_list, dtype=object) + basic_array_accel = np.array(basic_node_accel_list, dtype=object) force_array = np.array(basic_ele_force_list, dtype=object) ele_array = np.array(ele_nodes_list, dtype=object) @@ -2923,10 +3100,10 @@ def compile_data_array(self, local_force_option=True, main_ele_tags=None): force_array_beam = np.array(base_ele_force_list_beam) # create data array for each basic load case if any, else return - if basic_array.size: + if basic_array_disp.size: # displacement data array - basic_da = xr.DataArray( - data=basic_array, + basic_da_d = xr.DataArray( + data=basic_array_disp, dims=self.dim, coords={ self.dim[0]: basic_load_case_coord, @@ -2934,6 +3111,25 @@ def compile_data_array(self, local_force_option=True, main_ele_tags=None): self.dim[2]: self.displacement_component, }, ) + basic_da_v = xr.DataArray( + data=basic_array_vel, + dims=self.dim, + coords={ + self.dim[0]: basic_load_case_coord, + self.dim[1]: node, + self.dim[2]: self.vel_component, + }, + ) + + basic_da_a = xr.DataArray( + data=basic_array_accel, + dims=self.dim, + coords={ + self.dim[0]: basic_load_case_coord, + self.dim[1]: node, + self.dim[2]: self.acc_component, + }, + ) ele_nodes_beam = xr.DataArray( data=ele_array_beam, @@ -2967,7 +3163,9 @@ def compile_data_array(self, local_force_option=True, main_ele_tags=None): ) result = xr.Dataset( { - "displacements": basic_da, + "displacements": basic_da_d, + "velocity": basic_da_v, + "acceleration": basic_da_a, "forces_beam": force_da_beam, "forces_shell": force_da_shell, "ele_nodes_beam": ele_nodes_beam, @@ -2980,15 +3178,17 @@ def compile_data_array(self, local_force_option=True, main_ele_tags=None): dims=self.dim2, coords={ self.dim2[0]: basic_load_case_coord, - self.dim2[1]: ele_tag - if not local_force_option - else ele_tag_beam, + self.dim2[1]: ( + ele_tag if not local_force_option else ele_tag_beam + ), self.dim2[2]: self.force_component, }, ) result = xr.Dataset( { - "displacements": basic_da, + "displacements": basic_da_d, + "velocity": basic_da_v, + "acceleration": basic_da_a, "forces": force_da_beam, "ele_nodes": ele_nodes_beam, } @@ -3100,23 +3300,7 @@ def create_osp_model(self, pyfile=False): self.pyfile = pyfile if self.pyfile: - with open(self.filename, "w") as file_handle: - # create py file or overwrite existing - # writing headers and description at top of file - file_handle.write( - "# Grillage generator wizard\n# Model name: {}\n".format( - self.model_name - ) - ) - # time - now = datetime.now() - dt_string = now.strftime("%d/%m/%Y %H:%M:%S") - file_handle.write("# Constructed on:{}\n".format(dt_string)) - # write imports - file_handle.write( - "import numpy as np\nimport math\nimport openseespy.opensees as ops" - "\nimport vfo.vfo as opsplt\n" - ) + self._write_imports() # model() command self._write_op_model() # create grillage mesh object + beam element groups @@ -3129,11 +3313,11 @@ def create_osp_model(self, pyfile=False): file_handle.write(ele_str) else: eval(ele_str) + self.model_command_list.append(ele_str) # create rigid link command self._write_rigid_link() # create the result file for the Mesh object self.results = Results(self.Mesh_obj) - # flag # overwrites base class for beam element grillage - specific for Shell model def _create_standard_element_list(self): diff --git a/src/ospgrillage/postprocessing.py b/src/ospgrillage/postprocessing.py index 6e9a598f..fa648e9f 100644 --- a/src/ospgrillage/postprocessing.py +++ b/src/ospgrillage/postprocessing.py @@ -15,7 +15,7 @@ # if TYPE_CHECKING: from ospgrillage.load import ShapeFunction -from ospgrillage.static import solve_zeta_eta +from ospgrillage.utils import solve_zeta_eta import vfo.vfo as opsplt @@ -289,7 +289,7 @@ def plot_force( zz = [nodes[n]["coordinate"][2] for n in ele_node.values] # use ops_vis module to get force distribution on element s, al = opsv.section_force_distribution_3d( - ex=xx, ey=yy, ez=zz, pl=ele_components + ex=xx, ey=yy, ez=zz, pl=ele_components #TODO check inputs with opsv package ) # plot element force component @@ -358,7 +358,7 @@ def plot_defo( # check if component is provided, else default to dis_comp = component if component is None: - dis_comp = "dy" # default to dy + dis_comp = "y" # default to dy fig, ax = plt.subplots() # get all node information @@ -411,13 +411,16 @@ def __init__(self, grillage, result): # Union[xr.DataArray, xr.Dataset] self.shape_function_obj = ShapeFunction() def get_arbitrary_displacements( - self, point: list, shape_function_type: str = "linear" + self, point: list, shape_function_type: str = "linear", component: str = "y" ): - """Returns displacement values (translation and rotational) + """Returns displacement (translation and rotational) from an arbitrary point by interpolation param point: list of coordinate. Default three elements [x,y=0,z] type point: list - + param component: Displacement component. Default "y" + type component: str + param shape_function_type: The shape function for interpolation. Default "linear" + type shape_function_type: str """ node_displacements = [] node_coordinate = [] @@ -426,12 +429,14 @@ def get_arbitrary_displacements( nodes, grid_number = self.grillage._get_point_load_nodes( point=point ) # list of nodes + if nodes is None: + raise ValueError("Point is outside bridge mesh") # get results of each node of four nodes for node in nodes: node_displacements.append( - self.result.displacements.sel( - Component="dy", + self.result.sel( + Component=component, Node=node, ) .to_numpy() diff --git a/src/ospgrillage/static.py b/src/ospgrillage/utils.py similarity index 100% rename from src/ospgrillage/static.py rename to src/ospgrillage/utils.py diff --git a/tests/test_load.py b/tests/test_load.py index 58872f19..e7e8ac15 100644 --- a/tests/test_load.py +++ b/tests/test_load.py @@ -89,14 +89,21 @@ def test_point_load_getter( ULS_DL.add_load(Single) # ch example_bridge.add_load_case(ULS_DL) og.ops.wipe() + # + # assert example_bridge.load_case_list[0]["load_command"] == [ + # "ops.load(12, *[0, 0.6075807082987842, 0, 0.37389582049155967, 0, 0.34166943132701877])\n", + # "ops.load(17, *[0, 1.4724192917012129, 0, 0.9061041795084391, 0, -0.613915767971676])\n", + # "ops.load(18, *[0, 12.685458513118162, 0, -3.6244167180337583, 0, -5.289120462525217])\n", + # "ops.load(13, *[0, 5.234541486881841, 0, -1.4955832819662394, 0, 2.943613562202012])\n", + # ] + reference_command = example_bridge.load_case_list[0]["load_command"] - assert example_bridge.load_case_list[0]["load_command"] == [ - "ops.load(12, *[0, 0.6075807082987842, 0, 0.37389582049155967, 0, 0.34166943132701877])\n", - "ops.load(17, *[0, 1.4724192917012129, 0, 0.9061041795084391, 0, -0.613915767971676])\n", - "ops.load(18, *[0, 12.685458513118162, 0, -3.6244167180337583, 0, -5.289120462525217])\n", - "ops.load(13, *[0, 5.234541486881841, 0, -1.4955832819662394, 0, 2.943613562202012])\n", - ] + assert reference_command + assert example_bridge.load_case_list[0]["load_command"] == ['ops.load(12, *[0, np.float64(0.6075807082987873), 0, np.float64(0.3738958204915613), 0, np.float64(0.34166943132702027)])\n', + 'ops.load(17, *[0, np.float64(1.4724192917012136), 0, np.float64(0.9061041795084389), 0, np.float64(-0.6139157679716769)])\n', + 'ops.load(18, *[0, np.float64(12.685458513118142), 0, np.float64(-3.624416718033756), 0, np.float64(-5.289120462525214)])\n', + 'ops.load(13, *[0, np.float64(5.234541486881858), 0, np.float64(-1.4955832819662453), 0, np.float64(2.94361356220202)])\n'] # test point load returning None when point (loadpoint) is outside of mesh def test_point_load_outside_straight_mesh(bridge_model_42_negative): @@ -1111,5 +1118,48 @@ def test_compare_shell_beam_analysis(run_beam_model_point_load): og.plot_defo( shell_bridge, result_shell, member="interior_main_beam", option="nodes" ) - pass og.opsv.plot_defo() + + +def test_transient(beam_element_bridge): + """A test for the uncoupled transient analysis portion of VBI basic framework""" + positions = [5, 5, 10, 15] + P = 100 * kN + first_step = True + previous_state = None + + # start moving vehicle through all positions + for x in positions: + # create load case + mid_point_line_loadcase = og.create_load_case(name="VBI") + lp1 = og.create_load_vertex(x=x, y=0, z=3.5, p=P) + mid_point_line_load = og.create_load( + name="unit load", + point1=lp1, + ) + mid_point_line_loadcase.add_load(mid_point_line_load) + + # create bridge instance + beam_bridge = beam_element_bridge + beam_bridge.create_osp_model() + + print(og.ops.eigen(2)) + og.ops.rayleigh(0.0, 0.0, 0.0, 2 * 0.02 / 4) + + # set previous state + if not first_step: + beam_bridge.set_previous_state(previous_state) + first_step = False + + # add load case, run and store + beam_bridge.add_load_case(mid_point_line_loadcase) + beam_bridge.analyze(analysis_type="Transient", step=100, time_increment=0.1) + previous_state = beam_bridge.store_state() + result = beam_bridge.get_results() + print(result) + postprocessor = og.PostProcessor(beam_bridge, result.displacements) + + contactpoints = postprocessor.get_arbitrary_displacements(point=[5, 0, 3.5]) + # do something with results + print(result) + print("Next step") From 7e5e6fa07cc06d2ee91393eef6be9381131eed8c Mon Sep 17 00:00:00 2001 From: jngan Date: Tue, 22 Apr 2025 15:08:02 +1000 Subject: [PATCH 2/3] Add system dependency (BLAS) to yml workflow --- .github/workflows/codecov_test_suite.yml | 3 +++ .github/workflows/pytest.yml | 6 ++++++ 2 files changed, 9 insertions(+) diff --git a/.github/workflows/codecov_test_suite.yml b/.github/workflows/codecov_test_suite.yml index 8a44a119..16e21e97 100644 --- a/.github/workflows/codecov_test_suite.yml +++ b/.github/workflows/codecov_test_suite.yml @@ -42,6 +42,9 @@ jobs: - name: Install package run: python -m pip install -e .[test] + - name: Install system dependencies (BLAS) + run: sudo apt-get update && sudo apt-get install -y libblas3 liblapack3 + - name: Generate Report run: | pip install coverage diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 52407602..ea20602e 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -34,6 +34,9 @@ jobs: - name: Install package run: python -m pip install -e .[test] + - name: Install system dependencies (BLAS) + run: sudo apt-get update && sudo apt-get install -y libblas3 liblapack3 + - name: Test package run: python -m pytest @@ -63,6 +66,9 @@ jobs: - name: Install package run: python -m pip install -e .[test] + - name: Install system dependencies (BLAS) + run: sudo apt-get update && sudo apt-get install -y libblas3 liblapack3 + - name: Generate Report run: | pip install coverage From 78ab02c7b054558874bc3aab8163b5f914f28fd2 Mon Sep 17 00:00:00 2001 From: jngan Date: Tue, 22 Apr 2025 17:06:13 +1000 Subject: [PATCH 3/3] fix tests - improve assertion for strings with parsing and comparing parts --- src/ospgrillage/postprocessing.py | 7 +++-- tests/test_benchmark.py | 21 ++++++++----- tests/test_load.py | 45 +++++++++++++++++++++----- tests/test_osp_grillage.py | 52 ++++++++++++++++++++++++++++--- 4 files changed, 105 insertions(+), 20 deletions(-) diff --git a/src/ospgrillage/postprocessing.py b/src/ospgrillage/postprocessing.py index fa648e9f..a5fbb398 100644 --- a/src/ospgrillage/postprocessing.py +++ b/src/ospgrillage/postprocessing.py @@ -289,7 +289,10 @@ def plot_force( zz = [nodes[n]["coordinate"][2] for n in ele_node.values] # use ops_vis module to get force distribution on element s, al = opsv.section_force_distribution_3d( - ex=xx, ey=yy, ez=zz, pl=ele_components #TODO check inputs with opsv package + ex=xx, + ey=yy, + ez=zz, + pl=ele_components, # TODO check inputs with opsv package ) # plot element force component @@ -439,7 +442,7 @@ def get_arbitrary_displacements( Component=component, Node=node, ) - .to_numpy() + .displacements.to_numpy() .tolist()[0] ) node_coordinate.append(self.grillage.get_nodes(number=node)) diff --git a/tests/test_benchmark.py b/tests/test_benchmark.py index 9842becf..ddb8254d 100644 --- a/tests/test_benchmark.py +++ b/tests/test_benchmark.py @@ -849,11 +849,14 @@ def test_line_load_results(add_analysis_to_simple_grid): "3_Line_Test_Case.csv", ] print(Path.cwd()) - point_load_disp_lusas = pandas.read_csv(Path.cwd().joinpath(*point_lusas_disp_path)) + print(Path(__file__)) + point_load_disp_lusas = pandas.read_csv( + Path(__file__).resolve().parent.parent.joinpath(*point_lusas_disp_path) + ) # get ospg point load case point_load_disp_ospg = all_results["displacements"].sel( Loadcase="Points Test (Global)", - Component=["dx", "dy", "dz", "theta_x", "theta_y", "theta_z"], + Component=["x", "y", "z", "theta_x", "theta_y", "theta_z"], ) # lusas elements are 3 noded beam elemnt, this function extracts only the end nodes (first and third) of the model. @@ -863,23 +866,27 @@ def test_line_load_results(add_analysis_to_simple_grid): ) sorted_zip_ospg_node = sort_array_by_node_mapping( list_of_node=node_ospg, - data_of_node=point_load_disp_ospg.sel(Component="dy").values, + data_of_node=point_load_disp_ospg.sel(Component="y").values, ) - line_load_disp_lusas = pandas.read_csv(Path.cwd().joinpath(*line_lusas_disp_path)) + line_load_disp_lusas = pandas.read_csv( + Path(__file__).resolve().parent.parent.joinpath(*line_lusas_disp_path) + ) lusas_def = reduce_lusas_node_result( pd_data=line_load_disp_lusas["DZ[m]"], node_to_extract_list=node_lusas ) line_load_disp_ospg = all_results["displacements"].sel( Loadcase="Line Test", - Component=["dx", "dy", "dz", "theta_x", "theta_y", "theta_z"], + Component=["x", "y", "z", "theta_x", "theta_y", "theta_z"], ) sorted_zip_ospg_node = sort_array_by_node_mapping( list_of_node=node_ospg, - data_of_node=line_load_disp_ospg.sel(Component="dy").values, + data_of_node=line_load_disp_ospg.sel(Component="y").values, ) # lusas bending z - line_load_force_lusas = pandas.read_csv(Path.cwd().joinpath(*line_lusas_force_path)) + line_load_force_lusas = pandas.read_csv( + Path(__file__).resolve().parent.parent.joinpath(*line_lusas_force_path) + ) single_component_line_lusas = extract_lusas_ele_forces( list_of_ele=a, df_force=line_load_force_lusas, component="My[N.m]" ) diff --git a/tests/test_load.py b/tests/test_load.py index e7e8ac15..94bfb7a9 100644 --- a/tests/test_load.py +++ b/tests/test_load.py @@ -96,14 +96,45 @@ def test_point_load_getter( # "ops.load(18, *[0, 12.685458513118162, 0, -3.6244167180337583, 0, -5.289120462525217])\n", # "ops.load(13, *[0, 5.234541486881841, 0, -1.4955832819662394, 0, 2.943613562202012])\n", # ] - reference_command = example_bridge.load_case_list[0]["load_command"] + # reference_command = example_bridge.load_case_list[0]["load_command"] + # + # assert reference_command + import re + + def parse_load_cmd(cmd): + # Extract node ID and list of floats from string like: ops.load(12, *[0, np.float64(x), 0, ...]) + match = re.match(r"ops\.load\((\d+), \*\[([^\]]+)\]\)", cmd.strip()) + assert match, f"Failed to parse: {cmd}" + node = int(match.group(1)) + values = [ + float(x.split("(")[-1].rstrip(")")) + for x in match.group(2).split(",") + if "np.float64" in x + ] + return node, values + + expected_cmds = [ + (12, [0.6075807082987873, 0.3738958204915613, 0.34166943132702027]), + (17, [1.4724192917012136, 0.9061041795084389, -0.6139157679716769]), + (18, [12.685458513118142, -3.624416718033756, -5.289120462525214]), + (13, [5.234541486881858, -1.4955832819662453, 2.94361356220202]), + ] + + for cmd_str, (exp_node, exp_vals) in zip( + example_bridge.load_case_list[0]["load_command"], expected_cmds + ): + node, values = parse_load_cmd(cmd_str) + assert node == exp_node + for a, b in zip(values, exp_vals): + assert abs(a - b) < 1e-9 # tolerance for float comparison - assert reference_command + # assert example_bridge.load_case_list[0]["load_command"] == [ + # "ops.load(12, *[0, np.float64(0.6075807082987873), 0, np.float64(0.3738958204915613), 0, np.float64(0.34166943132702027)])\n", + # "ops.load(17, *[0, np.float64(1.4724192917012136), 0, np.float64(0.9061041795084389), 0, np.float64(-0.6139157679716769)])\n", + # "ops.load(18, *[0, np.float64(12.685458513118142), 0, np.float64(-3.624416718033756), 0, np.float64(-5.289120462525214)])\n", + # "ops.load(13, *[0, np.float64(5.234541486881858), 0, np.float64(-1.4955832819662453), 0, np.float64(2.94361356220202)])\n", + # ] - assert example_bridge.load_case_list[0]["load_command"] == ['ops.load(12, *[0, np.float64(0.6075807082987873), 0, np.float64(0.3738958204915613), 0, np.float64(0.34166943132702027)])\n', - 'ops.load(17, *[0, np.float64(1.4724192917012136), 0, np.float64(0.9061041795084389), 0, np.float64(-0.6139157679716769)])\n', - 'ops.load(18, *[0, np.float64(12.685458513118142), 0, np.float64(-3.624416718033756), 0, np.float64(-5.289120462525214)])\n', - 'ops.load(13, *[0, np.float64(5.234541486881858), 0, np.float64(-1.4955832819662453), 0, np.float64(2.94361356220202)])\n'] # test point load returning None when point (loadpoint) is outside of mesh def test_point_load_outside_straight_mesh(bridge_model_42_negative): @@ -1157,7 +1188,7 @@ def test_transient(beam_element_bridge): previous_state = beam_bridge.store_state() result = beam_bridge.get_results() print(result) - postprocessor = og.PostProcessor(beam_bridge, result.displacements) + postprocessor = og.PostProcessor(beam_bridge, result) contactpoints = postprocessor.get_arbitrary_displacements(point=[5, 0, 3.5]) # do something with results diff --git a/tests/test_osp_grillage.py b/tests/test_osp_grillage.py index 94402f11..5c35c0d2 100644 --- a/tests/test_osp_grillage.py +++ b/tests/test_osp_grillage.py @@ -370,10 +370,54 @@ def test_member_reassignment_feature(ref_bridge_properties): variant_one_model.create_osp_model(pyfile=False) # og.opsv.plot_model(element_labels=0, az_el=(-90, 0)) # plotting using ops_vis # og.plt.show() - assert ( - variant_one_model.element_command_list[2] - == 'ops.element("elasticBeamColumn", 2, *[2, 3], *[9.963e-02, 3.480e+10, 1.450e+10, 5.850e-04, 2.475e-04, 5.445e-04], 1, 106.272)\n' - ) + import re + + def parse_elastic_beam_cmd(cmd_str): + # Remove trailing newline + cmd_str = cmd_str.strip() + pattern = ( + r'ops\.element\(\s*"elasticBeamColumn"\s*,\s*(\d+)\s*,\s*\*\[(.*?)\]\s*,' + r'\s*\*\[(.*?)\]\s*,\s*(\d+)\s*,\s*"[-]?mass"\s*,\s*([\deE.+-]+)\s*\)' + ) + + # Match the command pattern + match = re.match( + pattern, + cmd_str, + ) + assert match, f"Failed to parse command: {cmd_str}" + + elem_id = int(match.group(1)) + nodes = list(map(int, match.group(2).split(","))) + section_vals = list(map(float, match.group(3).split(","))) + transf_tag = int(match.group(4)) + mass = float(match.group(5)) + + return elem_id, nodes, section_vals, transf_tag, mass + + # Expected values + expected_elem_id = 2 + expected_nodes = [2, 3] + expected_section_vals = [ + 9.963e-02, + 3.480e10, + 1.450e10, + 5.850e-04, + 2.475e-04, + 5.445e-04, + ] + expected_transf_tag = 1 + expected_mass = 106272.0 + + cmd = variant_one_model.element_command_list[2] + elem_id, nodes, section_vals, transf_tag, mass = parse_elastic_beam_cmd(cmd) + + assert elem_id == expected_elem_id + assert nodes == expected_nodes + for actual, expected in zip(section_vals, expected_section_vals): + assert abs(actual - expected) < 1e-9 + assert transf_tag == expected_transf_tag + assert abs(mass - expected_mass) < 1e-6 def test_create_offset_support(ref_bridge_properties):