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## 24## @deftypefn {Function File} @ 25## {[@var{A}]} = @ 26## bim2a_advection_upwind (@var{mesh}, @var{beta}) 27## 28## Build the UW stabilized stiffness matrix for an advection problem. 29## 30## The equation taken into account is: 31## 32## div (@var{beta} u) = f 33## 34## where @var{beta} is an element-wise constant vector function. 35## 36## Instead of passing the vector field @var{beta} directly one can pass 37## a piecewise linear conforming scalar function @var{phi} as the last 38## input. In such case @var{beta} = grad @var{phi} is assumed. 39## 40## If @var{phi} is a single scalar value @var{beta} is assumed to be 0 41## in the whole domain. 42## 43## @seealso{bim2a_rhs, bim2a_reaction, bim2c_mesh_properties} 44## @end deftypefn 45 46function A = bim2a_advection_upwind (mesh, beta) 47 48 ## Check input 49 if nargin != 2 50 error("bim2a_advection_upwind: wrong number of input parameters."); 51 elseif !(isstruct(mesh) && isfield(mesh,"p") && 52 isfield (mesh,"t") && isfield(mesh,"e")) 53 error("bim2a_advection_upwind: first input is not a valid mesh structure."); 54 endif 55 56 nnodes = columns(mesh.p); 57 nelem = columns(mesh.t); 58 59 alphaareak = reshape (mesh.area, 1, 1, nelem); 60 shg = mesh.shg(:,:,:); 61 62 ## Build local Laplacian matrix 63 Lloc = zeros(3,3,nelem); 64 65 for inode = 1:3 66 for jnode = 1:3 67 ginode(inode,jnode,:) = mesh.t(inode,:); 68 gjnode(inode,jnode,:) = mesh.t(jnode,:); 69 Lloc(inode,jnode,:) = sum( shg(:,inode,:) .* shg(:,jnode,:),1) .* alphaareak; 70 endfor 71 endfor 72 73 x = mesh.p(1,:); 74 x = x(mesh.t(1:3,:)); 75 y = mesh.p(2,:); 76 y = y(mesh.t(1:3,:)); 77 78 if all(size(beta)==1) 79 v12 = 0; 80 v23 = 0; 81 v31 = 0; 82 elseif all(size(beta)==[2,nelem]) 83 v12 = beta(1,:) .* (x(2,:)-x(1,:)) + beta(2,:) .* (y(2,:)-y(1,:)); 84 v23 = beta(1,:) .* (x(3,:)-x(2,:)) + beta(2,:) .* (y(3,:)-y(2,:)); 85 v31 = beta(1,:) .* (x(1,:)-x(3,:)) + beta(2,:) .* (y(1,:)-y(3,:)); 86 elseif all(size(beta)==[nnodes,1]) 87 betaloc = beta(mesh.t(1:3,:)); 88 v12 = betaloc(2,:)-betaloc(1,:); 89 v23 = betaloc(3,:)-betaloc(2,:); 90 v31 = betaloc(1,:)-betaloc(3,:); 91 else 92 error("bim2a_advection_upwind: coefficient beta has wrong dimensions."); 93 endif 94 95 [bp12, bm12] = deal (- (v12 - abs (v12))/2, (v12 + abs (v12))/2); 96 [bp23, bm23] = deal (- (v23 - abs (v23))/2, (v23 + abs (v23))/2); 97 [bp31, bm31] = deal (- (v31 - abs (v31))/2, (v31 + abs (v31))/2); 98 99 bp12 = reshape(bp12,1,1,nelem).*Lloc(1,2,:); 100 bm12 = reshape(bm12,1,1,nelem).*Lloc(1,2,:); 101 bp23 = reshape(bp23,1,1,nelem).*Lloc(2,3,:); 102 bm23 = reshape(bm23,1,1,nelem).*Lloc(2,3,:); 103 bp31 = reshape(bp31,1,1,nelem).*Lloc(3,1,:); 104 bm31 = reshape(bm31,1,1,nelem).*Lloc(3,1,:); 105 106 Sloc(1,1,:) = (-bm12-bp31); 107 Sloc(1,2,:) = bp12; 108 Sloc(1,3,:) = bm31; 109 110 Sloc(2,1,:) = bm12; 111 Sloc(2,2,:) = (-bp12-bm23); 112 Sloc(2,3,:) = bp23; 113 114 Sloc(3,1,:) = bp31; 115 Sloc(3,2,:) = bm23; 116 Sloc(3,3,:) = (-bm31-bp23); 117 118 A = sparse(ginode(:), gjnode(:), Sloc(:)); 119 120 121endfunction 122