* cleaned code up a little

* added documentation (partially)
* added support for interpolation in 2D

pdeutils/fespaces/lagrangeBasis.m

 function B = lagrangeBasis(dim, order)
+%lagrangeBasis Initializes a lagrange basis.
+%
+% Input:
+%   dim     Dimension of the domain over which the basis is defined.
+%   order   Polynomial order of the basis.
+%
+% Return:
+%   B       A structure representing the basis.
+%
+% Todo:
+%   Explain the basis struct.
 
-    if ((order==1) & (dim==1))
+    if ((order==1) && (dim==1))
         B = struct( ...
             'dim', 1, ...
             'evaluate', @D1P1Evaluate, ...
@@ -10,7 +21,7 @@ function B = lagrangeBasis(dim, order)
             'order', @(grid, elementIndex)(1),...
             'index', @P1Index,...
             'size', @(grid)(size(grid.nodes,2)));
-    elseif ((order==1) & (dim==2))
+    elseif ((order==1) && (dim==2))
         B = struct( ...
             'dim', 2, ...
             'evaluate', @D2P1Evaluate, ...
@@ -29,15 +40,28 @@ end
 
 
 
-function u = P1LocalInterpolation(grid, elementIndex, f, df, ddf, dddf)
-    element = grid.nodes(grid.elements(:,elementIndex));
-    n = size(element,2);
-    u = zeros(1,n);
-    for i=1:n
-        u(i) = f(element(:,i));
+function u = P1LocalInterpolation(grid, elementIndex, f, varargin)
+% IMPORTANT: f is expected to be given in global coordinates
+% Todo: docme
+
+    % Get the dimension that we are currently working in.
+    dim     = size(grid.nodes,1);
+    n       = dim+1;
+
+    % Retrieve the corners of the associated element.
+    % There are n=dim+1 corners per element.
+    cornerIds   = grid.elements(1:n, elementIndex);
+    corners     = grid.nodes(:,cornerIds);
+    
+    % Evaluate f in all those corners and return the resulting vector.
+    u   = zeros(n,1);
+    for i = 1:n
+        u(i)    = f(corners(:,i));
     end
+
 end
 
+
 function globalIndices = P1Index(grid, elementIndex, localIndices)
 % Compute global indices of basis functions with given localIndices in element.
     dim = size(grid.nodes,1);
@@ -49,9 +73,9 @@ end
 
 function y = D1P1Evaluate(grid, elementIndex, localIndex, x)
     switch localIndex
-        case 1,
+        case 1
             y = x;
-        case 2,
+        case 2
             y = 1-x;
     end
 end
@@ -59,9 +83,9 @@ end
 function Dx = D1P1EvaluateJacobian(grid, elementIndex, localIndex, x)
     Dx = x*0;
     switch localIndex
-        case 1,
+        case 1
             Dx(1,:) = 1;
-        case 2 ,
+        case 2
             Dx(1,:) = -1;
     end
 end
@@ -69,11 +93,11 @@ end
 
 function y = D2P1Evaluate(grid, elementIndex, localIndex, x)
     switch localIndex
-        case 1,
+        case 1
             y = 1 - x(1,:) - x(2,:);
-        case 2 ,
+        case 2
             y = x(1,:);
-        case 3,
+        case 3
             y = x(2,:);
     end
 end
@@ -81,13 +105,13 @@ end
 function Dx = D2P1EvaluateJacobian(grid, elementIndex, localIndex, x)
     Dx = x*0;
     switch localIndex
-        case 1,
+        case 1
             Dx(1,:) = -1;
             Dx(2,:) = -1;
-        case 2 ,
+        case 2
             Dx(1,:) = 1;
             Dx(2,:) = 0;
-        case 3,
+        case 3
             Dx(1,:) = 0;
             Dx(2,:) = 1;
     end