MuscleTendonMix.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 MuscleTendonMix
  :public models.muscle.AMuscleConfig {
    public:

        Variant;

        
    public:

        ylen;

        
    public:

        MuscleTendonMix(varargin) {
            this = this@models.muscle.AMuscleConfig(varargin[:]);
            this.addOption(" Variant ",1);
            this.init;

            /* this.NeumannCoordinateSystem = 'global'; */
        }


        function  configureModel(m) {
            configureModel@models.muscle.AMuscleConfig(this, m);
            m.T= 4;
            m.dt= .01;
            m.DefaultInput= 1;

            os = m.ODESolver;
            os.RelTol= .0001;
            os.AbsTol= .05;

            mu = m.DefaultMu;
            switch this.Options.Variant
                case [1 2]
                    mu(3) = 4;
                case 3
                    mu(3) = 4;
                    m.ODESolver.AbsTol= .1;
                case 4
                    mu(2) = m.T;
                    mu(3) = 0;
                    m.ODESolver.RelTol= .1;
                    m.ODESolver.AbsTol= .1;
                case 5
                    mu(2) = m.T;
                    m.ODESolver.RelTol= .1;
                    m.ODESolver.AbsTol= .1;
                case 6
                    mu(3) = 10;
                case 7
                    m.T= 30;
                    m.dt= .05;
                    mu(2) = m.T*2/3; /*  activate 2/3rds of the time */

                    mu(3) = 0;
            end
            mu(13) = .250; /*  Pmax [MPa] */

            m.DefaultMu= mu;
        }


        function u = getInputs() {
            u[1] = this.getAlphaRamp(1,1);
        }
        function tmr = getTendonMuscleRatio(points) {
            tmr = zeros(1,size(points,2));
            switch this.Options.Variant
                case 1
                    ratios = linspace(0,1,4);
                    for k=1:4
                        tmr(points(2,:) >= k-1 .* points(2,:) < k) = ratios(k);
                    end
                case [2 5]
                    tmr(points(3,:) > .66) = 1;
                case [3 4]
                    tmr = this.TMRFun(points(2,:),points(3,:));
                    /*  Take only half - too stiff otherwise */
                    tmr = tmr/2;
                case [6 7]
                    /*  Set the middle 20% percent to 100% muscle */
                    zeroperc = .2;

                    y = points(2,:);
                    cent = this.ylen/2;
                    left = y <= cent*(1-zeroperc/2);
                    tmr(left) = 1-y(left)/(cent*(1-zeroperc/2));
                    right = y > cent + this.ylen*zeroperc/2;
                    tmr(right) = (y(right)-(this.ylen*(.5+zeroperc/2)))/...
                        (this.ylen*(.5-zeroperc/2));
            end
            /*  Use max .2 tendon - works but needs to be checked for higher
             * tendon parts (and current default parameter mu!!!!) */
            if this.Options.Variant == 4
                tmr = 2*tmr/5;
            end
        }
        function tmr = TMRFun(z,matrix<double> y) {
            f = @(z,y)min(1,max(0,-.5+.5./((z/2+.1).*(1.5*y+.08*z+.1))));
            tmr = f(z,y) + f(6-z,2-y);
        }


        function  plotTMRFun() {
            [Y,Z] = meshgrid(0:.01:2.1,0:.01:6.1);
            tmr = this.TMRFun(Z,Y);
            pm = PlotManager;
            pm.LeaveOpen= true;
            ax = pm.nextPlot(" tmr_func "," Muscle/Tendon ratio function ");
            surfc(Z,Y,tmr," Parent ",ax);
            axis(ax," equal ");
        }


        function P = getBoundaryPressure(elemidx,faceidx) {
            P = [];
            switch this.Options.Variant
                /*  Pull on front face */
                case [1 2 3]
                    if any(elemidx == [1 2 7 8]) && faceidx == 3
                        P = 1;
                    end
                case 6
                    if any(elemidx == 1:4) && faceidx == 3
                        P = 1;
                    end
            end
        }
    protected:


        function geo = getGeometry() {
            switch this.Options.Variant
                case 1
                    /*  Four cubes in a row */
                    geo = fem.geometry.RegularHex8Grid(0:1,0:4,0:1);
                case [2 5]
                    /*  Single cube with same config as reference element */
                    geo = fem.geometry.RegularHex8Grid(0:1,0:1,0:1);
                case [3 4]
                    /*  2x4x2 geometry */
                    geo = fem.geometry.RegularHex8Grid(0:2,0:2:6,0:2);
            end
            if any(this.Options.Variant == [6 7])
                yl = 10;
                this.ylen= yl;
                geo = fem.geometry.Belly(6,yl," Radius ",1," InnerRadius ",.2," Gamma ",2);
            else
                geo = geo.toCube27Node;
            end
        }


        function displ_dir = setPositionDirichletBC(displ_dir) {
            geo = this.FEM.Geometry;

            switch this.Options.Variant
                case [1 2 5]
                    /*  Fix all on left and only the y,z directions of the back right
                     * nodes */
                    displ_dir(2,geo.Elements(geo.NumElements,geo.MasterFaces(4,:))) = true;
                case 3
                    displ_dir(2,geo.Elements([5 6 11 12],geo.MasterFaces(4,:))) = true;
                    displ_dir(:,geo.Elements(5,geo.MasterFaces(4,9))) = true;
                case 4
                    /* displ_dir(2,geo.Elements([5 6 11 12],geo.MasterFaces(4,:))) = true;
                     *displ_dir(2,geo.Elements([1 2 7 8],geo.MasterFaces(3,:))) = true; */
                    displ_dir(:,geo.Elements([5 6 11 12],geo.MasterFaces(4,:))) = true;
                    displ_dir(:,geo.Elements([1 2 7 8],geo.MasterFaces(3,:))) = true;
                case 6
                    displ_dir(:,geo.Elements(13:16,geo.MasterFaces(4,:))) = true;
                case 7
                    displ_dir(:,geo.Elements(1:4,geo.MasterFaces(3,:))) = true;
                    last = size(geo.Elements,1)-3:size(geo.Elements,1);
                    displ_dir(:,geo.Elements(last,geo.MasterFaces(4,:))) = true;
            end
        }
        function anull = seta0(anull) {
            switch this.Options.Variant
                case [1 2 5 6 7]
                    /*  Direction is y */
                    anull(2,:,:) = 1;
                case [3 4]
                    /* anull(2,:,:) = 1; */
                    anull([2 3],:,:) = -1;
            end
        }

    
    public: /* ( Static ) */

        static function  test_MuscleTendonMix(variant) {
            if nargin < 1
                variant = 1:7;
            end
            for k = variant
                fprintf(" Testing MuscleTendonMix Variant %d\n ",k);
                m = models.muscle.Model(models.muscle.examples.MuscleTendonMix(" Variant ",k));
                m.simulateAndPlot;
            end
        }

    

};
}
}
}