function index = get_edg_boundary(geometry,p1,p2,tol)
%this function finds the boundary index of a edge element between the
%points p1 = [x1 y1] and p2 =[x2 y2]
%index = get_edg_boundary(geometry,p1,p2)

 %set default tolerance=================================================================================================
 if nargin < 4 
    tol = eps;
 end
 
 % plot geometry =======================================================================================================
 %geomplot(geometry,'pointlabels','on','edgelabels','on');

 %get number of objects=================================================================================================
  gd = geominfo(geometry,'out',{'gd'});
 no = geominfo(geometry,'out',{'no'},'od',0:gd);
 
 %get coordinates of edges==============================================================================================
 X1=0;
 Y1=0;
 X2=0;
 Y2=0;

 for k = 1:no
  xy=flgeomed(geometry,k,0);
  X1(k)=xy(1);
  Y1(k)=xy(2);
  xy2=flgeomed(geometry,k,1);
  X2(k)=xy2(1);
  Y2(k)=xy2(2);  
 end


 %%calculate all boundary elements that may include the edg element in x-direction
 
  % getting elements with X1 <= x1, X1 <= x2 
  element_x_s1 = find(X1 <=  p1(1)+tol);
  element_x_s2 = find(X1 <=  p2(1)+tol);
  

  % getting elements with X2 >= x1, X2 >=x2 
  element_x_g1 = find(X2 >=  p1(1)-tol);
  element_x_g2 = find(X2 >=  p2(1)-tol);
  

		% creating matrix for start_x smaller x_vals
  coordinate_matrix_x_s1 = zeros(1,no(1));
  coordinate_matrix_x_s2 = zeros(1,no(1));
		
  % creating matrix for end_x greater x_vals
  coordinate_matrix_x_g1 = zeros(1,no(1));
  coordinate_matrix_x_g2 = zeros(1,no(1));
		
		% filling matrix with a one when element fullfills the smaller condition
  coordinate_matrix_x_s1(element_x_s1) = 1;
  coordinate_matrix_x_s2(element_x_s2) = 1;
  coordinate_matrix_x_s = coordinate_matrix_x_s1 & coordinate_matrix_x_s2;

  % filling matrix with a one when element fullfills the greater condition
  coordinate_matrix_x_g1(element_x_g1) = 1;
  coordinate_matrix_x_g2(element_x_g2) = 1;
  coordinate_matrix_x_g = coordinate_matrix_x_g1 & coordinate_matrix_x_g2;

  % Merge the x-coordinate_s and x-coordinate_g matrix  
  coordinate_matrix_x = coordinate_matrix_x_g & coordinate_matrix_x_s;
  
 %% calculate all boundary elements that may include the edg element in y-direction

  % getting elements with Y1 <= y1, Y1 <= y2 
  element_y_s1 = find(Y1 <=  p1(2)+tol);
  element_y_s2 = find(Y1 <=  p2(2)+tol);

  % getting elements with Y2 >= x1, Y2 >=x2
  element_y_g1 = find(Y2 >=  p1(2)-tol);
  element_y_g2 = find(Y2 >=  p2(2)-tol);

  % creating matrix for start_y smaller y_vals
  coordinate_matrix_y_s1 = zeros(1,no(1));
  coordinate_matrix_y_s2 = zeros(1,no(1));

  % creating matrix for end_x greater x_vals
  coordinate_matrix_y_g1 = zeros(1,no(1));
  coordinate_matrix_y_g2 = zeros(1,no(1));
		
		% filling matrix with a one when element fullfills the smaller condition
  coordinate_matrix_y_s1(element_y_s1) = 1;
  coordinate_matrix_y_s2(element_y_s2) = 1;

  % filling matrix with a one when element fullfills the greater condition
  coordinate_matrix_y_g1(element_y_g1) = 1;
  coordinate_matrix_y_g2(element_y_g2) = 1;

  % Merge the x-coordinate_s and x-coordinate_g matrix
  coordinate_matrix_y_g = coordinate_matrix_y_g1 & coordinate_matrix_y_g2;
  coordinate_matrix_y_s = coordinate_matrix_y_s1 & coordinate_matrix_y_s2;
  
  coordinate_matrix_y = coordinate_matrix_y_g & coordinate_matrix_y_s; 

 %% Calculating the Object which fullfills all the conditions
 
 % Addition of both matrices to get the element which is the searchobject
 index = find(coordinate_matrix_x & coordinate_matrix_y == 1);

 if length(index) > 1
  msgbox(['Error: More than one possible boundary found.' sprintf('\n') ...
           'for the edg element between p1 = (' num2str(p1(1),3) ', ' num2str(p1(2),3) ') and' sprintf('\n') ...
	   '        p2 = (' num2str(p2(1),3) ', ' num2str(p2(2),3) ') !'],'Error','error');
 end 
 if length(index) < 1
  msgbox(['Error: No boundary found for the edg element' sprintf('\n') ...
           'between p1 = (' num2str(p1(1),3) ', ' num2str(p1(2),3) ') and' sprintf('\n') ...
	   '        p2 = (' num2str(p2(1),3) ', ' num2str(p2(2),3) ') !'],'Error','error');
 end 
end