Dynamics.m

Go to the documentation of this file.

namespace models{
namespace muscle{


/* (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 Dynamics
  :public dscomponents.ACompEvalCoreFun {
    public:

                b1cf = 5.316372204148964;
d1cf = 0.014991843974911;
        function_handle alpha = @("t")0;
        ComputeUnassembled = false;
        Sigma;
        idx_uv_bc_glob_unass;
        fDim_unass;

        num_dv_unass;

        idx_v_elems_unass;


        ForceLengthFun;
        ForceLengthFunDeriv;

        
    public:

        APExp;

        AnisoPassiveMuscle;

        AnisoPassiveMuscleDeriv;

        AnisoPassiveTendon;

        AnisoPassiveTendonDeriv;

        
    public:

        lambda_dot_pos;
        lambda_dot;


        nfibres;

        
    public: /* ( Transient ) */

        LastBCResiduals;
        nfevals;

        nJevals;

        
    public: /* ( Transient ) */

        RampTime;
    private: /* ( Transient ) */

        usemassinv;
        crossfibres = false;
    private:

        fsys;
    public:

        Dynamics(sys) {
            this = this@dscomponents.ACompEvalCoreFun(sys);

            /* % Load AP expansion
             *             d = fileparts(which('models.muscle.Dynamics'));
             *             s = load(fullfile(d,'AP'));
             *             s.kexp.Ma = s.kexp.Ma(1,:);
             *             this.APExp = s.kexp; */
        }


        function  configUpdated() {
            sys = this.fsys;
            mc = sys.Model.Config;
            if ~isempty(mc.FibreTypeWeights)
                this.nfibres= size(mc.FibreTypeWeights,2);
            end
            if ~isempty(mc)
                this.xDim= sys.NumTotalDofs;
                this.fDim= sys.NumDerivativeDofs;
                this.JSparsityPattern= this.computeSparsityPattern;

                /* % Sigma assembly matrix */
                this.precomputeUnassembledData;
            end
        }


        function  prepareSimulation(colvec<double> mu) {
            prepareSimulation@dscomponents.ACompEvalCoreFun(this, mu);

            /*  Reset counters */
            this.nfevals= 0;
            this.nJevals= 0;

            sys = this.fsys;
            mc = sys.Model.Config;
            if ~isempty(mc.Pool)
                mc.Pool.prepareSimulation(sys.Model.T,mu(4));
            end
            /*  Returns an all zero function if mu(2) is less or equal to zero! */
            this.alpha= mc.getAlphaRamp(mu(2));

            /*  Prepare force-length fun */
            mc.setForceLengthFun(this);

            /*  Muscle anisotropic passive law */
            mlfg = models.muscle.functions.MarkertLawOriginal(mu(5),mu(6));
            [this.AnisoPassiveMuscle, this.AnisoPassiveMuscleDeriv] = mlfg.getFunction;

            /*  Tendon anisotropic passive law */
            mlfg = general.functions.CubicToLinear(mu(7),mu(8));
/*              mlfg = general.functions.MarkertLaw(1.3895e+07,11.1429,1.637893706954065e+05); */
            [this.AnisoPassiveTendon, this.AnisoPassiveTendonDeriv] = mlfg.getFunction;

            /*  Get the law function handles that also take b,d as arguments.
             * Needed due to possibly inhomogeneous material.
             *[~,~,this.MarkertLawFun,this.MarkertLawFunDeriv] = mlfg.getFunction; */

            /*  Cache stuff */
            this.usemassinv= sys.UseDirectMassInversion;
            this.crossfibres= sys.HasFibres && sys.UseCrossFibreStiffness;
        }


        function fx = evaluateCoreFun(colvec<double> x,double t) {
            error(" Do not call directly; have custom evaluate method. ");
        }


        function fx = evaluateComponentSet(integer nr,colvec<double> x,double t) {
            fx = this.evaluate(x,t);
            fx = this.V(this.PointSets[nr],:)*fx;
        }
        function fx = evaluateComponentSetMulti(integer nr,matrix<double> x,rowvec<double> t,matrix<double> mu) {
            fx = this.evaluateMulti(x,t,mu);
            fx = fx(this.PointSets[nr],:);
        }
        function res = test_Jacobian(matrix<double> y,double t,colvec<double> mu) {

            if nargin < 4
                mu = this.System.Model.DefaultMu;
                if nargin < 3
                    t = 1000;
                    if nargin < 2
                        y = this.System.getX0(mu);
                    end
                end
            end

            /*  Use nonzero t to have an effect */
            res = test_Jacobian@dscomponents.ACoreFun(this, y, t, mu);
        }
        function copy = clone() {
            copy = models.muscle.Dynamics(this.System);

            /*  Call superclass clone (for deep copy) */
            copy = clone@dscomponents.ACompEvalCoreFun(this, copy);

            /*  Copy local properties
             * No local properties are to be copied here, as so far everything is done in the
             * constructor. */
        }
        function  setSystem(sys) {
            setSystem@dscomponents.ACoreFun(this, sys);
            if isa(sys," models.ReducedSecondOrderSystem ")
                this.fsys= sys.Model.FullModel.System;
            else
                this.fsys= sys;
            end
        }


        function copy = project(V,W) {
             copy = project@dscomponents.ACompEvalCoreFun(this, V, W);
        }


        
#if 0 //mtoc++: 'set.ComputeUnassembled'
function  ComputeUnassembled(value) {
            if value && ~isempty(this.fDim_unass)/* #ok */

                this.fDim= this.fDim_unass;/* #ok */

            elseif ~value && ~isempty(this.System)
                this.fDim= this.System.NumDerivativeDofs;
            end
            this.ComputeUnassembled= value;
        }

#endif


        function [SPK , SPg , SPalpha , SPLamDot ] = computeSparsityPattern();
    
    protected:

        function fx = evaluateComponents(rowvec<integer> pts,rowvec<integer> ends,unused1,unused2,matrix<double> x,unused3) {
        }
        function dfx = evaluateComponentPartialDerivatives(rowvec<integer> pts,rowvec<integer> ends,rowvec<integer> idx,rowvec<integer> deriv,rowvec<integer> self,colvec<double> x,double t,dfxsel) {

        }
    private:

        function  precomputeUnassembledData() {
            sys = this.System;
            mc = sys.Model.Config;
            geo = mc.FEM.Geometry;
            this.num_dv_unass= 3*geo.NumElements * geo.DofsPerElement;

            /*  Velocity part: x,y,z velocities */
            [i, ~] = find(mc.FEM.Sigma);
            I = [3*(i" -1)+1; 3*(i "-1)+2; 3*(i^t-1)+3];

            if numel(I) ~= this.num_dv_unass
                error(" Size mismatch. Somethings wrong with FEM Sigma ");
            end

            n = numel(I);
            S = sparse(I,1:n,ones(n,1),3*geo.NumNodes,n);

            /*  Take out nodes with dirichlet BC on output side */
            S(sys.idx_v_bc_local,:) = [];
            /*  Find corresponding unassembled dofs that would be ignored
             * (due to dirichlet velocity values, pressure dirichlet not
             * implemented) */
            bc_unass = find(sum(S,1) == 0);
            /*  Remove them, too. The unassembled evaluation also removes the
             * corresponding entries of the unassembled vector. */
            S(:,bc_unass) = [];
            this.idx_uv_bc_glob_unass= bc_unass; /* [sys.idx_u_bc_glob' bc_unass]; */

            this.Sigma= S;
            this.fDim_unass= size(S,2);

            /*  Create boolean array that determines which unassembled dofs
             * belong to which element
             * dK dofs */
            hlp = repmat(1:geo.NumElements, 3*geo.DofsPerElement,1);
            hlp = hlp(:);
            hlp(bc_unass) = [];
            ass = false(geo.NumElements,length(hlp));
            for k = 1:geo.NumElements
                ass(k,:) = hlp == k;
            end
            this.idx_v_elems_unass= ass;
        }

};
}
}