% *************************************************************************
% MeshGen - gera malha de nós ligado a todos os demais nos 
% 
%   Gera Matriz incidencia
%   Gera Matriz coordenada
%   Gera matriz de deslocamentos impostos (cantilever)
%   Gera Vetor de carga (elemento central da ponta livre)
%   
%       Incidencia elemental
%       inci =[ (n. elemento) (tipo ) (n.tab mat) (n. tab geo)  noi  noj ]  
%
%       Coordenadas
%       coord = [x y]
%
%       Deslocamentos impostos
%       BC = [ nó   (u=1 ou v=2)  valor] 
%
%       Forças impostas
%       Loads = [ nó   (fx=1   ou fy=2)  magnitude(mag)]
%
%    struc_type -- type of Boundary conditions for Truss Structure - 
%     Ground  structure  model 
%           1 Cantilever Beam 
%           2 Simply supported beam
%      mesh_type 
%         1 full mesh - all nodes connected 
%         2 local mesh
% 
%
%   Author: Renato Pavanello
%   Copyright (c) 2023 by The DMC/FEM/UNICAMP.
%**************************************************************************

function [inci,coord,BC,Loads,nel,nnos] = meshGen(Lx,Ly,nodesx,nodesy,mag,kind,mesh_type)
%% Matriz Inicidencia e de Coordenadas
dx = (Lx/(nodesx-1));   
dy = (Ly/(nodesy-1));

nnos = nodesx*nodesy;  % total de nos
coord=zeros(nnos,2);

k = 1; %posicao da linha a ser atualizada
x_0 = 0; y_0 = 0; %coordenada do no 1

for i=1:nodesx
    for j=1:nodesy
        %matriz de coordenadas
        coord(k,1) = x_0;
        coord(k,2) = y_0;
        k = k+1;
        y_0 = y_0 +dy;
    end
    x_0 = x_0 + dx; y_0 = 0;
end

switch mesh_type
    case 1 
        %%  Full mesh - todos nos conectados entre sí. 
        nel = (nnos^2-nnos)/2;  % total de elementos
        v = 1:nnos; % vetor de nos
        b = 1:nel; % vetor de elementos
        c = nchoosek(v,2); % vetor combinacao
        inci= [b',c];
        one = ones(nel,2);
        inci = [inci(:,1) one inci(:,1) inci(:,2:end)]; %adicionando os outros elementos da matriz incidencia (material, geometria)

    case 2 
        %% malha local , sem superposições
         nel=1;    
        
         % barras verticais
        for i=1:nodesx   
        for j=1:nodesy-1
            inode=j+(i-1)*nodesy;
            inci(nel,:)=[ nel ,1 ,1 ,nel ,inode , inode+1]; nel=nel+1;
        end
        end        
        
        % barras horizontais
        for j=1:nodesy   
        for i=1:nodesx-1
            inode=j+(i-1)*nodesy;
            jnode=j+(i)*nodesy;
            inci(nel,:)=[ nel ,1 ,1 ,nel ,inode , jnode]; nel=nel+1;
        end
        end        
        
        % barras cruzadas
        for j=1:nodesy-1   
        for i=1:nodesx-1
            inode=j+(i-1)*nodesy;
            jnode=j+(i)*nodesy+1;
            inci(nel,:)=[ nel ,1 ,1 ,nel ,inode , jnode]; nel=nel+1;
            inci(nel,:)=[ nel ,1 ,1 ,nel ,inode+1 , jnode-1]; nel=nel+1;
         end
        end        
        
    otherwise               
        disp('*******************************************************');
        disp('*********  mesh_type not implemented  STOP program **********');
        disp('*******************************************************');
        pause;
end

nel=size(inci,1);


%% TIPO DE VIGA
no = 1;
u_v = 1; % u=1 v=2
switch kind
    case 1
        %% Displacements Constraints (cantilever beam )
        BC = zeros(nodesy*2,3);    %  Displacement constraint 

        for i = 1:2*nodesy
                if u_v <=2
                    BC(i,1) = no;
                    BC(i,2) = u_v;
                    u_v = u_v + 1;
                else
                    u_v = 1;
                    BC(i,2) = u_v;
                    no = no+1;
                    BC(i,1) = no;
                    u_v = u_v + 1;
                end

            end

        %% Load Case  (cantilever Beam)
        %  Loads = [ no   (fx=1   ou fy=2)  valor]
        %caso nao haja nó central, divide forca entre os dois nos centrais    
        if mod(nodesy,2)==0
            mag = mag/2;
            Loads=[nnos-nodesy/2+1 2 mag;
                   nnos-nodesy/2 2 mag]; 
        else
                Loads=[nnos-round(nodesy/2)+1 2 mag];      
        end
            
    case 2
        %% Deslocamentos impostos (viga biapoiada)
        BC = zeros(4,3);
       
           for i = 1:4

                if u_v <=2
                    BC(i,1) = no;
                    BC(i,2) = u_v;
                    u_v = u_v + 1;
                else
                    u_v = 1;
                    BC(i,2) = u_v;
                    no = (nodesx-1)*nodesy+1;
                    BC(i,1) = no;
                    u_v = u_v + 1;
                end

           end
         
        %% Carga imposta (viga biapoiada) 
        %caso nao haja nó central, divide forca entre os dois nos centrais    
        if mod(nodesx,2)==0
            mag = mag/2;
            Loads=[nnos/2+1 2 mag;
                   nnos/2-nodesy+1 2 mag]; 
        else
                Loads=[round((nnos-1)/2) 2 mag];      
        end
            
        
end


end