From cac70fb098c43dc935e613ee5e112d21696a8201 Mon Sep 17 00:00:00 2001 From: GBenedett Date: Mon, 20 Nov 2023 17:03:01 +0100 Subject: [PATCH] mesh size factor from geometry evaluation --- ceasiompy/CPACS2GMSH/__specs__.py | 63 ++++- ceasiompy/CPACS2GMSH/cpacs2gmsh.py | 13 +- ceasiompy/CPACS2GMSH/func/generategmesh.py | 21 +- ceasiompy/CPACS2GMSH/func/mesh_sizing.py | 308 +++++++++++++++++++++ ceasiompy/utils/commonxpath.py | 2 + 5 files changed, 387 insertions(+), 20 deletions(-) create mode 100644 ceasiompy/CPACS2GMSH/func/mesh_sizing.py diff --git a/ceasiompy/CPACS2GMSH/__specs__.py b/ceasiompy/CPACS2GMSH/__specs__.py index 18b06f3a4..a039ef4c5 100644 --- a/ceasiompy/CPACS2GMSH/__specs__.py +++ b/ceasiompy/CPACS2GMSH/__specs__.py @@ -18,6 +18,8 @@ GMSH_REFINE_FACTOR_XPATH, GMSH_SYMMETRY_XPATH, SU2MESH_XPATH, + GMSH_MESH_SIZE_FACTOR_FUSELAGE_XPATH, + GMSH_MESH_SIZE_FACTOR_WINGS_XPATH, ) from ceasiompy.utils.moduleinterfaces import CPACSInOut @@ -97,30 +99,55 @@ gui_group="Mesh size", ) +# cpacs_inout.add_input( +# var_name="fuselage_mesh_size", +# var_type=float, +# default_value=0.3, +# unit="[m]", +# descr="Value assigned for the fuselage surfaces mesh size", +# xpath=GMSH_MESH_SIZE_FUSELAGE_XPATH, +# gui=True, +# gui_name="Fuselage", +# gui_group="Mesh size", +# ) + cpacs_inout.add_input( - var_name="fuselage_mesh_size", + var_name="fuselage_mesh_size_factor", var_type=float, - default_value=0.3, - unit="[m]", + default_value=1, + unit="1", descr="Value assigned for the fuselage surfaces mesh size", - xpath=GMSH_MESH_SIZE_FUSELAGE_XPATH, + xpath=GMSH_MESH_SIZE_FACTOR_FUSELAGE_XPATH, gui=True, - gui_name="Fuselage", + gui_name="Fuselage factor", gui_group="Mesh size", ) +# cpacs_inout.add_input( +# var_name="wing_mesh_size", +# var_type=float, +# default_value=0.1, +# unit="[m]", +# descr="Value assigned for the wings surfaces mesh size", +# xpath=GMSH_MESH_SIZE_WINGS_XPATH, +# gui=True, +# gui_name="Wings", +# gui_group="Mesh size", +# ) + cpacs_inout.add_input( - var_name="wing_mesh_size", + var_name="wing_mesh_size_factor", var_type=float, - default_value=0.123, - unit="[m]", + default_value=1, + unit="1", descr="Value assigned for the wings surfaces mesh size", - xpath=GMSH_MESH_SIZE_WINGS_XPATH, + xpath=GMSH_MESH_SIZE_FACTOR_WINGS_XPATH, gui=True, - gui_name="Wings", + gui_name="Wings factor", gui_group="Mesh size", ) + cpacs_inout.add_input( var_name="engine_mesh_size", var_type=float, @@ -240,3 +267,19 @@ descr="Absolute path of the SU2 mesh", xpath=SU2MESH_XPATH, ) + +cpacs_inout.add_output( + var_name="fuselage_mesh_size", + var_type=float, + default_value=None, + descr="Value assigned for the fuselage surfaces mesh size", + xpath=GMSH_MESH_SIZE_FUSELAGE_XPATH, +) + +cpacs_inout.add_output( + var_name="wing_mesh_size", + var_type=float, + default_value=None, + descr="Value assigned for the wings surfaces mesh size", + xpath=GMSH_MESH_SIZE_WINGS_XPATH, +) \ No newline at end of file diff --git a/ceasiompy/CPACS2GMSH/cpacs2gmsh.py b/ceasiompy/CPACS2GMSH/cpacs2gmsh.py index ef15ccb81..046868395 100644 --- a/ceasiompy/CPACS2GMSH/cpacs2gmsh.py +++ b/ceasiompy/CPACS2GMSH/cpacs2gmsh.py @@ -37,8 +37,8 @@ GMSH_N_POWER_FIELD_XPATH, GMSH_INTAKE_PERCENT_XPATH, GMSH_MESH_SIZE_FARFIELD_XPATH, - GMSH_MESH_SIZE_FUSELAGE_XPATH, - GMSH_MESH_SIZE_WINGS_XPATH, + GMSH_MESH_SIZE_FACTOR_FUSELAGE_XPATH, + GMSH_MESH_SIZE_FACTOR_WINGS_XPATH, GMSH_MESH_SIZE_ENGINES_XPATH, GMSH_MESH_SIZE_PROPELLERS_XPATH, GMSH_OPEN_GUI_XPATH, @@ -82,8 +82,8 @@ def cpacs2gmsh(cpacs_path, cpacs_out_path): mesh_size_farfield = get_value_or_default(cpacs.tixi, GMSH_MESH_SIZE_FARFIELD_XPATH, 25) n_power_factor = get_value_or_default(cpacs.tixi, GMSH_N_POWER_FACTOR_XPATH, 2) n_power_field = get_value_or_default(cpacs.tixi, GMSH_N_POWER_FIELD_XPATH, 0.9) - mesh_size_fuselage = get_value_or_default(cpacs.tixi, GMSH_MESH_SIZE_FUSELAGE_XPATH, 0.4) - mesh_size_wings = get_value_or_default(cpacs.tixi, GMSH_MESH_SIZE_WINGS_XPATH, 0.23) + fuselage_mesh_size_factor = get_value_or_default(cpacs.tixi, GMSH_MESH_SIZE_FACTOR_FUSELAGE_XPATH,1) + wing_mesh_size_factor = get_value_or_default(cpacs.tixi, GMSH_MESH_SIZE_FACTOR_WINGS_XPATH,1) mesh_size_engines = get_value_or_default(cpacs.tixi, GMSH_MESH_SIZE_ENGINES_XPATH, 0.23) mesh_size_propellers = get_value_or_default(cpacs.tixi, GMSH_MESH_SIZE_PROPELLERS_XPATH, 0.23) refine_factor = get_value_or_default(cpacs.tixi, GMSH_REFINE_FACTOR_XPATH, 7.0) @@ -96,6 +96,7 @@ def cpacs2gmsh(cpacs_path, cpacs_out_path): export_brep(cpacs, brep_dir, (intake_percent, exhaust_percent)) mesh_path, _ = generate_gmsh( cpacs, + cpacs_path, brep_dir, results_dir, open_gmsh=open_gmsh, @@ -104,8 +105,8 @@ def cpacs2gmsh(cpacs_path, cpacs_out_path): mesh_size_farfield=mesh_size_farfield, n_power_factor=n_power_factor, n_power_field=n_power_field, - mesh_size_fuselage=mesh_size_fuselage, - mesh_size_wings=mesh_size_wings, + fuselage_mesh_size_factor=fuselage_mesh_size_factor, + wing_mesh_size_factor=wing_mesh_size_factor, mesh_size_engines=mesh_size_engines, mesh_size_propellers=mesh_size_propellers, refine_factor=refine_factor, diff --git a/ceasiompy/CPACS2GMSH/func/generategmesh.py b/ceasiompy/CPACS2GMSH/func/generategmesh.py index a07c18e16..47b6b5aa3 100644 --- a/ceasiompy/CPACS2GMSH/func/generategmesh.py +++ b/ceasiompy/CPACS2GMSH/func/generategmesh.py @@ -55,6 +55,8 @@ from ceasiompy.utils.ceasiomlogger import get_logger from ceasiompy.utils.ceasiompyutils import get_part_type +from ceasiompy.CPACS2GMSH.func.mesh_sizing import fuselage_size, wings_size + log = get_logger() @@ -557,6 +559,7 @@ def control_disk_actuator_normal(): def generate_gmsh( cpacs, + cpacs_path, brep_dir, results_dir, open_gmsh=False, @@ -565,8 +568,8 @@ def generate_gmsh( mesh_size_farfield=25, n_power_factor=2, n_power_field=0.9, - mesh_size_fuselage=0.4, - mesh_size_wings=0.23, + fuselage_mesh_size_factor=1, + wing_mesh_size_factor=1, mesh_size_engines=0.23, mesh_size_propellers=0.23, refine_factor=7.0, @@ -924,13 +927,23 @@ def generate_gmsh( # Thus be sure to define mesh size in a certain order to control # the size of the points on boundaries. + fuselage_maxlen, fuselage_minlen = fuselage_size(cpacs_path) + log.info(f"fuselage_minlen={fuselage_minlen}") + mesh_size_fuselage = fuselage_mesh_size_factor * fuselage_minlen + log.info(f"mesh_size_fuselage={mesh_size_fuselage}") + + wing_maxlen, wing_minlen = wings_size(cpacs_path) + log.info(f"wing_minlen={wing_minlen}") + mesh_size_wing = wing_mesh_size_factor * wing_minlen + log.info(f"mesh_size_wing={mesh_size_wing}") + for part in aircraft_parts: if part.part_type == "fuselage": part.mesh_size = mesh_size_fuselage gmsh.model.mesh.setSize(part.points, part.mesh_size) gmsh.model.setColor(part.surfaces, *MESH_COLORS[part.part_type], recursive=False) elif part.part_type in ["wing", "pylon"]: - part.mesh_size = mesh_size_wings + part.mesh_size = mesh_size_wing gmsh.model.mesh.setSize(part.points, part.mesh_size) gmsh.model.setColor(part.surfaces, *MESH_COLORS[part.part_type], recursive=False) elif part.part_type == "engine": @@ -967,7 +980,7 @@ def generate_gmsh( final_domain.volume_tag, aircraft, part, - mesh_size_wings, + mesh_size_wing, refine=refine_factor, refine_truncated=refine_truncated, ) diff --git a/ceasiompy/CPACS2GMSH/func/mesh_sizing.py b/ceasiompy/CPACS2GMSH/func/mesh_sizing.py new file mode 100644 index 000000000..2cabcb122 --- /dev/null +++ b/ceasiompy/CPACS2GMSH/func/mesh_sizing.py @@ -0,0 +1,308 @@ +""" +CEASIOMpy: Conceptual Aircraft Design Software + +Developed by CFS ENGINEERING, 1015 Lausanne, Switzerland + +This script contains different functions to classify and manipulate wing elements + +Python version: >=3.8 + +| Author: Giacomo Benedetti +| Creation: 2023-11-20 + + +""" + + +# ================================================================================================= +# IMPORTS +# ================================================================================================= + +import math +from ceasiompy.CPACS2SUMO.func.getprofile import get_profile_coord +from ceasiompy.utils.ceasiomlogger import get_logger + +from ceasiompy.utils.generalclasses import Transformation + +from ceasiompy.utils.commonxpath import FUSELAGES_XPATH, WINGS_XPATH +from cpacspy.cpacsfunctions import open_tixi + +log = get_logger() + +# ================================================================================================= +# FUNCTIONS +# ================================================================================================= + +def fuselage_size(cpacs_path): + + tixi = open_tixi(cpacs_path) + if tixi.checkElement(FUSELAGES_XPATH): + fus_cnt = tixi.getNamedChildrenCount(FUSELAGES_XPATH, "fuselage") + for i_fus in range(fus_cnt): + fus_xpath = FUSELAGES_XPATH + "/fuselage[" + str(i_fus + 1) + "]" + fus_uid = tixi.getTextAttribute(fus_xpath, "uID") + fus_transf = Transformation() + fus_transf.get_cpacs_transf(tixi, fus_xpath) + + # Positionings + if tixi.checkElement(fus_xpath + "/positionings"): + pos_cnt = tixi.getNamedChildrenCount(fus_xpath + "/positionings", "positioning") + + pos_x_list = [] + pos_y_list = [] + pos_z_list = [] + from_sec_list = [] + to_sec_list = [] + + for i_pos in range(pos_cnt): + pos_xpath = fus_xpath + "/positionings/positioning[" + str(i_pos + 1) + "]" + + length = tixi.getDoubleElement(pos_xpath + "/length") + sweep_deg = tixi.getDoubleElement(pos_xpath + "/sweepAngle") + sweep = math.radians(sweep_deg) + dihedral_deg = tixi.getDoubleElement(pos_xpath + "/dihedralAngle") + dihedral = math.radians(dihedral_deg) + + # Get the corresponding translation of each positioning + pos_x_list.append(length * math.sin(sweep)) + pos_y_list.append(length * math.cos(dihedral) * math.cos(sweep)) + pos_z_list.append(length * math.sin(dihedral) * math.cos(sweep)) + # Get which section are connected by the positioning + if tixi.checkElement(pos_xpath + "/fromSectionUID"): + from_sec = tixi.getTextElement(pos_xpath + "/fromSectionUID") + else: + from_sec = "" + from_sec_list.append(from_sec) + if tixi.checkElement(pos_xpath + "/toSectionUID"): + to_sec = tixi.getTextElement(pos_xpath + "/toSectionUID") + else: + to_sec = "" + to_sec_list.append(to_sec) + + for j_pos in range(pos_cnt): + if from_sec_list[j_pos] == "": + prev_pos_x = 0 + prev_pos_y = 0 + prev_pos_z = 0 + + elif from_sec_list[j_pos] == to_sec_list[j_pos - 1]: + prev_pos_x = pos_x_list[j_pos - 1] + prev_pos_y = pos_y_list[j_pos - 1] + prev_pos_z = pos_z_list[j_pos - 1] + + else: + index_prev = to_sec_list.index(from_sec_list[j_pos]) + prev_pos_x = pos_x_list[index_prev] + prev_pos_y = pos_y_list[index_prev] + prev_pos_z = pos_z_list[index_prev] + + pos_x_list[j_pos] += prev_pos_x + pos_y_list[j_pos] += prev_pos_y + pos_z_list[j_pos] += prev_pos_z + + else: + log.error('No "positionings" have been found!') + pos_cnt = 0 + + # Sections + sec_cnt = tixi.getNamedChildrenCount(fus_xpath + "/sections", "section") + + if pos_cnt == 0: + pos_x_list = [0.0] * sec_cnt + pos_y_list = [0.0] * sec_cnt + pos_z_list = [0.0] * sec_cnt + + # print(pos_cnt) + body_frm_height_values = [] + body_frm_width_values = [] + + for i_sec in range(sec_cnt): + sec_xpath = fus_xpath + "/sections/section[" + str(i_sec + 1) + "]" + sec_uid = tixi.getTextAttribute(sec_xpath, "uID") + + sec_transf = Transformation() + sec_transf.get_cpacs_transf(tixi, sec_xpath) + + # Elements + elem_cnt = tixi.getNamedChildrenCount(sec_xpath + "/elements", "element") + + for i_elem in range(elem_cnt): + elem_xpath = sec_xpath + "/elements/element[" + str(i_elem + 1) + "]" + elem_uid = tixi.getTextAttribute(elem_xpath, "uID") + elem_transf = Transformation() + elem_transf.get_cpacs_transf(tixi, elem_xpath) + + # Fuselage profiles + prof_uid = tixi.getTextElement(elem_xpath + "/profileUID") + prof_vect_x, prof_vect_y, prof_vect_z = get_profile_coord(tixi, prof_uid) + + prof_size_y = (max(prof_vect_y) - min(prof_vect_y)) / 2 + prof_size_z = (max(prof_vect_z) - min(prof_vect_z)) / 2 + + prof_vect_y[:] = [y / prof_size_y for y in prof_vect_y] + prof_vect_z[:] = [z / prof_size_z for z in prof_vect_z] + + prof_min_y = min(prof_vect_y) + prof_min_z = min(prof_vect_z) + + prof_vect_y[:] = [y - 1 - prof_min_y for y in prof_vect_y] + prof_vect_z[:] = [z - 1 - prof_min_z for z in prof_vect_z] + + pos_y_list[i_sec] += ((1 + prof_min_y) * prof_size_y) * elem_transf.scaling.y + pos_z_list[i_sec] += ((1 + prof_min_z) * prof_size_z) * elem_transf.scaling.z + + body_frm_height = ( + prof_size_z + * 2 + * elem_transf.scaling.z + * sec_transf.scaling.z + * fus_transf.scaling.z + ) + + body_frm_width = ( + prof_size_y + * 2 + * elem_transf.scaling.y + * sec_transf.scaling.y + * fus_transf.scaling.y + ) + + body_frm_height_values.append(body_frm_height) + body_frm_width_values.append(body_frm_width) + + # Extrapolate mesh values + circ_list = [] + min_radius = 10e6 + refine_factor = 1 + + for height, width in zip(body_frm_height_values, body_frm_width_values): + # Calculate the sum of squares for each element in the lists + circ_list.append(2 * math.pi * math.sqrt((height**2 + width**2) / 2)) + + # Get overall minimum radius (semi-minor axis for ellipse) + min_radius = min(min_radius, height, width) + + mean_circ = sum(circ_list) / len(circ_list) + + # Calculate mesh parameters from inputs and geometry + # In SUMO, fuselage_minlen is min_radius/2, but sometimes it leads to meshing errors + fuselage_maxlen = (0.08 * mean_circ) * refine_factor + fuselage_minlen = min(0.1 * fuselage_maxlen, min_radius / 2) * refine_factor + + log.info(f"maxlen={fuselage_maxlen:.3f}") + log.info(f"fuselage_minlen={fuselage_minlen:.3f}") + + return fuselage_maxlen, fuselage_minlen + + +def wings_size(cpacs_path): + + tixi = open_tixi(cpacs_path) + if tixi.checkElement(WINGS_XPATH): + wing_cnt = tixi.getNamedChildrenCount(WINGS_XPATH, "wing") + log.warning(str(wing_cnt) + " wings have been found.") + + chord_list = [] + + for i_wing in range(wing_cnt): + wing_xpath = WINGS_XPATH + "/wing[" + str(i_wing + 1) + "]" + wing_uid = tixi.getTextAttribute(wing_xpath, "uID") + wing_transf = Transformation() + wing_transf.get_cpacs_transf(tixi, wing_xpath) + + # Positionings + if tixi.checkElement(wing_xpath + "/positionings"): + pos_cnt = tixi.getNamedChildrenCount(wing_xpath + "/positionings", "positioning") + + pos_x_list = [] + pos_y_list = [] + pos_z_list = [] + from_sec_list = [] + to_sec_list = [] + + for i_pos in range(pos_cnt): + pos_xpath = wing_xpath + "/positionings/positioning[" + str(i_pos + 1) + "]" + + length = tixi.getDoubleElement(pos_xpath + "/length") + sweep_deg = tixi.getDoubleElement(pos_xpath + "/sweepAngle") + sweep = math.radians(sweep_deg) + dihedral_deg = tixi.getDoubleElement(pos_xpath + "/dihedralAngle") + dihedral = math.radians(dihedral_deg) + + # Get the corresponding translation of each positioning + pos_x_list.append(length * math.sin(sweep)) + pos_y_list.append(length * math.cos(dihedral) * math.cos(sweep)) + pos_z_list.append(length * math.sin(dihedral) * math.cos(sweep)) + # Get which section are connected by the positioning + if tixi.checkElement(pos_xpath + "/fromSectionUID"): + from_sec = tixi.getTextElement(pos_xpath + "/fromSectionUID") + else: + from_sec = "" + from_sec_list.append(from_sec) + if tixi.checkElement(pos_xpath + "/toSectionUID"): + to_sec = tixi.getTextElement(pos_xpath + "/toSectionUID") + else: + to_sec = "" + to_sec_list.append(to_sec) + + for j_pos in range(pos_cnt): + if from_sec_list[j_pos] == "": + prev_pos_x = 0 + prev_pos_y = 0 + prev_pos_z = 0 + + elif from_sec_list[j_pos] == to_sec_list[j_pos - 1]: + prev_pos_x = pos_x_list[j_pos - 1] + prev_pos_y = pos_y_list[j_pos - 1] + prev_pos_z = pos_z_list[j_pos - 1] + + else: + index_prev = to_sec_list.index(from_sec_list[j_pos]) + prev_pos_x = pos_x_list[index_prev] + prev_pos_y = pos_y_list[index_prev] + prev_pos_z = pos_z_list[index_prev] + + pos_x_list[j_pos] += prev_pos_x + pos_y_list[j_pos] += prev_pos_y + pos_z_list[j_pos] += prev_pos_z + + else: + log.warning('No "positionings" have been found!') + pos_cnt = 0 + + # Sections + sec_cnt = tixi.getNamedChildrenCount(wing_xpath + "/sections", "section") + + if pos_cnt == 0: + pos_x_list = [0.0] * sec_cnt + pos_y_list = [0.0] * sec_cnt + pos_z_list = [0.0] * sec_cnt + + for i_sec in range(sec_cnt): + sec_xpath = wing_xpath + "/sections/section[" + str(i_sec + 1) + "]" + sec_uid = tixi.getTextAttribute(sec_xpath, "uID") + + sec_transf = Transformation() + sec_transf.get_cpacs_transf(tixi, sec_xpath) + + # Salva la corda di ogni sezione nella lista + chord_list.append(sec_transf.scaling.x) + + refine_factor = 1 + refine_level = 0 + + ref_chord = sum(chord_list) / len(chord_list) + + # Calculate mesh parameter from inputs and geometry + wing_maxlen = (0.15 * ref_chord) * refine_factor + wing_minlen = (0.08 * wing_maxlen) * refine_factor + # in sumo it is 0.08*wing_maxlen or 0.7*min leading edge radius... + # Default value in SUMO + lerfactor = 1 / 4.0 + terfactor = 1 / 4. + if refine_level > 1: + lerfactor = 1 / (4.0 + 0.5 * (refine_level - 1)) + terfactor = 1 / (4.0 + 0.5 * (refine_level - 1)) + + return wing_maxlen, wing_minlen + diff --git a/ceasiompy/utils/commonxpath.py b/ceasiompy/utils/commonxpath.py index b69f3d441..2aa10a5f7 100644 --- a/ceasiompy/utils/commonxpath.py +++ b/ceasiompy/utils/commonxpath.py @@ -98,7 +98,9 @@ GMSH_N_POWER_FIELD_XPATH = GMSH_XPATH + "/n_power_field" GMSH_MESH_SIZE_FARFIELD_XPATH = GMSH_XPATH + "/mesh_size/farfield" GMSH_MESH_SIZE_FUSELAGE_XPATH = GMSH_XPATH + "/mesh_size/fuselage" +GMSH_MESH_SIZE_FACTOR_FUSELAGE_XPATH = GMSH_MESH_SIZE_FUSELAGE_XPATH + "/factor" GMSH_MESH_SIZE_WINGS_XPATH = GMSH_XPATH + "/mesh_size/wings" +GMSH_MESH_SIZE_FACTOR_WINGS_XPATH = GMSH_MESH_SIZE_WINGS_XPATH + '/factor' GMSH_MESH_SIZE_ENGINES_XPATH = GMSH_XPATH + "/mesh_size/engines" GMSH_MESH_SIZE_PROPELLERS_XPATH = GMSH_XPATH + "/mesh_size/propellers" GMSH_REFINE_FACTOR_XPATH = GMSH_XPATH + "/refine_factor"