1## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo 2## 3## This file is part of: 4## BIM - Diffusion Advection Reaction PDE Solver 5## 6## BIM is free software; you can redistribute it and/or modify 7## it under the terms of the GNU General Public License as published by 8## the Free Software Foundation; either version 2 of the License, or 9## (at your option) any later version. 10## 11## BIM is distributed in the hope that it will be useful, 12## but WITHOUT ANY WARRANTY; without even the implied warranty of 13## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14## GNU General Public License for more details. 15## 16## You should have received a copy of the GNU General Public License 17## along with BIM; If not, see <http://www.gnu.org/licenses/>. 18## 19## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> 20## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> 21 22## -*- texinfo -*- 23## @deftypefn {Function File} {[@var{jx},@var{jy}]} = @ 24## bim2c_global_flux(@var{mesh},@var{u},@var{alpha},@var{gamma},@var{eta},@var{beta}) 25## 26## Compute the flux associated with the Scharfetter-Gummel approximation 27## of the scalar field @var{u}. 28## 29## The vector field is defined as: 30## 31## J(@var{u}) = @var{alpha}* @var{gamma} * (@var{eta} * grad @var{u} - @var{beta} * @var{u})) 32## 33## where @var{alpha} is an element-wise constant scalar function, 34## @var{eta} and @var{gamma} are piecewise linear conforming scalar 35## functions, while @var{beta} is element-wise constant vector function. 36## 37## J(@var{u}) is an element-wise constant vector function. 38## 39## Instead of passing the vector field @var{beta} directly one can pass 40## a piecewise linear conforming scalar function @var{phi} as the last 41## input. In such case @var{beta} = grad @var{phi} is assumed. If 42## @var{phi} is a single scalar value @var{beta} is assumed to be 0 in 43## the whole domain. 44## 45## @seealso{bim2c_pde_gradient,bim2a_advection_diffusion} 46## @end deftypefn 47 48function [jx, jy] = bim2c_global_flux(mesh,u,alpha,gamma,eta,beta) 49 50 ## Check input 51 if nargin != 6 52 error("bim2c_global_flux: wrong number of input parameters."); 53 elseif !(isstruct(mesh) && isfield(mesh,"p") && 54 isfield (mesh,"t") && isfield(mesh,"e")) 55 error("bim2c_global_flux: first input is not a valid mesh structure."); 56 endif 57 58 nnodes = columns(mesh.p); 59 nelem = columns(mesh.t); 60 61 62 if !( isvector(u) && isvector(alpha) && isvector(gamma) && isvector(eta) ) 63 error("bim2c_global_flux: coefficients are not valid vectors."); 64 elseif (numel (u) != nnodes) 65 error("bim2c_global_flux: length of u is not equal to the number of nodes."); 66 elseif (numel (alpha) != nelem) 67 error("bim2c_global_flux: length of alpha is not equal to the number of elements."); 68 elseif (numel (gamma) != nnodes) 69 error("bim2c_global_flux: length of gamma is not equal to the number of nodes."); 70 elseif (numel (eta) != nnodes) 71 error("bim2c_global_flux: length of eta is not equal to the number of nodes."); 72 endif 73 74 nelem = columns(mesh.t); 75 nnodes = columns(mesh.p); 76 77 uloc = u(mesh.t(1:3,:)); 78 79 shgx = reshape(mesh.shg(1,:,:),3,nelem); 80 shgy = reshape(mesh.shg(2,:,:),3,nelem); 81 82 x = reshape(mesh.p(1,mesh.t(1:3,:)),3,[]); 83 dx = [ (x(3,:)-x(2,:)) ; 84 (x(1,:)-x(3,:)) ; 85 (x(2,:)-x(1,:)) ]; 86 87 y = reshape(mesh.p(2,mesh.t(1:3,:)),3,[]); 88 dy = [ (y(3,:)-y(2,:)) ; 89 (y(1,:) -y(3,:)) ; 90 (y(2,:) -y(1,:)) ]; 91 92 if all(size(beta)==1) 93 v12=0;v23=0;v31=0; 94 elseif all(size(beta)==[2,nelem]) 95 v23 = beta(1,:) .* dx(1,:) + beta(2,:) .* dy(1,:); 96 v31 = beta(1,:) .* dx(2,:) + beta(2,:) .* dy(2,:); 97 v12 = beta(1,:) .* dx(3,:) + beta(2,:) .* dy(3,:); 98 elseif all(size(beta)==[nnodes,1]) 99 betaloc = beta(mesh.t(1:3,:)); 100 v23 = betaloc(3,:)-betaloc(2,:); 101 v31 = betaloc(1,:)-betaloc(3,:); 102 v12 = betaloc(2,:)-betaloc(1,:); 103 else 104 error("bim2c_global_flux: coefficient beta has wrong dimensions."); 105 endif 106 107 etaloc = eta(mesh.t(1:3,:)); 108 109 eta23 = etaloc(3,:)-etaloc(2,:); 110 eta31 = etaloc(1,:)-etaloc(3,:); 111 eta12 = etaloc(2,:)-etaloc(1,:); 112 113 etalocm1 = bimu_logm(etaloc(2,:),etaloc(3,:)); 114 etalocm2 = bimu_logm(etaloc(3,:),etaloc(1,:)); 115 etalocm3 = bimu_logm(etaloc(1,:),etaloc(2,:)); 116 117 gammaloc = gamma(mesh.t(1:3,:)); 118 geloc = gammaloc.*etaloc; 119 120 gelocm1 = bimu_logm(geloc(2,:),geloc(3,:)); 121 gelocm2 = bimu_logm(geloc(3,:),geloc(1,:)); 122 gelocm3 = bimu_logm(geloc(1,:),geloc(2,:)); 123 124 [bp23,bm23] = bimu_bernoulli( (v23 - eta23)./etalocm1); 125 [bp31,bm31] = bimu_bernoulli( (v31 - eta31)./etalocm2); 126 [bp12,bm12] = bimu_bernoulli( (v12 - eta12)./etalocm3); 127 128 gfigfj = [ shgx(3,:) .* shgx(2,:) + shgy(3,:) .* shgy(2,:) ; 129 shgx(1,:) .* shgx(3,:) + shgy(1,:) .* shgy(3,:) ; 130 shgx(2,:) .* shgx(1,:) + shgy(2,:) .* shgy(1,:) ]; 131 132 jx = - alpha' .* ( gelocm1 .* etalocm1 .* dx(1,:) .* ... 133 gfigfj(1,:) .* ... 134 ( bp23 .* uloc(3,:)./etaloc(3,:) -... 135 bm23 .* uloc(2,:)./etaloc(2,:)) +... %% 1 136 gelocm2 .* etalocm2 .* dx(2,:) .* ... 137 gfigfj(2,:) .* ... 138 (bp31 .* uloc(1,:)./etaloc(1,:) -... 139 bm31 .* uloc(3,:)./etaloc(3,:)) +... %% 2 140 gelocm3 .* etalocm3 .* dx(3,:) .* ... 141 gfigfj(3,:) .* ... 142 (bp12 .* uloc(2,:)./etaloc(2,:) -... 143 bm12 .* uloc(1,:)./etaloc(1,:)) ... %% 3 144 ); 145 146 jy = - alpha' .* ( gelocm1 .* etalocm1 .* dy(1,:) .* ... 147 gfigfj(1,:) .* ... 148 ( bp23 .* uloc(3,:)./etaloc(3,:) -... 149 bm23 .* uloc(2,:)./etaloc(2,:)) +... %% 1 150 gelocm2 .* etalocm2 .* dy(2,:) .* ... 151 gfigfj(2,:) .* ... 152 (bp31 .* uloc(1,:)./etaloc(1,:) -... 153 bm31 .* uloc(3,:)./etaloc(3,:)) +... %% 2 154 gelocm3 .* etalocm3 .* dy(3,:) .* ... 155 gfigfj(3,:) .* ... 156 (bp12 .* uloc(2,:)./etaloc(2,:) -... 157 bm12 .* uloc(1,:)./etaloc(1,:)) ... %% 3 158 ); 159endfunction 160