evaluate.m

Go to the documentation of this file.

#include "Constraint.m"
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.
 */

function g = models.muscle.Constraint.evaluate(uvwdof,double t) {

    sys = this.fsys;
    isproj = ~isempty(this.V);
    /*  If we evaluate inside a projected (reduced) model, reconstruct */
    if isproj
        uvwdof = this.V*uvwdof;
    end

    /*  This should be more correct
 *     unassembled = ~isproj && this.ComputeUnassembled;
 *     unassembled = this.ComputeUnassembled; */

    /* % Include dirichlet values to state vector */
    uvwcomplete = sys.includeDirichletValues(t, uvwdof);

    /* % Evaluate g(u) */
    sys = this.fsys;
    mc = sys.Model.Config;
    fe_pos = mc.FEM;
    geo = fe_pos.Geometry;
    fe_press = mc.PressureFEM;
    pgeo = fe_press.Geometry;
    elem_idx_u_glob = sys.idx_u_elems_local;
    elem_idx_p_glob = sys.idx_p_elems_local;
    unassembled = this.ComputeUnassembled;
    dofsperelem_p = pgeo.DofsPerElement;
    num_gp = fe_pos.GaussPointsPerElem;
    num_elements = geo.NumElements;

    /*  Init result vector dvw */
    if unassembled
        g = zeros(this.fDim_unass,1);
    else
        g = zeros(pgeo.NumNodes,1);
    end

    for m = 1:num_elements
        elemidx_u = elem_idx_u_glob(:,:,m); /*  1:num_u_glob is all u */

        elemidx_p = elem_idx_p_glob(:,m);

        u = uvwcomplete(elemidx_u);

        integrand_p = zeros(dofsperelem_p,1);
        for gp = 1:num_gp
            pos = 3*(gp-1)+1:3*gp;
            dtn = fe_pos.transgrad(:,pos,m);

            /*  Deformation gradient */
            F = u * dtn;

            /* % Assembly part I - sum up contributions from gauss points */
            weight = fe_pos.GaussWeights(gp) * fe_pos.elem_detjac(m,gp);

            integrand_p = integrand_p + weight * (det(F)-1) * fe_press.Ngp(:,gp,m);
        end /*  end of gauss point loop */


        /* % Assembly part II - sum up contributions of elements
         * Unassembled or assembled? */
        if unassembled
            pos = (1:dofsperelem_p) + (m-1) * dofsperelem_p;
            g(pos) = integrand_p(:);
        else
            elemidx_p_out = elemidx_p;
            g(elemidx_p_out) = g(elemidx_p_out) + integrand_p;
        end
    end /*  end of element loop */

end
}
};
};