BcInfo.m

classdef BcInfo < handle
    %classdef BcInfo < handle
    % Class for boundary condition information. Allows more flexible usage.
    %
    % One can specify boundary parts via model.bnd_rect_corner1 and
    % model.bnd_rect_corner2. For boundary integrals one has to group the
    % rectangles with numbers < -1 in model.bnd_rect_index. For each of
    % these rectangles an integral kernel and coefficient functions can be
    % given in model. Then, model.matrix_boundary_int_kernel{i} or
    % model.rhs_boundary_int_kernel{i} is used on the boundary with index
    % -(i+1). If a term is zero, there should be an empty cell entry.
    %
    % Boundaries not assigned in this way are used for Dirichlet conditions
    % by default. It is also possible to define Dirichlet constraints for
    % chosen dimensions of the field variable in model.bnd_dimrange_index.
    % Then the rectangles have to be given explicitly (after the rectangles
    % for the integral terms).
    %
    % The Fem assembly requires a field named bc_info in model_data, which
    % is an object of this class.
    
    properties
        
        grid;
        df_info;
        nbint_parts;
        elinds;
        edgeinds;
        
        has_dirichlet_values;
        dirichlet_gids;
        dirichlet_dof_vector_components;
    end
    
    methods
        
        function bi = BcInfo(model, grid, df_info)
            %bi = BcInfo(model, grid, df_info)
            % Constructor preparing all necessary boundary condition
            % information for the Fem assembly.
            
            bi.grid = grid;
            bi.df_info = df_info;
            
            bi.nbint_parts = -min(min(grid.NBI)) - 1;
            bi.elinds = cell(bi.nbint_parts, 1);
            bi.edgeinds = cell(bi.nbint_parts, 1);
            
            for i = 1:bi.nbint_parts
                
                ind = find(grid.NBI == -(i+1));
                [bi.elinds{i}, bi.edgeinds{i}] = ind2sub(size(grid.NBI), ind);
            end
            
            bi.has_dirichlet_values = model.has_dirichlet_values;
            bi.dirichlet_gids = [];
            bi.dirichlet_dof_vector_components = {};
            
            if bi.has_dirichlet_values
                
                % find dirichlet indices               
                if (isfield(model, 'bnd_rect_index') && ...
                        size(model.bnd_rect_corner1, 2) > length(model.bnd_rect_index)) ...
                        || (~isfield(model, 'bnd_rect_index') && ...
                        isfield(model, 'bnd_dimrange_index'))
                    
                    bnd_ind = find(grid.NBI <= 0);
                    SX = grid.ECX(bnd_ind);
                    SY = grid.ECY(bnd_ind);
                    
                    % field bnd_dimrange_index of model is required
                    for i = 1:length(model.bnd_dimrange_index)
                        
                        dimrange_index = model.bnd_dimrange_index{i};
                        dir_corner_index = length(model.bnd_rect_index) + i;
                        dof_ids_local = cat(2, df_info.dof_ids_local{dimrange_index});
                        
                        inds = (SX > model.bnd_rect_corner1(1, dir_corner_index)) & ...
                            (SX < model.bnd_rect_corner2(1, dir_corner_index)) & ...
                            (SY > model.bnd_rect_corner1(2, dir_corner_index)) & ...
                            (SY < model.bnd_rect_corner2(2, dir_corner_index));
                        
                        bnd_inds_tmp = bnd_ind(inds);
                        [elids, edgeids] = ind2sub(size(grid.NBI), bnd_inds_tmp);
                        
                        for j = 1:grid.nneigh
                            
                            elids_tmp = elids(edgeids == j);
                            lids = find(df_info.dofs_edges_llcoord(:, j) > -1);
                            lids = intersect(lids, dof_ids_local);
                            
                            if ~isempty(lids) && ~isempty(elids_tmp)
                                
                                gids = get_global_dof_index(df_info, elids_tmp, lids);
                                bi.dirichlet_gids = union(bi.dirichlet_gids, gids(:));                             
                            end
                        end
                    end
                else
                    
                    if isfield(model, 'bnd_dimrange_index')
                        
                        dof_ids_local = cat(2, df_info.dof_ids_local{model.bnd_dimrange_index{1}});
                    end
                    
                    % default procedure
                    for j = 1:grid.nneigh
                        
                        elids = find(grid.NBI(:, j) == -1);
                        lids = find(df_info.dofs_edges_llcoord(:, j) > -1);
                        if isfield(model, 'bnd_dimrange_index')
                            lids = intersect(lids, dof_ids_local);
                        end
                        if ~isempty(lids) && ~isempty(elids)
                            
                            gids = get_global_dof_index(df_info, elids, lids);
                            bi.dirichlet_gids = union(bi.dirichlet_gids, gids(:));
                        end
                    end
                end
                
                old_decomp_mode = model.decomp_mode;
                model.decomp_mode = 1;
                bi.dirichlet_dof_vector_components = ...
                    Fem.Assembly.dirichlet_dof_vector(model, df_info, bi);
                model.decomp_mode = old_decomp_mode;
                
            end %has_dirichlet_values
            
        end %constructor
        
    end %methods
end %classdef