%%
% ---------------------------------------------------------------------------
%

classdef FEMElement < handle
    
  
    properties(GetAccess=private)
        v % v1=0.1;v2=0.4;v3=0.49;v4=0.499;
        x %Ortsvektor (4x2 Matrix)
         % Verschiebung x Richtung
         % & Verschiebung y Richtung
        Kraft        
        u 
        
       
    end
    %%
    
    methods 
        function el=FEMElement(Querkontraktionszahl,OrtsZeilenMatrix4x2)
        el.v=Querkontraktionszahl;
       
        el.x=OrtsZeilenMatrix4x2;
       %        4   1   %%  Knotennummern, Kraft setzt an Knoten 1 an
       %        3   2   %%  Knoten 2 und 3 = Fest, an Knoten 1 grift im
       %        45° Winkel eine Kraft nach oben rechts an.
       
       
        end
        
        %%
       
        
        function BerechnungVerschiebung(el)

            syms r s
           
           
            
         %  Summe der X Werte
           
                xdach=transpose(reshape(el.x(:,1),2,2));
                       xdachsumme=sum(sum(xdach.*[(1-r)*(1-s),(1-r)*(1+s);...
                                              (1+r)*(1+s),(1+r)*(1-s)]));
                                          
           
            
            % Summe der Y Werte
            ydach=transpose(reshape(el.x(:,2),2,2));
                        ydachsumme=sum(sum(ydach.*[(1-r)*(1-s),(1-r)*(1+s);...
                                                (1+r)*(1+s),(1+r)*(1-s)]));
            %JacobiMatrix
            JacMatrix=jacobian([xdachsumme, ydachsumme]);
         
            JacMatrixinv=inv(JacMatrix);
            % 2x4
            DNiDxi=0.25*JacMatrixinv*[1+s 1-s -1+s -1-s;...
                                      1+r 1-r -1+r -1-r]; %#ok<MINV> %
                                

            
            % B=3x8 Matrix Differentialoperatormatrix
            B=[ DNiDxi(1,1), 0         ,DNiDxi(1,2), 0            ,DNiDxi(1,3), 0         ,DNiDxi(1,4), 0         ;
                0          ,DNiDxi(2,1),   0       ,DNiDxi(2,2)   , 0         ,DNiDxi(2,3),   0       ,DNiDxi(2,4);
                DNiDxi(2,1),DNiDxi(1,1),DNiDxi(2,2),DNiDxi(1,2)   ,DNiDxi(2,3),DNiDxi(1,3),DNiDxi(2,4),DNiDxi(1,4)];
            
           
            % C = 3x3 Matrix
                C=[1 el.v 0;el.v 1 0;0 0 (1-el.v)/2]/(1-el.v^2);
               
            % Integration Element

                Integrant=B'*C*B;
           
            % Randbedingungen
            Integrant(3:6,:)=[];
            Integrant(:,3:6)=[];
          
        
        
            KElement=int(Integrant,r);
            a=subs(KElement,r,1);
            b=subs(KElement,r,-1);
            Kelementneu1=a-b;
            
            Kelementneu2=int(Kelementneu1,s);
            c=subs(Kelementneu2,s,1);
            d=subs(Kelementneu2,s,1);
            
            KElementend=c-d;
            
                      
            KElementinv=KElementend^-1;
        
            
            % Kraftvektor
            el.Kraft=[1/sqrt(2),1/sqrt(2),0,0]';
            
            % u = K^-1 * F
            el.u=KElementinv*el.Kraft;
            
            % el.x=el.x+el.u;
        end
    end
end
% -------------------------------------------------------------------------