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

        NeumannPressure(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;
        }


        function  configureModel(m) {
            configureModel@models.muscle.AMuscleConfig(this, m);
/*              sys = m.System;
             *f = sys.f; */
            m.T= 100;
            m.dt= 1;
            /*  Kleine Visko für schnelle einstellung des gleichgewichts */
            m.DefaultMu(1) = .1;/* 1e-4;
             * "Passende" Mooney-Rivlin parameter so dass das element nicht
             * allzu sehr degeneriert */

            m.DefaultMu([9 10]) = .1;
            /*  Genau lösen! */
            m.ODESolver.RelTol= 1e-7;
            m.ODESolver.AbsTol= 1e-7;
            m.DefaultInput= 1;
        }


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

    
    protected:


        function geo = getGeometry() {
            geo = fem.geometry.RegularHex8Grid([0 1],[0 1],[0 1]);
            geo = geo.toCube27Node;
        }
        function displ_dir = setPositionDirichletBC(displ_dir) {
            geo = this.FEM.Geometry;
            displ_dir(1,geo.Elements(1,geo.MasterFaces(1,:))) = true;
            displ_dir(:,geo.Elements(1,geo.MasterFaces(1,5))) = true;
        }


        function anull = seta0(anull) {
            if this.Options.Version == 2
                anull(1,:,:) = 1;
            end
        }

    
    public: /* ( Static ) */

        static function res = test_NeumannPressure() {
            for version = 1:2
                a = models.muscle.tests.NeumannPressure(version);
                m = a.createModel;
                if version == 1
                    m.setGaussIntegrationRule(4);
                end

                B = m.System.B.getMatrix(1);
                res = norm(sum(B) - 1) < 1e-10;
                extfidx = find(B);

                f1 = m.getPositionDirichletBCFaceIdx(1,1);
                mu = m.DefaultMu;
                maxf = [.1 rand(1,5)*.1];
                for k=1:length(maxf)
                   mu(3) = maxf(k); 
                   [t,y] = m.simulate(mu,1);
                   pdf = m.getResidualForces(t,y);
/*                     m.plot(t,y,'DF',pdf,'MR',true,'Pause',2); %,'Lambdas',true
 *                    m.plot(t,y,'F',2); */

                   force_l = sum(pdf(f1,:),1);
/*                     figure;
 *                    plot(t,force_l); */

                   /*  See that the forces at the pulled side match (at the end) */
                   fyend = m.System.f.evaluate(y(:,end),t(end));
                   res = res && norm(maxf(k)*B(extfidx) + fyend(extfidx)) < 1e-10;

                   fprintf(" Ext force on side %g[N] => Force left %gN\n ",...
                       maxf(k),force_l(end));
                   res = res && norm(maxf(k)-force_l(end)) < 1e-8;
                end
            end
        }


    

};
}
}
}