1## Copyright (C) 2006-2013  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: Matteo Porro       <meoo85 _AT_ users.sourceforge.net>
21
22## -*- texinfo -*-
23##
24## @deftypefn {Function File} {[@var{norm_u}]} = @
25## bim2c_norm(@var{mesh},@var{u},@var{norm_type})
26##
27## Compute the @var{norm_type}-norm of function @var{u} on the domain described
28## by the triangular grid @var{mesh}.
29##
30## The input function @var{u} can be either a piecewise linear conforming scalar
31## function or an elementwise constant scalar or vector function.
32##
33## The string parameter @var{norm_type} can be one among 'L2', 'H1' and 'inf'.
34##
35## Should the input function be piecewise constant, the H1 norm will not be
36## computed and the function will return an error message.
37##
38## For the numerical integration of the L2 norm the second order middle point
39## quadrature rule is used.
40##
41## @seealso{bim1c_norm, bim3c_norm}
42##
43## @end deftypefn
44
45function [norm_u] = bim2c_norm (m, u, norm_type)
46
47  ## Check input
48  if (nargin != 3)
49    error ("bim2c_norm: wrong number of input parameters.");
50  elseif (! (isstruct (m) && isfield (m,"p")
51          && isfield (m, "t")
52          && isfield (m, "e")))
53    error ("bim2c_norm: first input is not a valid mesh structure.");
54  endif
55
56  nnodes = columns (m.p);
57  nel    = columns (m.t);
58
59  if (isequal (size (u), [2, nel]))
60    u = u';
61  endif
62
63  if ((numel (u) != nnodes) && (rows (u) != nel))
64    error ("bim2c_norm: numel(u) != nnodes and rows(u) != nel.");
65  endif
66
67  if (! (strcmp (norm_type,'L2')
68         || strcmp (norm_type,'inf')
69         || strcmp (norm_type,'H1')))
70    error ("bim2c_norm: invalid norm type parameter.");
71  endif
72
73  if (strcmp (norm_type,'inf'))
74    norm_u = max (abs (u(:)));
75  else
76    if (numel (u) == nnodes)
77
78      M = __mass_matrix__ (m);
79
80      if (strcmp (norm_type, 'H1'))
81        A = bim2a_laplacian (m, 1, 1);
82        M += A;
83      endif
84
85      norm_u = sqrt(u' * M * u);
86
87    else
88
89      if (strcmp (norm_type, 'H1'))
90        error (["bim2c_norm: cannot compute the H1 norm ", ...
91                "of an elementwise constant function."]);
92      endif
93
94      norm_u = m.area' * (norm (u, 2, 'rows').^2);
95      norm_u = sqrt (norm_u);
96
97    endif
98  endif
99
100endfunction
101
102function M = __mass_matrix__ (mesh)
103
104  t      = mesh.t;
105  nnodes = columns (mesh.p);
106  nelem  = columns (t);
107
108  ## Local contributions
109
110  Mref = 1/12 * [2 1 1; 1 2 1; 1 1 2];
111  area = reshape (mesh.area, 1, 1, nelem);
112
113  ## Computation
114  for inode = 1:3
115    for jnode = 1:3
116      ginode(inode,jnode,:) = t(inode,:);
117      gjnode(inode,jnode,:) = t(jnode,:);
118    endfor
119  endfor
120  Mloc = area .* Mref;
121
122  ## assemble global matrix
123  M = sparse (ginode(:), gjnode(:), Mloc(:), nnodes, nnodes);
124
125endfunction
126
127%!test
128%!shared L, V, x, y, m
129%! L = rand (1); V = rand (1); x = linspace (0,L,4);  y = x;
130%! m = msh2m_structured_mesh (x,y,1,1:4);
131%! m.area = msh2m_geometrical_properties (m, 'area');
132%! m.shg  = msh2m_geometrical_properties (m, 'shg');
133%! u    = V * ones (columns(m.p),1);
134%! uinf = bim2c_norm (m, u, 'inf');
135%! uL2  = bim2c_norm (m, u, 'L2');
136%! uH1  = bim2c_norm (m, u, 'H1');
137%! assert ([uinf, uL2, uH1], [V, V*L, V*L], 1e-12);
138%!test
139%! u    = V * (m.p(1,:) + 2*m.p(2,:))';
140%! uinf = bim2c_norm (m, u, 'inf');
141%! uL2  = bim2c_norm (m, u, 'L2');
142%! uH1  = bim2c_norm (m, u, 'H1');
143%! assert ([uinf, uL2, uH1],
144%!         [3*L*V, V*L^2*sqrt(8/3), V*sqrt(8/3*L^4 + 5*L^2)],
145%!          1e-12);
146%!test
147%! u    = V * ones (columns(m.t),1);
148%! uinf = bim2c_norm (m, u, 'inf');
149%! uL2  = bim2c_norm (m, u, 'L2');
150%! assert ([uinf, uL2], [V, V*L], 1e-12);
151%!test
152%! u     = V * ones (columns(m.t),1);
153%! uvect = [u, 2*u];
154%! uinf  = bim2c_norm (m, uvect, 'inf');
155%! uL2   = bim2c_norm (m, uvect, 'L2');
156%! assert ([uinf, uL2], [2*V, V*L*sqrt(5)], 1e-12);
157