Update Visual Studio projects in dev-tc

.github/scripts/build_windows_executables.bat

@@ -41,15 +41,15 @@ echo on
 
 
 @REM Build MATLAB solution last
-echo on
-devenv vs-build/OpenFAST.sln /Build "Release_Matlab|x64"
-echo off
-if %ERRORLEVEL% NEQ 0 (
-    set "FailedSolutions=!FailedSolutions!Release_Matlab  "
-    set "OverallErrorLevel=1"
-    echo Build of OpenFAST.sln Release_Matlab failed!
-)
-echo on
+@REM echo on
+@REM devenv vs-build/OpenFAST.sln /Build "Release_Matlab|x64"
+@REM echo off
+@REM if %ERRORLEVEL% NEQ 0 (
+@REM     set "FailedSolutions=!FailedSolutions!Release_Matlab  "
+@REM     set "OverallErrorLevel=1"
+@REM     echo Build of OpenFAST.sln Release_Matlab failed!
+@REM )
+@REM echo on
 
 
 

.github/workflows/deploy.yml

@@ -163,39 +163,39 @@ jobs:
           rd /s/q "webimage_extracted"
           exit /b %installer_exit_code%
 
-      - name: Set up MATLAB
-        id: setup-matlab
-        uses: matlab-actions/setup-matlab@v2
-        with:
-          products: Simulink
+      # - name: Set up MATLAB
+      #   id: setup-matlab
+      #   uses: matlab-actions/setup-matlab@v2
+      #   with:
+      #     products: Simulink
 
       - name: Build Executables
         env:
           MATLAB_ROOT: ${{ steps.setup-matlab.outputs.matlabroot }}
         run: .github/scripts/build_windows_executables.bat
 
-      - name: Test version output
-        run: .github/scripts/test_windows_executables.bat
+      # - name: Test version output
+      #   run: .github/scripts/test_windows_executables.bat
 
-      - name: Build MATLAB Mex File
-        uses: matlab-actions/run-command@v2
-        with:
-          command: |
-            mexname = 'FAST_SFunc';
-            mex('-largeArrayDims', ...
-                '-v', ...
-                ['-L' fullfile('build','bin')], ...
-                ['-I' fullfile('modules','openfast-library','src')], ... % "FAST_Library.h"
-                ['-I' fullfile('modules','supercontroller','src')], ...  % "SuperController_Types.h"
-                ['-I' fullfile('modules','externalinflow','src')],  ...  % "ExternalInflow_Types.h"
-                ['-I' fullfile('modules','extloads','src')], ...         % "ExtLoadsDX_Types.h"
-                ['-I' fullfile(matlabroot,'simulink','include')], ...
-                ['-I' fullfile(matlabroot,'extern','include')], ...
-                ['COMPFLAGS=$COMPFLAGS -MT -DS_FUNCTION_NAME=' mexname], ...
-                '-lOpenFAST-Simulink_x64', ...
-                '-outdir', fullfile('build','bin'), ...
-                '-output', mexname, ...
-                fullfile('glue-codes','simulink','src','FAST_SFunc.c'));
+      # - name: Build MATLAB Mex File
+      #   uses: matlab-actions/run-command@v2
+      #   with:
+      #     command: |
+      #       mexname = 'FAST_SFunc';
+      #       mex('-largeArrayDims', ...
+      #           '-v', ...
+      #           ['-L' fullfile('build','bin')], ...
+      #           ['-I' fullfile('modules','openfast-library','src')], ... % "FAST_Library.h"
+      #           ['-I' fullfile('modules','supercontroller','src')], ...  % "SuperController_Types.h"
+      #           ['-I' fullfile('modules','externalinflow','src')],  ...  % "ExternalInflow_Types.h"
+      #           ['-I' fullfile('modules','extloads','src')], ...         % "ExtLoadsDX_Types.h"
+      #           ['-I' fullfile(matlabroot,'simulink','include')], ...
+      #           ['-I' fullfile(matlabroot,'extern','include')], ...
+      #           ['COMPFLAGS=$COMPFLAGS -MT -DS_FUNCTION_NAME=' mexname], ...
+      #           '-lOpenFAST-Simulink_x64', ...
+      #           '-outdir', fullfile('build','bin'), ...
+      #           '-output', mexname, ...
+      #           fullfile('glue-codes','simulink','src','FAST_SFunc.c'));
 
       - name: Upload executables
         uses: actions/upload-artifact@v4

modules/aerodyn/src/FVW_BiotSavart.f90

@@ -50,55 +50,61 @@ subroutine ui_seg_11(DeltaPa, DeltaPb, SegGamma, RegFunction, RegParam1, Uind)
    real(ReKi)              :: xa, ya, za, xb, yb, zb !< Coordinates of X-Xa and X-Xb
    real(ReKi)              :: exp_value      !< 
    ! 
-   Uind(1:3)=0.0_ReKi
+   Uind(1:3) = 0.0_ReKi
    xa=DeltaPa(1); ya=DeltaPa(2); za=DeltaPa(3)
    xb=DeltaPb(1); yb=DeltaPb(2); zb=DeltaPb(3)
    norm_a      = sqrt(xa*xa + ya*ya + za*za)
    norm_b      = sqrt(xb*xb + yb*yb + zb*zb)
    denominator  = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb) ! |r1|*|r2|*(|r1|*|r2| + r1.r2)
-   if (denominator>PRECISION_UI) then
-      crossprod(1) = ya*zb-za*yb; crossprod(2) = za*xb-xa*zb; crossprod(3) = xa*yb-ya*xb
-      norm2_orth   = crossprod(1)**2 + crossprod(2)**2 + crossprod(3)**2
-      if (norm2_orth>PRECISION_UI) then ! On the singularity, Uind(1:3)=0.0_ReKi
-         norm2_r0     = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) +(za-zb)*(za-zb) 
-         if (norm2_r0>PRECISION_UI) then ! segment of zero length
-            ! --- Far field TODO
-            ! --- Regularization (close field)
-            norm2_orth   = norm2_orth/norm2_r0 ! d = (r1xr2)/r0
-            select case (RegFunction) !
-            case ( idRegNone )         ! No vortex core model
-               Kv=1.0_ReKi
-            case ( idRegRankine )      ! Rankine
-               r_bar2    = norm2_orth/ RegParam1**2
-               if (r_bar2<1) then 
-                  Kv=r_bar2
-               else
-                  Kv=1.0_ReKi 
-               end if 
-            case ( idRegLambOseen )    ! Lamb-Oseen
-               r_bar2    = norm2_orth/ RegParam1**2
-               exp_value = -1.25643_ReKi*r_bar2
-               if(exp_value<MIN_EXP_VALUE) then ! Remove me when Far distance implemented
-                  Kv = 1.0_ReKi
-               else
-                  Kv = 1.0_ReKi-exp(exp_value)
-               endif
-            case ( idRegVatistas )    ! Vatistas n=2
-               r_bar2 = norm2_orth/ RegParam1**2
-               Kv     = r_bar2/sqrt(1.0_ReKi+r_bar2**2)
-            case ( idRegOffset )      ! Cut-off radius 
-               Kv        = 1.0_ReKi
-               denominator=denominator+RegParam1**2*norm2_r0
-            case default
-               print*,'Unknown SgmtReg', RegFunction
-               STOP ! Will never happen
-               Kv=1.0_ReKi !< Should be an error
-            end select 
-            Kv=SegGamma*fourpi_inv*Kv*(norm_a+norm_b)/denominator
-            Uind(1:3) = Kv*crossprod(1:3)
-         endif
+
+   if (denominator <= PRECISION_UI) return
+   crossprod(1) = ya*zb-za*yb; crossprod(2) = za*xb-xa*zb; crossprod(3) = xa*yb-ya*xb
+   norm2_orth = crossprod(1)**2 + crossprod(2)**2 + crossprod(3)**2
+
+   ! On the singularity, Uind(1:3)=0.0_ReKi
+   if (norm2_orth <= PRECISION_UI) return
+   norm2_r0 = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) + (za-zb)*(za-zb) 
+
+   ! segment of zero length
+   if (norm2_r0 <= PRECISION_UI) return 
+
+   ! --- Far field TODO
+   ! --- Regularization (close field)
+   norm2_orth = norm2_orth/norm2_r0 ! d = (r1xr2)/r0
+
+   select case (RegFunction) !
+   case ( idRegNone )         ! No vortex core model
+      Kv=1.0_ReKi
+   case ( idRegRankine )      ! Rankine
+      r_bar2    = norm2_orth/ RegParam1**2
+      if (r_bar2<1) then 
+         Kv=r_bar2
+      else
+         Kv=1.0_ReKi 
+      end if 
+   case ( idRegLambOseen )    ! Lamb-Oseen
+      r_bar2    = norm2_orth/ RegParam1**2
+      exp_value = -1.25643_ReKi*r_bar2
+      if(exp_value<MIN_EXP_VALUE) then ! Remove me when Far distance implemented
+         Kv = 1.0_ReKi
+      else
+         Kv = 1.0_ReKi-exp(exp_value)
       endif
-   endif
+   case ( idRegVatistas )    ! Vatistas n=2
+      r_bar2 = norm2_orth/ RegParam1**2
+      Kv     = r_bar2/sqrt(1.0_ReKi+r_bar2**2)
+   case ( idRegOffset )      ! Cut-off radius 
+      Kv        = 1.0_ReKi
+      denominator=denominator+RegParam1**2*norm2_r0
+   case default
+      print*,'Unknown SgmtReg', RegFunction
+      STOP ! Will never happen
+      Kv=1.0_ReKi !< Should be an error
+   end select 
+
+   Kv = SegGamma*fourpi_inv*Kv*(norm_a+norm_b)/denominator
+   Uind(1:3) = Kv*crossprod(1:3)
+
 end subroutine ui_seg_11
 
 
@@ -134,21 +140,23 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, &
    real(ReKi)              :: norm2_r0        !< Squared length of the segment d = (r1xr2)/r0
    real(ReKi)              :: xa, ya, za, xb, yb, zb !< Coordinates of X-Xa and X-Xb
    real(ReKi)              :: exp_value       !< 
+   real(ReKi)              :: CPs_icp(3)      !< 
 
    ! Branching based on regularization model
    ! NOTE: copy paste of code is done for optimization!
    !       The only thing changing is the part labelled "regularization"
    select case (RegFunction) 
    case ( idRegNone ) ! No vortex core 
       !$OMP PARALLEL default(shared)
-      !$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
+      !$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
       do icp=iCPStart,iCPEnd ! loop on CPs 
+         Uind = 0.0_ReKi
+         CPs_icp = CPs(:,icp)
          do is=iSegStart,iSegEnd ! loop on selected segments 
-            Uind = 0.0_ReKi
             P1   = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points
             P2   = SegPoints(1:3, SegConnct(2,is))
-            xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3);
-            xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3);
+            xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3);
+            xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3);
             norm_a      = sqrt(xa*xa + ya*ya + za*za)
             norm_b      = sqrt(xb*xb + yb*yb + zb*zb)
             denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
@@ -162,26 +170,27 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, &
                      ! --- Far field TODO
                      ! --- NO Regularization (close field)
                      Kv        = SegGamma(is)*fourpi_inv*(norm_a+norm_b)/(denominator + MINDENOM)
-                     Uind(1:3) = Kv*crossprod(1:3)
+                     Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3)
                   end if
                end if
             end if
-            Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
          end do ! Loop on segments
+         Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
       enddo ! Loop on control points
       !$OMP END DO 
       !$OMP END PARALLEL
 
    case ( idRegRankine )      ! Rankine
       !$OMP PARALLEL default(shared)
-      !$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
+      !$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
       do icp=iCPStart,iCPEnd ! loop on CPs 
+         Uind = 0.0_ReKi
+         CPs_icp = CPs(:,icp)
          do is=iSegStart,iSegEnd ! loop on selected segments 
-            Uind = 0.0_ReKi
             P1   = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points
             P2   = SegPoints(1:3, SegConnct(2,is))
-            xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3);
-            xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3);
+            xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3);
+            xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3);
             norm_a      = sqrt(xa*xa + ya*ya + za*za)
             norm_b      = sqrt(xb*xb + yb*yb + zb*zb)
             denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
@@ -201,26 +210,27 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, &
                         Kv=1.0_ReKi 
                      end if 
                      Kv        = SegGamma(is)*fourpi_inv*Kv*(norm_a+norm_b)/(denominator + MINDENOM)
-                     Uind(1:3) = Kv*crossprod(1:3)
+                     Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3)
                   end if
                end if ! denominator size or distances too small
             end if ! 
-            Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
          end do ! Loop on segments
+         Uind_out(1:3,icp) = Uind_out(1:3,icp) + Uind(1:3)
       enddo ! Loop on control points
       !$OMP END DO 
       !$OMP END PARALLEL
 
    case ( idRegLambOseen )      ! LambOseen
       !$OMP PARALLEL default(shared)
-      !$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb,exp_value) schedule(runtime)
+      !$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb,exp_value) schedule(runtime)
       do icp=iCPStart,iCPEnd ! loop on CPs 
+         Uind = 0.0_ReKi
+         CPs_icp = CPs(:,icp)
          do is=iSegStart,iSegEnd ! loop on selected segments 
-            Uind = 0.0_ReKi
             P1   = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points
             P2   = SegPoints(1:3, SegConnct(2,is))
-            xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3);
-            xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3);
+            xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3);
+            xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3);
             norm_a      = sqrt(xa*xa + ya*ya + za*za)
             norm_b      = sqrt(xb*xb + yb*yb + zb*zb)
             denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
@@ -241,61 +251,63 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, &
                         Kv = 1.0_ReKi-exp(exp_value)
                      endif
                      Kv        = SegGamma(is)*fourpi_inv*Kv*(norm_a+norm_b)/(denominator + MINDENOM)
-                     Uind(1:3) = Kv*crossprod(1:3)
+                     Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3)
                   endif 
                end if ! denominator size or distances too small
             end if ! 
-            Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
          end do ! Loop on segments
+         Uind_out(1:3,icp) = Uind_out(1:3,icp) + Uind(1:3)
       enddo ! Loop on control points
       !$OMP END DO 
       !$OMP END PARALLEL
 
    case ( idRegVatistas )      ! Vatistas n=2
       !$OMP PARALLEL default(shared)
-      !$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
+      !$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
       do icp=iCPStart,iCPEnd ! loop on CPs 
+         Uind = 0.0_ReKi
+         CPs_icp = CPs(:,icp)
          do is=iSegStart,iSegEnd ! loop on selected segments 
-            Uind = 0.0_ReKi
             P1   = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points
             P2   = SegPoints(1:3, SegConnct(2,is))
-            xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3);
-            xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3);
+            xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3);
+            xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3);
             norm_a      = sqrt(xa*xa + ya*ya + za*za)
             norm_b      = sqrt(xb*xb + yb*yb + zb*zb)
             denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
-            if (denominator>PRECISION_UI) then
-               crossprod(1) = ya*zb-za*yb; crossprod(2) = za*xb-xa*zb; crossprod(3) = xa*yb-ya*xb
-               norm2_orth   = crossprod(1)**2 + crossprod(2)**2 + crossprod(3)**2
-               if (norm2_orth>PRECISION_UI) then ! On the singularity, Uind(1:3)=0.0_ReKi
-                  norm2_r0     = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) +(za-zb)*(za-zb) 
-                  if (norm2_r0>PRECISION_UI) then
-                     ! --- Far field TODO
-                     ! --- Regularization (close field) --- Vatistas
-                     norm2_orth = norm2_orth/norm2_r0 ! d = (r1xr2)/r0
-                     r_bar2     = norm2_orth/RegParam(is)**2
-                     Kv         = r_bar2/sqrt(1+r_bar2**2)
-                     Kv         = SegGamma(is)*fourpi_inv*Kv*(norm_a+norm_b)/(denominator + MINDENOM)
-                     Uind(1:3)  = Kv*crossprod(1:3)
-                  end if
-               end if ! denominator size or distances too small
-            end if ! 
-            Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
+            ! denominator size or distances too small
+            if (denominator <= PRECISION_UI) cycle
+            crossprod(1) = ya*zb-za*yb; crossprod(2) = za*xb-xa*zb; crossprod(3) = xa*yb-ya*xb
+            norm2_orth   = crossprod(1)**2 + crossprod(2)**2 + crossprod(3)**2
+            ! On the singularity, cycle
+            if (norm2_orth <= PRECISION_UI) cycle
+            norm2_r0     = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) +(za-zb)*(za-zb) 
+            ! segment of zero length
+            if (norm2_r0 <= PRECISION_UI) cycle
+            ! --- Far field TODO
+            ! --- Regularization (close field) --- Vatistas
+            norm2_orth = norm2_orth/norm2_r0 ! d = (r1xr2)/r0
+            r_bar2     = norm2_orth/RegParam(is)**2
+            Kv         = r_bar2/sqrt(1.0_ReKi+r_bar2**2)
+            Kv         = SegGamma(is)*fourpi_inv*Kv*(norm_a+norm_b)/(denominator + MINDENOM)
+            Uind(1:3)  = Uind(1:3) + Kv*crossprod(1:3)
          end do ! Loop on segments
+         Uind_out(1:3,icp) = Uind_out(1:3,icp) + Uind(1:3)
       enddo ! Loop on control points
       !$OMP END DO 
       !$OMP END PARALLEL
 
    case ( idRegOffset )      ! Denominator offset
       !$OMP PARALLEL default(shared)
-      !$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
+      !$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
       do icp=iCPStart,iCPEnd ! loop on CPs 
+         Uind      = 0.0_ReKi
+         CPs_icp = CPs(:,icp)
          do is=iSegStart,iSegEnd ! loop on selected segments 
-            Uind      = 0.0_ReKi
             P1        = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points
             P2        = SegPoints(1:3, SegConnct(2,is))
-            xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3);
-            xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3);
+            xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3);
+            xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3);
             norm_a      = sqrt(xa*xa + ya*ya + za*za)
             norm_b      = sqrt(xb*xb + yb*yb + zb*zb)
             denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
@@ -309,12 +321,12 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, &
                      ! --- Regularization (close field) -- Offset
                      denominator = denominator+RegParam(is)**2*norm2_r0
                      Kv          = SegGamma(is)*fourpi_inv*(norm_a+norm_b)/(denominator + MINDENOM)
-                     Uind(1:3)   = Kv*crossprod(1:3)
+                     Uind(1:3)   = Uind(1:3) + Kv*crossprod(1:3)
                   end if
                end if ! denominator size or distances too small
             end if ! 
-            Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
          end do ! Loop on segments
+         Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
       enddo ! Loop on control points
       !$OMP END DO 
       !$OMP END PARALLEL