%%%%%%%%%%%%%%%%%%%% Easy DD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This routine refines the very coarse mesh only once
% to mark properly the nodes and edges
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Node,Edge,Face,...
       Node_flagx,Node_flagy,Edge_flagx,Edge_flagy,...
       Face_flag,Face_thick,Face_Node,Face_Parent,...
       nface_lvl,nnode_lvl] = ...
   my_Refine_hier(Node,Edge,Face,...
       Node_flagx,Node_flagy,Edge_flagx,Edge_flagy,...
       Face_flag,Face_thick,Face_Node,Face_Parent,...
       nface_lvl,nnode_lvl,lvl)

nnode = size(Node,2);
nedge = size(Edge,2);
nface = size(Face,2);

new_node = 0;
new_edge = 0;

% Make Node, Edge and Face immediately as large as needed
% Every face will be divided in 4
Face(1,4*nface)      = 0;
Face_flag(4*nface,1) = 0;
Face_thick(4*nface,1) = 1;
% every edge will produce a new node
Node(1,nnode+nedge)  = 0;
Node_flagx(4*nface,1) = 0;
Node_flagy(4*nface,1) = 0;
% Euler #F + #N - #E = 2  but outsideregion is also a Face
Edge(1,3*nface+2*nedge)       = 0;
Edge_flagx(3*nface+2*nedge,1) = 0;
Edge_flagy(3*nface+2*nedge,1) = 0;

% --------------------- halve the edges
for iedge = 1:nedge,
    inode1 = Edge(1,iedge); inode2 = Edge(2,iedge);
%disp(['Edge ',int2str(iedge)])
%disp(['End_1 ' int2str(inode1)])
%disp(['End_2 ' int2str(inode2)]) ']'])
%disp(['Midpoint coord: [' num2str(xM) ',' num2str(yM) ']'])
    new_node = new_node + 1;
    Node(:,nnode+new_node) = (Node(:,inode1)+Node(:,inode2))*0.5;
    new_edge = new_edge + 1;
    Edge(1,nedge+new_edge) = nnode+new_node; 
    Edge(2,nedge+new_edge) = Edge(2,iedge);
    Edge(2,iedge) = nnode+new_node;
    Edge_flagx(nedge+new_edge,1) = Edge_flagx(iedge,1);
    Edge_flagy(nedge+new_edge,1) = Edge_flagy(iedge,1);
    Node_flagx(nnode+new_node,1) = Edge_flagx(iedge,1);
    Node_flagy(nnode+new_node,1) = Edge_flagy(iedge,1);
end
% --------------------- add internal edges and new faces

new_face  = 0;
iedge     = zeros(1,6);
node_mid  = zeros(1,3);
Face2Edge = zeros(1,3);
CurEdge   = nedge*2+1;   % CurEdge=size(Edge,2)+1;

for iface = 1:nface,  % ------------- begin loop over faces 
    for j = 1:3,
        iedge(j) = Face(j,iface);
        iedge(j+3) = iedge(j) + nedge; % similar as Edge_temp
        node_mid(j)= Edge(2,iedge(j)); % these are the new points
        node_old(j)= Edge(1,iedge(j)); % these are the old points
    end

% ----------------------- add three more internal edges

    Edge(1,CurEdge)   = node_mid(3);
    Edge(2,CurEdge)   = node_mid(2);
    Edge(1,CurEdge+1) = node_mid(3);
    Edge(2,CurEdge+1) = node_mid(1);
    Edge(1,CurEdge+2) = node_mid(1);
    Edge(2,CurEdge+2) = node_mid(2);    

%disp(['newEdge_1:' int2str(CurEdge) ...
%    ' [' int2str(node_mid(3)) ',' int2str(node_mid(2)) ']'])
%disp(['newEdge_2:' int2str(CurEdge+1) ...
%    ' [' int2str(node_mid(3)) ',' int2str(node_mid(1)) ']'])
%disp(['newEdge_3:' int2str(CurEdge+2) ...
%    ' [' int2str(node_mid(1)) ',' int2str(node_mid(2)) ']']),wait

% ----------- put the middle triangle at the place of the old one
      for j = 1:3,
        Face(j, iface) = CurEdge + j-1; 
      end
      Face_Node(1:3,iface) =  node_mid';      
% ----------------------- update faces
      for j = 1:3,
        Face2Edge(1) = Face(j, iface);
        node_mid(1)  = Edge(1, Face2Edge(1));
        node_mid(2)  = Edge(2, Face2Edge(1));
        Exit = 0;
        for k = 1:6,
          for s = 1:2,
            if(Edge(s, iedge(k)) == node_mid(1)),
              for q = 1:6,
                for w = 1:2,
                  if( (Edge(w, iedge(q))    == node_mid(2))...
                     &(Edge(3 - w, iedge(q))== Edge(3 - s, iedge(k)))),
                    Face2Edge(2) = iedge(k);
                    Face2Edge(3) = iedge(q);
                    node_mid(3)  = Edge(3 - w, iedge(q)); % added
                    Exit = 1;
                  end
                  if(Exit == 1),
                    break;
                  end
                end              % end w-loop
                if(Exit == 1),
                  break;
                end
              end                % end q-loop
            end
            if(Exit == 1),
              break;
            end
          end                    % end s-loop
          if(Exit == 1),
            break;
          end
        end                      % end k-loop
        for k = 1:3,
          Face(k, nface + new_face + j)      = Face2Edge(k);
          Face_flag(nface + new_face + j, 1) = Face_flag(iface,1);
	      Face_thick(nface + new_face + j, 1) = Face_thick(iface,1);
          Face_Node(1:3,nface + new_face + j) =  node_mid';
        end
      end

%   Face_Parent(1:4,iface,lvl) = [iface;
%                                 nface+new_face+1;
%                                 nface+new_face+2;
%                                 nface+new_face+3];
   CurEdge = CurEdge + 3;
   new_face= new_face+ 3;

end   % ------------------------------ end loop over faces

return