BaseGeometry.m

Go to the documentation of this file.

namespace fem{
namespace geometry{


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

        Nodes;
        Elements;
        Edges;
        DofsPerElement;


        MasterFaces;


        PatchFacesIdx;
        PatchesPerFace;


        Faces;
        OrientationCheckIndices;

        ReverseAxesIndices;

        
    public: /* ( Dependent ) */

        NumElements;

        NumNodes;

        NumFaces;

        NodesPerFace;


        Width;
        Depth;
        Height;
    public:



        FaceNormals = "[-1 1  0 0  0 0\
                        0 0 -1 1  0 0\
                        0 0  0 0 -1 1]";
        FaceDims = logical("[0 0 1 1 1 1\
                            1 1 0 0 1 1\
                            1 1 1 1 0 0]");


        PatchFaces;

        
    public:


        BaseGeometry() {
        }


        function  plot(withdofnr,elem,pm) {
            
            if nargin < 4
                if nargin < 3
                    elem = 1:this.NumElements;
                    if nargin < 2
                        withdofnr = false;
                    end
                end
                pm = PlotManager;/* (false,1,2); */

                pm.LeaveOpen= true;
            end

            p = this.Nodes;
            h = pm.nextPlot(" geometry ",sprintf(" Geometry (%s) ",class(this)), " x [mm] ", " y [mm] ");
            plot3(h,p(1,:),p(2,:),p(3,:)," k. "," MarkerSize ",14);
            hold(h," on ");
            for k = 1:length(elem)
                el = this.Elements(elem(k),:);
                center = sum(p(:,el),2)/this.DofsPerElement;
                text(center(1),center(2),center(3),sprintf(" E_{%d} ",elem(k))," Parent ",h," Color "," m ");
                /*  Plot local numbering for first element */
                if k == 1
                    off = .04;
                    for i = 1:this.DofsPerElement
                        text(off+p(1,el(i)),off+p(2,el(i)),off+p(3,el(i)),sprintf(" %d ",i)," Parent ",h," Color "," r ");
                    end
                    for i = 1:6
                        facenodes = p(:,el(this.MasterFaces(i,:)));
                        mid = mean(facenodes,2);
                        text(mid(1),mid(2),mid(3),sprintf(" F_%d ",i)," Parent ",h," Color "," k ");
                    end
                end
            end
            eg = this.Edges;
            for i = 1:size(eg,1)
                plot3(h,[p(1,eg(i,1)) p(1,eg(i,2))],[p(2,eg(i,1)) p(2,eg(i,2))],[p(3,eg(i,1)) p(3,eg(i,2))]," r ");
            end
            if withdofnr
                for k = 1:this.NumNodes
                    text(p(1,k),p(2,k),p(3,k),sprintf(" %d ",k)," Parent ",h," Color "," k ");
                end
            end
            zlabel(" z [mm] ");
            view(h,[52 30]);
            daspect([1 1 1]);

            if nargin < 3
                pm.done;
            end
        }



        function bounds = getBoundingBox(marginfac) {
            if nargin < 2
                marginfac = 1;
            end
            m = min(this.Nodes,[],2)*marginfac;
            M = max(this.Nodes,[],2)*marginfac;
            bounds = [m(1) M(1) m(2) M(2) m(3) M(3)];
        }


        function commonidx = getCommonNodesWith(other) {
            commonidx = Utils.findVecInMatrix(this.Nodes,other.Nodes);
        }


        function  swapYZ() {
            error(" not implemented ");
        }


        function copy = scale(factor) {
            mc = metaclass(this);
            pts = this.Nodes * factor;
            cubes = this.Elements;
            eval(sprintf(" copy = %s(pts,cubes); ",mc.Name));
        }


        function  reverseAxis(dim) {
            this.reverseElementAxis(dim, 1:this.NumElements);
        }


        function [sub , nodeidx , faces ] = getSubMesh(elems,faces) {
            newelems = this.Elements(elems,:)^t;
            nodeidx = unique(newelems(:)," stable ");
            newnodes = this.Nodes(:,nodeidx);/* #ok
             * Fit element node reference indices */

            invidx(nodeidx) = 1:length(nodeidx);
            newelems = invidx(newelems)^t;/* #ok */

            sub = eval([class(this) " (newnodes, newelems) "]);

            /*  Auto-detect which faces are included for the submesh's
             * elements if none are specified */
            if nargin < 3
                isface = false(1,size(this.NumFaces,2));
                for k = elems
                    isface = isface | this.Faces(1,:) == k;
                end
                faces = find(isface);
            end
            sub.Faces= this.Faces(:,faces);

            /*  Fit face reference indices */
            invidx = [];
            invidx(elems) = 1:length(elems);
            sub.Faces(1,:) = invidx(sub.Faces(1,:));
            sub.faceComputations;
        }


        function  toCMGUI(folder,name) {
            if nargin < 3
                name = " musclegeom ";
                if nargin < 2
                    folder = pwd;
                end
            end
            /* % exnode file */
            exnode = fullfile(folder,[name " .exnode "]);
            f = fopen(exnode," w+ ");
            fprintf(f," Group name : %s\n ",name);
            fprintf(f," #Fields=1\n ");
            fprintf(f," 1) coordinates, coordinate, rectangular cartesian, #Components=3\n ");
            fprintf(f," x.  Value index= 1, #Derivatives= 0\n ");
            fprintf(f," y.  Value index= 2, #Derivatives= 0\n ");
            fprintf(f," z.  Value index= 3, #Derivatives= 0\n ");
            for n = 1:this.NumNodes
                fprintf(f," Node:\t%d\n ",n);
                fprintf(f," %E\n ",this.Nodes(:,n));
            end
            fclose(f);

            /* % exelem file */
            exelem = fullfile(folder,[name " .exelem "]);
            f = fopen(exelem," w+ ");
            fprintf(f," Group name : %s\n ",name);
            /* fprintf(f,'#Fields=0\n'); */

            /*  1D elements - lines */
            fprintf(f," Shape.  Dimension=1\n ");
            for e = 1:size(this.Edges,1)
                fprintf(f," Element: 0 0 %d\n ",e);
            end

            /*  2D elements - faces */
            fprintf(f," Shape.  Dimension=2\n ");
            for fa = 1:size(this.Faces,2)
                fprintf(f," Element: 0 %d 0\n ",fa);
                fprintf(f,"   Faces:\n ");
                elem = this.Faces(1,fa);
                /*  PatchFacesIdx returns the node positions in an order you can
                 * derive lines from it, MasterFaces doesn't */
                facenodeidx = this.Elements(elem,this.PatchFacesIdx(this.Faces(2,fa),:));
                /*  Insert last edge closing the face */
                facenodeidx = sort([facenodeidx; facenodeidx(2:end) facenodeidx(1)]);
                edgenumbers = Utils.findVecInMatrix(this.Edges^t,facenodeidx);
                fprintf(f,"     0 0 %d\n ",edgenumbers);
            end

            /*  3D elements */
            fprintf(f," Shape.  Dimension=3\n ");
            fprintf(f," #Scale factor sets=0\n ");
            fprintf(f," #Nodes=0\n ");
/*              fprintf(f,'#Nodes=%d\n',size(this.Elements,2)); */
            fprintf(f," #Fields=0\n ");
/*              fprintf(f,'1) coordinates, coordinate, rectangular cartesian, #Components=3\n');
 *             fprintf(f,'x.  l.Lagrange*l.Lagrange*l.Lagrange, no modify, standard node based.\n');
 *             fprintf(f,'y.  l.Lagrange*l.Lagrange*l.Lagrange, no modify, standard node based.\n');
 *             fprintf(f,'z.  l.Lagrange*l.Lagrange*l.Lagrange, no modify, standard node based.\n'); */
            for m = 1:this.NumElements
                fprintf(f," Element: %d 0 0\n ",m);
                /*  See if this element has any faces */
                facesofelem = find(this.Faces(1,:) == m);
                if ~isempty(facesofelem)
                    fprintf(f,"   Faces:\n ");
                    fprintf(f,"     0 %d 0\n ",facesofelem);
                end
                fprintf(f,"   Nodes: %s\n ",sprintf(" %d\t ",this.Elements(m,:)));
            end

            fclose(f);
        }


    
    protected:

        function faces = computeFaces() {
            if this.NumElements > 1
                inner = zeros(this.NumElements,6);
                for Afaceidx = 1:6
                    for Bfaceidx = 1:6
                        if Afaceidx ~= Bfaceidx
                            /*  Also sort the element node indices in case of
                             * nonidentical numbering within each element */
                            A = sort(this.Elements(:,this.MasterFaces(Afaceidx,:)),2);
                            B = sort(this.Elements(:,this.MasterFaces(Bfaceidx,:)),2);
                            for l = 1:this.NumElements
                                opposing = sum(abs(A - circshift(B,l-1)),2) == 0;
                                if any(opposing)
                                    opposedto = circshift(opposing,-(l-1));
                                    inner(opposing, Afaceidx) = find(opposedto);
                                    inner(opposedto, Bfaceidx) = find(opposing);
                                end
                            end
                        end
                    end
                end
                [elemnr, facenr] = find(inner == 0);
                faces = [elemnr facenr]^t;
            else
                mf = size(this.MasterFaces,1);
                faces = [ones(1,mf); 1:mf];
            end
            this.faceComputations(faces);
        }
        function  faceComputations(faces) {
            if ~isempty(this.PatchFacesIdx)
                if nargin < 2
                    faces = this.Faces;
                end
                nf = size(faces,2);
                ppf = this.PatchesPerFace;
                /* np = size(this.PatchFacesIdx,1); */
                pf = zeros(nf*ppf,size(this.PatchFacesIdx,2));
                for k=1:nf
                    elem = faces(1,k);
                    face = faces(2,k);
                    pos = (face-1)*ppf + (1:ppf);
                    patchidx = this.PatchFacesIdx(pos,:);
                    pf((k-1)*ppf+(1:ppf),:) = reshape(this.Elements(elem,patchidx),ppf,[]);
                end
                this.PatchFaces= pf;
            end
        }
        function  checkOrientation() {
            return;
            oi = this.OrientationCheckIndices;
            pi = ProcessIndicator(" %s: Checking orientation of %d elements in x,y,z directions ",3*this.NumElements,false,class(this),this.NumElements);
            for dim = 1:3
                for e = 1:this.NumElements
                    chk = diff(reshape(this.Nodes(dim,this.Elements(e,oi(:,:,dim)" ) "),size(oi,2),[]),[],1);
                    if any(chk <= 0)
                        if all(chk < 0)
                            this.reverseElementAxis(dim, e);
                            /* fprintf('Axis %d has negative orientation. Reversing element %d\n',dim,e); */
                        else
                            [ip,jp] = find(chk > 0);
                            [in,jn] = find(chk < 0);
                            [is,js] = find(chk == 0);
                            if ~isempty(is)
                                errel = this.Elements(e,oi(js,[is is+1]^t,dim));
                                fprintf(" \n%s: Degenerate element nr %3d in dimension %d: %3d equal! Nodes %s (%s)\n ",...
                                    class(this),e,dim,length(is),int2str(errel),num2str(this.Nodes(dim,errel)));
                            else
                                if length(ip) > length(in)
                                    i = in; j = jn;
                                    str = " positive ";
                                    reverse = false;
                                else
                                    i = ip; j = jp;
                                    str = " negative ";
                                    reverse = true;
                                end
                                errel = this.Elements(e,oi(j,[i i+1],dim));
                                fprintf(" \n%s: Opposite orientation in dimension %d in mainly %s oriented element nr %3d (%3d pos/%3d neg): Nodes %s (%s)\n ",...
                                    class(this),dim,str,e,length(ip),length(in),int2str(errel),num2str(this.Nodes(dim,errel)));
                                if reverse
                                    this.reverseElementAxis(dim, e);
                                end
                            end
                        end
                    end
                    pi.step;
                end
            end
            pi.stop;
        }

    
    private:

        function  reverseElementAxis(dim,e) {
            this.Elements(e,:) = this.Elements(e,this.ReverseAxesIndices(dim,:));
        }

    
    public:

        
#if 0 //mtoc++: 'get.NumElements'
function nc = NumElements() {
            nc = size(this.Elements,1);
        }

#endif


        
#if 0 //mtoc++: 'get.NumNodes'
function nc = NumNodes() {
            nc = size(this.Nodes,2);
        }

#endif


        
#if 0 //mtoc++: 'get.NumFaces'
function nf = NumFaces() {
            nf = size(this.Faces,2);
        }

#endif



        
#if 0 //mtoc++: 'get.NodesPerFace'
function npf = NodesPerFace() {
            npf = size(this.MasterFaces,2);
        }

#endif


        
#if 0 //mtoc++: 'get.Width'
function w = Width() {
            w = max(this.Nodes(1,:))-min(this.Nodes(1,:));
        }

#endif


        
#if 0 //mtoc++: 'get.Depth'
function w = Depth() {
            w = max(this.Nodes(2,:))-min(this.Nodes(2,:));
        }

#endif


        
#if 0 //mtoc++: 'get.Height'
function w = Height() {
            w = max(this.Nodes(3,:))-min(this.Nodes(3,:));
        }

#endif

    
    public: /* ( Static ) */

        static function res = test_DemoGrids() {
            fem.geometry.RegularHex8Grid(1:3);
            fem.geometry.RegularHex8Grid(1:3,1:4);
            fem.geometry.RegularHex8Grid(1:3,1:4,-1:3);
            fem.geometry.RegularHex8Grid(1:3,1:4,-1:3,.2);
            fem.geometry.RegularHex20Grid(1:2);
            fem.geometry.RegularHex20Grid(1:2,0:2);
            fem.geometry.RegularHex20Grid(-1:1,1:3,-1:1);
            fem.geometry.RegularHex20Grid(-1:1,1:2,-1:1,.2);
            res = true;
        }


        static function  test_subMesh() {
            g = fem.geometry.RegularHex27Grid(-1:1,1:2,-1:1,.2);
            [sg, node] = g.getSubMesh(1:3:g.NumElements);
            sg.plot;
            g = fem.geometry.RegularHex8Grid(-20:1,1:4,-1:1,.2);
            [sg, node] = g.getSubMesh([1 5 9 10 15 56 79 99]);
            sg.plot;
        }

    
    
};
}
}