function resultfem=remove_internal_boundaries(fem)
%this function removes internal boundaries from a fem strukture and keeps the mesh
%fem=remove_internal_boundaries(fem)

 %delete interior boundaries of geometry--------------------------------------------------------------------------------
 new_geometry=geomdel(fem.geom);
 %geomplot(new_geometry);
  
 %read mesh information from fem----------------------------------------------------------------------------------------
 el=get(fem.mesh,'el');
 p=get(fem.mesh,'p');

 %write triangle element information to new el structure----------------------------------------------------------------
 for i=1:length(el)
 	if isequal(el{i}.type,'tri')
   tri_elem=el{i}.elem;
		end
	end
 el_structure{1}=struct('type','tri','elem',tri_elem);

 %create new mesh ------------------------------------------------------------------------------------------------------
 fem.mesh=femmesh(p,el_structure);
 fem.mesh=meshenrich(fem); 
 
 %the new mesh has to be modified to fit the already existing geometry==================================================
 %read new mesh information from fem------------------------------------------------------------------------------------
 clear el; clear p;
 el=get(fem.mesh,'el');
 p=get(fem.mesh,'p');
 tri=el{1};
 vtx=el{2};
 edg=el{3}
 
 %correct vtx information-----------------------------------------------------------------------------------------------
 
 %get elem vector
 elem=vtx.elem;
 
 %get according coordinates
 elem_coord=p(:,elem)';
 [n m] = size(elem_coord);
 
 %get according dom numbers
 for k = 1:n
		vtx.dom(k)=get_point(new_geometry, elem_coord(k,1),elem_coord(k,2));
 end
 
 %correct edg information-----------------------------------------------------------------------------------------------
 
 %get elem vector
 elem=edg.elem;
 
 %get according coordinates
 elem_coord_1=p(:,elem(1,:))';
 elem_coord_2=p(:,elem(2,:))';
 [n m] = size(elem_coord_1);
 
 for k = 1:n
		%get according dom numbers
		x1=elem_coord_1(k,1);
		y1=elem_coord_1(k,2);
		xy1=[x1 y1];
		x2=elem_coord_2(k,1);
		y2=elem_coord_2(k,2);
		xy2=[x2 y2];
		edg.dom(k)=get_edg_boundary(new_geometry, xy1, xy2);
		
		%calculate up-down subdomains  
			%calculate test point that lies left of boundary
			epsi=3*eps;
			if ~(abs(x2-x1)<epsi)
				test_point_x=x1+(x2-x1)/2;
				if x2-x1 > 0
					test_point_y=y1+(y2-y1)/2+epsi;
				else
					test_point_y=y1+(y2-y1)/2-epsi;
				end
			elseif ~(abs(y2-y1)<epsi)
				test_point_y=y1+(y2-y1)/2;
				if y2-y1 > 0
					test_point_x=x1+(x2-x1)/2-epsi;
				else
					test_point_x=x1+(x2-x1)/2+epsi;
				end
			else
				msgbox('Error in mesh correction.');
			end
			%test if test point lies inside the subdomain polygon and set ud numbers
				%create convex polygon 
				poly_ind=convhull(p(1,:),p(2,:));
				poly_p_x=p(1,poly_ind);
				poly_p_y=p(2,poly_ind);
			
				%test and set ud number
				edg.ud(1,k)=inpolygon(test_point_x,test_point_y,poly_p_x,poly_p_y);
				edg.ud(2,k)=~edg.ud(1,k);			
		
		%transform param values from absolute distances to relative distances	
			%get length of boundary
			XY1=flgeomed(new_geometry,edg.dom(k),0);
			XY2=flgeomed(new_geometry,edg.dom(k),1);
			bound_length=sqrt((XY2(2)-XY1(2))^2+(XY2(1)-XY1(1))^2);
			%set param values
			edg.param(:,k)=edg.param(:,k)/bound_length; 
						
		%check orientation and flip, if edg orientation is not equal boundary orientation
		if ~and(isequal(XY2(1)-XY1(1)>epsi,x2-x1>epsi),isequal(XY2(2)-XY1(2)>epsi,y2-y1>epsi))
			%flip elem
			saveelem=edg.elem(1,k);
			edg.elem(1,k)=edg.elem(2,k);
			edg.elem(2,k)=saveelem;
			%flip ud
			saveud=edg.ud(1,k);
			edg.ud(1,k)=edg.ud(2,k);
			edg.ud(2,k)=saveud;
			%calculate new param
			edg.param(1,k)=1-edg.param(1,k);
			edg.param(2,k)=1-edg.param(2,k);		
		end 
		
 end
		
 %assign corrected mesh information-------------------------------------------------------------------------------------
 el{1}=tri;
 el{2}=edg;
 el{3}=vtx;
 fem.mesh=femmesh(p,el);
 
 %assign new geometry---------------------------------------------------------------------------------------------------
 fem.geom=new_geometry;
 
 resultfem=fem;
end