CoreFun3D.m

Go to the documentation of this file.

namespace models{
namespace pcd{


/* (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 CoreFun3D
  :public models.pcd.BaseCoreFun {
    public: /* ( setObservable ) */

        tDecaySecs = 2;
    private:

        A;


        nodes;


        g;

        
    public:


        CoreFun3D(dynsys) {
            this = this@models.pcd.BaseCoreFun(dynsys);
        }


        function copy = clone() {
            copy = models.pcd.CoreFun3D(this.sys);

            /*  Call superclass method */
            copy = clone@models.pcd.BaseCoreFun(this, copy);

            /*  copy reference!
             *copy.sys = this.sys; % already done in constructor */
            copy.A= this.A;
            copy.dim= this.dim;
            copy.g= this.g;
            copy.nodes= this.nodes;
        }


        function  newSysDimension() {

            /*  Create diffusion matrix */
            n = prod(this.System.Dims);
            this.nodes= this.System.geo.Points;
            this.fDim= 4*n;
            this.xDim= 4*n;
        }
        function fx = evaluateCoreFun(colvec<double> x,double t) {
            fx = zeros(size(x));

            m = this.nodes;
            s = this.System;
            n = this.xDim;
            g = s.geo;

            /*  Uncomment if reaction coeffs become real params again */
            mu = [s.ReacCoeff; this.mu];

            /*  Extract single functions */
            xa = x(1:m);
            xan = xa.^n;
            ya = x(m+1:2*m);
            xi = x(2*m+1:3*m);
            yi = x(3*m+1:end);

            /* % Compile boundary conditions
             * Determine the points on each side that are affected of
             * boundary conditions (constant over each simulation as mu is
             * fixed then) */
            h = s.h;
            a = s.Omega;
            xr = a(1,2)-a(1,1);
            yr = a(2,2)-a(2,1);
            zr = a(3,2)-a(3,1);
            rb = zeros(m,1);
            ud = (t < this.tDecaySecs) + (t >= this.tDecaySecs)*max(0,2-.5*t);

            /* % Front & back (only x/y coords relevant) */
            [i,j] = ind2sub(g.Dims,g.F); /* #ok<*PROP> */

            xd = abs((i-1)*h-.5*xr);
            yd = abs((j-1)*h-.5*yr);
            rbp = xd < xr*mu(9)/2 & yd < yr*mu(9)/2;
            rbF = sub2ind(g.Dims, i(rbp), j(rbp), ones(size(i(rbp))));
            rbp = xd < xr*mu(10)/2 & yd < yr*mu(10)/2;
            rbBa = sub2ind(g.Dims, i(rbp), j(rbp), ones(size(i(rbp)))*g.Dims(3));
            rb(rbF) = (xi(rbF)*mu(15)*ud)/h;
            rb(rbBa) = rb(rbBa) + (xi(rbBa)*mu(16)*ud)/h;

            /* % Left & Right (only x/z coords relevant) */
            [i,~,k] = ind2sub(g.Dims,g.L);
            xd = abs((i-1)*h-.5*xr);
            zd = abs((k-1)*h-.5*zr);
            rbp = xd < xr*mu(11)/2 & zd < zr*mu(11)/2;
            rbL = sub2ind(g.Dims, i(rbp), ones(size(i(rbp))), k(rbp));
            rbp = xd < xr*mu(12)/2 & zd < zr*mu(12)/2;
            rbR = sub2ind(g.Dims, i(rbp), ones(size(i(rbp)))*g.Dims(2), k(rbp));
            rb(rbL) = rb(rbL) + (xi(rbL)*mu(17)*ud)/h;
            rb(rbR) = rb(rbR) + (xi(rbR)*mu(18)*ud)/h;

            /* % Top & Bottom (only y/z coords relevant) */
            [i,j,k] = ind2sub(g.Dims,g.T);
            yd = abs((j-1)*h-.5*yr);
            zd = abs((k-1)*h-.5*zr);
            rbp = yd < yr*mu(13)/2 & zd < zr*mu(13)/2;
            rbT = sub2ind(g.Dims, ones(size(i(rbp))), j(rbp), k(rbp));
            rbp = yd < yr*mu(14)/2 & zd < zr*mu(14)/2;
            rbBo = sub2ind(g.Dims, ones(size(i(rbp)))*g.Dims(1), j(rbp), k(rbp));
            rb(rbT) = rb(rbT) + (xi(rbT)*mu(19)*ud)/h;
            rb(rbBo) = rb(rbBo) + (xi(rbBo)*mu(20)*ud)/h;

            /* % Actual evaluation
             * Xa */
            fx(1:m) = mu(1)*xi.*ya - mu(3)*xa + rb;
            /*  Ya */
            fx(m+1:2*m) = mu(2)*yi.*xan - mu(4)*ya;
            /*  Xi */
            fx(2*m+1:3*m) = -mu(1)*xi.*ya - mu(5)*xi + mu(7) - rb;
            /*  Yi */
            fx(3*m+1:end) = -mu(2)*yi.*xan - mu(6)*yi + mu(8);
        }
    protected:

        function fxj = evaluateComponents(pts,ends,unused1,unused2,X,double t) {
            error(" Not yet implemented ");
        }


};
}
}