Activation.m

Go to the documentation of this file.

namespace models{
namespace muscle{
namespace tests{


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

        Activation(version) {
            if nargin < 1 
                version = 1;
            end
            /*  Creates a Debug simple muscle model configuration. */
            this = this@models.muscle.AMuscleConfig;
            this.addOption(" Version ",version);
            this.init;

            /*  Use classical FL function with fl(1)=1 */
            this.Options.FL= 3;
        }


        function  configureModel(m) {
            configureModel@models.muscle.AMuscleConfig(this, m);
/*              sys = m.System;
             *f = sys.f; */
            m.T= 15;
            m.dt= 1;
            /*  Activate over 50ms */
            m.DefaultMu(2) = 10;
            m.ODESolver.RelTol= 1e-10;
            m.ODESolver.AbsTol= 1e-10;
        }

    
    protected:


        function geo = getGeometry() {
            geo = fem.geometry.RegularHex8Grid([0 2],[0 1],[0 1]);
            if this.Options.Version == 2
                geo = fem.geometry.RegularHex8Grid([0:.5:1],[0:1:2],[0:.5:1]);
            end
            geo = geo.toCube27Node;
        }
        function displ_dir = setPositionDirichletBC(displ_dir) {
            geo = this.FEM.Geometry;
            /*  Always fix back side */
            switch this.Options.Version
                case 1
                    displ_dir(:,geo.Elements(1,geo.MasterFaces(1:2,:))) = true;
                case 2
                    displ_dir(:,geo.Elements(1:4,geo.MasterFaces(1,:))) = true;
                    displ_dir(:,geo.Elements(5:8,geo.MasterFaces(2,:))) = true;
/*                  otherwise
 *                     displ_dir(:,geo.Elements(1,geo.MasterFaces(1,:))) = true;
                     *displ_dir(:,geo.Elements(1:4,geo.MasterFaces(1,:))) = true; */
            end
        }


        function anull = seta0(anull) {
            anull(1,:,:) = 1;
        }

    
    public: /* ( Static ) */

        static function res = test_Activation() {
            for version = 1:2
                a = models.muscle.tests.Activation(version);
                m = a.createModel;
                if version == 1
                    f1 = m.getPositionDirichletBCFaceIdx(1,1);
                    f2 = m.getPositionDirichletBCFaceIdx(1,2);
                    m.setGaussIntegrationRule(4);
                else
                    f1 = m.getPositionDirichletBCFaceIdx(1:4,1);
                    f2 = m.getPositionDirichletBCFaceIdx(5:8,2);
                end

                mu = m.DefaultMu;
                maxp = [1 rand(1,5)*10];
                res = true;
                for k=1:length(maxp)
                   mu(13) = maxp(k); 
                   [t,y] = m.simulate(mu);
                   pdf = m.getResidualForces(t,y);
                   /*  Check overall zero force
                    *max(abs(sum(pdf))) */
                   res = res && all(abs(sum(pdf)) < 1e-8);

                   force_l = sum(pdf(f1,end),1);
                   force_r = sum(pdf(f2,end),1);
                   fprintf(" P_{max}=%g => Force left %gN, right %gN\n ",...
                       maxp(k),force_l,force_r);
                   /* norm(maxp(k)-force_l) */
                   res = res && norm(version*maxp(k)-force_l) < 1e-8;
                end
            end
        }
};
}
}
}