Debug.m

Go to the documentation of this file.

namespace models{
namespace muscle{
namespace examples{


/* (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 Debug
  :public models.muscle.AMuscleConfig {
    public:

        Debug(varargin) {
            if nargin == 1 && isscalar(varargin[1])
                varargin = [" Version ",varargin[1]];
            end
            /*  Creates a Debug simple muscle model configuration. */
            this = this@models.muscle.AMuscleConfig(varargin[:]);
            this.addOption(" Version ",3);
            this.init;

            switch this.Options.Version
                case [4,5,6]
                    this.VelocityBCTimeFun= general.functions.ConstantUntil(15);
                    this.Options.FL= 2;
                case [7,8,9]
                    this.VelocityBCTimeFun= general.functions.ConstantUntil(15);
                case 14
                    /*  Use classical FL function with fl(1)=1 */
                    this.Options.FL= 3;
            end
        }
        function  configureModel(m) {
            configureModel@models.muscle.AMuscleConfig(this, m);
/*              sys = m.System;
 *             f = sys.f; */
            m.T= 50;
            m.dt= m.T / 30;
            m.ODESolver.RelTol= 1e-5;
            m.ODESolver.AbsTol= 1e-3;
            switch this.Options.Version
            case 1
                m.DefaultMu(2) = -1;
            case 2
                m.DefaultMu(1:2) = [1; m.T/2];
            case 3
                m.T= 400;
                m.dt= 1;
                types = [0 .2 .4 .6 .8 1];
                fe = this.FEM;
                geo = fe.Geometry;
                ftw = zeros(fe.GaussPointsPerElem,length(types),geo.NumElements);
                /*  Test: Use only slow-twitch muscles */
                ftw(:,1,:) = .4;
                ftw(:,2,:) = .05;
                ftw(:,3,:) = .05;
                ftw(:,4,:) = .1;
                ftw(:,5,:) = .2;
                ftw(:,6,:) = .2;
                this.FibreTypeWeights= ftw;
                p = models.motorunit.Pool;
                p.FibreTypes= types;
                this.Pool= p;
/*                  m.ODESolver.RelTol = .01;
 *                 m.ODESolver.AbsTol = .1; */
                m.Plotter.DefaultArgs= [" Pool ",true];
                m.DefaultMu(4) = 4;
            case [4,5,6]
                /* % Material configuration from CMISS/3Elem_sprenger.xml                
                 * c3M = 3.05907e-10 MPa */
                m.DefaultMu(5) = 3.05907e-10; /*  b1 [MPa]
                 * c4M = 47.270456264135881 */

                m.DefaultMu(6) = 47.270456264135881; /*  d1 [-]
                 * c1M = 3.56463903963e-02 MPa */

                m.DefaultMu(9) = 35.6463903963; /*  c10 [MPa]
                 * c2M = 3.859558659683e-3 MPa */

                m.DefaultMu(10) = 3.859558659683; /*  c01 [MPa] */


                /*  sigma_max_calculation  */
                m.DefaultMu(13) = .300; /*  [MPa]
                 * lambda_ofl_calculation */

                m.DefaultMu(14) = 1.4; /*  [-] */

                m.T= 20;
                m.dt= .25;
            case [7,8,9]
                /* % Using the Material configuration of Heidlauf */
                m.DefaultMu(5) = 0.00355439810963035e-3; /*  b1 [MPa] */

                m.DefaultMu(6) = 12.660539325481963; /*  d1 [-] */

                m.DefaultMu(9) = 6.352e-13; /*  c10 [MPa] */

                m.DefaultMu(10) = 3.627e-3; /*  c01 [MPa] */

                m.T= 50;
                m.dt= .5;
            case 10
                m.ODESolver.RelTol= .001;
                m.ODESolver.AbsTol= .03;
                m.DefaultMu(1:3) = [0;0;1];
                m.DefaultInput= 1;
                m.T= 1;
                m.dt= .01;
            case 11
                m.ODESolver.RelTol= .001;
                m.ODESolver.AbsTol= .05;
                m.T= 10;
                m.dt= .05;
            case 12
                m.DefaultMu(1:2) = [.1; 50];
                m.System.UseCrossFibreStiffness= true;
            end
        }


        function  configureModelFinal() {
            configureModelFinal@models.muscle.AMuscleConfig(this)
            m = this.Model;
            /*  Have the PODReducer always generate the full reduced model */
            m.SpaceReducer.Value= length(m.SpaceReducer.TargetDimensions);
        }


        function P = getBoundaryPressure(elemidx,faceidx) {
            P = [];
            if this.Options.Version == 10
                if elemidx == 1 && faceidx == 2
                   P = 1;
                end
            end
        }
        function u = getInputs() {
            u = [];
            if this.Options.Version == 10
                u = [this.getAlphaRamp(this.Model.T/2,.1)];
            end
        }

    
    protected:


        function geo = getGeometry() {
            geo = fem.geometry.RegularHex8Grid([0 1],[0 1],[0 1]);
            /* [pts, cubes] = fem.geometry.RegularHex8Grid([0:.5:1],[0:.5:1],[0:.5:1]); */
            if this.Options.Version == 11 || this.Options.Version == 12
                geo = fem.geometry.RegularHex8Grid([-1 1],[-1 1],[-1 1]);
                pts = geo.Nodes;
                cubes = geo.Elements;
                pts(1,[6 8]) = 2;
                theta = .3;
                R = [cos(theta) -sin(theta) 0
                     sin(theta) cos(theta)  0
                     0          0           1];
                pts = circshift(R,[1 1])*R*pts;
                geo = fem.geometry.Cube8Node(pts, cubes);
            end
            geo = geo.toCube27Node;
        }
        function displ_dir = setPositionDirichletBC(displ_dir) {
            geo = this.FEM.Geometry;
            /*  Always fix back side */
            switch this.Options.Version
                case 4
                    displ_dir(:,geo.Elements(1,geo.MasterFaces(4,:))) = true;
                case 10
                    displ_dir(1,geo.Elements(1,geo.MasterFaces(1,:))) = true;
                case 14
                    displ_dir(:,geo.Elements(1,geo.MasterFaces(1:2,:))) = true;
                otherwise
                    displ_dir(:,geo.Elements(1,geo.MasterFaces(1,:))) = true;
            end
        }


        function [velo_dir , velo_dir_val ] = setVelocityDirichletBC(velo_dir,velo_dir_val) {
            geo = this.FEM.Geometry;
            switch this.Options.Version

            case 4
                velo_dir(2,geo.Elements(1,geo.MasterFaces(3,:))) = true;
                velo_dir_val(velo_dir) = -.01;
            case [5,6]
                velo_dir(1,geo.Elements(1,geo.MasterFaces(2,:))) = true;
                velo_dir_val(velo_dir) = .01;
            case [7,8,9]
                velo_dir(1,geo.Elements(1,geo.MasterFaces(2,:))) = true;
                velo_dir_val(velo_dir) = .04;
            case 11
                velo_dir(:,geo.Elements(1,geo.MasterFaces(2,:))) = true;
                hlp = velo_dir;
                hlp([1 3],:) = 0;
                velo_dir_val(hlp) = .01;
            end
        }
        function anull = seta0(anull) {
            switch this.Options.Version
            case [2,3,5,8,12,14]
                anull(1,:,:) = 1;
            case [4,7]
                /*  No fibres */
            case [6,9]
                /*  Stretch perpendicular to fibre direction */
                anull(2,:,:) = 1;
            case [10,11]
                anull(:,:,:) = 0;
            end
        }

    
    public: /* ( Static ) */

        static function  test_DebugConfigV4_9(version) {
            if nargin < 1
                version = 4;
            end
            m = models.muscle.Model(models.muscle.examples.Debug(version));
            [t,y] = m.simulate;
            pdf = m.getResidualForces(t,y);

            /*  1..27 are residual forces from fixed position nodes */
            force_on_fixed_side = sum(pdf(1:27,:),1);
            /*  28..36 (=9) are residual forces from moved nodes */
            force_on_moved_side = sum(pdf(28:end,:),1);
            figure;
            plot(t,force_on_fixed_side," r ",t,force_on_moved_side," b ");
            xlabel(" time "); ylabel(" force [N] ");

            /*  Node 15 is center -> dof indices (15-1)*3 + [1 2 3] */
            yall = m.System.includeDirichletValues(t,y);
            xpos = yall((15-1)*3+1,end);
            fprintf(" Force at T=%g, x-position (center node) %g: %gN\n ",t(end),xpos,force_on_moved_side(end));

            m.plot(t,y," DF ",pdf," Velo ",true," Pause ",.02);
        }
        static function  test_DebugConfigs(versions) {
            if nargin < 1
                versions = [1 2 3 10 11];
            end
            for k = versions
                fprintf(" Testing DebugConfig Version %d\n ",k);
                m = models.muscle.Model(models.muscle.examples.Debug(k));
                [t,y] = m.simulate;
                df = m.getResidualForces(t,y);
                m.plot(t,y," DF ",df," Velo ",true);
            end
        }

    

};
}
}
}