Add body forces to api and strip forces fixes

Makefile

@@ -21,10 +21,6 @@ clean:
 	rm -fr src/build/*.o
 	rm -fr src/build/*.a
 	rm -fr src/build/*.so
-	rm -fr src_cs/build/*.mod
-	rm -fr src_cs/build/*.o
-	rm -fr src_cs/build/*.a
-	rm -fr src_cs/build/*.so
 	rm -fr optvl/*.so
 	rm -f *~ config.mk;
 

docs/optvl_api.md

@@ -50,13 +50,13 @@ The commands from the oper and mode menus are available
 | get shear moment distribution | VM | ovl.get_strip_forces() |
 | get control surface derivatives (e.g. dCL/dElevator)| ST | ovl.get_control_derivs() |
 | get stability derivatives | ST | ovl.get_stab_derivs()|
-| get stability derivatives in the body axis| SB | - |
+| get stability derivatives in the body axis| SB | ovl.get_body_axis_derivs() |
 | get/set reference data | RE | ovl.get_reference_data/set_reference_data()|
-| get/set  design variables specified in AVL file | DE | -|
+| get/set twist design variables specified in AVL file | DE | not supported (set twist directly)|
 | get surface forces | FN | ovl.get_surface_forces() |
 | get force distribution on strips| FS| ovl.get_strip_forces() |
 | get force distribution on elements | FE | - |
-| get body forces| FB | -|
+| get body forces| FB | ovl.get_body_forces()|
 | get high moments| HM | ovl.get_hinge_moments() |
 
 

examples/run_analysis_body.py

@@ -1,16 +1,19 @@
 from optvl import OVLSolver
 import numpy as np
 
-ovl_solver = OVLSolver(geo_file="../geom_files/supra.avl", debug=False)
+ovl = OVLSolver(geo_file="../geom_files/supra.avl", debug=False)
 # set the angle of attack
-ovl_solver.set_parameter("Mach", 0.0)
+ovl.set_parameter("Mach", 0.0)
+ovl.execute_run()
 
+body_names = ovl.get_body_names()
 print("----------------- alpha sweep ----------------")
-print("   Angle        Cl           Cd          Cdi          Cdv          Cm")
+print("   Angle        CL           Cd        Cm        Cl")
 for alpha in range(10):
-    ovl_solver.set_variable("alpha", alpha)
-    ovl_solver.execute_run()
-    run_data = ovl_solver.get_total_forces()
-    print(
-        f' {alpha:10.6f}   {run_data["CL"]:10.6f}   {run_data["CD"]:10.6f}   {run_data["CDi"]:10.6f}   {run_data["CDv"]:10.6f}   {run_data["Cm"]:10.6f}'
-    )
+    ovl.set_variable("alpha", alpha)
+    ovl.execute_run()
+    run_data = ovl.get_body_forces()
+    for body_name in body_names:
+        print(
+            f' {alpha:10.2f}   {run_data[body_name]["CL"]:10.6e}   {run_data[body_name]["CD"]:10.6e}   {run_data[body_name]["Cm"]:10.6e}  {run_data[body_name]["Cl"]:10.6e}'
+        )

geom_files/supra.avl

@@ -18,7 +18,7 @@ TRANSLATE
 0.0  0.0  0.0
 #
 BFIL
-fuseSupra.dat
+../geom_files/fuseSupra.dat
 
 #==============================================================
 SURFACE

meson.build

@@ -1,7 +1,8 @@
 project(
   'optvl',
   'c',
-  version: '2.1.1',
+  # don't forget to change the value in pyproject.toml
+  version: '2.2.0',
   license: 'GPL-3.0',
   meson_version: '>= 0.64.0',
   default_options: [

optvl/optvl_class.py

@@ -114,6 +114,19 @@ class OVLSolver(object):
         # spanwise efficiency
         "e": ["CASE_R", "SPANEF"],
     }
+
+    body_forces_var_to_fort_var = {
+        'length':        ['BODY_R', 'ELBDY'],
+        'surface area':  ['BODY_R', 'SRFBDY'],
+        'volume':        ['BODY_R', 'VOLBDY'],
+        'CL':            ['BODY_R', 'CLBDY'],
+        'CD':            ['BODY_R', 'CDBDY'],
+        'Cm':            ['BODY_R', 'CMBDY', 1],
+        'CY':            ['BODY_R', 'CFBDY', 1],
+        'Cn':            ['BODY_R', 'CMBDYBAX', 2],
+        'Cl':            ['BODY_R', 'CMBDYBAX', 0]
+    }
+
     ref_var_to_fort_var = {
         "Sref": ["CASE_R", "SREF"],
         "Cref": ["CASE_R", "CREF"],
@@ -159,14 +172,15 @@ class OVLSolver(object):
         "CZ": ["SURF_R", "CFSURF", 2],
 
         # non-dimensionalized moments (body frame)
-        "Cl": ["SURF_R", "CMSURF", 0],
-        "Cm": ["SURF_R", "CMSURF", 1],
-        "Cn": ["SURF_R", "CMSURF", 2],
+        "Cl": ["SURF_R", "CMSURFBAX", 0],
+        "Cm": ["SURF_R", "CMSURFBAX", 1],
+        "Cn": ["SURF_R", "CMSURFBAX", 2],
 
         # forces non-dimentionalized by surface quantities
         # uses surface area instead of sref and takes moments about leading edge
         "CL surf" : ["SURF_R", "CL_LSRF"],
         "CD surf" : ["SURF_R", "CD_LSRF"],
+        "CDv surf" : ["SURF_R", "CDVSURF"],
         "CMLE_LSTRP surf" : ["SURF_R", "CMLE_LSRF"],
 
         #TODO: add CF_LSRF(3,NFMAX), CM_LSRF(3,NFMAX)
@@ -1096,7 +1110,7 @@ def set_section_naca(self, isec: int, isurf: int, nasec: int, naca: str, xfminma
         zf = np.zeros_like(xf)
         zf[xf < pos] = cam * (2.0 * pos * xf[xf < pos] - 1.0) * xf[xf < pos] / (pos**2)
         zf[xf > pos] = cam * ((1 - 2.0 * pos) + (2.0 * pos - xf[xf > pos]) * xf[xf > pos]) / (1 - pos) ** 2
-        theta = np.atan(slopes)
+        theta = np.arctan(slopes)
 
         # Set airfoil section
         self.set_avl_fort_arr("SURF_GEOM_I", "NASEC", nasec, slicer=(isurf, isec))
@@ -1503,6 +1517,29 @@ def get_total_forces(self) -> Dict[str, float]:
 
         return total_data
 
+    def get_body_forces(self) -> Dict[str, float]:
+        """Get the aerodynamic data for the last run case and return it as a dictionary.
+
+        Returns:
+            Dict[str, float]: Dictionary of aerodynamic data. The keys the aerodyanmic coefficients.
+        """
+
+        body_data = {}
+        body_names = self.get_body_names()
+
+        for body_name in body_names:
+            body_data[body_name] = {}
+
+        for key, avl_key in self.body_forces_var_to_fort_var.items():
+            val = self.get_avl_fort_arr(avl_key[0], avl_key[1])
+            for idx_body, body_name in enumerate(body_names):
+                if (len(avl_key) > 2):
+                    body_data[body_name][key] = val[idx_body,avl_key[2]]
+                else:
+                    body_data[body_name][key] = val[idx_body]
+
+        return body_data
+
     def get_control_stab_derivs(self) -> Dict[str, float]:
         """Get the control surface derivative data, i.e. dCL/dElevator,
         for the current analysis run
@@ -1897,7 +1934,7 @@ def get_strip_forces(self) -> Dict[str, Dict[str, np.ndarray]]:
             strip_data[surf_key]["area"] = strip_data[surf_key]["width"] * strip_data[surf_key]["chord"]
 
             # formula is directly from AVL
-            xcp = 0.25 - strip_data[surf_key]["Cm c/4"] / strip_data[surf_key]["CL perp"]
+            xcp = 0.25 - strip_data[surf_key]["Cm c/4"] / strip_data[surf_key]["CL strip"]
             strip_data[surf_key]["CP x/c"] = xcp
 
         # get length of along the surface of each strip
@@ -4190,11 +4227,121 @@ def add_mesh_plot(
                         }
                         axis.plot(pts[xaxis], pts[yaxis], mesh_style, color=color, alpha=0.7, linewidth=mesh_linewidth)
 
-    def plot_geom(self, axes=None):
+    def add_body_plot(
+        self,
+        axis,
+        xaxis: str = "x",
+        yaxis: str = "y",
+        color: str = "magenta",
+    ):
+        """Adds a plot of the body geometry to the axis
+
+        Args:
+            axis: axis to add the plot to
+            xaxis: what variable should be plotted on the x axis. Options are ['x', 'y', 'z']
+            yaxis: what variable should be plotted on the y-axis. Options are ['x', 'y', 'z']
+            color: what color should the body be plotted in
+        """
+        import numpy as np
+
+        # Check if any bodies exist
+        num_bodies = self.get_avl_fort_arr("CASE_I", "NBODY")
+        if num_bodies == 0:
+            return
+
+        # Get body data
+        nl = self.get_avl_fort_arr("BODY_I", "NL", slicer=slice(0, num_bodies))
+        lfrst = self.get_avl_fort_arr("BODY_I", "LFRST", slicer=slice(0, num_bodies))
+        print(lfrst)
+
+        # Calculate total number of body nodes needed
+        # We need lfrst[i] + nl[i] to include all nodes we might access
+        max_node_idx = max(lfrst[i] + nl[i] for i in range(num_bodies))
+
+        # Get body line coordinates and radii
+        rl = self.get_avl_fort_arr("VRTX_R", "RL", slicer=(slice(None), slice(0, max_node_idx)))
+        radl = self.get_avl_fort_arr("VRTX_R", "RADL", slicer=slice(0, max_node_idx))
+
+        # Number of points around each cross-section circle
+        kk = 51
+
+        # Plot each body
+        for n in range(num_bodies):
+            # First and last node indices for this body (Fortran uses 1-based indexing)
+            l1 = lfrst[n] - 1  # Convert to 0-based
+            ln = l1 + nl[n]
+
+            # Plot centerline
+            centerline_pts = {
+                "x": rl[l1:ln, 0],
+                "y": rl[l1:ln, 1],
+                "z": rl[l1:ln, 2],
+            }
+            axis.plot(centerline_pts[xaxis], centerline_pts[yaxis], "-", color=color, linewidth=1.5)
+
+            # Plot cross-section circles at each node
+            for l in range(l1, ln):
+                # Calculate tangent vector from neighboring nodes
+                lm = max(l - 1, l1)
+                lp = min(l + 1, ln - 1)
+
+                dxl = rl[lp,0] - rl[lm, 0]
+                dyl = rl[lp,1] - rl[lm, 1]
+                dzl = rl[lp,2] - rl[lm, 2]
+                dsl = np.sqrt(dxl**2 + dyl**2 + dzl**2)
+
+                if dsl != 0.0:
+                    dxl /= dsl
+                    dyl /= dsl
+                    dzl /= dsl
+
+                # Calculate perpendicular unit vectors
+                # Horizontal vector (in x-y plane, perpendicular to tangent)
+                uhx = -dyl
+                uhy = dxl
+                uhz = 0.0
+
+                # Vertical vector (perpendicular to both tangent and horizontal)
+                uvx = dyl * uhz - dzl * uhy
+                uvy = dzl * uhx - dxl * uhz
+                uvz = dxl * uhy - dyl * uhx
+
+                # Generate circle points
+                circle_x = []
+                circle_y = []
+                circle_z = []
+
+                for k in range(kk + 1):  # +1 to close the circle
+                    theta = 2.0 * np.pi * k / kk
+
+                    # Components of the circle point
+                    hx = uhx * np.cos(theta)
+                    hy = uhy * np.cos(theta)
+                    hz = uhz * np.cos(theta)
+
+                    vx = uvx * np.sin(theta)
+                    vy = uvy * np.sin(theta)
+                    vz = uvz * np.sin(theta)
+
+                    # Circle point at radius RADL
+                    circle_x.append(rl[l, 0] + (hx + vx) * radl[l])
+                    circle_y.append(rl[l, 1] + (hy + vy) * radl[l])
+                    circle_z.append(rl[l, 2] + (hz + vz) * radl[l])
+
+                # Plot the circle
+                circle_pts = {
+                    "x": np.array(circle_x),
+                    "y": np.array(circle_y),
+                    "z": np.array(circle_z),
+                }
+                axis.plot(circle_pts[xaxis], circle_pts[yaxis], "-", color=color, linewidth=0.5)
+
+    def plot_geom(self, axes=None, body_color="magenta"):
         """Generate a matplotlib plot of geometry
 
         Args:
             axes: Matplotlib axis object to add the plots too. If none are given, the axes will be generated.
+            body_color: Color to use for plotting body geometry (default: 'magenta')
         """
 
         if axes == None:
@@ -4213,8 +4360,10 @@ def plot_geom(self, axes=None):
             ax1, ax2 = axes
 
         self.add_mesh_plot(ax1, xaxis="y", yaxis="x")
+        self.add_body_plot(ax1, xaxis="y", yaxis="x", color=body_color)
 
         self.add_mesh_plot(ax2, xaxis="y", yaxis="z")
+        self.add_body_plot(ax2, xaxis="y", yaxis="z", color=body_color)
 
         if axes == None:
             # assume that if we don't provide axes that we want to see the plot

pyproject.toml

@@ -28,7 +28,7 @@ dependencies = [
     "numpy>=1.19",
 ]
 readme = "README.md"
-version = "2.1.1" # this automatically updates __init__.py and setup_deprecated.py
+version = "2.2.0" # this automatically updates __init__.py and setup_deprecated.py
 
 
 [tool.cibuildwheel]

src/aero.f

@@ -193,6 +193,12 @@ SUBROUTINE AERO
       CRBAX = DIR*CMTOT(1)  
       CMBAX =     CMTOT(2)
       CNBAX = DIR*CMTOT(3)
+
+      do IS = 1,NSURF
+        CMSURFBAX(1,IS) = DIR*CMSURF(1,IS) 
+        CMSURFBAX(2,IS) =     CMSURF(2,IS)
+        CMSURFBAX(3,IS) = DIR*CMSURF(3,IS)
+      enddo
 
 C
       ! compute the stability derivatives every time (it's quite cheap)
@@ -1137,6 +1143,12 @@ SUBROUTINE SFFORC
         CA_LSTRP(J) = CAXL0*COSAINC - CNRM0*SINAINC
         CN_LSTRP(J) = CNRM0*COSAINC + CAXL0*SINAINC
 C
+C------ vector at chord reference point from rotation axes
+        RROT(1) = XSREF(J) - XYZREF(1)
+        RROT(2) = YSREF(J) - XYZREF(2)
+        RROT(3) = ZSREF(J) - XYZREF(3)
+c        print *,"WROT ",WROT
+C
 C------ set total effective velocity = freestream + rotation
         CALL CROSS(RROT,WROT,VROT)
         VEFF(1) = VINF(1) + VROT(1)
@@ -1470,6 +1482,8 @@ SUBROUTINE BDFORC
       REAL CDBDY_U(NUMAX), CYBDY_U(NUMAX), CLBDY_U(NUMAX), 
      &     CFBDY_U(3,NUMAX), 
      &     CMBDY_U(3,NUMAX)
+      CHARACTER*50 SATYPE
+
 C
 C
       BETM = SQRT(1.0 - MACH**2)
@@ -1635,6 +1649,13 @@ SUBROUTINE BDFORC
         ENDDO
  200  CONTINUE
 C
+      ! compute the forces on the body in the body axis
+      CALL GETSA(LNASA_SA,SATYPE,DIR)
+
+      do IB = 1,NBODY
+            CMBDYBAX(3,IB) = DIR*CMBDY(3,IB)
+            CMBDYBAX(1,IB) = DIR*CMBDY(1,IB)
+      enddo
       RETURN
       END ! BDFORC
 

src/build/Makefile

@@ -31,7 +31,7 @@ f77FilesNoDir=$(notdir $(f77Files))
 # Generate two separate list of .F90 and .f90 files using the filter command
 
 # Finally convert all source files to .o
-OFILES=$(f90Files:%.f90=%.o) $(f77FilesNoDir:%.f=%.o) $(cFilesNoDir:%.c=%.o)
+OFILES=$(f90FilesNoDir:%.f90=%.o) $(f77FilesNoDir:%.f=%.o) $(cFilesNoDir:%.c=%.o)
 
 # Compile sources
 
@@ -81,10 +81,16 @@ default: lib ../f2py/libavl.pyf
 
 # build the module first
 avl_heap_inc.o avl_heap_inc.mod: ../avl_heap_inc.f90
-	$(FF90) $(FF90_ALL_FLAGS) -cpp -I ../includes -c $<
+	$(FF90) $(FF90_ALL_FLAGS) -cpp -I ../includes -c $<  -o $*.o
+	@echo
+	@echo "        --- Compiled $*.f90 successfully ---"
+	@echo
 
 avl_heap_diff_inc.o avl_heap_diff_inc.mod: ../avl_heap_diff_inc.f90
-	$(FF90) $(FF90_ALL_FLAGS) -cpp -I ../includes -c $<
+	$(FF90) $(FF90_ALL_FLAGS) -cpp -I ../includes -c  $< -o $*.o
+	@echo
+	@echo "        --- Compiled $*.f90 successfully ---"
+	@echo
 
 sources: $(OFILES)