MusclePlotter.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 MusclePlotter
  :public handle {
    public:

        GeoView = "[46 30]";

        DefaultArgs = {""};

    
    public:

        System;

        Config;

        
    protected:

        plotdata;

        
    public:

        MusclePlotter(sys) {
            this.System= sys;
            this.Config= sys.Model.Config;
            this.plotdata= this.initPlotData;
        }


        function [pm , h ] = plot(double t,y_dofs,varargin) {
            opts = this.parsePlotArgs(varargin);

            mc = this.Config;

            if isempty(opts.PM)
                if ~isempty(mc.Pool) && opts.Pool
                    pm = PlotManager(false,2,1);
                else
                    pm = PlotManager;
                end
                pm.LeaveOpen= true;
            else
                pm = opts.PM;
            end

            [pd, t] = this.updatePlotData(this.plotdata, opts, t, y_dofs);
            this.plotdata= pd;

            if opts.Vid
                [p,n,e] = fileparts(opts.Vid);
                if ~isempty(p) && exist(p," file ") ~= 7
                    error(" Directory %s not found ",p);
                end
                if isempty(p)
                    if isa(mc," models.muscle.AExperimentModelConfig ")
                        p = mc.OutputDir;
                    else
                        p = pwd;
                    end
                end
                if isempty(n)
                    n = " output ";
                end
                if isempty(e)
                    e = " .avi ";
                end
                avifile = fullfile(p,[n e]);
                vw = VideoWriter(avifile);
                vw.FrameRate= 25;
                vw.open;
            end

            /* % Loop over time */
            if ~isempty(mc.Pool) && opts.Pool
                hf = pm.nextPlot(" force "," Activation force "," t [ms] "," alpha ");
                axis(hf,[0 t(end) 0 1]);
                hold(hf," on ");
            end

            h = pm.nextPlot(" geo ",sprintf(" Deformation at t=%g ",t(end))," x [mm] "," y [mm] ");
            zlabel(h," z [mm] ");
            axis(h, pd.geo_plotbox);
            daspect([1 1 1]);
            view(h, this.GeoView);
            hold(h," on ");

            for ts = 1:length(t)
                /*  Quit if figure has been closed */
                if ~ishandle(h)
                    if opts.Vid
                        vw.close;
                    end
                    break;
                end
                this.plotGeometry(h, t(ts), pd.yfull(:,ts), ts, opts);

                if ~isempty(mc.Pool) && opts.Pool
                    times = t(1:ts);
                    alpha = mc.Pool.getActivation(times);
                    walpha = mc.FibreTypeWeights(1,:,1) * alpha;
                    cla(hf);
                    plot(hf,times,alpha);
                    plot(hf,times,walpha," LineWidth ",2);
                end

/*                  % Quit if figure has been closed (extra check - runtime
 *                 % thing)
 *                 if ~ishandle(h)
 *                     if opts.Vid
 *                         vw.close;
 *                     end
 *                     break;
 *                 end */

                if ~isempty(opts.Vid) && ishandle(h)
                    vw.writeVideo(getframe(gcf));
                else
                    pause(opts.Pause);
                end
            end

            if opts.Vid
                vw.close;
            end

            if isempty(opts.PM)
                pm.done;
            end
        }

    
    protected:

        function  plotGeometry(h,double t,y_dofs,ts,opts) {
            casecap = " Deformation ";
            if opts.Stretch
                casecap = " Lambda-stretches ";
            end
            title(h,sprintf(" %s at t=%g ",casecap,t));

            pd = this.plotdata;
            sys = this.System;
            nelem = length(opts.Elems);

            u = reshape(y_dofs(1:pd.vstart-1),3,[]);
            v = reshape(y_dofs(pd.vstart:pd.pstart-1),3,[]);
            cla(h);

            geo = this.Config.FEM.Geometry;
            if opts.Skel
                e = pd.e;
                plot3(h,u(1,pd.no_bc),u(2,pd.no_bc),u(3,pd.no_bc)," r. "," MarkerSize ",14);
                for k=1:size(e,1)
                    plot3(h,u(1,[e(k,1) e(k,2)]),u(2,[e(k,1) e(k,2)]),u(3,[e(k,1) e(k,2)])," r ");
                end 
            elseif ~opts.Stretch
                p = patch(" Faces ",geo.PatchFaces," Vertices ",u^t,...
                    " Parent ",h," FaceAlpha ",.3);
                if sys.HasTendons
                    set(p," FaceVertexCData ",pd.vertexcol," FaceColor ",...
                        " interp "," FaceAlpha ",.5," EdgeColor "," interp ");
                else
                    set(p," EdgeColor ",.8*pd.musclecol," FaceColor ",pd.musclecol);
                end
            end

            /*  Velocities */
            if opts.Velo
                quiver3(h,u(1,:),u(2,:),u(3,:),v(1,:),v(2,:),v(3,:),1," b ", " MarkerSize ",14);
            end

            /* % Dirichlet conditions
             * Displacement */
            u3dir = u(:,pd.bc_dir_3pos_applies);
            plot3(h,u3dir(1,:),u3dir(2,:),u3dir(3,:)," . "," MarkerSize ",20," Color ",[0 0 0]);
            u2dir = u(:,pd.bc_dir_2pos_applies);
            plot3(h,u2dir(1,:),u2dir(2,:),u2dir(3,:)," . "," MarkerSize ",20," Color ",[.5 .5 .5]);
            u1dir = u(:,pd.bc_dir_1pos_applies);
            plot3(h,u1dir(1,:),u1dir(2,:),u1dir(3,:)," . "," MarkerSize ",20," Color ",[.7 .7 .7]);

            /*  Velocity */
            uvdir = u(:,pd.bool_expl_v_bc_nodes_applies);
            plot3(h,uvdir(1,:),uvdir(2,:),uvdir(3,:)," g. "," MarkerSize ",20);

            /* % Dirichlet Forces */
            if ~isempty(pd.DF)
                residuals = pd.residuals;
                have_residuals = pd.have_residuals;
                udir = u(:,have_residuals);
/*                  residuals(pd.residuals_pos) = pd.DF(pd.sortidx,ts)/pd.maxdfval; */
                residuals(pd.residuals_pos) = pd.DF(:,ts)/pd.maxdfval;
                quiver3(h,udir(1,:),udir(2,:),udir(3,:),...
                    residuals(1,have_residuals),residuals(2,have_residuals),...
                    residuals(3,have_residuals),1," k ", " MarkerSize ",14);
            end

            /* % Neumann condition forces */
            if opts.Forces
                /*  Plot force vectors at nodes
 *                 uforce = u(:,pd.forces_apply);
 *                 quiver3(h,uforce(1,:),uforce(2,:),uforce(3,:),...
 *                     pd.forces(1,:),pd.forces(2,:),pd.forces(3,:),0,'Color',[.8 .8 1]);
                 * Plot force vectors at face centers */
                for k=1:pd.numfaceswithforce
                    masterfacenodeidx = geo.MasterFaces(pd.force_elem_face_idx(2,k),:);
                    facenodeidx = geo.Elements(pd.force_elem_face_idx(1,k),masterfacenodeidx);

                    facecenter = mean(u(:,facenodeidx),2);
                    mf = pd.meanforces(:,k)*this.System.mu(3);
                    quiver3(h,facecenter(1),facecenter(2),facecenter(3),...
                        mf(1),mf(2),mf(3),0," LineWidth ",2," Color "," b "," MaxHeadSize ",1);

                    if ~isempty(pd.NF)
                        pd.residual_neumann_forces(sys.idx_neumann_bc_glob) = pd.NF(:,ts);
                        meanforce = mean(pd.residual_neumann_forces(:,facenodeidx),2);
                        quiver3(h,facecenter(1),facecenter(2),facecenter(3),...
                        meanforce(1),meanforce(2),meanforce(3),0," LineWidth ",2," Color "," k "," MaxHeadSize ",1);
                    end
                end
            end

            /* % Pressure */
            if opts.Pressure
                p = y_dofs(pd.pstart:sys.num_uvp_glob);
                pn = sys.idx_p_to_u_nodes;
                pneg = p<0;
                /*  Negative pressures */
                if any(pneg)
                    scatter3(h,u(1,pn(pneg)),u(2,pn(pneg)),u(3,pn(pneg)),...
                        -p(pneg)*10+1,[188 188 255]/255);/* 
                end
                 * Positive pressures */

                if any(~pneg)
                    scatter3(h,u(1,pn(~pneg)),u(2,pn(~pneg)),u(3,pn(~pneg)),...
                        p(~pneg)*10+1,[255 188 188]/255); /* [244 149 25]/255 */

                end
            end

            /* % a0 fibres & gp values */
            fibres = sys.HasFibres && opts.Fibres;
            doI1 = any(opts.Invariants == 1);
            doLam = opts.Lambdas;
            doMR = opts.MR;
            doTMR = opts.MTRatio;
            doStr = opts.Stretch;
            if fibres || ~isempty(opts.Invariants) || doLam || doMR || doTMR || doStr
                for elemidx = 1:nelem
                    m = opts.Elems(elemidx);
                    u = y_dofs(1:pd.vstart-1);
                    u = u(sys.idx_u_elems_local(:,:,m));
                    gps = u*pd.Ngp;

                    if doStr
                        scatter3(h,gps(1,:),gps(2,:),gps(3,:),1,pd.stretchcol(:,:,elemidx,ts)" , "." , "SizeData^t,30);
                    else
                        if fibres
                            anull = u*sys.dNa0(:,:,m);
                            quiver3(gps(1,:),gps(2,:),gps(3,:),anull(1,:),anull(2,:),anull(3,:),.5," . "," Color "," w ");
                        end

                        /* % tendon ratio */
                        if ~isempty(pd.gpcol)
                            scatter3(h,gps(1,:),gps(2,:),gps(3,:),1,pd.gpcol(:,:,m)," . "," SizeData ",30);
                        end

                        /* % Invariants */
                        if doI1
                            values = sprintfc(" I_1=%3.3g ",pd.I1(:,elemidx,ts));
                            text(gps(1,:),gps(2,:),gps(3,:),values," Parent ",h)
                        end
                        /* % Lambdas */
                        if doLam
                            /* values = sprintfc('\\lambda=%3.3g',pd.lambdas(:,m,ts)); */
                            values = sprintfc(" %3.4g ",pd.lambdas(:,elemidx,ts));
                            text(gps(1,:),gps(2,:),gps(3,:),values," Parent ",h)
                        end
                        /* % Mooney-Rivlin tensors */
                        if doMR
                            values = cellfun(@(v)sprintf(" %g  ",diag(v)), pd.MR(:,elemidx,ts)," UniformOutput ",false);
                            text(gps(1,:),gps(2,:),gps(3,:),values," Parent ",h)
                        end

                        /* % Tendon-Muscle ratio */
                        if doTMR
                            values = sprintfc(" %3.4g ",sys.MuscleTendonRatioGP(:,m,ts));
                            /* values = sprintfc('%3.4g',diags{:}); */
                            text(gps(1,:),gps(2,:),gps(3,:),values," Parent ",h)
                        end
                    end
                end
                if doStr
                    colormap(pd.stretchcmap);
                    colorbar(" YTickLabel ",pd.stretchcmaplbl);
                end
            end
        }


        function [pd , doublet , matrix<double>y ] = updatePlotData(pd,opts,double t,matrix<double> y) {
            mc = this.Config;
            sys = this.System;
            dfem = mc.FEM;
            geo = dfem.Geometry;

            pd.DF= opts.DF;
            pd.NF= opts.NF;
            /* % "Speedup" factor for faster plots */
            if ~isempty(opts.F)
                t = t(1:opts.F:end);
                y = y(:,1:opts.F:end);
                if ~isempty(pd.DF)
                    pd.DF= pd.DF(:,1:opts.F:end);
                end
                if ~isempty(pd.NF)
                    pd.NF= pd.NF(:,1:opts.F:end);
                end
            end

            /* % Re-add the dirichlet nodes for geometry plotting */
            pd.yfull= sys.includeDirichletValues(t, y);
            pd.geo_plotbox= this.getPlotBox(pd.yfull);

            /* % Dirichlet plotting */
            if ~isempty(opts.DF)
                pd.maxdfval= max(abs(opts.DF(:)))/10;
            end

            /* % Forces plotting */
            if opts.Forces
                /*  Get forces on each node in x,y,z directions
                 * This is where the connection between plane index and
                 * x,y,z coordinate is "restored" */
                forces = zeros(size(sys.bool_u_bc_nodes));
                forces(sys.idx_neumann_bc_glob) = sys.bc_neum_forces_val;

                /*  Forces at face centers */
                force_elem_face_idx = geo.Faces(:,sys.FacesWithForce);
                pd.numfaceswithforce= size(force_elem_face_idx,2);
                meanforces = zeros(3,pd.numfaceswithforce);
                for k=1:pd.numfaceswithforce
                    masterfacenodeidx = geo.MasterFaces(force_elem_face_idx(2,k),:);
                    facenodeidx = geo.Elements(force_elem_face_idx(1,k),masterfacenodeidx);
                    meanforces(:,k) = mean(forces(:,facenodeidx),2);
                end
                pd.meanforces= meanforces;
                pd.force_elem_face_idx= force_elem_face_idx;

                /*  Also save forces at x,y,z nodes for plotting */
                pd.forces_apply= sum(abs(forces),1) ~= 0;
                pd.forces= forces(:,pd.forces_apply);

                if ~isempty(opts.NF)
                    pd.residual_neumann_forces= zeros(size(sys.bool_u_bc_nodes));
                end
            end

            /* % Skeleton plotting */
            if ~opts.Skel
        /*                  light('Position',[1 1 1],'Style','infinite','Parent',h); */
                pd.musclecol= [0.854688, 0.201563, 0.217188];
            end

            /* % Show invariants or lambda stretch values */
            if ~isempty(opts.Invariants) || opts.Lambdas || opts.MR || opts.Stretch
                nt = length(t);
                nelem = length(opts.Elems);
                doI1 = any(opts.Invariants == 1);
                doLam = opts.Lambdas;
                doMR = opts.MR;
                doStr = opts.Stretch;
                if doI1 || doMR
                    pd.I1= zeros(dfem.GaussPointsPerElem,nelem,nt);
                end
                if doLam || doStr
                    pd.lambdas= zeros(dfem.GaussPointsPerElem,nelem,nt);
                    if doStr
                        pd.stretchcol= zeros(3,dfem.GaussPointsPerElem,nelem,nt);
                    end
                end
                if doMR
                    pd.MR= cell(dfem.GaussPointsPerElem,nelem,nt);
                    c10 = sys.MuscleTendonParamc10;
                    c01 = sys.MuscleTendonParamc01;
                end
                num_gp = dfem.GaussPointsPerElem;
                for elemidx = 1:nelem
                    m = opts.Elems(elemidx);
                    for gp = 1:num_gp
                        pos = 3*(gp-1)+1:3*gp;
                        dtn = dfem.transgrad(:,pos,m);        
                        elemidx_u = sys.idx_u_elems_local(:,:,m); 
                        for ts = 1:nt
                            /*  Deformation gradient */
                            yts = pd.yfull(:,ts) ;
                            F = yts(elemidx_u) * dtn;
                            C = F^t*F;
                            if doI1 || doMR
                                pd.I1(gp,elemidx,ts) = C(1,1)+C(2,2)+C(3,3);
                            end
                            if doLam || doStr
                                fibres = sys.a0Base(:,:,(m-1)*num_gp + gp);
                                pd.lambdas(gp,elemidx,ts) = norm(F*fibres(:,1));
                            end
                            if doMR
                                pd.MR[gp,elemidx,ts] = sys.MooneyRivlinICConst(gp,m)*eye(3) ...
                                    + 2*(c10(gp,m) + pd.I1(gp,elemidx,ts)*c01(gp,m))*F ...
                                    - 2*c01(gp,m)*F*C;
                            end
                        end
                    end
                end
                if doStr
                    Mlam = max(pd.lambdas(:));
                    mlam = min(pd.lambdas(:));
                    center = (1-mlam)/(Mlam-mlam);
                    if abs(Mlam-mlam) < 1e-8
                        /*  Green for lambda=1 everywhere */
                        pd.stretchcol(2,:,:,:) = 1;
                    else
                        colfun = @(x)[(x>center).*((x-center)/(1-center))
                            (x<center).*x/center + (x>=center).*(1-(x-center)/(1-center));
                            (x<center).*(1-x/center);];
                        for ts = 1:nt
                            for elemidx = 1:nelem
                                normlam = (pd.lambdas(:,elemidx,ts)-mlam) / (Mlam-mlam);
                                pd.stretchcol(:,:,elemidx,ts) = colfun(normlam^t);
                            end
                        end
                    end
                    pd.stretchcmap= colfun(linspace(0,1,125))^t;
                    pd.stretchcmaplbl= sprintfc(" %3.3g ",linspace(mlam,Mlam,12));
                end
            end
        }


        function opts = parsePlotArgs(args) {
            
            i.KeepUnmatched= true;
            i.addParamValue('Vid',[],@(v)~isempty(v));
            i.addParamValue('Forces',false,@(v)islogical(v));
            i.addParamValue('Velo',false,@(v)islogical(v));
            i.addParamValue('Pressure',false,@(v)islogical(v));
            i.addParamValue('Fibres',true,@(v)islogical(v));
            i.addParamValue('Skel',false,@(v)islogical(v));
            i.addParamValue('Pool',false,@(v)islogical(v));
            i.addParamValue('PM',[],@(v)isa(v,'PlotManager'));
            i.addParamValue('DF',[]);
            i.addParamValue('NF',[]);
            i.addParamValue('F',[]);
            i.addParamValue('Elems',1:this.Config.FEM.Geometry.NumElements);
            i.addParamValue('Invariants',[]);
            i.addParamValue('Lambdas',false,@(v)islogical(v));
            i.addParamValue('MR',false);
            i.addParamValue('MTRatio',false,@(v)islogical(v));
            i.addParamValue('Pause',.01);
            i.addParamValue('Stretch',false,@(v)islogical(v));

            args = [this.DefaultArgs args];
            i.parse(args[:]);
            opts = i.Results;
            if ~isempty(opts.NF)
                opts.Forces= true;
            end
        }


        function pd = initPlotData() {
            sys = this.System;
            pd = struct;
            hlp = sum(sys.bool_u_bc_nodes,1);
            pd.bc_dir_3pos_applies= hlp == 3;
            pd.bc_dir_2pos_applies= hlp == 2;
            pd.bc_dir_1pos_applies= hlp == 1;
            pd.bc_dir_pos_applies= hlp >= 1; 
            pd.bool_expl_v_bc_nodes_applies= sum(sys.bool_expl_v_bc_nodes,1) >= 1;
            pd.no_bc= ~pd.bc_dir_pos_applies & ~pd.bool_expl_v_bc_nodes_applies;

            mc = this.Config;
            dfem = mc.FEM;
            geo = dfem.Geometry;
            num_u_glob = geo.NumNodes * 3;
            pd.vstart= num_u_glob+1;
            pd.pstart= 2*num_u_glob+1;
            pd.e= geo.Edges;

            /* % Dirichlet forces */
            pd.have_residuals= pd.bc_dir_pos_applies | pd.bool_expl_v_bc_nodes_applies;
            /*  By sorting and combining the pos/velo Dir BC, the
             * plotting of mixed BCs on one node is plotted correctly. */
            pd.residuals_pos= sys.bool_u_bc_nodes | sys.bool_expl_v_bc_nodes;
/*              [~, pd.sortidx] = sort([sys.idx_u_bc_glob; sys.idx_expl_v_bc_glob]);
             * Preallocate the residuals matrix */
            pd.residuals= zeros(size(pd.residuals_pos));

            /* % Fibres
             *if sys.HasFibres || ~isempty(opts.Invariants) */
                pd.Ngp= dfem.N(dfem.GaussPoints);
            /* end */

            /*  Set muscle color either way */
            pd.musclecol= [0.854688, 0.201563, 0.217188];
            pd.gpcol= [];
            /* % Muscle/Tendon patch vertex colors */
            if sys.HasTendons
                tendoncol = [1, .8, .6]; /* Micha [.9 .7 .5]
                 * Vertex coloring */

                pd.vertexcol= ones(geo.NumNodes,1)*pd.musclecol ...
                    + sys.MuscleTendonRatioNodes^t*(tendoncol-pd.musclecol);
                /*  Gauss point coloring */
                tmr = sys.MuscleTendonRatioGP;
                pd.gpcol(:,:,1) = pd.musclecol(1) + tmr*(tendoncol(1)-pd.musclecol(1));
                pd.gpcol(:,:,2) = pd.musclecol(2) + tmr*(tendoncol(2)-pd.musclecol(2));
                pd.gpcol(:,:,3) = pd.musclecol(3) + tmr*(tendoncol(3)-pd.musclecol(3));
                pd.gpcol= permute(pd.gpcol,[1 3 2]);/* *.8; */

            end
        }

        function box = getPlotBox(uvw) {
            xpos = 1:3:this.Config.FEM.Geometry.NumNodes*3;
            box = [min(min(uvw(xpos,:))) max(max(uvw(xpos,:)))...
                min(min(uvw(xpos+1,:))) max(max(uvw(xpos+1,:)))...
                min(min(uvw(xpos+2,:))) max(max(uvw(xpos+2,:)))];
            diam = diff(box)*.1;
            box = box + [-diam(1) diam(1) -diam(3) diam(3) -diam(5) diam(5)];
        }
};
}
}