% Класс задачи.
% Copyright 2017 Daniil S. Denisov
classdef StructFEProblem < handle
    properties (Access = public)
        % Имя кейс-файла.
        filename
        % Матрица жесткости.
        K
        % Матрица масс.
        M
        % Вектор правой части.
        F
        % Сетка с ГУ.
        mesh
        % Шаг по времени (def=0).
        ts=0
        % Число шагов по времени (def=0).
        tsNum=0
        % Время задачи (def=0).
        tDur=0
        % Вывод диагностической информации в командное окно.
        verbose=false
        % Построение сетки и графиков результатов.
        plotting=false
    end
    methods
        % Конструктор с аргументом.
        function obj = StructFEProblem(filename, options)
            if nargin < 2
                options = struct();
            end
            if ~isstruct(options)
                error('StructFEProblem:InvalidOptions', ...
                    'Options must be provided as a struct.');
            end
            if isfield(options, 'verbose')
                if ~isscalar(options.verbose) || ...
                        ~(islogical(options.verbose) || isnumeric(options.verbose)) || ...
                        ~isfinite(options.verbose) || ...
                        ~ismember(double(options.verbose), [0 1])
                    error('StructFEProblem:InvalidOptions', ...
                        'verbose must be a scalar logical value.');
                end
                obj.verbose = logical(options.verbose);
            end
            if isfield(options, 'plotting')
                if ~isscalar(options.plotting) || ...
                        ~(islogical(options.plotting) || isnumeric(options.plotting)) || ...
                        ~isfinite(options.plotting) || ...
                        ~ismember(double(options.plotting), [0 1])
                    error('StructFEProblem:InvalidOptions', ...
                        'plotting must be a scalar logical value.');
                end
                obj.plotting = logical(options.plotting);
            end
            obj.filename = filename;
            % Чтение кейса, создание сетки.
            resourceLimits = struct();
            if isfield(options, 'resourceLimits')
                resourceLimits = options.resourceLimits;
            end
            obj.mesh = FEMesh(obj.filename, resourceLimits);
            % Вывод сетки на график.
            if obj.plotting
                obj.mesh.Plot2DMesh();
            end
            % Печать номеров узлов в элементах.
            if obj.verbose
                obj.mesh.DispNN();
            end
            % Печать матрицы соответствия.
            if obj.verbose
                obj.mesh.DispIM();
            end
            % Создание глобальной матрицы жесткости, вектора F и матрицы
            % масс в зависимости от кол-ва СС на узел.
            % Определение общего числа степеней свободы в системе
            % и сборка разреженных глобальных матриц из элементных триплетов.
            systemDOF = obj.mesh.numberOfDOFs;
            [GlobK, GlobM] = assembleGlobalMatrices( ...
                obj.mesh.allMeshElems, systemDOF);
            obj.K = GlobK;
            obj.M = GlobM;
            % Вектор правой части (сил). Преаллокация без ГУ.
            obj.F = zeros(systemDOF,1);

        end
        % Совместимый метод возвращает разбиение СС, не изменяя МЖ и ММ.
        function [fixedDOFs, freeDOFs] = ApplyFixBC(this)
            model = this.GetAnalysisModel();
            [fixedDOFs, freeDOFs] = partitionDOFs(model);
        end
        % Совместимый метод собирает свежий вектор или историю нагрузок.
        function ApplyForceBC(this)
            model = this.GetAnalysisModel();
            if this.ts == 0
                this.F = buildStaticLoad(model);
            else
                if this.tsNum < 1
                    error('MKEF:InvalidTransientOptions', ...
                        'A positive tsNum is required to build a load history.');
                end
                this.F = buildTransientLoad(model, this.ts, this.tsNum);
            end
        end
        % Возвращает копию данных модели для чистого численного ядра.
        function model = GetAnalysisModel(this)
            model = createAnalysisModel(this.K, this.M, this.mesh);
        end
        % Метод запуска расчета статического нагружения.
        function result = RunStatic(this)
            model = this.GetAnalysisModel();
            result = solveStatic(model);
            if this.verbose
                disp('DOFs (displ. components):')
                disp(result.displacements);
                disp('Reactions (components):')
                disp(result.reactions);
            end
        end
        % Метод запуска расчета собственных колебаний.
        function result = RunModal(this)
            model = this.GetAnalysisModel();
            result = solveModal(model);
            if this.verbose
                disp('Natural frequencies (Hz):')
                disp(result.frequenciesHz);
            end
        end
        % Запуск расчёта, выбранного в необязательной секции analysis.
        function result = RunSelected(this)
            configuration = this.mesh.analysisConfiguration;
            if isempty(fieldnames(configuration))
                error('MKEF:MissingAnalysisConfiguration', ...
                    ['The input file has no analysis section. Use RunStatic, ' ...
                     'RunModal, or RunTransient explicitly.']);
            end
            switch configuration.type
                case 'static'
                    result = this.RunStatic();
                case 'modal'
                    result = this.RunModal();
                case 'transient'
                    result = this.RunTransient(configuration.timeStep, ...
                        configuration.duration, ...
                        configuration.monitorNode, ...
                        configuration.monitorDOF);
                otherwise
                    error('MKEF:UnsupportedAnalysisType', ...
                        'Unsupported configured analysis type %s.', ...
                        configuration.type);
            end
        end

        % Метод запуска анализа динамики со внешними силами.
        function result = RunTransient(this,tStep,tDur,node,dofToPlot,options)
            if nargin < 5 || ~isnumeric(node) || ~isscalar(node) || ...
                    ~isfinite(node) || node ~= fix(node) || node < 1 || ...
                    node > this.mesh.numberOfNodes || ...
                    ~isnumeric(dofToPlot) || ~isscalar(dofToPlot) || ...
                    ~isfinite(dofToPlot) || dofToPlot ~= fix(dofToPlot) || ...
                    dofToPlot < 1 || dofToPlot > this.mesh.dofPerNode
                error('MKEF:InvalidPlotSelection', ...
                    ['Transient node and DOF must be integer indices within ' ...
                     'the model, even when plotting is disabled.']);
            end
            if nargin < 6
                options = struct();
            end
            if ~isstruct(options)
                error('StructFEProblem:InvalidOptions', ...
                    'Transient options must be provided as a struct.');
            end
            options.timeStep = tStep;
            options.duration = tDur;
            model = this.GetAnalysisModel();
            result = solveTransient(model, options);
            if this.plotting
                globalDOF = this.mesh.iMnod(node, dofToPlot);
                if globalDOF == 0
                    error('MKEF:InactiveRotationMonitor','Node %d has no shared thetaZ.',node);
                end
                plotTransientResult(result, globalDOF);
            end
        end
    end
end