%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
% --------------------------
% Discoef(1,k) = nju in region k
% Discoef(2,k) = E   in region k
%
function ...
[Node,Edge,Face,Node_flagx,Node_flagy,Edge_flagx,Edge_flagy,...
           Face_flag,Face_thick,Face_Node,...
           Disco,Discoef,wh] = Square_disco_hier(jump)
wh='jp';
Disco(1,1,1) = 0.0;     %form x
Disco(2,1,1) = 0.5;     %to   x
Disco(1,2,1) = 0.0;     %from y
Disco(2,2,1) = 1.0;     %to   y
%
Disco(1,1,2) = 0.5;     %form x
Disco(2,1,2) = 1.0;     %to   x
Disco(1,2,2) = 0.0;     %from y
Disco(2,2,2) = 0.5;     %to   y
%
Disco(1,1,3) = 0.5;     %form x
Disco(2,1,3) = 0.75;    %to   x
Disco(1,2,3) = 0.5;     %from y
Disco(2,2,3) = 0.75;    %to   y
%
Disco(1,1,4) = 0.5;     %form x
Disco(2,1,4) = 0.75;    %to   x
Disco(1,2,4) = 0.5;     %from y
Disco(2,2,4) = 1.0;     %to   y

%
Disco(1,1,5) = 0.75;    %form x
Disco(2,1,5) = 1.0;     %to   x
Disco(1,2,5) = 0.5;     %from y
Disco(2,2,5) = 1.0;     %to   y

Discoef(1,1:5) = 1;     % coefficient in the bilinear form
Discoef(1,3)   = jump;  % coefficient in the bilinear form
Discoef(2,1:5) = 1;     % not used in the scalar case

% ---------------- nodes -----
Node(1, 1) = 0;
Node(2, 1) = 0;
Node(1, 2) = 0.5;
Node(2, 2) = 0;
Node(1, 3) = 1;
Node(2, 3) = 0;
Node(1, 4) = 0.75;
Node(2, 4) = 0.25;
Node(1, 5) = 0;
Node(2, 5) = 0.5;
Node(1, 6) = 0.5;
Node(2, 6) = 0.5;
Node(1, 7) = 0.75;
Node(2, 7) = 0.5;
Node(1, 8) = 1;
Node(2, 8) = 0.5;
Node(1, 9) = 0.25;
Node(2, 9) = 0.75;
Node(1,10) = 0.5;
Node(2,10) = 0.75;
Node(1,11) = 0.75;
Node(2,11) = 0.75;
Node(1,12) = 0;
Node(2,12) = 1;
Node(1,13) = 0.5;
Node(2,13) = 1;
Node(1,14) = 1;
Node(2,14) = 1;
% ---------------- edges -----
Edge(1, 1) = 1;
Edge(2, 1) = 2;
Edge(1, 2) = 2;
Edge(2, 2) = 3;
Edge(1, 3) = 1;
Edge(2, 3) = 5;
Edge(1, 4) = 2;
Edge(2, 4) = 5;
Edge(1, 5) = 2;
Edge(2, 5) = 6;
Edge(1, 6) = 2;
Edge(2, 6) = 4;
Edge(1, 7) = 4;
Edge(2, 7) = 6;
Edge(1, 8) = 4;
Edge(2, 8) = 7;
Edge(1, 9) = 4;
Edge(2, 9) = 8;
Edge(1,10) = 3;
Edge(2,10) = 4;
Edge(1,11) = 3;
Edge(2,11) = 8;
Edge(1,12) = 5;
Edge(2,12) = 6;
Edge(1,13) = 6;
Edge(2,13) = 7;
Edge(1,14) = 7;
Edge(2,14) = 8;
Edge(1,15) = 5;
Edge(2,15) = 12;
Edge(1,16) = 5;
Edge(2,16) = 9;
Edge(1,17) = 6;
Edge(2,17) = 9;
Edge(1,18) = 6;
Edge(2,18) = 10;
Edge(1,19) = 7;
Edge(2,19) = 10;
Edge(1,20) = 7;
Edge(2,20) = 11;
Edge(1,21) = 8;
Edge(2,21) = 11;
Edge(1,22) = 8;
Edge(2,22) = 14;
Edge(1,23) = 9;
Edge(2,23) = 10;
Edge(1,24) = 10;
Edge(2,24) = 11;
Edge(1,25) = 9;
Edge(2,25) = 12;
Edge(1,26) = 9;
Edge(2,26) = 13;
Edge(1,27) = 10;
Edge(2,27) = 13;
Edge(1,28) = 11;
Edge(2,28) = 13;
Edge(1,29) = 11;
Edge(2,29) = 14;
Edge(1,30) = 12;
Edge(2,30) = 13;
Edge(1,31) = 13;
Edge(2,31) = 14;

% ---------------- Face -----
Face(1, 1) = 1;
Face(2, 1) = 3;
Face(3, 1) = 4;
Face(1, 2) = 4;
Face(2, 2) = 5;
Face(3, 2) = 12;
Face(1, 3) = 5;
Face(2, 3) = 6;
Face(3, 3) = 7;
Face(1, 4) = 2;
Face(2, 4) = 6;
Face(3, 4) = 10;
Face(1, 5) = 9;
Face(2, 5) = 10;
Face(3, 5) = 11;
Face(1, 6) = 7;
Face(2, 6) = 8;
Face(3, 6) = 13;
Face(1, 7) = 8;
Face(2, 7) = 9;
Face(3, 7) = 14;
Face(1, 8) = 12;
Face(2, 8) = 17;
Face(3, 8) = 16;
Face(1, 9) = 17;
Face(2, 9) = 18;
Face(3, 9) = 23;
Face(1,10) = 13;
Face(2,10) = 19;
Face(3,10) = 18;
Face(1,11) = 19;
Face(2,11) = 20;
Face(3,11) = 24;
Face(1,12) = 14;
Face(2,12) = 21;
Face(3,12) = 20;
Face(1,13) = 21;
Face(2,13) = 22;
Face(3,13) = 29;
Face(1,14) = 15;
Face(2,14) = 16;
Face(3,14) = 25;
Face(1,15) = 23;
Face(2,15) = 27;
Face(3,15) = 26;
Face(1,16) = 24;
Face(2,16) = 28;
Face(3,16) = 27;
Face(1,17) = 25;
Face(2,17) = 26;
Face(3,17) = 30;
Face(1,18) = 28;
Face(2,18) = 29;
Face(3,18) = 31;

Node_flagx(1:14,1) = 0;
Node_flagx(1, 1) = Flags('Dirichlet');
Node_flagx(2, 1) = Flags('Dirichlet');
Node_flagx(3, 1) = Flags('Dirichlet');
Node_flagx(5, 1) = Flags('Dirichlet');
Node_flagx(8, 1) = Flags('Dirichlet');
Node_flagx(12,1) = Flags('Dirichlet');
Node_flagx(14,1) = Flags('Dirichlet');
Node_flagx(13,1) = Flags('Neumann');

Edge_flagx(1:31,1) = 0;
Edge_flagx(1, 1) = Flags('Dirichlet');
Edge_flagx(2, 1) = Flags('Dirichlet');
Edge_flagx(3, 1) = Flags('Dirichlet');
Edge_flagx(11,1) = Flags('Dirichlet');
Edge_flagx(15,1) = Flags('Dirichlet');
Edge_flagx(22,1) = Flags('Dirichlet');
Edge_flagx(30,1) = Flags('Neumann');
Edge_flagx(31,1) = Flags('Neumann');

Node_flagy(1:14,1) = 0;
Node_flagy(1, 1) = Flags('Dirichlet');
Node_flagy(2, 1) = Flags('Dirichlet');
Node_flagy(3, 1) = Flags('Dirichlet');
Node_flagy(5, 1) = Flags('Neumann');
Node_flagy(8, 1) = Flags('Neumann');
Node_flagy(12,1) = Flags('Neumann');
Node_flagy(13,1) = Flags('Neumann');
Node_flagy(14,1) = Flags('Neumann');

Edge_flagy(1:31,1) = 0;
Edge_flagy(1, 1) = Flags('Dirichlet');
Edge_flagy(2, 1) = Flags('Dirichlet');
Edge_flagy(3, 1) = Flags('Neumann');
Edge_flagy(11,1) = Flags('Neumann');
Edge_flagy(15,1) = Flags('Neumann');
Edge_flagy(22,1) = Flags('Neumann');
Edge_flagy(30,1) = Flags('Neumann');
Edge_flagy(31,1) = Flags('Neumann');

Face_flag(1:18,1)  = 1;
Face_flag(10,1)    = 3;
Face_flag(11,1)    = 3;
Face_thick(1:18,1) = 1;

Face_Node(1:3, 1,1) = [1;2;5];
Face_Node(1:3, 2,1) = [2;5;6];
Face_Node(1:3, 3,1) = [2;4;6];
Face_Node(1:3, 4,1) = [2;3;4];
Face_Node(1:3, 5,1) = [4;3;8];
Face_Node(1:3, 6,1) = [4;6;7];
Face_Node(1:3, 7,1) = [4;7;8];
Face_Node(1:3, 8,1) = [5;6;9];
Face_Node(1:3, 9,1) = [6;9;10];
Face_Node(1:3,10,1) = [6;7;10];
Face_Node(1:3,11,1) = [7;10;11];
Face_Node(1:3,12,1) = [7;8;11];
Face_Node(1:3,13,1) = [8;11;14];
Face_Node(1:3,14,1) = [5;9;12];
Face_Node(1:3,15,1) = [9;10;13];
Face_Node(1:3,16,1) = [10;11;13];
Face_Node(1:3,17,1) = [9;13;12];
Face_Node(1:3,18,1) = [11;14;13];

return