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