ReducedSystem.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 ReducedSystem
  :public models.BaseFirstOrderSystem {
    public:

        R;
    private:

        plotPtr;

        
    public:

        ReducedSystem(rmodel) {
            this = this@models.BaseFirstOrderSystem(rmodel);
        }
        function  build() {
            disp(" Start building reduced system... ");

            rmodel = this.Model;
            fullmodel = rmodel.FullModel;
            fullsys = fullmodel.System;

            V = rmodel.V;

            if ~isempty(fullmodel.SpaceReducer) && isempty(V)
                error([" Model has a SpaceReducer but no projection "...
                        " matrix V is given. Forgot to call offlineGenerations? "]);
            end

            W = rmodel.W;
            /*  Use V if W is empty (Galerkin projection) */
            if ~isempty(V) && isempty(W)
                W = V;
            end

            this.updateDimensions;

            /*  The state scaling for the reduced system is one, as the scaling matrices have been
             * incorporated into the V,W matrices via V := SV and W = S^-1W inside the ReducedModel
             * setModel method. */
            this.StateScaling= 1;
            SV = 1; SW = 1; /*  Default: no scaling */

            s = fullmodel.System.StateScaling;
            if s ~= 1
                if isscalar(s)
                    dim = fullsys.NumStateDofs;
                    s(1:dim,1) = s;
                else
                    dim = length(s);
                end
                SV = spdiags(s,0,dim,dim);
                SW = spdiags(1./s,0,dim,dim);
            end
            /*  Incorporate the state scaling into the projection matrices.
             * They are used to project the x0 initial values and output C. */
            if ~isempty(V)
                SV = SV * V;
                SW = SW * W;
            end

            /*  Clones the full system's basic (all but functions)
             * properties */
            this.Params= fullsys.Params;
            this.Inputs= fullsys.Inputs;
            this.MaxTimestep= fullsys.MaxTimestep;

            /*  Copy component handles (it's NOT cloning them!)
             * This is the default as if no reduction methods are
             * applied (sensless but yet allowed) an identical system
             * will be the result. */
            this.f= fullsys.f;
            this.B= fullsys.B;

            /*  SV ~= 1 means that projection or nontrivial scaling is used. */
            if ~isequal(SV,1) 
                if ~isempty(fullsys.x0)  /*  otherwise error for static models */

                    this.x0= fullsys.x0.project(SV,SW);
                end
                SR = this.compileReconstructionMatrix(SV);
                this.C= fullsys.C.project(SR,SW);
            else
                if ~isempty(fullsys.x0)  /*  otherwise error for static models */

                    this.x0= fullsys.x0.clone;
                end
                this.C= fullsys.C.clone;
            end

            /*  Check whether projection was setup for this system */
            if ~isempty(V)
                this.R= this.compileReconstructionMatrix(V);

                /*  Project input */
                if ~isempty(fullsys.B)
                    this.B= fullsys.B.project(V,W);
                end
                /*  Project linear part */
                if ~isempty(fullsys.A)
                    this.A= fullsys.A.project(V,W);
                end
                /*  Project the approximated CoreFun of the full model if exists */
                if ~isempty(fullmodel.Approx)
                    this.f= fullmodel.Approx.project(V,W);
                elseif ~isempty(fullsys.f)
                    /*  Otherwise project the models' full function. */
                    this.f= fullsys.f.project(this.R,W);
                end
                /*  Project mass matrix */
                if ~isempty(fullsys.M)
                    this.M= fullsys.M.project(V,W);
                end
                /*  Project algebraic constraints function
                 * This needs only affect the getStateJacobian function of
                 * ACoreFun */
                if ~isempty(fullsys.g)
                    this.g= fullsys.g.project(this.R,1);
                    this.g.setSystem(this);
                end
            else
                /*  Only use approximated version if set */
                if ~isempty(fullmodel.Approx)
                    this.f= fullmodel.Approx;
                end
            end
            /*  Set the reduced system as base system for the projected
             * ACoreFun - it will need to change references to the full
             * system itself if required */
            this.f.setSystem(this);

            /*  Set the plot-wrapper (uses the plot method from the full
             * system) */
            this.plotPtr= @fullsys.plot;
        }
        function  updateSparsityPattern() {
            this.SparsityPattern= [];
        }
        function  plot(models.BaseFullModel model,double t,matrix<double> y) {
            this.plotPtr(model, t, y);
        }
        function dx = ODEFun(double t,colvec<double> x) {
            est = this.Model.ErrorEstimator;
            haveest = ~isempty(est) && est.Enabled;
            if haveest
                xall = x;
                x = x(1:end-est.ExtraODEDims,:);    
            end

            dx = ODEFun@models.BaseFirstOrderSystem(this,t,x);

            if haveest
                dx = [dx; est.evalODEPart(xall, t, this.u(t))];
            end
        }


        function J = getJacobian(double t,xc) {
            est = this.Model.ErrorEstimator;
            haveest = ~isempty(est) && est.Enabled;
            if haveest
                xc = xc(1:end-est.ExtraODEDims,:);    
            end
            J = getJacobian@models.BaseFirstOrderSystem(this, t, xc);
            if haveest
                J = blkdiag(J,0);
            end
        }


        function x0 = getX0(colvec<double> mu) {

            x0 = getX0@models.BaseFirstOrderSystem(this, mu);
            m = this.Model;
            if ~isempty(m.ErrorEstimator) && m.ErrorEstimator.Enabled
                x0 = [x0; m.ErrorEstimator.getE0(mu)];
            end
        }
    protected:


        function R = compileReconstructionMatrix(V) {
            R = V;
        }


        function  updateDimensions() {
            rm = this.Model;
            if ~isempty(rm.V)
                /*  This is the projection case */
                this.NumStateDofs= size(rm.V,2);
            else
                /*  No projection case */
                this.NumStateDofs= rm.FullModel.System.NumStateDofs;
            end
            this.NumAlgebraicDofs= rm.FullModel.System.NumAlgebraicDofs;
            this.NumTotalDofs= this.NumStateDofs + this.NumAlgebraicDofs;
        }


        function  validateModel(models.BaseFullModel model) {
            if ~isa(model, " models.ReducedModel ")
                error(" The Model property has to be a child of models.ReducedModel ");
            end
        }
};
}