1## Copyright (C) 2006-2014 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## author: Matteo porro <meoo85 _AT_ users.sourceforge.net> 22## author: Emanuela Abbate <emanuela.abbate _AT_ mail.polimi.it> 23 24## -*- texinfo -*- 25## 26## @deftypefn {Function File} @ 27## {[@var{A}]} = @ 28## bim2a_axisymmetric_advection_upwind (@var{mesh}, @var{beta}) 29## 30## Build the Upwind stabilized stiffness matrix for an advection problem 31## in cylindrical coordinates with axisymmetric configuration. 32## 33## The equation taken into account is: 34## 35## 1/r * d/dr (r * @var{beta}_r u) + d/dz (@var{beta}_z u) = f 36## 37## where @var{beta} 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. 42## 43## If @var{phi} is a single scalar value @var{beta} is assumed to be 0 44## in the whole domain. 45## 46## @seealso{bim2a_axisymmetric_rhs, bim2a_axisymmetric_reaction, 47## bim2a_axisymmetric_advection_diffusion, bim2c_mesh_properties} 48## @end deftypefn 49 50function A = bim2a_axisymmetric_advection_upwind (mesh, beta) 51 52 ## Check input 53 if nargin != 2 54 error("bim2a_axisymmetric_advection_upwind: wrong number of input parameters."); 55 elseif !(isstruct(mesh) && isfield(mesh,"p") && 56 isfield (mesh,"t") && isfield(mesh,"e")) 57 error("bim2a_axisymmetric_advection_upwind: first input is not a valid mesh structure."); 58 elseif !(all(mesh.p(1,:) >= 0) || all(mesh.p(1,:) <= 0)) 59 error("bim2a_axisymmetric_advection_upwind: the input mesh cannot intersect the rotation axis r=0."); 60 endif 61 62 nnodes = columns(mesh.p); 63 nelem = columns(mesh.t); 64 65 x = abs (mesh.p(1,:)); 66 x = x(mesh.t(1:3,:)); 67 y = mesh.p(2,:); 68 y = y(mesh.t(1:3,:)); 69 70 alphaareak = reshape (mesh.area, 1, 1, nelem); 71 shg = mesh.shg(:,:,:); 72 73 ## Build local Laplacian matrix 74 Lloc = zeros(3,3,nelem); 75 76 for inode = 1:3 77 for jnode = 1:3 78 ginode(inode,jnode,:) = mesh.t(inode,:); 79 gjnode(inode,jnode,:) = mesh.t(jnode,:); 80 Lloc(inode,jnode,:) = sum( shg(:,inode,:) .* shg(:,jnode,:),1) .* alphaareak; 81 endfor 82 endfor 83 84 if all(size(beta)==1) 85 v12 = 0; 86 v23 = 0; 87 v31 = 0; 88 elseif all(size(beta)==[2,nelem]) 89 v12 = beta(1,:) .* (x(2,:)-x(1,:)) + beta(2,:) .* (y(2,:)-y(1,:)); 90 v23 = beta(1,:) .* (x(3,:)-x(2,:)) + beta(2,:) .* (y(3,:)-y(2,:)); 91 v31 = beta(1,:) .* (x(1,:)-x(3,:)) + beta(2,:) .* (y(1,:)-y(3,:)); 92 elseif all(size(beta)==[nnodes,1]) 93 betaloc = beta(mesh.t(1:3,:)); 94 v12 = betaloc(2,:)-betaloc(1,:); 95 v23 = betaloc(3,:)-betaloc(2,:); 96 v31 = betaloc(1,:)-betaloc(3,:); 97 else 98 error("bim2a_axisymmetric_advection_upwind: coefficient beta has wrong dimensions."); 99 endif 100 101 [bp12, bm12] = deal (- (v12 - abs (v12))/2, (v12 + abs (v12))/2); 102 [bp23, bm23] = deal (- (v23 - abs (v23))/2, (v23 + abs (v23))/2); 103 [bp31, bm31] = deal (- (v31 - abs (v31))/2, (v31 + abs (v31))/2); 104 105 r12 = (x(2,:) + x(1,:)) / 2; 106 r23 = (x(3,:) + x(2,:)) / 2; 107 r31 = (x(1,:) + x(3,:)) / 2; 108 109 bp12 = reshape(r12 .* bp12,1,1,nelem).*Lloc(1,2,:); 110 bm12 = reshape(r12 .* bm12,1,1,nelem).*Lloc(1,2,:); 111 bp23 = reshape(r23 .* bp23,1,1,nelem).*Lloc(2,3,:); 112 bm23 = reshape(r23 .* bm23,1,1,nelem).*Lloc(2,3,:); 113 bp31 = reshape(r31 .* bp31,1,1,nelem).*Lloc(3,1,:); 114 bm31 = reshape(r31 .* bm31,1,1,nelem).*Lloc(3,1,:); 115 116 Sloc(1,1,:) = (-bm12-bp31); 117 Sloc(1,2,:) = bp12; 118 Sloc(1,3,:) = bm31; 119 120 Sloc(2,1,:) = bm12; 121 Sloc(2,2,:) = (-bp12-bm23); 122 Sloc(2,3,:) = bp23; 123 124 Sloc(3,1,:) = bp31; 125 Sloc(3,2,:) = bm23; 126 Sloc(3,3,:) = (-bm31-bp23); 127 128 A = sparse(ginode(:), gjnode(:), Sloc(:)); 129 130endfunction 131