function BMRSMES2Dtest
close all
obsz = 1;
nh = 6;
a = 1;
b = 1;
E = 20;
v = 0.1;
q = 1;
P = 1;
obc = 1; % 1 - ci¹g³e, 2 - si³a skupiona
eps = 10^-12;
koniec = 0;
siat = 1;
adap = 2; % 1 - regularne, 2 - nieregularne adaptacyjne
siatmax = 25;
if obc==2
    siatmax = 27;
end
if obsz==3
    siatmax = 1;
end
if adap==1
    siatmax = 10;
end
Rsr = zeros(100,3);
Rsrx = zeros(100,3);
Rsry = zeros(100,3);
Rmax = zeros(100,3);
NSS = zeros(100,1);
Ex = zeros(100,3);
Ey = zeros(100,3);
oppop = 1;
if obsz==3
    obc = 1;
end
if obc==1
    P = 0;
end
if obc==2
    q = 0;
end
while (koniec == 0)
    if obsz==3 & siat==1
        p = fopen('pp.mmg','rb');
        n = fscanf(p,'%d',[1 1]);
        X = zeros(n,1);
        Y = zeros(n,1);
        WB = zeros(n,1);
        WS = zeros(n,1);
        POP = zeros(n,1);
        NAS = zeros(n,1);
        for i=1:n
            X(i) = fscanf(p,'%f',[1 1]);
            Y(i) = fscanf(p,'%f',[1 1]);
            WB(i) = fscanf(p,'%d',[1 1]);
            z = fscanf(p,'%f',[1 1]);
            if WB(i)==1
                n1 = fscanf(p,'%d',[1 1]);
                n2 = fscanf(p,'%d',[1 1]);
                POP(i) = n1 + 1;
                NAS(i) = n2 + 1;
                if Y(i)<=3
                    WB(i) = 2;
                end
                if Y(i)<=-3
                    WB(i) = 2;
                end
                if Y(i)>=3
                    WB(i) = 3;
                end
            end
        end
        iltr = fscanf(p,'%d',[1 1]);
        T = fscanf(p,'%d',[3 iltr]);
        T = T' + 1;
        fclose(p);
    end
    if obsz<=2 & (siat==1 | adap==1)
        n = nh*nh;
        hx = a/(nh-1);
        hy = b/(nh-1);
        X = zeros(n,1);
        Y = zeros(n,1);
        WB = zeros(n,1);
        WS = zeros(n,1);
        il = 0;
        x = 0;
        y = 0;
        amp = 0.5;
        for i=1:nh
            for j=1:nh
                il=il+1;
                X(il) = x;
                Y(il) = y;
                if abs(x)<eps | abs(y)<eps | abs(x-a)<eps | abs(y-b)<eps
                    if WB(il)==3
                        WS(il) = 1;
                    end
                    WB(il) = 1; % brzeg swobodny
                end
                if abs(y-b)<eps
                    if WB(il)>0
                        WS(il) = 1;
                    end
                    WB(il) = 3; % obci¹¿enie
                    if obc==2
                        WB(il) = 1;
                    end
                end
                if abs(x)<eps
                    if WB(il)==3
                        WS(il) = 1;
                    end
                    WB(il) = 2; % zamocowanie
                end
                x = x + hx;
                if obsz==2 & WB(il)==0
                    X(il) = X(il) + (amp/2-amp*rand(1,1))*hx;
                    Y(il) = Y(il) + (amp/2-amp*rand(1,1))*hy;
                end
            end
            x = 0;
            y = y + hy;
        end
    end
    
    Nss = 2*n;
    il = 25;
    if il>n
        il = n;
    end
    %il = n; % !!!!!!!!!!!!!!!!!!!
    
    ilHO = 0;
    HO = 0;
    edop = 10^-3;
    HOsr = 1;
    WEK = zeros(n,2);
    WEZ = 0;
    dmin = 6;
    for i=1:n
        d = sqrt( (X(i)-0)^2 + (Y(i)-3)^2 );
        if d<dmin
            d = dmin;
            WEZ = i;
        end
    end
    while (HO<=ilHO & HOsr>edop)
        A = zeros(Nss,Nss);
        B = zeros(Nss,1);
        num = 0;
        for i=1:n
            if WB(i)==0
                G = gwiazda(X,Y,X(i),Y(i));
                ild = il;
                if obsz==4
                    NR = 0;
                    for j=1:iltr                       
                        for k=1:3
                            if i==T(j,k)
                                NR = j;
                            end
                        end                       
                    end
                    G = [T(NR,1) 0; T(NR,2) 0; T(NR,3) 0];
                    ild = 3;  
                end
                [M,K] = mwls(X,Y,X(i),Y(i),G,ild,2);   
                num = num + 1;
                U = zeros(il,1);
                for j=1:ild
                    nr = G(j,1);
                    A(num,2*nr-1) = M(4,j) + M(6,j) * (1-v)/2;
                    A(num,2*nr) = v * M(5,j) + M(5,j) * (1-v)/2;
                    A(num+1,2*nr-1) = v * M(5,j) + M(5,j) * (1-v)/2;
                    A(num+1,2*nr) = M(6,j) + M(4,j) * (1-v)/2;
                end
                if HO>0
                    delx = K*Wx1(i,7:15)';
                    dely = K*Wy1(i,7:15)';
                    B(num) = delx(4) + delx(6) * (1-v)*0.5 + v*dely(5) + dely(5)*(1-v)*0.5;
                    B(num+1) = v*delx(5) + delx(5)*(1-v)*0.5 + dely(6) + dely(4)*(1-v)*0.5;
                end
                num = num + 1;
            end
            if WB(i)>0
                Gw = zeros(n,2);
                ile = 0;
                for j=1:n
                    if (i~=j & (WB(j)>0 | WS(j)==1))
                        ile = ile + 1;
                        d = sqrt((X(j) - X(i))^2 + (Y(j) - Y(i))^2);
                        Gw(ile,:) = [j d];
                    end
                end          
                Gw(ile+1:end,:) = [];
                Gw = sortrows(Gw,2);
                hs = (Gw(1,2) + Gw(2,2))/2;
                if obsz<=2
                    if abs(X(i)-a)<eps | abs(X(i)<eps);
                        hs = Gw(1,2)/2;
                    end
                end
            end
            if WB(i)==3 | WB(i)==1
                ilb = il;
                Gd = gwiazda(X,Y,X(i),Y(i));  
                G = Gd;       
                if obsz==4
                    NR = 0;
                    for j=1:iltr
                        for k=1:3
                            if i==T(j,k)
                                NR = j;
                            end
                        end
                    end
                    G = [T(NR,1) 0; T(NR,2) 0; T(NR,3) 0];
                    ilb = 3;  
                end                               
                [M,K] = mwls(X,Y,X(i),Y(i),G,ilb,2);
                num = num + 1;
                x = X(i);
                y = Y(i);
                if abs(y)<eps
                    Wek = [0 -1];
                end  
                if abs(x-a)<eps
                    Wek = [1 0];
                end
                if abs(y-b)<eps
                    Wek = [0 1];
                end                   
                if abs(x)<eps
                    Wek = [-1 0];
                end
                if obc<=2    
                    nr1 = Gw(1,1);
                    nr2 = Gw(2,1);
                    Wek = [abs(Y(nr1)-Y(nr2)) abs(X(nr1)-X(nr2))];
                    Wek = Wek/norm(Wek);
                    if X(i)<eps
                        Wek(1) = -Wek(1);
                    end
                    if Y(i)<eps
                        Wek(2) = -Wek(2);
                    end
                    if obsz<=3
                        if obsz==3 & siat==1
                            nr1 = POP(i);
                            nr2 = NAS(i);
                        end
                        sk = 0.4;
                        Wekx = [X(nr1)-X(nr2) Y(nr1)-Y(nr2)];
                        Weky1 = [Wekx(2) -Wekx(1)];
                        Weky2 = [-Wekx(2) Wekx(1)];
                        Weky1 = Weky1/norm(Weky1);
                        Weky2 = Weky2/norm(Weky2);
                        Weky1s = sk*Weky1/norm(Weky1);
                        Weky2s = sk*Weky2/norm(Weky2);
                        P1 = [X(nr1) Y(nr1); X(nr1)+Weky1s(1) Y(nr1)+Weky1s(2); X(nr2)+Weky1s(1) Y(nr2)+Weky1s(2); X(nr2) Y(nr2)];
                        P2 = [X(nr1) Y(nr1); X(nr1)+Weky2s(1) Y(nr1)+Weky2s(2); X(nr2)+Weky2s(1) Y(nr2)+Weky2s(2); X(nr2) Y(nr2)];
                        il1 = prostokat(P1,X,Y,Gd(1:9,:));
                        il2 = prostokat(P2,X,Y,Gd(1:9,:));
                        Wek = Weky1;
                        if il1>il2
                            Wek = Weky2;
                        end                    
                    end
                    WEK(i,:) = Wek;                    
                end
                for j=1:ilb
                    nr = G(j,1);
                    A(num,2*nr-1) = (E/(1-v^2))*M(2,j)*Wek(1) + (E/(1-v^2))*((1-v)/2)*M(3,j)*Wek(2);
                    A(num,2*nr) = (E/(1-v^2))*M(3,j)*v*Wek(1) + (E/(1-v^2))*((1-v)/2)*M(2,j)*Wek(2);
                    A(num+1,2*nr-1) = (E/(1-v^2))*((1-v)/2)*M(3,j)*Wek(1) + (E/(1-v^2))*v*M(2,j)*Wek(2);
                    A(num+1,2*nr) = (E/(1-v^2))*((1-v)/2)*M(2,j)*Wek(1) + (E/(1-v^2))*M(3,j)*Wek(2);
                end
                if HO>0
                    delx = K*Wx1(i,7:15)';
                    dely = K*Wy1(i,7:15)';
                    B(num) = (E/(1-v^2))*delx(2)*Wek(1) + (E/(1-v^2))*((1-v)/2)*delx(3)*Wek(2) + (E/(1-v^2))*dely(3)*v*Wek(1) + (E/(1-v^2))*((1-v)/2)*dely(2)*Wek(2);
                    B(num+1) = (E/(1-v^2))*((1-v)/2)*delx(3)*Wek(1) + (E/(1-v^2))*v*delx(2)*Wek(2) + (E/(1-v^2))*((1-v)/2)*dely(2)*Wek(1) + (E/(1-v^2))*dely(3)*Wek(2);
                end
                if WB(i)==3
                    B(num) = B(num) + 0;
                    B(num+1) = B(num+1) - q;
                end                
                if WB(i)>0 & abs(X(i) - a)<eps & abs(Y(i)-b)<eps                  
                    B(num+1) = B(num+1)-P/hs;
                end
                num = num + 1;
            end
            if WB(i)==2
                num = num + 1;
                A(num,2*i-1) = 1;
                A(num+1,2*i) = 1;
                num = num + 1;
            end
        end
        cond(A)
        det(A)
        U1 = A^-1 * B;        
        if HO>0
            Ux1p = Ux1;
            Uy1p = Uy1;
        end
        Ux1 = U1(1:2:Nss);
        Uy1 = U1(2:2:Nss);
        
        if HO<ilHO
            Wx1 = zeros(n,15);
            Wy1 = zeros(n,15);
            ilp = n;
            stp = 4;
            Bx = zeros(ilp,1);
            By = zeros(ilp,1);
            n1 = [2 3 2 3 4 6 4 6 5];
            n2 = [4 4 5 5 4 4 5 5 5];
            for i=1:n
                G = gwiazda(X,Y,X(i),Y(i));
                [M,K] = mwls(X,Y,X(i),Y(i),G,ilp,stp);
                for j=1:ilp
                    nr = G(j,1);
                    Bx(j) = Ux1(nr);
                    By(j) = Uy1(nr);
                end
                wx = M*Bx;
                wy = M*By;
                delx = zeros(6,1);
                dely = zeros(6,1);
                if HO>0 & stp==2
                    delx = K*Wx1(i,7:15)';
                    dely = K*Wy1(i,7:15)';
                end
                Wx1(i,1) = Ux1(i);
                Wy1(i,1) = Uy1(i);
                Wx1(i,[2:6]) = wx(2:6)-delx(2:6);
                Wy1(i,[2:6]) = wy(2:6)-dely(2:6);
            end
            if stp<4
                for i=1:n
                    G = gwiazda(X,Y,X(i),Y(i));
                    M = mwls(X,Y,X(i),Y(i),G,il,2);
                    for k=1:9
                        for j=1:il
                            nr = G(j,1);
                            Wx1(i,k+6) = Wx1(i,k+6) + M(n1(k),j) * Wx1(nr,n2(k));
                            Wy1(i,k+6) = Wy1(i,k+6) + M(n1(k),j) * Wy1(nr,n2(k));
                        end
                    end
                end
            end
            if HO>0
                HOsr = norm(Ux1p - Ux1)/norm(Ux1) + norm(Uy1p - Uy1)/norm(Uy1);
            end
        end
        HO = HO + 1;        
    end    
    
    if obsz<=2
        T = delaunay(X,Y);
        iltr = length(T);      
    end
    N = iltr;
        
    K = zeros(Nss,Nss);
    F = zeros(Nss,1);
    D = (E/(1-v^2))*[1 v 0;v 1 0; 0 0 (1-v)/2];
    for i=1:N
        nr = [T(i,1) T(i,2) T(i,3)];
        Xt = [X(nr(1)) X(nr(2)) X(nr(3))];
        Yt = [Y(nr(1)) Y(nr(2)) Y(nr(3))];
        M = zeros(3,3);
        for j=1:3
            M(j,:) = [1 Xt(j) Yt(j)];
        end
        w = M^-1 * eye(3,3);
        w = w';
        B = [w(1,2) 0 w(2,2) 0 w(3,2) 0; 0 w(1,3) 0 w(2,3) 0 w(3,3); w(1,3) w(1,2) w(2,3) w(2,2) w(3,3) w(3,2)];
        Bo = zeros(6,Nss);
        for j=1:3
            Bo(2*j-1,2*nr(j)-1) = 1;
            Bo(2*j,2*nr(j)) = 1;
        end
        wek = [Xt(1) - Xt(2), Yt(1) - Yt(2); Xt(2) - Xt(3), Yt(2) - Yt(3); Xt(3)-Xt(1), Yt(3)-Yt(1)];
        dl = [sqrt(wek(1,1)^2 + wek(1,2)^2), sqrt(wek(2,1)^2 + wek(2,2)^2), sqrt(wek(3,1)^2 + wek(3,2)^2)];
        p = 0.5*(dl(1) + dl(2) + dl(3));
        J = sqrt(p*(p-dl(1))*(p-dl(2))*(p-dl(3)));
        Kel = J * B' * D * B;
        K = K + Bo' * Kel * Bo;
    end
    for i=1:n
        if WB(i)>0
            G = zeros(n,2);
            ile = 0;
            for j=1:n
                if (i~=j & (WB(j)==3 | WS(j)==1))
                    ile = ile + 1;
                    d = sqrt((X(j) - X(i))^2 + (Y(j) - Y(i))^2);
                    G(ile,:) = [j d];
                end
            end          
            G(ile+1:end,:) = [];
            G = sortrows(G,2);
            hs = (G(1,2) + G(2,2))/2;
            if obsz<=2
                if abs(X(i)-a)<eps | abs(X(i)<eps);
                    hs = G(1,2)/2;
                end
            end
        end
        if WB(i)==3              
            F(2*i) = -q*hs;            
        end
        if WB(i)>0 & abs(X(i) - a)<eps & abs(Y(i)-b)<eps
            F(2*i) = F(2*i) - P;
        end
        if WB(i)==2
            F(2*i) = 0;
            F(2*i-1) = 0;
            K(2*i,:) = 0;
            K(:,2*i) = 0;
            K(2*i-1,:) = 0;
            K(:,2*i-1) = 0;
            K(2*i,2*i) = 1;
            K(2*i-1,2*i-1) = 1;
        end        
    end
    U2 = K^-1 * F;
    Ux2 = U2(1:2:Nss);
    Uy2 = U2(2:2:Nss);
    
    A = zeros(Nss,Nss);
    B = zeros(Nss,1);
    Ng = 4;
    for i=1:N
        nr = [T(i,1) T(i,2) T(i,3)];
        Xt = [X(nr(1)) X(nr(2)) X(nr(3))];
        Yt = [Y(nr(1)) Y(nr(2)) Y(nr(3))];
        M = zeros(3,3);
        for j=1:3
            M(j,:) = [1 Xt(j) Yt(j)];
        end
        w = M^-1 * eye(3,3);
        vx = w(2,:);
        vy = w(3,:);
        wek = [Xt(1) - Xt(2), Yt(1) - Yt(2); Xt(2) - Xt(3), Yt(2) - Yt(3); Xt(3)-Xt(1), Yt(3)-Yt(1)];
        dl = [sqrt(wek(1,1)^2 + wek(1,2)^2), sqrt(wek(2,1)^2 + wek(2,2)^2), sqrt(wek(3,1)^2 + wek(3,2)^2)];
        p = 0.5*(dl(1) + dl(2) + dl(3));
        J = sqrt(p*(p-dl(1))*(p-dl(2))*(p-dl(3)));
        [PGx,PGy,WAG] = cal_troj(Ng,X,Y,nr);
        w1 = E/(1-v^2);
        w2 = (1-v)/2;
        for k=1:Ng
            pg = [PGx(k) PGy(k)];
            G = gwiazda(X,Y,pg(1),pg(2));
            M = mwls(X,Y,pg(1),pg(2),G,il,1);
            for p = 1:3
                for j=1:il
                    nr1 = G(j,1);
                    A(2*nr(p)-1,2*nr1-1) = A(2*nr(p)-1,2*nr1-1) + J*WAG(k)*(w1*M(2,j)*vx(p) + w1*w2*M(3,j)*vy(p));
                    A(2*nr(p)-1,2*nr1) = A(2*nr(p)-1,2*nr1) + J*WAG(k)*(w1*v*M(3,j)*vx(p) + w1*w2*M(2,j)*vy(p));
                    A(2*nr(p),2*nr1-1) = A(2*nr(p),2*nr1-1) + J*WAG(k)*(w1*w2*M(3,j)*vx(p) + w1*v*M(2,j)*vy(p));
                    A(2*nr(p),2*nr1) = A(2*nr(p),2*nr1) + J*WAG(k)*(w1*w2*M(2,j)*vx(p) + w1*M(3,j)*vy(p));
                end
            end
        end
    end
    for i=1:n
        if WB(i)>0
            G = zeros(n,2);
            ile = 0;
            for j=1:n
                if (i~=j & (WB(j)==3 | WS(j)==1))
                    ile = ile + 1;
                    d = sqrt((X(j) - X(i))^2 + (Y(j) - Y(i))^2);
                    G(ile,:) = [j d];
                end
            end          
            G(ile+1:end,:) = [];
            G = sortrows(G,2);
            hs = (G(1,2) + G(2,2))/2;
            if obsz<=2
                if abs(X(i)-a)<eps | abs(X(i)<eps);
                    hs = G(1,2)/2;
                end
            end
        end
        if WB(i)==3           
            B(2*i) = -q*hs;
        end
        if WB(i)==2
            B(2*i) = 0;
            B(2*i-1) = 0;
            A(2*i,:) = 0;
            A(:,2*i) = 0;
            A(2*i-1,:) = 0;
            A(:,2*i-1) = 0;
            A(2*i,2*i) = 1;
            A(2*i-1,2*i-1) = 1;
        end
        if WB(i)>0 & abs(X(i) - a)<eps & abs(Y(i)-b)<eps
            B(2*i) = B(2*i)-P;
        end
    end
    U3 = A^-1 * B;
    Ux3 = U3(1:2:Nss);
    Uy3 = U3(2:2:Nss);
    
    
    iladap = 0;
    ADAP = zeros(3*iltr,4);
    eps = 10^-12;
    Rsr(siat,:) = 0;
    Rmax(siat,:) = 0;
    NSS(siat) = Nss;
    for i=1:N
        nr = [T(i,1) T(i,2) T(i,3)];
        Xt = [X(nr(1)) X(nr(2)) X(nr(3)) X(nr(1))];
        Yt = [Y(nr(1)) Y(nr(2)) Y(nr(3)) Y(nr(1))];
        for j=1:3
            xt = 0.5*(Xt(j) + Xt(j+1));
            yt = 0.5*(Yt(j) + Yt(j+1));
            ile = 0;
            for k=1:iladap
                dx = abs(ADAP(k,1) - xt);
                dy = abs(ADAP(k,2) - yt);
                if (dx<eps & dy<eps)
                    ile = ile + 1;
                end
            end
            if ile==0
                iladap = iladap + 1;
                ADAP(iladap,1) = xt;
                ADAP(iladap,2) = yt;
                if obsz<=2
                    if abs(xt)<eps | abs(yt)<eps | abs(xt-a)<eps | abs(yt-b)<eps
                        ADAP(iladap,4) = 1; % brzeg swobodny
                    end
                    if abs(yt-b)<eps & obc==1
                        ADAP(iladap,4) = 3; % obci¹¿enie
                    end
                    if abs(xt)<eps
                        ADAP(iladap,4) = 2; % zamocowanie
                    end
                end
                if obsz==3
                    if yt<=-3
                        ADAP(iladap,4) = 2;
                    end
                    if yt>=3
                        ADAP(iladap,4) = 3;
                    end
                end
                G = gwiazda(X,Y,xt,yt);
                ilr = 25;
                if n<ilr
                    ilr = n;
                end
                M = mwls(X,Y,xt,yt,G,ilr,2);
                Bx = zeros(ilr,1);
                By = zeros(ilr,1);
                for p=1:3
                    if p==1
                        Ux = Ux1;
                        Uy = Uy1;
                    end
                    if p==2
                        Ux = Ux2;
                        Uy = Uy2;
                    end
                    if p==3
                        Ux = Ux3;
                        Uy = Uy3;
                    end
                    for k=1:ilr
                        nr1 = G(k,1);
                        Bx(k) = Ux(nr1);
                        By(k) = Uy(nr1);
                    end
                    wx = M*Bx;
                    wy = M*By;
                    r1 = (wx(4) + v*wy(5)) + 0.5*(1-v)*(wx(6) + wy(5));
                    r2 = 0.5*(1-v)*(wx(5) + wy(4)) + (v*wx(5) + wy(6));
                    r = sqrt(r1^2 + r2^2);
                    if ADAP(iladap,4)==0
                        Rsr(siat,p) = Rsr(siat,p) + r^2;
                        if r>Rmax(siat,p)
                            Rmax(siat,p) = r;
                        end
                    end
                    if p==1
                        ADAP(iladap,3) = r;
                    end
                end
            end
        end
    end
    Rsr(siat,:) = (Rsr(siat,:)./iladap).^0.5;
    
    ADAP(iladap+1:end,:) = [];
    
    Rsr(siat,:) = 0;
    Rsrx(siat,:) = 0;
    Rsry(siat,:) = 0;
    Rmax(siat,:) = 0;
    for k=1:iltr
        nr = [T(k,1) T(k,2) T(k,3)];
        Xt = [X(nr(1)) X(nr(2)) X(nr(3))];
        Yt = [Y(nr(1)) Y(nr(2)) Y(nr(3))];
        wek = [Xt(1) - Xt(2), Yt(1) - Yt(2); Xt(2) - Xt(3), Yt(2) - Yt(3); Xt(3)-Xt(1), Yt(3)-Yt(1)];
        dl = [sqrt(wek(1,1)^2 + wek(1,2)^2), sqrt(wek(2,1)^2 + wek(2,2)^2), sqrt(wek(3,1)^2 + wek(3,2)^2)];
        p = 0.5*(dl(1) + dl(2) + dl(3));
        J = sqrt(p*(p-dl(1))*(p-dl(2))*(p-dl(3)));
        xs = (1/3)*(sum(Xt));
        ys = (1/3)*(sum(Yt));
        G = gwiazda(X,Y,xs,ys);
        ilr = 16;
        if n<ilr
            ilr = n;
        end
        % ilr = n;
        M = mwls(X,Y,xs,ys,G,ilr,2);
        Bx = zeros(ilr,1);
        By = zeros(ilr,1);
        for p=1:3
            if p==1
                Ux = Ux1;
                Uy = Uy1;
            end
            if p==2
                Ux = Ux2;
                Uy = Uy2;
            end
            if p==3
                Ux = Ux3;
                Uy = Uy3;
            end
            for j=1:ilr
                nr1 = G(j,1);
                Bx(j) = Ux(nr1);
                By(j) = Uy(nr1);
            end
            wx = M*Bx;
            wy = M*By;
            r1 = (wx(4) + v*wy(5)) + 0.5*(1-v)*(wx(6) + wy(5));
            r2 = 0.5*(1-v)*(wx(5) + wy(4)) + (v*wx(5) + wy(6));
            Rsrx(siat,p) = Rsrx(siat,p) + r1^2;
            Rsry(siat,p) = Rsry(siat,p) + r2^2;
            r = sqrt(r1^2 + r2^2);
            Rsr(siat,p) = Rsr(siat,p) + r^2;
            if r>Rmax(siat,p)
                Rmax(siat,p) = r;
            end
        end
    end
    for i=1:3
        Rsr(siat,i) = sqrt(Rsr(siat,i)/iltr);
        Rsrx(siat,i) = sqrt(Rsrx(siat,i)/iltr);
        Rsry(siat,i) = sqrt(Rsry(siat,i)/iltr);
    end
    %Rsr(siat,:) = (Rsr(siat,:)./iltr).^0.5;
    %Rsrx(siat,:) = (Rsrx(siat,:)./iltr).^0.5;
    %Rsry(siat,:) = (Rsry(siat,:)./iltr).^0.5;
    %iladap
    %ADAP
    
    %     for i=1:n
    %         G = gwiazda(X,Y,X(i),Y(i));
    %         M = mwls(X,Y,X(i),Y(i),G,il,2);
    %         Bx = zeros(il,1);
    %         By = zeros(il,1);
    %         for k=1:il
    %             nr1 = G(k,1);
    %             Bx(k) = Ux1(nr1);
    %             By(k) = Uy1(nr1);
    %         end
    %         wx = M*Bx;
    %         wy = M*By;
    %         r1 = (wx(4) + v*wy(5)) + 0.5*(1-v)*(wx(6) + wy(5));
    %         r2 = 0.5*(1-v)*(wx(5) + wy(4)) + (v*wx(5) + wy(6));
    %         r = sqrt(r1^2 + r2^2);
    %         if WB(i)==0
    %             r
    %         end
    %     end
    
    if siat>1
        ux1maxp = ux1max;
        ux2maxp = ux2max;
        ux3maxp = ux3max;
        uy1maxp = uy1max;
        uy2maxp = uy2max;
        uy3maxp = uy3max;
    end
    ux1max = max(abs(Ux1));
    ux2max = max(abs(Ux2));
    ux3max = max(abs(Ux3));
    uy1max = max(abs(Uy1));
    uy2max = max(abs(Uy2));
    uy3max = max(abs(Uy3));
    if siat>1
        Ex(siat-1,1) = abs(ux1maxp - ux1max)/abs(ux1max);
        Ex(siat-1,2) = abs(ux2maxp - ux2max)/abs(ux2max);
        Ex(siat-1,3) = abs(ux3maxp - ux3max)/abs(ux3max);
        Ey(siat-1,1) = abs(uy1maxp - uy1max)/abs(uy1max);
        Ey(siat-1,2) = abs(uy2maxp - uy2max)/abs(uy2max);
        Ey(siat-1,3) = abs(uy3maxp - uy3max)/abs(uy3max);
    end
    
    f = figure;
    ilwyk = 8;
    subplot(1,1,1);
    hold on
    grid on
    kol = 'rbgk';
    sym = '..od';
    sk = 0.05;
    for i=1:n
        p = plot(X(i),Y(i),[kol(WB(i)+1) sym(WB(i)+1)]);
        set(p,'MarkerFaceColor',kol(WB(i)+1));
        if WB(i)>0
            plot([X(i) X(i)+WEK(i,1)*sk],[Y(i) Y(i)+WEK(i,2)*sk],'g-');
        end
    end
    if obsz<=2
        axis([0 1.1*a 0 1.1*b]);
    end
    axis equal;
    title(['WÊZ£Y = ' num2str(n) ', ELEMENTY = ' num2str(N) ', SS = ' num2str(Nss) ', PKT.ADAP. = ' num2str(iladap)]);
%     for i=1:n
%         if WB(i)>0
%             plot([X(i) X(i)+WEK(i,1)*sk],[Y(i) Y(i)+WEK(i,2)*sk],'g-');
%             n1 = POP(i);
%             n2 = NAS(i);
%             A = [X(n1) 1;X(n2) 1];
%             B = [Y(n1) Y(n2)]';
%             a1 = A^-1*B;
%             plot([X(n1) X(n2)],[a1(1)*X(n1)+a1(2) a1(1)*X(n2)+a1(2)],'k-');
%             a2(1) = -1/a1(1);
%             a2(2) = Y(i) - a2(1)*X(i);
%             plot([X(n1) X(n2)],[a2(1)*X(n1)+a2(2) a2(1)*X(n2)+a2(2)],'k-');
%             pause
%         end
%     end
    pause(0.1);
    x = zeros(6,iltr);
    y = zeros(6,iltr);
    c = zeros(6,iltr);
    op = oppop;
    ilop = ilwyk + 4;
    while op<ilop
        if oppop~=ilwyk+2
            op = menu('WYBIERZ WIELKOŒÆ DO RYSOWANIA','DEFORMACJA','ODKSZTA£CENIA (Exx)','ODKSZTA£CENIA (Eyy)','ODKSZTA£CENIA (Exy)',...
                'NAPRÊ¯ENIA (Sxx)','NAPRÊ¯ENIA (Syy)','NAPRÊ¯ENIA (Sxy)','WYTÊ¯ENIE H-M-H','PUNKTY ADAPTACJI','ADAPTACJA','WYKRES ZBIE¯NOŒCI','KONIEC');
            oppop = op;
        end
        if op==ilop
            koniec = 1;
        end
        if op==ilwyk+1
            for i=1:iladap
                p = plot(ADAP(i,1),ADAP(i,2),[kol(ADAP(i,4)+1) '+']);
                set(p,'MarkerFaceColor',kol(ADAP(i,4)+1));
            end
        end
        if op==ilwyk+3
            Rmax(:,1)
            Rsr(:,1)
            figure;
            hold on
            grid on
            title('WYKRES ZBIE¯NOŒCI ROZWI¥ZANIA');
            zn = 'o*d';
            for i=1:3
                plot(log10(NSS(2:siat,1)),log10(Ex(1:siat-1,i)),['b' zn(i) '-']);
                plot(log10(NSS(2:siat,1)),log10(Ey(1:siat-1,i)),['g' zn(i) '-']);
              %  plot(log10(NSS(:,1)),log10(Rsrx(:,i)),['b' zn(i) '-']);
               % plot(log10(NSS(:,1)),log10(Rsry(:,i)),['g' zn(i) '-']);
                %        plot(log10(NSS(:,1)),log10(Rmax(:,i)),['r' zn(i) '-']);
            end
        end
        if op==ilwyk+2
            siat = siat + 1;
            if siat>=siatmax
                oppop = 1;
            end
            if adap==1
                nh = nh + 1;
            end
            if adap==2
                prog = 0.8;
                ADAP = sortrows(ADAP,3);
                rmax = max(abs(ADAP(:,3)));
                rprog = prog * rmax;
                ilmax = 2;
                if ilmax>iladap
                    ilmax = iladap;
                end
                ile = 0;
                r = rmax;
                Xd = zeros(ilmax,1);
                Yd = zeros(ilmax,1);
                WBd = zeros(ilmax,1);
                nr = iladap;
                while (r>=rprog & ile<=ilmax)
                    if ADAP(nr,4) == 0
                        ile = ile + 1;
                        r = ADAP(nr,3);
                        Xd(ile) = ADAP(nr,1);
                        Yd(ile) = ADAP(nr,2);
                        WBd(ile) = ADAP(nr,4);
                    end
                    nr = nr - 1;
                end
                n = n + ile;
                Xd(ile+1:end) = [];
                Yd(ile+1:end) = [];
                WBd(ile+1:end) = [];
                X = [X;Xd];
                Y = [Y;Yd];
                WB = [WB;WBd];
                WS = [WS;zeros(ile,1)];
                ile1 = 0;
                for k=nr:-1:1
                    G = gwiazda(X,Y,ADAP(k,1),ADAP(k,2));
                    liczba = 16;
                    wspd = 0.4;
                    wspb = 0.3;
                    if obc==1
                        liczba = 9;
                        wspd = 0.45;
                        wspb = 0.3;
                    end
                    dmin = min(G(1:liczba,2));
                    dmax = max(G(1:liczba,2));
                    wsp = wspd;
                    if ADAP(k,4)>0
                        wsp = wspb;
                    end
                    if abs(dmin/dmax)>wsp
                        n = n + 1;
                        X(n,1) = ADAP(k,1);
                        Y(n,1) = ADAP(k,2);
                        WB(n,1) = ADAP(k,4);
                        WS(n,1) = 0;
                        ile1 = ile1 + 1;
                    end
                end
                disp(['MAKSYMALNE RESIDUUM = ' num2str(ADAP(iladap,3),6)]);
                disp(['STARA LICZBA WEZLOW = ' num2str(n-ile1-ile)]);
                disp(['WEZLY Z KRYTERIUM BLEDU = ' num2str(ile)]);
                disp(['WEZLY Z KRYTERIUM GLADKOSCI = ' num2str(ile1)]);
                disp(['NOWA LICZBA WEZLOW = ' num2str(n)]);
            end
            break;
        end
        if op<=ilwyk
            wyk = op;            
            xmin = 0;
            xmax = a;
            ymin = 0;
            ymax = b;
            for k=1:4
                subplot(2,2,k);
                plot(0,0,'k.');
                hold on
                grid on
                Ux = zeros(n,3);
                Uy = zeros(n,3);                
                if k==2
                    Ux(:,1) = Ux1;
                    Uy(:,1) = Uy1;
                end
                if k==3
                    Ux(:,1) = Ux2;
                    Uy(:,1) = Uy2;
                end
                if k==4
                    Ux(:,1) = Ux3;
                    Uy(:,1) = Uy3;
                end
                wsp = 0.2/max(sqrt(Ux(:,1).^2 + Uy(:,1).^2));
                if isinf(wsp)==1
                    wsp = 1;
                end
                if wyk>1
                    if (k==2 | k==4)
                        st = 1;
                        ilp = il;
                        if k==2
                            st = 1;
                        end
                        for i=1:n
                            G = gwiazda(X,Y,X(i),Y(i));
                            M = mwls(X,Y,X(i),Y(i),G,ilp,st);
                            Bx = zeros(ilp,1);
                            By = zeros(ilp,1);
                            for j=1:ilp
                                nr = G(j,1);
                                Bx(j) = Ux(nr,1);
                                By(j) = Uy(nr,1);
                            end
                            wx = M*Bx;
                            wy = M*By;
                            Ux(i,2:3) = wx(2:3);
                            Uy(i,2:3) = wy(2:3);
                        end                        
                    end
                    if k==3
                        ILE = zeros(n,1);
                        for i=1:N
                            nr = [T(i,1) T(i,2) T(i,3)];
                            Xt = zeros(1,3);
                            Yt = zeros(1,3);
                            M = zeros(3,3);
                            for j=1:3
                                Xt(j) = X(nr(j));
                                Yt(j) = Y(nr(j));
                                M(j,:) = [1 Xt(j) Yt(j)];
                            end
                            w = M^-1 * eye(3,3);
                            px = w*[Ux(nr(1)) Ux(nr(2)) Ux(nr(3))]';
                            py = w*[Uy(nr(1)) Uy(nr(2)) Uy(nr(3))]';             
                            for j=1:3
                                ILE(nr(j)) = ILE(nr(j))+1;
                                Ux(nr(j),2) = Ux(nr(j),2) + px(2);
                                Ux(nr(j),3) = Ux(nr(j),3) + px(3);
                                Uy(nr(j),2) = Uy(nr(j),2) + py(2);
                                Uy(nr(j),3) = Uy(nr(j),3) + py(3);
                            end
                        end
                        for i=1:n
                            if ILE(i)>0
                                Ux(i,2) = Ux(i,2)/ILE(i);
                                Ux(i,3) = Ux(i,3)/ILE(i);
                                Uy(i,2) = Uy(i,2)/ILE(i);
                                Uy(i,3) = Uy(i,3)/ILE(i);                                
                            end
                        end                        
                    end
                end
                if wyk==1 | k==1
                    for i=1:n
                        p = plot(X(i)+Ux(i,1)*wsp,Y(i)+Uy(i,1)*wsp,[kol(WB(i)+1) sym(WB(i)+1)]);
                        set(p,'MarkerFaceColor',kol(WB(i)+1));
                    end
                end
                for i=1:iltr
                    nr = [T(i,1) T(i,2) T(i,3)];
                    Xt = zeros(1,3);
                    Yt = zeros(1,3);
                    for j=1:3
                        Xt(j) = X(nr(j)) + Ux(nr(j),1)*wsp;
                        Yt(j) = Y(nr(j)) + Uy(nr(j),1)*wsp;
                        if Xt(j)<xmin
                            xmin = Xt(j);
                        end
                        if Xt(j)>xmax
                            xmax = Xt(j);
                        end
                        if Yt(j)<ymin
                            ymin = Yt(j);
                        end
                        if Yt(j)>ymax
                            ymax = Yt(j);
                        end
                    end
                    if wyk==1 | k==1
                        plot([Xt Xt(1)],[Yt Yt(1)],'b--');
                    end
                    if (wyk>1 & k>1)
                        for j=1:3
                            m = nr(j);
                            x(2*j-1,i) = X(m);%+Ux(m,1);
                            y(2*j-1,i) = Y(m);%+Uy(m,1);
                            if wyk==2
                                c(2*j-1,i) = Ux(m,2);
                            end
                            if wyk==3
                                c(2*j-1,i) = Uy(m,3);
                            end
                            if wyk==4
                                c(2*j-1,i) = Ux(m,3) + Uy(m,2);
                            end
                            if wyk==5
                                c(2*j-1,i) = (E/(1-v^2))*(Ux(m,2) + v*Uy(m,3));
                            end
                            if wyk==6
                                c(2*j-1,i) = (E/(1-v^2))*(v*Ux(m,2) + Uy(m,3));
                            end
                            if wyk==7
                                c(2*j-1,i) = (E/(1-v^2))*((1-v)/2)*(Ux(m,3) + Uy(m,2));
                            end
                            if wyk==8
                                s1 = (E/(1-v^2))*(Ux(m,2) + v*Uy(m,3));
                                s2 = (E/(1-v^2))*(v*Ux(m,2) + Uy(m,3));
                                s3 = (E/(1-v^2))*((1-v)/2)*(Ux(m,3) + Uy(m,2));
                                %c(2*j-1,i) = (1/sqrt(2)) * sqrt( (s1 - s2)^2 + s1^2 + s2^2 +  6 * s3^2 );
                                lam = eig([s1 s3;s3 s2]);
                                sg1 = lam(1);
                                sg2 = lam(2);
                                c(2*j-1,i) = sqrt(sg1^2 - sg1*sg2 + sg2^2);
                            end
                        end                        
                        for nod=2:2:6
                            n1 = nod-1;
                            n2 = nod+1;
                            if n2>6
                                n2=1;
                            end
                            x(nod,i) = (x(n1,i) + x(n2,i))/2;
                            y(nod,i) = (y(n1,i) + y(n2,i))/2;
                            c(nod,i) = (c(n1,i) + c(n2,i))/2;
                        end
                        t = fill(x(:,i),y(:,i),c(:,i));
                        set(t,'Edgecolor','interp');
                    end
                end
                if k==1
                    title(['WÊZ£Y (n=' num2str(n) ') DLA BMRS i ELEMENTY (N=' num2str(N) ') DLA MES']);
                end
                Nap = ['E_x_x^m^a^x';'E_y_y^m^a^x';'E_x_y^m^a^x';'S_x_x^m^a^x';'S_y_y^m^a^x';'S_x_y^m^a^x';'S_h_m^r^e^d'];
                if k==2
                    if wyk==1
                        title(['BMRS_L_O_K: U_x^m^a^x = ' num2str(ux1max,4) ', U_y^m^a^x = ' num2str(uy1max,4)]);
                    end
                    if wyk>1
                        title(['BMRS_L_O_K: ' Nap(wyk-1,:) ' = ' num2str(max(max(abs(c))),4)]);
                    end
                end
                if k==3
                    if wyk==1
                        title(['MES: U_x^m^a^x = ' num2str(ux2max,4) ', U_y^m^a^x = ' num2str(uy2max,4)]);
                    end
                    if wyk>1
                        title(['MES: ' Nap(wyk-1,:) ' = ' num2str(max(max(abs(c))),4)]);
                    end
                end
                if k==4
                    if wyk==1
                        title(['BMRS_W_A_R: U_x^m^a^x = ' num2str(ux3max,4) ', U_y^m^a^x = ' num2str(uy3max,4)]);
                    end
                    if wyk>1
                        title(['BMRS_W_A_R: ' Nap(wyk-1,:) ' = ' num2str(max(max(abs(c))),4)]);
                    end
                end
            end
            
            for i=1:4
                subplot(2,2,i);
                axis([xmin xmax ymin ymax]);
            end
            s = '';
            if wyk==1
                s = 'PRZEMIESZCZENIA (U_x i U_y)';
            end
            if wyk==2
                s = 'ODKSZTA£CENIA (E_x_x)';
            end
            if wyk==3
                s = 'ODKSZTA£CENIA (E_y_y)';
            end
            if wyk==4
                s = 'ODKSZTA£CENIA (E_x_y)';
            end
            if wyk==5
                s = 'NAPRÊ¯ENIA (S_x_x)';
            end
            if wyk==6
                s = 'NAPRÊ¯ENIA (S_y_y)';
            end
            if wyk==7
                s = 'NAPRÊ¯ENIA (S_x_y)';
            end
            if wyk==8
                s = 'NAPRÊ¯ENIA HMH (S_H_M_H)';
            end
            subplot(2,2,1);
            t = text(1*xmax,1.1*ymax,s);
            set(t,'FontSize',[12]);
            set(f,'Position',[1 41 1440 790]);
            
            for i=1:4
                subplot(2,2,i);
                hold off
            end
        end
    end
    close all
end

function [M,K] = mwls(X,Y,x,y,G,il,st)
licz = liczba(st);
m = licz(1);
m1 = licz(2);
A = zeros(il,m);
W = zeros(il,il);
A1 = zeros(il,m1);
for i=1:il
    nr = G(i,1);
    hx = X(nr) - x;
    hy = Y(nr) - y;
    d = sqrt(hx*hx + hy*hy);
    inter = [1 hx hy 0.5*hx^2 hx*hy 0.5*hy^2 (1/6)*hx^3 (1/2)*hx^2*hy (1/2)*hx*hy^2 (1/6)*hy^3 (1/24)*hx^4 (1/6)*hx^3*hy (1/4)*hx^2*hy^2 (1/6)*hx*hy^3 (1/24)*hy^4];
    A(i,:) = inter(1:m);
    W(i,i) = 1/(d^(2*(st+1)) + 10^-12);
    A1(i,:) = inter((m+1):(m+m1));
end
M = (A'*W*A)^-1 * A'*W;
K = M*A1;

function [licz] = liczba(st)
if st==1
    m = 3;
    m1 = 3;
end
if st==2
    m = 6;
    m1 = 9;
end
if st==3
    m = 10;
    m1 = 5;
end
if st==4
    m = 15;
    m1 = 0;
end
licz = [m m1];


function [G] = gwiazda(X,Y,x,y)
n = length(X);
G = zeros(n,2);
for i=1:n
    d = sqrt((X(i) - x)^2 + (Y(i) - y)^2);
    G(i,:) = [i d];
end
G = sortrows(G,2);

function [pg,WAG] = calka(Ng)
switch Ng
    case 1
        pg = [0];
        WAG = [2];
    case 2
        pg = [-1./sqrt(3.) 1./sqrt(3.)];
        WAG = [1 1];
    case 3
        pg = [-1/sqrt(0.6) 0 1/sqrt(0.6)];
        WAG = [ 5/9 8/9 5/9 ];
    case 4
        pg = [-0.861136311594053 -0.339981043584856 0.339981043584856 0.861136311594053];
        WAG = [0.347854845137454 0.652145154862546 0.652145154862546 0.347854845137454];
    case 5
        mz = [0.906179845939 0.538469310106];
        pg = [-mz 0 mz(2) mz(1)];
        wg = [0.236926885056 0.478628670499 0.568888888889];
        WAG = [wg wg(2) wg(1)];
    case 6
        mz = [0.932469514203 0.661209386466 0.238619186083];
        wg = [0.171324492379 0.360761573048 0.467913934573];
        pg = [-mz mz(3) mz(2) mz(1)];
        WAG = [wg wg(3) wg(2) wg(1)];
    case 7
        mz = [0.949107912343 0.741531185599 0.405845151377];
        wg = [0.129484966169 0.279705391489 0.381830050505 0.417959183673];
        pg = [-mz 0 mz(3) mz(2) mz(1)];
        WAG = [wg wg(3) wg(2) wg(1)];
    otherwise
        pg = zeros(Ng,1);
        WAG = zeros(Ng,1);
        for i=1:Ng
            pg(i) = cos(pi*(2*i-1)/(2*Ng));
            WAG(i) = pi/Ng;
        end
end

function [PGx,PGy,WAG1] =  cal_troj(Ng1,X,Y,w)
PGx = zeros(1,Ng1);
PGy = zeros(1,Ng1);
WAG1 = zeros(1,Ng1);
switch Ng1
    case 1
        WAG1 = [1];
    case 2
        WAG1 = [1/3 1/3 1/3];
    case 4
        WAG1 = [-27/48 25/48 25/48 25/48];
    case 7
        WAG1 = [0.225 0.1323941527 0.1323941527 0.1323941527 0.1259391805 0.1259391805 0.1259391805];
end
switch Ng1
    case 1
        PGx(1) = (1./3.) * (X(w(1)) + X(w(2)) + X(w(3)));
        PGy(1) = (1./3.) * (Y(w(1)) + Y(w(2)) + Y(w(3)));
    case 2
        PGx(1) = 0.5*(X(w(1)) + X(w(2)));
        PGy(1) = 0.5*(Y(w(1)) + Y(w(2)));
        PGx(2) = 0.5*(X(w(2)) + X(w(3)));
        PGy(2) = 0.5*(Y(w(2)) + Y(w(3)));
        PGx(3) = 0.5*(X(w(1)) + X(w(3)));
        PGy(3) = 0.5*(Y(w(1)) + Y(w(3)));
    case 4
        PGx(1) = (1./3.) * (X(w(1)) + X(w(2)) + X(w(3)));
        PGy(1) = (1./3.) * (Y(w(1)) + Y(w(2)) + Y(w(3)));
        PGx(2) = 0.6*X(w(1)) + 0.2*X(w(2)) + 0.2*X(w(3));
        PGy(2) = 0.6*Y(w(1)) + 0.2*Y(w(2)) + 0.2*Y(w(3));
        PGx(3) = 0.2*X(w(1)) + 0.6*X(w(2)) + 0.2*X(w(3));
        PGy(3) = 0.2*Y(w(1)) + 0.6*Y(w(2)) + 0.2*Y(w(3));
        PGx(4) = 0.2*X(w(1)) + 0.2*X(w(2)) + 0.6*X(w(3));
        PGy(4) = 0.2*Y(w(1)) + 0.2*Y(w(2)) + 0.6*Y(w(3));
    case 7
        a1 = 0.0597158717;
        b1 = 0.4701420641;
        a2 = 0.7974269853;
        b2 = 0.1012865073;
        PGx(1) = (1/3)*(X(w(1)) + X(w(2)) + X(w(3)));
        PGy(1) = (1/3)*(Y(w(1)) + Y(w(2)) + Y(w(3)));
        PGx(2) = a1 * X(w(1)) + b1 * X(w(2)) + b1 * X(w(3));
        PGy(2) = a1 * Y(w(1)) + b1 * Y(w(2)) + b1 * Y(w(3));
        PGx(3) = b1 * X(w(1)) + a1 * X(w(2)) + b1 * X(w(3));
        PGy(3) = b1 * Y(w(1)) + a1 * Y(w(2)) + b1 * Y(w(3));
        PGx(4) = b1 * X(w(1)) + b1 * X(w(2)) + a1 * X(w(3));
        PGy(4) = b1 * Y(w(1)) + b1 * Y(w(2)) + a1 * Y(w(3));
        PGx(5) = a2 * X(w(1)) + b2 * X(w(2)) + b2 * X(w(3));
        PGy(5) = a2 * Y(w(1)) + b2 * Y(w(2)) + b2 * Y(w(3));
        PGx(6) = b2 * X(w(1)) + a2 * X(w(2)) + b2 * X(w(3));
        PGy(6) = b2 * Y(w(1)) + a2 * Y(w(2)) + b2 * Y(w(3));
        PGx(7) = b2 * X(w(1)) + b2 * X(w(2)) + a2 * X(w(3));
        PGy(7) = b2 * Y(w(1)) + b2 * Y(w(2)) + a2 * Y(w(3));
end

function [ilwew] = prostokat(P,Xw,Yw,G)
%P = [0 0; 1 2; 0 3; -1 1];
%P = [0 0; 2 0; 2 3; 0 3];
%plot([P(:,1)' P(1,1)],[P(:,2)' P(1,2)],'ro-');
%hold on
%grid on

P1(1,:) = [min(P(:,1)) min(P(:,2))];
P1(2,:) = [max(P(:,1)) min(P(:,2))];
P1(3,:) = [max(P(:,1)) max(P(:,2))];
P1(4,:) = [min(P(:,1)) max(P(:,2))];
P = P1;

ilwew = 0;
Xp = [P(1,1) P(2,1) P(3,1) P(4,1) P(1,1)];
Yp = [P(1,2) P(2,2) P(3,2) P(4,2) P(1,2)];

eps = 10^-12;
n = length(G);
for i=1:n+1
    if i==1
        X(i) = 0.5*(Xp(1)+Xp(3));
        Y(i) = 0.5*(Yp(1)+Yp(3));
    end
    if i>1
        nr = G(i-1,1);
        X(i) = Xw(nr);
        Y(i) = Yw(nr);
    end    
    sum  = 0;
    for j=1:4
        n1 = j;
        n2 = j+1;
        w1 = [Xp(n1)-X(i) Yp(n1)-Y(i)];
        w2 = [Xp(n2)-X(i) Yp(n2)-Y(i)];
        w3 = [Xp(n1)-Xp(n2) Yp(n1)-Yp(n2)];
        d1 = norm(w1);
        d2 = norm(w2);
        d3 = norm(w3);
        p = 0.5*(d1 + d2 + d3);
        J = sqrt(p*(p-d1)*(p-d2)*(p-d3));
        sum = sum + J;
    end
    if i==1
        Pole = sum;
    else
        if abs(Pole - sum)<eps
            ilwew = ilwew +1;
        end       
    end
end