DefaultInfo.m

classdef DefaultInfo < Fem.IFemInfo & DataTree.ILeafNode
    %classdef LagrangeInfo < Fem.IFemInfo & DataTree.ILeafNode
    % scalar Lagrange FE functions on triangular grid
    
    properties (SetAccess = private)
        
        pdeg;
        dimrange = 1;
        ndofs;
        ndofs_per_element;
        grid;
        
        basis_weight_matrix;
        
        % index vector of global DOFs
        global_dof_index;
        
        dofs_lcoord;
        dofs_edges_llcoord;
        
        % weight matrix for `L^2` inner product
        l2_inner_product_matrix;
        
        % weight matrix for `H^1_0` inner product
        h10_inner_product_matrix;
        
        dof_ids_local;
    end
    
    properties (SetAccess = private, Dependent)
        
        detDF;
    end
    
    methods
        
        function df_info = DefaultInfo(pdeg, grid)
            %function df_info = DefaultInfo(pdeg, grid)
            %
            df_info.pdeg = pdeg;
            df_info.dimrange = 1;
            df_info.grid = grid;
            df_info.ndofs = lagrange_ndofs(df_info);
            df_info.ndofs_per_element = (pdeg+1) * (pdeg+2) / 2;
            
            df_info.dof_ids_local = {1:df_info.ndofs_per_element};
            
            df_info.global_dof_index = Fem.Lagrange.global_dof_index(df_info, grid);
            
            df_info.basis_weight_matrix = lagrange_basis_weight_matrix(df_info);
            
            df_info.dofs_lcoord = Fem.Lagrange.nodes_lcoord(pdeg);
            df_info.dofs_edges_llcoord = Fem.Lagrange.nodes_edges_llcoord(df_info.dofs_lcoord);
            
            % reaction term == 1 then reaction matrix is l2matrix
            model = [];
            model.disable_caching = 1;
            model.reaction = ...
                @(grid, eindices, loc, params) ones(length(eindices), 1);
            model.qdeg = 2*df_info.pdeg;
            df_info.l2_inner_product_matrix = Fem.Assembly.matrix_part(...
                @Fem.IntegralKernels.matrix_reaction, model, df_info);
            
            % H10_inner_product_matrix: identity diffusion
            model.diffusivity_tensor = @(grid, eindices, loc, params) ...
                [ones(length(eindices), 1), zeros(length(eindices), 1), ...
                zeros(length(eindices), 1), ones(length(eindices), 1)];
            df_info.h10_inner_product_matrix = Fem.Assembly.matrix_part(...
                @Fem.IntegralKernels.matrix_diffusion, model, df_info);        
        end
        
        % copy constructor
        function obj2 = copy(obj1)
            obj2 = Fem.Lagrange.DefaultInfo(obj1.pdeg, obj1.grid);
        end
        
        function detDF = get.detDF(this)
            detDF = this.grid.A(:)*2;
        end
        
        %function dof_ids_local = get.dof_ids_local(this)
        %    dof_ids_local = {1:this.ndofs_per_element};
        %end
        
        function res = evaluate_basis(this, lcoord)
            res = this.basis_weight_matrix * power_vector2(lcoord, this.pdeg);
        end
        
        function res = evaluate_scalar_basis_derivative(this, lcoord, ~)
            res = this.basis_weight_matrix * power_vector2_derivative(lcoord, this.pdeg);
        end
        
        function res = evaluate_basis_function(this, lcoord, i)
            res = this.basis_weight_matrix(i, :) * power_vector2(lcoord, this.pdeg);
        end
        
        function res = evaluate_basis_function_derivative(this, lcoord, i)
            res = this.basis_weight_matrix(i,:) * power_vector2_derivative(lcoord, this.pdeg);
        end
        
        function res = get_global_dof_index(this, varargin)
            if isempty(varargin)
                res = this.global_dof_index;
            else
                res = this.global_dof_index(varargin{:});
            end
        end
        
        function res = get_dof_ids_global(this, ~)            
            res = 1:this.ndofs;
        end
        
        function df = interpol_local(this, func, params)
            
            if nargin<3
                params = [];
            end
            
            df = Fem.DiscFunc([], this);            
            for i = 1:this.ndofs_per_element               
                df.dofs(this.global_dof_index(:, i)) = func(this.grid, 1:this.grid.nelements, ...
                    this.dofs_lcoord(i, :), params);
            end
        end
        
        p = plot_dofmap(this,params);
        
    end
    
    methods (Access = private, Hidden)
        
        function res = lagrange_ndofs(this)
            ndof_nodes = this.grid.nvertices;
            ndof_edges = (this.grid.nedges_interior + this.grid.nedges_boundary) * (this.pdeg-1);
            ndof_elements = this.grid.nelements * (this.pdeg-1) * (this.pdeg-2) / 2;
            res = ndof_nodes + ndof_edges + ndof_elements;
        end
        
        function V = lagrange_basis_weight_matrix(this)
            
            switch this.pdeg
                case 1
                    V = [1,-1,-1 ; ...
                        0, 1, 0; ...
                        0, 0, 1];
                case 2
                    V = [1    -3    -3     2     4     2
                        0     4     0    -4    -4     0
                        0    -1     0     2     0     0
                        0     0     4     0    -4    -4
                        0     0     0     0     4     0
                        0     0    -1     0     0     2
                        ];
                case 3
                    V = [
                        2   -11   -11    18    36    18    -9   -27   -27    -9
                        0    18     0   -45   -45     0    27    54    27     0
                        0    -9     0    36     9     0   -27   -27     0     0
                        0     2     0    -9     0     0     9     0     0     0
                        0     0    18     0   -45   -45     0    27    54    27
                        0     0     0     0    54     0     0   -54   -54     0
                        0     0     0     0    -9     0     0    27     0     0
                        0     0    -9     0     9    36     0     0   -27   -27
                        0     0     0     0    -9     0     0     0    27     0
                        0     0     2     0     0    -9     0     0     0     9
                        ] / 2;
                case 4
                    V = [
                        3   -25   -25    70   140    70   -80  -240  -240   -80    32   128   192   128    32
                        0    48     0  -208  -208     0   288   576   288     0  -128  -384  -384  -128     0
                        0   -36     0   228    84     0  -384  -432   -48     0   192   384   192     0     0
                        0    16     0  -112   -16     0   224    96     0     0  -128  -128     0     0     0
                        0    -3     0    22     0     0   -48     0     0     0    32     0     0     0     0
                        0     0    48     0  -208  -208     0   288   576   288     0  -128  -384  -384  -128
                        0     0     0     0   288     0     0  -672  -672     0     0   384   768   384     0
                        0     0     0     0   -96     0     0   480    96     0     0  -384  -384     0     0
                        0     0     0     0    16     0     0   -96     0     0     0   128     0     0     0
                        0     0   -36     0    84   228     0   -48  -432  -384     0     0   192   384   192
                        0     0     0     0   -96     0     0    96   480     0     0     0  -384  -384     0
                        0     0     0     0    12     0     0   -48   -48     0     0     0   192     0     0
                        0     0    16     0   -16  -112     0     0    96   224     0     0     0  -128  -128
                        0     0     0     0    16     0     0     0   -96     0     0     0     0   128     0
                        0     0    -3     0     0    22     0     0     0   -48     0     0     0     0    32
                        ]/3;
                otherwise
                    % these computations should not be done everytime!
                    % for higher pdeg, perhaps also use higher accuracy!!!
                    disp(['please hardcode the following matrix into' ...
                        ' Fem.Lagrange.DefaultInfo.lagrange_basis_weight_matrix:']);
                    % V = inv(P) where P = (p(l1),...p(lm)) powervector of lagrange-nodes
                    lagrange_nodes = this.dofs_lcoord;
                    P = zeros(size(lagrange_nodes,1));
                    for i = 1: size(lagrange_nodes,1)
                        P(:,i) = power_vector2(lagrange_nodes(i,:), this.pdeg);
                    end
                    V = inv(P);
            end
        end
        
    end
    
end