BaseSecondOrderSystem.m

Go to the documentation of this file.

namespace models{


/* (Autoinserted by mtoc++)
 * This source code has been filtered by the mtoc++ executable,
 * which generates code that can be processed by the doxygen documentation tool.
 *
 * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
 * Except for the comments, the function bodies of your M-file functions are untouched.
 * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
 * attached source files that are highly readable by humans.
 *
 * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
 * the correct locations in the source code browser.
 * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
 */

class BaseSecondOrderSystem
  :public models.BaseFirstOrderSystem {
    public: /* ( setObservable ) */

        dscomponents.AffLinCoreFun,dscomponents.LinearCoreFun D;
    public:

        NumDerivativeDofs;


        DerivativeDirichletPosInStateDofs;
    protected:

        JDX;
    public:

        BaseSecondOrderSystem(models.BaseFullModel model) {
            this = this@models.BaseFirstOrderSystem(model);

            /*  Register default properties */
            this.registerProps(" D ");
        }
        function rsys = getReducedSystemInstance(models.ReducedModel rmodel) {
            rsys = models.ReducedSecondOrderSystem(rmodel);
        }
        function  prepareSimulation(colvec<double> mu,integer inputidx) {
            prepareSimulation@models.BaseFirstOrderSystem(this, mu, inputidx);

            /*  Forward preparation call to D, if present */
            if ~isempty(this.D)
                this.D.prepareSimulation(mu);
            end
        }


        function dx_xdot_c = ODEFun(double t,x_xdot_c) {
            dx_xdot_c = zeros(size(x_xdot_c));

            num_x_dof = this.NumStateDofs;
            num_xdot_dof = this.NumDerivativeDofs;
            x = x_xdot_c(1:num_x_dof);
            xdot = x_xdot_c(num_x_dof+(1:num_xdot_dof));

            /* % Set x'=xdot
             * Check for derivative dirichlet conditions - insert them here
             * as direct values */
            if num_x_dof ~= num_xdot_dof /*  Quickest check */

                dx = zeros(num_x_dof,1);
                /*  Execute callback - subclasses must provide the (possibly
                 * time-dependent) valuess */
                val = this.getDerivativeDirichletValues(t);
                pos = this.DerivativeDirichletPosInStateDofs;
                /*  Dirichlet values */
                dx(pos) = val;
                /*  Dof values */
                dx(~pos) = xdot;
                /*  Set in global dof vector */
                dx_xdot_c(1:num_x_dof) = dx;
            else
                /*  Direct forwarding - no derivative dirichlet conditions */
                dx_xdot_c(1:num_x_dof) = xdot;
            end

            /* % Set x''=xdot' = A(x)+D(y)+f(x,y)+Bu */
            dy = zeros(num_xdot_dof,1);
            if ~isempty(this.A)
                dy = dy + this.A.evaluate(x, t);
            end
            if ~isempty(this.D)
                dy = dy + this.D.evaluate(xdot, t);
            end
            if ~isempty(this.f)
                dy = dy + this.f.evaluate(x_xdot_c, t);
            end
            if ~isempty(this.B) && ~isempty(this.inputidx)
                B = this.B.evaluate(t, this.mu)*this.u(t);
                dy = dy + B;
            end
            dx_xdot_c(num_x_dof+(1:num_xdot_dof)) = dy;
            /*  Append algebraic constraint evaluation */
            if ~isempty(this.g)
                dx_xdot_c(num_x_dof+num_xdot_dof+1:end) = this.g.evaluate(x_xdot_c,t);
            end
        }
        function J = getJacobian(double t,x_xdot_c) {
            td = this.NumTotalDofs;
            sd = this.NumStateDofs;
            dd = this.NumDerivativeDofs;
            ad = this.NumAlgebraicDofs;
            xdotpos = sd+(1:dd);
            xpos = 1:sd;
            cpos = sd+dd+(1:ad);
            x = x_xdot_c(xpos);
            xdot = x_xdot_c(xdotpos);
            J = sparse(td,td);
            J(xpos,:) = this.JDX;
            if ~isempty(this.A)
                J(xdotpos,xpos) = J(xdotpos,xpos) + this.A.getStateJacobian(x, t);
            end
            if ~isempty(this.D)
                J(xdotpos,xdotpos) = J(xdotpos,xdotpos) + this.D.getStateJacobian(xdot, t);
            end
            if ~isempty(this.f)
                J(xdotpos,:) = J(xdotpos,:) + this.f.getStateJacobian(x_xdot_c, t);
            end
            if ~isempty(this.g)
                J(cpos,:) = this.g.getStateJacobian(x_xdot_c,t);
            end
        }
        function  updateSparsityPattern() {
            td = this.NumTotalDofs;
            sd = this.NumStateDofs;
            dd = this.NumDerivativeDofs;
            ad = this.NumAlgebraicDofs;
            xdotpos = sd+(1:dd);
            JP = logical(sparse(td,td));
            I = speye(sd);
            I(:,this.DerivativeDirichletPosInStateDofs) = [];
            this.JDX= [sparse(sd,sd) I sparse(sd,ad)];
            JP(1:sd,:) = this.JDX;
            if ~isempty(this.A) && ~isempty(this.A.JSparsityPattern)
                JP(xdotpos,1:sd) = JP(xdotpos,1:sd) | this.A.JSparsityPattern;
            end
            if ~isempty(this.D) && ~isempty(this.D.JSparsityPattern)
                JP(xdotpos,xdotpos) = JP(xdotpos,xdotpos) | this.D.JSparsityPattern;
            end
            if ~isempty(this.f) && ~isempty(this.f.JSparsityPattern)
                JP(xdotpos,:) = JP(xdotpos,:) | this.f.JSparsityPattern;
            end
            if ~isempty(this.g)
                JP(sd+dd+1:end,:) = this.g.JSparsityPattern;
            end
            this.SparsityPattern= JP;
        }
        function x_xdot_c0 = getX0(colvec<double> mu) {

            /*  Gets the initial state variable at `t=0`. */
            x_c0 = getX0@models.BaseFirstOrderSystem(this, mu);
            num_x_dof = this.NumStateDofs;
            num_xdot_dof = this.NumDerivativeDofs;
            x_xdot_c0 = zeros(this.NumTotalDofs,1);

            /*  Insert state initial values */
            x_xdot_c0(1:num_x_dof) = x_c0(1:num_x_dof);

            /*  Insert derivative initial values (xdot0) between */
            x_xdot_c0(num_x_dof+(1:num_xdot_dof)) = this.getDerivativeInitialValues(mu);

            /*  Insert alg cond initial values */
            x_xdot_c0(num_x_dof+num_xdot_dof+1:end) = x_c0(num_x_dof+1:end);
        }
        function M = getMassMatrix() {
            M = this.M;
            if isempty(M)
                M = speye(this.NumDerivativeDofs);
            else
                if ~isa(M," dscomponents.ConstMassMatrix ")
                    error(" Non-constant mass matrices not yet implemented for second order systems ");
                end
                M = this.M.M;
            end    
            sd = this.NumStateDofs;
            ad = this.NumAlgebraicDofs;
            M = blkdiag(speye(sd),M,sparse(ad,ad));
            /*  Tell the mass matrix which components are algebraic
             * constraint dofs - those wont be used determinig a suitable
             * initial slope in e.g. solvers.MLWrapper */
            algdofs = sd+this.NumDerivativeDofs+(1:ad);
            M = dscomponents.ConstMassMatrix(M,algdofs);
        }

    
    protected:


        function  updateDimensions() {
            this.NumTotalDofs= this.NumStateDofs ...
                + this.NumDerivativeDofs + this.NumAlgebraicDofs;
        }


        function  validateModel(models.BaseFullModel model) {
            if ~isa(model, " models.BaseModel ")
                error(" The Model property has to be a child of models.BaseModel ");
            end
        }
        function val = getDerivativeInitialValues(unused1) {
            val = zeros(this.NumDerivativeDofs,1);
        }
    protected: /* ( Abstract ) */

        virtual function val = getDerivativeDirichletValues(double t) = 0;
};
}