DLTNonlinearCoreFun.m

Go to the documentation of this file.

namespace models{
namespace beam{


/* (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 DLTNonlinearCoreFun
  :public dscomponents.ACoreFun {
    private:

        Kx;

        currentx = "[]";

        dKx;

        
    public:

        DLTNonlinearCoreFun(system) {
            this = this@dscomponents.ACoreFun(system);

            free = this.System.Model.free;
            K0 = this.System.K0(free,free);
            this.fDim= size(K0,1);
            this.xDim= size(K0,2);
            JP = logical(K0);
            this.JSparsityPattern= [JP 0*JP];
        }


        function fx = evaluateCoreFun(colvec<double> x,double t) {
            this.updateB(x);
            /*  Funktion f auswerten */
            fx = -this.Kx;
        }
        function J = getStateJacobian(colvec<double> x,double t,colvec<double> mu) {
            this.updateB(x);
            /*  Jacobi-Matrix ausgeben */
            J = -this.dKx;
        }
        function pr = project(V,W) {
            error(" not implemented ")
        }

    
    private:

        function  updateB(colvec<double> x) {
            m = this.System.Model;
            /*  Updates the x, J and R big versions (if needed) */
            if isempty(this.currentx) || (norm(this.currentx - x) > 1e-8)

                /*  Only pass the spatial part of x to the weak form! */
                [Kx_full, dKx_full] = this.weak_form(x(1:length(m.free)));

                /*  K(x) */
                this.Kx= Kx_full(m.free);

                /*  \partial K|\partial x
                 * dKx = [0 0; dK 0] */
                this.dKx= dKx_full(m.free, m.free);
                /*  Augment to represent zero velocity dof dependency */
                this.dKx= [this.dKx zeros(size(this.dKx))];

                this.currentx= x;
            end
        }


        function [R , K ] = weak_form(u) {
            m = this.System.Model;
            data = m.data;
            /*  Full dimension */
            fdim = 7*data.num_knots;

            /*  Tangentiale Steifigkeitsmatrix aufstellen und schwache Form mit der aktuellen Verschiebung auswerten */
            el = m.Elements;

            nel = length(el);
            nK = zeros(nel,1); nR = nK;
            for k=1:nel
                nR(k) = length(el[k].getGlobalIndices);
            end

            i = zeros(sum(nR.^2),1); j = i; K = i;
            R = zeros(sum(nR),1); ir = R;

            u_full = zeros(fdim,1);
            /*  load dirichlet values */
            u_full(m.dir_u) = m.u_dir;
            /*  set current values, only use x */
            u_full(m.free) = u;
            Roff = 0;
            Koff = 0;
            for k=1:length(el)
                index_glob = el[k].getGlobalIndices;
                [K_lok, R_lok] = el[k].getLocalTangentials(u_full(index_glob));

                /*  Unklar welche alternative schneller ist; muss man in
                 * sp�teren tests mal profilen. Idealerweise m�ssen die i,j
                 * etc vektoren auch vorinitialisiert werden.. */
                [li,lj] = meshgrid(index_glob);
                pos = Koff + (1:nR(k)^2);
                i(pos) = li(:);
                j(pos) = lj(:);
                K(pos) = K_lok(:);

                pos = Roff + (1:nR(k));
                ir(pos) = index_glob;
                R(pos) = R_lok;

                Roff = Roff + nR(k);
                Koff = Koff + nR(k)^2;
            end
            /*  Steifigkeitsmatrix für u und T */
            K = sparse(i,j,K,fdim,fdim);
            /*  Residuumsvektor */
            R = sparse(ir,ones(size(ir)),R,fdim,1);
        }
};
}
}