integration test with new framat

Author: LeonDeligny

ceasiompy/AeroFrame_new/aeroframe_run.py

@@ -27,6 +27,7 @@
 
 from ceasiompy.PyAVL.pyavl import main as run_avl
 from ceasiompy.utils.ceasiompyutils import call_main
+from ceasiompy.PyAVL.func.utils import create_case_dir
 from ceasiompy.PyAVL.func.plot import convert_ps_to_pdf
 from ceasiompy.utils.ceasiompyutils import get_aeromap_conditions
 from cpacspy.cpacsfunctions import (
@@ -43,7 +44,7 @@
     plot_deformed_wing
 )
 from ceasiompy.AeroFrame_new.func.config import (
-    read_AVL_fe_file,
+    read_avl_fe_file,
     create_framat_model,
     get_material_properties,
     get_section_properties,
@@ -71,7 +72,12 @@
 
 
 # TODO: Reduce complexity
-def aeroelastic_loop(cpacs_path, CASE_PATH, q, xyz, fxyz):
+def aeroelastic_loop(
+    cpacs,
+    results_dir,
+    case_dir_path,
+    q, xyz, fxyz
+):
     """Function to execute the aeroelastic-loop.
 
     Function 'aeroelastic_loop' executes the aeroelastic-loop,
@@ -90,12 +96,10 @@ def aeroelastic_loop(cpacs_path, CASE_PATH, q, xyz, fxyz):
         res (list): residual of the mid-chord wing tip deflection.
     """
 
-    cpacs = CPACS(cpacs_path)
-
-    AVL_ITER1_PATH = Path(CASE_PATH, "Iteration_1", "AVL")
+    AVL_ITER1_PATH = Path(case_dir_path, "Iteration_1", "AVL")
 
     # Get the path to the undeformed/initial AVL geometry
-    for path in AVL_ITER1_PATH.glob('*.avl'):
+    for path in results_dir.glob('*.avl'):
         AVL_UNDEFORMED_PATH = path
     AVL_UNDEFORMED_COMMAND = Path(AVL_ITER1_PATH, "avl_commands.txt")
 
@@ -111,10 +115,10 @@ def aeroelastic_loop(cpacs_path, CASE_PATH, q, xyz, fxyz):
         wg_center_z_list,
         wg_chord_list,
         wg_scaling
-    ) = compute_cross_section(cpacs_path)
+    ) = compute_cross_section(cpacs)
 
     # Get the properties of the wing material and the number of beam nodes to use
-    young_modulus, shear_modulus, material_density = get_material_properties(cpacs_path)
+    young_modulus, shear_modulus, material_density = get_material_properties(cpacs)
     N_beam = get_value_or_default(cpacs.tixi, FRAMAT_MESH_XPATH + "/NumberNodes", 8)
 
     # Define the coordinates of the wing root and tip
@@ -155,7 +159,7 @@ def aeroelastic_loop(cpacs_path, CASE_PATH, q, xyz, fxyz):
                         Zle_array[-1]])
 
     # Get cross-section properties from CPACS file (if constants)
-    area_const, Ix_const, Iy_const = get_section_properties(cpacs_path)
+    area_const, Ix_const, Iy_const = get_section_properties(cpacs)
 
     if area_const > 0:
         area_list = area_const * np.ones(len(wg_center_y_list))
@@ -189,10 +193,10 @@ def aeroelastic_loop(cpacs_path, CASE_PATH, q, xyz, fxyz):
         log.info("")
         log.info(f"----- FramAT: Deformation {n_iter} -----")
 
-        Path(CASE_PATH, f"Iteration_{n_iter + 1}", "AVL").mkdir(parents=True, exist_ok=True)
-        Path(CASE_PATH, f"Iteration_{n_iter}", "FramAT").mkdir(parents=True, exist_ok=True)
-        AVL_ITER_PATH = Path(CASE_PATH, f"Iteration_{n_iter + 1}", "AVL")
-        FRAMAT_ITER_PATH = Path(CASE_PATH, f"Iteration_{n_iter}", "FramAT")
+        Path(case_dir_path, f"Iteration_{n_iter + 1}", "AVL").mkdir(parents=True, exist_ok=True)
+        Path(case_dir_path, f"Iteration_{n_iter}", "FramAT").mkdir(parents=True, exist_ok=True)
+        AVL_ITER_PATH = Path(case_dir_path, f"Iteration_{n_iter + 1}", "AVL")
+        FRAMAT_ITER_PATH = Path(case_dir_path, f"Iteration_{n_iter}", "FramAT")
         AVL_DEFORMED_PATH = Path(AVL_ITER_PATH, "deformed.avl")
         AVL_DEFORMED_COMMAND = Path(AVL_ITER_PATH, "avl_commands.txt")
 
@@ -217,7 +221,7 @@ def aeroelastic_loop(cpacs_path, CASE_PATH, q, xyz, fxyz):
                                    Iy_profile,
                                    chord_profile,
                                    twist_profile,
-                                   CASE_PATH,
+                                   case_dir_path,
                                    AVL_UNDEFORMED_PATH,
                                    wg_scaling)
 
@@ -283,7 +287,7 @@ def aeroelastic_loop(cpacs_path, CASE_PATH, q, xyz, fxyz):
 
         # Read the new forces file to extract the loads
         FE_PATH = Path(AVL_ITER_PATH, "fe.txt")
-        _, _, _, xyz_list, p_xyz_list, _ = read_AVL_fe_file(FE_PATH, plot=False)
+        _, _, _, xyz_list, p_xyz_list, _ = read_avl_fe_file(FE_PATH, plot=False)
 
         xyz = xyz_list[0]
         fxyz = np.array(p_xyz_list[0]) * q
@@ -368,7 +372,7 @@ def main(cpacs: CPACS, results_dir: Path) -> None:
 
     # First AVL run
     # You need to first load the default values of AVL
-    # Then you add back the ones that you wanted to specify...
+    # Then you add back the ones that you wanted to specify.
     run_avl(cpacs, results_dir)
 
     for i_case, _ in enumerate(alt_list):
@@ -380,38 +384,47 @@ def main(cpacs: CPACS, results_dir: Path) -> None:
         Atm = Atmosphere(alt)
         fluid_density = Atm.density[0]
         fluid_velocity = Atm.speed_of_sound[0] * mach
-        q = 0.5 * fluid_density * fluid_velocity**2
-
-        case_dir_name = (
-            f"Case{str(i_case).zfill(2)}_alt{alt}_mach{round(mach, 2)}"
-            f"_aoa{round(aoa, 1)}_aos{round(aos, 1)}"
+        rho = 0.5 * fluid_density * fluid_velocity**2
+
+        case_dir_path = create_case_dir(
+            results_dir,
+            i_case,
+            alt,
+            mach=mach,
+            aoa=aoa,
+            aos=aos,
+            q=0.0,
+            p=0.0,
+            r=0.0,
         )
+        Path(case_dir_path, "Iteration_1", "AVL").mkdir(parents=True, exist_ok=True)
 
-        CASE_PATH = Path(results_dir, case_dir_name)
-        Path(CASE_PATH, "Iteration_1", "AVL").mkdir(parents=True, exist_ok=True)
-        AVL_ITER_PATH = Path(CASE_PATH, "Iteration_1", "AVL")
+        avl_iter_path = Path(case_dir_path, "Iteration_1", "AVL")
 
-        for file_path in CASE_PATH.iterdir():
+        for file_path in case_dir_path.iterdir():
             if file_path.is_file():
-                shutil.move(str(file_path), str(AVL_ITER_PATH / file_path.name))
+                shutil.move(str(file_path), str(avl_iter_path / file_path.name))
 
-        FE_PATH = Path(AVL_ITER_PATH, "fe.txt")
-        _, _, _, xyz_list, p_xyz_list, _ = read_AVL_fe_file(FE_PATH, plot=False)
+        fe_path = Path(avl_iter_path, "fe.txt")
+        _, _, _, xyz_list, p_xyz_list, _ = read_avl_fe_file(fe_path, plot=False)
 
-        f_xyz_array = np.array(p_xyz_list) * q
+        f_xyz_array = np.array(p_xyz_list) * rho
 
         # Start the aeroelastic loop
-        tip_deflection, residuals = aeroelastic_loop(cpacs_path,
-                                                     CASE_PATH,
-                                                     q,
-                                                     xyz_list[0],
-                                                     f_xyz_array[0])
+        tip_deflection, residuals = aeroelastic_loop(
+            cpacs,
+            results_dir,
+            case_dir_path,
+            rho,
+            xyz_list[0],
+            f_xyz_array[0],
+        )
 
         # Write results in CPACS out
         create_branch(tixi, FRAMAT_RESULTS_XPATH + "/TipDeflection")
         tixi.updateDoubleElement(FRAMAT_RESULTS_XPATH + "/TipDeflection", tip_deflection[-1], "%g")
 
-        plot_convergence(tip_deflection, residuals, wkdir=CASE_PATH)
+        plot_convergence(tip_deflection, residuals, wkdir=case_dir_path)
 
 
 # =================================================================================================

ceasiompy/AeroFrame_new/func/config.py

@@ -16,7 +16,6 @@
 # ==============================================================================
 
 import re
-import math
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
@@ -25,6 +24,10 @@
 from cpacspy.cpacsfunctions import open_tixi
 from cpacspy.cpacsfunctions import get_value_or_default
 from ceasiompy.CPACS2SUMO.func.getprofile import get_profile_coord
+from ceasiompy.utils.geometryfunctions import (
+    sum_points,
+    get_positionings,
+)
 from ceasiompy.utils.mathsfunctions import (
     euler2fix,
     rotate_points,
@@ -39,8 +42,8 @@
 from scipy import interpolate
 from cpacspy.cpacspy import CPACS
 from ceasiompy.utils.generalclasses import (
+    Point,
     Transformation,
-    SimpleNamespace,
 )
 
 from ceasiompy import log
@@ -173,7 +176,7 @@ def parse_AVL_surface(string):
     return surface_name, nspanwise, nchord_strip, xyz, p_xyz, slope_list
 
 
-def read_AVL_fe_file(FE_PATH, plot=False):
+def read_avl_fe_file(FE_PATH, plot=False):
     """Function to read AVL 'fe.txt' element force file,
     and extract the aerodynamic loads calling the function
     'parse_AVL_surface'.
@@ -355,7 +358,7 @@ def create_framat_model(young_modulus, shear_modulus, material_density,
     return model
 
 
-def get_material_properties(cpacs_path):
+def get_material_properties(cpacs):
     """Function to read the material properties for structural
     calculations.
 
@@ -372,20 +375,19 @@ def get_material_properties(cpacs_path):
         material_density (float): Density of the material [kg/m^3].
 
     """
-    cpacs = CPACS(cpacs_path)
-
-    young_modulus = get_value_or_default(cpacs.tixi,
+    tixi = cpacs.tixi
+    young_modulus = get_value_or_default(tixi,
                                          FRAMAT_MATERIAL_XPATH + "/YoungModulus", 70)
 
-    shear_modulus = get_value_or_default(cpacs.tixi,
+    shear_modulus = get_value_or_default(tixi,
                                          FRAMAT_MATERIAL_XPATH + "/ShearModulus", 27)
 
-    material_density = get_value_or_default(cpacs.tixi, FRAMAT_MATERIAL_XPATH + "/Density", 1960)
+    material_density = get_value_or_default(tixi, FRAMAT_MATERIAL_XPATH + "/Density", 1960)
 
     return young_modulus, shear_modulus, material_density
 
 
-def get_section_properties(cpacs_path):
+def get_section_properties(cpacs):
     """Function reads the cross-section properties for structural
     calculations.
 
@@ -402,11 +404,11 @@ def get_section_properties(cpacs_path):
         Iy (float): second moment of area about the x-axis [m^4].
 
     """
-    cpacs = CPACS(cpacs_path)
+    tixi = cpacs.tixi
 
-    area = get_value_or_default(cpacs.tixi, FRAMAT_SECTION_XPATH + "/Area", -1)
-    Ix = get_value_or_default(cpacs.tixi, FRAMAT_SECTION_XPATH + "/Ix", -1)
-    Iy = get_value_or_default(cpacs.tixi, FRAMAT_SECTION_XPATH + "/Iy", -1)
+    area = get_value_or_default(tixi, FRAMAT_SECTION_XPATH + "/Area", -1)
+    Ix = get_value_or_default(tixi, FRAMAT_SECTION_XPATH + "/Ix", -1)
+    Iy = get_value_or_default(tixi, FRAMAT_SECTION_XPATH + "/Iy", -1)
 
     return area, Ix, Iy
 
@@ -613,7 +615,7 @@ def compute_distance_and_moment(row):
 
 
 # TODO: Reduce complexity
-def compute_cross_section(cpacs_path):
+def compute_cross_section(cpacs):
     """
     Function 'compute_cross_section' computes the area, the second moments of area,
     and additional geometric properties of each cross-section of the wing.
@@ -638,7 +640,7 @@ def compute_cross_section(cpacs_path):
         wg_chord_list (list): list of the chord length of each cross-section [m].
 
     """
-    tixi = open_tixi(cpacs_path)
+    tixi = cpacs.tixi
 
     # Wing(s) ------------------------------------------------------------------
     if tixi.checkElement(WINGS_XPATH):
@@ -681,76 +683,7 @@ def compute_cross_section(cpacs_path):
                       round(wing_transf.scaling.z, 3)]
 
         # Positionings
-        if tixi.checkElement(wing_xpath + "/positionings"):
-            pos_cnt = tixi.getNamedChildrenCount(wing_xpath + "/positionings", "positioning")
-            log.info(str(pos_cnt) + ' "positioning" has been found : ')
-
-            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)
-
-            # Re-loop though the positioning to re-order them
-            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")
-        log.info("    -" + str(sec_cnt) + " wing sections have been found")
-
-        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
-
+        sec_cnt, pos_x_list, pos_y_list, pos_z_list = get_positionings(tixi, wing_xpath, "wing")
         for i_sec in range(sec_cnt):
             sec_xpath = wing_xpath + "/sections/section[" + str(i_sec + 1) + "]"
             sec_uid = tixi.getTextAttribute(sec_xpath, "uID")
@@ -802,13 +735,12 @@ def compute_cross_section(cpacs_path):
 
                 # Add rotation from element and sections
                 # Adding the two angles: Maybe not work in every case!!!
-                add_rotation = SimpleNamespace()
-                add_rotation.x = elem_transf.rotation.x + \
-                    sec_transf.rotation.x + wg_sk_transf.rotation.x
-                add_rotation.y = elem_transf.rotation.y + \
-                    sec_transf.rotation.y + wg_sk_transf.rotation.y
-                add_rotation.z = elem_transf.rotation.z + \
-                    sec_transf.rotation.z + wg_sk_transf.rotation.z
+                x, y, z = sum_points(
+                    elem_transf.rotation,
+                    sec_transf.rotation,
+                    sec_transf.rotation
+                )
+                add_rotation = Point(x, y, z)
 
                 # Get Section rotation
                 wg_sec_rot = euler2fix(add_rotation)