Merge origin/main: make_dual_shadow builds every wing and keeps the reference point

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: 1-Bort-1

CHANGELOG.md

@@ -12,6 +12,8 @@
   doi:10.1088/1742-6596/2767/2/022068): each section gets a drag increment and a force
   along its span from the flow across it, in `solve!`, `solve` and `linearize`. Opt-in
   via `is_with_viscous_drag_correction` (default `false`) on the solver settings.
+- `linearize` takes a `BodyAerodynamics` with more than one wing; `theta_idxs` and
+  `delta_idxs` then run over the unrefined sections of all wings in order.
 
 ### Changed
 
@@ -26,6 +28,9 @@
   body z axis; it was read from the settings file and ignored.
 - The `VSMSolution` docstring gives `lift_dist`, `drag_dist` and `panel_moment_dist` in
   the per-unit-span units they hold, [N/m] and [Nm/m], instead of [N] and [Nm].
+- Inside its vortex core, `velocity_3D_trailing_vortex!` induces an azimuthal velocity
+  instead of a radial one. Only points within the millimetre-scale Oseen core of a
+  panel's chordwise trailing segment were affected.
 - With `artificial_damping` on, an iteration whose circulation is already smooth no longer
   re-applies the previous iteration's damping correction.
 

docs/src/private_functions.md

@@ -63,6 +63,7 @@ panel_loads
 ```@docs
 velocity_3D_bound_vortex!
 velocity_3D_trailing_vortex!
+velocity_3D_vortex_segment!
 velocity_3D_trailing_vortex_semiinfinite!
 calculate_velocity_induced_bound_2D!
 calculate_velocity_induced_single_ring_semiinfinite!
@@ -83,6 +84,7 @@ calculate_filaments_for_plotting
 ### Mesh refinement and billowing
 ```@docs
 unrefined_deform!
+unrefined_section_range
 deform!
 compute_refined_panel_mapping!
 compute_refined_section_interpolation!

src/body_aerodynamics.jl

@@ -256,6 +256,37 @@ function calculate_stall_angle_list!(stall_angles::AbstractVector,
     return nothing
 end
 
+"""
+    unrefined_section_range(body_aero::BodyAerodynamics, wing_idx)
+
+Indices of the unrefined sections of wing `wing_idx` in a distribution that runs over the
+unrefined sections of all wings in order, such as `moment_unrefined_dist`.
+"""
+function unrefined_section_range(body_aero::BodyAerodynamics, wing_idx)
+    offset = 0
+    for i in 1:wing_idx-1
+        offset += body_aero.wings[i].n_unrefined_sections
+    end
+    return offset .+ (1:body_aero.wings[wing_idx].n_unrefined_sections)
+end
+
+"""
+    unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles)
+
+Deform each wing of `body_aero` by its entries of `theta_angles` and `delta_angles` [rad],
+which run over the unrefined sections of all wings in order; `nothing` leaves that angle
+unchanged. Call [`reinit!`](@ref) afterwards to update the panels.
+"""
+function unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles)
+    for (wing_idx, wing) in enumerate(body_aero.wings)
+        section_range = unrefined_section_range(body_aero, wing_idx)
+        unrefined_deform!(wing,
+            isnothing(theta_angles) ? nothing : view(theta_angles, section_range),
+            isnothing(delta_angles) ? nothing : view(delta_angles, section_range))
+    end
+    return nothing
+end
+
 """
     reinit!(body_aero::BodyAerodynamics; init_aero, va, omega, refine_mesh, recompute_mapping, sort_sections)
 

src/filament.jl

@@ -39,10 +39,11 @@ function reinit!(filament::BoundFilament{T}, x1, x2, vec=zeros(MVector{3, T})) w
 end
 
 """
-    velocity_3D_bound_vortex(vel, filament::BoundFilament, XVP, 
-                           gamma, core_radius_fraction, work_vectors)
+    velocity_3D_bound_vortex!(vel, filament::BoundFilament, XVP,
+                              gamma, core_radius_fraction, work_vectors)
 
-Calculate induced velocity by a bound vortex filament at a point in space.
+Calculate induced velocity by a bound vortex filament at a point in space, with a
+core radius of `core_radius_fraction` times the filament length.
 """
 function velocity_3D_bound_vortex!(
     vel,
@@ -52,87 +53,24 @@ function velocity_3D_bound_vortex!(
     core_radius_fraction,
     work_vectors
 )
-    r1, r2, r1Xr2, r1Xr0, r2Xr0, r1r2norm, r1_proj, r2_proj,
-        r1_projXr2_proj, vel_ind_proj = work_vectors
-    r0 = filament.r0
-    nr0 = filament.length
-    r1 .= XVP .- filament.x1
-    r2 .= XVP .- filament.x2
-
-    epsilon = core_radius_fraction * nr0
-
-    cross3!(r1Xr0, r1, r0)
-
-    # Check point location relative to filament
-    nr1Xr0 = norm3(r1Xr0)
-    if nr1Xr0 / nr0 > epsilon
-        cross3!(r1Xr2, r1, r2)
-        nr1 = norm3(r1)
-        nr2 = norm3(r2)
-        @inbounds for k in 1:3
-            r1r2norm[k] = r1[k]/nr1 - r2[k]/nr2
-        end
-        nr1Xr2 = norm3(r1Xr2)
-        coeff = (gamma / (4π)) / (nr1Xr2^2) * dot3(r0, r1r2norm)
-        @inbounds for k in 1:3
-            vel[k] = coeff * r1Xr2[k]
-        end
-    elseif nr1Xr0 / nr0 < 1e-12 * epsilon
-        vel .= 0.0
-    else
-        @debug "inside core radius"
-        @debug "distance from control point to filament: $(nr1Xr0 / nr0)"
-
-        nr0sq = nr0 * nr0
-        d_r1_r0 = dot3(r1, r0)
-        d_r2_r0 = dot3(r2, r0)
-        r_rad = r1Xr0
-        @inbounds for k in 1:3
-            r_rad[k] = r1[k] - d_r1_r0 * r0[k] / nr0sq
-        end
-        nr_rad = norm3(r_rad)
-        @inbounds for k in 1:3
-            r1_proj[k] = d_r1_r0 * r0[k] / nr0sq +
-                         epsilon * r_rad[k] / nr_rad
-            r2_proj[k] = d_r2_r0 * r0[k] / nr0sq +
-                         epsilon * r_rad[k] / nr_rad
-        end
-        cross3!(r1_projXr2_proj, r1_proj, r2_proj)
-
-        nr1pXr2p = norm3(r1_projXr2_proj)
-        nr1_proj = norm3(r1_proj)
-        nr2_proj = norm3(r2_proj)
-        d_sum = 0.0
-        @inbounds for k in 1:3
-            d_sum += r0[k] * (r1_proj[k]/nr1_proj -
-                              r2_proj[k]/nr2_proj)
-        end
-        coeff = (gamma / (4π)) / (nr1pXr2p^2) * d_sum
-        @inbounds for k in 1:3
-            vel_ind_proj[k] = coeff * r1_projXr2_proj[k]
-        end
-
-        scale = nr1Xr0 / (nr0 * epsilon)
-        @inbounds for k in 1:3
-            vel[k] = scale * vel_ind_proj[k]
-        end
-    end
-    nothing
+    epsilon = core_radius_fraction * filament.length
+    velocity_3D_vortex_segment!(vel, filament, XVP, gamma, epsilon, work_vectors)
 end
 
 """
-    velocity_3D_trailing_vortex(vel, filament::BoundFilament, 
-                              XVP, gamma, va, work_vectors)
+    velocity_3D_trailing_vortex!(vel, filament::BoundFilament,
+                                 XVP, gamma, va, work_vectors)
 
-Calculate induced velocity by a trailing vortex filament.
+Calculate induced velocity by a trailing vortex filament, with a Lamb–Oseen core
+radius grown over the axial distance of `XVP` from the filament start.
 
 # Arguments
 - `XVP`: Control point coordinates
 - `gamma`: Vortex strength
 - `va`: Inflow velocity magnitude
 - work_vectors: preallocated array of intermediate variables
 
-Reference: Rick Damiani et al. "A vortex step method for nonlinear airfoil polar data 
+Reference: Rick Damiani et al. "A vortex step method for nonlinear airfoil polar data
 as implemented in KiteAeroDyn".
 """
 @inline function velocity_3D_trailing_vortex!(
@@ -143,71 +81,78 @@ as implemented in KiteAeroDyn".
     va,
     work_vectors
 )
-    r1 = work_vectors[2]
-    r2 = work_vectors[3]
-    r_perp = work_vectors[4]
-    r1Xr2 = work_vectors[5]
-    r1Xr0 = work_vectors[6]
-    r2Xr0 = work_vectors[7]
-    normr1r2 = work_vectors[8]
+    r1 = work_vectors[1]
+    r1 .= XVP .- filament.x1
+    axial_distance = abs(dot3(r1, filament.r0)) / filament.length
+    epsilon = sqrt(4 * ALPHA0 * NU * axial_distance / va)
+    velocity_3D_vortex_segment!(vel, filament, XVP, gamma, epsilon, work_vectors)
+end
 
+"""
+    velocity_3D_vortex_segment!(vel, filament::BoundFilament, XVP,
+                                gamma, epsilon, work_vectors)
+
+Calculate the Biot–Savart velocity induced by a straight vortex segment at `XVP`.
+Inside the core radius `epsilon` the velocity is evaluated on the core boundary and
+scaled linearly with the distance to the axis.
+"""
+@inline function velocity_3D_vortex_segment!(
+    vel,
+    filament::BoundFilament,
+    XVP,
+    gamma,
+    epsilon,
+    work_vectors
+)
+    r1, r2, r1Xr2, r1Xr0, r1r2norm, r1_proj, r2_proj = work_vectors
     r0 = filament.r0
     nr0 = filament.length
     r1 .= XVP .- filament.x1
     r2 .= XVP .- filament.x2
 
-    nr0sq = nr0 * nr0
-    d_r1_r0 = dot3(r1, r0)
-
-    # Cut-off radius. The perpendicular component has length |r1.r0|/|r0|, so the
-    # vector itself is only needed inside the core.
-    epsilon = sqrt(4 * ALPHA0 * NU * abs(d_r1_r0) / nr0 / va)
-
     cross3!(r1Xr0, r1, r0)
-
-    # Check point location relative to filament
     nr1Xr0 = norm3(r1Xr0)
     if nr1Xr0 / nr0 > epsilon
         cross3!(r1Xr2, r1, r2)
         nr1 = norm3(r1)
         nr2 = norm3(r2)
         @inbounds for k in 1:3
-            normr1r2[k] = r1[k]/nr1 - r2[k]/nr2
+            r1r2norm[k] = r1[k]/nr1 - r2[k]/nr2
         end
         nr1Xr2 = norm3(r1Xr2)
-        coeff = (gamma / (4π)) / (nr1Xr2^2) * dot3(r0, normr1r2)
+        coeff = (gamma / (4π)) / (nr1Xr2^2) * dot3(r0, r1r2norm)
         @inbounds for k in 1:3
             vel[k] = coeff * r1Xr2[k]
         end
     elseif nr1Xr0 / nr0 < 1e-12 * epsilon
         vel .= 0.0
     else
-        # Project onto core radius — reuse r_perp, normr1r2
-        r1_proj = r_perp
-        r2_proj = normr1r2
-        cross3!(r2Xr0, r2, r0)
-        nr2Xr0 = norm3(r2Xr0)
+        nr0sq = nr0 * nr0
+        d_r1_r0 = dot3(r1, r0)
         d_r2_r0 = dot3(r2, r0)
+        r_rad = r1Xr0
+        @inbounds for k in 1:3
+            r_rad[k] = r1[k] - d_r1_r0 * r0[k] / nr0sq
+        end
+        nr_rad = norm3(r_rad)
         @inbounds for k in 1:3
             r1_proj[k] = d_r1_r0 * r0[k] / nr0sq +
-                         epsilon * r1Xr0[k] / nr1Xr0
+                         epsilon * r_rad[k] / nr_rad
             r2_proj[k] = d_r2_r0 * r0[k] / nr0sq +
-                         epsilon * r2Xr0[k] / nr2Xr0
+                         epsilon * r_rad[k] / nr_rad
         end
-
         cross3!(r1Xr2, r1_proj, r2_proj)
-        nr1Xr2_val = norm3(r1Xr2)
         nr1_proj = norm3(r1_proj)
         nr2_proj = norm3(r2_proj)
         d_sum = 0.0
         @inbounds for k in 1:3
             d_sum += r0[k] * (r1_proj[k]/nr1_proj -
                               r2_proj[k]/nr2_proj)
         end
-        coeff = (gamma / (4π)) / (nr1Xr2_val^2) * d_sum
         scale = nr1Xr0 / (nr0 * epsilon)
+        coeff = scale * (gamma / (4π)) / (norm3(r1Xr2)^2) * d_sum
         @inbounds for k in 1:3
-            vel[k] = scale * coeff * r1Xr2[k]
+            vel[k] = coeff * r1Xr2[k]
         end
     end
     nothing
@@ -263,8 +208,7 @@ function velocity_3D_trailing_vortex_semiinfinite!(
     GAMMA = -GAMMA * filament.filament_direction
     r1 .= XVP .- filament.x1
 
-    # Core radius. `r_perp` is `(r1.Vf) Vf`, so its length is `|r1.Vf| |Vf|` and
-    # the vector itself is only needed inside the core.
+    # Core radius, grown with the axial distance of `XVP` along `Vf`.
     d_r1_Vf = dot3(r1, Vf)
     nVf = norm3(Vf)
     epsilon = sqrt(4 * ALPHA0 * NU * abs(d_r1_Vf) * nVf / va)