%Initial triangulation
% the array 'Disco' describes the discontinuity regions
% Disco(*,*,k) is the kth region
% Disco(*,1,k) x-interval
% Disco(*,2,k) y-interval
% --------------------------
%
function ...
[Node,Edge,Face,Node_flagx,Node_flagy,Edge_flagx,Edge_flagy,...
           Face_flag,Face_thick,Face_Node,...
           Disco,Discoef,wh] = Square_ConvDiff
wh='cd';
Disco(1,1,1) = 0;    %form x
Disco(2,1,1) = 1;    %to   x
Disco(1,2,1) = 0;    %from y
Disco(2,2,1) = 1;    %to   y
Discoef(1,1) = 1;  % 
Discoef(2,1) = 1;  % 

% ---------------- nodes -----
Node(1,1) = 0;
Node(2,1) = 0;
Node(1,2) = 1;
Node(2,2) = 0;
Node(1,3) = 0;
Node(2,3) = 1;
Node(1,4) = 1;
Node(2,4) = 1;
% ---------------- edges -----
Edge(1, 1) = 1;
Edge(2, 1) = 2;
Edge(1, 2) = 3;
Edge(2, 2) = 1;
Edge(1, 3) = 2;
Edge(2, 3) = 3;
Edge(1, 4) = 4;
Edge(2, 4) = 2;
Edge(1, 5) = 3;
Edge(2, 5) = 4;
% ---------------- Face -----
Face(1,1) = 1;
Face(2,1) = 3;
Face(3,1) = 2;
Face(1,2) = 3;
Face(2,2) = 5;
Face(3,2) = 4;

Node_flagx(1,1) = Flags('Dirichlet');
Node_flagx(2,1) = Flags('Dirichlet');
Node_flagx(3,1) = Flags('Dirichlet');
Node_flagx(4,1) = Flags('Neumann');
Node_flagy = Node_flagx;  %dummy

Edge_flagx(1,1) = Flags('Dirichlet');
Edge_flagx(2,1) = Flags('Dirichlet');
Edge_flagx(4,1) = Flags('Neumann');
Edge_flagx(5,1) = Flags('Neumann');
Edge_flagx(3,1) = 0;
Edge_flagy = Edge_flagx;  %dummy

Face_flag(1:2,1) = 1;
Face_thick(1:2,1) = 1;

Face_Node(1:3,1,1)=[1;2;3];
Face_Node(1:3,2,1)=[2;3;4];


return