## [Octave-cvsupdate] SF.net SVN: octave:[8193] trunk/octave-forge/extra/bim/inst

 [Octave-cvsupdate] SF.net SVN: octave:[8193] trunk/octave-forge/extra/bim/inst From: - 2011-03-31 12:23:40 Revision: 8193 http://octave.svn.sourceforge.net/octave/?rev=8193&view=rev Author: cdf Date: 2011-03-31 12:23:34 +0000 (Thu, 31 Mar 2011) Log Message: ----------- 3d bugs solved Modified Paths: -------------- trunk/octave-forge/extra/bim/inst/bim3a_advection_diffusion.m trunk/octave-forge/extra/bim/inst/bim3a_rhs.m Modified: trunk/octave-forge/extra/bim/inst/bim3a_advection_diffusion.m =================================================================== --- trunk/octave-forge/extra/bim/inst/bim3a_advection_diffusion.m 2011-03-29 21:00:46 UTC (rev 8192) +++ trunk/octave-forge/extra/bim/inst/bim3a_advection_diffusion.m 2011-03-31 12:23:34 UTC (rev 8193) @@ -22,7 +22,7 @@ ## ## @deftypefn {Function File} @ ## {[@var{A}]} = @ -## bim3a_advection_diffusion(@var{mesh},@var{alpha}, @var{v}) +## bim3a_advection_diffusion (@var{mesh}, @var{alpha}, @var{v}) ## ## Build the Scharfetter-Gummel stabilized stiffness matrix for a ## diffusion-advection problem. @@ -70,7 +70,7 @@ endfor endfor - vloc = v(t(:,iel)); + 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,:)); Modified: trunk/octave-forge/extra/bim/inst/bim3a_rhs.m =================================================================== --- trunk/octave-forge/extra/bim/inst/bim3a_rhs.m 2011-03-29 21:00:46 UTC (rev 8192) +++ trunk/octave-forge/extra/bim/inst/bim3a_rhs.m 2011-03-31 12:23:34 UTC (rev 8193) @@ -21,7 +21,7 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{b}]} = @ -## bim3a_rhs(@var{mesh},@var{f},@var{g}) +## bim3a_rhs (@var{mesh}, @var{f}, @var{g}) ## ## Build the finite element right-hand side of a diffusion problem ## employing mass-lumping. @@ -57,8 +57,8 @@ error("bim3a_rhs: first input is not a valid mesh structure."); endif - nnodes = length(p); - nelem = length(t); + nnodes = columns (mesh.p); + nelem = length (mesh.t); ## Turn scalar input to a vector of appropriate size if isscalar(f) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 

 [Octave-cvsupdate] SF.net SVN: octave:[7461] trunk/octave-forge/extra/bim/inst From: - 2010-07-01 14:44:34 Revision: 7461 http://octave.svn.sourceforge.net/octave/?rev=7461&view=rev Author: cdf Date: 2010-07-01 14:44:27 +0000 (Thu, 01 Jul 2010) Log Message: ----------- pure advection Added Paths: ----------- trunk/octave-forge/extra/bim/inst/bim1a_advection_upwind.m trunk/octave-forge/extra/bim/inst/bim2a_advection_upwind.m Added: trunk/octave-forge/extra/bim/inst/bim1a_advection_upwind.m =================================================================== --- trunk/octave-forge/extra/bim/inst/bim1a_advection_upwind.m (rev 0) +++ trunk/octave-forge/extra/bim/inst/bim1a_advection_upwind.m 2010-07-01 14:44:27 UTC (rev 7461) @@ -0,0 +1,87 @@ +## Copyright (C) 2010 Carlo de Falco +## +## This file is part of: +## BIM - Diffusion Advection Reaction PDE Solver +## +## BIM is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## BIM is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with BIM; If not, see ;. +## +## author: Carlo de Falco +## author: Massimiliano Culpo + +## -*- texinfo -*- +## +## @deftypefn {Function File} @ +## {[@var{A}]} = bim1a_advection_upwind (@var{mesh}, @var{beta}) +## +## Build the UW stabilized stiffness matrix for an advection problem. +## +## The equation taken into account is: +## +## @iftex +## @tex +## $(\vect{beta} u)' = f$ +## @end tex +## @end iftex +## @ifinfo +## (@var{beta} u)' = f +## @end ifinfo +## @ifhtml +## (@var{beta} u ))' = f +## @end ifhtml +## +## where @var{beta} is an element-wise constant. +## +## Instead of passing the vector field @var{beta} directly one can pass +## a piecewise linear conforming scalar function @var{phi} as the last +## input. In such case @var{beta} = grad @var{phi} is assumed. +## +## If @var{phi} is a single scalar value @var{beta} is assumed to be 0 +## in the whole domain. +## +## @seealso{bim1a_rhs, bim1a_reaction, bim1a_laplacian, +## bim2a_advection_diffusion} +## @end deftypefn + +function A = bim1a_advection_upwind (x, beta) + + ## Check input + if nargin != 2 + error("bim1a_advection_upwind: wrong number of input parameters."); + endif + + nnodes = length(x); + nelem = nnodes-1; + + areak = reshape(diff(x),[],1); + + if (length(beta) == 1) + vk = 0; + elseif (length(beta) == nelem) + vk = beta .* areak; + elseif (length(beta) == nnodes) + vk = diff(beta); + else + error("bim1a_advection_upwind: coefficient beta has wrong dimensions."); + endif + + ck = 2./([0; areak(2:end)]+[areak(1:end-1); 0]); + bmk = (vk+abs(vk))/2; + bpk = -(vk-abs(vk))/2; + + dm1 = [-(ck.*bmk); NaN]; + dp1 = [NaN; -(ck.*bpk)]; + d0 = [(ck(1).*bmk(1)); ((ck.*bmk)(2:end) + (ck.*bpk)(1:end-1)); (ck(end).*bpk(end))]; + A = spdiags([dm1, d0, dp1],-1:1,nnodes,nnodes); + +endfunction Added: trunk/octave-forge/extra/bim/inst/bim2a_advection_upwind.m =================================================================== --- trunk/octave-forge/extra/bim/inst/bim2a_advection_upwind.m (rev 0) +++ trunk/octave-forge/extra/bim/inst/bim2a_advection_upwind.m 2010-07-01 14:44:27 UTC (rev 7461) @@ -0,0 +1,131 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## BIM - Diffusion Advection Reaction PDE Solver +## +## BIM is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## BIM is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with BIM; If not, see ;. +## +## author: Carlo de Falco +## author: Massimiliano Culpo + +## -*- texinfo -*- +## +## @deftypefn {Function File} @ +## {[@var{A}]} = @ +## bim2a_advection_upwind (@var{mesh}, @var{beta}) +## +## Build the UW stabilized stiffness matrix for an advection problem. +## +## The equation taken into account is: +## +## @iftex +## @tex +## $(\vect{beta} u )' = f$ +## @end tex +## @end iftex +## @ifinfo +## (@var{beta} u )' = f +## @end ifinfo +## @ifhtml +## (@var{beta} u )' = f +## @end ifhtml +## +## where @var{beta} is an element-wise constant vector function. +## +## Instead of passing the vector field @var{beta} directly one can pass +## a piecewise linear conforming scalar function @var{phi} as the last +## input. In such case @var{beta} = grad @var{phi} is assumed. +## +## If @var{phi} is a single scalar value @var{beta} is assumed to be 0 +## in the whole domain. +## +## @seealso{bim2a_rhs, bim2a_reaction, bim2c_mesh_properties} +## @end deftypefn + +function A = bim2a_advection_upwind (mesh, beta) + + ## Check input + if nargin != 2 + error("bim2a_advection_upwind: wrong number of input parameters."); + elseif !(isstruct(mesh) && isfield(mesh,"p") && + isfield (mesh,"t") && isfield(mesh,"e")) + error("bim2a_advection_upwind: first input is not a valid mesh structure."); + endif + + nnodes = columns(mesh.p); + nelem = columns(mesh.t); + + alphaareak = reshape (mesh.area, 1, 1, nelem); + shg = mesh.shg(:,:,:); + + ## Build local Laplacian matrix + Lloc = zeros(3,3,nelem); + + for inode = 1:3 + for jnode = 1:3 + ginode(inode,jnode,:) = mesh.t(inode,:); + gjnode(inode,jnode,:) = mesh.t(jnode,:); + Lloc(inode,jnode,:) = sum( shg(:,inode,:) .* shg(:,jnode,:),1) .* alphaareak; + endfor + endfor + + x = mesh.p(1,:); + x = x(mesh.t(1:3,:)); + y = mesh.p(2,:); + y = y(mesh.t(1:3,:)); + + if all(size(beta)==1) + v12 = 0; + v23 = 0; + v31 = 0; + elseif all(size(beta)==[2,nelem]) + v12 = beta(1,:) .* (x(2,:)-x(1,:)) + beta(2,:) .* (y(2,:)-y(1,:)); + v23 = beta(1,:) .* (x(3,:)-x(2,:)) + beta(2,:) .* (y(3,:)-y(2,:)); + v31 = beta(1,:) .* (x(1,:)-x(3,:)) + beta(2,:) .* (y(1,:)-y(3,:)); + elseif all(size(beta)==[nnodes,1]) + betaloc = beta(mesh.t(1:3,:)); + v12 = betaloc(2,:)-betaloc(1,:); + v23 = betaloc(3,:)-betaloc(2,:); + v31 = betaloc(1,:)-betaloc(3,:); + else + error("bim2a_advection_upwind: coefficient beta has wrong dimensions."); + endif + + [bp12, bm12] = deal (- (v12 - abs (v12))/2, (v12 + abs (v12))/2); + [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,:); + + Sloc(1,1,:) = (-bm12-bp31); + Sloc(1,2,:) = bp12; + Sloc(1,3,:) = bm31; + + Sloc(2,1,:) = bm12; + Sloc(2,2,:) = (-bp12-bm23); + Sloc(2,3,:) = bp23; + + Sloc(3,1,:) = bp31; + Sloc(3,2,:) = bm23; + Sloc(3,3,:) = (-bm31-bp23); + + A = sparse(ginode(:), gjnode(:), Sloc(:)); + + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
 [Octave-cvsupdate] SF.net SVN: octave:[8016] trunk/octave-forge/extra/bim/inst From: - 2010-12-12 09:41:38 Revision: 8016 http://octave.svn.sourceforge.net/octave/?rev=8016&view=rev Author: cdf Date: 2010-12-12 09:41:31 +0000 (Sun, 12 Dec 2010) Log Message: ----------- 3d advdiff Modified Paths: -------------- trunk/octave-forge/extra/bim/inst/bim2a_advection_diffusion.m Added Paths: ----------- trunk/octave-forge/extra/bim/inst/bim3a_advection_diffusion.m Modified: trunk/octave-forge/extra/bim/inst/bim2a_advection_diffusion.m =================================================================== --- trunk/octave-forge/extra/bim/inst/bim2a_advection_diffusion.m 2010-12-10 16:00:16 UTC (rev 8015) +++ trunk/octave-forge/extra/bim/inst/bim2a_advection_diffusion.m 2010-12-12 09:41:31 UTC (rev 8016) @@ -128,8 +128,7 @@ for jnode = 1:3 ginode(inode,jnode,:) = mesh.t(inode,:); gjnode(inode,jnode,:) = mesh.t(jnode,:); - Lloc(inode,jnode,:) = \ - sum( shg(:,inode,:) .* shg(:,jnode,:),1) .* alphaareak; + Lloc(inode,jnode,:) = sum( shg(:,inode,:) .* shg(:,jnode,:),1) .* alphaareak; endfor endfor Added: trunk/octave-forge/extra/bim/inst/bim3a_advection_diffusion.m =================================================================== --- trunk/octave-forge/extra/bim/inst/bim3a_advection_diffusion.m (rev 0) +++ trunk/octave-forge/extra/bim/inst/bim3a_advection_diffusion.m 2010-12-12 09:41:31 UTC (rev 8016) @@ -0,0 +1,125 @@ +## Copyright (C) 2010 Carlo de Falco +## +## This file is part of: +## BIM - Diffusion Advection Reaction PDE Solver +## +## BIM is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## BIM is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with BIM; If not, see ;. +## +## author: Carlo de Falco + +## -*- texinfo -*- +## +## @deftypefn {Function File} @ +## {[@var{A}]} = @ +## bim3a_advection_diffusion(@var{mesh},@var{alpha}, @var{v}) +## +## Build the Scharfetter-Gummel stabilized stiffness matrix for a +## diffusion-advection problem. +## +## The equation taken into account is: +## +## @iftex +## @tex +## $- ( \alpha ( u' - \nabla v u ))' = f$ +## @end tex +## @end iftex +## @ifinfo +## - (@var{alpha} ( u' - @var{v}' u ))' = f +## @end ifinfo +## @ifhtml +## - (@var{alpha} ( u' - @var{v}' u ))' = f +## @end ifhtml +## +## where @var{v} is a piecewise linear continuous scalar +## functions and @var{alpha} is a piecewise constant scalar function. +## +## @seealso{bim3a_rhs, bim3a_reaction, bim3a_laplacian, bim3c_mesh_properties} +## @end deftypefn + + +function SG = bim3a_advection_diffusion (mesh, acoeff, v) + + t = mesh.t; + + nelem = columns(mesh.t); + + ## Local contributions + Lloc = zeros(4,4,nelem); + + epsilonareak = reshape (acoeff .* mesh.area', 1, 1, nelem); + shg = mesh.shg(:,:,:); + + ## 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 + + vloc = v(t(:,iel)); + [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; + + Sloc(2,1,:) = bm12; + Sloc(2,2,:) = -bp12-bm23-bm24; + Sloc(2,3,:) = bp23; + Sloc(2,4,:) = bp24; + + 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; + + ## assemble global matrix + SG = sparse(ginode(:),gjnode(:),Sloc(:)); + +endfunction + This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
 [Octave-cvsupdate] SF.net SVN: octave:[8193] trunk/octave-forge/extra/bim/inst From: - 2011-03-31 12:23:40 Revision: 8193 http://octave.svn.sourceforge.net/octave/?rev=8193&view=rev Author: cdf Date: 2011-03-31 12:23:34 +0000 (Thu, 31 Mar 2011) Log Message: ----------- 3d bugs solved Modified Paths: -------------- trunk/octave-forge/extra/bim/inst/bim3a_advection_diffusion.m trunk/octave-forge/extra/bim/inst/bim3a_rhs.m Modified: trunk/octave-forge/extra/bim/inst/bim3a_advection_diffusion.m =================================================================== --- trunk/octave-forge/extra/bim/inst/bim3a_advection_diffusion.m 2011-03-29 21:00:46 UTC (rev 8192) +++ trunk/octave-forge/extra/bim/inst/bim3a_advection_diffusion.m 2011-03-31 12:23:34 UTC (rev 8193) @@ -22,7 +22,7 @@ ## ## @deftypefn {Function File} @ ## {[@var{A}]} = @ -## bim3a_advection_diffusion(@var{mesh},@var{alpha}, @var{v}) +## bim3a_advection_diffusion (@var{mesh}, @var{alpha}, @var{v}) ## ## Build the Scharfetter-Gummel stabilized stiffness matrix for a ## diffusion-advection problem. @@ -70,7 +70,7 @@ endfor endfor - vloc = v(t(:,iel)); + 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,:)); Modified: trunk/octave-forge/extra/bim/inst/bim3a_rhs.m =================================================================== --- trunk/octave-forge/extra/bim/inst/bim3a_rhs.m 2011-03-29 21:00:46 UTC (rev 8192) +++ trunk/octave-forge/extra/bim/inst/bim3a_rhs.m 2011-03-31 12:23:34 UTC (rev 8193) @@ -21,7 +21,7 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{b}]} = @ -## bim3a_rhs(@var{mesh},@var{f},@var{g}) +## bim3a_rhs (@var{mesh}, @var{f}, @var{g}) ## ## Build the finite element right-hand side of a diffusion problem ## employing mass-lumping. @@ -57,8 +57,8 @@ error("bim3a_rhs: first input is not a valid mesh structure."); endif - nnodes = length(p); - nelem = length(t); + nnodes = columns (mesh.p); + nelem = length (mesh.t); ## Turn scalar input to a vector of appropriate size if isscalar(f) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.