Add full directional wave examples and test

Full_Directional_Waves/MHKiT_NDBC_Waves/userDefinedFunctions.m

@@ -0,0 +1,22 @@
+%Example of user input MATLAB file for post processing
+
+%Plot waves
+waves.plotElevation(simu.rampTime);
+try 
+    waves.plotSpectrum();
+catch
+end
+
+% Plot RY forces for body 1
+output.plotForces(1,5)
+
+%Plot RY response for body 1
+output.plotResponse(1,5);
+
+% Plot x forces for body 2
+output.plotForces(2,1)
+
+% Save waves and response as video
+% output.saveViz(simu,body,waves,...
+%     'timesPerFrame',5,'axisLimits',[-50 50 -50 50 -12 20],...
+%     'startEndTime',[100 125]);

Full_Directional_Waves/MHKiT_NDBC_Waves/wecSimInputFile.m

@@ -0,0 +1,41 @@
+%% Simulation Data
+simu = simulationClass();                    % Initialize Simulation Class
+simu.simMechanicsFile = 'OSWEC.slx';         % Specify Simulink Model File
+simu.mode      = 'normal';                   % Specify Simulation Mode ('normal','accelerator','rapid-accelerator')
+simu.explorer  = 'on';                       % Turn SimMechanics Explorer (on/off)
+simu.startTime = 0;                          % Simulation Start Time [s]
+simu.endTime   = 400;                        % Simulation End Time [s]
+simu.rampTime  = 100;                        % Wave Ramp Time [s]
+simu.solver    = 'ode4';                     % simu.solver = 'ode4' for fixed step & simu.solver = 'ode45' for variable step 
+simu.dt        = 0.05; 					  	 % Simulation time-step [s]
+simu.cicEndTime = 30;
+
+%% Wave Information 
+% Full directional waves
+waves = waveClass('spectrumImportFullDir');
+waves.spectrumFile = ('fullDirSpectrum.mat');
+waves.phaseSeed = 1;
+
+%% Body Data
+% Flap
+body(1) = bodyClass('../../_Common_Input_Files/OSWEC/hydroData/oswec.h5');      % Initialize bodyClass for Flap
+body(1).geometryFile = '../../_Common_Input_Files/OSWEC/geometry/flap.stl';     % Geometry File
+body(1).mass = 127000;                          % User-Defined mass [kg]
+body(1).inertia = [1.85e6 1.85e6 1.85e6];       % Moment of Inertia [kg-m^2]
+
+% Base
+body(2) = bodyClass('../../_Common_Input_Files/OSWEC/hydroData/oswec.h5');      % Initialize bodyClass for Base
+body(2).geometryFile = '../../_Common_Input_Files/OSWEC/geometry/base.stl';     % Geometry File
+body(2).mass = 999;                             % Placeholder mass for a fixed body
+body(2).inertia = [999 999 999];                % Placeholder inertia for a fixed body
+
+%% PTO and Constraint Parameters
+% Fixed
+constraint(1)= constraintClass('Constraint1');  % Initialize ConstraintClass for Constraint1
+constraint(1).location = [0 0 -10];             % Constraint Location [m]
+
+% Rotational PTO
+pto(1) = ptoClass('PTO1');                      % Initialize ptoClass for PTO1
+pto(1).stiffness = 0;                           % PTO Stiffness Coeff [Nm/rad]
+pto(1).damping = 12000;                         % PTO Damping Coeff [Nsm/rad]
+pto(1).location = [0 0 -8.9];                   % PTO Location [m]

Full_Directional_Waves/OOI_Waves/convertOOIToIEC.m

@@ -0,0 +1,51 @@
+function [frequenciesIEC, spectrumIEC, spreadIEC, directionsIEC] = convertOOIToIEC(spectrumDataOOI, directionBins, plotFlag)
+
+    directionBins = wrapTo180(directionBins);
+    if any(directionBins == 180) && any(directionBins == -180)
+        warning('Directions include both -180 and 180. Removing -180 to avoid duplicate directions.')
+        directionBins(directionBins == -180) = [];
+    end
+    if length(unique(mod(directionBins + 180, 360) - 180)) < length(directionBins)
+        error('Duplicate directions were found.')
+    end
+
+    frequenciesOOI = spectrumDataOOI(:,1);
+    spectrumOOI = spectrumDataOOI(:,2);
+    directionsOOI = spectrumDataOOI(:,3);
+    spreadOOI = spectrumDataOOI(:,4);
+
+    % convert to IEC: 
+    directionsIEC = directionBins(:);
+    frequenciesIEC = frequenciesOOI;
+    spectrumIEC = spectrumOOI;
+    spreadIEC = zeros(length(frequenciesIEC),length(directionsIEC));
+
+    % for spread, we need to apply gaussian distribution
+    for ii = 1:length(frequenciesOOI)
+        % convert to IEC spread calculation
+        directions = deg2rad(wrapTo180(directionBins - wrapTo180(directionsOOI(ii))));
+        energySpread = (1./(deg2rad(spreadOOI(ii)).*sqrt(2*pi))) .* exp(-(directions.^2) ./ (2.*deg2rad(spreadOOI(ii)).^2));
+        checkSum = trapz(deg2rad(directionBins),energySpread);
+        if checkSum < 0.95 % if this is true, then less than 95% of the initial energy at this frequency is contained over the considered directions.
+            warning('Number of directional bins inadequate at frequency number %i. Directional approximation weak. \n \r', ii);
+        end
+        energySpread = energySpread./checkSum;
+        spreadIEC(ii,:) = energySpread;
+    end
+
+    if plotFlag == 1
+        spectrumFullDir = spreadIEC.*spectrumIEC;
+        meanDirection = sum(wrapTo180(directionsOOI).*spectrumOOI)/sum(spectrumOOI);
+
+        [plotDirs,plotFreqs] = meshgrid(directionsIEC,frequenciesIEC);
+
+        figure()
+        polarscatter(deg2rad(plotDirs(:)),plotFreqs(:),4,spectrumFullDir(:),'filled')
+        hold on
+        polarplot(deg2rad([meanDirection meanDirection]), [0 max(plotFreqs(:))], 'k--'); % 
+        c = colorbar;
+        c.Label.String = 'Spectrum (m^2/Hz/deg)';
+        title('Elevation variance spectrum');
+        legend('spectrum','approx. mean direction','Location','northwest')
+    end
+end

Full_Directional_Waves/OOI_Waves/userDefinedFunctions.m

@@ -0,0 +1,22 @@
+%Example of user input MATLAB file for post processing
+
+%Plot waves
+waves.plotElevation(simu.rampTime);
+try 
+    waves.plotSpectrum();
+catch
+end
+
+% Plot RY forces for body 1
+output.plotForces(1,5)
+
+%Plot RY response for body 1
+output.plotResponse(1,5);
+
+% Plot x forces for body 2
+output.plotForces(2,1)
+
+% Save waves and response as video
+% output.saveViz(simu,body,waves,...
+%     'timesPerFrame',5,'axisLimits',[-50 50 -50 50 -12 20],...
+%     'startEndTime',[100 125]);

Full_Directional_Waves/OOI_Waves/wecSimInputFile.m

@@ -0,0 +1,48 @@
+%% Simulation Data
+simu = simulationClass();                    % Initialize Simulation Class
+simu.simMechanicsFile = 'OSWEC.slx';         % Specify Simulink Model File
+simu.mode      = 'normal';                   % Specify Simulation Mode ('normal','accelerator','rapid-accelerator')
+simu.explorer  = 'on';                       % Turn SimMechanics Explorer (on/off)
+simu.startTime = 0;                          % Simulation Start Time [s]
+simu.endTime   = 400;                        % Simulation End Time [s]
+simu.rampTime  = 100;                        % Wave Ramp Time [s]
+simu.solver    = 'ode4';                     % simu.solver = 'ode4' for fixed step & simu.solver = 'ode45' for variable step 
+simu.dt        = 0.05; 					  	 % Simulation time-step [s]
+simu.cicEndTime = 30;
+
+%% Wave Information 
+load dirSpectrumOOI.mat
+spectrumDataOOI = dataWaveSnip;
+directions = -179:2:179;
+
+[frequencies, spectrum, spread, directions] = convertOOIToIEC(spectrumDataOOI, directions, 1);
+
+save 'fullDirSpectrum.mat' spectrum spread frequencies directions
+
+waves = waveClass('spectrumImportFullDir');
+waves.spectrumFile = ('fullDirSpectrum.mat');
+waves.phaseSeed = 1;
+
+%% Body Data
+% Flap
+body(1) = bodyClass('../../_Common_Input_Files/OSWEC/hydroData/oswec.h5');      % Initialize bodyClass for Flap
+body(1).geometryFile = '../../_Common_Input_Files/OSWEC/geometry/flap.stl';     % Geometry File
+body(1).mass = 127000;                          % User-Defined mass [kg]
+body(1).inertia = [1.85e6 1.85e6 1.85e6];       % Moment of Inertia [kg-m^2]
+
+% Base
+body(2) = bodyClass('../../_Common_Input_Files/OSWEC/hydroData/oswec.h5');      % Initialize bodyClass for Base
+body(2).geometryFile = '../../_Common_Input_Files/OSWEC/geometry/base.stl';     % Geometry File
+body(2).mass = 999;                             % Placeholder mass for a fixed body
+body(2).inertia = [999 999 999];                % Placeholder inertia for a fixed body
+
+%% PTO and Constraint Parameters
+% Fixed
+constraint(1)= constraintClass('Constraint1');  % Initialize ConstraintClass for Constraint1
+constraint(1).location = [0 0 -10];             % Constraint Location [m]
+
+% Rotational PTO
+pto(1) = ptoClass('PTO1');                      % Initialize ptoClass for PTO1
+pto(1).stiffness = 0;                           % PTO Stiffness Coeff [Nm/rad]
+pto(1).damping = 12000;                         % PTO Damping Coeff [Nsm/rad]
+pto(1).location = [0 0 -8.9];                   % PTO Location [m]

Full_Directional_Waves/TestFullDirectionalWaves.m

@@ -0,0 +1,64 @@
+classdef TestFullDirectionalWaves < matlab.unittest.TestCase
+
+    properties
+        OriginalDefault
+        testDir
+        h5Dir = '../_Common_Input_Files/OSWEC/hydroData/'
+        h5Name = 'oswec.h5'
+        outName = 'oswec.out'
+    end    
+
+    methods (Access = 'public')        
+        function obj = TestFullDirectionalWaves
+            obj.testDir = fileparts(mfilename('fullpath'));
+        end    
+    end
+
+    methods (TestMethodSetup)
+        function killPlots (~)
+            set(0,'DefaultFigureVisible','off');
+        end
+    end
+
+    methods(TestClassSetup)        
+        function captureVisibility(testCase)
+            testCase.OriginalDefault = get(0,'DefaultFigureVisible');
+        end        
+        function runBemio(testCase)            
+            cd(testCase.h5Dir);
+            if isfile(testCase.h5Name)
+                fprintf('runBemio skipped, *.h5 already exists\n')
+            else
+                bemio
+            end
+            cd(testCase.testDir)            
+        end        
+    end
+
+    methods(TestMethodTeardown)
+        function returnHome(testCase)
+            cd(testCase.testDir)
+        end
+    end
+
+    methods(TestClassTeardown)        
+        function checkVisibilityRestored(testCase)
+            set(0,'DefaultFigureVisible',testCase.OriginalDefault);
+            testCase.assertEqual(get(0,'DefaultFigureVisible'),     ...
+                                 testCase.OriginalDefault);
+        end
+    end
+
+    methods(Test)        
+        function testMHKiT_NDBC_Waves(testCase)
+            cd MHKiT_NDBC_Waves
+            wecSim
+            close_system('OSWEC',0)
+        end        
+        function testOOI_Waves(testCase)
+            cd OOI_Waves
+            wecSim
+            close_system('OSWEC',0)
+        end        
+    end    
+end