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

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

            this.VelocityBCTimeFun= general.functions.ConstantUntil(20);
        }
        function  configureModel(m) {
            configureModel@models.muscle.AMuscleConfig(this, m);

            m.T= 40;
            m.dt= .5;
            m.ODESolver.RelTol= 1e-6;
            m.ODESolver.AbsTol= 1e-3;
            m.DefaultMu(2) = 0;
        }


        function tmr = getTendonMuscleRatio(points) {
            tmr = zeros(1,size(points,2));
/*              tmr(:)=1; */
            switch this.Options.Flip
                case 3
                    tmr(points(3,:)<.5 & points(1,:)>.5) = 1;
                    tmr(points(3,:)>2.5 & points(1,:)<.5) = 1;
                case 2
                    tmr(points(2,:)<.5 & points(3,:)>.5) = 1;
                    tmr(points(2,:)>2.5 & points(3,:)<.5) = 1;
                otherwise
                    tmr(points(1,:)<.5 & points(3,:)>.5) = 1;
                    tmr(points(1,:)>2.5 & points(3,:)<.5) = 1;
            end
        }
    protected:


        function geo = getGeometry() {
            xr = linspace(0,3,this.Options.GeoNr+1);
            switch this.Options.Flip
                case 3
                    geo = fem.geometry.RegularHex8Grid([0 1],[0 1], xr);
                case 2
                    geo = fem.geometry.RegularHex8Grid([0 1],xr,[0 1]);
                otherwise
                    geo = fem.geometry.RegularHex8Grid(xr,[0 1],[0 1]);
            end
            geo = geo.toCube27Node;
        }
        function displ_dir = setPositionDirichletBC(displ_dir) {
            geo = this.FEM.Geometry;
            switch this.Options.Flip
                case 3
                    displ_dir(:,geo.Elements(1,geo.MasterFaces(5,:))) = true;
                    /* displ_dir([2 3],geo.Elements(this.Options.GeoNr,geo.MasterFaces(6,:))) = true; */
                case 2
                    displ_dir(:,geo.Elements(1,geo.MasterFaces(3,:))) = true;
                    /* displ_dir([2 3],geo.Elements(this.Options.GeoNr,geo.MasterFaces(4,:))) = true; */
                otherwise
                    displ_dir(:,geo.Elements(1,geo.MasterFaces(1,:))) = true;
                    /* displ_dir([2 3],geo.Elements(this.Options.GeoNr,geo.MasterFaces(2,:))) = true; */
            end

        }


        function [velo_dir , velo_dir_val ] = setVelocityDirichletBC(velo_dir,velo_dir_val) {
            geo = this.FEM.Geometry;
            switch this.Options.Flip
                case 3
                    velo_dir(3,geo.Elements(this.Options.GeoNr,geo.MasterFaces(6,:))) = true;
                case 2
                    velo_dir(2,geo.Elements(this.Options.GeoNr,geo.MasterFaces(4,:))) = true;
                otherwise
                    velo_dir(1,geo.Elements(this.Options.GeoNr,geo.MasterFaces(2,:))) = true;
            end
            velo_dir_val(velo_dir) = .02;
        }
        function anull = seta0(anull) {
            anull(this.Options.Flip,:,:) = 1;
        }

    
    public: /* ( Static ) */

        static function res = test_SubElem(gn) {
            if nargin < 1
                gn = 3;
            end
            /* %  */
            v = zeros(3,1);
            o = ones(3,1);
            for fl = 1:3
                c = models.muscle.examples.SubElemInhomogMaterial(" GeoNr ",gn," Flip ",fl);
                m = c.createModel;
                m.setGaussIntegrationRule(6);
                [t,y] = m.simulate;
                df = m.getResidualForces(t,y);
                idx = m.getPositionDirichletBCFaceIdx(1,1+(fl-1)*2);
                force = sum(df(idx,end));
                fprintf(" Flip %d - Forces at pos dirichlet side at end: %g\n ",fl,force);
                v(fl) = force;
                m.plot(t,y," DF ",df," F ",2);
            end
            err = tril(v*o" -(v*o ")^t)
            rel1 = err ./ (v*o^t)
            rel2 = err ./ (v*o" ) "
            res = all(abs(rel1(:)) < 1e-4) && all(abs(rel2(:)) < 1e-4)
        }

    

};
}
}
}