diff --git a/scripts/area_on_point b/scripts/area_on_point
new file mode 100644
index 0000000000000000000000000000000000000000..817e4978c17680cfffae140ddb2815cf545222aa
--- /dev/null
+++ b/scripts/area_on_point
@@ -0,0 +1,56 @@
+function [ a ] = area_on_point( polygon,Polygon,Point )
+% Calculates intersection areas of convex polygons
+% For a small polygon surrounding a point in a much coarser triangulation,
+% this function calculates the intersection area with the coarse polygon
+% Pol. Point is one of the points defining Pol. Note that both size
+% difference and convexity assumptions must hold!
+% f1 er point 1 og point 2, f2 er p2 og p3 ... fn er pn og p1
+N       = size(Polygon,1);
+n       = size(polygon,1);
+inPol   = find(ismember(Polygon,Point,'rows'));
+
+% Find the faces of Pol adjacent to Point
+if inPol == 1
+    FirstFace = [Polygon(N,:);Polygon(1,:)];
+else
+    FirstFace = [Polygon(inPol-1,:);Polygon(inPol,:)];
+end
+if inPol == N
+    SecondFace = [Polygon(N,:);Polygon(1,:)];
+else
+    SecondFace = [Polygon(inPol,:);Polygon(inPol+1,:)];
+end
+
+pol2 = [polygon;polygon(1,:)];
+for i = 1:n
+    face          = [pol2(i,:);pol2(i+1,:)];
+    intersection1 = intersection_between_inclusive(FirstFace, ...
+                FirstFace(1)==FirstFace(2), face, face(1)==face(2),true);
+    if ~isempty(intersection1)
+        break
+    end
+end
+for j = 1:n
+    face          = [pol2(j,:);pol2(j+1,:)];
+    intersection2 = intersection_between_inclusive(SecondFace, ...
+                SecondFace(1)==SecondFace(2), face, face(1)==face(2),true);
+    if ~isempty(intersection2)
+        break
+    end
+end
+if i ~= j
+    smallPointsInside = inpolygon(polygon(:,1),polygon(:,2),Polygon(:,1),Polygon(:,2));
+    smallPointsInside = polygon(find(smallPointsInside),:);
+else
+    smallPointsInside = [];
+end
+newPoints   = [Point; intersection1;intersection2;smallPointsInside];
+
+
+if size( unique( [ newPoints(:,1) newPoints(:,2) ], 'rows' ), 1 ) > 2
+    [~, a]      = convhull(newPoints(:,1),newPoints(:,2));
+else
+    a = 0;
+end
+
+end