AFEMConfig.m

Go to the documentation of this file.

namespace fem{


/* (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 AFEMConfig
  :public handle {
    public:

        FEM;


        Geometry;


        Model;

        
    public:

        char NeumannCoordinateSystem = "local";
        Options;

        
    public:

        VelocityBCTimeFun;
    private:

        iP;

        optArgs = {""};

        
    public:

        AFEMConfig(varargin) {
            this.optArgs= varargin;

            
            i.KeepUnmatched= true;
            this.iP= i;

            this.addOption(" GeoNr ",1);
            mc = metaclass(this);
            this.addOption(" Tag ",mc.Name);
        }


        function  configureModel(models.BaseFullModel model) {
            this.Model= model;
        }
        function  configureModelFinal() {
        }


        function  prepareSimulation(colvec<double> mu,integer inputidx) {

            /*  do nothing by default */
        }
        function P = getBoundaryPressure(elemidx,faceidx) {
            P = [];
        }
        function str = getOptionStr(withtag) {
            if nargin < 2
                withtag = true;
            end
            o = this.Options;
            fieldnames = fields(o);
            strs = [];
            for k=1:length(fieldnames)
                if (withtag || ~strcmp(fieldnames[k]," Tag ")) && ~isempty(o.(fieldnames[k]))
                    fmt = " %s-%g ";
                    if isa(o.(fieldnames[k])," char ")
                        fmt = " %s-%s ";
                    end
                    strs[end+1] = sprintf(fmt,fieldnames[k],o.(fieldnames[k]));/* #ok */

                end
            end
            str = Utils.implode(strs," _ ");
        }


        function  plotGeometryInfo(allnode,elemnr) {
            if nargin < 3
                elemnr = 1;
                if nargin < 2
                    allnode = false;
                end
            end
            g = this.FEM.Geometry;
            g.plot(allnode,elemnr);
        }

    
    protected:


        function  init() {
            this.iP.parse(this.optArgs[:]);
            this.Options= this.iP.Results;

            /* % Get the geometry */
            geo = this.getGeometry;
            if isa(geo," fem.geometry.Cube27Node ")
                this.FEM= fem.HexahedronTriquadratic(geo);
            elseif isa(geo," fem.geometry.Cube20Node ")
                this.FEM= fem.HexahedronSerendipity(geo);
            elseif isa(geo," fem.geometry.Cube8Node ")
                this.FEM= fem.HexahedronTrilinear(geo);
            end
            this.Geometry= geo;
        }
        function [velo_dir , velo_dir_val ] = setVelocityDirichletBC(velo_dir,velo_dir_val) {
        }
    protected: /* ( Sealed ) */

        function  addOption(name,default,varargin) {
            this.iP.addParamValue(name, default, varargin[:]);
        }

    
    protected: /* ( Abstract ) */

        virtual function displ_dir = setPositionDirichletBC(displ_dir) = 0;


        virtual function geo = getGeometry() = 0;
    public: /* ( Sealed ) */

        function [displ_dir , velo_dir , velo_dir_val ] = getBC() {
            N = this.Geometry.NumNodes;
            displ_dir = false(3,N);
            displ_dir = this.setPositionDirichletBC(displ_dir);
            velo_dir = false(3,N);
            velo_dir_val = zeros(3,N);
            [velo_dir, velo_dir_val] = this.setVelocityDirichletBC(velo_dir, velo_dir_val);

            if any(any(displ_dir & velo_dir)) 
                error(" Cannot impose displacement and velocity dirichlet conditions on same DoF ");
            end
        }


        function [force , nodeidx , faceswithforce ] = getSpatialExternalForces() {
            fe = this.FEM;
            geo = fe.Geometry;
            ngp = fe.GaussPointsPerElemFace;
            force = zeros(geo.NumNodes * 3,1);
            faceswithforce = false(1,geo.NumFaces);

            globalcoord = strcmp(this.NeumannCoordinateSystem," global ");
            for fn = 1:geo.NumFaces
                elemidx = geo.Faces(1,fn);
                faceidx = geo.Faces(2,fn);
                masterfacenodeidx = geo.MasterFaces(faceidx,:);
                /*  So far: Constant pressure on all gauss points! */
                P = this.getBoundaryPressure(elemidx, faceidx);
                if ~isempty(P)
                    faceswithforce(fn) = true;
                    integrand = zeros(3,geo.NodesPerFace);
                    if globalcoord
                        N = geo.FaceNormals(:,faceidx);
                    end
                    for gi = 1:ngp
                        if ~globalcoord
                            N = fe.NormalsOnFaceGP(:,gi,fn);
                        end
                        PN = (P * N) * fe.Ngpface(:,gi,fn)^t;
                        integrand = integrand + fe.FaceGaussWeights(gi)*PN*fe.face_detjac(fn,gi);
                    end
                    facenodeidx = geo.Elements(elemidx,masterfacenodeidx);
                    facenodeidx = (facenodeidx-1)*3+1;
                    facenodeidx = [facenodeidx; facenodeidx+1; facenodeidx+2];/* #ok */

                    force(facenodeidx(:)) = force(facenodeidx(:)) + integrand(:);
                end
            end
            /*  Augment to u,v,w vector */
            nodeidx = find(force);
        }

    

};
}