Model.m

Go to the documentation of this file.

namespace models{
namespace muscle{


/* (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 Model
  :public models.BaseFullModel {
    public:

        MuscleDensity = 1.1e-3;
        Plotter;
    public:

        .int RandSeed = 1;
        Config;


            Gravity = 9.80665e-3;
    public: /* ( Dependent ) */

        Geo;

        
    public: /* ( Dependent ) */

        Model(conf,basedir) {
            if nargin < 2
                basedir = KerMor.App.DataDirectory;
                if nargin < 1
                    conf = Debug;
                end
            end
            /*  Creates a new muscle model */
            name = sprintf(" FEM Muscle model %s ",class(conf));
            this = this@models.BaseFullModel(name);

            this.SaveTag= sprintf(" musclemodel_%s ",class(conf));
            this.Data= data.ModelData(this,basedir);

            this.System= models.muscle.System(this);
            /*  Sets DefaultMu */
            this.initDefaultParameter;

            this.TrainingParams= [1 2];

            this.T= 10; /*  [ms] */

            this.dt= .01; /*  [ms] */

            s = solvers.MLode15i;
            /*  Use relatively coarse precision settings. This skips fine
             * oscillations but yields the correct long time results. */
            s.RelTol= 1e-1;
            s.AbsTol= 1e-1;
            this.ODESolver= s;
            this.System.MaxTimestep= []; /* model.dt; */


            /*  MOR pre-setup.
             * If you assign a new SpaceReducer instance, dont forget to set
             * the TargetDimensions property accordingly (is now
             * automatically set within the configureModel base routine)
             *s = spacereduction.PODGreedy;
             *s.Eps = 1e-9;
             *s.MaxSubspaceSize = 500; */
            s = spacereduction.PODReducer;
            s.IncludeInitialSpace= true;
            s.Mode= " abs ";
            this.SpaceReducer= s;

            /*  Call the config-specific model configuration */
            conf.configureModel(this);

            /*  Set the config to the model, triggering geometry related
             * pre-computations */
            this.setConfig(conf);

            conf.configureModelFinal;

            /* % Health tests
             * Propagate the current default param */
            this.System.prepareSimulation(this.DefaultMu, this.DefaultInput);

/*              fprintf('Running Jacobian health check..');
 *             res = this.System.f.test_Jacobian;
 *             fprintf('Done. Success=%d\n',res);
 *             chk = this.Config.FEM.test_JacobiansDefaultGeo;
 * %             chk = chk && this.Config.FEM.test_QuadraticBasisFun;
 *             chk = chk && this.Config.PressureFEM.test_JacobiansDefaultGeo;
 *             if ~chk
 *                 error('Health tests failed!');
 *             end */
        }


        function [doublet , matrix<double>y ] = simulateAndPlot(withResForce,varargin) {
            if nargin < 2
                withResForce = true;
            end
            [t,y] = this.simulate(varargin[:]);
            xargs = [];
            if (withResForce)
                [df,nf] = this.getResidualForces(t,y);
                xargs = [" NF ",nf," DF ",df];
            end
            this.plot(t,y,xargs[:]);
        }


        function [doublet , colvec<double>x , time , cache ] = computeTrajectory(colvec<double> mu,integer inputidx) {
            this.Config.prepareSimulation(mu, inputidx);
            [t, x, time, cache] = computeTrajectory@models.BaseFullModel(this, mu, inputidx);
            f = this.System.f;
            fprintf(" Finished after %gs (fevals:%d, Jacobians: %d)\n ",time,f.nfevals,f.nJevals);
        }
        function doublet = getConfigTable(colvec<double> mu) {
            if nargin < 2
                mu = this.DefaultMu;
            end
            f = this.System.f;
            t = PrintTable(" Configuration of Model %s ",this.Name);
            t.HasHeader= true;
            t.addRow(" \rho_0 ", " c_{10} [MPa] "," c_{01} [MPa] "," b_1 [MPa] "," d_1 [-] "," P_{max} [MPa] "," \lambda_f^{opt} ");
            t.addRow(this.MuscleDensity, mu(9),mu(10),mu(5),mu(6),mu(13),mu(14));
            t.Format= " tex ";
        }


        function  plotForceLengthCurve(colvec<double> mu,pm) {
            if nargin < 3
                pm = PlotManager(false,2,2);
                pm.LeaveOpen= true;
                if nargin < 2
                    mu = this.DefaultMu;
                end
            end
            sys = this.System;
            f = sys.f;
            sys.prepareSimulation(mu,this.DefaultInput);
            this.Config.setForceLengthFun(f);

            lambda = .2:.005:2;

            /* % Plain Force-length function */
            fl = f.ForceLengthFun(lambda);
            h = pm.nextPlot(" force_length_plain "," Direct force-length curve of muscle/sarcomere ",...
                sprintf(" \\lambda [-], fl_p1=%g ",mu(14)),...
                " force-length relation [-] ");
            plot(h,lambda,fl," r ");

            /* % Effective force-length function for muscle tissue
             * effective signal from active part */
            fl_eff = (mu(13)./lambda) .* fl;
            /*  Passive markert laws */
            aniso_passive_muscle = f.AnisoPassiveMuscle(lambda);

            /*  Find a suitable position to stop plotting (otherwise the
             * passive part will steal the show) */
            pos = find(aniso_passive_muscle > max(fl_eff)*1.4,1," first ");
            if ~isempty(pos)
                lambda = lambda(1:pos);
                fl_eff = fl_eff(1:pos);
                aniso_passive_muscle = aniso_passive_muscle(1:pos);
            end

            h = pm.nextPlot(" force_length_eff ",...
                " Effective force-length curve of muscle material ",...
                " \lambda [-] "," pressure [MPa] ");
            plot(h,lambda,fl_eff," r ",lambda,aniso_passive_muscle," g ",lambda,fl_eff+aniso_passive_muscle," b ");
            legend(h," Active "," Passive "," Total "," Location "," NorthWest ");

/*              %% Effective force-length function derivative for muscle tissue
 *             dfl = (mu(13)./lambda) .* f.ForceLengthFunDeriv(lambda);
 *             dmarkertf = (lambda>=1).*(b1./lambda.^3).*((d1-2)*lambda.^d1 + 2);
 *             h = pm.nextPlot('force_length_eff_deriv',...
 *                 'Effective force-length curve derivative of muscle material',...
 *                 '\lambda','deriv [MPa/ms]');
 *             plot(h,lambda,dfl,'r',lambda,dmarkertf,'g',lambda,dfl + dmarkertf,'b'); */

/*              %% Cross fibre stuff
 *             if this.System.UseCrossFibreStiffness
 *                 error('fixme');
 *                 markertf = max(0,(f.b1cf./lambda.^2).*(lambda.^f.d1cf-1));
 *                 h = pm.nextPlot('force_length_xfibre',sprintf('Force-Length curve in cross-fibre direction for model %s',this.Name),'\lambda [-]','pressure [MPa]');
 *                 plot(h,lambda,markertf,'r');
 *                 axis(h,[0 2 0 150]);
 *                 legend(h,'Passive cross-fibre pressure','Location','NorthWest');
 *                 
 *                 dmarkertf = (lambda >= 1).*(f.b1cf./lambda.^3).*((f.d1cf-2)*lambda.^f.d1cf + 2);
 *                 h = pm.nextPlot('force_length_xfibre_deriv',sprintf('Derivative of Force-Length curve in cross-fibre direction for model %s',this.Name),'\lambda [-]','pressure [MPa]');
 *                 plot(h,lambda,dmarkertf,'r');
 *             end */

            /* % Passive force-length function for 100% tendon tissue
             * Passive markert law */
            if sys.HasTendons
                aniso_passive_tendon = f.AnisoPassiveTendon(lambda);
                h = pm.nextPlot(" force_length_tendon ",...
                    " Effective force-length curve of tendon material (=passive) ",...
                    " \lambda [-] "," pressure [MPa] ");
                plot(h,lambda,aniso_passive_tendon," g ");

                /* % Effective force-length surface for muscle-tendon tissue
                 * Sampled ratios */
                tmr = 0:.03:1;
                /* tmrlog = [0 logspace(-4,0,length(tmr)-1)]; */
                tmrlog = tmr;

                /* % Passive part */
                [LAM,TMR] = meshgrid(lambda,tmr);
                passive = repmat(aniso_passive_muscle,length(tmr),1) ...
                    + tmrlog^t*(aniso_passive_tendon-aniso_passive_muscle);

                /* % Active part */
                FL = (1-tmr)^t*fl_eff;

                /* % Plot dat stuff! */
                h = pm.nextPlot(" force_length_muscle_tendon ",...
                    " Effective force-length curve between muscle/tendon material ",...
                    " \lambda [-] "," muscle/tendon ratio [m=0,t=1] ");
                surfc(LAM,TMR,passive+FL," Parent ",h," EdgeColor "," interp ");
                zlabel(" pressure [MPa] ");
                zlim([0, 3*max(fl_eff(:))]);
            end

            if nargin < 2
                pm.done;
            end
        }


        function  plotAnisotropicPressure(colvec<double> mu) {
            if nargin < 2
                mu = this.DefaultMu;
            end
/*              pm = PlotManager(false,1,2); */
            pm = PlotManager;
            pm.LeaveOpen= true;
            f = this.System.f;

            /*  Simply use muscle parameter values */
            b1 = mu(1);
            d1 = mu(2);
            warning(" Using muscle parameters only (ignoring tendon) ");

            [lambda, alpha] = meshgrid(.02:.02:1.5,0:.01:1);

            fl = f.ForceLengthFun(lambda/mu(14));
            active = mu(13)./lambda.*fl.*alpha;/*  + 1*max(0,(lambda-1)).*alpha; */

            passive = max(0,(b1./lambda.^2).*(lambda.^d1-1));

            h = pm.nextPlot(" aniso_pressure ",sprintf(" Pressure in fibre direction for model %s ",this.Name)," Stretch \lambda "," Activation \alpha ");
            surf(h,lambda,alpha,active+passive," EdgeColor "," k "," FaceColor "," interp ");

            pm.done;
        }


        function  plotActivation() {
            pm = PlotManager;
            pm.LeaveOpen= true;
            f = this.System.f;

            h = pm.nextPlot(" activation ",sprintf(" Activation curve for model %s ",this.Name)," time [ms] "," alpha [-] ");
            plot(h,this.Times,f.alpha(this.scaledTimes)," r ");
            pm.done;
        }


        function varargout = plotGeometrySetup(varargin) {
            if isempty(this.System.mu)
                this.System.prepareSimulation(this.DefaultMu, this.DefaultInput);
            end
            if ~isempty(varargin) && isa(varargin[1]," PlotManager ")
                varargin = [[" PM "] varargin];
            end
            x0 = this.System.getX0(this.System.mu);
            [~, nf] = this.getResidualForces(0, x0);
            if ~isempty(nf)
                varargin(end+1:end+2) = [" NF ",nf];
            end
            varargin(end+1:end+2) = [" Velo ",true];
            /*  Plot without default args (time-plotting might want to
             * suppress fibres for speed, but we want them here) */
            old = this.Plotter.DefaultArgs;
            this.Plotter.DefaultArgs= [];
            [varargout[1:nargout]] = this.plot(0,x0,varargin[:]);
            this.Plotter.DefaultArgs= old;
        }


        function  plotGeometryInfo(varargin) {
            this.Config.plotGeometryInfo(varargin[:]);
        }


        function [residuals_dirichlet , residuals_neumann ] = getResidualForces(double t,uvw) {
            sys = this.System;
            num_bc = length(sys.idx_u_bc_glob)+length(sys.idx_expl_v_bc_glob);
            residuals_dirichlet = zeros(num_bc,length(t));
            residuals_neumann = zeros(length(sys.idx_neumann_bc_glob),length(t));
            for k=1:length(t)
                dy = sys.f.evaluate(uvw(:,k),t(k));
                if ~isempty(residuals_neumann)
                    residuals_neumann(:,k) = dy(sys.idx_neumann_bc_dof);
                end
                residuals_dirichlet(:,k) = sys.f.LastBCResiduals;
            end
        }


        function colvec<integer>idx = getFaceDofsGlobal(integer elem,rowvec<integer> faces,rowvec<integer> dim) {
            if nargin < 4
                dim = 1:3;
            end
            geo = this.Config.FEM.Geometry;
            idxXYZ = false(3,geo.NumNodes);
            for k=1:length(faces)
                idxXYZ(dim,geo.Elements(elem,geo.MasterFaces(faces(k),:))) = true;
            end
            idx = find(idxXYZ(:));
        }
        function idx = getPositionDirichletBCFaceIdx(integer elem,integer face,rowvec<integer> dim) {
            if nargin < 4
                dim = 1:3;
            end
            geo = this.Config.FEM.Geometry;
            idx_face = false(size(this.System.bool_u_bc_nodes));
            idx_face(dim,geo.Elements(elem,geo.MasterFaces(face,:))) = true;
            fidx = find(this.System.bool_u_bc_nodes & idx_face);
            [~, idx] = intersect(this.System.idx_u_bc_glob, fidx);
        }
        function idx = getVelocityDirichletBCFaceIdx(integer elem,integer face,rowvec<integer> dim) {
            if nargin < 4
                dim = 1:3;
            end
            geo = this.Config.FEM.Geometry;
            idx_face = false(size(this.System.bool_u_bc_nodes));
            idx_face(dim,geo.Elements(elem,geo.MasterFaces(face,:))) = true;
            fidx = find(this.System.bool_expl_v_bc_nodes & idx_face);
            /*  Include the offset to the indices for velocity dofs */
            fidx = fidx + geo.NumNodes*3;
            [~, idx] = intersect(this.System.idx_v_bc_glob, fidx);
            /*  Also include the offset of velocity components within the
             * dirichlet force vector */
            idx = idx + length(this.System.val_u_bc);
        }
        function  setConfig(value) {
            if ~isa(value, " models.muscle.AMuscleConfig ")
                error(" Config must be a models.muscle.AMuscleConfig instance ");
            end
            this.Config= value;
            this.System.configUpdated;
            this.Plotter= models.muscle.MusclePlotter(this.System);
        }


        function  setGaussIntegrationRule(value) {
            mc = this.Config;
            mc.FEM.GaussPointRule= value;
            mc.PressureFEM.GaussPointRule= value;
            this.setConfig(mc);
            s = this.System;
            mu = s.mu;
            if isempty(mu)
                mu = this.DefaultMu;
            end
            in = s.inputidx;
            if isempty(in)
                in = this.DefaultInput;
            end
            s.prepareSimulation(mu,in);
        }
        function varargout = plot(varargin) {
            [varargout[1:nargout]] = this.Plotter.plot(varargin[:]);
        }


        
#if 0 //mtoc++: 'get.Geo'
function value = Geo() {
            value = this.Config.Geometry;
        }

#endif


    
    public: /* ( Static ) */ /* ( Dependent ) */

        static function res = test_ModelVersions() {
            res = true;
            try
                m = models.muscle.Model(models.muscle.examples.Debug);
                mu = m.getRandomParam;
                [t,y] = m.simulate(mu);
                m.System.UseDirectMassInversion= true;
                [t,y] = m.simulate(mu);

                /*  Version with "constant" fibre activation forces */
                m = models.muscle.Model(models.muscle.examples.Debug(2));
                [t,y] = m.simulate(mu);
                m.System.UseDirectMassInversion= true;
                [t,y] = m.simulate(mu);

                /*  Version with "neurophysiological" fibre activation forces */
                m = models.muscle.Model(models.muscle.examples.Debug(3));
                [t,y] = m.simulate(mu);
                m.System.UseDirectMassInversion= true;
                [t,y] = m.simulate(mu);
                m.System.UseDirectMassInversion= false;

                /*  "Disable" viscosity */
                m.DefaultMu(1) = 0;
                [t,y] = m.simulate(mu);
                m.System.UseDirectMassInversion= true;
                [t,y] = m.simulate(mu);
            catch ME
                ME.getReport(" extended ");
                res = false;
            end
        }


        static function res = test_JacobianApproxGaussRules() {
            m = models.muscle.Model(models.muscle.examples.Debug(2));
            m.simulate;
            f = m.System.f;
            m.dt= .2;
            m.T= 1;
            res = f.test_Jacobian;
            m.setGaussIntegrationRule(4);
            res = res && f.test_Jacobian;
            m.setGaussIntegrationRule(5);
            res = res && f.test_Jacobian;
        }
    protected: /* ( Static ) */ /* ( Dependent ) */

        static function this = loadobj() {
            if ~isa(this, " models.muscle.Model ")
                sobj = this;
                this = models.muscle.Model;
                this.RandSeed= sobj.RandSeed;
                this.Config= sobj.Config;
                this = loadobj@models.BaseFullModel(this, sobj);
            else
                this = loadobj@models.BaseFullModel(this);
            end
        }

};
}
}