Mesh.m

classdef Mesh
    % Mesh Class
    % Works as a wrapper for the node coordinate matrix (coords) and the
    % element connectivity matrix (connect). Provides a series of helper
    % methods to assemble functions along the mesh.
    % Its main methods are assemble_matrix and assemble_vector.
    % Currently contains specific shell methods  and properties such as 
    % nodal_system and normals.
    properties
        coords
        connect
        normals
        thickness
        laminate_ids
        laminates
        ele_type
        nodes_with_mass
        mass_values
    end
    properties (Dependent)
        n_nodes
        n_ele
        nodes_per_ele
    end
    methods
        function obj = Mesh(ele_type,laminates,laminate_ids,coords,connect,thickness_in)
            % obj = Mesh(coords,connect,node_normals)
            % Sanity Checks
            require(isnumeric(coords) && isnumeric(connect), ...
                'ArgumentError: all arguments should be numeric')
            require(all(all(~mod(connect,1))), ...
                'ArgumentError: connect should be all integers')
            require(max(max(connect))==size(coords,1), ...
                'ArgumentError: biggest node in connect should be size coords')
            % Setting variables
            obj.ele_type = ele_type;
            obj.coords = coords;
            obj.connect = connect;
            obj.laminates = laminates;
            if isempty(laminate_ids)
                obj.laminate_ids = ones(size(connect,1),1);
            else
                require(length(laminate_ids)==size(connect,1), ...
                    'ArgumentError: length(laminate_ids) should = size(connect,2)');
                obj.laminate_ids = laminate_ids;
            end
            if EleType.is_shell(ele_type)
                obj.normals = obj.nodal_systems();
                obj.thickness = thickness_in;
            end
        end
        function nodal_systems = nodal_systems(mesh)
            % nodal_systems = nodal_systems(mesh)
            % nodal_systems [3x3xn_nodes]: One coordinate system per node
            % Shell specific
            switch (mesh.ele_type)
                case {EleType.AHMAD4}
                    type = EleType.Q4;
                case {EleType.AHMAD8}
                    type = EleType.Q8;
                case {EleType.AHMAD9}
                    type = EleType.Q9;
            end
            ksi = [ -1  1  1 -1  0  1  0 -1  0 ];
            eta = [ -1 -1  1  1 -1  0  1  0  0 ];
            tol = 1e-9;     % Used for length comparisson
            nodal_systems = zeros(3,3,mesh.n_nodes);
            for n = 1:mesh.n_nodes
                % eles:           Find the elements where the node is used
                % local_index:    Find which node it represents to the element
                [eles, local_index] = find(mesh.connect == n);
                n_eles = length(eles);     % Number of connect the node belongs to.
                % Loop through those connect to find v
                v = zeros(n_eles,3);      % Vector directions (v3) of the node in each element
                for e = 1:n_eles
                    % Create the element
                    ele_nodes = mesh.connect(eles(e),:);
                    ele_coords = mesh.coords(ele_nodes,:);
                    element = Element(1,type,[],ele_coords,[],[]);
                    
                    % Get 2D Jacobian of the element
                    % dN = element.dN(ksi(localNode(e)),eta(localNode(e)));
                    jac = element.jacobian(ksi(local_index(e)),eta(local_index(e)),0);
                    
                    % Use the jacobian to get v3
                    v3 = cross(jac(1,:),jac(2,:));
                    v(e,:) = v3/norm(v3);
                end
                aux_v3 = mean(v,1)';
                % V1 and V2 are totally dependet on V3
                % It is a choice to calculate them here and store them.
                % They could easily be dependent properties form v3
                v3 = aux_v3/norm(aux_v3);
                v1_aux = cross([0 1 0]',v3);
                v1_norm = norm(v1_aux);
                if v1_norm > tol
                    v1 = v1_aux/v1_norm;
                else
                    v1 = [1 0 0]';
                end
                % v2 is orthogonal both to v1 and v3
                v2_aux = cross(v3,v1);
                v2 = v2_aux/norm(v2_aux);
                % save v1 v2 and v3 in the global coords' directions
                nodal_systems(:,:,n) = [v1 v2 v3];
            end
        end
        function L = integral_along_surface(mesh,dofs_per_node, ...
                dofs_per_ele,s_nodes,fun_in)
            % L = integral_along_surface(mesh,s_nodes,fun_in)
            % s_nodes [n x 1]: nodes in the surface
            % fun_in [FHandle](element,surface_c,surface_value)
            ele_ids = mesh.unique_eles(s_nodes);
            ele_s = mesh.ele_surfaces(s_nodes,ele_ids);
            [~, s_coord, s_val] = mesh.ele_type.surfaces();
            function L_out = aux_fun(element)
                surfaces = ele_s{element.id == ele_ids};    % Element's surface #
                L_out = [];
                for s = surfaces
                    if isempty(L_out)
                        L_out = fun_in(element,s_coord(s),s_val(s));
                    else
                        L_out = L_out + fun_in(element,s_coord(s),s_val(s));
                    end
                end
            end
            L = mesh.assembly_vector_along(dofs_per_node,dofs_per_ele, ...
                                                        ele_ids,@aux_fun);
        end
        function L = assembly_vector_along(mesh,dofs_per_node,dofs_per_ele, ...
                ele_ids,fun_in)
            % L = assembly_vector_along(mesh,dofs_per_node,dofs_per_ele, ...
            %                                    ele_ids,fun_in)
            % K [n_dofs x n_dofs]
            % dofs_per_node [Int]
            % dofs_per_ele [Int]
            % fun_in [Fhandle] f(element) -> [dofs_per_ele x 1]
            % fun_in follows convention [node_dofs ele_dofs]
            require(all(~mod([dofs_per_node,dofs_per_ele],1)), ...
                'ArgumentError: dofs should be integers');
            L = sparse(mesh.n_dofs(dofs_per_node,dofs_per_ele),1);
            % Loop through elements
            for e = ele_ids
                ele = mesh.ele(e);
                dofs = mesh.all_eles_dofs(dofs_per_node,dofs_per_ele,e);
                L(dofs) = L(dofs) + fun_in(ele);
            end
        end
        function new_mesh = add_point_mass(mesh,nodes_with_mass,mass_values)
            % add_point_mass(mesh,nodes_with_mass,mass_val)
            % Getter
            require(all(size(nodes_with_mass) == size(mass_values)), ...
                'ArgumentError: nodes_with_mass and mass_val should be the same size');
            mesh.nodes_with_mass = nodes_with_mass;
            mesh.mass_values     = mass_values;
            new_mesh = mesh;
        end
        function L = assembly_vector(mesh,dofs_per_node,dofs_per_ele,fun_in)
            % L = assembly_vector(mesh,dofs_per_node,dofs_per_ele,fun_in)
            % L [n_dofs x 1]
            % dofs_per_node [Int]
            % dofs_per_ele [Int]
            % fun_in [Fhandle] f(element) -> [dofs_per_ele x 1]
            % fun_in follows convention [node_dofs ele_dofs]
            L = assembly_vector_along(mesh,dofs_per_node,dofs_per_ele, ...
                                                    1:mesh.n_ele,fun_in);
        end
        function K = assembly_matrix(mesh,dofs_per_node,dofs_per_ele,fun_in)
            % K = assembly_matrix(mesh,dofs_per_node,dofs_per_ele,fun_in)
            % K [n_dofs x n_dofs]
            % dofs_per_node [Int]
            % dofs_per_ele [Int]
            % fun_in [Fhandle] f(element) -> [dofs_per_ele x dofs_per_ele]
            % fun_in follows convention [node_dofs ele_dofs]
            require(all(~mod([dofs_per_node,dofs_per_ele],1)), ...
                'ArgumentError: dofs should be integers');
            n_dofs = mesh.n_dofs(dofs_per_node,dofs_per_ele);
            K = sparse(n_dofs,n_dofs);
            % Loop through elements
            for e = 1:mesh.n_ele
                ele = mesh.ele(e);
                dofs = mesh.all_eles_dofs(dofs_per_node,dofs_per_ele,e);
                K(dofs,dofs) = K(dofs,dofs) + fun_in(ele);
            end
        end
        function nodes = find_nodes(mesh,condition)
            % nodes = find_nodes(fun_in)
            % nodes [n_nodes x 1][Bool]
            % condition [FHandle][(x,y,z) -> Bool]
            % Returns a logical matrix where the true-valued nodes satisfy
            % the condition
            nodes = false(mesh.n_nodes,1);
            for n = 1:mesh.n_nodes
                r = mesh.coords(n,:);   % r = [x,y,z]
                nodes(n) = condition(r(1),r(2),r(3));
            end
        end
        %% DOF helpers
        function s_dofs = slave_dofs(mesh,dofs_per_node,dofs_per_ele,master_dofs)
            % s_dofs = slave_dofs(mesh,master_dofs)
            % Returns all the dofs that are not master_dofs
            all_dofs = (1:mesh.n_dofs(dofs_per_node,dofs_per_ele))';
            s_dofs = setdiff(all_dofs,master_dofs);
        end
        function dofs = master_ele_dofs(mesh,dofs_per_node,dofs_per_ele, ...
                                                        ele_ids,dofs_type)
            % dofs = master_ele_dofs(mesh,dofs_per_node,dofs_per_ele, ...
            %                        ele_ids,dofs_type)
            % dofs [Int]: Master Dof ids
            % dofs_per_node [Int]
            % dofs_per_ele [Int]
            % ele_ids: [E x 1][Int]: Vector containing the element ids
            %   accordint to the mesh.
            % dofs_type [e_types x 1][Int]: Vector containg the type of
            %   dofs that are neede
            % Finds the dofs corresponding to the elements in ele_ids, and
            % then filters out by dof_types
            require(all(dofs_type <= dofs_per_ele), ...
                'ArgumentError: all values in dofs_type should be smaller than dofs_per_ele');
            require(all(dofs_type == unique(dofs_type)), ...
                'ArgumentError: dofs_type should have unique values');
            require(dofs_per_ele > 0, ...
                'ArgumentError: dof_per_ele should be greater than zero');
            E = length(ele_ids);                % Number of elements to consider
            e_types = length(dofs_type);    % Number of dofs_per_ele to consider
            dofs = zeros(E*e_types,1);
            l = mesh.last_node_dof(dofs_per_node);
            % Loop through the selected elements
            for e = 1:E
                % Find all the element's dofs
                all_ele_dofs = l + index_range(dofs_per_ele,ele_ids(e));
                % Keep only the interesting ones, specified by dofs_type
                dofs(index_range(e_types,e)) = all_ele_dofs(dofs_type);
            end
        end
        function dofs = master_node_dofs(mesh,dofs_per_node,node_ids, ...
                                                        dofs_type)
            % dofs = master_node_dofs(mesh,dofs_per_node,dofs_per_ele,node_ids,dof_type)
            % Finds the dofs corresponding to the nodes in node_ids, and
            % then filters out the dof_types
            % dofs [Int]: Master Dof ids
            % dofs_per_node [Int]
            % node_ids: [N x 1][Int]: Vector containing the node ids
            %   accordint to the mesh.
            % dofs_type [n_types x 1][Int]: Vector containg the type of
            %   dofs that are neede
            % Finds the dofs corresponding to the nodes in node_ids, and
            % then filters out by dof_types
            require(all(dofs_type <= dofs_per_node), ...
                'ArgumentError: all values in dofs_type should be smaller than dofs_per_ele');
            require(all(dofs_type == unique(dofs_type)), ...
                'ArgumentError: dofs_type should have unique values');
            require(dofs_per_node > 0, ...
                'ArgumentError: dof_per_ele should be greater than zero');
            N = length(node_ids);           % Number of nodes to consider
            n_types = length(dofs_type);    % Number of dofs_per_node to consider
            dofs = zeros(N*n_types,1);      
            % Loop through selected nodes
            for n = 1:N
                % Find all the node's dofs
                all_node_dofs = index_range(dofs_per_node,node_ids(n));
                % Keep only the interesting ones, specified by dofs_type
                dofs(index_range(n_types,n)) = all_node_dofs(dofs_type);
            end
        end
        function dofs = all_eles_dofs(mesh,dofs_per_node,dofs_per_ele,ele_id)
            % eles_dofs(mesh,dofs_per_node,dofs_per_ele,ele_id)
            % Computes the complete Dof list for a certain element
            node_dofs = mesh.ele_node_dofs(dofs_per_node,ele_id);
            ele_dofs = mesh.eles_dofs(dofs_per_node,dofs_per_ele,ele_id);
            dofs = [node_dofs; ele_dofs];
        end
        function dofs = eles_dofs(mesh,dofs_per_node,dofs_per_ele,ele_id)
            % dofs(mesh,dofs_per_node,dofs_per_ele,ele_id)
            % Finds the dofs that correspond to the element but not to the 
            % nodes of the element
            dofs = mesh.last_node_dof(dofs_per_node) + ...
                index_range(dofs_per_ele,ele_id)';
        end
        function dofs = ele_node_dofs(mesh,dofs_per_node,ele_id)
            % dofs = ele_node_dofs(mesh,dofs_per_node,ele_id)
            % Finds the dofs that correspond only to the nodes of the
            % element
            nodes = mesh.ele_nodes(ele_id);
            % Loop through nodes and add their dofs
            dofs = zeros(mesh.nodes_per_ele*dofs_per_node,1);
            for i = 1:mesh.nodes_per_ele
                dofs(index_range(dofs_per_node,i)) = ... 
                    mesh.node_dofs(dofs_per_node,nodes(i));
            end
        end
        function dofs = node_dofs(mesh,dofs_per_node,node_id)
            % dofs = node_dofs(mesh,dofs_per_node,node_id)
            % Finds the dofs that correspond only to the node_id
            % Could be Static, it's not for consistency to other similar
            % methods
            dofs = index_range(dofs_per_node,node_id);
        end
        function out = n_dofs(mesh,dofs_per_node,dofs_per_ele)
            % out = n_dofs(mesh,dofs_per_node,dofs_per_ele)
            % Total number of DOFs
            out = mesh.last_node_dof(dofs_per_node) + ...
                mesh.n_ele * dofs_per_ele;
        end
        function out = last_node_dof(mesh,dofs_per_node)
            % out = last_node_dof(mesh,dofs_per_node)
            out = mesh.n_nodes*dofs_per_node;
        end
        function out = all_node_dofs(mesh,dofs_per_node)
            out = 1:mesh.last_node_dof(dofs_per_node);
        end
        function out = all_element_dofs(mesh,dofs_per_node,dofs_per_ele)
            last_ele_dofs = mesh.all_eles_dofs(dofs_per_node,dofs_per_ele,mesh.n_ele);
            out = (mesh.last_node_dof(dofs_per_node)+1):last_ele_dofs(end);
        end
        %% Element Helpers
        function ele = ele(mesh,ele_id)
            % ele = ele(mesh,ele_id)
            % ele [Element]
            % Create element with ID ele_id from the mesh
            nodes = mesh.ele_nodes(ele_id);
            if EleType.is_shell(mesh.ele_type)
                laminate = mesh.laminates(mesh.laminate_ids(ele_id));
                ele = Element(ele_id,mesh.ele_type,laminate, ...
                    mesh.coords(nodes,:),mesh.normals(:,:,nodes),mesh.thickness(nodes));
            else
                ele = Element(ele_id,mesh.ele_type,mesh.coords(nodes,:), ...
                                [],[]);
            end
        end
        function lam = ele_laminate(mesh,ele_id)
            if length(mesh.laminate_list) == 1
                lam = mesh.laminate_list;
            else
                lam = mesh.laminate_list(ele_id);
            end
        end
        function node_ids = ele_nodes(mesh,ele_id)
            % nodes_ids = ele_nodes(mesh,ele_id)
            % nodes_ids [1xnodes_per_element]
            % Returns all the IDs corresponding to ele_id
            node_ids = mesh.connect(ele_id,:);
        end
        function ele_ids = node_eles(mesh,node_id)
            % ele_ids = node_eles(mesh,node_id)
            % node_id [n_nodes x 1][Int]: Ids of the nodes.
            % ele_ids {n_nodes x 1}[Int]: Ids of the elements that contain
            % those nodes.
            % A cell array is used because each node might return a
            % different number of elements
            n_nodes = length(node_id);
            ele_ids = cell(n_nodes,1);
            for n = 1:n_nodes
                ele_ids{n} = find(any((mesh.connect == node_id(n))')')';
            end
        end
        function ele_ids = unique_eles(mesh,node_ids)
            % ele_ids = unique_elements(mesh,node_ids)
            % node_id [n_nodes x 1][Int]: Ids of the nodes.
            % ele_ids [? x 1][Int]: Ids of the elements that contain those nodes.
            node_eles = mesh.node_eles(node_ids);
            ele_ids = [];
            for n = 1:length(node_eles)
                ele_ids = [ele_ids node_eles{n}];
            end
            ele_ids = unique(ele_ids);
        end
        function ele_surfaces = ele_surfaces(mesh,nodes,ele_ids)
            % ele_surfaces = ele_surfaces(mesh,nodes,ele_ids)
            % ele_surfaces {n_ele_ids x 1}[Int] the surfaces (between 1
            % and 6) of each element that contain the nodes.
            % nodes: [n x 1][Int] node_ids of the nodes in the global surface
            n_ele = length(ele_ids);    % Number of elements
            ele_surfaces = cell(n_ele,1);
            % surfaces contains data from the element type
            surfaces = mesh.ele_type.surfaces();
            % Loop through all elements
            for e = 1:n_ele
                % Get all of its nodes
                ele_nodes = mesh.ele_nodes(ele_ids(e));
                % Loop through the nodes
                ele_surface_nodes = false(1,length(ele_nodes));
                for n = 1:length(ele_nodes)
                    % Check if the node is in the global surface
                    ele_surface_nodes(n) = any(nodes == ele_nodes(n));
                end
                % Those who were, are stored in ele_surface_nodes
                ele_surface_nodes = sort(find(ele_surface_nodes));
                ele_surface = [];
                % Check which surfaces contain ele_surface_nodes
                for s = 1:size(surfaces,2)
                    if all(ele_surface_nodes == surfaces(s,:))
                        ele_surface = [ele_surface s];
                    end
                end
                ele_surfaces{e} = ele_surface;
            end
        end
        %% Plotting
        function plot(mesh)
            % Nodes per element (to be plotted)
            nodes_per_ele = 4;  % HARDCODED FOR RECTANGULAR SHELL ELEMENTS                              
            % Initialization of the required matrices
            X = zeros(nodes_per_ele,mesh.n_ele); 
            Y = zeros(nodes_per_ele,mesh.n_ele); 
            Z = zeros(nodes_per_ele,mesh.n_ele);
            C = char(mesh.n_ele,1);
            cc = ['w','r'];
            for e = 1:mesh.n_ele
                e_nodes = mesh.ele_nodes(e);    % nodes for ele with id = e
                C(e) = cc(mesh.laminate_ids(e));
                for n = 1:nodes_per_ele
                    X(n,e) = mesh.coords(e_nodes(n),1); % extract x value of the node
                    Y(n,e) = mesh.coords(e_nodes(n),2); % extract y value of the node
                    Z(n,e) = mesh.coords(e_nodes(n),3); % extract z value of the node
                end
            end
            % Plotting the FEM mesh, diaplay Node numbers and Element numbers
            f1 = figure ;
            set(f1,'name','Mesh','numbertitle','off') ;
            hold on
            for e = 1:mesh.n_ele
                patch(X(:,e),Y(:,e),Z(:,e),C(e));
            end
            e_nodes = mesh.connect(:,1:end)';
%             for n = 1:mesh.n_ele
%                 text(X(:,n),Y(:,n),Z(:,n),int2str(e_nodes(1:4,n)), ...
%                                                 'fontsize',8,'color','k');
%                 text(sum(X(:,n))/4,sum(Y(:,n))/4,sum(Z(:,n))/4,int2str(n), ...
%                                                 'fontsize',10,'color','r') ;
%             end
            if ~isempty(mesh.nodes_with_mass)
                nodes = mesh.nodes_with_mass;
                plot3(mesh.coords(nodes,1),mesh.coords(nodes,2), ...
                                        mesh.coords(nodes,3),'b.')
            end
            hold off
        end
        function plot_displacements(mesh,scale,displacements)
            % plot_displacements(mesh,scale,displacements)
            % displacements [n_nodes x 3]
            % scale [Float]: the bigger scale, the bigger displacements look
            aux_mesh = mesh;
            aux_mesh.coords = mesh.coords + scale*displacements;
            aux_mesh.plot();
        end
        function plot_vector_field(mesh,vector_field,options)
            % plot_vector_field(mesh,vector_field,option)  
            % Wrapper for quiver3
            % vector_field [n_nodes x 3][Float]
            % options: MATLAB plotting options, i.e. 'k-'
            quiver3(mesh.coords(:,1),mesh.coords(:,2),mesh.coords(:,3), ...
            	vector_field(:,1),vector_field(:,2),vector_field(:,3), ...
                options);
        end
        %% Dependent Properties
        function out = get.n_nodes(mesh)
            % out = n_nodes(mesh)
            % Number of Nodes
            out = size(mesh.coords,1);
        end
        function out = get.n_ele(mesh)
            % out = n_nodes(mesh)
            % Number of Elements
            out = size(mesh.connect,1);
        end
        function out = get.nodes_per_ele(mesh)
            % out = nodes_per_ele(mesh)
            out = size(mesh.connect,2);
        end
    end
end