% Esempio di uso del programma di reticolazione EasyMesh
% per la soluzione di un problema elettromagnetico 2D
% (supponendo che 'easymesh.exe' e 'dos4gw.exe' si trovino
% nel direttorio corrente)
%
% Versione 12/7/2005 - M. Politi

clear all
format compact

%
% Prepara file dati per reticolazione
% (Questo è solo l'esempio generato automaticamente da EasyMesh)
%

!easymesh +example

%
% Lancia programma di reticolazione
% (EasyMesh da 'EXAMPLE.D' genera 'EXAMPLE.N', 'EXAMPLE.S' e 'EXAMPLE.E')
% 

fname = 'EXAMPLE';
command = ['!easymesh ' fname];
eval(command)

%
% Leggi dati reticolazione
% ('Get_em_data' provvede anche a convertire la numerazione di 
% nodi, lati e facce da 0-(N-1) a 1-N; provvede inoltre a correggere 
% l'orientamento dei lati al contorno se non congruenti) 
%
%   'Node': matrice Nn x 4 (da file EXAMPLE.N)
%   'Side': matrice Ns x 6 (da file EXAMPLE.S)
%   'Elem': matrice Nn x 13 (da file EXAMPLE.E)
% 

[Node, Nn, Side, Ns, Elem, Ne] = Get_em_data(fname);

%
% Visualizza reticolo Delaunay & Voronoi
%

trimesh(Elem(:,2:4), Node(:,2), Node(:,3), 'Color',[0 0 1])
hold on
    X = Node(:,2);
    Y = Node(:,3);
    Tri = Elem(:,2:4);
    [Xv, Yv] = voronoi(X, Y, Tri);
    plot(Xv,Yv,'-r')
hold off
axis equal

%
% Calcola parametri geometrici
%
%   'L': vettore lunghezze lati (Ns x 1)
%   'Ld': vettore lunghezze lati duali (Ns x 1)
%   'S': vettore aree elementi (Ne x 1)
%   'Sd': vettore aree elementi duali (Nn x 1)
%

[L, Ld, S, Sd] = Calc_em_par(Node, Nn, Side, Ns, Elem, Ne);

%
% Elabora ...
% (seguono esempi di operazioni elementari)
%

%           matrice di incidenza lati-nodi (matrice gradiente)

n1 = Side(:,2);         % indici 1.o nodo (1 per lato)
n2 = Side(:,3);         % indici 2.o nodo (idem)
GG = sparse(Side(:,1), n1, -1, Ns, Nn) + sparse(Side(:,1), n2, 1, Ns, Nn);

%           matrice lunghezze lati duali

LLd = spdiags(Ld, 0, Ns, Ns);

%           nodi al contorno con marker = 3

ind = find(Node(:,4) == 3);

%           ...
