# HG changeset patch
# User cdf
# Date 1433753466 0
# Node ID 4cacfa5f9470ba0e37e456f70c842fa68b22b492
# Parent  a5790043b6a9701c67e00d729094b07e6f960f38
push changes made for release

diff -r a5790043b6a9 -r 4cacfa5f9470 extra/bim/DESCRIPTION
--- a/extra/bim/DESCRIPTION	Mon Jun 01 07:38:49 2015 +0000
+++ b/extra/bim/DESCRIPTION	Mon Jun 08 08:51:06 2015 +0000
@@ -1,10 +1,10 @@
 Name: bim
-Version: 1.1.4
-Date: 2014-03-30
+Version: 1.1.5
+Date: 2014-10-17
 Author: Carlo de Falco, Culpo Massimiliano, Matteo Porro, Emanuela Abbate
 Maintainer: Carlo de Falco
 Title: PDE Solver using a Finite Element/Finite Volume approach
 Description: Package for solving Diffusion Advection Reaction (DAR) Partial Differential Equations 
-Depends: octave (>= 3.6.0), fpl, msh 
+Depends: octave (>= 3.8.0), fpl, msh 
 Autoload: no
 License: GPLv2+
diff -r a5790043b6a9 -r 4cacfa5f9470 extra/bim/NEWS
--- a/extra/bim/NEWS	Mon Jun 01 07:38:49 2015 +0000
+++ b/extra/bim/NEWS	Mon Jun 08 08:51:06 2015 +0000
@@ -1,3 +1,9 @@
+Summary of important user-visible changes for bim 1.1.5:
+-------------------------------------------------------------------
+
+ ** Improvement of the functions for stiffness matrix assembly in 2D 
+    axisymmetric configuration.
+
 Summary of important user-visible changes for bim 1.1.4:
 -------------------------------------------------------------------
 
diff -r a5790043b6a9 -r 4cacfa5f9470 extra/bim/inst/bim2a_axisymmetric_advection_diffusion.m
--- a/extra/bim/inst/bim2a_axisymmetric_advection_diffusion.m	Mon Jun 01 07:38:49 2015 +0000
+++ b/extra/bim/inst/bim2a_axisymmetric_advection_diffusion.m	Mon Jun 08 08:51:06 2015 +0000
@@ -101,14 +101,12 @@
     error("bim2a_axisymmetric_advection_diffusion: length of eta is not equal to the number of nodes.");
   endif
   
-  x = abs( mesh.p(1,:));
+  x = abs(mesh.p(1,:));
   x = x(mesh.t(1:3,:));
   y = mesh.p(2,:);
   y = y(mesh.t(1:3,:));
-
-  rcm = sum (x, 1) / (rows (mesh.t) - 1);
    
-  alphaareak = reshape (alpha.*mesh.area.*rcm',1,1,nelem);
+  alphaareak = reshape (alpha.*mesh.area,1,1,nelem);
   shg        = mesh.shg(:,:,:);
   
   ## Build local Laplacian matrix
@@ -145,6 +143,10 @@
   eta23 = etaloc(3,:) - etaloc(2,:);
   eta31 = etaloc(1,:) - etaloc(3,:);
   
+  r12 = (x(2,:) + x(1,:)) / 2;
+  r23 = (x(3,:) + x(2,:)) / 2;
+  r31 = (x(1,:) + x(3,:)) / 2;
+  
   etalocm1 = bimu_logm(etaloc(2,:),etaloc(3,:));
   etalocm2 = bimu_logm(etaloc(3,:),etaloc(1,:));
   etalocm3 = bimu_logm(etaloc(1,:),etaloc(2,:));
@@ -160,12 +162,12 @@
   [bp23,bm23] = bimu_bernoulli ((v23 - eta23) ./ etalocm1);
   [bp31,bm31] = bimu_bernoulli ((v31 - eta31) ./ etalocm2);
   
-  bp12 = reshape(gelocm3.*etalocm3.*bp12,1,1,nelem).*Lloc(1,2,:);
-  bm12 = reshape(gelocm3.*etalocm3.*bm12,1,1,nelem).*Lloc(1,2,:);
-  bp23 = reshape(gelocm1.*etalocm1.*bp23,1,1,nelem).*Lloc(2,3,:);
-  bm23 = reshape(gelocm1.*etalocm1.*bm23,1,1,nelem).*Lloc(2,3,:);
-  bp31 = reshape(gelocm2.*etalocm2.*bp31,1,1,nelem).*Lloc(3,1,:);
-  bm31 = reshape(gelocm2.*etalocm2.*bm31,1,1,nelem).*Lloc(3,1,:);
+  bp12 = reshape(r12.*gelocm3.*etalocm3.*bp12,1,1,nelem).*Lloc(1,2,:);
+  bm12 = reshape(r12.*gelocm3.*etalocm3.*bm12,1,1,nelem).*Lloc(1,2,:);
+  bp23 = reshape(r23.*gelocm1.*etalocm1.*bp23,1,1,nelem).*Lloc(2,3,:);
+  bm23 = reshape(r23.*gelocm1.*etalocm1.*bm23,1,1,nelem).*Lloc(2,3,:);
+  bp31 = reshape(r31.*gelocm2.*etalocm2.*bp31,1,1,nelem).*Lloc(3,1,:);
+  bm31 = reshape(r31.*gelocm2.*etalocm2.*bm31,1,1,nelem).*Lloc(3,1,:);
   
   Sloc(1,1,:) = (-bm12-bp31)./reshape(etaloc(1,:),1,1,nelem);
   Sloc(1,2,:) = bp12./reshape(etaloc(2,:),1,1,nelem);
@@ -310,11 +312,11 @@
 %! f = @(r,z) -D./r.*duexdr(r,z) - D.*d2uexdr2(r,z) ...
 %!           + vr./r .* uex(r,z) + vr * duexdr(r,z) ...
 %!           - D.*d2uexdz2(r,z) + vz * duexdz(r,z);
-%! rhs = bim2a_axisymmetric_rhs(mesh, ones(Nelements,1), f(mesh.p(1,:), mesh.p(2,:)));
+%! rhs = bim2a_axisymmetric_rhs(mesh, ones(Nelements,1), f(abs(mesh.p(1,:)), mesh.p(2,:)));
 %! S   = bim2a_axisymmetric_advection_diffusion(mesh,alpha,gamma,eta,beta);
-%! u   = zeros(Nnodes,1); u(Dnodes) = uex(mesh.p(1,Dnodes), mesh.p(2,Dnodes));
+%! u   = zeros(Nnodes,1); u(Dnodes) = uex(abs(mesh.p(1,Dnodes)), mesh.p(2,Dnodes));
 %! u(Varnodes) = S(Varnodes,Varnodes)\(rhs(Varnodes) - S(Varnodes,Dnodes)*u(Dnodes));
-%! assert(u,uex(mesh.p(1,:), mesh.p(2,:))',1e-7)
+%! assert(u,uex(abs(mesh.p(1,:)), mesh.p(2,:))',1e-7)
 
 %!test
 %! n         = 10;
@@ -335,11 +337,11 @@
 %! eta       = ones(Nnodes,1);
 %! beta      = 1/D * 1/2*(mesh.p(1,:)').^2;
 %! f         = @(r,z) 1./r.*(1+r).*(r-D) .* uex(r,z);
-%! rhs = bim2a_axisymmetric_rhs(mesh, ones(Nelements,1), f(mesh.p(1,:), mesh.p(2,:)));
+%! rhs = bim2a_axisymmetric_rhs(mesh, ones(Nelements,1), f(abs(mesh.p(1,:)), mesh.p(2,:)));
 %! S   = bim2a_axisymmetric_advection_diffusion(mesh,alpha,gamma,eta,beta);
-%! u   = zeros(Nnodes,1); u(Dnodes) = uex(mesh.p(1,Dnodes), mesh.p(2,Dnodes));
+%! u   = zeros(Nnodes,1); u(Dnodes) = uex(abs(mesh.p(1,Dnodes)), mesh.p(2,Dnodes));
 %! u(Varnodes) = S(Varnodes,Varnodes)\(rhs(Varnodes) - S(Varnodes,Dnodes)*u(Dnodes));
-%! assert(u,uex(mesh.p(1,:), mesh.p(2,:))',1e-3)
+%! assert(u,uex(abs(mesh.p(1,:)), mesh.p(2,:))',1e-3)
 
 %!test
 %! [mesh] = msh2m_structured_mesh([0:.1:1],[0:.1:1],1,1:4);
diff -r a5790043b6a9 -r 4cacfa5f9470 extra/bim/inst/bim2a_axisymmetric_advection_upwind.m
--- a/extra/bim/inst/bim2a_axisymmetric_advection_upwind.m	Mon Jun 01 07:38:49 2015 +0000
+++ b/extra/bim/inst/bim2a_axisymmetric_advection_upwind.m	Mon Jun 08 08:51:06 2015 +0000
@@ -67,9 +67,7 @@
   y = mesh.p(2,:);
   y = y(mesh.t(1:3,:));
   
-  rcm = sum (x, 1) / (rows (mesh.t) - 1);
-
-  alphaareak = reshape (mesh.area .* rcm', 1, 1, nelem);
+  alphaareak = reshape (mesh.area, 1, 1, nelem);
   shg        = mesh.shg(:,:,:);
   
   ## Build local Laplacian matrix	
@@ -104,12 +102,16 @@
   [bp23, bm23] = deal (- (v23 - abs (v23))/2, (v23 + abs (v23))/2);
   [bp31, bm31] = deal (- (v31 - abs (v31))/2, (v31 + abs (v31))/2);
   
-  bp12 = reshape(bp12,1,1,nelem).*Lloc(1,2,:);
-  bm12 = reshape(bm12,1,1,nelem).*Lloc(1,2,:);
-  bp23 = reshape(bp23,1,1,nelem).*Lloc(2,3,:);
-  bm23 = reshape(bm23,1,1,nelem).*Lloc(2,3,:);
-  bp31 = reshape(bp31,1,1,nelem).*Lloc(3,1,:);
-  bm31 = reshape(bm31,1,1,nelem).*Lloc(3,1,:);
+  r12 = (x(2,:) + x(1,:)) / 2;
+  r23 = (x(3,:) + x(2,:)) / 2;
+  r31 = (x(1,:) + x(3,:)) / 2;
+
+  bp12 = reshape(r12 .* bp12,1,1,nelem).*Lloc(1,2,:);
+  bm12 = reshape(r12 .* bm12,1,1,nelem).*Lloc(1,2,:);
+  bp23 = reshape(r23 .* bp23,1,1,nelem).*Lloc(2,3,:);
+  bm23 = reshape(r23 .* bm23,1,1,nelem).*Lloc(2,3,:);
+  bp31 = reshape(r31 .* bp31,1,1,nelem).*Lloc(3,1,:);
+  bm31 = reshape(r31 .* bm31,1,1,nelem).*Lloc(3,1,:);
   
   Sloc(1,1,:) = (-bm12-bp31);
   Sloc(1,2,:) = bp12;
diff -r a5790043b6a9 -r 4cacfa5f9470 extra/bim/inst/bim2a_axisymmetric_rhs.m
--- a/extra/bim/inst/bim2a_axisymmetric_rhs.m	Mon Jun 01 07:38:49 2015 +0000
+++ b/extra/bim/inst/bim2a_axisymmetric_rhs.m	Mon Jun 08 08:51:06 2015 +0000
@@ -45,7 +45,7 @@
   ## Check input
   if (nargin != 3)
     error("bim2a_axisymmetric_rhs: wrong number of input parameters.");
-  elseif !(isstruct(mesh)     && isfield(mesh,"p") &&
+  elseif !(isstruct(mesh) && isfield(mesh,"p") &&
 	  isfield (mesh,"t") && isfield(mesh,"e"))
     error("bim2a_axisymmetric_rhs: first input is not a valid mesh structure.");
   elseif !(all(mesh.p(1,:) >= 0) || all(mesh.p(1,:) <= 0))
diff -r a5790043b6a9 -r 4cacfa5f9470 extra/bim/inst/bim2c_intrp.m
--- a/extra/bim/inst/bim2c_intrp.m	Mon Jun 01 07:38:49 2015 +0000
+++ b/extra/bim/inst/bim2c_intrp.m	Mon Jun 08 08:51:06 2015 +0000
@@ -18,9 +18,9 @@
 ##
 ## @deftypefn {Function File} {@var{data}} = bim2c_intrp (@var{msh}, @var{n_data}, @var{e_data}, @var{points}) 
 ##
-## Compute interpolated values of node centered multicomponent node centered field @var{n_data} and 
-## cell centered field @var{n_data} at an arbitrary set of points whos coordinates are given in the 
-## n_by_2 matrix  @var{points}. 
+## Compute interpolated values of multicomponent node centered field @var{n_data} and/or
+## cell centered field @var{n_data} at an arbitrary set of points whose coordinates are given in the 
+## n_by_2 matrix @var{points}. 
 ##
 ## @end deftypefn
 
@@ -59,4 +59,4 @@
 %! x = y = linspace (0, 1, 100).';
 %! u = msh.p(1, :).';
 %! ui = bim2c_intrp (msh, u, [], [x, y]); 
-%! assert (ui, linspace (0, 1, 100), 10*eps);
\ No newline at end of file
+%! assert (ui, linspace (0, 1, 100), 10*eps);