--- a/extra/bim/INDEX	Fri Sep 13 13:57:58 2013 +0000
+++ b/extra/bim/INDEX	Mon Sep 16 10:40:48 2013 +0000
@@ -17,6 +17,7 @@
   bim2a_rhs
   bim3a_rhs
   bim2a_boundary_mass
+  bim3a_boundary_mass
 Pre-processing and Post-processing computations
   bim2c_mesh_properties
   bim3c_mesh_properties

--- a/extra/bim/inst/bim2c_tri_to_nodes.m	Fri Sep 13 13:57:58 2013 +0000
+++ b/extra/bim/inst/bim2c_tri_to_nodes.m	Mon Sep 16 10:40:48 2013 +0000
@@ -34,13 +34,13 @@
 
 function [u_nod, m_tri] = bim2c_tri_to_nodes (m, u_tri)
 
-  if (nargout >1)
+  if (nargout > 1)
     if (isstruct (m))
-      nel  = columns (m.t);
-      nnod = columns (m.p);
-      ii = m.t(1:3, :);
-      jj = repmat (1:nel, 3, 1);
-      vv = repmat (m.area(:)', 3, 1) / 3;
+      nel   = columns (m.t);
+      nnod  = columns (m.p);
+      ii    = m.t(1:3, :);
+      jj    = repmat (1:nel, 3, 1);
+      vv    = repmat (m.area(:)', 3, 1) / 3;
       m_tri = bim2a_reaction (m, 1, 1) \ sparse (ii, jj, vv, nnod, nel); 
     elseif (ismatrix (m))
       m_tri = m;
@@ -49,12 +49,17 @@
     endif
     u_nod = m_tri * u_tri;
   else
-    rhs  = bim2a_rhs (m, u_tri, 1);
-    mass = bim2a_reaction (m, 1, 1);
-    u_nod = full (mass \ rhs);
+    if (isstruct (m))
+      rhs   = bim2a_rhs (m, u_tri, 1);
+      mass  = bim2a_reaction (m, 1, 1);
+      u_nod = full (mass \ rhs);
+    elseif (ismatrix (m))
+      u_nod = m * u_tri;
+    else
+      error ("bim2c_tri_to_nodes: first input parameter is of incorrect type");
+    endif
   endif
 
-
 endfunction
 
 %!test

--- a/extra/bim/inst/bim3a_osc_laplacian.m	Fri Sep 13 13:57:58 2013 +0000
+++ b/extra/bim/inst/bim3a_osc_laplacian.m	Mon Sep 16 10:40:48 2013 +0000
@@ -32,8 +32,7 @@
 ##
 ## - div (@var{epsilon} grad (u)) = f
 ## 
-## where @var{epsilon} is an element-wise constant scalar function,
-## while @var{kappa} is a piecewise linear conforming scalar function.
+## where @var{epsilon} is an element-wise constant scalar function.
 ##
 ## @seealso{bim3a_rhs, bim3a_reaction, bim2a_laplacian, bim3a_laplacian}
 ## @end deftypefn
@@ -47,7 +46,7 @@
              && isfield (msh, "p") 
              && isfield (msh, "t") 
              && isfield (msh, "e")))
-    error (["bim3a_laplacian: first input ", ...
+    error (["bim3a_osc_laplacian: first input ", ...
             "is not a valid msh structure"]);
   endif
 

--- a/extra/bim/inst/bim3c_global_flux.m	Fri Sep 13 13:57:58 2013 +0000
+++ b/extra/bim/inst/bim3c_global_flux.m	Mon Sep 16 10:40:48 2013 +0000
@@ -29,7 +29,7 @@
 ##
 ## The vector field is defined as:
 ##
-## F =- @var{alpha} ( grad (u) - gard (@var{v}) u ) 
+## F =- @var{alpha} ( grad (u) - grad (@var{v}) u ) 
 ##
 ## where @var{v} is a piecewise linear continuous scalar
 ## functions and @var{alpha} is a piecewise constant scalar function.
@@ -40,50 +40,84 @@
 
 function F = bim3c_global_flux (mesh, u, acoeff, v)
 
-  t = mesh.t;
-  nelem = columns(mesh.t);
-  F = zeros (3, nelem);
+  t     = mesh.t;
+  nelem = columns (mesh.t);
+  F     = zeros (3, nelem);
 
-  for iel = 1:nelem
-    ## Local matrices
-    Lloc = zeros(4,4); ## diffusion
-    Sloc = zeros(4,4); ## advection-diffusion
+  ## Local contributions
+  Lloc = zeros (4,4,nelem);
+
+  epsilonareak = reshape (acoeff .* mesh.area', 1, 1, nelem);
+  shg = mesh.shg(:,:,:);
 
-    epsilonareak = acoeff(iel) .* mesh.area(iel);
-    shg = mesh.shg(:,:,iel);    
-    vloc = v(t(1:4, iel));
-    uloc = u(t(1:4, iel));
-
-    ## Computation
-    for inode = 1:4
-      for jnode = 1:4
-        Lloc(inode, jnode)   = sum (shg(:, inode) .* shg(:, jnode)) .* epsilonareak;
-      endfor
-    endfor	
+  ## Computation
+  for inode = 1:4
+    for jnode = 1:4
+      ginode(inode,jnode,:) = t(inode,:);
+      gjnode(inode,jnode,:) = t(jnode,:);
+      Lloc(inode,jnode,:)   = sum (shg(:,inode,:) .* shg(:,jnode,:), 1) ...
+                              .* epsilonareak;
+    endfor
+  endfor
 
-    ## SGloc=[...
-    ##        -bm12-bm13-bm14,bp12            ,bp13           ,bp14     
-    ##        bm12           ,-bp12-bm23-bm24 ,bp23           ,bp24
-    ##        bm13           ,bm23            ,-bp13-bp23-bm34,bp34
-    ##        bm14           ,bm24            ,bm34           ,-bp14-bp24-bp34...
-    ##        ];
-    
-    for inode = 1:4
-      for jnode = inode+1:4
-        [Sloc(inode, jnode), Sloc(jnode, inode)] = bimu_bernoulli (vloc(jnode)-vloc(inode));
-        Sloc(inode, jnode) *= Lloc(inode, jnode);
-        Sloc(jnode, inode) *= Lloc(inode, jnode);
-      endfor
-    endfor
+  uloc = u(t(1:4, :));
+  vloc = v(t(1:4, :));
+  [bp12,bm12] = bimu_bernoulli (vloc(2,:)-vloc(1,:));
+  [bp13,bm13] = bimu_bernoulli (vloc(3,:)-vloc(1,:));
+  [bp14,bm14] = bimu_bernoulli (vloc(4,:)-vloc(1,:));
+  [bp23,bm23] = bimu_bernoulli (vloc(3,:)-vloc(2,:));
+  [bp24,bm24] = bimu_bernoulli (vloc(4,:)-vloc(2,:));
+  [bp34,bm34] = bimu_bernoulli (vloc(4,:)-vloc(3,:));
+  
+  bp12 = reshape (bp12, 1, 1, nelem) .* Lloc(1,2,:);
+  bm12 = reshape (bm12, 1, 1, nelem) .* Lloc(1,2,:);
+  bp13 = reshape (bp13, 1, 1, nelem) .* Lloc(1,3,:);
+  bm13 = reshape (bm13, 1, 1, nelem) .* Lloc(1,3,:);
+  bp14 = reshape (bp14, 1, 1, nelem) .* Lloc(1,4,:);
+  bm14 = reshape (bm14, 1, 1, nelem) .* Lloc(1,4,:);
+  bp23 = reshape (bp23, 1, 1, nelem) .* Lloc(2,3,:);
+  bm23 = reshape (bm23, 1, 1, nelem) .* Lloc(2,3,:);
+  bp24 = reshape (bp24, 1, 1, nelem) .* Lloc(2,4,:);
+  bm24 = reshape (bm24, 1, 1, nelem) .* Lloc(2,4,:);
+  bp34 = reshape (bp34, 1, 1, nelem) .* Lloc(3,4,:);
+  bm34 = reshape (bm34, 1, 1, nelem) .* Lloc(3,4,:);
+  
+  ## SGloc=[...
+  ##        -bm12-bm13-bm14,bp12            ,bp13           ,bp14     
+  ##        bm12           ,-bp12-bm23-bm24 ,bp23           ,bp24
+  ##        bm13           ,bm23            ,-bp13-bp23-bm34,bp34
+  ##        bm14           ,bm24            ,bm34           ,-bp14-bp24-bp34
+  ##        ];
+  
+  Sloc(1,1,:) = -bm12-bm13-bm14;
+  Sloc(1,2,:) = bp12;
+  Sloc(1,3,:) = bp13;
+  Sloc(1,4,:) = bp14;
 
-    for inode = 1:4
-      Sloc(inode, inode) = - sum (Sloc (:, inode));
-    endfor
+  Sloc(2,1,:) = bm12;
+  Sloc(2,2,:) = -bp12-bm23-bm24; 
+  Sloc(2,3,:) = bp23;
+  Sloc(2,4,:) = bp24;
 
-   r = Sloc * uloc;
-   f = Lloc(1:3, 1:3) \ r(1:3);
+  Sloc(3,1,:) = bm13;
+  Sloc(3,2,:) = bm23;
+  Sloc(3,3,:) = -bp13-bp23-bm34;
+  Sloc(3,4,:) = bp34;
+  
+  Sloc(4,1,:) = bm14;
+  Sloc(4,2,:) = bm24;
+  Sloc(4,3,:) = bm34;
+  Sloc(4,4,:) = -bp14-bp24-bp34;
 
-   F(:,iel) = shg(:,1:3) * f;
+  r = zeros (4, nelem);
+  f = zeros (3, nelem);
+  
+  for iel = 1:nelem
+
+   r(:,iel) = Sloc(:,:,iel) * uloc(:,iel);
+   f(:,iel) = Lloc(1:3, 1:3, iel) \ r(1:3, iel);
+
+   F(:,iel) = shg(:,1:3, iel) * f(:, iel);
 
   endfor
 
@@ -95,4 +129,4 @@
 %! v = ones (N^3, 1);
 %! alpha = ones (columns (msh.t), 1);
 %! F =  bim3c_global_flux (msh, u, alpha, v);
-%! assert (norm (F(:), inf), 0, 100*eps);
\ No newline at end of file
+%! assert (norm (F(:), inf), 0, 100*eps);

--- a/extra/bim/inst/bim3c_tri_to_nodes.m	Fri Sep 13 13:57:58 2013 +0000
+++ b/extra/bim/inst/bim3c_tri_to_nodes.m	Mon Sep 16 10:40:48 2013 +0000
@@ -34,7 +34,7 @@
 
 function [u_nod, m_tri] = bim3c_tri_to_nodes (m, u_tri)
 
-  if (nargout >1)
+  if (nargout > 1)
     if (isstruct (m))
       nel  = columns (m.t);
       nnod = columns (m.p);
@@ -49,12 +49,17 @@
     endif
     u_nod = m_tri * u_tri;
   else
-    rhs  = bim3a_rhs (m, u_tri, 1);
-    mass = bim3a_reaction (m, 1, 1);
-    u_nod = full (mass \ rhs);
+    if (isstruct (m))
+      rhs   = bim3a_rhs (m, u_tri, 1);
+      mass  = bim3a_reaction (m, 1, 1);
+      u_nod = full (mass \ rhs);
+    elseif (ismatrix (m))
+      u_nod = m * u_tri;
+    else
+      error ("bim3c_tri_to_nodes: first input parameter is of incorrect type");
+    endif
   endif
 
-
 endfunction
 
 %!test