From 32bcdc6543860eac70a8298960c479486bad1935 Mon Sep 17 00:00:00 2001 From: Ernesto Camarena Date: Thu, 19 Oct 2023 09:17:41 -0600 Subject: [PATCH 1/9] fixing ansys example --- src/pynumad/analysis/ansys/run.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pynumad/analysis/ansys/run.py b/src/pynumad/analysis/ansys/run.py index 43396eb..2b7faa1 100644 --- a/src/pynumad/analysis/ansys/run.py +++ b/src/pynumad/analysis/ansys/run.py @@ -5,7 +5,7 @@ def call_ansys(script_name,log=False,script_out='ansys.out',ncpus=1): for f in glob.glob("*.lock"): os.remove(f) - ansys_path=SOFTWARE_PATHS['ansys_path'] + ansys_path=SOFTWARE_PATHS['ansys'] MAXnLicenceTries=100 try: From 628fb7fe7167701b92334801cbccabb2757e36ad Mon Sep 17 00:00:00 2001 From: Ernesto Camarena Date: Thu, 19 Oct 2023 10:30:22 -0600 Subject: [PATCH 2/9] making sure ansys and cubit example work. adding to gitignore. adding to software_paths --- .gitignore | 19 +++++++++++++++++++ docs/user-guide/meshing.rst | 2 +- examples/cubit_beam.py | 19 +++++++++++++++---- examples/cubit_solid.py | 4 ++-- src/pynumad/analysis/make_models.py | 12 ++++++------ src/pynumad/software_paths.json | 3 ++- 6 files changed, 45 insertions(+), 14 deletions(-) diff --git a/.gitignore b/.gitignore index 82fa191..034f892 100644 --- a/.gitignore +++ b/.gitignore @@ -183,3 +183,22 @@ cython_debug/ # pynumad src/pynumad/software_paths.json + +#VABS +*.in +*.in.* + +#Cubit +*.cub +*.g +*.jou + +#BeamDyn +*BeamDyn*dat + +#Sierra +sd.i +sm.i +euler +mat_ori.py + diff --git a/docs/user-guide/meshing.rst b/docs/user-guide/meshing.rst index c42d9eb..b387f62 100644 --- a/docs/user-guide/meshing.rst +++ b/docs/user-guide/meshing.rst @@ -76,7 +76,7 @@ Then add Cubit to the path import sys sys.path.append(pynumad.SOFTWARE_PATHS['cubit']) - sys.path.append(pynumad.SOFTWARE_PATHS['cubitEnhancements']) + sys.path.append(pynumad.SOFTWARE_PATHS['cubit_enhancements']) import cubit diff --git a/examples/cubit_beam.py b/examples/cubit_beam.py index bcd27c6..cd24876 100644 --- a/examples/cubit_beam.py +++ b/examples/cubit_beam.py @@ -2,14 +2,14 @@ import pynumad sys.path.append(pynumad.SOFTWARE_PATHS['cubit']) -sys.path.append(pynumad.SOFTWARE_PATHS['cubitEnhancements']) +sys.path.append(pynumad.SOFTWARE_PATHS['cubit_enhancements']) import cubit from pynumad.analysis.cubit.make_blade import * import numpy as np from pynumad.analysis.make_models import write_beam_model - +import logging def get_cs_params(): @@ -75,12 +75,23 @@ def get_cs_params(): cs_params=get_cs_params() settings={} -settings['make_input_for']='VABS' #SM, VABS, ANBA, or None +settings['make_input_for']='VABSbeamdyn' #SM, VABS, ANBA, or None settings['export']='cubg' #cub, g, or None cubit_make_cross_sections(blade,wt_name,settings,cs_params,'2D',stationList=[],directory=dirName) #Note an empty list for stationList will make all cross sections. - +#Proportional damping values for BeamDyn file. mu=[0.00257593, 0.0017469, 0.0017469, 0.0017469, 0.00257593, 0.0017469] + + +#Set up logging +log =logging.getLogger(__name__) +log.setLevel(logging.DEBUG) +fh=logging.FileHandler(wt_name+'driver.log',mode='w') +log.addHandler(fh) + +#Read in a fresh new blade +blade=pynumad.Blade() +blade.read_yaml('example_data/'+yamlName+'.yaml') file_names=write_beam_model(wt_name,settings,blade,mu,log,directory=dirName) \ No newline at end of file diff --git a/examples/cubit_solid.py b/examples/cubit_solid.py index e387a3c..ff2d482 100644 --- a/examples/cubit_solid.py +++ b/examples/cubit_solid.py @@ -2,7 +2,7 @@ import pynumad sys.path.append(pynumad.SOFTWARE_PATHS['cubit']) -sys.path.append(pynumad.SOFTWARE_PATHS['cubitEnhancements']) +sys.path.append(pynumad.SOFTWARE_PATHS['cubit_enhancements']) import cubit from pynumad.analysis.cubit.make_blade import * @@ -75,7 +75,7 @@ def get_cs_params(): cs_params=get_cs_params() settings={} -settings['make_input_for']='Sd' #SM, VABS, ANBA, or None +settings['make_input_for']='SmSd' #SM, VABS, ANBA, or None settings['export']='cubg' #cub, g, or None materials_used=cubit_make_solid_blade(blade, wt_name, settings, cs_params, stationList=[2,3]) diff --git a/src/pynumad/analysis/make_models.py b/src/pynumad/analysis/make_models.py index 1b27ccb..b250a78 100644 --- a/src/pynumad/analysis/make_models.py +++ b/src/pynumad/analysis/make_models.py @@ -4,7 +4,7 @@ import glob import numpy as np from pynumad.utils.misc_utils import copy_and_replace - +from pynumad.paths import SOFTWARE_PATHS def write_beam_model(wt_name,settings,blade,mu,log,directory='.'): import pynumad.analysis.beam_utils as beam_utils @@ -31,7 +31,7 @@ def write_beam_model(wt_name,settings,blade,mu,log,directory='.'): for filePath in glob.glob(directory+'/'+wt_name+'*.in'): fileCount+=1 try: - this_cmd = 'VABS ' +filePath + this_cmd = SOFTWARE_PATHS['vabs']+' ' +filePath log.info(f' running: {this_cmd}') licenseAvailable=False @@ -100,7 +100,7 @@ def write_sierra_model(wt_name,settings,blade,materials_used,directory='.'): # #Runs VABS or OpenSG to homogenize # #Makes beamDyn or GEBT files - + template_path=SOFTWARE_PATHS['pynumad']+'src/data/templates/' # radial_stations=blade.ispan/blade.ispan[-1] # nStations=len(radial_stations) @@ -109,7 +109,7 @@ def write_sierra_model(wt_name,settings,blade,materials_used,directory='.'): materials = blade.definition.materials solver_string=settings['make_input_for'].lower() if 'sm' in solver_string or 'sd' in solver_string: - templateFileName='mat_ori.py.template' + templateFileName=template_path+'mat_ori.py.template' mat_ori_file_name='mat_ori.py' copy_and_replace(templateFileName, mat_ori_file_name, @@ -118,7 +118,7 @@ def write_sierra_model(wt_name,settings,blade,materials_used,directory='.'): }) if 'sm' in settings['make_input_for'].lower(): - templateFileName='sm.i.template' + templateFileName=template_path+'sm.i.template' adagioFileName='sm.i' materialLines=f'' @@ -167,7 +167,7 @@ def write_sierra_model(wt_name,settings,blade,materials_used,directory='.'): 'BLADE_BLOCKS': blockLines, }) if 'sd' in settings['make_input_for'].lower(): - templateFileName='sd.i.template' + templateFileName=template_path+'sd.i.template' adagioFileName='sd.i' materialLines=f'' diff --git a/src/pynumad/software_paths.json b/src/pynumad/software_paths.json index 755693e..8a18e88 100644 --- a/src/pynumad/software_paths.json +++ b/src/pynumad/software_paths.json @@ -1,6 +1,7 @@ { - "numad": "", + "pynumad": "", "ansys": "", + "vabs": "", "cubit": "", "cubit_enhancements": "" } From e920839cc75bde869d3a464a934cae507aa13a5d Mon Sep 17 00:00:00 2001 From: Camarena Date: Mon, 17 Feb 2025 17:27:40 -0500 Subject: [PATCH 3/9] adding layup angle read from yaml --- src/pynumad/objects/bom.py | 7 ++++++- src/pynumad/objects/stackdb.py | 4 ++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/pynumad/objects/bom.py b/src/pynumad/objects/bom.py index ba9a856..fa321c7 100644 --- a/src/pynumad/objects/bom.py +++ b/src/pynumad/objects/bom.py @@ -96,6 +96,8 @@ def generate(self, definition: Definition, keypoints: KeyPoints): cur_bom.area = regionarea cur_bom.thickness = mat.layerthickness cur_bom.weight = mat.drydensity * regionarea + cur_bom.angle = comp.fabricangle + self.indices["hp"].append( [begin_station[ks], end_station[ks], *hpRegion] ) @@ -127,6 +129,7 @@ def generate(self, definition: Definition, keypoints: KeyPoints): cur_bom.area = regionarea cur_bom.thickness = mat.layerthickness cur_bom.weight = mat.drydensity * regionarea + cur_bom.angle = comp.fabricangle self.indices["lp"].append( [begin_station[ks], end_station[ks], *lpRegion] ) @@ -187,6 +190,7 @@ def sorter(e): cur_bom.area = regionarea cur_bom.thickness = mat.layerthickness cur_bom.weight = mat.drydensity * regionarea + cur_bom.angle = comp.fabricangle self["sw"][swnum].append(cur_bom) self.indices["sw"][swnum].append( [begin_station[ks], end_station[ks]] @@ -271,4 +275,5 @@ def __init__(self): self.avgwidth: float = None self.area: float = None self.thickness: float = None - self.weight: float = None \ No newline at end of file + self.weight: float = None + self.angle: float = None \ No newline at end of file diff --git a/src/pynumad/objects/stackdb.py b/src/pynumad/objects/stackdb.py index a2a9f3d..da84ed9 100644 --- a/src/pynumad/objects/stackdb.py +++ b/src/pynumad/objects/stackdb.py @@ -72,7 +72,7 @@ def generate(self, keypoints: KeyPoints, bom: BillOfMaterials): cur_ply.component = bom["hp"][k].name cur_ply.materialid = bom["hp"][k].materialid cur_ply.thickness = bom["hp"][k].thickness - cur_ply.angle = 0 # TODO, set to 0 for now, bom['lp'](k, ); + cur_ply.angle = bom["hp"][k].angle # TODO, set to 0 for now, bom['lp'](k, ); cur_ply.nPlies = 1 # default to 1, modified in addply() if necessary # ... and add the ply to every area that is part of the region @@ -88,7 +88,7 @@ def generate(self, keypoints: KeyPoints, bom: BillOfMaterials): cur_ply.component = bom["lp"][k].name cur_ply.materialid = bom["lp"][k].materialid cur_ply.thickness = bom["lp"][k].thickness - cur_ply.angle = 0 # TODO, set to 0 for now, bom['lp'](k, ); + cur_ply.angle = bom["lp"][k].angle # TODO, set to 0 for now, bom['lp'](k, ); cur_ply.nPlies = 1 # default to 1, modified in addply() if necessary # ... and add the ply to every area that is part of the region From 4e471c720151c2e473715ad25de699de98f27360 Mon Sep 17 00:00:00 2001 From: Camarena Date: Thu, 6 Mar 2025 09:05:18 -0500 Subject: [PATCH 4/9] re-adding tet mesh capability to mat_orientations. Removed parallel processing from material orientations due to multi argment functions needed for mixed hex or tet capability. --- src/pynumad/analysis/cubit/make_blade.py | 117 +++++++++-------------- 1 file changed, 45 insertions(+), 72 deletions(-) diff --git a/src/pynumad/analysis/cubit/make_blade.py b/src/pynumad/analysis/cubit/make_blade.py index b90a426..abff448 100644 --- a/src/pynumad/analysis/cubit/make_blade.py +++ b/src/pynumad/analysis/cubit/make_blade.py @@ -110,7 +110,7 @@ def write_path_node_angles_to_file(set_verts,prepend,directory='.'): file.write(f'R_32 {" ".join(map(str,dcms[7]))}\n') file.write(f'R_33 {" ".join(map(str,dcms[8]))}\n') file.close() -def get_hex_orientations_euler(volume_id): +def get_orientations_euler(volume_id,element_shape_string): global_el_ids_in_vol=[] theta1s_in_vol=[] theta2s_in_vol=[] @@ -122,8 +122,14 @@ def get_hex_orientations_euler(volume_id): #volume_name = cubit.get_entity_name("volume", volume_id) #t0 = time.time() - for el_id in get_volume_hexes(volume_id): - coords = cubit.get_center_point("hex", el_id) + if 'hex' in element_shape_string: + element_ids = get_volume_hexes(volume_id) + elif 'tet' in element_shape_string: + element_ids = get_volume_tets(volume_id) + + + for el_id in element_ids: + coords = cubit.get_center_point(element_shape_string, el_id) surface_normal = vectNorm( list(sign*np.array(get_surface_normal_at_coord(surf_id_for_mat_ori, coords)))) @@ -138,7 +144,7 @@ def get_hex_orientations_euler(volume_id): globalAxisBasisVectors = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] - global_id=get_global_element_id('hex',el_id) + global_id=get_global_element_id(element_shape_string,el_id) dcm = getDCM(globalAxisBasisVectors, newCoordinateSystemVectors) @@ -152,7 +158,7 @@ def get_hex_orientations_euler(volume_id): return global_el_ids_in_vol,theta1s_in_vol,theta2s_in_vol,theta3s_in_vol -def get_hex_orientations_two_points(volume_id): +def get_orientations_two_points(volume_id,element_shape_string): global_el_ids_in_vol=[] @@ -164,8 +170,14 @@ def get_hex_orientations_two_points(volume_id): #volume_name = cubit.get_entity_name("volume", volume_id) #t0 = time.time() - for el_id in get_volume_hexes(volume_id): - coords = cubit.get_center_point("hex", el_id) + if 'hex' in element_shape_string: + element_ids = get_volume_hexes(volume_id) + elif 'tet' in element_shape_string: + element_ids = get_volume_tets(volume_id) + + + for el_id in element_ids: + coords = cubit.get_center_point(element_shape_string, el_id) surface_normal = vectNorm( list(sign*np.array(get_surface_normal_at_coord(surf_id_for_mat_ori, coords)))) @@ -176,7 +188,7 @@ def get_hex_orientations_two_points(volume_id): perimeter_direction = vectNorm(np.cross(surface_normal, spanwise_direction)) - global_id=get_global_element_id('hex',el_id) + global_id=get_global_element_id(element_shape_string,el_id) global_el_ids_in_vol.append(global_id) @@ -185,7 +197,7 @@ def get_hex_orientations_two_points(volume_id): return global_el_ids_in_vol,spanwise_directions_in_vol,perimeter_directions_in_vol -def get_hex_orientations_vectors(volume_id): +def get_orientations_vectors(volume_id,element_shape_string): global_el_ids_in_vol=[] @@ -197,9 +209,14 @@ def get_hex_orientations_vectors(volume_id): surf_id_for_mat_ori,sign = get_mat_ori_surface(volume_id) #volume_name = cubit.get_entity_name("volume", volume_id) #t0 = time.time() + if 'hex' in element_shape_string: + element_ids = get_volume_hexes(volume_id) + elif 'tet' in element_shape_string: + element_ids = get_volume_tets(volume_id) - for el_id in get_volume_hexes(volume_id): - coords = cubit.get_center_point("hex", el_id) + + for el_id in element_ids: + coords = cubit.get_center_point(element_shape_string, el_id) surface_normal = vectNorm( list(sign*np.array(get_surface_normal_at_coord(surf_id_for_mat_ori, coords)))) @@ -210,7 +227,7 @@ def get_hex_orientations_vectors(volume_id): perimeter_direction = vectNorm(np.cross(surface_normal, spanwise_direction)) - global_id=get_global_element_id('hex',el_id) + global_id=get_global_element_id(element_shape_string,el_id) global_el_ids_in_vol.append(global_id) @@ -220,52 +237,6 @@ def get_hex_orientations_vectors(volume_id): surface_normal_directions_in_vol.append(surface_normal) return global_el_ids_in_vol,spanwise_directions_in_vol,perimeter_directions_in_vol,surface_normal_directions_in_vol - -def get_tet_orientations(volume_id): - global_el_ids_in_vol=[] - theta1s_in_vol=[] - theta2s_in_vol=[] - theta3s_in_vol=[] - - - spanwise_directions_in_vol = [] - perimeter_directions_in_vol = [] - - surf_id_for_mat_ori,sign = get_mat_ori_surface(volume_id) - #volume_name = cubit.get_entity_name("volume", volume_id) - #t0 = time.time() - - for el_id in get_volume_tets(volume_id): - coords = cubit.get_center_point("tet", el_id) - - surface_normal = vectNorm( - list(sign*np.array(get_surface_normal_at_coord(surf_id_for_mat_ori, coords)))) - - ref_line_direction = [0,0,1] - #https://www.maplesoft.com/support/help/maple/view.aspx?path=MathApps%2FProjectionOfVectorOntoPlane - spanwise_direction = vectNorm(np.array(ref_line_direction)-np.dot(ref_line_direction,surface_normal)*np.array(surface_normal)) - - perimeter_direction = vectNorm(np.cross(surface_normal, spanwise_direction)) - - newCoordinateSystemVectors = [spanwise_direction,perimeter_direction,surface_normal] - - globalAxisBasisVectors = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] - - global_id=get_global_element_id('tet',el_id) - - dcm = getDCM(globalAxisBasisVectors, newCoordinateSystemVectors) - - temp1, temp2, temp3 = dcmToEulerAngles(dcm) - - global_el_ids_in_vol.append(global_id) - theta1s_in_vol.append(-1*temp1) - theta2s_in_vol.append(-1*temp2) - theta3s_in_vol.append(-1*temp3) - - spanwise_directions_in_vol.append(spanwise_direction) - perimeter_directions_in_vol.append(perimeter_direction) - - return global_el_ids_in_vol,theta1s_in_vol,theta2s_in_vol,theta3s_in_vol,spanwise_directions_in_vol,perimeter_directions_in_vol def assign_material_orientations(orientation_data,output_format = 'euler'): #Apply Material Orientation @@ -309,38 +280,40 @@ def compute_material_orientations(element_shape,output_format = 'euler',ncpus = # # ### Assign material orientations ### # # #################################### + element_shape_string = element_shape.lower() parse_string = f'with name "*volume*"' all_volume_ids = parse_cubit_list("volume", parse_string) t0 = time.time() print(f'Calculating material orientations with {ncpus} CPU(s)...') - if 'hex' in element_shape: + if 'hex' in element_shape_string or 'tet' in element_shape.lower(): if ncpus==1: ans = [] if 'euler' in output_format: for vol_id in all_volume_ids: - ans.append(get_hex_orientations_euler(vol_id)) + ans.append(get_orientations_euler(vol_id,element_shape_string)) elif 'two_points' in output_format: for vol_id in all_volume_ids: - ans.append(get_hex_orientations_two_points(vol_id)) + ans.append(get_orientations_two_points(vol_id,element_shape_string)) elif 'vectors' in output_format: for vol_id in all_volume_ids: - ans.append(get_hex_orientations_vectors(vol_id)) + ans.append(get_orientations_vectors(vol_id,element_shape_string)) else: raise NameError(f'Material Orientation output format: {output_format} is not supported') else: - pool_obj = multiprocessing.Pool(ncpus) - if 'euler' in output_format: - ans = pool_obj.map(get_hex_orientations_euler,all_volume_ids) - elif 'two_points' in output_format: - ans = pool_obj.map(get_hex_orientations_two_points,all_volume_ids) - elif 'vectors' in output_format: - ans = pool_obj.map(get_hex_orientations_vectors,all_volume_ids) - else: - raise NameError(f'Material Orientation output format: {output_format} is not supported') + raise ValueError(f'Number of CPUs has to be 1 for material orientations due to multiple args in get_orientations_XXXXXX() functions') + # pool_obj = multiprocessing.Pool(ncpus) + # if 'euler' in output_format: + # ans = pool_obj.map(get_hex_orientations_euler,all_volume_ids) + # elif 'two_points' in output_format: + # ans = pool_obj.map(get_hex_orientations_two_points,all_volume_ids) + # elif 'vectors' in output_format: + # ans = pool_obj.map(get_hex_orientations_vectors,all_volume_ids) + # else: + # raise NameError(f'Material Orientation output format: {output_format} is not supported') - pool_obj.close() + # pool_obj.close() else: raise NameError(f' element shape {element_shape} unsupported.') t1 = time.time() From d537b26c78719f2fea5d36940222e55be03bff20 Mon Sep 17 00:00:00 2001 From: Camarena Date: Wed, 12 Mar 2025 13:06:07 -0400 Subject: [PATCH 5/9] converting radians to degrees for use in pynumad --- src/pynumad/objects/bom.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pynumad/objects/bom.py b/src/pynumad/objects/bom.py index fa321c7..ff73fe9 100644 --- a/src/pynumad/objects/bom.py +++ b/src/pynumad/objects/bom.py @@ -96,7 +96,7 @@ def generate(self, definition: Definition, keypoints: KeyPoints): cur_bom.area = regionarea cur_bom.thickness = mat.layerthickness cur_bom.weight = mat.drydensity * regionarea - cur_bom.angle = comp.fabricangle + cur_bom.angle = comp.fabricangle*180/3.141592653589793 self.indices["hp"].append( [begin_station[ks], end_station[ks], *hpRegion] @@ -129,7 +129,7 @@ def generate(self, definition: Definition, keypoints: KeyPoints): cur_bom.area = regionarea cur_bom.thickness = mat.layerthickness cur_bom.weight = mat.drydensity * regionarea - cur_bom.angle = comp.fabricangle + cur_bom.angle = -1*comp.fabricangle*180/3.141592653589793 self.indices["lp"].append( [begin_station[ks], end_station[ks], *lpRegion] ) @@ -190,7 +190,7 @@ def sorter(e): cur_bom.area = regionarea cur_bom.thickness = mat.layerthickness cur_bom.weight = mat.drydensity * regionarea - cur_bom.angle = comp.fabricangle + cur_bom.angle = comp.fabricangle*180/3.141592653589793 self["sw"][swnum].append(cur_bom) self.indices["sw"][swnum].append( [begin_station[ks], end_station[ks]] From 2605b116e18d0a1e08e38e1305dec7e63108f792 Mon Sep 17 00:00:00 2001 From: Camarena Date: Thu, 13 Mar 2025 17:19:07 -0400 Subject: [PATCH 6/9] adding nonzero layup angle capability to cubit. unifying how material orietation vectors are computed. --- .../analysis/cubit/connect_cross_sections.py | 12 +- src/pynumad/analysis/cubit/make_blade.py | 179 ++++++++++-------- 2 files changed, 105 insertions(+), 86 deletions(-) diff --git a/src/pynumad/analysis/cubit/connect_cross_sections.py b/src/pynumad/analysis/cubit/connect_cross_sections.py index e6ee717..ac094b6 100644 --- a/src/pynumad/analysis/cubit/connect_cross_sections.py +++ b/src/pynumad/analysis/cubit/connect_cross_sections.py @@ -367,7 +367,7 @@ def get_spanwise_splines_for_a_volume(current_surface_id,next_surface_id,spanwis -def make_all_volumes_for_a_part(surface_dict, ordered_list, i_station_end,spanwise_splines): +def make_all_volumes_for_a_part(surface_dict, volume_dict,ordered_list, i_station_end,spanwise_splines): # nIntervals=3 vol_list=[] i_start = len(spanwise_splines) @@ -384,8 +384,14 @@ def make_all_volumes_for_a_part(surface_dict, ordered_list, i_station_end,spanwi spanwise_splines_for_a_volume = get_spanwise_splines_for_a_volume(current_surface_id,next_surface_id,spanwise_splines[part_surface_ids], surface_dict[current_surface_id]["verts"],surface_dict[next_surface_id]["verts"]) make_a_volume(i_span,current_surface_id,next_surface_id,spanwise_splines_for_a_volume,surface_dict,i_station_end) - vol_list.append(get_last_id("volume")) - # assign_intervals(get_last_id("volume"),nIntervals) + vol_id = get_last_id("volume") + vol_list.append(vol_id) + + volume_dict[vol_id] = {} + # volume_dict[vol_id]["material_name"] = material_name + volume_dict[vol_id]["ply_angle"] = surface_dict[current_surface_id]["ply_angle"] + del surface_dict[current_surface_id] + else: raise ValueError("Can't make volumes with only one cross section.") diff --git a/src/pynumad/analysis/cubit/make_blade.py b/src/pynumad/analysis/cubit/make_blade.py index abff448..e4ca3f7 100644 --- a/src/pynumad/analysis/cubit/make_blade.py +++ b/src/pynumad/analysis/cubit/make_blade.py @@ -197,26 +197,21 @@ def get_orientations_two_points(volume_id,element_shape_string): return global_el_ids_in_vol,spanwise_directions_in_vol,perimeter_directions_in_vol -def get_orientations_vectors(volume_id,element_shape_string): - global_el_ids_in_vol=[] - - - - spanwise_directions_in_vol = [] - perimeter_directions_in_vol = [] - surface_normal_directions_in_vol = [] - - surf_id_for_mat_ori,sign = get_mat_ori_surface(volume_id) - #volume_name = cubit.get_entity_name("volume", volume_id) - #t0 = time.time() - if 'hex' in element_shape_string: - element_ids = get_volume_hexes(volume_id) - elif 'tet' in element_shape_string: - element_ids = get_volume_tets(volume_id) +def get_orientations_vectors(element_ids,volume_dict): + + spanwise_directions = [] + perimeter_directions = [] + surface_normal_directions = [] for el_id in element_ids: - coords = cubit.get_center_point(element_shape_string, el_id) + volume_id = parse_cubit_list('volume', f'in element {el_id}')[0] + surf_id_for_mat_ori,sign = get_mat_ori_surface(volume_id) + + node_ids = parse_cubit_list('node',f'in element {el_id}') + coords = [list(get_nodal_coordinates(node_id)) for node_id in node_ids] + coords = np.array(coords) + coords = np.mean(coords, 0) surface_normal = vectNorm( list(sign*np.array(get_surface_normal_at_coord(surf_id_for_mat_ori, coords)))) @@ -226,17 +221,20 @@ def get_orientations_vectors(volume_id,element_shape_string): spanwise_direction = vectNorm(np.array(ref_line_direction)-np.dot(ref_line_direction,surface_normal)*np.array(surface_normal)) perimeter_direction = vectNorm(np.cross(surface_normal, spanwise_direction)) + + #Additional rotation about surface normal + theta = math.radians(volume_dict[volume_id]['ply_angle']) + c = math.cos(theta) + s = math.sin(theta) + beta = np.array([[c,s, 0],[-s, c, 0],[0,0,1]]) - global_id=get_global_element_id(element_shape_string,el_id) - - global_el_ids_in_vol.append(global_id) + new_directions = beta @ [spanwise_direction,perimeter_direction,surface_normal] - - spanwise_directions_in_vol.append(spanwise_direction) - perimeter_directions_in_vol.append(perimeter_direction) - surface_normal_directions_in_vol.append(surface_normal) + spanwise_directions.append(list(new_directions[0])) + perimeter_directions.append(list(new_directions[1])) + surface_normal_directions.append(list(new_directions[2])) - return global_el_ids_in_vol,spanwise_directions_in_vol,perimeter_directions_in_vol,surface_normal_directions_in_vol + return spanwise_directions,perimeter_directions,surface_normal_directions def assign_material_orientations(orientation_data,output_format = 'euler'): #Apply Material Orientation @@ -275,34 +273,48 @@ def assign_material_orientations(orientation_data,output_format = 'euler'): return -def compute_material_orientations(element_shape,output_format = 'euler',ncpus = 1): +def compute_material_orientations(element_shape,volume_dict, output_format = 'euler',ncpus = 1): # # #################################### # # ### Assign material orientations ### # # #################################### - element_shape_string = element_shape.lower() - parse_string = f'with name "*volume*"' - all_volume_ids = parse_cubit_list("volume", parse_string) + # element_shape_string = element_shape.lower() + parse_string = f'in volume with name "*volume*"' + global_element_ids = parse_cubit_list("element", parse_string) t0 = time.time() print(f'Calculating material orientations with {ncpus} CPU(s)...') - if 'hex' in element_shape_string or 'tet' in element_shape.lower(): - if ncpus==1: - ans = [] - if 'euler' in output_format: - for vol_id in all_volume_ids: - ans.append(get_orientations_euler(vol_id,element_shape_string)) - elif 'two_points' in output_format: - for vol_id in all_volume_ids: - ans.append(get_orientations_two_points(vol_id,element_shape_string)) - elif 'vectors' in output_format: - for vol_id in all_volume_ids: - ans.append(get_orientations_vectors(vol_id,element_shape_string)) - else: - raise NameError(f'Material Orientation output format: {output_format} is not supported') - else: - raise ValueError(f'Number of CPUs has to be 1 for material orientations due to multiple args in get_orientations_XXXXXX() functions') + ans = [] + if 'euler' in output_format: + for vol_id in all_volume_ids: + ans.append(get_orientations_euler(vol_id,element_shape_string)) + elif 'two_points' in output_format: + for vol_id in all_volume_ids: + ans.append(get_orientations_two_points(vol_id,element_shape_string)) + elif 'vectors' in output_format: + spanwise_directions,perimeter_directions,surface_normal_directions = get_orientations_vectors(global_element_ids,volume_dict) + # for global_element_id in global_element_ids: + # ans.append() + else: + raise NameError(f'Material Orientation output format: {output_format} is not supported') + + # if 'hex' in element_shape_string or 'tet' in element_shape.lower(): + # if ncpus==1: + # ans = [] + # if 'euler' in output_format: + # for vol_id in all_volume_ids: + # ans.append(get_orientations_euler(vol_id,element_shape_string)) + # elif 'two_points' in output_format: + # for vol_id in all_volume_ids: + # ans.append(get_orientations_two_points(vol_id,element_shape_string)) + # elif 'vectors' in output_format: + # for vol_id in all_volume_ids: + # ans.append(get_orientations_vectors(vol_id,element_shape_string)) + # else: + # raise NameError(f'Material Orientation output format: {output_format} is not supported') + # else: + # raise ValueError(f'Number of CPUs has to be 1 for material orientations due to multiple args in get_orientations_XXXXXX() functions') # pool_obj = multiprocessing.Pool(ncpus) # if 'euler' in output_format: # ans = pool_obj.map(get_hex_orientations_euler,all_volume_ids) @@ -314,50 +326,50 @@ def compute_material_orientations(element_shape,output_format = 'euler',ncpus = # raise NameError(f'Material Orientation output format: {output_format} is not supported') # pool_obj.close() - else: - raise NameError(f' element shape {element_shape} unsupported.') + # else: + # raise NameError(f' element shape {element_shape} unsupported.') t1 = time.time() print(f'Total time for material orientations: {t1-t0}') - ans=np.array(ans,dtype=object) - global_ids=[] + # ans=np.array(ans,dtype=object) + # global_ids=[] - if 'euler' in output_format: - theta1s=[] - theta2s=[] - theta3s=[] - for i in range(len(all_volume_ids)): - global_ids+=list(ans[i][0]) - theta1s+=list(ans[i][1]) - theta2s+=list(ans[i][2]) - theta3s+=list(ans[i][3]) - - return [global_ids,theta1s,theta2s,theta3s] + # if 'euler' in output_format: + # theta1s=[] + # theta2s=[] + # theta3s=[] + # for i in range(len(all_volume_ids)): + # global_ids+=list(ans[i][0]) + # theta1s+=list(ans[i][1]) + # theta2s+=list(ans[i][2]) + # theta3s+=list(ans[i][3]) + + # return [global_ids,theta1s,theta2s,theta3s] - elif 'two_points' in output_format: - spanwise_directions = [] - perimiter_directions = [] + # elif 'two_points' in output_format: + # spanwise_directions = [] + # perimiter_directions = [] - for i in range(len(all_volume_ids)): - global_ids+=list(ans[i][0]) - spanwise_directions+=list(ans[i][1]) - perimiter_directions+=list(ans[i][2]) + # for i in range(len(all_volume_ids)): + # global_ids+=list(ans[i][0]) + # spanwise_directions+=list(ans[i][1]) + # perimiter_directions+=list(ans[i][2]) - return [global_ids,spanwise_directions,perimiter_directions] + # return [global_ids,spanwise_directions,perimiter_directions] - elif 'vectors' in output_format: - spanwise_directions = [] - perimiter_directions = [] - normal_directions = [] + # elif 'vectors' in output_format: + # spanwise_directions = [] + # perimiter_directions = [] + # normal_directions = [] - for i in range(len(all_volume_ids)): - global_ids+=list(ans[i][0]) - spanwise_directions+=list(ans[i][1]) - perimiter_directions+=list(ans[i][2]) - normal_directions+=list(ans[i][3]) + # for i in range(len(all_volume_ids)): + # global_ids+=list(ans[i][0]) + # spanwise_directions+=list(ans[i][1]) + # perimiter_directions+=list(ans[i][2]) + # normal_directions+=list(ans[i][3]) - return [global_ids,spanwise_directions,perimiter_directions,normal_directions] + return [global_element_ids,spanwise_directions,perimeter_directions,surface_normal_directions] def order_path_points(points, ind): @@ -1054,12 +1066,13 @@ def cubit_make_solid_blade( # Make volumes along the span. ### ### ### ### mesh_vol_list = [] + volume_dict ={} part_name = "shell" ordered_list = get_ordered_list(part_name) spanwise_splines=[] if len(ordered_list) > 0: - shell_vol_list,spanwise_splines = make_all_volumes_for_a_part(surface_dict, ordered_list, i_station_end,spanwise_splines) + shell_vol_list,spanwise_splines = make_all_volumes_for_a_part(surface_dict, volume_dict, ordered_list, i_station_end,spanwise_splines) else: shell_vol_list=[] @@ -1067,14 +1080,14 @@ def cubit_make_solid_blade( ordered_list = get_ordered_list(part_name) ordered_list_web = ordered_list.copy() if ordered_list and len(ordered_list[0]) > 1: - web_vol_list,spanwise_splines = make_all_volumes_for_a_part(surface_dict, ordered_list, i_station_end,spanwise_splines) + web_vol_list,spanwise_splines = make_all_volumes_for_a_part(surface_dict, volume_dict, ordered_list, i_station_end,spanwise_splines) else: web_vol_list=[] part_name = "roundTEadhesive" ordered_list = get_ordered_list(part_name) if ordered_list and len(ordered_list[0]) > 1: - roundTEadhesive_vol_list,spanwise_splines = make_all_volumes_for_a_part(surface_dict, ordered_list, i_station_end,spanwise_splines) + roundTEadhesive_vol_list,spanwise_splines = make_all_volumes_for_a_part(surface_dict, volume_dict, ordered_list, i_station_end,spanwise_splines) else: roundTEadhesive_vol_list=[] @@ -1083,7 +1096,7 @@ def cubit_make_solid_blade( ordered_list = get_ordered_list(part_name) if ordered_list and len(ordered_list[0]) > 1: - flatTEadhesive_vol_list,spanwise_splines = make_all_volumes_for_a_part(surface_dict, ordered_list, i_station_end,spanwise_splines) + flatTEadhesive_vol_list,spanwise_splines = make_all_volumes_for_a_part(surface_dict, volume_dict, ordered_list, i_station_end,spanwise_splines) else: flatTEadhesive_vol_list=[] @@ -1461,4 +1474,4 @@ def cubit_make_solid_blade( # cubit.cmd(f'save as "{wt_name}.cub" overwrite') - return materials_used + return materials_used, volume_dict From 7f77b473e16a0b97fad5422bae72e51c650b6f48 Mon Sep 17 00:00:00 2001 From: Camarena Date: Fri, 14 Mar 2025 12:28:48 -0400 Subject: [PATCH 7/9] updated how material orientations are assigned to genisis mesh with new way of computing orientations. --- .../analysis/cubit/connect_cross_sections.py | 88 ------- src/pynumad/analysis/cubit/make_blade.py | 237 ++++++------------ 2 files changed, 76 insertions(+), 249 deletions(-) diff --git a/src/pynumad/analysis/cubit/connect_cross_sections.py b/src/pynumad/analysis/cubit/connect_cross_sections.py index ac094b6..f5fb3a2 100644 --- a/src/pynumad/analysis/cubit/connect_cross_sections.py +++ b/src/pynumad/analysis/cubit/connect_cross_sections.py @@ -625,92 +625,4 @@ def get_mat_ori_surface(volume_id): else: sign = -1.0 return mat_ori_surface,sign -def old_get_mat_ori_surface(volume_id, spanwise_mat_ori_curve): - # This function is used when assigning material orientations - # This gets returns the surface within a volume that will be used to get surface normals. - # The sign +-1 is also returned since some of the surfaces are oriented the wrong way - - approximate_thickness_direction = get_approximate_thickness_direction_for_volume(volume_id) - - - surface_ids = [] - volumeSurfaces = cubit.volume(volume_id).surfaces() - for current_surface in volumeSurfaces: - - parse_string = f'with name "*thickness*" in surface {current_surface.id()}' - thickness_curve_ids = parse_cubit_list("curve", parse_string) - - if len(thickness_curve_ids)==0: - surface_ids.append(current_surface.id()) - - # # Create a list of surface IDs in the given volume - # surface_ids = [] - # volumeSurfaces = cubit.volume(volume_id).surfaces() - # for current_surface in volumeSurfaces: - # surface_ids.append(current_surface.id()) - - # # Eliminate surfaces that have two curves named thickness: - # surface_ct = 0 - # for current_surface in volumeSurfaces: - # curve_ct = ( - # 0 # Counts the number of curves in the surface with name 'layer_thickness' - # ) - # for current_curve in current_surface.curves(): - # curve_name = cubit.get_entity_name("curve", current_curve.id()) - # if "layer_thickness" in curve_name: - # curve_ct += 1 - - # if curve_ct >= 2: - # surface_ct += 1 - # surface_ids.remove(current_surface.id()) - - # surface_ids now has the list of surfaces w/o thickness curves - if len(surface_ids) == 3 or len(surface_ids) == 2 or len(surface_ids) == 1: - if len(surface_ids) == 2: - surface_name = cubit.get_entity_name("surface", surface_ids[0]) - if "topFace" in surface_name: - surface_id = surface_ids[0] - else: - surface_id = surface_ids[-1] - elif len(surface_ids) == 1: # Web overwrap - surface_id = surface_ids[0] - elif len(surface_ids) == 3: - if 'web' in cubit.get_entity_name("volume", volume_id).lower(): - #This is for when birdsmouth is made and it does - #not cut into the next station - max_area=0 - for i_surf in surface_ids: - area =cubit.surface(i_surf).area() - if area > max_area: - max_area=area - surface_id=i_surf - else: - raise ValueError( - "The number of thickness curves in volume is unexpected. Cannot assign material orientation" - ) - - coords = cubit.get_center_point("surface", surface_id) - surface_normal = cubit.surface(surface_id).normal_at(coords) - - if dotProd(surface_normal, approximate_thickness_direction) > 0: - sign = 1.0 - else: - sign = -1.0 - - - - - elif len(surface_ids) == 0: - - # LE adhesive and/or TE adhesive for round cross-sections - # print(f'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~volume_id {volume_id}') - surface_id = False - sign = 1.0 - #parse_string = f'with name "*layer_thickness*" in volume {volume_id}' - #thickness_curve_ids = parse_cubit_list("curve", parse_string) - else: - raise ValueError( - "The number of thickness curves in volume is unexpected. Cannot assign material orientation" - ) - return surface_id, sign diff --git a/src/pynumad/analysis/cubit/make_blade.py b/src/pynumad/analysis/cubit/make_blade.py index e4ca3f7..3825f97 100644 --- a/src/pynumad/analysis/cubit/make_blade.py +++ b/src/pynumad/analysis/cubit/make_blade.py @@ -80,9 +80,9 @@ def write_path_node_angles_to_file(set_verts,prepend,directory='.'): #https://www.maplesoft.com/support/help/maple/view.aspx?path=MathApps%2FProjectionOfVectorOntoPlane spanwise_direction = vectNorm(np.array(ref_line_direction)-np.dot(ref_line_direction,surface_normal)*np.array(surface_normal)) - perimeter_direction = vectNorm(np.cross(surface_normal, spanwise_direction)) + hoop_direction = vectNorm(np.cross(surface_normal, spanwise_direction)) - newCoordinateSystemVectors = [spanwise_direction,perimeter_direction,surface_normal] + newCoordinateSystemVectors = [spanwise_direction,hoop_direction,surface_normal] globalAxisBasisVectors = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] @@ -138,9 +138,9 @@ def get_orientations_euler(volume_id,element_shape_string): #https://www.maplesoft.com/support/help/maple/view.aspx?path=MathApps%2FProjectionOfVectorOntoPlane spanwise_direction = vectNorm(np.array(ref_line_direction)-np.dot(ref_line_direction,surface_normal)*np.array(surface_normal)) - perimeter_direction = vectNorm(np.cross(surface_normal, spanwise_direction)) + hoop_direction = vectNorm(np.cross(surface_normal, spanwise_direction)) - newCoordinateSystemVectors = [spanwise_direction,perimeter_direction,surface_normal] + newCoordinateSystemVectors = [spanwise_direction,hoop_direction,surface_normal] globalAxisBasisVectors = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] @@ -158,50 +158,11 @@ def get_orientations_euler(volume_id,element_shape_string): return global_el_ids_in_vol,theta1s_in_vol,theta2s_in_vol,theta3s_in_vol -def get_orientations_two_points(volume_id,element_shape_string): - global_el_ids_in_vol=[] - - - - spanwise_directions_in_vol = [] - perimeter_directions_in_vol = [] - - surf_id_for_mat_ori,sign = get_mat_ori_surface(volume_id) - #volume_name = cubit.get_entity_name("volume", volume_id) - #t0 = time.time() - - if 'hex' in element_shape_string: - element_ids = get_volume_hexes(volume_id) - elif 'tet' in element_shape_string: - element_ids = get_volume_tets(volume_id) - - - for el_id in element_ids: - coords = cubit.get_center_point(element_shape_string, el_id) - - surface_normal = vectNorm( - list(sign*np.array(get_surface_normal_at_coord(surf_id_for_mat_ori, coords)))) - - ref_line_direction = [0,0,1] - #https://www.maplesoft.com/support/help/maple/view.aspx?path=MathApps%2FProjectionOfVectorOntoPlane - spanwise_direction = vectNorm(np.array(ref_line_direction)-np.dot(ref_line_direction,surface_normal)*np.array(surface_normal)) - - perimeter_direction = vectNorm(np.cross(surface_normal, spanwise_direction)) - - global_id=get_global_element_id(element_shape_string,el_id) - - global_el_ids_in_vol.append(global_id) - - spanwise_directions_in_vol.append(spanwise_direction) - perimeter_directions_in_vol.append(perimeter_direction) - - return global_el_ids_in_vol,spanwise_directions_in_vol,perimeter_directions_in_vol - -def get_orientations_vectors(element_ids,volume_dict): +def get_element_orientations_vectors(element_ids,volume_dict): spanwise_directions = [] - perimeter_directions = [] + hoop_directions = [] surface_normal_directions = [] for el_id in element_ids: @@ -220,7 +181,7 @@ def get_orientations_vectors(element_ids,volume_dict): #https://www.maplesoft.com/support/help/maple/view.aspx?path=MathApps%2FProjectionOfVectorOntoPlane spanwise_direction = vectNorm(np.array(ref_line_direction)-np.dot(ref_line_direction,surface_normal)*np.array(surface_normal)) - perimeter_direction = vectNorm(np.cross(surface_normal, spanwise_direction)) + hoop_direction = vectNorm(np.cross(surface_normal, spanwise_direction)) #Additional rotation about surface normal theta = math.radians(volume_dict[volume_id]['ply_angle']) @@ -228,148 +189,102 @@ def get_orientations_vectors(element_ids,volume_dict): s = math.sin(theta) beta = np.array([[c,s, 0],[-s, c, 0],[0,0,1]]) - new_directions = beta @ [spanwise_direction,perimeter_direction,surface_normal] + new_directions = beta @ [spanwise_direction,hoop_direction,surface_normal] spanwise_directions.append(list(new_directions[0])) - perimeter_directions.append(list(new_directions[1])) + hoop_directions.append(list(new_directions[1])) surface_normal_directions.append(list(new_directions[2])) - return spanwise_directions,perimeter_directions,surface_normal_directions - -def assign_material_orientations(orientation_data,output_format = 'euler'): - #Apply Material Orientation + return spanwise_directions,hoop_directions,surface_normal_directions +def assign_material_orientation_angles(orientation_data): + #Assigne Euler angles to exodus mesh variables global_ids=orientation_data[0] n_el = len(global_ids) - if output_format =='euler': - theta1s=orientation_data[1] - theta2s=orientation_data[2] - theta3s=orientation_data[3] - cubit.set_element_variable(global_ids, 'rotation_angle_one', theta1s) - cubit.set_element_variable(global_ids, 'rotation_angle_two', theta2s) - cubit.set_element_variable(global_ids, 'rotation_angle_three', theta3s) + theta1s=orientation_data[1] + theta2s=orientation_data[2] + theta3s=orientation_data[3] - cubit.set_element_variable(global_ids, 'rotation_axis_one', 1*np.ones(n_el)) - cubit.set_element_variable(global_ids, 'rotation_axis_two', 2*np.ones(n_el)) - cubit.set_element_variable(global_ids, 'rotation_axis_three', 3*np.ones(n_el)) - elif output_format =='vectors': + cubit.set_element_variable(global_ids, 'rotation_angle_one', theta1s) + cubit.set_element_variable(global_ids, 'rotation_angle_two', theta2s) + cubit.set_element_variable(global_ids, 'rotation_angle_three', theta3s) - one_axis=np.array(orientation_data[1]) - two_axis=np.array(orientation_data[2]) - three_axis=np.array(orientation_data[3]) + cubit.set_element_variable(global_ids, 'rotation_axis_one', 1*np.ones(n_el)) + cubit.set_element_variable(global_ids, 'rotation_axis_two', 2*np.ones(n_el)) + cubit.set_element_variable(global_ids, 'rotation_axis_three', 3*np.ones(n_el)) - cubit.set_element_variable(global_ids, 'matCoord_1_x', one_axis[:,0]) - cubit.set_element_variable(global_ids, 'matCoord_1_y', one_axis[:,1]) - cubit.set_element_variable(global_ids, 'matCoord_1_z', one_axis[:,2]) + return +def assign_material_orientation_vectors(orientation_data): + #Assign material orientation vectors to exodus mesh variables + global_ids=orientation_data[0] + n_el = len(global_ids) + + - cubit.set_element_variable(global_ids, 'matCoord_2_x', two_axis[:,0]) - cubit.set_element_variable(global_ids, 'matCoord_2_y', two_axis[:,1]) - cubit.set_element_variable(global_ids, 'matCoord_2_z', two_axis[:,2]) + one_axis=np.array(orientation_data[1]) + two_axis=np.array(orientation_data[2]) + three_axis=np.array(orientation_data[3]) - cubit.set_element_variable(global_ids, 'matCoord_3_x', three_axis[:,0]) - cubit.set_element_variable(global_ids, 'matCoord_3_y', three_axis[:,1]) - cubit.set_element_variable(global_ids, 'matCoord_3_z', three_axis[:,2]) + cubit.set_element_variable(global_ids, 'matCoord_1_x', one_axis[:,0]) + cubit.set_element_variable(global_ids, 'matCoord_1_y', one_axis[:,1]) + cubit.set_element_variable(global_ids, 'matCoord_1_z', one_axis[:,2]) + + cubit.set_element_variable(global_ids, 'matCoord_2_x', two_axis[:,0]) + cubit.set_element_variable(global_ids, 'matCoord_2_y', two_axis[:,1]) + cubit.set_element_variable(global_ids, 'matCoord_2_z', two_axis[:,2]) + + cubit.set_element_variable(global_ids, 'matCoord_3_x', three_axis[:,0]) + cubit.set_element_variable(global_ids, 'matCoord_3_y', three_axis[:,1]) + cubit.set_element_variable(global_ids, 'matCoord_3_z', three_axis[:,2]) return -def compute_material_orientations(element_shape,volume_dict, output_format = 'euler',ncpus = 1): - # # #################################### - # # ### Assign material orientations ### - # # #################################### - - # element_shape_string = element_shape.lower() - parse_string = f'in volume with name "*volume*"' - global_element_ids = parse_cubit_list("element", parse_string) - +def get_material_orientation_angles(orientation_vectors): + # orientation_vectors is the output of assign_material_orientation_vectors() + t0 = time.time() - print(f'Calculating material orientations with {ncpus} CPU(s)...') + print(f'Converting orientation vectors to Euler Angles...') - ans = [] - if 'euler' in output_format: - for vol_id in all_volume_ids: - ans.append(get_orientations_euler(vol_id,element_shape_string)) - elif 'two_points' in output_format: - for vol_id in all_volume_ids: - ans.append(get_orientations_two_points(vol_id,element_shape_string)) - elif 'vectors' in output_format: - spanwise_directions,perimeter_directions,surface_normal_directions = get_orientations_vectors(global_element_ids,volume_dict) - # for global_element_id in global_element_ids: - # ans.append() - else: - raise NameError(f'Material Orientation output format: {output_format} is not supported') - - # if 'hex' in element_shape_string or 'tet' in element_shape.lower(): - # if ncpus==1: - # ans = [] - # if 'euler' in output_format: - # for vol_id in all_volume_ids: - # ans.append(get_orientations_euler(vol_id,element_shape_string)) - # elif 'two_points' in output_format: - # for vol_id in all_volume_ids: - # ans.append(get_orientations_two_points(vol_id,element_shape_string)) - # elif 'vectors' in output_format: - # for vol_id in all_volume_ids: - # ans.append(get_orientations_vectors(vol_id,element_shape_string)) - # else: - # raise NameError(f'Material Orientation output format: {output_format} is not supported') - # else: - # raise ValueError(f'Number of CPUs has to be 1 for material orientations due to multiple args in get_orientations_XXXXXX() functions') - # pool_obj = multiprocessing.Pool(ncpus) - # if 'euler' in output_format: - # ans = pool_obj.map(get_hex_orientations_euler,all_volume_ids) - # elif 'two_points' in output_format: - # ans = pool_obj.map(get_hex_orientations_two_points,all_volume_ids) - # elif 'vectors' in output_format: - # ans = pool_obj.map(get_hex_orientations_vectors,all_volume_ids) - # else: - # raise NameError(f'Material Orientation output format: {output_format} is not supported') - - # pool_obj.close() - # else: - # raise NameError(f' element shape {element_shape} unsupported.') - t1 = time.time() - print(f'Total time for material orientations: {t1-t0}') + globalAxisBasisVectors = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] + R_1_angles = [] + R_2_angles = [] + R_3_angles = [] + for i_el, element_id in enumerate(orientation_vectors[0]): + spanwise_direction = orientation_vectors[1][i_el] + hoop_direction = orientation_vectors[2][i_el] + surface_normal = orientation_vectors[3][i_el] - # ans=np.array(ans,dtype=object) - # global_ids=[] + newCoordinateSystemVectors = [spanwise_direction,hoop_direction,surface_normal] - # if 'euler' in output_format: - # theta1s=[] - # theta2s=[] - # theta3s=[] - # for i in range(len(all_volume_ids)): - # global_ids+=list(ans[i][0]) - # theta1s+=list(ans[i][1]) - # theta2s+=list(ans[i][2]) - # theta3s+=list(ans[i][3]) + dcm = getDCM(globalAxisBasisVectors, newCoordinateSystemVectors) - # return [global_ids,theta1s,theta2s,theta3s] - - # elif 'two_points' in output_format: - # spanwise_directions = [] - # perimiter_directions = [] + temp1, temp2, temp3 = dcmToEulerAngles(dcm) - # for i in range(len(all_volume_ids)): - # global_ids+=list(ans[i][0]) - # spanwise_directions+=list(ans[i][1]) - # perimiter_directions+=list(ans[i][2]) + R_1_angles.append(-1*temp1) + R_2_angles.append(-1*temp2) + R_3_angles.append(-1*temp3) - # return [global_ids,spanwise_directions,perimiter_directions] + t1 = time.time() + print(f'Total time for Euler angles: {t1-t0}') + return orientation_vectors[0],R_1_angles,R_2_angles,R_3_angles - # elif 'vectors' in output_format: - # spanwise_directions = [] - # perimiter_directions = [] - # normal_directions = [] +def get_material_orientation_vectors(volume_dict,ncpus = 1): + # # #################################### + # # ### Get material orientations ### + # # #################################### - # for i in range(len(all_volume_ids)): - # global_ids+=list(ans[i][0]) - # spanwise_directions+=list(ans[i][1]) - # perimiter_directions+=list(ans[i][2]) - # normal_directions+=list(ans[i][3]) + parse_string = f'in volume with name "*volume*"' + global_element_ids = parse_cubit_list("element", parse_string) + + t0 = time.time() + print(f'Calculating material orientations with {ncpus} CPU(s)...') + spanwise_directions,hoop_directions,surface_normal_directions = get_element_orientations_vectors(global_element_ids,volume_dict) + t1 = time.time() + print(f'Total time for material orientations: {t1-t0}') - return [global_element_ids,spanwise_directions,perimeter_directions,surface_normal_directions] + return [global_element_ids,spanwise_directions,hoop_directions,surface_normal_directions] def order_path_points(points, ind): From 6f883e338a07202fac3fd358387c09eb99de36a6 Mon Sep 17 00:00:00 2001 From: Camarena Date: Fri, 14 Mar 2025 14:20:04 -0400 Subject: [PATCH 8/9] updating sierra examples to work with new material orientations. updating sierra template files. adding new function yaml_mesh_to_cubit in make_blade.py --- examples/cubit_solid.py | 10 +-- .../analysis/cubit/connect_cross_sections.py | 55 +--------------- src/pynumad/analysis/cubit/make_blade.py | 62 +++++++++++++++++++ src/pynumad/data/templates/sd.i.template | 17 +++-- src/pynumad/data/templates/sm.i.template | 10 +-- 5 files changed, 86 insertions(+), 68 deletions(-) diff --git a/examples/cubit_solid.py b/examples/cubit_solid.py index b7a2437..26c08be 100644 --- a/examples/cubit_solid.py +++ b/examples/cubit_solid.py @@ -80,13 +80,15 @@ def get_cs_params(): #Make Cubit Geometry station_list = [2,3] -materials_used=cubit_make_solid_blade(blade, wt_name, settings, cs_params, stationList=station_list) +materials_used, volume_dict=cubit_make_solid_blade(blade, wt_name, settings, cs_params, stationList=station_list) #Compute material orientation -orientation_data=compute_material_orientations(cs_params['element_shape'],ncpus = 1) +orientation_vectors=get_material_orientation_vectors(volume_dict,ncpus = 1) +orientation_angles=get_material_orientation_angles(orientation_vectors) #assign material orientation in Cubit -assign_material_orientations(orientation_data) +assign_material_orientation_vectors(orientation_vectors) +assign_material_orientation_angles(orientation_angles) #Export mesh in Genisis format cubit.cmd(f'export mesh "{wt_name}.g" overwrite') @@ -94,7 +96,7 @@ def get_cs_params(): #Write Sierra input file from pynumad.paths import SOFTWARE_PATHS -template_path=SOFTWARE_PATHS['pynumad']+'src/data/templates/' +template_path=SOFTWARE_PATHS['pynumad']+'src/pynumad/data/templates/' write_sierra_sm_model(template_path+'sm.i.template',wt_name,station_list,blade,materials_used,'.') diff --git a/src/pynumad/analysis/cubit/connect_cross_sections.py b/src/pynumad/analysis/cubit/connect_cross_sections.py index f5fb3a2..e78fa84 100644 --- a/src/pynumad/analysis/cubit/connect_cross_sections.py +++ b/src/pynumad/analysis/cubit/connect_cross_sections.py @@ -535,60 +535,7 @@ def make_birds_mouth( return list(parse_cubit_list("volume", parse_string)) #update the web volumes -# cubit.cmd('open "/home/ecamare/myprojects/bar/cubitDev/python/python0.cub"') - - -# def get_approximate_thickness_direction_for_volume(volume_id): -# # This function is used when assigning material orientations - -# # Get thickness direction tangents - -# #Get list of curves with name layer_thickness -# parse_string = f'with name "*thickness*" in volume {volume_id}' -# thickness_curve_ids = parse_cubit_list("curve", parse_string) - -# approximate_thickness_direction = [] -# for i_curve in thickness_curve_ids: -# current_curve=cubit.curve(i_curve) -# coords = current_curve.position_from_fraction(0.5) -# approximate_thickness_direction.append(current_curve.tangent(coords)) -# approximate_thickness_direction = np.array(approximate_thickness_direction) -# n_thickness_curves=len(thickness_curve_ids) -# # approximate_thickness_direction = [] -# # for current_curve in cubit.volume(volume_id).curves(): -# # curve_name = cubit.get_entity_name("curve", current_curve.id()) -# # if "layer_thickness" in curve_name: -# # coords = current_curve.position_from_fraction(0.5) -# # approximate_thickness_direction.append(current_curve.tangent(coords)) -# # approximate_thickness_direction = np.array(approximate_thickness_direction) -# # n_thickness_curves, _ = approximate_thickness_direction.shape - -# if n_thickness_curves == 4: # All other cases -# return np.mean(approximate_thickness_direction, 0) -# elif n_thickness_curves == 8: # LE adhesive case and round TE adhesive -# return 0 -# elif n_thickness_curves == 6: # Web overwrap -# # Take the mean of all curves with name 'layer_thickness' -# mean = np.mean(approximate_thickness_direction, 0) - -# errorList = [] -# for i in range(n_thickness_curves): -# diff = approximate_thickness_direction[i] - mean - -# errorList.append(sqrt(dotProd(diff, diff))) -# sortIndex = np.argsort(errorList)[:4] # Take the first four. -# # This discards the two directions -# # with the largest deviation from the -# # average - -# return np.mean(approximate_thickness_direction[sortIndex, :], 0) -# else: -# cubit.cmd(f'save as "Debug.cub" overwrite') -# raise ValueError( -# f"The number of thickness curves in volume is unexpected. Cannot assign material orientation. n_thickness_curves: {n_thickness_curves}" -# ) - -# return + def get_mat_ori_surface(volume_id): diff --git a/src/pynumad/analysis/cubit/make_blade.py b/src/pynumad/analysis/cubit/make_blade.py index 3825f97..4ec5abd 100644 --- a/src/pynumad/analysis/cubit/make_blade.py +++ b/src/pynumad/analysis/cubit/make_blade.py @@ -1390,3 +1390,65 @@ def cubit_make_solid_blade( return materials_used, volume_dict + +def yaml_mesh_to_cubit(yaml_file_base,element_type,plot_mat_ori = True): + import yaml + + print(f'Importing {yaml_file_base} to cubit ...') + with open(f'{yaml_file_base}.yaml', 'r') as file: + mesh_data = yaml.load(file, Loader=yaml.CLoader) + + print(' Making nodes ...') + for node_coords in mesh_data['nodes']: + cubit.silent_cmd(f'create node location {node_coords[0]}') + + print(' Making elements ...') + for element_conn in mesh_data['elements']: + + cubit.silent_cmd(f'create {element_type} node {element_conn[0]}') + + print(' Making blocks ...') + for i_mat, mat_data in enumerate(mesh_data['materials']): + mat_name = mat_data['name'] + for el_set in mesh_data['sets']['element']: + if mat_name.lower() == el_set['name'].lower(): + cubit.silent_cmd(f'block {i_mat+1} add {element_type} {str(el_set["labels"]).replace("[", "").replace("]", "")}') + cubit.silent_cmd(f'block {i_mat+1} name "{mat_name}"') + + + if plot_mat_ori: + print(' Material orientation lines ...') + dir_strings = ['xdir','ydir','zdir',] + for i_el, element_ori in enumerate(mesh_data['elementOrientations']): + hex_id = i_el+1 + + node_ids = cubit.get_expanded_connectivity(element_type, hex_id) + coords = [] + for iNd, node_id in enumerate(node_ids): + coords.append(list(get_nodal_coordinates(node_id))) + coords = np.array(coords) + + # #######For Plotting - find the largest element side length ####### + distances=[] + for iNd,node_id in enumerate(node_ids): + for jNd,node_idj in enumerate(node_ids): + distances.append(norm(vectSub(coords[iNd],coords[jNd]))) + length=max(distances) + + coords = np.mean(coords, 0) + + cubit.create_vertex(coords[0],coords[1],coords[2]) + iVert1=get_last_id("vertex") + for i_dir in range(3): + index = 3*i_dir + cubit.create_vertex(coords[0]+length*element_ori[index],coords[1]+length*element_ori[index+1],coords[2]+length*element_ori[index+2]) + iVert2=get_last_id("vertex") + cubit.silent_cmd(f'create curve vertex {iVert1} {iVert2}') + cubit.silent_cmd(f'curve {get_last_id("curve")} name "{dir_strings[i_dir]}"') + + + cubit.silent_cmd(f'color curve with name "xdir*" geometry seagreen') + cubit.silent_cmd(f'color curve with name "ydir*" geometry blue') + cubit.silent_cmd(f'color curve with name "zdir*" geometry red') + + print(f'Done importing! \n') diff --git a/src/pynumad/data/templates/sd.i.template b/src/pynumad/data/templates/sd.i.template index a838418..6f44084 100644 --- a/src/pynumad/data/templates/sd.i.template +++ b/src/pynumad/data/templates/sd.i.template @@ -7,13 +7,20 @@ SOLUTION END file - geometry_file 'IN_MESH' + geometry_file 'IN_MESH' + initialize variable name = material_direction_1 + read variable = matCoord_1_ + variable type = element + initialize variable name = material_direction_2 + read variable = matCoord_2_ + variable type = element + initialize variable name = material_direction_3 + read variable = matCoord_3_ + variable type = element end BLADE_BLOCKS - - outputs disp stress @@ -26,8 +33,6 @@ end BLADE_MATERIALS boundary - nodeset root + nodeset ROOT_STATION_ns fixed - nodeset tip - z=-10.5 end diff --git a/src/pynumad/data/templates/sm.i.template b/src/pynumad/data/templates/sm.i.template index 1a5f887..32dbed1 100644 --- a/src/pynumad/data/templates/sm.i.template +++ b/src/pynumad/data/templates/sm.i.template @@ -180,6 +180,8 @@ BLADE_BLOCKS function = ramp scale factor = -750 end traction + + ### ------------------### ### Solver definition ### ### ------------------### @@ -188,11 +190,11 @@ BLADE_BLOCKS type = secant end begin cg - target relative residual = 2.0e-5 - acceptable relative residual = 4.0e-5 + target relative residual = 2.0e-3 + acceptable relative residual = 4.0e-3 target residual = 1.0e-7 acceptable residual = 2.0e-7 - maximum iterations = 3000 + maximum iterations = 4000 iteration print = 100 reference = belytschko preconditioner = probe @@ -208,7 +210,7 @@ BLADE_BLOCKS cutback factor = 0.5 growth factor = 1.0 maximum failure cutbacks = 10 - target iterations = 3000 + target iterations = 4000 iteration window = 200 end adaptive time stepping end adagio region adagio From d8d91d4dab202c796e4e91ea5a71943938ad54b4 Mon Sep 17 00:00:00 2001 From: Camarena Date: Fri, 14 Mar 2025 14:41:26 -0400 Subject: [PATCH 9/9] updating documentation on how to use material orientations --- docs/user-guide/solid_models.rst | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/docs/user-guide/solid_models.rst b/docs/user-guide/solid_models.rst index 8b046d2..59bdc0e 100644 --- a/docs/user-guide/solid_models.rst +++ b/docs/user-guide/solid_models.rst @@ -7,23 +7,25 @@ Pure solid model (explicitly descritezed sandwich panels) ========================================================= Pure solid models can be made in Sierra SM, Sierra SD or Abaqus. -Currenly only SD is partially supported since it does not allow for -spatially varying material orientations. + Continuous meshes in Sierra ---------------------------- -#. The first step is to make a Genesis mesh file and the associated files with - :py:func:`~pynumad.analysis.cubit.make_blade.cubit_make_solid_blade`. This - will create the following file +#. The first step is to make a solid element Cubit model with + :py:func:`~pynumad.analysis.cubit.make_blade.cubit_make_solid_blade`. - * {wt_name}.g. Genesis mesh file. +#. The next step is to compute material orientation vectors for each element with + :py:func:`~pynumad.analysis.cubit.make_blade.get_material_orientation_vectors`. -#. The next step is to Compute material orientation with - :py:func:`~pynumad.analysis.cubit.make_blade.compute_material_orientations` +#. If Sierra SD is a target model, these orientation vectors can be assigned directly with + :py:func:`~pynumad.analysis.cubit.make_blade.assign_material_orientation_vectors`. -#. Next, the orientation data needs to be assigned with - :py:func:`~pynumad.analysis.cubit.make_blade.assign_material_orientations` +#. If Sierra SM is a target model, the orientation vectors need to be converted to Euler Angles with + :py:func:`~pynumad.analysis.cubit.make_blade.get_material_orientation_angles`. + +#. The Euler Angles are then assigned to the mesh with + :py:func:`~pynumad.analysis.cubit.make_blade.assign_material_orientation_angles`. #. Finally, the mesh needs to be exported in Genisis format .. code-block:: python @@ -62,7 +64,7 @@ Continuous meshes in Sierra #. AND/OR run Sierra SD with: launch -n 10 salinas -i sd.i, where n is the number of CPUs. -An example called `cubit_solid.py` exists in the examples folder. +An example called `cubit_solid.py` exists in the examples folder illustrates these steps. Discontinuous meshes in Abaqus ------------------------------