BubbleInfo.m

classdef BubbleInfo < Fem.IFemInfo & DataTree.ILeafNode
    %classdef LagrangeBubbleInfo < Fem.IFemInfo & DataTree.ILeafNode
    % scalar Lagrange FE functions (degree 1) with additional bubble function
    % on triangular grid
    
    % could possibly be generalized to arbitrary pdeg of Lagrange basis
    
    % I. Maier 14.08.2012
    
    properties (SetAccess = private)
        
        pdeg = 3;
        
        dimrange = 1;
        
        ndofs;
        
        ndofs_per_element = 4;
        
        grid;
        
        basis_weight_matrix;
        
        dofs_lcoord;
        
        dofs_edges_llcoord;
        
        l2_inner_product_matrix;
        
        h10_inner_product_matrix;
        
        dof_ids_local;
    end
    
    properties (SetAccess = private, Dependent)
        
        detDF;
    end
    
    methods
        
        function lbi = BubbleInfo(~, grid)
            %function velocity_df_info = BubbleInfo(~, grid)
            % constructor
            
            lbi.dof_ids_local = {1:lbi.ndofs_per_element};
            
            lbi.grid = grid;
            lbi.ndofs = grid.nvertices + grid.nelements;
            
            % basis weight matrix
            lbi.basis_weight_matrix = ...
                [1,-1,-1, 0, 0, 0, 0, 0, 0, 0; ...
                0, 1, 0, 0, 0, 0, 0, 0, 0, 0; ...
                0, 0, 1, 0, 0, 0, 0, 0, 0, 0; ...
                0, 0, 0, 0,27, 0, 0,-27,-27, 0];
            
            % local coordinates of dofs
            lbi.dofs_lcoord = [0, 0; 1, 0; 0, 1; 1/3 1/3];
            
            % edge-local coordinates of dofs
            lbi.dofs_edges_llcoord = Fem.Lagrange.nodes_edges_llcoord(lbi.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*lbi.pdeg;
            lbi.l2_inner_product_matrix = Fem.Assembly.matrix_part(...
                @Fem.IntegralKernels.matrix_reaction, model, lbi);
            
            % 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)];
            lbi.h10_inner_product_matrix = Fem.Assembly.matrix_part(...
                @Fem.IntegralKernels.matrix_diffusion, model, lbi);
        end
        
        % copy constructor
        function obj2 = copy(obj1)
            obj2 = Fem.Lagrange.BubbleInfo([], obj1.grid);
        end
        
        %function dof_ids_local = get.dof_ids_local(this)
        %    dof_ids_local = {1:this.ndofs_per_element};
        %end
        
        function detDF = get.detDF(this)
            detDF = 2 * this.grid.A;
        end
        
        function res = get_global_dof_index(this, einds, gids)
            % returns entries of global-dof-map for several elements and several
            % local dof ids
            global_dof_index = ...
                [this.grid.VI, this.grid.nvertices + (1:this.grid.nelements)'];
            res = global_dof_index(einds, gids);
        end
        
        function res = get_dof_ids_global(this, ~)
            % returns global dof ids belonging to one component (indices into
            % DOF-vector)
            res = 1:this.ndofs;
        end
        
        function res = evaluate_basis(this, lcoord)
            % evaluation of all basis functions in local coordinate
            res = this.basis_weight_matrix * power_vector2(lcoord, this.pdeg);
        end
        
        function res = evaluate_scalar_basis_derivative(this, lcoord, ~)
            % evaluation of the gradient of a scalar basis in local coordinate
            res = this.basis_weight_matrix * power_vector2_derivative(lcoord, this.pdeg);
        end
        
        function res = evaluate_basis_function(this, lcoord, i)
            % evaluation of i-th basis function in local coordinate
            res = this.basis_weight_matrix(i, :) * power_vector2(lcoord, this.pdeg);
        end
        
        function res = evaluate_basis_function_derivative(this, lcoord, i)
            % evaluation of the gradient of i-th basis function in local
            % coordinate
            res = this.basis_weight_matrix(i,:) * power_vector2_derivative(lcoord, this.pdeg);
        end
        
        function df = interpol_local(this, func, params)
            % interpolating a given analytical function
            % input:
            % - func: function which can be evaluated at
            % @(grid, elinds, lcoord, params)
            
            if nargin<3
                params = [];
            end
            
            df = Fem.DiscFunc([], this);
            
            for i = 1:(this.ndofs_per_element - 1)
                
                gids = this.get_global_dof_index(1:this.grid.nelements, i);
                df.dofs(gids) = func(this.grid, 1:this.grid.nelements, ...
                    this.dofs_lcoord(i, :), params);
            end
            
            interp_error = triaquadrature(3, @(lcoord) ...
                func(this.grid, 1:this.grid.nelements, lcoord, params) - ...
                df(1:this.grid.nelements, lcoord));
            
            bubble_integr = triaquadrature(3, ...
                @(lcoord) this.evaluate_basis_function(lcoord, this.ndofs_per_element));
            
            gids = this.get_global_dof_index(1:this.grid.nelements, this.ndofs_per_element);
            df.dofs(gids) = interp_error / bubble_integr(1); % detDF cancels out
        end
        
    end
    
end